Heavy-Trafﬁc Comparison of a Discrete-Time Generalized Processor Sharing Queue and a Pure Randomly Alternating Service Queue

: This paper compares two discrete-time single-server queueing models with two queues. In both models, the server is available to a queue with probability 1/2 at each service opportunity. Since obtaining easy-to-evaluate expressions for the joint moments is not feasible, we rely on a heavy-trafﬁc limit approach. The correlation coefﬁcient of the queue-contents is computed via the solution of a two-dimensional functional equation obtained by reducing it to a boundary value problem on a hyperbola. In most server-sharing models, it is assumed that the system is work-conserving in the sense that if one of the queues is empty, a customer of the other queue is served with probability 1. In our second model, we omit this work-conserving rule such that the server can be idle in case of a non-empty queue. Contrary to what we would expect, the resulting heavy-trafﬁc approximations reveal that both models remain different for critically loaded queues.


Introduction
The present article is motivated by mathematical challenges that arise in the study of queueing systems that involve two types of customers that compete for the same server. Such two-class queueing systems give rise to Markov processes with two-dimensional state spaces. When the storage capacity for both types of customers is infinite, the state space is unbounded in two dimensions which causes the analysis to be substantially different from that of systems that involve but one queue. The earliest and most well-known approach to tackling two-class queueing systems is the boundary value method [1][2][3]. This method relies on a generating function approach to analyze the steady-state distribution of the twodimensional Markov process. An expression for the joint generating function of the system contents is found by solving a boundary problem for analytic functions. Unfortunately, this approach is notorious for its cumbersome (exact) integral formulas. While these formulas are exact, it is far from easy to extract performance measures from them, see for example Part IV in [1], and [4,5]. While the calculations for the general case are cumbersome, in some fortunate cases, usually with a simple arrival process, there is little to no numerical work needed; this is for example possible in the discrete-time 2 × 2 switch model [6,7].
In this paper, we study two simple and conceptual discrete-time queueing models involving two different queues that get an equal share of the single-server capacity. The server is shared in the following way: At each service opportunity, the server chooses a queue with equal probability. Based on what happens when an empty queue is chosen, we consider two different models. In model I, no one gets service if an empty queue is chosen even if the other queue is non-empty. In contrast, in model II, the server will always serve a non-empty queue in order to keep the whole system work-conserving. From a queueing theoretical point of view, we are interested in the influence of this different rule on the (joint) performance measures of the queueing system. Clearly, some capacity of the server is lost in model I. The rationale for model I is that the server has either no information about the queue lengths, or that obtaining this information is a costly operation. The latter can, for example, be the case in applications where energy efficiency is an issue. The service discipline in model I does not depend on queue lengths, and hence it is not work-conserving. One of the first papers having this property that is similar to model I is [8], where a single server visits a finite number of stations for fixed time periods. More recent queueing models with a single sever visiting multiple queues in a non-work-conserving matter have been studied in [9][10][11][12][13]. In our earlier paper [14], model I was introduced in a more general setting than in this paper. We used a numerical method, based on tail probability asymptotics, to gain some insights in the performance measures of the system. Model II is what is often called the discrete-time Generalized Processor Sharing (GPS) model [15]. In the continuous-time queueing literature, this is perhaps better known as a coupled processor model [16,17]. Model II, in a more general setting, has been approximately analyzed in [18], using the power series approximation for generating functions.
In the case of a work-conserving service discipline like in model II, the distribution of the total number of customers is easy to calculate since the total system is identical to a single-server single-queue model. In model I, this property is inevitably lost. However, in contrast, the marginal distributions of the number of customers present in both queues are easy to calculate for model I. Indeed, each queue is identical to a single-server queueing model with server interruptions (we kindly refer the reader to Chapter 3.2 in [19] for a comprehensive treatment of discrete-time models with server interruptions). These two observations allow for comparing both models by means of a simple mean value analysis. We are however also interested in joint performance measures, in particular the coefficient of correlation between the two queue lengths. Obtaining joint performance measures brings us back to the machinery of the boundary value method, of which it was already argued that it rarely provides explicit, let alone easy-to-calculate, performance measures.
In this paper, we consider the heavy-traffic limit of the joint system content distribution of both models. In our case, this means that we let the mean arrival rate of the models go to its critical value of instability such that the queues are nearly saturated. The reason for considering the heavy-traffic limit is two-fold. The first reason is that by considering the heavy-traffic limit, the boundary value method can be applied in a simpler manner compared to the non-heavy-traffic case, such that the numerical work is limited (model II) or even not necessary (model I). Moreover, since it is typical for heavy-traffic results to be rather insensitive to the exact form of the arrival (and service) process [20], we can assume a general batch arrival process for our queueing models. The second reason is that in heavy traffic, we have the interesting question relating to if both models converge to the same model. It is reasonable to think that the answer to this question is affirmative since it is expected that the queues empty less and less near saturation, and hence the only rule that sets the two models apart applies less frequently. However, we show in this paper that this is not true.
Earlier studies that consider the heavy-traffic limit of similar, continuous-time, queueing models are [21][22][23]. It is worth mentioning that these papers use a heavy-traffic diffusion approximation (see for example also [24,25] for a more detailed explanation of this method). Broadly speaking, with the diffusion approximation one replaces the balance equations by a diffusion equation with appropriate boundary conditions. Recently, in the doctoral thesis Chapter 4 in [26], a heavy-traffic limit is obtained using transforms, similar as in the classic work of [20], but with the added difficulty that a boundary value problem problem for analytic functions has to be solved. However, a nice boundary value problem is obtained, where the boundary consist of a parabola, for which an explicit solution is obtained. We note that the model in Chapter 4 in [26] is a continuous-time queueing model, referred to as the random-time limited polling model with two queues, and is in a sense somewhat similar to the (discrete-time) model I in the current paper.
The structure of the remainder of this paper is as follows. In Section 2, the model under investigation is described in detail. Section 3 presents the problem statement and lists the main results of this paper. In Section 4, the analysis for model I is carried out. Finally, Section 5 is devoted to the analysis of model II. In both Sections 4 and 5, we also discuss numerical results.

Mathematical Model and Preliminary Results
We consider a discrete-time single-server queueing system with two infinitely sized queues, say queue 1 and queue 2, and two independent input lines. Customers arriving at queue 1 and queue 2 are referred to as type-1 and type-2 customers, respectively. Time is assumed to be slotted. The numbers of type-1 and type-2 arrivals during slot k are denoted by a 1,k and a 2,k , respectively. The sequences of discrete random variables (a 1,k ) k∈N and (a 2,k ) k∈N are assumed to be i.i.d. and are specified by a common probability generating function (pgf) A(z), i.e., We thus assume that both types of customers have the same arrival process, i.e., a 1,k and a 2,k have the same probability distribution. Moreover, we assume that these two arrival processes are independent, i.e., the random variables a 1,k and a 2,k are independent. For further use, let λ denote the mean numbers of type-1 and type-2 arrivals per slot, At the beginning of every slot, the single-server selects either queue with probability 1 2 . In case a non-empty queue is selected, a customer of this queue gets served. We assume that the service of each customer type requires exactly one slot, regardless of whether the customer is of type 1 or type 2. Based on what happens when an empty queue is selected, we make the distinction between the following two service disciplines.
Service discipline 1 (Non-work-conserving policy). We assume that when an empty queue is chosen, no service occurs in that slot, even when the other queue is non-empty. We will refer to this service discipline as the non-work-conserving policy.
Service discipline 2 (Work-conserving policy). When only one of both queues is non-empty and the other is empty, the non-empty queue is served during that slot, even when initially the empty queue was selected. We will refer to this service discipline as the work-conserving policy.
We now write down the key functional equations of the joint PGF of the system contents.

The Non-Work-Conserving Policy
Let u 1 and u 2 indicate the number of type-1 and type-2 customers respectively, in steady-state. We define their joint PGF as: As result of the non-work-conserving property, the queues can be analyzed separately. Both queues are then equivalent to simple single-server queues with geometric service times with mean 2. The arrival load offered by type-j customers is given by 2λ. The stability condition of a single queue is given by: which is therefore also the stability condition of the complete queueing system. Section 2.1 in [14] contains the expression for a functional equation for the joint PGF of the number of type-1 and type-2 customers, for a more general setting with different arrival distributions for type-1 and type-2 customers and the probability that the server selects queue 1 is parametrized. The functional equation in our case reads: with We will refer to K(z 1 , z 2 ) as the kernel of Equation (4). Equation (4) is a functional equation for U(z 1 , z 2 ) as the boundary functions U(z 1 , 0) and U(0, z 2 ) are present in the RHS of the equation. These boundary functions can be found by constructing a boundary value problem for analytic functions. As already mentioned in the introduction, this will require a numerical approach. Therefore, we propose a heavy-traffic approximation that requires only a minimum amount of (or no) numerical computations. We postpone this analysis to Sections 3 and 4. First we give some, more elementary, results of the joint PGF U(z 1 , z 2 ).
It is not difficult to see that the joint PGF U(z 1 , z 2 ) of the system contents is symmetric in z 1 and z 2 , i.e., In particular, this implies that there is only one boundary function present in the functional equation for U(z 1 , z 2 ), since: The PGF S(z) of the total number of customers can be expressed in terms of this boundary function. Substituting z 1 = z 2 = z into (4) yields: On the other hand, the marginal PGFs U(z, 1) and U(1, z) of the number of type-1 and type-2 customers do not depend on U(z, 0). Moreover U(z, 1) and U(1, z) are identical and given by: A first important performance measure that can be derived from this PGF, is the probability of an empty queue at the beginning of a random slot: Further, all the marginal moments of interest of the random variables u 1 and u 2 can thus be computed. For example, the mean is given by:

The Work-Conserving Policy
Let v 1 and v 2 be the number of type-1 and type-2 customers in steady-state. Their joint PGF is defined by: This model is a special case of the model studied in [18]. Introducing the assumptions of the present paper (i.e., single slot service times, the same arrival distributions and weights), the following functional equation for V(z 1 , z 2 ) is obtained: where K is given by (5). Note that both functional equations have the same kernel.
In case of a work-conserving service discipline and single-slot service times, the distribution of the total number of customers is easy to calculate since the total system is identical to a single-server model with a single queue. Indeed, set z 1 = z 2 = z in (12) to obtain the relation: and then take the limit z = 1 to obtain (using l'Hôpital's rule): All the moments of the total system content v 1 + v 2 can therefore be computed. For example, we obtain that the mean total number of customers in the system is given by: It is not difficult to see that the joint pgf V(z 1 , z 2 ) is symmetric in z 1 , z 2 , i.e., As a consequence, there is only one boundary function present in (12). Moreover, the marginal pgfs V(z, 1) and V(1, z) of the system contents of type-1 and type-2 customers are identical and given by: where we used (14). It is worth noting that (14) indicates the probability that the total system is empty at the beginning of a random slot. Due to the work-conserving service discipline, the stability condition of the system is naturally given by 2λ < 1, i.e., (3). This is the same stability condition as for the non-work-conserving policy. This is not a coincidence, and it is precisely the reason that it was assumed a symmetrical system. We emphasize that in the case of a non-symmetrical system, both systems have a different stability condition (cf. Section 2 in [14] and Section 4.2 in [15]).

Problem Statement and Main Results
We are interested in the influence of the service discipline on the correlation structure between the number of type-1 and type-2 customers, when λ is near the critical value 1 2 . The simplest possible performance measure to quantify the correlation structure is the correlation coefficient between the number of type-1 and type-2 customers. The objective of our analysis is to obtain expressions for lim Our approach is to obtain the joint Laplace-Stieltjes transform (LST) of the scaled system contents (1 − 2λ)u 1 and (1 − 2λ)u 2 as λ ↑ 1 2 and the joint LST of the scaled system contents Notice that for two random variables X and Y: A functional equation for the two limiting joint LSTs can easily be obtained from (4) and (12). The kernel in these new functional equations is much simpler than K(z 1 , z 2 ), as will be shown in Section 4.
In this paper, we assume that A (1) has a finite limit as λ ↑ 1 2 . We define: This mathematical quantity can be expressed in terms of the physical quantity var[a j,k ], since var[a j,k ] = A (1) + λ − λ 2 . Therefore, we define σ 2 h as the variance of the number of type-j arrivals per slot, as λ ↑ 1 2 , i.e., Notice that σ 2 h is the same for both customer types, because of the symmetry of the arrival process. Since A(z) is a power series in z with positive coefficients, we have that λ 11 ≥ 0. Consequently, we have the following important lower bound for σ 2 h : When λ approaches its critical value 1 2 in the non-work-conserving model, it is expected that the probability to select an empty server tends to zero because Pr[u 1 In that perspective, the non-work-conserving model resembles the work-conserving model when λ is close to 1 2 , because the number of wasted slots in the non-work-conserving model decreases. However, a simple mean value analysis already reveals that even in the limit λ ↑ 1 2 , the work-conserving model is still strictly more efficient. Indeed, from (10) we have that when λ ↑ 1 2 , the mean total scaled system content (1 − 2λ)(u 1 + u 2 ) tends to: while from (15), we easily obtain that: Hence, we have that for λ ↑ 1 2 : where f (x) ∼ g(x) means that f (x) g(x) tends to 1 as x → x 0 . Consequently, we have shown that the mean total system contents are asymptotically not equivalent.
In Sections 4 and 5, we obtain the limiting LSTs of the scaled system contents. From these expressions, the correlation coefficient for λ ↑ 1 2 under the different service discipline is obtained. The final, complete expressions are given in the theorem below. Theorem 1. The correlation coefficient between the number of type-1 and type-2 customers in the system with the non-work-conserving policy, for λ ↑ 1 2 , is given by: The correlation coefficient between the number of type-1 and type-2 customers in the system with the work-conserving policy, for λ ↑ 1 2 , is given by: The result of Theorem 1 is discussed in Sections 4.4 and 5.4, but it has already been mentioned that both correlation coefficients are significantly different when compared to each other.

The Non Work-Conserving Policy in Heavy-Traffic
The purpose of the analysis in this section is to obtain the joint LST of the scaled system contents (1 − 2λ)u 1 and (1 − 2λ)u 2 as λ ↑ 1 2 . The joint LST of (1 − 2λ)u 1 and (1 − 2λ)u 2 is: Hence, from (4) we obtain an equation for U λ (s 1 , s 2 ): We assume that the following limit exists: The subscript h is to indicate that the LST corresponds to the heavy-traffic limit. We also define: as the limiting LST of the total scaled system content.
The boundary function U(e −s(1−2λ) , 0), as λ ↑ 1 2 , can be written in terms of S h (s). From (7), we obtain that: We have to use the change of variables z = e −s (1−2λ) in this expression and take the limit for λ ↑ 1 2 . However, we have to be careful with derivatives of A(e −s(1−2λ) ) with respect to λ since A(·) itself also depends on λ. For ease of presentation, we therefore make this dependency on λ explicit and write A(λ, z), instead of A(z). Therefore, if we substitute z = e −s(1−2λ) in (31) and take the limit for λ ↑ 1 2 we obtain that: Note that A(λ, 1) = 1 for all λ, because of the normalization condition of a PGF. Hence, e −s(1−2λ) − A 2 (λ, e −s(1−2λ) ) goes to zero as λ ↑ 1 2 . Therefore, we will have to apply l'Hôpital's rule at least once to the limit in (32). As we will show further, the following partial derivatives will occur: with whereby k n represents a series consisting of n consecutive k's. Perhaps the easiest way to obtain (33)-(37), is to look at a series expansion of A(λ, z) about z = 1. Since the first two moments of A(z) exist, we can write the following series expansion of A(λ, z) about z = 1 for fixed λ: If we differentiate (39) with respect to λ once or twice, all terms have a common factor (z − 1). Hence, by then evaluating at z = 1, we get (33) and (34). Equations (35) and (36) follow from the first two moments of A(z), where we used definition (1) to simplify (35). Finally, Equation (37) can be found by first taking the derivative of (39) with respect to z, yielding λ + O(z − 1), and then taking the derivative with respect to λ. Thereafter, by evaluating at z = 1 we obtain (37).
Returning to the limit in (32), we have that: In the second equality we used l'Hôpital's rule again, since (33) and (35) imply that the numerator in the first equality goes to zero. In the third equality, we used (33)-(37) and (18). The final equality follows from (19). Substituting this result into (32) gives us: We can now proceed with the determination of a functional equation for U h (s 1 , s 2 ). Dividing (28) by e −s 1 (1−2λ) − 1 e −s 2 (1−2λ) − 1 , taking the limit λ ↑ 1 2 and using (41), we obtain that: We will refer to: as the kernel of the functional Equation (42). It is easily seen that the marginal LSTs U h (s, 0) and U h (0, s) can be obtained by choosing either {s 1 = s, s 2 = 0} or {s 1 = 0, s 2 = s} in Equation (42). We obtain that: which is the LST of the distribution of an exponential random variable with mean σ 2 h + 1 4 . Clearly, we have the same result for U h (0, s).

Areas of Convergence
In this section, we examine in which region(s) the LSTs U h (s 1 , s 2 ) and S h (s) are analytic. From Laplace transform theory, we have that U h (s, 0) is analytic for Re[s] > − 4 1+4σ 2 h and so we may conclude that U h (s 1 , s 2 ) is at least analytic in the region: Since U h (s 1 , s 2 ) is joint analytic inside this region, this also holds for K h (s 1 , s 2 )U h (s 1 , s 2 ). Finally, using functional Equation (42)

Solution of the Functional Equation
In this section, we solve the functional Equation (42) for U h (s 1 , s 2 ). In order to solve this equation, we have to determine the function S h (s). This can be done by exploiting the fact that when K h (s 1 , s 2 ) vanishes for a certain pair (s 1 , s 2 ) for which U h (s 1 , s 2 ) is finite, it must be that the RHS of (42) vanishes for that pair (s 1 , s 2 ).
Observe that K h (s 1 , s 2 ) is a quadratic polynomial both in s 1 and s 2 , where s 1 and s 2 may be complex-valued. For any given value s 1 , there exist exactly two values, says 2 and s 2 , such that K h (s 1 ,s 2 ) = 0, K h (s 1 ,ŝ 2 ) = 0. Using the well-known relations for the product and sum of the roots of a quadratic equation, we always have that: since K h (s 1 , s 2 ) = (1 + 4σ 2 h )s 2 2 + (4 − 2s 1 )s 2 + 4s 1 + (1 + 4σ 2 h )s 2 1 . If U h (s 1 ,s 2 ) and U h (s 1 ,ŝ 2 ) are finite, we obtain from (42) that: so that by eliminating S h (s 1 ): We consider zeros such thatŝ 2 =s 2 . As we will see further, this choice of zeros is sufficient to determine S h (s). Letting: we will first prove the existence and the exact location of such zeros. Subsitutings 2 = x + iy andŝ 2 = x − iy into (46) yields: and Solving (50) for s 1 yields, and substituting this into (49) gives us: Completing the square and dividing by 2 yields the equation of a hyperbola: This hyperbola is shown in Figure 1 for a specific choice of the parameter σ 2 h . Notice that this hyperbola is symmetric with respect to the x-axis and symmetric with respect to the vertical axis y = − 1 2σ 2 h . Moreover, it holds that: , y = 0 are the two solutions of Equation (52) intersecting the x-axis. Hence, for values s 2 = x + iy on the right-branch of this hyperbola, it holds that: Since σ 2 h ≥ 1 4 , a part of this branch is always located in the left half-plane. In particular, s 2 = 0 is located in the interior of the region bounded by the right-branch of the hyperbola. We have thus shown that for all complex values z located on the hyperbola with Equation (52), there exists a unique, positive real s 1 given by (51) such that K h (s 1 , z) = K h (s 1 ,z) = 0. However, in order to guarantee that U h (s 1 , z) and U h (s 1 ,z) are finite, we restrict ourselves to the right-branch of the hyperbola. As a result of (20), it holds that − 4 Consequently, by restricting to the right-branch of the hyperbola, U h (s 1 , s 2 ) is certainly finite and we can safely use Equation (47). Let us denote the rightbranch of the hyperbola by Σ. More precisely, we define: The curve Σ divides the complex plane into two parts. The region onto the right of it, containing the origin, will in the remainder be referred to as the interior region of Σ. Fors 2 = s,ŝ 2 =s, s ∈ Σ, Equation (47) implies that: or, using that S h (s) = S h (s), Observe that σ 2 h + 1 s S h (s) is analytic in the interior region of Σ, except at s = 0 where it has a simple pole. We thus have reduced our problem to that of determining a function that is analytic inside the interior region of Σ, except for a simple pole at s = 0, with prescribed boundary values of its imaginary part (54). The solution to this problem in case the boundary is the unit circle, is given in Section I.3.3 in [1]. Hence, our boundary value problem can be solved using a conformal mapping. Therefore let f 0 be a conformal mapping from the unit disk onto the interior of Σ. Since the real axis is an axis of symmetry of the interior of Σ, it is natural to chose (as we may) f 0 symmetric with respect to the real axis, i.e., f 0 (s) = f 0 (s). We can then rewrite (54) as: Finally, let us denote s 0 as the unique value in the open interval [−1, 1] such that f 0 (s 0 ) = 0. We are now facing the following problem: Find a function σ 2 h + 1 f 0 (s) S h ( f 0 (s)) analytic inside the unit disk, except for a simple pole at s = s 0 , satisfying the boundary condition (55). The solution of this problem is given by Section I.3.3 in [1]: with c 0 , c 1 two unknown constants that still have to be determined. Let f = f −1 0 be the conformal mapping from the interior region of Σ to the unit disk. Then, or Since S h (s) is the LST of a random variable, it is required that S h (0) = 1. Taking the limit s → 0 into the expression above yields: We thus see that ic 1 Further, from the initial value theorem for LSTs, we must have that lim s→∞ sS h (s) is finite, which can only be if c 0 = − 2 f (0) 1− f 2 (0) . As a result, The function f can be explictly obtained. The linear function z 1 = z + 1 2σ 2 h maps the region described by: z = x + iy, to the region: In [27] (page 186, Equation (13)), it is shown that the interior of the right branch of such a hyperbola is mapped to the upper half-plane by the function z 2 = i √ 2 cosh π 2ϕ Arccosh z 1 k , with: Finally, the Möbius function , maps the upper half-plane to the unit disk such that i √ 2 is mapped to 0 and −i √ 2 to infinity (see for example [28], page 175, Equation (17)). The composition of these mappings gives us the mapping f : where we used that: Note that due to the particular choice of the Möbius function z 3 , the hyperbolic tangent tanh shows up. One can verify that the function f (s) is indeed symmetric with respect to the x-axis. Further it is not difficult to obtain that: Substituting f (s) and f (0) into (59) and using some trigonometric identities gives us the following final result for S h (s): The joint LST of the random variables (1 − 2λ)u 1 and (1 − 2λ)u 2 , as λ ↑ 1 2 is then obtained by solving (42) for U h (s 1 , s 2 ) and substituting (65).
We note that, despite the presence of the function Arccosh, S h (s) is indeed analytic in the right half-plane (as it should be for the LST of a continuous random variable). The function Arccosh is the analytic inverse of cosh : {z ∈ C : 0 < Im z < π} → C \ {a : a ∈ R, |a| ≥ 1}. Clearly, Arccosh is not continuous on the branch cut [1, +∞[, since it holds that for x ∈ [1, +∞[, lim z−>x,Im z<0 Arccosh z = − lim z−>x,Im z>0 Arccosh z = −Arccosh x. However since cosh(−z) = cosh(z), the limit of the composite transformation yields in both cases the same result. By the Riemann continuation theorem, the function S h is indeed analytic in the right half-plane.
Finally, remark that S h (s) is the product of the LST of an exponential variable with mean σ 2 h and the function: which satisfies also the normalization property. It is not known to us if this function is the LST of a known distribution.

Calculation of Moments
The determination of all (mixed) moments of the scaled system contents, as λ ↑ 1 2 , can be deduced from the LST U h (s 1 , s 2 ). In this section, we are particularly interested in the coefficient of correlation. First recall that we obtained the marginal LST of the scaled system contents already from the start (44), which was the LST of an exponential distribution. In particular, we have that: By symmetry the mean and variance of (1 − 2λ)u 2 are equal to the mean and variance of (1 − 2λ)u 1 . Therefore, we also obtain the mean of the total scaled system content: Note that we already obtained this result as (21) in Section 3. In addition, note that this expression can be found by taking the first derivative of (65), substituting s = 0 and multiply by −1. Making use of (65), we can compute higher moments of the scaled total system content using the moment-generating property of LSTs. For the second moment, we get: Since for a random variable X it holds that var[X] = E[X 2 ] − E[X] 2 , it easily follows that: Since also for two random variables X, Y it holds that var[X The correlation between the system contents can then be found as: (72)

Examples and Discussions
The results obtained in the previous subsection depend only on the system parameter σ 2 h . Note that ϕ is given by: Hence, if σ 2 h ranges from 1 4 to +∞, it follows that ϕ ranges from: We emphasize that the results from the previous subsection are valid for any arrival distribution with pgf A(z), provided that A (1) has a finite limit for λ ↑ 1 2 . Before discussing the results in general, we first treat the two extreme cases for σ 2 h , i.e., σ 2 h = 1 4 and σ 2 h → +∞.

Bernoulli Arrivals
Let us assume that the number of type-j arrivals within a slot is Bernoulli distributed, i.e., Since A (1) = 0, the lower bound (20) for σ 2 h is an equality, i.e., The parameter ϕ defined in (61) is then given by: Using the fact that cosh(3z) = 4 cosh 3 (z) + 3 cosh(z) for every z, we obtain that the LST S h (s) simplifies to: We thus observe that, in case of Bernoulli arrivals, S h (s) is the product of three LSTs of exponential random variables. This simple expression is completely in accordance with the results obtained in [29]. In [29], we studied the non-work-conserving model under the assumption of (not necessarily symmetric) Bernoulli arrivals. For this particular case of arrivals, we obtained the joint probability distribution of the number of type-1 and type-2 customers, in steady state. The corresponding joint probability generating function in [29] turned out to be a rational function.
The coefficient of correlation between the type-1 and type-2 customers, using (72), is equal to:

Arrivals with Infinite Asymptotic Variance
We have made the restriction that A (1) has a finite limit for λ ↑ 1 2 . As a consequence, this implies that σ 2 h remains finite. However, in the correlation coefficient (72) we still can take the limit as σ 2 h → +∞. For this limit, we have that ϕ → π 2 . Using (72), we obtain that: This result suggests that in case the mean and the variance of the number of arrivals is very high in both queues, the two queues become uncorrelated.

Other Arrival Processes
For other well-known arrival processes, such as the geometrically distribution, binomial distribution, and the Poisson distribution, the LST S h (s) does not simplify to a rational function. Nevertheless, for every value of σ 2 h ∈ 1 4 , +∞ , we can easily compute the correlation coefficient using (72). In Figure 2, we show the correlation coefficient for σ 2 h ∈ 1 4 , 10 , which quantifies the correlation between the number of type-1 and type-2 customers in the system at the beginning of a slot for λ ↑ 1 2 . Looking at Figure 2, we observe that the correlation coefficient is always negative. It is worth noting that in [29] it has been shown that for the special case of two independent Bernoulli arrival processes, the correlation coefficient is negative for all allowed values of λ. The reason for this negative correlation coefficient is the service policy. If a queue is getting longer, this is because either there were a lot of arrivals in the previous slots or because the queue was not being served the previous slots. So in the latter case, if queue 1 is large, it is likely that queue 2 is small. In the heavy-traffic limit λ ↑ 1 2 , the lengths of both queues go to infinity. However, according to Figure 2, the correlation effect of the service policy, still present for σ 2 h , is not too high. Finally, we observe that the correlation coefficient is a strictly increasing function of σ 2 h and goes from − 1 4 to 0 as σ 2 h increases from 1 4 to +∞.

The Work-Conserving Policy in Heavy-Traffic
We now analyze the joint distribution of v 1 and v 2 in case of heavy traffic. More concretely, the purpose is to obtain the joint LST of the scaled random variables (1 − 2λ)v 1 and (1 − 2λ)v 2 as λ ↑ 1 2 . We assume that this limit distribution exists and we define the corresponding LST of this distribution as: The subscript h is again to indicate that the LST corresponds to the heavy-traffic limit. Similarly we define: as the limiting LST of the number of type-1 customers in the system. Obviously because of the symmetry we also have: as the limiting LST of the number of type-2 customers in the system. To obtain a functional equation for V h (s 1 , s 2 ), we can substitute z 1 = e −s 1 (1−2λ) and z 2 = e −s 2 (1−2λ) into (12) and take the limit λ ↑ 1/2. Let us first rewrite the boundary function V(e −(1−2λ)s , 0) as a function of Q h (s). To that end, we solve Equation (16) for V(z, 0) − V(0, 1). This gives us: Next, we divide both sides by 1 − 2λ and use the change of variable z = e −s(1−2λ) . Finally, taking the limit λ ↑ 1/2 we obtain that: The detailed computations to obtain the limit above are omitted. The computations are very similar to the ones in Section 4 to obtain Equation (41), i.e., by carefully taking into account the dependency of A(z) on the arrival rate λ and applying l'Hôpital's rule several times.
A functional equation for V h (s 1 , s 2 ) can now be obtained, similarly as in Section 4. Rewrite (12) as: Next, we substitute z 1 = e −s 1 (1/2−λ) and z 2 = e −s 2 (1/2−λ) into the above. Finally, taking the limit λ ↑ 1 2 , we obtain the following functional equation for V h (s 1 , s 2 ): The LST of the total scaled system content can easily be obtained by substituting s 1 = s 2 = s in Equation (77), yielding: The above LST is the LST of an exponentially distributed random variable with mean σ 2 h .

Solution of the Functional Equation
In this section, we will determine the boundary function Q h (s), and hence the solution V h (s 1 , s 2 ) of the functional Equation (77). We solve Equation (77) in the same manner as we did in Section 4, i.e., using the boundedness of the function V h (s 1 , s 2 ) and the zeros of K h (s 1 , s 2 ) in order to obtain a boundary value problem for the remaining unknown function Q h (s). Regarding the boundedness, we observe from (78) also analytic for: In particular, the marginal LST Q h (s) = V h (s, 0) is analytic for at least Re[s] > − 1 We now proceed as in Section 4. For any given value s 1 there exist exactly two values, says 2 andŝ 2 , such that K h (s 1 ,s 2 ) = 0, K h (s 1 ,ŝ 2 ) = 0. If moreover V h (s 1 ,s 2 ) and V h (s 1 ,ŝ 2 ) are finite, we obtain from (77) that: so that by eliminating Q h (s 1 ): We emphasize that the kernel K h (s 1 , s 2 ) in Equation (77) is the same as in Section 4. The analysis of the kernel in Section 4 can thus be applied in this section as well. For every s 2 ∈ Σ, defined by (53), there exists a unique, positive real s 1 such that K h (s 1 , s 2 ) = K h (s 1 ,s 2 ) = 0. Since for every s 2 ∈ Σ it holds that Re[ bounded whenever s 2 ∈ Σ and when s 1 is given by (51). Then, Equation (80) implies that fors 2 = s,ŝ 2 =s, s ∈ Σ, (with x = Re[s], y = Im[s]): where s 1 is given by (51). Substituting (51) and (49) into the above gives us: Since Q h (s) = Q h (s), we have obtained that: Let f 0 again be a conformal map from the unit disk to interior region of Σ, such that f 0 is symmetric with respect to the real axis. We can thus write that: The solution of this boundary value problem is given by the Schwarz integral formula (Theorem 7.38 in [30]), yielding: where C 0 is the positively oriented unit circle, D a normalization constant constant, and Finally, let f = f −1 0 be the inverse mapping of f 0 . Then Q h (s) is given by: for all s in the interior of Σ. The function f is given by (62). We emphasize that s = 0 is located in this region, hence we can obtain integral formulas from (85) for all derivatives of Q h (s) evaluated at s = 0. The joint LST V h (s 1 , s 2 ) is fully determined by substituting (85) into (77).
In the remainder of this section, we will rewrite the contour integral in (85) as a real integral that is suitable for numerical integration.
Note that the first term of the right-hand-side is an even function in t ∈ R, while the second term is an odd function in t ∈ R. Therefore the first term does not contribute in the integral and we obtain that: We have thus obtained that Q h (s) is given by: The constant D can be obtained from the normalization condition Q h (0) = 1, yielding: where we have used (64).

Calculation of Moments
In this section we will compute the mean and variance of the number of type-j customers (j = 1, 2) and the mean and variance of the total number of customers, as λ ↑ 1 2 . From these performance measures, the covariance and correlation between the number of type-1 and type-2 customers can be deduced.

Conclusions
We compared two similar, albeit slightly different, discrete-time two-class queueing models that fall within a class of queueing models that are known to be hard to analyze. In this paper, we have combined two well-known methods, namely the heavy-traffic limit method and the boundary value method. By combining these two approaches, we have succeeded in deriving easy-to-use formulas for the correlation coefficient between the number of type-1 and type-2 customers, when the queues are close to saturation. We emphasize that we have assumed a symmetrical arrival process, but that it can have any distribution. The results reveal that the two queueing models are in fact very different from each other in heavy traffic, which goes against our prior intuition.