Analysis of the queue lengths in a priority retrial queue with constant retrial policy

In this paper, we analyze a priority queueing system with a regular queue and an orbit. Customers in the regular queue have priority over the customers in the orbit. Only the first customer in the orbit (if any) retries to get access to the server, if the queue and server are empty (constant retrial policy). In contrast with existing literature, we assume different service time distributions for the high- and low-priority customers. Closed-form expressions are obtained for the probability generating functions of the number of customers in the queue and orbit, in steady-state. Another contribution is the extensive singularity analysis of these probability generating functions to obtain the stationary (asymptotic) probability mass functions of the queue and orbit lengths. Influences of the service times and the retrial policy are illustrated by means of some numerical examples.


1.
Introduction. In queueing systems, it is not uncommon that some customers who find the server busy upon arrival leave the system, but return after some random time to request for service again. This class of customers can be visualised as being located in a virtual waiting room, hereafter called the orbit [11]. In queueing theory, such queueing systems are called retrial queues. There exists a vast literature on the theory of retrial queues [3]. Perhaps the most well known books are those of [4,10], which can serve as a solid introduction to retrial queues. In the literature, two types of retrial policies are typically identified. The first one, called the classical retrial policy, assumes that all the customers in the orbit retry to get access to the server. For analytical feasibility reasons, these retrials are typically modelled as independent Poisson processes with a common arrival rate. An excellent survey of these models can be found in [8]. The second policy, referred to as the constant retrial policy, assumes that only one customer of the orbit can retry for service. Therefore, in the constant retrial policy, the orbit is commonly presented as a FCFS queue. In contrast to the classical retrial policy, the retrial times in the constant retrial policy are independent of the orbit length. This property may be important when modeling practical applications, as illustrated in [20,21] where the principle of call blending in call centres is modelled as queues with a constant retrial policy.
The constant retrial policy was first introduced in [12], where a telephone exchange model with only one line is examined. In [12], blocked customers (incoming calls) form a feedback queue (this is the orbit) and only the customer at the head of the feedback queue can retry for service after an exponentially distributed retrial time. The model presented in [12], but with a general retrial time distribution is studied in [7]. Finally, in [16] the constant retrial policy with general retrial times and general service times is extensively studied. The author has obtained the probability generating function of the orbit length and the Laplace-Stieltjes transform of the waiting time, in steady-state. The most important feature of a retrial queue with constant retrial policy and general retrial times is that, after each service completion, 1 a competition between an exponential law (the arrivals) and a general retrial time determines the next customer who receives service. Retrial queues with constant retrial policy and general retrial times are therefore related to many queueing systems where the server spends some random amount of time to find the next customer to serve after each service completion.
If retrial queues also possess a regular queue, priority occurs naturally. Customers in the regular queue have priority over the customers in the orbit. This means that the customer at the head of the orbit can only start retrying when both the regular queue and the server becomes empty. If another customer arrives during a retrial time, this customer is served and the retrial has to start over when the regular queue becomes empty again. The priority can in principal be nonpreemptive or preemptive, the non-preemptive again being the most natural one. In [15], a preemptive constant retrial policy is studied. The author has obtained the probability generating functions of the stationary queue and orbit length, and the Laplace-Stieltjes transforms of the waiting time of both types of customers. In [5], a non-preemptive constant retrial policy is studied. The probability generating functions of the stationary queue and orbit lengths are obtained. An interesting result from this paper, is that the distribution of the priority queue length is independent of the retrial time distribution. While in [15] it is assumed that the service time distributions of the customers in the priority queue and customers in the orbit are different, this is not the case in [5]. However, there might be situations where the service time distributions are not necessarily the same. It might even be the reason why a customer joins either the queue or the orbit. For instance, customers with a small (non-urgent) service time who face a queue upon arrival, are (in certain situations) more likely to join the orbit. In [9], we generalized the model of [5] in two ways. We assumed different service time distributions between both classes of customers. Further, in contrast to [5], we analyzed the waiting times of both classes of customers.
Our present paper is an extended version of [9]. In [9], we obtained explicit expressions for the probability generating functions of the stationary queue and orbit length. Theoretically, these contain all information of the queue and orbit lengths. For example, we can compute relatively easy the mean and the variance of the queue and orbit lengths from their probability generating functions. However, obtaining the distribution from the probability generating function is a much more delicate task. If we want to know the probability that the orbit length equals n customers, we have to compute the n-th derivative of the probability generating function and evaluate at zero. This becomes quickly problematic for large n. An elegent way to approximate these probabilities is to asymptotically invert the probability generating function, resulting into explicit expressions for the probability mass function [14], which are accurate for large n. Asymptotics for the M/G/1 retrial queue, for the classical retrial policy, is studied for light-tailed service time distributions in [17] and for a subclass of heavy-tailed service time distributions in [22]. The asymptotics of the marginal and joint queue lengths in the (non-preemptive) M/M/1 priority queue is studied in [18]. We extend our analysis of [9] to include asymptotics of the queue and orbit lengths, for several different types of service time distributions (e.g. light-tailed, fat-tailed, . . .). Since we assume different service time distributions between customer classes, we have to examine several possible combinations of types of distributions. By obtaining these asymptotics, we now demonstrate a complete characterisation of the queue and orbit lengths of this queueing model. Our methodology is systematic in the sense that we use the same singularity analysis of the probability generating functions to obtain the asymptotics for all combinations of the service time distributions. This is the novelty and difference in comparison with previous studies in which adhoc methods were developed for specific types of service time distributions [17,22].
The paper is organized as follows. In Section 2, we provide a detailed description of the queueing model under consideration. Next, in Section 3 we discuss the stability condition and analyze the steady-state distributions using the supplementary variable technique. In Section 4 we give an overview of singularity analysis, which is the method we apply in Section 5 and Section 6 to analyze the asymptotic distributions of the queue and orbit lengths respectively. Section 7 discusses the most important performance measures of the queueing system. Finally, we discuss some numerical examples in Section 8. 2. Mathematical model. We assume a queueing system with an infinite-sized queue, an infinite-sized orbit, and one server. Class-1 and class-2 customers arrive to the system according to two independent Poisson processes with rate λ 1 and λ 2 respectively. If a new customer arrives to the system and the server is free, he gets service immediately. If the server is busy upon arrival, a class-1 customer waits in the queue, while a class-2 customer joins the orbit. We assume that only the first customer in the orbit retries to get access to the server when the queue and server are idle. Retrial times are characterized by means of i.i.d. continuous random variables with common cumulative distribution (cdf) R(x) and common Laplace-Stieltjes transform (LST) R * (s). Service times of class-j customers, j = 1, 2, are specified by i.i.d. continuous random variables with cdf B j (x), LST B * j (s) and mean β j . We define the arrival load as ρ λ 1 β 1 + λ 2 β 2 . Further, we assume that service times, retrial times and inter-arrival times are all mutually independent.
In the rest of this section we fix notations for the rest of the paper. We use the notation P for the probability measure, E for the expectation operator and 1 for the indicator function.
Let N 1 (t) be the number of class-1 customers in the priority queue at time t, not counting the one in the server. Similarly, let N 2 (t) the number of class-2 customers in the orbit at time t. Further, we add four auxiliary processes. First, C(t) summarizes the state of the system, more precisely C(t) = 0 when the server is idle at time t, C(t) = 1 when the server is serving a class-1 customer at time t and C(t) = 2 when the server is serving a class-2 customer at time t. When C(t) = 0 and N 2 (t) > 0, ξ 0 (t) represents the elapsed retrial time at time t; when C(t) = 1, ξ 1 (t) denotes the elapsed class-1 service time at time t; when C(t) = 2, ξ 2 (t) denotes the elapsed class-2 service time at time t. The stochastic vector {C(t), N 1 (t), N 2 (t), ξ 0 (t), ξ 1 (t), ξ 2 (t)} is then a Markov Process. We assume that the stationary distribution of this stochastic vector exists. Therefore we define the following limiting probabilities: We also define the following probability generating functions (pgf) Next, we definer(x) as the probability density function (pdf) of the retrial time R, given that the retrial time is at least x. We defineb j (x) as the pdf of the service time of a class-j customer, given that the service time is at least x. We thus havē Let λ λ 1 +λ 2 the total arrival rate to the system. Let p 1 and p 2 the probability that a customer is of class-1 and of class-2 respectively. These probabilities are given by p j = λ j /λ, j = 1, 2.
For further use, we define A j (z 1 , z 2 ) as the joint pgf of the class-1 and class-2 arrivals during a class-j service time and A(z 1 , z 2 ) as the joint pgf of the class-1 and class-2 arrivals during a randomly chosen service time. They are calculated as In particular, we are interested in the steady-state distributions of the orbit and queue lengths. We define as the partial pgf of the number of customers in the orbit when the server is idle, and as the partial joint pgf of the number of customers in queue and orbit when the server is serving a type-j customer.
3. Queueing analysis. In this section, we use the supplementary variable technique to derive expressions for the steady-state (partials) pgfs of the the number of customers in the orbit when the server is idle and of the number of customers in queue and orbit when the server is serving a type-j customer, j = 1, 2. First, we indicate the stability condition of the queueing system.

Stability condition.
A stability condition could be derived by studying the embedded Markov Chain at departure epochs, similarly as in [5]. It can be shown that this yields the same condition as in [5], namely Due to the competition of two exponentials with rate λ j , j = 1, 2 and a general retrial time, R * (λ) is the probability that there are no new arrivals during a retrial time and is therefore also the probability that a single retrial is successful. We can also rewrite (9) as to which we will now give an intuitive explanation. From (6) it follows that p 2 ρ is the mean number of class-2 arrivals during a randomly chosen service time. In stability, the probability that the regular queue is empty after a departure instant is equal to 1 − p 1 ρ. Notice now that between two departure instants, a customer from the orbit (if any) can only have entered the server if at the first departure instant the regular queue was left empty and his retrial was successful. Since only the first customer from the orbit can try to get access to the server, R * (λ)(1 − p 1 ρ) is the expected number of customers from the orbit entering the server between two departure instants. Condition (10) can therefore be interpret as follows: the average number of customers entering the orbit during a randomly chosen service time must be strictly less then the average number of customers from the orbit entering the server between two departure instants.
3.2. Steady-state distributions. If the regular queue is non-empty, a customer from that queue is served next. Otherwise, a customer from the orbit (if any) retries to access the server. If his retrial time finishes before a new arrival occurs, he is served. Otherwise, the newly arriving customer starts service, irrespective of his type, and the customer at the head of the orbit stops retrying until the server and the priority queue are empty again. This leads to the following equilibrium equations The boundary conditions are p (2) i,n (0) = 1(i = 0, n = 0)λ 2 p 0 + 1(i = 0, n > 0)λ 2 ∞ 0 p (0) 0,n (x)dx Assuming stationarity, the equilibrium Equations (12), (13) and (14) are transformed to the following partial differential equations for the pgfs P 0 (x, z), P 1 (x, z 1 , z 2 ) and P 2 (x, z 1 , z 2 ): , whose solutions are given by In terms of the pgfs, the boundary conditions can be written as where we also used (11) in the first equation.
If we substitute (18), (19) and (20) in the previous equations and by using (5), we get We are left with five unknowns: p 0 , P 0 (0, z 2 ), P 1 (0, z 1 , z 2 ), P 1 (0, 0, z 2 ) and P 2 (0, 0, z 2 ). 2 First, we substitute (26) in (24) and solve for P 0 (0, z 2 ), Substituting (26) in (25) and solving to P 1 (0, z 1 , z 2 ), we can write It is well-known that f (z 1 ) := z 1 − A 1 (z 1 , z 2 ) has one zero inside the complex unit disk for each z 2 with |z 2 | < 1, see for instance [25]. Denote this zero by Y (z 2 ). It can be proven that this is an analytic function inside the complex unit disk and that Y (1) = 1, at least when the queueing system is stable. For more details, we refer to Appendix A in [23]. By substituting z 1 by Y (z 2 ) in (28), we find Substituting the above equation into (27) yields Substituting Equation (30) back into (27) and solving for P 1 (0, 0, z 2 ) yields Inserting (30) and (31) into (28) gives The forth unknown P 2 (0, z 1 , z 2 ) is given as Finally, the last unknown p 0 can be found from the normalisation condition Q(1) + U 1 (1, 1) + U 2 (1, 1) = 1. Using the Definitions (7), (8) and the Expressions (18), (19) and (20), one easily finds It now follows that Remark that the stability condition implies that p 0 > 0, i.e. the total system will occasionally be idle. Summarized, we have obtained partial pgfs of the number of customers in the queue and orbit when the server is either busy or idle. The complete expressions are given by: From these three important pgfs, the marginal pgfs of the numbers of customers in the queue and orbit are easily deduced. The pgf N 1 (z) of N 1 is given by The pgf of N 1 is thus independent of the distribution of the retrial times. As mentioned before, this was also observed in [5]. The pgf N 2 (z) of N 2 is given by 4. Singularity analysis of a probability generating function. Now we have obtained the stationary pgf of the queue and orbit length, we are interested in deriving asymptotic expressions for the stationary queue and orbit length distribution. This section gives a brief overview of the method we apply to obtain such asymptotics in Section 5 of the queue length distribution and in Section 6 of the orbit length distribution. Denote by X a generic discrete random variable with pgf X(z). By definition, the Taylor expansion of X(z) at the origin is equal to Let ζ be the radius of convergence of X(z). Because X(z) is a pgf, we have that ζ ≥ 1. According to Pringsheim's theorem [14, Th. IV.6], ζ is a singularity of X(z). Throughout the rest of this paper, we will assume that the encountered pgfs have a single dominant singularity 3 . If the dominant singularity of a pgf X(z) is a simple pole, it then relatively easily follows, from the Laurent series of X(z) at ζ, that P[X = n] has an asymptotic geometric behaviour. In priority queueing systems, other type of singularities can occur (c.f. [18,19]), even if the service time distributions are well-behaved functions. Therefore, we need a more general method. The method we apply, is a singularity analysis of the pgf X(z) of interest. For a detailed explanation of this method we refer the reader to [14,Chap. 6] or [24,Sect. 3]. We provide Corollary VI.1 in [14] as a theorem, which we will make use of frequently.
Assume the following two conditions on an analytic function at the origin f (z), with Taylor coefficients f n : 1. There exists numbers r, φ, r > 1 and 0 < φ < π 2 , such that the function f (z) is analytic in ∆(φ, r), where we write f n ∼ g n for n → ∞ if lim n→∞ g n /f n = 1 and Γ(·) stands for the Gamma function.
1. The first assumption is rather technical. It means that the function can be analytic continued in a domain, larger than its disk of convergence. This domain has a specific form, defined by (41), and is called a ∆-domain. 4 . A function that is analytic in a ∆-domain is then called a ∆-analytic function. 2. By the rescaling rule, for z → ζ in some ∆-domain, we have that Our strategy to derive tail asymptotics for a discrete random variable X consists of determining ζ, c and α, such that we can write Applying Theorem 4.1 then yields Notice that the analytic terms d j (ζ − z) j , j = 0, . . . , −α do not contribute to the asymptotic coefficients. Therefore, unless needed, we never explicitly compute the constants d i . Usually, we will write such unspecified constants as x i ,x i orx i . We will assume that B * 1 (s) and B * 2 (s) are ∆-analytic to ensure N 1 (z) and N 2 (z) are ∆-analytic. For notational purposes, we will not write explicitly that the limit is taken in a ∆-domain in the remainder.
The singularities of N 1 (z) and N 2 (z) depend on the pgfs A 1 (z 1 , z 2 ), A 2 (z 1 , z 2 ) and Y (z). These pfgs are defined through the LST of the service time distributions B * 1 (s) and B * 2 (s). Therefore, it is clear that the asymptotics strongly depend on the service time distributions. Let us define −τ j as the dominant singularity of B * j (s) for j = 1, 2. Recall that for an LST, the dominant singularity is the rightmost singularity on the negative real axis; hence, τ j ≥ 0. To proceed carefully, we classify service time distributions according to the behaviour of their LST B * j (s) in their rightmost singularity −τ j : In case the distribution has no singularities, we say it is of type-1 (τ j = ∞). This classification for probability distributions on the positive halfline was first introduced in [2]. Type-1 distributions are well-behaved functions and are also referred to as light-tailed distributions. It is worth noting that the dominant singularity of these distributions is not necessary a pole, for instance the Gamma distribution with LST (1 + θs) −k , with parameters θ and k real and strictly positive. Type-2 distributions are called semi-exponential distributions, because the tail is dominated by an exponential but is not a pure exponential one, like type-1 distributions. Lastly, type-3 distributions are the heavy-tailed distributions, because their tail is not dominated by an exponential one. Often encountered type-3 distributions in practice are the Pareto distribution and the lognormal distribution. For more discussion, we refer the reader to [1]. In our analysis, we will consider (almost) every possible distributions types for B * 1 (s) and B * 2 (s). In case of type-3 distributions, we restrict ourself to the important subclass of distributions of the form where −α j > 1. We refer to these distributions as fat-tailed distributions.
From this expression, we see that these are distributions where the first m moments exist, m = −α j , and that have an infinite (m + 1) st moment. The assumption −α j > 1 is to ensure that the first moment β j of B j exists. We want to emphasize that the singularity analysis of regularly varying distributions 5 and fat-tailed distributions does not differ significantly, since the slowly varying function L(t) can in fact be treated as a constant. The analysis can be easily extended to regularly varying tails, using Theorem 5 in [13]. However we prefer not to, since this only leads to more complicated notations, without making the analysis different.
For type-2 distributions, we restrict ourself to type-2 distributions with an LST similar as (46), namely as follows where α j < 0 and k j is a real constant such that k 2 Γ(α j ) > 0.

5.
Asymptotics of the queue length. In order to obtain asymptotics for the queue length, we perform a singularity analysis of the pgf N 1 (z), given by (38), as described in Section 4. The singularities of N 1 (z) are those of A 1 (z, 1), A 2 (z, 1) and the possible zeros of the denominator. Let us define ξ 1 and ξ 2 as the dominant singularity of A 1 (z, 1) and A 2 (z, 1) respectively. Using Definition (5) it follows that Lastly, we define ξ as the dominant singularity of N 1 (z).

5.1.
Type-1 service time distributions for class-1 customers. Because in this case, 1 (1, 1) < 1 and A 1 (z, 1) is a strictly increasing, convex function on the positive real axis, the denominator of (38) has a simple zero z 0 , such that 1 < z 0 < ξ 1 . The dominant singularity ξ of N 1 (z) is then either z 0 or ξ 2 . 5.1.1. Type-1 service time distributions for class-2 customers. We have two possible cases depending on the minimum of z 0 and ξ 2 . Case 1: z 0 < ξ 2 . The dominant singularity of N 1 (z) is z 0 . The singular expansion of N 1 (z) in the neighbourhood of z 0 singularity can be obtained from its Laurent series and is given by Transferring to coefficients yields Remark that the asymptotic distribution of N 1 consists exclusively of an exponential term z −n 0 , besides a constant factor. We say that N 1 has, in this case, a geometric behaviour. Case 2: ξ 2 < z 0 . The dominant singularity of N 1 (z) is now ξ 2 . The singularity depends on the LST B * 2 (s). Without loss of generality, we write that B * 2 (s) has a singular expansion at s = τ 2 of the form where α 2 > 0 and k 2 is a real constant wherefore k 2 Γ(α 2 ) > 0. Because A 2 (z, 1) = B * 2 (λ 1 (1 − z)), we find that Notice that A 1 (z, 1) is analytic at z = ξ 2 . It is sufficient to approximate 1/(A 1 (z, 1)− z) as 1 where we write f (s) = O(g(s)) as s → s 0 if lim s→s0 sup |f (s)/g(s)| < ∞.
Using the expansion of A 1 (z, 1) and A 2 (z, 1) at z = ξ 2 , we obtain that Using Theorem 4.1, we conclude that If α 2 = 1, the distribution of N 1 is, in this case, is the product of a power law n α2−1 and an exponential ξ −n 2 . We say that N 1 has, in this case, a semi-geometric behaviour. Note that α 2 = 1 means that −τ 2 is a simple pole of B * 2 (s). If this is the case, N 1 has a geometric behaviour.

5.1.2.
Type-2 service time distributions for class-2 customers. The singularity analysis yields the same result as the previous section. If z 0 < ξ 2 , the analysis is completely the same. If ξ 2 < z 0 , we make use of (47), it follows that where now −α 2 > 0. This yields a singular expansion for N 1 (z) of the form Where x j are some constants. Applying Theorem 4.1 results again into (54).

5.1.3.
Type-3 service time distributions for class-2 customers. Since B 2 is a type-3 distribution, we have that ξ 2 = 1, hence ξ = 1. Further, by making use of (46) we find that The numerator of (38) can therefore be estimated for z → 1 as Since the denominator of (38) vanishes for z = 1, we now estimate the denominator of (38) as p 1 (1 − λ 1 β 1 )(1 − z) for z → 1. If we now cancel the common factor (1 − z) in the estimated numerator and estimated denominator of (38) near z = 1, we find the following singular expansion for N 1 (z) at z = 1 Transferring to coefficients yields The distribution of N 1 consists in this case solely of a power term n α2 .

5.2.1.
Type-1 and type-2 service time distributions for class-2 customers. Case 1: z 0 does not exist, and ξ 1 < ξ 2 . We deduce that For z → ξ 1 , we estimate the factor 1/(A 1 (z, 1) − z) as follows where we in the second step expanded the denominator and kept only the terms of order −α 1 or lower (remark that −α 1 > 0). Plugging this expansion into (38) yields Transferring to coefficients then yields Case 2: z 0 does not exist, and ξ 2 < ξ 1 . This case yields the same singularity analysis as in Case 2 of Section 5.1.1. Case 3: z 0 exists. This case yields the same singularity analysis as in Section 5.1.1 for both z 0 < ξ 2 and ξ 2 < z 0 .

5.2.2.
Type-3 service time distributions for class-2 customers. The singularity analysis and asymptotics are the same as in Section 5.1.3.

5.3.
Type-3 service time distributions for class-1 customers. If the class-1 service time distributions are of type-3, the denominator no longer possesses a simple zero z 0 , wherefore z 0 > 1. Because ξ 1 = 1, the dominant singularity of N 1 (z) is that of A 1 (z, 1), i.e. ξ = ξ 1 . If we assume that the LST of B 1 is of the form (46), it follows that Consequently, for z → 1.

Dividing this expansion by Expansion
Hence, applying Theorem 4.1, we obtain that The asymptotic distribution of N 1 is again a pure power-law.

5.3.2.
Type-3 service time distributions for class-2 customers. We have that ξ 2 = 1 and the singular expansion of A 2 (z, 1) at z = 1 is given by (57). The expansion for the numerator is the same as (58), while the expansion for the denominator is given by (63). We divide (58) by (63) and cancel the common factor (1 − z).
If we now expand the denominator and keep only the dominant terms, the result depends on the minimum between −α 1 and −α 2 . Namely, if −α 1 < −α 2 , we obtain the singular expansion (64). If −α 2 < −α 1 , we obtain the singular Expansion (59). Suppose now the boundary case α 1 = α 2 . In this case, the singular expansion at z = 1 yields which leads to Notice that constant factor in this expression is the sum of the constants in Expressions (60) and (65).
6. Asymptotics of the orbit length. In this section we derive the asymptotic distribution of N 2 . Expression (39) contains the implicit defined function Y (z). Therefore, the singularity analysis of N 2 (z) is more involved than the singularity analysis of N 1 (z). It will be convenient to define the following function so that N 2 (z) can be rewritten as The possible singularities for N 2 (z) arise from the functions Y (z) and A 2 (Y (z), z) and the zeros of z − D(z). In order to perform a singularity analysis on N 2 (z), we first have to determine the dominant singularities of Y (z) and A 2 (Y (z), z). We emphasize that both Y (z) and A 2 (Y (z), z) are in fact pgfs (c.f. [23]). To make the analysis more clear, we introduce some more new variables, ζ := Dominant singularity of N 2 (z) .
We have to determine successively ζ 1 , ζ 2 , ζ 3 and ζ. First, we will show that the dominant singularity of the implicit defined function Y (z) is of square-root type, if B 1 is a light-tailed distribution. If B 1 is a fat-tailed distribution, the dominant singularity of Y (z) will be that of B * 1 (s). In case of type-2 distributions, both situations can occur. Secondly, to determine ζ 2 , we introduce the auxiliary function ). Notice that the function h(z) is a decreasing function for z ∈ [0, ζ 1 [ with h(1) = 0. There are three possible cases for the singularity of A 2 (Y (z), z), depending on the location of −τ 2 compared to h(ζ 1 ) [14, Ch. VI.9]: 1. Supercritical case: h(ζ 1 ) < −τ 2 . By the intermediate value theorem, there exists a θ, 1 < θ < ζ 1 , such that h(θ) = −τ 2 . The dominant singularity of A 2 (Y (z), z) is θ, thus ζ 2 = θ. Around this point, Y (z) is analytic because θ < ζ 1 . A singular expansion of A 2 (Y (z), z) around θ is obtained by combining the singular expansion of B * 2 (s) and the Taylor expansion of Y (z) at θ.

Subcritical case
When z increases from 0, the singularity of h(z) is achieved before the function h(z) achieves the critical value −τ 2 . Therefore the singularity of A 2 (Y (z), z) is now driven by the singularity of h(z). Hence, ζ 2 = ζ 1 and the singular expansion of A 2 (Y (z), z) at ζ 1 is obtained by combining the Taylor expansion of B * 2 (s) and the singular expansion of Y (z) at ζ 1 . 3. Critical case: h(ζ 1 ) = −τ 2 .
We have that ζ 2 = ζ 1 and A 2 (Y (ζ 2 ), ζ 2 ) = B * 2 (−τ 2 ). The singular expansion is obtained by composition rules of the singular expansions. Remark that not all cases are possible if ζ 1 = 1 or τ 2 = 0. Indeed, if we have that ζ 1 = 1, the supercritical case can never occur, because h(1) = 0 and τ 2 ≥ 0. If we have that τ 2 = 0, the subcritical case never occurs, because h(ζ 1 ) ≤ 0. Except when both service time distributions are of type-3, we will never treat the critical case because this case only occurs for a very specific choice of parameters and service time distributions.
Thirdly, from the definition of D(z), see (68), it is clear that ζ 3 = min(ζ 1 , ζ 2 ). But from the above considerations we already know that ζ 2 ≤ ζ 1 . Finally, if there exists a z 0 ∈]1, ζ 3 [ such that z 0 − D(z 0 ) = 0, it necessary holds that ζ = z 0 . If this is not the case, we have that ζ = ζ 3 because of the previous reasonings.
6.1. Type-1 service time distributions for class-1 customers. In this section we assume that the service time distribution of class-1 customers are type-1 distributions. We start with the singularity analysis of Y (z). This is done by applying the Smooth implicit-function scheme [14, Th. VII.3, pg 468]. Lemma 6.1. We have that ζ 1 = ω, where ω is a square-root singularity. The tuple (ω, Y (ω)) is the positive solution (z, w) of the following characteristic system of equations, The expansion of Y (z) in the neighbourhood of ω is equal to the expansion being valid in a ∆-domain. The constant c Y is given by .
1. If a solution of the characteristic system exists within the positive quadrant, where A 1 (z, w) is analytic, then it will be necessary unique, because the Taylor coefficients of Y (z) cannot admit two asymptotic expressions with different parameters [14].
In the remainder of this section, we use notation ζ 1 instead of ω as the dominant singularity of Y (z). 6.1.1. Type-1 service time distributions for class-2 customers. Supercritical case. In the supercritical case and since B 2 is a type-1 distribution, we have that A 2 (Y (ζ 2 ), ζ 2 ) = B 2 (−τ 2 ) = ∞. Therefore, D(ζ 3 ) = ∞. This property, together with the fact that D(z) is an increasing, convex function on the positive real-axis, D(1) = 1 and implies that z − D(z) has a simple zero z 0 wherefore 1 < z 0 < ζ 3 . The singular expansion of N 2 (z) in the neighbourhood of ζ = z 0 is given by Transferring to coefficients yields Subcritical case. In this case we have that the dominant singularity of Y (z) is transferred to A 2 (Y (z), z) and D(z). It follows that D(ζ 3 ) < ∞. Hence, unlike in the supercritical case, the simple zero z 0 of z − D(z) does not necessarily exist. Therefore, there are two possibilities for ζ. Either ζ = z 0 , if z 0 exists, or ζ = ζ 3 . The first case will lead to the same coefficient asymptotics as (80). In order to obtain the subexponential factor in the latter case, we consecutively write down singular expansions of the pgfs A 2 (Y (z), z) and D(z) at z = ζ 1 . Because we are in the subcritical case, B * 2 (s) is analytic in h(ζ 1 ). Hence, we can write A singular expansion of D(z) is simply obtained by using those of (81) and (76). This yields The constant c D is given by We can now make use of (76) and (82) to obtain a singular expansion of N 2 (z). For z → ζ 2 , we have that The constant c 1 is given by Using Theorem 4.1, we obtain In case that both service time distributions are of type-1, we can conclude that the asymptotic distribution of N 2 either has a geometric behaviour (z −n 0 ) or a semigeometric behaviour (n −3/2 ζ −n 1 ), depending on the existence of z 0 .
6.1.2. Type-2 service time distributions for class-2 customers. Supercritical case. In the supercritical case, we have that ζ 3 = ζ 2 < ζ 1 and h(ζ 2 ) = −τ 2 . Hence, ζ is either the simple zero in ]1, ζ 2 [ of z − D(z), if any, or ζ = ζ 2 . The first case yields the same result as (80). In the rest of this paragraph, we will assume that ζ = ζ 2 . We assume that B * 2 (s) has a singular expansion near −τ 2 like (47). Using this expansion, we can write for h(z) → −τ 2 Notice now that h(z) is analytic in ζ 2 and admits therefore a Taylor expansion around ζ 2 . Inserting the Taylor expansion of h(z) around z = ζ 2 into (87) yields where the x i are some constants and where we used that h (z) = −(λ 1 Y (z) + λ 2 ). Further, we can write Using a similar procedure as in the previous subsection, it follows that The constant c 2 is given as (90) Finally we obtain the asymptotic distribution of N 2 using the transform theorem, Subcritical case. The analysis for the subcritical case is completely the same as for the case of a type-1 distribution for the class-2 service time. This is due to the fact that in the subcritical case, the singularity is driven by that of Y (z).
6.1.3. Type-3 service time distributions for class-2 customers. As already pointed out, the subcritical case can never occur if τ 2 = 0. Moreover, because ζ 1 > 1, the scheme is always supercritical. We assume that B 2 is a fat-tailed distribution. Then we can write for z → 1 and (93) From the singular expansion of D(z) near z = 1, we find that with c 3 defined as To obtain this singular expansion of N 2 (z), we cancelled first a common factor (1 − z) in the expansions of the numerator and denominator. Then we expanded the denominator again. We used thatx 1 = −D (1). We have now obtained the following asymptotic: Remark that R * (λ) is incorporated in the constant p 0 , so the retrial time distribution does have an influence on the constant c 3 . Similar as for the priority queue, we observe that in this case the distribution of N 2 follows a power-law n α2 .
6.2. Type-2 service time distributions for class-1 customers. Because B * 1 (s) is no longer infinite in −τ 1 , a solution of the characteristic system (75) in the positive quadrant no longer necessarily exists. See also Appendix A in [24] for a detailed discussion. If a solution exists, the previous results apply. Therefore, in this section we will only discuss the case that Lemma 6.1 does not hold, i.e. there exists a θ, wherefore h(θ) = −τ 1 . Moreover, because Y (θ) = B * 1 (−τ 1 ), we obtain an explicit expression for θ: As before, we will write the dominant singularity of Y (z) as ζ 1 , in fact ζ 1 = θ. Note that A 1 (w, z) is still bivariate analytic in (1, 1), because of τ 1 > 0, and A 1 (Y (1), 1) = 1. Hence, the implicit function theorem ensures us that ζ 1 > 1. We assume the expansion of B * 1 (s) in −τ 1 is like (47). The singular expansion of Y (z) near ζ 1 can then be written as 6.2.1. Type-1 and type-2 service time distributions for class-2 customers. Supercritical case. The supercritical scheme is the same as in Section 6.1.1 and Section 6.1.2. Subcritical case. We have two possible instances for the dominant singularity ζ of N 2 (z). Either ζ = ζ 1 , or ζ = z 0 , where z 0 is a simple zero of z − D(z) contained in the interval ]1, ζ 1 [, if any. The latter case leads to the same asymptotics as (80).
In the following analysis, we assume that such a z 0 does not exist. A singular expansion for A 2 (Y (z), z) at z = ζ 1 is obtained by combining the Taylor expansion of B * 2 (s) at s = −τ 1 and the singular expansion of Y (z), given by (97). This yields for z → ζ 1 . Consequently, where the constant c D is given by With similar techniques as before, we obtain a singular expansion of N 2 (z), The constant c 4 is given as We obtain the following asymptotic Remark that k 1 and Γ(α 1 ) have the same sign by assumption. Further, based on the properties of an LST, it can easily be seen that c D < 0. Finally, ζ 1 − D(ζ 1 ) > 0, exactly because we assumed that z 0 does not exist. From these three remarks one can verify, as a minor check, that the obtained tail asymptotic is positive.

6.2.2.
Type-3 service time distributions for class-2 customers. The subcritical case cannot occur if B 2 is a type-3 distribution. The supercritical case is the same as in Section 6.1.3.
6.3. Type-3 service time distributions for class-1 customers. As explained in Appendix A [24], Lemma 6.1 never holds, hence ζ 1 = 1. As a consequence, the scheme is always subcritical if the class-2 distribution is type-1 or type-2. The scheme is critical in case that τ 2 = 0. Using the Expansion (46), we find the following singular expansion of Y (z) for z → 1 6.3.1. Type-1 and type-2 service time distributions for class-2 customers. In these two cases, B * 2 (s) is analytic in s = 0. Hence, for the singular expansion of A 2 (Y (z), z) for z → 1 we get We also obtain the singular expansion of D(z) for z → 1, with c D defined as Using the expansion of D(z) and Y (z) for z → 1, yields The constant c 5 is equal to where we used thatx 1 = −D (1). Finally, using the transfer theorem, we find that Again, similar as for the priority queue, we conclude that the asymptotic distribution of N 2 follows a power-law n α1 in this case.
6.3.2. Type-3 service time distributions for class-2 customers. In this particular case, both service time distributions are of type-3, i.e. their dominant singularities are zero. The scheme is now critical and the expansion of A 2 (Y (z), z) for z → 1 will depend on the minimum of −α 1 and −α 2 , namely The first case −α 1 < −α 2 corresponds to (105). The asymptotics will yield the same result as (110). The second case −α 2 < −α 1 corresponds to (92). The corresponding asymptotics will yield the same result as (95). The last case −α 1 = −α 2 is a boundary case. However, since it is an important boundary case (both service time distributions are equally "fat"), we treat it explicitly. It can be shown that, if −α 1 = −α 2 , the asymptotic expansion of D(z) in z = 1 is given by where the auxiliary constant c D is given by Using similar mathematically manipulations as before, we find that where the constant c 6 is given by Notice that c 6 = c 3 + c 5 . Finally, we obtain that 7. Performance measures. In this section we list some important performance measures of the queueing system. The probability that the server is empty is given by Q(1), It is worth noting that this is not equal to the probability that the system is empty, which is given by p 0 , see (34). Also, this performance measure indicates that the probability that a new arrival will be served immediately is independent of the retrial time distribution. The mean length of the priority queue is found by taking the first derivative of (38), yielding .
The mean orbit length is found by taking the first derivative of (39), which leads to From the derived asymptotics in Sections 5 and 6, we can derive tail asymptotics, i.e. P[N j > n] as n → ∞. It holds that If the dominant singularity is strictly greater than one, we have that In case the dominant singularity is equal to one, the singular expansion of N j (z) is of the form 1 + −α −1 j=1 x j (1 − z) j + cΓ(1 + α)(1 − z) −α−1 . Therefore, after cancelling (1 − z), the singular expansion of (120) is of the form Using Theorem 4.1, we obtain that where we used that Γ(α + 2) = (α + 1)Γ(α + 1). Notice that −α > 1 by assumption. Hence (122) is always positive. 8. Numerical examples. To conclude this paper, we show and discuss some numerical examples. We assume the service times of type-1 customers to be exponentially distributed with mean β 1 , i.e.
In order to study the influence of different service times, we assume the service times of type-2 customers to be of the gamma type with shape parameter 1/2 and mean β 2 , i.e. B * 2 (s) = (1 + 2β 2 s) −1/2 . The reason for this distribution is twofold. The first reason is to have two different service time distributions. The second reason is to have a coefficient of variation 6 that is bigger than one; it is equal to √ 2. The coefficient of variation for B * 1 (s) obviously equals one. Next, we suppose the retrial times follow an Erlang distribution of order two and rate ν > 0, i.e. (125) Finally, we define α as the fraction of class-1 customers load in the overall traffic mix, i.e., with ρ = λ 1 β 1 + λ 2 β 2 as before.  Figure 1 shows the mean queue length and mean orbit length versus the mean service time of class-2 customers for ρ = 0.7, β 1 = 1, ν = 3 and α = 0.25, 0.5 and 0.75. From Figure 1, it can be seen that the mean queue length increases linearly with β 2 , while the mean orbit length decreases with increasing β 2 . The decrease of the mean orbit length is caused by the decrease of λ 2 when β 2 increases (ρ is kept constant), hence less class-2 customers arrive to the system. The increase in the mean queue length can be explained as follows. When β 2 increases, more class-1 customers will arrive during the time periods that a class-2 customer is being served. Due to the non-preemptive priority policy, arriving class-1 customers cannot interrupt the service of class-2 customers. Hence, they have to wait until this class-2 customer leaves the system, which results in a larger mean queue length.
In Figure 2 the mean queue and orbit lengths are shown as a function of the (mean) class-1 service time with ρ = 0.7, β 1 = 1.5, ν = 3 and α = 0.25, 0.5 and 0.75. On the one hand, we observe that the mean queue length decreases with increasing β 1 . This is because less class-1 customers arrive when β 1 increases (since λ 1 β 1 is kept constant). The mean orbit length on the other hand first decreases rapidly and then increases almost linearly, for increasing β 1 . The latter effect is due to the priority policy, the former to the retrial policy. To motivate this, we present in Figure 3 the same model as in Figure 2, but where the retrial times are zero, i.e. R * (s) = 1. In this way, we obtained the ordinary non-preemptive priority queueing model. The mean queue lengths are the same in both figures, because the mean queue length is independent of the retrial time distribution. The decrease for small β 1 observed in Figure 2 can be explained as follows. The probability of a successful retrial R * (λ 1 + λ 2 ) increases, for decreasing λ 1 (for fixed λ 2 ). Intuitively, this is because less class-1 customers arrive and therefore there is less competition between the retrial and the arriving customers. Mathematically, this can be seen by taking the first derivative of (125). Consequently, if β 1 increases (hence λ 1 decreases), the effective time until a successful retrial decreases. This initially results in a decrease of the orbit length because the retrials gets more successful. This positive effect continues until β 1 achieves a critical value. When the (mean) class-1 service time gets longer, the variance in the lengths of the periods that both queue and server are empty grows, which implies higher (mean) orbit lengths [23]. The variances of the queue and orbit lengths can be derived from the pgfs (38) and (39). In Figure 4, the variances of the queue and orbit lengths are shown as functions of the load when α = 0.25, 0.5 and 0.75 respectively.  We now take a look at the tail probabilities of the orbit length. We assume that both service time distributions are of type-1, with τ 1 = 1/β 1 and τ 2 = 1/(2β 2 ). Figure 5 shows which scheme applies, according to the analysis in Section 6.1.1, for λ 1 = λ 2 = 0.3 and ν = 3. The line β 1 +β 2 = c, for a certain constant c, indicates the stability border for the queueing system. The parameterspace (β 1 , β 2 ) is first divided into two parts, the supercritical and the subcritical scheme. In the supercritical scheme, the tail of the orbit has a geometric behaviour. In the subcritical scheme, there are two cases, geometric or the product of a power law with a geometric (semi-geometric). This yields the second partition of the parameterspace (β 1 , β 2 ). For a large set of the parameters β 1 and β 2 , the tail behaviour is thus geometric. Which scheme applies, i.e. the supercritical or the subcritical one, is independent of the retrial time distribution. However, which case applies within the subcritical scheme does depend on the retrial time distribution, because this is determined by the existence of a positive real zero greater than one of the denominator of (39).
To illustrate this, we show in Figure 6 which scheme applies if the retrial times are zero. Notice, that in this case the stability condition is different, which explains the move upwards of the stability border in Figure 6. We indeed see that the area that indicates semi-geometric behaviour is changed compared to Figure 5. Finally, Figure 7 shows the tail probabilities of the orbit length for λ 1 = λ 2 = 0.3, β 1 = 1.2, β 2 = 1 and ν = 3, when one of the service time distributions is changed to a fat-tailed distribution, for α 1 = α 2 = −2.5 (see also (46)). The straight line in Figure 7 is the case when both service time distributions are of type-1. From Figure  5 we know that the tail behaviour is then geometric. As expected, the tail of the orbit length is much heavier if (one of) the service time distributions is fat-tailed.  9. Conclusion. In this paper, We have investigated a priority retrial queue with constant retrial policy and different service time distributions between class-1 and class-2 customers. We were able to derive expressions for the pgfs of the stationary queue and orbit lengths. The main contribution of this paper is the singularity analysis of these pgfs, for several different types of service time distributions (including a large subclass of heavy-tailed distributions). The result of this analysis is that, it characterizes the effects of the service time distributions and the retrial time distributions on the queue and orbit lengths distributions. It turns out that the sort of asymptotic behaviour for the orbit length is strongly determined by the service time distributions. However, if both class-1 and class-2 service time distributions are not heavy-tailed, then also the retrial time distribution has an influence on the asymptotic behaviour of the orbit length, since the pgf of the orbit length may then possess a simple pole as dominant singularity. The existence of this simple pole is determined by the service time distributions, but also of the retrial time distribution if the scheme is subcritical. In case one of the service time distributions is a fat-tailed distribution, both queue and orbit length follow a power-law distribution.