VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA302 10.15559/26-VMSTA302 Research Article On mixed time-changed Erlang queue https://orcid.org/0009-0004-5178-3305 PoteRohini Bhagwanraorohinib@iitbhilai.ac.in∗ KatariaKuldeep Kumarkuldeepk@iitbhilai.ac.in Indian Institute of Technology Bhilai, Durg 491002, India Corresponding author. 2026 96202600127218202527320262352026 © 2026 The Author(s). Published by VTeX2026 Open access article under the CC BY license.

We study a time-changed variant of the Erlang queue by taking the first hitting time of a mixed stable subordinator as the time-changing component. We call it the mixed time-changed Erlang queue. We derive the system of fractional differential equations that governs its state probabilities. The explicit expressions for the state probabilities of mixed time-changed Erlang queue and their Laplace transform are derived. Also, an equivalent representation of this time-changed queue in terms of phases is provided, and its mean queue length is obtained. Some of its distributional properties such as the distribution of its inter-arrival times, inter-phase times, service times and busy period are derived. Later, its conditional waiting time is discussed and some plots of sample paths simulation are presented.

 Erlang service distribution state probabilities inverse mixed stable subordinator Mittag-Leffler function 60K15 60K20 60K25 26A33 The research of first author was supported by a UGC fellowship, NTA reference no. 231610158041, Govt. of India. Introduction

The Erlang queue is a queueing system in which the arrival of customers is modeled according to a Poisson process with rate λ and its service system has Erlang distribution with shape parameter k and mean 1 / μ. Its inter-arrival times are exponentially distributed with parameter λ and its service system consists of k phases, where the service time of each phase is exponentially distributed with parameter k μ. In the field of telecommunication, the Erlang queues are quite useful. In call centers, the Erlang queue is used to predict number of agents required to handle incoming calls (see [29]).

In [21, 22] a single channel queue with Poisson arrivals and a general class of service-time distributions is studied. The transient solution of Erlang queue is obtained in [16]. Cahoy et al. [7] proposed the first fractional generalization of M / M / 1 queue where they obtained its state probabilities and an algorithm to simulate fractional M / M / 1 queue. Di Crescenzo et al. [9] studied M / M / 1 queue with catastrophes and its fractional variant is discussed in [3]. Giorno et al. [14] studied a single-server queueing system with Poisson arrivals and state-dependent service mechanism characterized by logarithmic steady-state distribution. The Erlang queue time-changed with inverse stable subordinator is studied by Ascione et al. [4], where its state probabilities, mean queue length along with the distributions its of busy period, inter-arrival, inter-phase and service times are derived.

Let N denote the set of positive integers. The Erlang queue { Q ( t ) } t ≥ 0 with state space H 0 = H ∪ { ( 0 , 0 ) }, where H = { ( n , s ) ∈ N × N : s ≤ k } can be described as follows: Q ( t ) = ( N ( t ) , S ( t ) ). Here, N ( t ) denotes the number of customers in the system at time t ≥ 0 and S ( t ) is the phase of customer being served at time t provided N ( t ) > 0. For N ( t ) = 0, we take S ( t ) = 0.

For t ≥ 0, let us denote the transient state probabilities of Erlang queue as follows: p 0 ( t ) = Pr ( Q ( t ) = ( 0 , 0 ) | Q ( 0 ) = ( 0 , 0 ) ) and p n , s ( t ) = Pr ( Q ( t ) = ( n , s ) | Q ( 0 ) = ( 0 , 0 ) ), ( n , s ) ∈ H. That is, p 0 ( t ) is the probability that there are no customers in the system at time t and p n , s ( t ) is the probability that there are n customers in the system at time t and the customer being served is at phase s.

Equivalently, the Erlang queue can be represented as a queue length process in terms of phase count (see [4], p. 3252). Consider a bijective map g k : H 0 → N 0 defined as follows: g k ( n , s ) : = k ( n − 1 ) + s , ( n , s ) ∈ H , 0 , ( n , s ) = ( 0 , 0 ) , where N 0 = N ∪ { 0 }. Its inverse map ( a k ( m ) , b k ( m ) ) is such that b k ( m ) = min { s > 0 : s ≡ m ( mod k ) } , m > 0 , 0 , m = 0 , and a k ( m ) = m − b k ( m ) k + 1 , m > 0 , 0 , m = 0 . Then, L ( t ) = g k ( Q ( t ) ) is the length of Erlang queue in terms of number of phases at time t ≥ 0. Its state probabilities are denoted by p n ( t ) = Pr ( L ( t ) = n | L ( 0 ) = 0 ), n ≥ 0. Note that p n ( t ) = p a k ( n ) , b k ( n ) ( t ).

It is well known that the occurrence of catastrophic events such as earthquakes, tsunami, etc. exhibits long memory. It is empirically proven that these real life situations can be modeled by time-changed processes as they exhibits long range behaviour. Recently, the time-changed processes are extensively studied by many researchers, for example, the time fractional Poisson process (see [6, 23]), space fractional Poisson process (see [27]), time-changed birth-death processes (see [19, 26]), etc. Beghin [5] studied random-time processes governed by differential equations of fractional distributed order, Kataria and Khandakar [18] studied mixed fractional risk process which are based on the inverse mixed stable subordinator. Recently, Pote and Kataria [28] studied a time-changed variant of the Erlang queue with multiple arrivals where the time-changing component used is the inverse stable subordinator. For some more recent literature in this direction, we refer the reader to the works on multiparameter generalized counting process and its time-changed variants (see [10]), tempered space-time fractional negative binomial process (see [11]), bivariate gamma subordination for a Poisson shock model (see [25]), non-homogeneous generalized fractional Skellam process (see [30]), etc.

In a queue time-changed with inverse stable subordinator, the waiting times depends on a single fractional parameter, and hence has a fixed memory scale (see [4]). However, whenever a process is time-changed via an inverse mixed stable subordinator it mixes different waiting times. This makes it potentially applicable for modeling systems with distributed-order fractional dynamics, which cannot be captured by an inverse stable subordinator (see [5]).

In this paper, we introduce and study a time-changed variant of the Erlang queue where the time changing component used is the first hitting time of a mixed stable subordinator. We call it the mixed time-changed Erlang queue. The paper is organized as follows:

First, we give some preliminary results and known definitions of special functions, fractional derivative, subordinators and their inverse processes, etc. Then, we define the mixed time-changed Erlang queue by time-changing it with an inverse mixed stable subordinator which is the first hitting time of a mixed stable subordinator. The system of fractional differential equations that governs its state probabilities is derived. The Laplace transform of its state probabilities is obtained whose inversion yields its distribution. Also, the mixed time-changed Erlang queue is represented as a queue length process which characterizes it in terms of number of phases. The fractional differential equation that governs its mean queue length is derived, and an explicit expression for the mean queue length is obtained. Moreover, we derive the distribution of inter-arrival times, inter-phase times, service times and busy period of mixed time-changed Erlang queue, and discuss its conditional waiting time. Some plots of sample paths simulation are given.

 Preliminaries

Here, we collect some known definitions and results that will be used.

Mittag-Leffler function

The three-parameter Mittag-Leffler function also known as the Prabhakar function is defined as (see [20]) E α , β γ ( t ) : = ∑ r = 0 ∞ Γ ( γ + r ) t r r ! Γ ( γ ) Γ ( α r + β ) , t ∈ R , α , β , γ > 0 , where R denotes the set of real numbers.

The following result is useful in the theory of stable subordinators (see [17]): L ( t β − 1 E α , β γ ( ω t α ) ) ( z ) = z α γ − β ( z α − ω ) γ , ω ∈ R , z > 0 , | ω z − α | < 1 .

The following result will be used (see [17], Eq. (17.6)): L − 1 ( z ρ − 1 z α + a z β + b ; t ) = t α − ρ ∑ h = 0 ∞ ( − a ) h t ( α − β ) h E α , α + ( α − β ) h − ρ + 1 h + 1 ( − b t α ) , where α > β > 0, α − ρ + 1 > 0 and | a z β / ( z α + b ) | < 1.

 Fractional derivative

The Caputo fractional derivative of an absolutely continuous function g ( · ) is defined as follows (see [20]): d α d t α g ( t ) = 1 Γ ( 1 − α ) ∫ 0 t ( t − y ) − α g ′ ( y ) d y , 0 < α < 1 , g ′ ( t ) , α = 1 . Its Laplace transform is given by L ( d α d t α g ( t ) ) ( z ) = z α g ˜ ( z ) − z α − 1 g ( 0 ) , where g ˜ ( z ) denotes the Laplace transform of the function g ( t ).

 Stable subordinators and their inverse

An α-stable subordinator is a one dimensional Lévy process { D α ( t ) } t ≥ 0 , 0 < α < 1 whose Laplace transform is given by E ( e − z D α ( t ) ) = e − t z α , z > 0 (see [2]). Its first hitting time process { Y α ( t ) } t ≥ 0 , where Y α ( t ) = inf { s ≥ 0 : D α ( s ) > t } is known as the inverse α-stable subordinator. The Laplace transform of its probability density function is given by (see [24]) L ( Pr ( Y α ( t ) ∈ d y ) ) ( z ) = z α − 1 e − y z α d y , z > 0 .

The mixed stable subordinator { D α 1 , α 2 ( t ) } t ≥ 0 is a Lévy process whose Laplace transform is given by E ( e − z D α 1 , α 2 ( t ) ) = e − t ( c 1 z α 1 + c 2 z α 2 ) , z > 0 , where c 1 ≥ 0, c 2 ≥ 0, c 1 + c 2 = 1, 0 < α 2 < α 1 < 1. Its first hitting time process { Y α 1 , α 2 ( t ) } t ≥ 0 is known as the inverse mixed stable subordinator, where Y α 1 , α 2 ( t ) = inf { s ≥ 0 : D α 1 , α 2 ( s ) > t } , t > 0 with Y α 1 , α 2 ( 0 ) = 0. For c 1 , c 2 ∈ { 0 , 1 }, it reduces to inverse stable subordinator. The probability density function f α 1 , α 2 ( t , u ) = d d u Pr ( Y α 1 , α 2 ( t ) ≤ u ) of inverse mixed stable subordinator has the following Laplace transform (see [1]): ∫ 0 ∞ e − z t f α 1 , α 2 ( t , u ) d t = c 1 z α 1 + c 2 z α 2 z e − u ( c 1 z α 1 + c 2 z α 2 ) , z > 0 .

Mixed time-changed Erlang queue

In this section, we introduce a time-changed Erlang queue where the time changing component is the first hitting time of mixed stable subordinator. We call it the mixed time-changed Erlang queue.

Let { Q ( t ) } t ≥ 0 be an Erlang queue, that is, Q ( t ) = ( N ( t ) , S ( t ) ), t ≥ 0, where N ( t ) denotes the number of customers in the system at time t and S ( t ) denotes the phase of customer being served at time t. We define the mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 , 0 < α 2 < α 1 < 1 as follows: Q α 1 , α 2 ( t ) : = Q ( Y α 1 , α 2 ( t ) ) , where the Erlang queue { Q ( t ) } t ≥ 0 is independent of the inverse mixed stable subordinator { Y α 1 , α 2 ( t ) } t ≥ 0 . For t ≥ 0, let us denote its state probabilities by p 0 α 1 , α 2 ( t ) = Pr ( Q α 1 , α 2 ( t ) = ( 0 , 0 ) | Q α 1 , α 2 ( 0 ) = ( 0 , 0 ) ) , p n , s α 1 , α 2 ( t ) = Pr ( Q α 1 , α 2 ( t ) = ( n , s ) | Q α 1 , α 2 ( 0 ) = ( 0 , 0 ) ) , ( n , s ) ∈ H , where H = { ( n , s ) ∈ N × N : s ≤ k }.

In the next result, we obtain the system of fractional differential equations whose solution gives the state probabilities of mixed time-changed Erlang queue.

Let c 1 ≥ 0, c 2 ≥ 0 such that c 1 + c 2 = 1 and d α d t α be the Caputo fractional derivative of order α as defined in (7). Then, the state probabilities of mixed time-changed Erlang queue solve the following system of differential equations: ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p 0 α 1 , α 2 ( t ) = − λ p 0 α 1 , α 2 ( t ) + k μ p 1 , 1 α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p 1 , s α 1 , α 2 ( t ) = − ( λ + k μ ) p 1 , s α 1 , α 2 ( t ) + k μ p 1 , s + 1 α 1 , α 2 ( t ) , 1 ≤ s ≤ k − 1 , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p 1 , k α 1 , α 2 ( t ) = − ( λ + k μ ) p 1 , k α 1 , α 2 ( t ) + k μ p 2 , 1 α 1 , α 2 ( t ) + λ p 0 α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p n , s α 1 , α 2 ( t ) = − ( λ + k μ ) p n , s α 1 , α 2 ( t ) + k μ p n , s + 1 α 1 , α 2 ( t ) + λ p n − 1 , s α 1 , α 2 ( t ) , 1 ≤ s ≤ k − 1 , n ≥ 2 , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p n , k α 1 , α 2 ( t ) = − ( λ + k μ ) p n , k α 1 , α 2 ( t ) + k μ p n + 1 , 1 α 1 , α 2 ( t ) + λ p n − 1 , k α 1 , α 2 ( t ) , n ≥ 2 , with p 0 α 1 , α 2 ( 0 ) = 1 and p n , s α 1 , α 2 ( 0 ) = 0 for 1 ≤ s ≤ k, n ≥ 1.

For ( n , s ) ∈ H, we have p n , s α 1 , α 2 ( t ) = Pr ( Q ( Y α 1 , α 2 ( t ) ) = ( n , s ) | Q ( 0 ) = ( 0 , 0 ) ) = ∫ 0 ∞ Pr ( Q ( y ) = ( n , s ) | Q ( 0 ) = ( 0 , 0 ) ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) = ∫ 0 ∞ p n , s ( y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) , t ≥ 0 and p 0 α 1 , α 2 ( t ) = ∫ 0 ∞ p 0 ( y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) , t ≥ 0 where p n , s ( y ) and p 0 ( y ) are the state probability and zero state probability of classical Erlang queue, respectively.

Let the generating function of { Q α 1 , α 2 ( t ) } t ≥ 0 be given by G α 1 , α 2 ( x , t ) = p 0 α 1 , α 2 ( t ) + ∑ n = 1 ∞ ∑ s = 1 k p n , s α 1 , α 2 ( t ) x ( n − 1 ) k + s , t ≥ 0 , | x | ≤ 1 . On substituting (11) and (12) in (13), we get G α 1 , α 2 ( x , t ) = ∫ 0 ∞ G ( x , y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) , t ≥ 0 , where G ( x , t ) is the generating function of { Q ( t ) } t ≥ 0 . On taking the Laplace transform of (12) and (14), and using (10), we get p 0 ˜ α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ p 0 ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 and G ˜ α 1 , α 2 ( x , z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ G ( x , y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 .

On multiplying the first equation by x, the second by x s + 1 , the third by x k + 1 , the fourth by x k ( n − 1 ) + s + 1 and the fifth by x n k + 1 in the system of differential equations given in Theorem 1, and then summing these equations over the range of n and s, it can be established that the state probabilities of { Q α 1 , α 2 ( t ) } t ≥ 0 are solution of the system of differential equations given in Theorem 1 if and only if its generating function is the solution of following equation: x ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) G α 1 , α 2 ( x , t ) = ( 1 − x ) ( ( k μ − λ ( x + ⋯ + x k ) ) G α 1 , α 2 ( x , t ) − k μ p 0 α 1 , α 2 ( t ) ) , with initial condition G α 1 , α 2 ( x , 0 ) = 1. On taking the Laplace transform on both sides of (17) and using (8), we obtain x ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( z G ˜ α 1 , α 2 ( x , z ) − 1 ) = ( 1 − x ) ( ( k μ − λ ( x + ⋯ + x k ) ) G ˜ α 1 , α 2 ( x , z ) − k μ p 0 ˜ α 1 , α 2 ( z ) ) , z > 0 . By using (15) and (16) in (18), we have x ( ( c 1 z α 1 + c 2 z α 2 ) ∫ 0 ∞ G ( x , y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y − 1 ) = ∫ 0 ∞ ( 1 − x ) ( ( k μ − λ ( x + ⋯ + x k ) ) G ( x , y ) − k μ p 0 ( y ) ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 . Now, on substituting Eq. (6) of [16] in (19), we get the following identity: x ( ( c 1 z α 1 + c 2 z α 2 ) ∫ 0 ∞ G ( x , y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y − 1 ) = x ∫ 0 ∞ ∂ ∂ y G ( x , y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 . This proves the required result.  □

In Theorem 1, the Caputo fractional derivatives can be written as follows: ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) f ( t ) = ∫ 0 1 d α d t α f ( t ) c ( α ) d α , where c ( α ) = c 1 δ ( α − α 1 ) + c 2 δ ( α − α 2 ), c 1 ≥ 0, c 2 ≥ 0, c 1 + c 2 = 1 and 0 < α 2 < α 1 < 1. So, the fractional derivative depends on the entire time span [ 0 , t ] via the weight function involved in the definition of Caputo derivative given in (7). Also, the memory depends on the two fractional parameters α 1 and α 2 . For instance, with the increase or decrease in the values of α 1 and α 2 , the memory effects can increase or decrease.

In the following result, we obtain the Laplace transform of zero state probability of mixed time-changed Erlang queue which we use to derive its explicit expression:

The Laplace transform of zero state probability of { Q α 1 , α 2 ( t ) } t ≥ 0 is given by p ˜ 0 α 1 , α 2 ( z ) = ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) a m , r 0 , z > 0 where a m , r 0 = m + r ( k + 1 ) and A m , r 0 = m λ r ( k μ ) m + r k − 1 r ! ( m + r k ) ! ( m + r ( k + 1 ) − 1 ) ! .

By using Eq. (5) of [4] in (15), we obtain p 0 ˜ α 1 , α 2 ( z ) = ∑ m = 1 ∞ ∑ r = 0 ∞ m λ r ( k μ ) m + r k − 1 r ! ( m + r k ) ! ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ y m + r ( k + 1 ) − 1 · e − ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) y d y = ∑ m = 1 ∞ ∑ r = 0 ∞ m λ r ( k μ ) m + r k − 1 r ! ( m + r k ) ! ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( m + r ( k + 1 ) − 1 ) ! ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) m + r ( k + 1 ) , z > 0 . On substituting (21) and (22) in (23), we get the required result.  □

For a positive integer N, let us define f N ( t ) = t α 1 − ( α 1 − 1 N ) − 1 ∑ h = 0 ∞ ( − c 2 c 1 ) h t ( α 1 − α 2 ) h E α 1 , α 1 + ( α 1 − α 2 ) h + 1 − α 1 N h + 1 ( − ( λ + k μ c 1 ) t α 1 ) and g N ( t ) = t α 1 − ( α 2 − 1 N ) − 1 ∑ h = 0 ∞ ( − c 2 c 1 ) h t ( α 1 − α 2 ) h E α 1 , α 1 + ( α 1 − α 2 ) h + 1 − α 2 N h + 1 ( − ( λ + k μ c 1 ) t α 1 ) .

The next result gives the probability of the system being empty at time t.

Let f ∗ ( n ) denote the n-fold convolution of any function f ( · ) and E α , β γ ( · ) be the three parameter Mittag-Leffler function defined in (4). Then, the zero state probability of { Q α 1 , α 2 ( t ) } t ≥ 0 is given by p 0 α 1 , α 2 ( t ) = ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 c 1 a m , r 0 ( c 1 f a m , r 0 ∗ ( a m , r 0 ) ( t ) + c 2 g a m , r 0 ∗ ( a m , r 0 ) ( t ) ) , t ≥ 0 , where a m , r 0 , A m , r 0 , f a m , r 0 ( · ) and g a m , r 0 ( · ) are given in (21), (22), (24), and (25) respectively.

On taking the inverse Laplace transform of (20), we obtain p 0 α 1 , α 2 ( t ) = ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 ( 1 c 1 a m , r 0 − 1 L − 1 ( ( z α 1 − 1 a m , r 0 z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) a m , r 0 ; t ) + c 2 c 1 a m , r 0 L − 1 ( ( z α 2 − 1 a m , r 0 z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) a m , r 0 ; t ) ) . On using (6) in (27), we get the required result.  □

In (26), if c 1 = 1 then c 2 = 0. So, p 0 α 1 , α 2 ( t ) = ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 f a m , r 0 ∗ ( a m , r 0 ) ( t ) , t ≥ 0 and f a m , r 0 ( t ) = t α 1 − ( α 1 − 1 a m , r 0 ) − 1 E α 1 , α 1 + 1 − α 1 a m , r 0 1 ( − ( λ + k μ ) t α 1 ) , t ≥ 0 . On taking the Laplace transform of (29) and using (5), we have f ˜ a m , r 0 ( z ) = z α 1 − 1 a m , r 0 z α 1 + λ + k μ , z > 0 . Thus, the Laplace transform of a 0 m , r -fold convolution of (30) is given by f ˜ a m , r 0 ∗ ( a m , r 0 ) ( z ) = z α 1 − 1 ( z α 1 + λ + k μ ) a m , r 0 , z > 0 . In (28), on taking the inverse Laplace transform of (31) and by using (5), we get the zero state probability of fractional Erlang queue (see [4], Eq. (36)). Similar result holds if c 2 = 1.

In the subsequent results, we will be using the following notations: B i n , s = λ n + i ( k μ ) k ( i + 1 ) − s ( k ( i + 1 ) − s ) ! ( n + i ) ! ( n + k − s + i ( k + 1 ) ) ! , δ i n , s = n − s + ( i + 1 ) ( k + 1 ) , C i , m , r n , s = k μ B i n , s A m , r 0 , ν i , m , r n , s = δ i n , s + a m , r 0 , D i , m , r n , s = k μ B i n , s + 1 A m , r 0 , s = 1 , 2 , … , k − 1 , k μ B i n + 1 , 1 A m , r 0 , s = k , π i , m , r n , s = δ i n , s + 1 + a m , r 0 , s = 1 , 2 , … , k − 1 , δ i n + 1 , 1 + a m , r 0 , s = k , where a m , r 0 and A m , r 0 are as given in (21) and (22), respectively.

In the following result, we obtain the Laplace transform of state probabilities of mixed time-changed Erlang queue which we use to derive its explicit expressions.

The Laplace transform of the state probabilities of { Q α 1 , α 2 ( t ) } t ≥ 0 is given by p ˜ n , s α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( ∑ i = 0 ∞ B i n , s ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) δ i n , s + ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ C i , m , r n , s ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) ν i , m , r n , s − ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ D i , m , r n , s ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) π i , m , r n , s ) , n ≥ 1 , 1 ≤ s ≤ k , z > 0 .

On taking the Laplace transform of (11) and using (10), we get p ˜ n , s α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ p n , s ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 . For n ≥ 1 and 1 ≤ s ≤ k − 1, by using Eq. (6) of [4] in (34), we obtain p ˜ n , s α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( ∑ i = 0 ∞ λ n + i ( k μ ) k ( i + 1 ) − s ( k ( i + 1 ) − s ) ! ( n + i ) ! ∫ 0 ∞ y n + k − s + i ( k + 1 ) · e − ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) y d y + ∑ i = 0 ∞ λ n + i ( k μ ) k ( i + 1 ) − s + 1 ( k ( i + 1 ) − s ) ! ( n + i ) ! · ∫ 0 ∞ ∫ 0 y p 0 ( u ) ( y − u ) n + k − s + i ( k + 1 ) e − ( λ + k μ ) ( y − u ) e − ( c 1 z α 1 + c 2 z α 2 ) y d u d y − ∑ i = 0 ∞ λ n + i ( k μ ) k ( i + 1 ) − s ( k ( i + 1 ) − s − 1 ) ! ( n + i ) ! ∫ 0 ∞ ∫ 0 y p 0 ( u ) ( y − u ) n + k − s − 1 + i ( k + 1 ) · e − ( λ + k μ ) ( y − u ) e − ( c 1 z α 1 + c 2 z α 2 ) y d u d y ) , z > 0 . Also, by using Eq. (7) of [4] in (34), we have p ˜ n , k α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( ∑ i = 0 ∞ λ n + i ( k μ ) k i ( k i ) ! ( n + i ) ! ∫ 0 ∞ y n + i ( k + 1 ) e − ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) y d y + ∑ i = 0 ∞ λ n + i ( k μ ) k i + 1 ( k i ) ! ( n + i ) ! ∫ 0 ∞ ∫ 0 y p 0 ( u ) ( y − u ) n + i ( k + 1 ) e − ( λ + k μ ) ( y − u ) · e − ( c 1 z α 1 + c 2 z α 2 ) y d u d y − ∑ i = 0 ∞ λ n + i + 1 ( k μ ) k ( i + 1 ) ( k ( i + 1 ) − 1 ) ! ( n + i + 1 ) ! · ∫ 0 ∞ ∫ 0 y p 0 ( u ) ( y − u ) n + k + i ( k + 1 ) e − ( λ + k μ ) ( y − u ) e − ( c 1 z α 1 + c 2 z α 2 ) y d u d y ) , z > 0 . By using Lemma 5.4 of [4] and (32) in (35) and (36), we get the required result.  □

For c 1 = 1 or c 2 = 1 in (33), the Laplace transform of state probabilities of mixed time-changed Erlang queue reduces to that of fractional Erlang queue (see [4], Theorem 5.5). Observe that the state probabilities of fractional Erlang queue have dependence only on one fractional parameter, whereas that of mixed time-changed Erlang queue have dependence on two distinct fractional parameters. That is, past memory effect is due to two distinct parameters in mixed time-changed Erlang queue.

Next, we obtain the probability that there are n ≥ 1 customers in the system and the customer being served is in phase s ∈ { 1 , 2 , … , k } at time t.

For 1 ≤ s ≤ k, the state probabilities of { Q α 1 , α 2 ( t ) } t ≥ 0 are given by p n , s α 1 , α 2 ( t ) = ∑ i = 0 ∞ B i n , s c 1 δ i n , s ( c 1 f δ i n , s ∗ ( δ i n , s ) ( t ) + c 2 g δ i n , s ∗ ( δ i n , s ) ( t ) ) + ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ C i , m , r n , s c 1 ν i , m , r n , s ( c 1 f ν i , m , r n , s ∗ ( ν i , m , r n , s ) ( t ) + c 2 g ν i , m , r n , s ∗ ( ν i , m , r n , s ) ( t ) ) − ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ D i , m , r n , s c 1 π i , m , r n , s ( c 1 f π i , m , r n , s ∗ ( π i , m , r n , s ) ( t ) + c 2 g π i , m , r n , s ∗ ( π i , m , r n , s ) ( t ) ) , n ≥ 1 , t ≥ 0 , where f N ( · ) and g N ( · ) are given in (24) and (25), respectively.

On taking the inverse Laplace transform of (33), we obtain p n , s α 1 , α 2 ( t ) = ∑ i = 0 ∞ B i n , s ( 1 c 1 δ i n , s − 1 L − 1 ( ( z α 1 − 1 δ i n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) δ i n , s ; t ) + c 2 c 1 δ i n , s L − 1 ( ( z α 2 − 1 δ i n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) δ i n , s ; t ) ) + ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ C i , m , r n , s ( 1 c 1 ν i , m , r n , s − 1 L − 1 ( ( z α 1 − 1 ν i , m , r n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) ν i , m , r n , s ; t ) + c 2 c 1 ν i , m , r n , s L − 1 ( ( z α 2 − 1 ν i , m , r n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) ν i , m , r n , s ; t ) ) − ∑ i = 0 ∞ ∑ m = 1 ∞ ∑ r = 0 ∞ D i , m , r n , s ( 1 c 1 π i , m , r n , s − 1 L − 1 ( ( z α 1 − 1 π i , m , r n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) π i , m , r n , s ; t ) + c 2 c 1 π i , m , r n , s L − 1 ( ( z α 2 − 1 π i , m , r n , s z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) π i , m , r n , s ; t ) ) , t ≥ 0 . On using (6) in (37), we get the required result.  □

Queue length process

Here, we define the queue length process of mixed time-changed Erlang queue as follows: L α 1 , α 2 ( t ) : = g k ( Q α 1 , α 2 ( t ) ) , t ≥ 0 , where g k ( · ) is a bijective map as defined in (1). Let p n α 1 , α 2 ( t ) = Pr ( L α 1 , α 2 ( t ) = n | L α 1 , α 2 ( 0 ) = 0 ), n ≥ 0 be its state probabilities such that p − n α 1 , α 2 ( t ) = 0 for all n ≥ 1.

The mean queue length of mixed time-changed Erlang queue is defined as M α 1 , α 2 ( t ) : = E ( L α 1 , α 2 ( t ) | L α 1 , α 2 ( 0 ) = 0 ) , t ≥ 0 .

As g k ( · ) is a bijective map, the mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 is equivalently represented by the queue length process { L α 1 , α 2 ( t ) } t ≥ 0 . Thus, p n , s α 1 , α 2 ( t ) = p g k ( n , s ) α 1 , α 2 ( t ) , t ≥ 0 and p n α 1 , α 2 ( t ) = p a k ( n ) , b k ( n ) α 1 , α 2 ( t ) , t ≥ 0 , where a k ( · ) and b k ( · ) are as given in (3) and (2), respectively.

Now, by using (39) in Theorem 1, we get the following system of differential equations that governs the state probabilities of { L α 1 , α 2 ( t ) } t ≥ 0 : ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p 0 α 1 , α 2 ( t ) = − λ p 0 α 1 , α 2 ( t ) + k μ p 1 α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p n α 1 , α 2 ( t ) = − ( λ + k μ ) p n α 1 , α 2 ( t ) + k μ p n + 1 α 1 , α 2 ( t ) , 1 ≤ n ≤ k − 1 , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p n α 1 , α 2 ( t ) = − ( λ + k μ ) p n α 1 , α 2 ( t ) + k μ p n + 1 α 1 , α 2 ( t ) + λ p n − k α 1 , α 2 ( t ) , n ≥ k with initial condition p 0 α 1 , α 2 ( 0 ) = 1. Here, d α i d t α i is a Caputo fractional derivative of order α i defined in (7).

Next, we obtain the expected number of remaining phases to be served.

The mean queue length of mixed time-changed Erlang queue solves the following fractional differential equation: ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) M α 1 , α 2 ( t ) = k ( λ − μ ) + k μ p 0 α 1 , α 2 ( t ) , t ≥ 0 with initial condition M α 1 , α 2 ( 0 ) = 0.

From (38), we have M α 1 , α 2 ( t ) = ∑ n = 0 ∞ n p n α 1 , α 2 ( t ) , t ≥ 0 . Note that M α 1 , α 2 ( 0 ) = 0 as p 0 α 1 , α 2 ( 0 ) = 1. By using (41) and (43), we get ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) M α 1 , α 2 ( t ) = − ( λ + k μ ) M α 1 , α 2 ( t ) + k μ ∑ n = 1 ∞ n p n + 1 α 1 , α 2 ( t ) + λ ∑ n = k ∞ n p n − k α 1 , α 2 ( t ) = − ( λ + k μ ) M α 1 , α 2 ( t ) + k μ ( M α 1 , α 2 ( t ) + p 0 α 1 , α 2 ( t ) − 1 ) + λ ( M α 1 , α 2 ( t ) + k ) = k ( λ − μ ) + k μ p 0 α 1 , α 2 ( t ) . This proves the result.  □

For any positive integer N, let us define the following function: h N ( t ) = t α 1 − ( N − 1 N ) ∑ j = 0 ∞ ( − c 2 c 1 ) j t ( α 1 − α 2 ) j E α 1 , α 1 + ( α 1 − α 2 ) j + 1 N j + 1 ( − ( λ + k μ c 1 ) t α 1 ) .

Let E α , β γ ( · ) be the three parameter Mittag-Leffler function as defined in (4) and h a m , r 0 ( · ) be as defined in (45). Then, the mean queue length of mixed time-changed Erlang queue is given by M α 1 , α 2 ( t ) = k ( λ − μ ) c 1 t α 1 E α 1 − α 2 , α 1 + 1 1 ( − c 2 c 1 t α 1 − α 2 ) + k μ ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 c 1 a m , r 0 h a m , r 0 ∗ ( a m , r 0 ) ( t ) , z > 0 where a m , r 0 and A m , r 0 are given in (21) and (22), respectively.

By using (11) and (40) in (43), we get M α 1 , α 2 ( t ) = ∫ 0 ∞ M ( y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) = k ( λ − μ ) ∫ 0 ∞ y Pr ( Y α 1 , α 2 ( t ) ∈ d y ) + k μ ∫ 0 ∞ ∫ 0 y p 0 ( u ) d u Pr ( Y α 1 , α 2 ( t ) ∈ d y ) , t ≥ 0 , where the second equality follows from Eq. (21) of [22]. Here, p 0 ( u ) is the zero state probability of classical Erlang queue. On taking the Laplace transform on both sides of (47) and using (10), we get M ˜ α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( k ( λ − μ ) ∫ 0 ∞ y e − y ( c 1 z α 1 + c 2 z α 2 ) d y + k μ ∫ 0 ∞ ∫ 0 y p 0 ( u ) e − y ( c 1 z α 1 + c 2 z α 2 ) d u d y ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( k ( λ − μ ) ( c 1 z α 1 + c 2 z α 2 ) 2 + k μ ∫ 0 ∞ p 0 ( u ) · ∫ u ∞ e − y ( c 1 z α 1 + c 2 z α 2 ) d y d u ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( k ( λ − μ ) ( c 1 z α 1 + c 2 z α 2 ) 2 + k μ ( c 1 z α 1 + c 2 z α 2 ) · ∫ 0 ∞ p 0 ( u ) e − u ( c 1 z α 1 + c 2 z α 2 ) d u ) = k ( λ − μ ) c 1 z α 1 + 1 + c 2 z α 2 + 1 + k μ z ∫ 0 ∞ p 0 ( u ) e − u ( c 1 z α 1 + c 2 z α 2 ) d u = k ( λ − μ ) c 1 z α 1 + 1 + c 2 z α 2 + 1 + k μ z p 0 ˜ α 1 , α 2 ( z ) ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( by using (15) ) = k ( λ − μ ) c 1 z α 1 + 1 + c 2 z α 2 + 1 + k μ z ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) a m , r 0 , z > 0 , where we have used (20) in the last step.

On taking the inverse Laplace transform on both sides of (48), we get M α 1 , α 2 ( t ) = k ( λ − μ ) c 1 L − 1 ( z − α 2 − 1 z α 1 − α 2 + c 2 c 1 ; t ) + k μ ∑ m = 1 ∞ ∑ r = 0 ∞ A m , r 0 c 1 a m , r 0 L − 1 ( ( z − 1 a m , r 0 z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) a m , r 0 ; t ) . From (5) and (6), the required result follows.  □

 Inter-arrival and inter-phase times of mixed time changed Erlang queue

Here, we view the mixed time-changed Erlang queue in terms of its inter-arrival and service times.

Let us denote the following: q n α 1 , α 2 ( t ) = Pr ( N α 1 , α 2 ( t ) = n | Q α 1 , α 2 ( 0 ) = ( 0 , 0 ) ) , n ≥ 1 , t ≥ 0 r s α 1 , α 2 ( t ) = Pr ( S α 1 , α 2 ( t ) = s | Q α 1 , α 2 ( 0 ) = ( 0 , 0 ) ) , 1 ≤ s ≤ k , t ≥ 0 . Then, q n α 1 , α 2 ( t ) = ∑ s = 1 k p n , s α 1 , α 2 ( t ) and r s α 1 , α 2 ( t ) = ∑ n = 1 ∞ p n , s α 1 , α 2 ( t ) , t ≥ 0 . For 1 ≤ s ≤ k, we add the differential equations given in Theorem 1 to obtain ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) p 0 α 1 , α 2 ( t ) = − λ p 0 α 1 , α 2 ( t ) + k μ p 1 , 1 α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) q 1 α 1 , α 2 ( t ) = − λ q 1 α 1 , α 2 ( t ) + k μ ( p 2 , 1 α 1 , α 2 ( t ) − p 1 , 1 α 1 , α 2 ( t ) ) + λ p 0 α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) q n α 1 , α 2 ( t ) = − λ q n α 1 , α 2 ( t ) + k μ ( p n + 1 , 1 α 1 , α 2 ( t ) − p n , 1 α 1 , α 2 ( t ) ) + λ q n − 1 α 1 , α 2 ( t ) , n ≥ 2 , with p 0 α 1 , α 2 ( 0 ) = 1 and q n α 1 , α 2 ( 0 ) = 0, n ≥ 1.

Similarly, for the classical Erlang queue, we have (see [4], p. 3257) d d t p 0 ( t ) = − λ p 0 ( t ) + k μ p 1 , 1 ( t ) , d d t q 1 ( t ) = − λ q 1 ( t ) + k μ ( p 2 , 1 ( t ) − p 1 , 1 ( t ) ) + λ p 0 ( t ) , d d t q n ( t ) = − λ q n ( t ) + k μ ( p n + 1 , 1 ( t ) − p n , 1 ( t ) ) + λ q n − 1 ( t ) , n ≥ 2 , with p 0 ( 0 ) = 1 and q n ( 0 ) = 0, n ≥ 1. Here, q n ( t ) = Pr ( N ( t ) = n | Q ( 0 ) = ( 0 , 0 ) ) , n ≥ 1 , t ≥ 0 r s ( t ) = Pr ( S ( t ) = s | Q ( 0 ) = ( 0 , 0 ) ) , 1 ≤ s ≤ k , t ≥ 0 . Note that q n ( t ) = ∑ s = 1 k p n , s ( t ) and r s ( t ) = ∑ n = 1 ∞ p n , s ( t ).

Recall that the first arrival time is the random time at which the first arrival occurs, and the inter-arrival time is defined as the time between two successive arrivals.

In the following results, we obtain the distributions of inter-arrival, inter-phase and sojourn times of mixed time-changed Erlang queue. These results are useful in understanding the customers arrival, phase service and overall system pattern.

The inter-arrival times T α 1 , α 2 of the mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 are independent and identically distributed (iid) with the following distribution: Pr ( T α 1 , α 2 > t ) = ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , ( α 1 − α 2 ) r + 1 r + 1 ( − λ c 1 t α 1 ) − ∑ r = 0 ∞ ( − c 2 c 1 ) r + 1 t ( α 1 − α 2 ) ( r + 1 ) E α 1 , ( α 1 − α 2 ) ( r + 1 ) + 1 r + 1 ( − λ c 1 t α 1 ) , t ≥ 0 where E α , β γ ( · ) is the three parameter Mittag-Leffler function as defined in (4).

Note that the inter-arrival times of mixed time-changed Erlang queue are independent. To obtain the distribution of inter-arrival times, let us define a process { N ∗ ( t ) } t ≥ 0 whose state probabilities q n ( t ) = Pr ( N ∗ ( t ) = n ), n ≥ 0 solve the following Cauchy problem (see [4], p. 3258): d d t q m ( t ) = − λ q m ( t ) , d d t q m + 1 ( t ) = λ q m ( t ) , d d t q n ( t ) = 0 , n ≥ 0 , n ≠ m , m + 1 , with initial conditions q m ( 0 ) = 1 and q n ( 0 ) = 0 for n ≥ 0, n ≠ m. Here, { N ∗ ( t ) } t ≥ 0 is a counting process with the same arrival rate λ as that of { N ( t ) } t ≥ 0 . It starts at m and m + 1 is its absorbing state. Observe that q m ( t ) + q m + 1 ( t ) = 1, t ≥ 0.

Let us denote T as the arrival time of first new customer in the classical Erlang queue { Q ( t ) } t ≥ 0 starting from N ∗ ( 0 ) = m. Then, its distribution function F T ( · ) is given by F T ( t ) : = Pr ( T ≤ t ) = 1 − Pr ( T > t ) = 1 − Pr ( N ∗ ( t ) = m ) = Pr ( N ∗ ( t ) = m + 1 ) = q m + 1 ( t ) , t ≥ 0 .

Let N ∗ α 1 , α 2 ( t ) = N ∗ ( Y α 1 , α 2 ( t ) ), t ≥ 0 be the time-changed process whose state probabilities are denoted by q m α 1 , α 2 ( t ) = Pr ( N ∗ α 1 , α 2 ( t ) = m ) and q m + 1 α 1 , α 2 ( t ) = Pr ( N ∗ α 1 , α 2 ( t ) = m + 1 ). Here, { N ∗ ( t ) } t ≥ 0 and { Y α 1 , α 2 ( t ) } t ≥ 0 are independent. So, { N ∗ α 1 , α 2 ( t ) } t ≥ 0 is a process which starts at m with m + 1 being its absorbing state. Let T α 1 , α 2 be the arrival time of the first new customer in the mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 starting from N ∗ α 1 , α 2 ( 0 ) = m. Thus, F T α 1 , α 2 ( t ) = Pr ( T α 1 , α 2 ≤ t ) = q m + 1 α 1 , α 2 ( t ).

Let { K ( t ) } t ≥ 0 be the counting process which counts the number of arrivals in the Erlang queue { Q ( t ) } t ≥ 0 . This counting process is obtained by assuming service times to be degenerate to ∞, that is, by setting μ = 0 in the Erlang queue. As { K ( t ) } t ≥ 0 is a Markov process, its time-changed variant { K α 1 , α 2 ( t ) } t ≥ 0 , where K α 1 , α 2 ( t ) = K ( Y α 1 , α 2 ( t ) ), t ≥ 0 is a semi-Markov process. Let T n be its nth jump time. Then, discrete time process K α 1 , α 2 ( T n ) is a time-homogeneous Markov chain (see [13]). Thus, we can say that the inter-arrival times of mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 are distributed according to T α 1 , α 2 .

Now, we show that the state probabilities q m α 1 , α 2 ( t ) and q m + 1 α 1 , α 2 ( t ) of { N ∗ α 1 , α 2 ( t ) } t ≥ 0 solve the following differential equations: ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) q m α 1 , α 2 ( t ) = − λ q m α 1 , α 2 ( t ) , ( c 1 d α 1 d t α 1 + c 2 d α 2 d t α 2 ) q m + 1 α 1 , α 2 ( t ) = λ q m α 1 , α 2 ( t ) , with q m α 1 , α 2 ( 0 ) = 1 and q m + 1 α 1 , α 2 ( 0 ) = 0.

The Eq. (51) holds true if and only if the Laplace transform q ˜ m α 1 , α 2 ( z ) of q m α 1 , α 2 ( t ) solves: ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ( z q ˜ m α 1 , α 2 ( z ) − 1 ) = − λ q ˜ m α 1 , α 2 ( z ) , z > 0 , where we have used (8).

On taking the Laplace transform on both sides of q m α 1 , α 2 ( t ) = Pr ( N ∗ ( Y α 1 , α 2 ( t ) ) = m ) = ∫ 0 ∞ q m ( y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) , t ≥ 0 and by using (10), we get q ˜ m α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ q m ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y , z > 0 . On substituting (54) in (53), we obtain ( c 1 z α 1 + c 2 z α 2 ) ∫ 0 ∞ q m ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y − 1 = − λ ∫ 0 ∞ q m ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y . By using (50) in (55) and on solving it, we get an identity. This establishes that (51) holds true. Similarly, it can be shown that (52) holds true.

On using Eq. (2.41) of [5], the Eq. (51) along with the initial condition q m α 1 , α 2 ( 0 ) = 1 can be solved to get the required result. This completes the proof.  □

From (53), it follows that the Laplace transform of q m α 1 , α 2 ( t ) = Pr ( T α 1 , α 2 > t ) is q ˜ m α 1 , α 2 ( z ) = c 1 z α 1 − 1 + c 2 z α 2 − 1 λ + c 1 z α 1 + c 2 z α 2 , z > 0 .

The probability density function of inter-arrival times of mixed time-changed Erlang queue is given by (see [5], Eq. (2.38)) f T α 1 , α 2 ( t ) = λ c 1 t α 1 − 1 ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , α 1 + ( α 1 − α 2 ) r r + 1 ( − λ t α 1 c 1 ) , t ≥ 0 , which has the following asymptotic relation (see [5], Eq. (2.47)): f T α 1 , α 2 ( t ) ∼ λ t α 1 − 1 c 1 Γ ( α 1 ) as t → 0. Here, f ( t ) ∼ g ( t ) as t → a means f and g are asymptotically equal as t → a, that is, lim t → a f ( t ) g ( t ) = 1. So, the asymptotic behavior of inter-arrival times depends on fractional parameter α 1 as t → 0. Also, f T α 1 , α 2 ( t ) ∼ c 2 α 2 t − α 2 − 1 λ Γ ( 1 − α 2 ) as t → ∞ (see [5], (2.48)). Hence, the asymptotic behavior of inter-arrival times depends on the fractional parameter α 2 as t → ∞.

The proof of the following result follows similar lines to that of Theorem 6. Thus, it is omitted.

The inter-phase times P α 1 , α 2 of mixed time-changed Erlang queue are iid with the following distribution: Pr ( P α 1 , α 2 > t ) = ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , ( α 1 − α 2 ) r + 1 r + 1 ( − k μ c 1 t α 1 ) − ∑ r = 0 ∞ ( − c 2 c 1 ) r + 1 t ( α 1 − α 2 ) ( r + 1 ) E α 1 , ( α 1 − α 2 ) ( r + 1 ) + 1 r + 1 ( − k μ c 1 t α 1 ) , t ≥ 0 where E α , β γ ( · ) is the three parameter Mittag-Leffler function as defined in (4).

The probability density function of inter-phase times of mixed time-changed Erlang queue is given by (see [5], Eq. (2.38)) f P α 1 , α 2 ( t ) = k μ c 1 t α 1 − 1 ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , α 1 + ( α 1 − α 2 ) r r + 1 ( − k μ t α 1 c 1 ) , t ≥ 0 , which has the following asymptotic relation (see [5], Eq. (2.47)): f P α 1 , α 2 ( t ) ∼ k μ t α 1 − 1 c 1 Γ ( α 1 ) as t → 0. So, the asymptotic behavior of inter-phase times depends on fractional parameter α 1 as t → 0. Also, f P α 1 , α 2 ( t ) ∼ c 2 α 2 t − α 2 − 1 k μ Γ ( 1 − α 2 ) as t → ∞ (see [5], (2.48)). Hence, the asymptotic behavior of inter-phase times depends on the fractional parameter α 2 as t → ∞.

Let Y 1 , Y 2 , … , Y k be the inter-phase times of mixed time changed Erlang queue and f Y 1 ( · ) be the probability density function of Y 1 , that is, f Y 1 ( t ) = d d t Pr ( Y 1 ≤ t ). It is given by f Y 1 ( t ) = f P α 1 , α 2 ( t ) t ≥ 0 . Let X = Y 1 + Y 2 + ⋯ + Y k . As Y i ’s are iid, we have f X ( t ) = d d t Pr ( X ≤ t ) = f Y 1 ∗ ( k ) ( t ) , t ≥ 0 . So, its Laplace transform is (see [5], Eq. (2.37)) f ˜ X ( z ) = ( f ˜ Y 1 ( z ) ) k = ( k μ k μ + c 1 z α 1 + c 2 z α 2 ) k , z > 0 . Hence, the service times of mixed time-changed Erlang queue are iid with the probability density function as given in (57).

For c 1 = 1 or c 2 = 1 in (58), it reduces to the Laplace transform of service times of fractional Erlang queue (see [4], p. 3260).

Recall that the sojourn time of a process is the time spent in a particular state before its transition to another state.

The proof of the following result follows along the similar lines to that of Theorem 6. Thus, it is omitted. Here, we need to work with the state probabilities of queue length process { L α 1 , α 2 ( t ) } t ≥ 0 .

The sojourn times S α 1 , α 2 of queue length process { L α 1 , α 2 ( t ) } t ≥ 0 in a state n > 0 are iid with the following distribution function: Pr ( S α 1 , α 2 > t ) = ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , ( α 1 − α 2 ) r + 1 r + 1 ( − ( λ + k μ ) c 1 t α 1 ) − ∑ r = 0 ∞ ( − c 2 c 1 ) r + 1 t ( α 1 − α 2 ) ( r + 1 ) E α 1 , ( α 1 − α 2 ) ( r + 1 ) + 1 r + 1 ( − ( λ + k μ ) c 1 t α 1 ) , t ≥ 0 where E α , β γ ( · ) is the three parameter Mittag-Leffler function defined in (4).

The probability density function of sojourn times times for the queue length process of mixed time-changed Erlang queue is given by (see [5], Eq. (2.38)) f S α 1 , α 2 ( t ) = ( λ + k μ ) c 1 t α 1 − 1 ∑ r = 0 ∞ ( − c 2 c 1 ) r t ( α 1 − α 2 ) r E α 1 , α 1 + ( α 1 − α 2 ) r r + 1 ( − ( λ + k μ ) t α 1 c 1 ) , t ≥ 0 , which has the following asymptotic relation (see [5], Eq. (2.47)): f S α 1 , α 2 ( t ) ∼ ( λ + k μ ) t α 1 − 1 c 1 Γ ( α 1 ) as t → 0. So, the asymptotic behavior of sojourn times depends on fractional parameter α 1 as t → 0. Also, f S α 1 , α 2 ( t ) ∼ c 2 α 2 t − α 2 − 1 ( λ + k μ ) Γ ( 1 − α 2 ) as t → ∞ (see [5], (2.48)). Hence, the asymptotic behavior of sojourn times depends on the fractional parameter α 2 as t → ∞.

Recall that the mixed time-changed Erlang queue is a generalization of fractional Erlang queue, for c 1 = 1 or c 2 = 1 it reduces to the fractional Erlang queue. As the inter-arrival and inter-phase times of fractional Erlang queue are not independent (see [4], Corollary 4.6), the inter-arrival and inter-phase times of mixed time-changed Erlang queue are not necessarily independent.

 Distribution of busy period

The busy period of a queue is defined as a duration of time since the first customer enters the empty system till the system becomes empty again.

In the next result, we obtain the distribution of a busy period. It is useful in understanding for how much time system stays busy when customers arrive in an empty system.

Let B α 1 , α 2 be the busy period of mixed time-changed Erlang queue, and let F B α 1 , α 2 ( t ) = Pr ( B α 1 , α 2 ≤ t ) be its distribution function. Then F B α 1 , α 2 ( t ) = ∑ h = 0 ∞ k μ A k , h 0 c 1 a k , h 0 h a k , h 0 ∗ ( a k , h 0 ) ( t ) , t ≥ 0 where a k , h 0 , A k , h 0 and h a k , h 0 ( · ) are as given in (21), (22) and (45), respectively.

Let us consider a process { L ∗ ( t ) } t ≥ 0 whose state probabilities p n ( t ) = Pr ( L ∗ ( t ) = n | L ∗ ( 0 ) = k ), n ≥ 0 solve the following system of differential equations (see [22], Eq. (36)-Eq. (38)): d d t p 0 ( t ) = k μ p 1 ( t ) , d d t p n ( t ) = − ( λ + k μ ) p n ( t ) + k μ p n + 1 ( t ) , 1 ≤ n ≤ k , d d t p n ( t ) = − ( λ + k μ ) p n ( t ) + k μ p n + 1 ( t ) + λ p n − k ( t ) , n ≥ k + 1 , with the initial condition p k ( 0 ) = 1.

Note that the process { L ∗ ( t ) } t ≥ 0 behaves similar to { L ( t ) } t ≥ 0 except that it starts from k instead of 0, which indicates that the first customer entered the system, and 0 is its absorbing state. Thus, by construction, the distribution of busy period of classical Erlang queue is F B ( t ) = p 0 ( t ), where F B ( t ) = Pr ( B ≤ t ) is given in [22].

Let { L ∗ ( t ) } t ≥ 0 and { Y α 1 , α 2 ( t ) } t ≥ 0 be independent. Also, let us define L ∗ α 1 , α 2 ( t ) : = L ∗ ( Y α 1 , α 2 ( t ) ), t ≥ 0 and denote the nth jump time of { L α 1 , α 2 ( t ) } t ≥ 0 by T n . Then, L α 1 , α 2 ( T n ) is time-homogeneous Markov chain. So, F B α 1 , α 2 ( t ) = p 0 α 1 , α 2 ( t ). Thus, F B α 1 , α 2 ( t ) = ∫ 0 ∞ F B ( y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) . On taking the Laplace transform of (60), we get F ˜ B α 1 , α 2 ( z ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∫ 0 ∞ F B ( y ) e − y ( c 1 z α 1 + c 2 z α 2 ) d y = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∑ h = 0 ∞ k λ h ( k μ ) k ( h + 1 ) h ! ( h k + k ) ! ∫ 0 ∞ ∫ 0 y u k + h ( k + 1 ) − 1 e − ( λ + k μ ) u · e − y ( c 1 z α 1 + c 2 z α 2 ) d u d y ( see [4], Eq. (13) ) = ( c 1 z α 1 − 1 + c 2 z α 2 − 1 ) ∑ h = 0 ∞ k λ h ( k μ ) k ( h + 1 ) h ! ( h k + k ) ! ∫ 0 ∞ u k + h ( k + 1 ) − 1 e − ( λ + k μ ) u · ∫ u ∞ e − y ( c 1 z α 1 + c 2 z α 2 ) d y d u = 1 z ∑ h = 0 ∞ k λ h ( k μ ) k ( h + 1 ) h ! ( h k + k ) ! ∫ 0 ∞ u k + h ( k + 1 ) − 1 e − ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) u d u = ∑ h = 0 ∞ k λ h ( k μ ) k ( h + 1 ) ( k + h ( k + 1 ) − 1 ) ! h ! ( h k + k ) ! z − 1 ( λ + k μ + c 1 z α 1 + c 2 z α 2 ) k + h ( k + 1 ) . On taking the inverse Laplace transform of (61), we obtain F B α 1 , α 2 ( t ) = ∑ h = 0 ∞ k μ A k , h 0 c 1 a k , h 0 L − 1 ( ( z − 1 a k , h 0 z α 1 + c 2 c 1 z α 2 + λ + k μ c 1 ) a k , h 0 ; t ) . Finally, on using (6) in (62), we get the required result.  □

 Conditional waiting time distribution

Methods used for deriving the distribution of waiting time in classical single-server queue depends on its Markovian structure (see [12]). However, these methods are not applicable in the case of mixed time-changed Erlang queue as it is a semi-Markov process. But a partial information about the distribution of its waiting time can be discussed by using some sort of conditioning arguments, as done in [4].

The following distribution function will be used: ψ α 1 , α 2 ( t ) = Pr ( P α 1 , α 2 > t ). It is given in (56). Let Y be a random variable with the following distribution function: F Y ( t ) = Pr ( P α 1 , α 2 ≤ t 0 + t | P α 1 , α 2 ≥ t 0 ) = 1 − ψ α 1 , α 2 ( t 0 + t ) ψ α 1 , α 2 ( t 0 ) , t ≥ 0 with parameters t 0 > 0 and k μ.

On taking the Laplace transform of f ψ ( t ) = ( t 0 + t ) r , r > 0, we obtain f ˜ ψ ( z ) = ∫ 0 ∞ e − z t ( t 0 + t ) r d t = e z t 0 ∫ t 0 ∞ e − z y y r d y = e z t 0 z r + 1 Γ ( r + 1 , z t 0 ) , z > 0 , where Γ ( y , z ) is the analytic extension of the upper incomplete gamma function on C 2 (see [15]). Also, let f Y ( t ) be the probability density function of Y. Then, f Y ( t ) = − 1 ψ α 1 , α 2 ( t 0 ) d d t ψ α 1 , α 2 ( t 0 + t ) . On taking the Laplace transform of (65) and by using (64), we get f ˜ Y ( z ) = 1 − ( ∑ r = 0 ∞ ∑ m = 0 ∞ ( − c 2 c 1 ) r ( − k μ c 1 ) m ( m + r ) ! Γ ( ( α 1 − α 2 ) r + α 1 m + 1 , z t 0 ) r ! m ! Γ ( α 1 m + ( α 1 − α 2 ) r + 1 ) · z ( α 2 − α 1 ) r − α 1 m − ∑ r = 0 ∞ ∑ m = 0 ∞ ( − c 2 c 1 ) r + 1 ( − k μ c 1 ) m ( m + r ) ! r ! m ! · Γ ( ( α 1 − α 2 ) ( r + 1 ) + α 1 m + 1 , z t 0 ) Γ ( α 1 m + ( α 1 − α 2 ) ( r + 1 ) + 1 ) z ( α 2 − α 1 ) ( r + 1 ) − α 1 m ) e z t 0 ψ α 1 , α 2 ( t 0 ) , z > 0 , where we have used (4) and (56).

Let W t α 1 , α 2 be the total time spend by the customer who enters the system at time t > 0. Consider the following random set: E t α 1 , α 2 ( ω ) = { 0 ≤ s ≤ t : L α 1 , α 2 ( s − ) ( ω ) = L α 1 , α 2 ( s ) ( ω ) + 1 } , which denotes the set of time instants before t at which a phase completion occurs. Now, we define the following random variable: N t α 1 , α 2 = | E t α 1 , α 2 | , t ≥ 0 which gives the count of phases completed before time t. Also, consider a random variable τ t α 1 , α 2 = sup E t α 1 , α 2 , N t α 1 , α 2 > 0 , t , N t α 1 , α 2 = 0 , namely, the last time instant before t at which a phase is completed. For a fixed t, the random time τ t α 1 , α 2 is a Markov time which takes values among the jump times of { S α 1 , α 2 ( t ) } t ≥ 0 , except for the case when τ t α 1 , α 2 = t. Here, { S α 1 , α 2 ( t ) } t ≥ 0 is a semi-Markov process and its semi-regenerative set coincides with its discontinuity points (see [8]). Moreover, if τ t α 1 , α 2 = t 0 for some t 0 < t, then the past and future of { S α 1 , α 2 ( t ) } t ≥ 0 are independent. Therefore, for the waiting time, we condition on the events { L α 1 , α 2 ( t ) = n } and { τ t α 1 , α 2 = t 0 }.

Let w α 1 , α 2 ( x ; t , t 0 , n ) d x = Pr ( W t α 1 , α 2 ∈ d x | τ t α 1 , α 2 = t 0 , L α 1 , α 2 ( t ) = n ) . A customer has to wait for the completion of n independent phases (as we have conditioned on { L α 1 , α 2 ( t ) = n }), one of the n phases started at t 0 that has not been completed till time t. Let T j , 1 ≤ j ≤ n be independent random variables, such that for 1 ≤ j ≤ n − 1, T j denotes the duration of jth phase whose service is not yet started and T n be the remaining time duration of the phase which is in service at time t. Thus, we have E ( W t α 1 , α 2 | τ t α 1 , α 2 = t 0 , L α 1 , α 2 ( t ) = n ) = ∑ j = 1 n T j . Note that T j , 1 ≤ j ≤ n − 1 are iid from Theorem 7. Thus, W = ∑ j = 1 n − 1 T j has the distribution given by (57). Also, T n is distributed as (63) with parameters t − t 0 and k μ. This is because we have conditioned on { τ t α 1 , α 2 = t 0 } which indicates that the last phase was completed at t 0 . Then, we obtain w α 1 , α 2 ( x ; t , t 0 , n ) = f W ∗ f T n ( x ) . On taking the Laplace transform of (67), and using (58) and (66), we get w ˜ ( z ; t , t 0 , n ) = ( k μ k μ + c 1 z α 1 + c 2 z α 2 ) n − 1 ( 1 − e z ( t − t 0 ) ψ α 1 , α 2 ( t − t 0 ) ( ∑ r = 0 ∞ ∑ m = 0 ∞ ( − c 2 c 1 ) r ( − k μ c 1 ) m · ( m + r ) ! Γ ( ( α 1 − α 2 ) r + α 1 m + 1 , z ( t − t 0 ) ) r ! m ! Γ ( α 1 m + ( α 1 − α 2 ) r + 1 ) z ( α 2 − α 1 ) r − α 1 m − ∑ r = 0 ∞ ∑ m = 0 ∞ ( − c 2 c 1 ) r + 1 ( − k μ c 1 ) m ( m + r ) ! r ! m ! · Γ ( ( α 1 − α 2 ) ( r + 1 ) + α 1 m + 1 , z ( t − t 0 ) ) Γ ( α 1 m + ( α 1 − α 2 ) ( r + 1 ) + 1 ) z ( α 2 − α 1 ) ( r + 1 ) − α 1 m ) ) , z > 0 .

Sample paths simulation

Here, we give plots of sample path simulations of the mixed time-changed Erlang queue. The following results give the distributional properties of exponential random variable and mixed stable subordinator:

Let X be an exponential random variable with rate λ and { D α 1 , α 2 ( t ) } t ≥ 0 be mixed stable subordinator which is independent of X. Then, D α 1 , α 2 ( X ) = d T α 1 , α 2 , where T α 1 , α 2 is the inter-arrival time of the mixed time-changed Erlang queue.

Consider the following: Pr ( D α 1 , α 2 ( X ) > t ) = Pr ( X > Y α 1 , α 2 ( t ) ) = ∫ 0 ∞ Pr ( X > y ) Pr ( Y α 1 , α 2 ( t ) ∈ d y ) = ∫ 0 ∞ e − λ y Pr ( Y α 1 , α 2 ( t ) ∈ d y ) . On taking the Laplace transform of (69) and by using (10), we obtain L ( Pr ( D α 1 , α 2 ( X ) > t ) ) ( z ) = c 1 z α 1 − 1 + c 2 z α 2 − 1 λ + c 1 z α 1 + c 2 z α 2 , z > 0 whose inversion yields the required result.  □

Let X be a non-negative random variable independent of the mixed stable subordinator { D α 1 , α 2 ( t ) } t ≥ 0 . Then, D α 1 , α 2 ( X ) = d c 1 1 / α 1 X 1 / α 1 D α 1 ( 1 ) + c 2 1 / α 2 X 1 / α 2 D α 2 ( 1 ) , where D α 1 ( 1 ) and D α 2 ( 1 ) are independent random variables.

Note that, Pr ( D α 1 , α 2 ( X ) > t ) = ∫ 0 ∞ Pr ( D α 1 , α 2 ( y ) > t ) Pr ( X ∈ d y ) = ∫ 0 ∞ Pr ( c 1 1 / α 1 D α 1 ( y ) + c 2 1 / α 2 D α 2 ( y ) > t ) Pr ( X ∈ d y ) = ∫ 0 ∞ Pr ( c 1 1 / α 1 y 1 / α 1 D α 1 ( 1 ) + c 2 1 / α 2 y 1 / α 2 D α 2 ( 1 ) > t ) Pr ( X ∈ d y ) = Pr ( c 1 1 / α 1 X 1 / α 1 D α 1 ( 1 ) + c 2 1 / α 2 X 1 / α 2 D α 2 ( 1 ) > t ) , t ≥ 0 , where we have used D α 1 , α 2 ( t ) = d c 1 1 / α 1 D α 1 ( t ) + c 2 1 / α 2 D α 2 ( t ), t ≥ 0 (see [1], p. 704) such that { D α 1 ( t ) } t ≥ 0 and { D α 2 ( t ) } t ≥ 0 are independent stable subordinators. This proves the required result.  □

Plot (a) represents the sample path simulation of N α 1 , α 2 ( t ) and Plot (b) represents the sample path simulation of N α 1 , α 2 ( T n ) such that T n are the jump times of { N α 1 , α 2 ( t ) } t ≥ 0 for parameters c 1 = 0.4, c 2 = 0.6, α 1 = 0.5, α 2 = 0.3, k = 4, λ = 6 and μ = 5

Two step function plots: (a) vs t shows peaks near 0 and 1.0; (b) vs n shows peaks near n=4 and rising trend after n=12.

By using Lemma 1 and Lemma 2, we have the following result:

Let X be an exponential random variable with rate λ > 0, { D α 1 , α 2 ( t ) } t ≥ 0 be mixed stable subordinator independent of X, and Y be a random variable such that Y = d T α 1 , α 2 . Then, Y = d c 1 1 / α 1 X 1 / α 1 D α 1 ( 1 ) + c 2 1 / α 2 X 1 / α 2 D α 2 ( 1 ) , where D α 1 ( 1 ) and D α 2 ( 1 ) are independent.

Plot (a) represents the sample path simulation of N ( t ) with parameters k = 4, λ = 6 and μ = 5, and Plot (b) represents the sample path simulation of N α 1 , α 2 ( t ) with parameters c 1 = 0.4, c 2 = 0.6, α 1 = 0.5, α 2 = 0.3, k = 4, λ = 6 and μ = 5

Two step function plots vs t (0–1.8). (a) values rise to ~7 then decline; (b) values rise to 4, stabilize at 3 after t=0.8.

Plot (a) represents the sample path simulation of N α 1 ( t ) with parameters α 1 = 0.5, k = 4, λ = 6 and μ = 5, and Plot (b) represents the sample path simulation of N α 1 , α 2 ( t ) with parameters c 1 = 0.4, c 2 = 0.6, α 1 = 0.5, α 2 = 0.3, k = 4, λ = 6 and μ = 5

Two step function plots vs t. (a) rises from 0 to ~7 with jumps near t=0.5 and t=2.6. (b) oscillates between 0 and 2 with jumps near t=0.1 and t=2.6.

We use the modified version of Gillespie algorithm (see [7]) for simulating the mixed time-changed Erlang queue. It is based on the fact that Q α 1 , α 2 ( T n ) is again a Markov chain with the same transition probabilities as that of Erlang queue { Q ( t ) } t ≥ 0 , where T n is the nth jump time of { Q α 1 , α 2 ( t ) } t ≥ 0 . This algorithm is used in [4] for simulating the fractional Erlang queue. We use a similar algorithm to simulate the mixed time-changed Erlang queue in which a different distribution is utilized for the random variable I that appears in the algorithm used in [4].

The simulation algorithm is described as follows.

Step 1. Initialize the process by setting Q α 1 , α 2 ( 0 ) = ( 0 , 0 ) , T 0 = 0 . Step 2. Let T n be the nth jump time of the mixed time-changed Erlang queue { Q α 1 , α 2 ( t ) } t ≥ 0 . Assume that the queue has been simulated till the nth event. If Q ( T n ) = 0, generate a random variable I = d T α 1 , α 2 ; otherwise, generate I = d S α 1 , α 2 , where the distributions of T α 1 , α 2 and S α 1 , α 2 are given in (49) and (59), respectively.

Step 3. Define the next jump time as T n + 1 = T n + I .

Step 4. Whenever Q α 1 , α 2 ( T n ) = ( 0 , 0 ), set Q α 1 , α 2 ( T n + 1 ) = ( 1 , k ) . If Q α 1 , α 2 ( T n ) ∈ H, generate a random variable U such that it is uniformly distributed on ( 0 , 1 ). For U < λ λ + k μ , set N ( T n + 1 ) = N ( T n ) + 1 , S ( T n + 1 ) = S ( T n ) . For U ≥ λ λ + k μ , phase is served and the state is updated according to the following cases:

If S α 1 , α 2 ( T n ) > 1, set N α 1 , α 2 ( T n + 1 ) = N α 1 , α 2 ( T n ) , S α 1 , α 2 ( T n + 1 ) = S α 1 , α 2 ( T n ) − 1 .

If S α 1 , α 2 ( T n ) = 1 and N α 1 , α 2 ( T n ) > 1, set N α 1 , α 2 ( T n + 1 ) = N α 1 , α 2 ( T n ) − 1 , S α 1 , α 2 ( T n + 1 ) = k .

If S α 1 , α 2 ( T n ) = 1 and N α 1 , α 2 ( T n ) = 1, set N α 1 , α 2 ( T n + 1 ) = S α 1 , α 2 ( T n + 1 ) = 0 .

Step 5. The sample path of the whole queue can be obtained by setting Q α 1 , α 2 ( t ) = Q α 1 , α 2 ( T n − 1 ) , t ∈ [ T n − 1 , T n ) .

On substituting c 1 = 1 or c 2 = 1 in (20), (42), Theorem 1, Lemma 2 and Proposition 3, the results for mixed time-changed Erlang queue reduce to the corresponding results for fractional Erlang queue. On substituting c 1 = 1 in (46), (68) and Theorem 3, Theorem 6 - Theorem 9, the results for mixed time-changed Erlang queue reduce to the corresponding results for fractional Erlang queue (see [4]).

 Acknowledgments

The authors would like to thank the anonymous reviewers for carefully reading the manuscript. Reviewer’s comments were useful in improving the quality of paper.

 References Aletti, G., Leonenko, N., Merzbach, E.: Fractional Poisson fields. J. Stat. Phys. 170(4), 700–730 (2018). MR3764004. https://doi.org/10.1007/s10955-018-1951-y Applebaum, D.: Lévy Processes and Stochastic Calculus. Cambridge University Press, Cambridge (2009). MR2512800. https://doi.org/10.1017/CBO9780511809781 Ascione, G., Leonenko, N., Pirozzi, E.: Fractional queues with catastrophes and their transient behaviour. Math. 6(9), 159 (2018). MR4092404. https://doi.org/10.3390/math6090159 Ascione, G., Leonenko, N., Pirozzi, E.: Fractional Erlang queues. Stoch. Process. Appl. 130(6), 3249–3276 (2020). MR4092404. https://doi.org/10.1016/j.spa.2019.09.012 Beghin, L.: Random-time processes governed by differential equations of fractional distributed order. Chaos Solitons Fractals 45(11), 1314–1327 (2012). MR2990245. https://doi.org/10.1016/j.chaos.2012.07.001 Beghin, L., Orsingher, E.: Fractional Poisson processes and related planar random motions. Electron. J. Probab. 14(61), 1790–1827 (2009). MR2535014. https://doi.org/10.1214/EJP.v14-675 Cahoy, D.O., Polito, F., Phoha, V.: Transient behavior of fractional queues and related processes. Methodol. Comput. Appl. Probab. 17(3), 739–759 (2015). MR3377858. https://doi.org/10.1007/s11009-013-9391-2 Cinlar, E.: Markov additive processes and semi-regeneration (1974). Discussion Paper No. 118. Northwestern University. Di Crescenzo, A., Giorno, V., Kumar, B.K., Nobile, A.G.: On the M / M / 1 queue with catastrophes and its continuous approximation. Queueing Syst. 43(4), 329–347 (2003). MR1976263. https://doi.org/10.1023/A:1023261830362 Dhillon, M., Kataria, K.K.: On multiparameter generalized counting process and its time-changed variants. J. Theor. Probab. 38(4), 81 (2025). MR4963663. https://doi.org/10.1007/s10959-025-01451-8 Garg, S., Pathak, A.K., Maheshwari, A.: Tempered space-time fractional negative binomial process. Methodol. Comput. Appl. Probab. 27, 51 (2025). MR4920025. https://doi.org/10.1007/s11009-025-10179-1 Gaver, D.P.: The influence of servicing times in queuing processes. J. Oper. Res. Soc. Am. 2(2), 139–149 (1954). https://doi.org/10.1287/opre.2.2.139 Gikhman, I.I., Skorokhod, A.V.: The Theory of Stochastic Processes. II. Grundlehren Math. Wiss., vol. 218 (1975). MR0375463 Giorno, V., Nobile, A.G., Pirozzi, E.: A state-dependent queueing system with asymptotic logarithmic distribution. J. Math. Anal. Appl. 458(2), 949–966 (2018). MR3724709. https://doi.org/10.1016/j.jmaa.2017.10.004 Gradshteyn, I.S., Ryzhik, I.M., Jeffrey, A., Zwillinger, D.: Table of Integrals, Series, and Products, 7th edn. Academic Press, Boston (2007). MR2360010 Griffiths, J.D., Leonenko, G.M., Williams, J.E.: The transient solution to M / E k / 1 queue. Oper. Res. Lett. 34, 349–354 (2006). MR2204349. https://doi.org/10.1016/j.orl.2005.05.010 Haubold, H.J., Mathai, A.M., Saxena, R.K.: Mittag-Leffler functions and their applications. J. Appl. Math. 298628, 51 pp. (2011). MR2800586. https://doi.org/10.1155/2011/298628 Kataria, K.K., Khandakar, M.: Mixed fractional risk process. J. Math. Anal. Appl. 504(1) (2021). MR4269613. https://doi.org/10.1016/j.jmaa.2021.125379 Kataria, K.K., Vishwakarma, P.: On time-changed linear birth-death-immigration process. J. Theoret. Probab. 38(1), Article 21 (2025). MR4838514. https://doi.org/10.1007/s10959-024-01387-5 Kilbas, A.A., Srivastava, H.M., Trujillo, J.J.: Theory and Applications of Fractional Differential Equations. Elsevier Science B. V., Amsterdam (2006). MR2218073. https://doi.org/10.1016/S0304-0208(06)80001-0 Luchak, G.: The solution of the single-channel queuing equations characterized by a time-dependent Poisson-distributed arrival rate and a general class of holding times. Oper. Res. 4(6), 711–732 (1956). MR0083415. https://doi.org/10.1287/opre.4.6.711 Luchak, G.: The continuous time solution of the equations of the single channel queue with a general class of service-time distributions by the method of generating functions. J. R. Stat. Soc. Ser. B, Methodol. 20(1), 176–181 (1958). MR0096312. https://doi.org/10.1111/j.2517-6161.1958.tb00287.x Meerschaert, M.M., Nane, E., Vellaisamy, P.: The fractional Poisson process and the inverse stable subordinator. Electron. J. Probab. 16(59), 1600–1620 (2011). MR2835248. https://doi.org/10.1214/EJP.v16-920 Meerschaert, M.M., Straka, P.: Inverse stable subordinators. Math. Model. Nat. Phenom. 8(2), 1–16 (2013). MR3049524. https://doi.org/10.1051/mmnp/20138201 Meoli, A.: Bivariate gamma subordination for a Poisson shock model. ALEA Lat. Am. J. Probab. Math. Stat. 22, 73–91 (2025). MR4860602. https://doi.org/10.30757/alea.v22-02 Orsingher, E., Polito, F.: On a fractional linear birth-death process. Bernoulli 17(1), 114–137 (2011). MR2797984. https://doi.org/10.3150/10-BEJ263 Orsingher, E., Polito, F.: The space-fractional Poisson process. Stat. Probab. Lett. 82(4), 852–858 (2012). MR2899530. https://doi.org/10.1016/j.spl.2011.12.018 Pote, R.B., Kataria, K.K.: On Erlang Queue with Multiple Arrivals and its Time-Changed Variant. Methodol. Comput. Appl. Probab. 27, 69 (2025). MR4951585. https://doi.org/10.1007/s11009-025-10200-7 Spath, D., Fähnrich, K.P.: Advances in Services Innovations. Springer, Berlin, Heidelberg (2006) https://doi.org/10.1007/978-3-540-29860-1 Tathe, K., Ghosh, S.: Non-homogeneous generalized fractional Skellam process. Fract. Calc. Appl. Anal. 28(6), 2687–2747 (2025). MR5009328. https://doi.org/10.1007/s13540-025-00450-0