Asymptotic behavior of a system of two coupled queues when the content of one queue is very high

We consider a system of two parallel discrete-time single-server queues, queue 1 and queue 2. The service time of any customer in either queue is equal to 1 time slot. Arrivals during consecutive slots occur independently from slot to slot. However, the arrival streams into both queues are possibly mutually interdependent, i.e., during any slot, the numbers of arrivals in queue 1 and queue 2 need not be statistically independent. Their joint probability generating function (pgf) A(x, y) fully characterizes the queueing model. As a consequence of the possible intra-slot correlation in the arrival process, the numbers of customers present (“system contents”) in queues 1 and 2, at any given slot boundary, are not necessarily independent either. In a previous paper, we have already discussed the mathematical difficulty of computing their steady-state joint pgf U(z1,z2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U(z_1,z_2)$$\end{document}; explicit closed-form results can only be obtained for specific choices of A(x, y). In this paper, we therefore look at the problem from an other angle. Specifically, we study the (asymptotic) conditional steady-state behavior of the system under the condition that the content of queue 1 is (temporarily) very high (goes to infinity). For ease of terminology, we refer to the system as the “asymptotic system” in these circumstances. We prove that the asymptotic system is nearly identical to the original (unconditional) system, but with a modified joint arrival pgf A∗(x,y)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^*(x,y)$$\end{document} that can be computed explicitly from A(x, y). This fundamental result allows us to determine the stability condition of queue 2 in the asymptotic system, and explicitly compute the classical queueing performance metrics of queue 2, such as the pgf, the moments and the approximate tail distribution of its system content, when this condition is fulfilled. It also leads to accurate approximative closed-form expressions for the joint tail distribution of the system contents in both queues, in the original (unconditional) system. We extensively illustrate our methodology by means of various specific (popular) choices of A(x, y). In some examples, where an explicit solution for U(z1,z2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U(z_1,z_2)$$\end{document} or for the (approximative) joint tail distribution is known, we retrieve the known results easily. In other cases, new results are found for arrival pgfs A(x, y) for which no explicit results were known until now.


Introduction
The purpose of this paper is to contribute to a subdomain of queueing theory concerned with the joint analysis of multiple queues.Whereas a multitude of papers exists in which mathematical techniques are developed for the analysis of one single, isolated queue, analyses where several queues are considered simultaneously are rather scarce, unless, of course, these queues operate completely independently of each other.A joint analysis of multiple queues becomes (nearly) unavoidable as soon as these queues influence each other; such queues are sometimes referred to as "coupled queues".There may be various reasons for the coupling between queues: the arrival streams into the queues may be mutually dependent, the queues may have to share the same service facilities, (part of) the input into one queue may be (part of) the output of other queues, etc.It is well-known that the joint analysis of coupled queues is, in general, quite hard, even when restricted to the situation where only two queues are at play.
Various mathematical techniques have been proposed to study the behavior of two coupled queues, where the purpose of the analysis basically comes down to an attempt to compute the steady-state joint distribution or joint probability generating function (pgf) of the numbers of customers present (the "system contents") in these queues, or, at least, some (mixed) moments, tail approximations, etc. of the joint distribution.With no claim on completeness, we briefly summarize some of those techniques here.
In many cases, a kernel-type functional equation for the joint pgf of the two system contents pops up in the analysis of two coupled queues.A particularly well-known general principal method, from the field of complex analysis, to find exact solutions of such equations is the so-called "boundary-value approach", which is described in great detail in the classical books by Cohen and Boxma [12] and by Fayolle et al. [17].Although the boundary-value approach can deal with various kinds of kernel-type functional equations, it has the disadvantage that it involves singular integrals and conformal mappings, which may be very complicated, and also requires quite some additional numerical work.It does not usually lead to "simple", easy-to-use formulas.The theory behind the boundary-value approach was developed several decades ago; yet the technique has been used frequently ever since; some examples of (relatively recent) applications are in the context of tandem queues (see, e.g., Resing and Ormeci [26], van Leeuwaarden and Resing [28]), polling systems (see, e.g., Saxena et al. [27]) and retrial queues (see, e.g., various papers by Dimitriou [14][15][16]), and there are many more.
Another method that has been used on various occasions to study the behavior of two coupled queues is the "compensation approach", a rather versatile technique for the analysis of two-dimensional Markov chains satisfying certain conditions, without transforms (pgfs), initiated by Adan and various co-authors, and very well documented in the paper Adan et al. [4].Basically, in this method, the desired unknown joint distribution is expressed as a sum of bivariate geometric product forms satisfying the inner balance equations of the Markov chain, and the coefficients of the individual terms in the sum are determined in a clever way.The compensation approach has been used in various applications and contexts; in particular, we mention the study of switching elements in communications networks (Boxma and van Houtum [6]) and the analysis of multiprogramming queues (Adan et al. [2]); a more general application of the compensation approach for random walks in the quarter plane is discussed in Adan, van Leeuwaarden and Raschel [3].
Other researchers have successfully applied the so-called "uniformization method", a complex-function-based technique, for various models; see, e.g., Kingman [24] and Flatto and McKean [19] for the study of two parallel queues and Jaffe [22,23] for the analysis of clocked switching elements.In the uniformization method, a parameter representation is introduced for the algebraic curve that represents the set of zerotuples of the kernel, the unknown boundary functions in the functional equation are also expressed in terms of this parameter, referred to as the "uniformization variable", and equations are derived and subsequently solved for these boundary functions.In our view, the concept of analytic (or meromorphic) continuation plays a more important role in these analyses than the "uniformization" itself, which does not really reflect the essence of the technique, but the terminology seems, by now, widely accepted.Often, analytic continuation can also be applied directly to the functional equation, without the need of parametrization (or "uniformization"); see, e.g., Cohen [10,11] for the analysis of the so-called clocked buffered switch, and Devos et al. [13] in the study of a two-class randomly alternating service model.
We also mention the so-called "power-series technique", mainly devised by Hooghiemstra, Keane and Van de Ree [20] and Blanc (see, e.g., [5]), and further developed by, among others, Hooghiemstra and Koole [21], as a useful method to approximately analyze the steady-state joint distribution of some two-queue systems.This technique consists of expanding the desired distribution (or the associated pgf) as a power series in terms of some system parameter (usually, the load in some queue), and then finding the coefficients of the consecutive powers in that power series as the solutions of "easier" equations than the original functional equation.An adapted version of the power-series technique was successfully applied by Walraevens et al. [30] for the approximate analysis of generalized processor sharing (GPS) queues.
Of course, there are other (ad-hoc) methods as well.An excellent overview of methods that have been used to find the joint distribution of the system contents in two-queue systems can be found in the paper by Adan et al. [1].
In the current paper, we look at the problem of simultaneously analyzing two coupled queues from yet another perspective.Instead of trying to solve the associated balance equations or functional equation, we focus on the asymptotic conditional behavior of the two-queue system when one queue is (temporarily) overloaded, i.e., when the system content of one queue is (temporarily) very high (goes to infinity).
We do this for a very simple example system of two coupled parallel queues, where the coupling is caused solely by the fact that the two queues receive external arrival streams that are possibly mutually interdependent.In order to keep the discussion as "pure" as possible, we intentionally allow no other sources of coupling between the two queues; in particular, we assume that arrivals in either queue are not dependent on the states of the queues (as could be the case, for instance, when customers, upon arrival, always join the shortest queue, or when arrivals in one queue are departures from the other queue) and that both queues have their own dedicated server which is permanently available for the service of customers of that queue only.Systems of two parallel queues have been investigated frequently in the past, as illustrated by the numerous references mentioned earlier.We emphasize in particular the early (seminal) works by Kingman [24] and Flatto and McKean [19] who consider continuous-time models for the join-the-shorter-line scheduling scheme, and Konheim et al. [25] who consider one shared server for both queues, also through a continuous-time model.But the specific system considered in the current paper, i.e., two parallel queues with interdependent arrivals (and no other sources of coupling) has been dealt with in just a few papers, mainly in a discrete-time setting, especially in the context of the analysis of a clocked buffered switch, in papers by Jaffe [22,23], Cohen [10,11], Boxma and van Houtum [6] and Adan et al. [1], as mentioned earlier.The paper Bruneel [7] also falls within this category, albeit that the joint arrival distribution can be more general than "binomial of order two" (as is the case in the clocked buffered switch) there.In the current paper, we essentially investigate the same general model as Bruneel [7], but from the perspective of analyzing it under the specific condition that one of both queues contains a large number of customers.The approach taken in this paper does not require any specific assumptions with respect to the joint arrival distribution, and, to the best of our knowledge, is new.
The rest of the paper is organized as follows.Section 2 describes the full details of the mathematical model we consider.In Sect.3, we establish the basic functional equation for the steady-state joint pgf of the two system contents, and derive explicit expressions for the two marginal pgfs (and mean values) of the system contents in the individual queues.In Sect.4, we define the "asymptotic system" as a convenient short terminology for the original two-queue system, considered under the explicit condition that the system content in queue 1 is (temporarily) very large (goes to infinity).We then formally prove that the asymptotic system is nearly identical to the original (unconditional) system, i.e., it has two parallel discrete-time single-server queues, deterministic 1-slot service times, independent arrivals from slot to slot, but a modified joint arrival distribution, that can be computed explicitly from the original joint arrival distribution, regardless of its nature.We then apply this fundamental result to examine the stability condition and the steady-state queueing performance of queue 2 in the asymptotic system; we also derive accurate approximative closed-form expressions for the joint tail distribution of the two system contents in the original (unconditional) system.Section 5 illustrates our methodology extensively by considering a number of (well-known) example arrival distributions in great detail.Here, we easily retrieve previously known results, when these are available, but also various new results for systems that have not been solved so far.In Sect.6, we compare the performance of queue 2 in the asymptotic system with its performance in the original system, in terms Fig. 1 System of two parallel single-server queues with joint arrival pgf A(x, y).A 1 (x) and A 2 (y) denote the marginal arrival pgfs in queue 1 and queue 2, respectively of the arrival rate, the mean system content and the (approximate) tail distribution of the system content, and try to explain the observed differences in terms of the nature of the intra-slot correlation in the arrival process.Some conclusions are stated in Sect.7. The paper also has two appendices; "Appendix A" simply proves a mathematical lemma, required in Sect.4, but "Appendix B" is important in its own right and shows how to solve the functional equation explicitly for a very specific joint arrival distribution that was also considered in one of the examples of Sect.5; we had this interesting (non-trivial) result for a number of years but it was never published before.

Mathematical model
In this paper, we investigate a discrete-time queueing system, consisting of two parallel queues, named queue 1 and queue 2, with infinite waiting room and one server each.A graphical representation of the system is shown in Fig. 1.
As in all discrete-time models, the time axis is divided into fixed-length intervals referred to as (time) slots.New customers may enter the system at any given (continuous) point on the time axis, but services are synchronized to (i.e., can only start and end at) slot boundaries.Customers arriving to queue 1 and to queue 2 are referred to as type-1 customers and type-2 customers, respectively.We assume that the service of each customer requires exactly one time slot, regardless of whether the customer is of type 1 or type 2.
The arrival process of new customers in the system is characterized by means of a sequence of i.i.d.nonnegative discrete random vectors (a 1,k , a 2,k ) with common joint probability mass function (pmf) a(i, j) and common joint probability generating function (pgf) A(x, y), respectively.Here a 1,k and a 2,k indicate the numbers of type-1 and type-2 arrivals, respectively, during slot k, and, hence, a(i, j) and A(x, y) are defined as which are independent of k.The (marginal) pgfs of a 1,k and a 2,k are given by respectively.The mean number of arrivals of type i per slot is denoted as Note that, just as the joint pgf A(x, y), the marginal pgfs A 1 (x) and A 2 (y) and the mean arrival rates λ i do not depend on the slot index k.
3 Steady-state queueing analysis

Stability conditions of the system
The main purpose of our paper is to analyze some aspects of the steady-state behavior of the queueing system under study.It is not difficult to see that the system is stable, i.e., a steady state exists, if and only if both queue 1 and queue 2 are stable.As both queues are simple single-server queues with deterministic 1-slot service times, the stability conditions thus read

System evolution as a function of time
Let u 1,k and u 2,k indicate the system contents, i.e., the numbers of customers present in queue 1 and queue 2, respectively, at the beginning of slot k.We indicate their joint pgf as We note that, in this paper, the term system content always refers to the total number of customers in the queue, including the customer in service, if any.As each queue of the system has exactly one server that can serve exactly one customer per time slot, the following system equations are valid between the random variables u 1,k and u 2,k on the one hand, and the random variables u 1,k+1 and u 2,k+1 , one slot later, on the other hand: Here we have introduced the notation (x) + to indicate the quantity max(0, x).Equations (6) can be translated into one corresponding equation between the joint pgfs U k (z 1 , z 2 ) and U k+1 (z 1 , z 2 ), by using definition (5), which results in

Steady state
We now let the time parameter k go to infinity.Assuming the system reaches a steady state, i.e., assuming the stability conditions (4) are met, then both functions U k (•) and U k+1 (•) converge to a common limit function U (•), which denotes the joint pgf of the two system contents at the beginning of an arbitrary slot in the steady state.As a result, Eq. ( 7) translates into a linear equation for U (z 1 , z 2 ), which can be written as where the unknown function L(z 1 , z 2 ) is defined as and the "kernel" As has been discussed in e.g.Bruneel [7], the functional equation ( 8) is very hard to solve analytically, for general arrival pgfs A(z 1 , z 2 ).Only for specific classes of arrival pgfs, the explicit determination of the boundary functions U (z 1 , 0) and U (0, z 2 ), and, from this, the unknown function L(z 1 , z 2 ), and the whole pgf U (z 1 , z 2 ), in closed analytic form, seems feasible.We will discuss some rare examples in Sect. 5 (notably in Sects.5.1 and 5.2.2).The marginal pgf U i (z) of the system content in queue i, however, can be easily deducted from (8) by choosing either {z 1 = z, z 2 = 1} or {z 1 = 1, z 2 = z}, which results in where the marginal pgf A i (z) and mean value λ i of the number of type-i arrivals per slot were defined in Eqs. ( 2) and (3).The mean system content in queue i is given by These simple formulas for the individual queues are not surprising, as queue i is a single-server queue with 1-slot service times and arrival pgf A i (z).These formulas can be found in many papers and books on discrete-time queueing theory; see, e.g., chapter 4 of Bruneel and Kim [8].

Asymptotic behavior
From now on, let us assume that the system has reached a steady state.We then define the random variables u 1 and u 2 as the system contents in queues 1 and 2, respectively, at the beginning of an arbitrary slot in the steady state.We now shift our focus from the analysis of the full joint distribution (or pgf) of the two system contents to the study of the conditional asymptotic behavior of one queue, say queue 2, given that the system content in the other queue, i.e., in queue 1, is very high.Specifically, we are interested in determining such quantities as the conditional pgf, mean value, approximate tail distribution,... of u 2 , given that u 1 = n, for large values of n.In order to do so, we first make some very mild assumptions on the marginal arrival pgfs A 1 (z) and A 2 (z), as stated in the "arrival assumption" below.

Arrival assumption
The pgf A i (z) has a radius of convergence i > 1 and lim z→ i A i (z) = +∞.
Pgfs that satisfy the above condition are sometimes called "well-behaved".The arrival assumption implies, among other things, that A 1 (z) and A 2 (z) do not have singularities on the unit circle in the complex z-plane, because their radius of convergence is strictly larger than 1.Note that most commonly used discrete arrival distributions have pgfs that satisfy these conditions.

Tail distribution of the system content in one queue
It is well known (see, e.g., Flajolet and Odlyzko [18], Bruneel and Kim [8], Bruneel et al. [9], Van Mieghem [29]) that the tail distribution of a discrete random variable can be very well approximated based on the so-called "dominant singularity" of the pgf of that random variable.Specifically, if the pgf X (z) of a discrete random variable X has a radius of convergence τ > 1, then τ (a positive real number) is a singularity of X (z) and X (z) cannot possess singularities with a smaller absolute value than τ , i.e., τ is a "dominant" singularity of X (z).Moreover, in the very common case where X (z) has no other (complex) singularities with the same modulus as τ (i.e., τ is the only dominant singularity of X (z)), and τ is a pole of multiplicity 1, the tail distribution of X can be very well approximated by a geometric form with decay rate τ −1 : for large values of n.The constant c is the negative residue of the pgf X (z) at z = τ , i.e., This approximation is usually referred to as the "dominant-pole approximation" of the tail distribution of X .The terminology "dominant" stems from the fact that the exact value of the inverse z-transform of X (z), i.e., the pmf Prob[X = n], can be expressed as a sum of (geometric or similar) contributions of all the singularities of X (z), each with a decay rate inversely proportional to the singularity at hand; for large n, the contribution of the singularity with smallest modulus therefore "dominates" all the others, and the sum can be well approximated by the contribution of the dominant singularity alone.We now apply the above approximation technique to the marginal pgf U i (z) of the queue-i system content.In order to do so, we first notice that the singularities of U i (z), as given by Eq. ( 9), are the zeroes of the expression z − A i (z).We show in "Appendix A" that, if A i (z) is "well-behaved" with radius of convergence i > 1, as defined above in the "arrival assumption", then z − A i (z) has exactly one real positive zero, say τ i , with 1 < τ i < i , which has multiplicity 1; obviously τ i satisfies the equation furthermore, in "Appendix A", we also show that and that z − A i (z) has no other zeroes outside the unit disk in the complex z-plane with modulus smaller than τ i .This means that τ i is the dominant pole of the pgf U i (z) and an approximation such as (11) can be used.As a result, we obtain where, according to ( 9) and (12), or, upon use of de l'Hospital's rule and Eq. ( 13), In Fig. 2, we compare the dominant-pole tail approximation (15), (16) for Prob[u 1 = n] with simulation results for three different well-known arrival distributions with mean arrival rate λ 1 : 2 (binomial arrivals of order two), with λ 1 = 0.9.
The graphs show that the approximation is extremely accurate, even for small values of n.In Sect.5.2.1, we shall see that, for geometric arrivals, u 1 happens to be geometrically distributed and the "tail approximation" is even exact for all values of n; see Eq. (60).

Conditional arrival process
Our main goal being the study of the conditional behavior of the system content in queue 2, given that the system content in queue 1 is (temporarily) very large, it would seem natural to focus immediately on the (existence and-if so-the computation of the) conditional pgf of the queue-2 system content u 2 given the (high) value (n) of u 1 at the beginning of the same (arbitrary) slot in the steady state, i.e., the conditional pgf and its limit as n → ∞, which we denote as Here, however, we take a different approach.As the system content in queue 2 at the beginning of any slot in the steady state, say slot S, is merely the result of the arrivals and the services in queue 2 before slot S, and the services are, in fact, very deterministic, we first concentrate on the conditional arrival process in queue 2 before slot S, given that the system content in queue 1 at the start of slot S is equal to n, for large values of n.It turns out that it is not more difficult to analyze the conditional joint arrival process in both queues 1 and 2 than in just queue 2.Moreover, we shall see that this yields some interesting additional information.Our approach is based on this observation.Formally, we proceed as follows.Let S be an arbitrary slot in the steady state.Let u 1,S denote the system content at the beginning of slot S in queue 1.As S is an arbitrary slot, u 1,S is distributed according to the known pgf U 1 (z) defined in (9).Let a 1,S− and a 2,S− ( ≥ 1) denote the numbers of arriving customers of type 1 and 2, respectively, slots earlier than slot S, i.e., during slot S− .For any integer m ≥ 1, the (unconditional) joint distribution of the random variables {a 1,S− , a 2,S− , 1 ≤ ≤ m} is well-known; as the arrivals in our system are independent from slot to slot, the associated joint pgf is given by a simple product form, i.e., Note that this result is valid for any choice of the slot S, and therefore the notation (B m ) has no dependence on S.
Next, assume that, for the chosen slot S, u 1,S is equal to a given value n; if desired, the probability that this is the case, i.e., can be computed from the known pgf U 1 (z) (see Eq. ( 9)) by inverse z-transformation.
Of course, if it is known that u 1,S = n, this has an impact on the possible values of the numbers of arrivals in the slots before S (and their probabilities).In general, it is to be expected that, under this condition, the random variables {a 1,S− , a 2,S− , 1 ≤ ≤ m} have a different joint pgf than the above defined B m , and are even not necessarily mutually independent anymore from slot to slot.We are therefore interested in the computation of the conditional joint pgf of the random variables {a 1,S− , a 2,S− , 1 ≤ ≤ m} given that u 1,S = n.Formally, we define this conditional joint pgf as where we may again suppress the dependence on S in the notation.We now go through a series of consecutive steps to derive an explicit expression for the function F m|n , for large values of n.

Step 1: the joint pgf V m
We start by defining the joint pgf of {a 1,S− , a 2,S− , 1 ≤ ≤ m} and u 1,S as This definition can be easily extended to include the boundary case m = 0, by agreeing that an empty product is equal to 1, i.e., where we have used (9).Note that, although the function V m , as defined in (21), does not explicitly depend on the choice of the slot S, we will (temporarily) add a subscript S in our notation (i.e., V m,S = V m ), in order to improve the clarity of the derivations that follow.

Step 2: a recurrence relation for V m
Using system Eq.( 6) for k = S − 1 in ( 21), we can express V m+1,S in terms of V m,S−1 : Since the arrivals during slot S − 1 are independent of all the earlier arrivals (and, hence, also of u 1,S−1 ), the above expected value can be expressed as a product of two expected values, i.e., Here the second factor can be expressed in terms of the function V m,S−1 , and so Suppressing the dependence on S (or S−1) again, we obtain the following (recurrence) relation between V m+1 and V m :

Step 3: z = 1 is the dominant pole of V m
For m = 0, (23) leads to or in view of ( 22), Similarly, for m = 1, we get from ( 23) and ( 25) Analogously, we could derive an explicit expression for the function V m , for any finite value of m, from the recurrence relation (23).However, such an expression would be complicated (for general m), and, more importantly, is not really needed for the purpose of computing the function F m|n , defined in (20), which is our true goal.The expressions ( 22), ( 25) and ( 26) for the functions V 0 , V 1 , V 2 do make clear that these functions, considered as functions of the variable z, all have a single pole at z = τ 1 , which is a zero of their common denominator (z − A 1 (z)).It is readily seen that this statement holds true for all the functions V m , m ≥ 0.Moreover, the pole at z = τ 1 is the dominant pole of the functions V m if -as usual -we restrict the values of the other variables (x and y ) to the closed unit disks in their respective complex planes, i.e., |x | ≤ 1, |y | ≤ 1.One easy intuitive way to see this is by observing that factors such as A(x z, y ), appearing in the numerators of the formulas for V m remain finite when |z| < τ 1 : which is certainly finite for |z| < τ 1 , since the radius of convergence 1 of the marginal arrival pgf A 1 (z) is larger than τ 1 , as discussed in Sect.4.1 and "Appendix A".

Step 4: deriving F m|n from V m
From the definitions ( 20) and ( 21), it easily follows (by an application of the law of total expectation) that Hence, the product u 1 (n)F m|n is nothing but the coefficient of the power z n in the power series expansion of the function V m , considered as a function of the variable z, i.e., for given values of the variables {(x , y ), 1 ≤ ≤ m}.In principle, this product can thus be retrieved from V m through inverse z-transformation.Here, we are interested in large values of n, and, therefore, the inverse z-transform of V m is dominated by the term associated to its dominant pole, which is z = τ 1 .For large n, we can therefore use the following accurate dominant-pole approximation: and, hence, where W m denotes the negative residue of V m at z = τ 1 , i.e.,

Step 5: determining W m
Multiplying both sides of the recurrence relation ( 23) by a factor −(z − τ 1 ) and then taking the limit for z → τ 1 , we easily get Equation ( 33) is a (very simple) recursive equation for the functions W m ; applying (33) repeatedly for consecutive values of m, we can express W m as The function W 1 can be easily derived from (25) as which, eventually, leads to (36)

Step 6: approximate expression for F m|n
We have shown in Sect.4.1 (Eqs.( 15) and ( 16)) that, for large n, the probability u 1 (n) can be well approximated as Using ( 36) and ( 37) in (31), we find that, for large values of n, the function F m|n can be accurately approximated as which turns out to be independent of n.Introducing the notation we can thus conclude that, for large values of n, As a consequence of the properties of dominant-pole approximations, as used in the above derivation (in both the numerator and the denominator of (31)), the above approximative equality becomes a true equality when n → ∞.An easy way to see that this statement is true, without complicating the notations too much, is as follows.
First, replace the approximative equality ( 30) by an exact equality by introducing a correction term o 1 (τ −n−1

1
) that summarizes the contributions of all the non-dominant singularities in the inverse z-transform of the function V m : where, for notational convenience we have (temporarily) omitted the arguments of the functions F m|n and W m , and ) is a function of n that goes to zero more rapidly than ). (42) From ( 41) and (42), we then get the following exact expression for F m|n : Taking the limit for n → ∞ in (43), and in view of the properties of the auxiliary functions o 1 and o 2 , we then easily get Finally, restoring the full notations again, and using (36), we obtain the following exact result: We note that the "asymptotic" conditional joint pgf F m has a similar product form as the unconditional joint pgf B m , defined in (18), the only difference being that the unconditional joint pgf A(x, y) of the numbers of type-1 and type-2 arrivals in one slot has been replaced by the "asymptotic" conditional joint pgf A * (x, y), defined in (39).

The "asymptotic" system
In this paper, we basically study a system of two parallel single-server queues with 1slot service times and independent arrivals from slot to slot, characterized by the joint pgf A(x, y), as defined in (1).Intuitively speaking, we have shown in the previous section that, conditional on the fact that the system content u 1 = n and n is large at the beginning of an arbitrary slot in the steady state, the arrivals into the twoqueue system during the preceding slots are (approximately) independent from slot to slot, and distributed according to the joint pgf A * (x, y), as defined in (39); the approximation becomes exact when n → ∞.Let us define the "asymptotic system" as the (conditional) "limit version" of the (original) system where the condition "u 1 → ∞" is imposed.We have shown in the previous section that the asymptotic system is almost identical to the original system, the only difference being that the arrivals are characterized by the joint pgf A * (x, y) instead of A(x, y).We emphasize that the asymptotic system is not a new, different system than the original system; it is just a convenient short terminology to refer to the original system, whenever the latter is observed under the condition of u 1 approaching infinity.In the asymptotic system, according to Eq. (39), the marginal pgfs of the numbers of arrivals (per slot) of type 1 and type 2 are given by The mean arrival rates (per slot) in queues 1 and 2 are According to (14), A 1 (τ 1 ) > 1, and we can thus conclude that i.e., in the asymptotic system, the arrival rate in queue 1 is strictly greater than 1.This matches perfectly with the defining condition of the asymptotic system: u 1 → ∞, i.e., queue 1 builds up continuously in the asymptotic system.On the other hand, the stability of queue 2 in the asymptotic system depends on the value of λ * 2 .Specifically, in the asymptotic system, the second queue is stable if and only if The correct interpretation of these results deserves some discussion.Recall that we are studying the original two-queue system under the assumption that it has reached a steady state, globally speaking.To that end, we have assumed that both λ 1 and λ 2 are strictly smaller than 1, i.e., on average both queue 1 and queue 2 receive less than 1 customer per slot.However, this does not imply that the (conditional) mean arrival rate in either queue cannot be temporarily greater than 1.As far as queue 1 is concerned, according to (48), this is exactly what happens in the asymptotic system.We assume there that, at some moment in time, the system content u 1 is very large (theoretically: goes to infinity).Of course, this can only happen if the arrivals during (infinitely) many preceding slots have (on average) exceeded the service capacity of the queue.We emphasize that the probability of u 1 being very large may not be very high in a system that is globally stable, but yet the event is possible (as illustrated in equation ( 15)) and-in the asymptotic system-we study exactly this case.Our results reveal that, in the asymptotic system, i.e., the original system conditioned on the event that, at some time instant in the steady state, u 1 → ∞, the (conditional) mean arrival rate λ * 1 in queue 1 (in the preceding slots) is indeed larger than 1, as could be expected.More interestingly, they also show that, if u 1 takes very large values, then queue 2 may be either stable or instable, depending on λ * 2 being smaller or greater than 1.
The essential reason why the system content in queue 1 has an impact on the arrivals in queue 2 is, of course, that the arrival streams in both queues are possibly nonindependent (within a slot).It is therefore to be expected that, in the special case where the arrivals in queues 1 and 2 are strictly mutually independent, this phenomenon disappears.Indeed, assume that then, according to (46) and ( 13), i.e., the arrivals in queue 2 in the asymptotic system are identical to their counterparts in the original, unconditioned system.

Conditional system content in queue 2
We now take a closer look at the behavior of queue 2 in the asymptotic system.From our study of the conditional arrival process in Sect.4.2, it follows that queue 2 behaves as a regular discrete-time single-server system with 1-slot service times, independent arrivals from slot to slot, and arrival pgf A * 2 (y), as defined in Eq. ( 46).Queue 2 is stable if and only if λ * 2 , as defined in (48), is smaller than 1.If this is the case, then the steady-state pgf of the system content in queue 2, i.e., the pgf U * 2 (z), defined in (17), exists and is given by a formula similar to the ones in Eq. ( 9), but with a modified arrival process: Various (conditional) moments of u 2 , given that u 1 → ∞, can be easily computed from (49), in the usual way.For instance, the mean value can be obtained from It is also possible to derive an accurate approximation for the (conditional) tail distribution of u 2 , given that u 1 → ∞, by using a dominant-pole approximation as summarized in the Eqs.(11) and (12), with X (z) replaced by U * 2 (z), i.e., for large values of , where and where σ denotes the dominant pole of U * 2 (z), i.e., σ is the smallest real positive solution, larger than 1, of the equation The existence of a dominant pole σ for the pgf U * 2 (z) basically follows from the same kind of argument as the existence of the dominant pole τ i of the pgf U i (z), as discussed in "Appendix A", and the observation (in all the examples that we considered) that

Some specific examples
In this section, we explore some specific instances of the general two-queue model discussed in the previous sections.This amounts to considering specific choices for A(x, y).

Identical arrivals in both queues
Assume that the two queues receive exactly the same number of arrivals in each slot, i.e., for all values of the slot index k, a 1,k = a 2,k .In this case, A(x, y) can be expressed as a function of the product of its arguments x and y, i.e., A(x, y) = B(x y), where B(z) is a well-behaved pgf with radius of convergence > 1.The marginal arrival pgfs and mean values are It is very clear that, in the steady state (possibly after some finite transient period), the two queues contain exactly the same number of customers (at the beginning of each slot), and, consequently, the joint pgf U (z 1 , z 2 ) of the two system contents can also be expressed as a function of the product of its two arguments.Specifically, where U 1 (•) and U 2 (•) are the known functions defined in (9).We note that, according to Bruneel [7], this is one of the few known cases in which a fully explicit analytic expression for the joint pgf U (z 1 , z 2 ) can be derived, which is, of course, due to the very specific and simple nature of the system.From our formulas derived in Sect.4, we easily find that the dominant poles τ i (as defined in ( 13)) of the marginal system-content pgfs are τ 1 = τ 2 τ , where τ is the only (real positive) solution of the equation τ = B(τ ), 1 < τ < , and from (14), we also know that B (τ ) > 1.According to (39), in the asymptotic system, the joint arrival pgf is given by α C(z) Fig. 3 A global arrival stream with pgf C(z) is split probabilistically into two substreams with pgfs A 1 (x) and A 2 (y), where α and 1 − α indicate the fractions of customers routed to queues 1 and 2, respectively which is also a function of the product x • y.It then easily follows that the marginal arrival pgfs and mean values in the asymptotic system are As expected, we find that λ * 1 > 1, signifying that, in the (conditional) asymptotic system, queue 1 is unstable, in accordance with its defining condition u 1 → ∞.However, since in this particular case, λ * 1 = λ * 2 , we find that queue 2 is unstable too, i.e., if u 1 → ∞ then also u 2 → ∞, which is in complete agreement with the fact that the system contents in both queues, at the beginning of any slot in the steady state, are identical.In conclusion, our methodology clearly reveals that, in the case of identical arrivals in both queues, a steady-state pgf U * 2 (z), as defined in (49), does not exist.

One global arrival stream probabilistically routed to both queues
In our model, the total number of arrivals into the system (i.e., both queues together) during slot k is given by c k = a 1,k + a 2,k .It is easily seen that the c k s also form a sequence of i.i.d.nonnegative discrete random variables with common pgf C(z) A(z, z).We call C(z) the "pgf of the global arrival stream".In this subsection, we assume that a global arrival stream of customers with corresponding pgf C(z) presents itself at the entrance of our two-queue system.Consecutive customers are then routed probabilistically, i.e., independently from customer to customer, and with fixed probabilities, to queue 1 or queue 2, respectively.Let α and 1 − α denote the probabilities that an arriving customer is routed to queue 1 or queue 2, respectively.Figure 3 represents this setting graphically.
It is readily seen that, in this situation, the joint arrival pgf A(x, y) and the marginal arrival pgfs A 1 (x) and A 2 (y) are given by Let μ indicate the total mean arrival rate, i.e., μ C (1).Then, the mean arrival rates in the individual queues (see (3)) are and, of course, λ 1 + λ 2 = μ.The two-queue system is stable if and only if λ 1 < 1 and λ 2 < 1 (see ( 4)), or, in this setting, if and only if Note that this implies μ < 2.
It is not possible to obtain more detailed results on the behavior of the asymptotic system for general C(z).In the following subsections, we therefore consider some specific choices for C(z).

Global geometric arrivals
Assume a geometric distribution with mean μ for the global arrival stream, i.e., The joint and marginal arrival pgfs A(x, y), A 1 (x) and A 2 (y) can be obtained from (55) as Note the remarkable property that the marginal arrival distributions in queues 1 and 2 are also geometric, with means λ 1 = μα and λ 2 = μ(1 − α), respectively.The marginal pgf and mean value of the system content in queue i can be found from ( 9) and (10) as We thus observe that the system content in queue i is also geometrically distributed, notably with parameter λ i .The corresponding pmf is well-known: The dominant pole τ i is the (only) zero of the denominator of U i (z), as given in (59): According to (46), (58), the marginal arrival pgfs and mean values in the asymptotic system are and Some first conclusions can be drawn from these results.In the asymptotic system, the marginal arrival distributions are geometric, just as in the original system, but the mean arrival rates are λ * 1 and λ * 2 instead of λ 1 and λ 2 .From (62), we learn that λ * 1 is strictly greater than one, because it is the inverse of λ 1 < 1, which again illustrates the instability of queue 1 in the asymptotic system.On the other hand, (63) shows that the stability of queue 2 in the asymptotic system is determined solely by the value of the routing probability α, or, equivalently, the ratio of λ 1 and λ 2 .Queue 2 is stable if and only if This result is intuitively appealing: in order for queue 2 to be stable when queue 1 goes into saturation (i.e., u 1 → ∞), it should receive less customers per slot than queue 1.Moreover, the (conditional) arrival process in queue 2 does not depend on the total arrival rate μ of the original system, but is fully determined by the value of α, a most intriguing result.Now let us assume that λ 2 < λ 1 .Then, in the asymptotic system, queue 2 behaves as a stable single-server queue with 1-slot service times and geometric arrivals with mean λ * 2 , implying that the pgf and the mean value of the queue-2 system content exist and are given by The dominant pole σ of U * 2 (z), as defined in (53), is the (only) zero of its denominator, i.e., As discussed in Sect.4.4, the (conditional) tail distribution of u 2 , given that u 1 → ∞, can be well approximated by using the dominant-pole technique, i.e., for large , Hence, combining (66) and (67), we get, for large , The joint pmf of (u 1 , u 2 ) can be expressed, in general, as Here the first factor is given by the exact formula (60), valid for any n, and the second factor can be well approximated by the dominant-pole form (68), valid for large n and .As a result, we obtain a closed-form approximation for the joint tail distribution of (u 1 , u 2 ), i.e., for large values of both n and , The high accuracy of formula (70) for large n (say n ≥ 10), is illustrated by one specific example, in Fig. 4, where a comparison is made with simulation results.Note that we did not even select high values of here.We obtained a similar accuracy when keeping n fixed and varying .Producing simulation results for high n and at the same time, is extremely time-consuming, since it concerns very rare events; evaluating the closed formula (70), on the contrary, requires virtually no time at all.
An additional result can be obtained by exploiting the symmetry of our model.Indeed, we can interchange the roles of queue 1 and queue 2, in the sense that we define the asymptotic system by the condition "u 2 → ∞" instead of "u 1 → ∞", and then study the behavior of queue 1 under this condition.This implies that the roles of the routing probabilities α and 1 − α must also be interchanged.As a result, we find, We note that, to the best of our knowledge, a fully explicit analytic expression for the whole joint pgf U (z 1 , z 2 ) has never been derived for the specific joint arrival distribution considered here, i.e., for the joint arrival pgf A(x, y) defined in Eq. (58).We strongly suspect that the fact that the approximative expressions (70) and (71) for the joint tail distribution of (u 1 , u 2 ) are profoundly different depending on the routing probability α being greater or smaller than 1  2 implies that such an analytic expression of U (z 1 , z 2 ), valid for all possible values of α, may not exist.We also note that, if α is exactly equal to 1  2 , then queue 2 is unstable when u 1 → ∞ and queue 1 is unstable when u 2 → ∞, and we cannot find a joint tail approximation at all.

Global shifted geometric arrivals
Assume a shifted geometric distribution with mean μ ≥ 1 for the global arrival stream, i.e., The joint and marginal arrival pgfs A(x, y), A 1 (x) and A 2 (y) can be obtained from (55) as .
The marginal pgfs of the system contents in queues 1 and 2 can be found from (9) as The mean system content E[u 2 ] = U 2 (1) follows from (73) as The dominant pole τ i is the (only) zero of the denominator of U i (z); from (73) we get According to (46), the marginal arrival pgfs in the asymptotic system are The corresponding mean arrival rates (per slot) in queues 1 and 2 are Again, just as in the case of (unshifted) geometric arrivals, it is clear that λ * 1 is strictly greater than one, because it is the inverse of λ 1 < 1, reflecting the instability of queue 1 in the asymptotic system.On the other hand, in this case, λ * 2 turns out to be always strictly smaller than one, regardless of the value of the routing probability α.This follows easily from the stability condition (57).This means that queue 2 is always stable when u 1 → ∞, in this case.Hence, the asymptotic pgf and mean value of the system content in queue 2 always exist and can be found from (49), ( 76), ( 77) and (50): It is remarkable that U 2 (z) and U * 2 (z) have the same denominator (upon a multiplicative constant).The dominant pole As discussed in Sect.4.4, the (conditional) tail distribution of u 2 , given that u 1 → ∞, can be well approximated, for large values of , as Hence, combining (79) and (80), we get Again, this can be combined with a tail approximation for the system-content distribution in queue 1, as per Eqs.( 15) and ( 16), i.e., for large values of n, where Hence Combining ( 81) and ( 84), we find that, for large values of both n and , We note that, in this case, we have found an approximation which is valid for all values of the parameter α.Moreover, formula (85) is invariant to the substitutions α ↔ 1 − α, n ↔ .In part, this is due to the remarkable property that the dominant poles σ and τ 2 are identical here, i.e., the decay rate of the tail distribution of u 2 is not influenced by u 1 , and vice versa.The case of global shifted geometric arrivals is also one of the very few instances of our model where an explicit analytic expression for the whole joint pgf U (z 1 , z 2 ) of the two system contents can be found, as opposed, for instance, to the seemingly very similar case of (unshifted) geometric arrivals, where no solution has been reported so far.Indeed, using an adapted version of a technique briefly discussed in Bruneel [7], we have been able to derive such an expression.The details of this derivation-which is interesting in its own right-are discussed in "Appendix B"; here we merely state the final result: where τ 1 and τ 2 are the dominant poles of U 1 (z) and U 2 (z), respectively, given in (75).Note that, apparently, τ 1 and τ 2 are also the only singularities of the joint pgf Having an explicit expression for U (z 1 , z 2 ) at our disposal, gives a unique opportunity to check the correctness of the asymptotic results obtained through the methodology developed in this paper.First, consider the series expansion of U (z 1 , z 2 ) in powers of z 1 , for given z 2 : If we consider U (z 1 , z 2 ) as a function of the variable z 1 , then (86) shows that U (z 1 , z 2 ) has a single dominant pole at z 1 = τ 1 , and hence, for large n, S n (z 2 ) can be well approximated as where c S (z 2 ) is given by The conditional marginal pgf of the system content in queue 2, given that u 1 = n, can be derived from this through For large n, according to (87), ( 83), ( 82) and ( 88), this can be well approximated by This is exactly the same expression as we obtained earlier in Eq. ( 78).
Next, we can also invert the known pgf U (z 1 , z 2 ) completely, i.e., compute the joint probabilities Prob[u 1 = n, u 2 = ] such that Departing from expression (86), we easily get Identifying equal powers of z 1 and z 2 in ( 89) and ( 90), we can compute Prob[u 1 = n, u 2 = ] explicitly, for any n, ≥ 0. It is not difficult to see that this leads precisely to expression (85), as soon as n ≥ 1, ≥ 1, which illustrates, once again, the accuracy, and, in this case even the exact nature of the asymptotic dominant-pole results derived via the methodology developed in this paper.

Global binomial arrivals of order two
Assume a binomial distribution of order two with mean μ for the global arrival stream, i.e., The mean arrival rates are λ 1 = μα and λ 2 = μ(1 − α) (see ( 56)) and the joint and marginal arrival pgfs can be expressed in terms of these as (see (55)) Here, again, we observe the remarkable property that the marginal arrival processes in queues 1 and 2 are also binomial of order two, with means λ 1 and λ 2 , respectively.The marginal pgf and mean value of the system content in queue i are (see (9)) The dominant pole τ i is the (only) zero of the denominator of U i (z); from (92), we get According to (46), the marginal arrival pgfs and mean arrival rates in the asymptotic system are Some conclusions can be drawn from these expressions.First, we observe that the two marginal arrival streams in the asymptotic system are also of type "binomial of order two", but with different mean arrival rates, i.e., λ * i instead of λ i .Again, as expected, it is clear that λ * 1 is strictly greater than one, because λ 1 < 1.And, just as in the case of global shifted geometric arrivals, λ * 2 turns out to be always strictly smaller than one, regardless of the value of the routing probability α, because λ 1 < 1 and λ 2 < 1.So, queue 2 is always stable when u 1 → ∞.Hence, the asymptotic marginal pgf and mean value of the system content in queue 2 always exist.Exploiting the formal resemblance between A * 2 (y) and A 2 (y), we easily find similar expressions as in the unconditional case (see (92)): The dominant pole σ of U * 2 (z) is the (only) zero of its denominator, i.e., i.e., a similar formula as for τ i in (93).As before, the (conditional) tail distribution of u 2 , given that u 1 → ∞, can be well approximated by using the dominant-pole technique, i.e., for large , where Hence, combining (96), ( 97) and (98), we get, for large , Again, this can be combined with a tail approximation for the system-content distribution in queue 1, as per Eqs.( 15) and ( 16), i.e., for large values of n, Combining (99) and (100), we find that, for large values of both n and , In terms of the original parameters μ and α of our model, u 1 (n, ) can be expressed as Interchanging the roles of queues 1 and 2, we can similarly obtain a dual approximation: It may seem surprising that expressions (102) and (103) are not identical, i.e., are not invariant to the substitutions α ↔ 1 − α, n ↔ .In order to clarify this observation, let us consider the special case α = 1 2 .Our model then reduces to the well-known model of the "(symmetric) clocked buffered switch", that has been analyzed by various authors in the past (Jaffe [22,23], Boxma and van Houtum [6], Cohen [10], Adan et al. [1]), by using various kinds of techniques.For α = 1 2 , we get the following results for the dominant poles τ 1 , τ 2 and σ : The joint tail approximations u 1 (n, ) and u 2 (n, ), given in (102) and (103), reduce to Jaffe [22,23] also provides an approximative expression for the joint tail distribution of the two "queue contents" (q 1 , q 2 ), where the queue content q i in queue i is defined as the number of customers in queue i, excluding the customer in service, if any.It is clear that the "queue content" q i is related to the "system content" u i in queue i (as we use in this paper) by the simple equation q i = (u i − 1) + , and, therefore, for large values of n and , Our two tail approximations u 1 (n, ) and u 2 (n, ) for Prob[u 1 = n, u 2 = ] can thus be translated to two corresponding tail approximations q 1 (n, ) and q 2 (n, ) for Prob[q 1 = n, q 2 = ] as follows: We illustrate the accuracy of formula (104) for large n (say n ≥ 10) by one specific example, in Fig. 5, where a comparison is made with simulation results.Producing simulation results for high n and at the same time, is extremely time-consuming, since it concerns very rare events; evaluating the closed formula (104), on the contrary, requires virtually no time at all.
The corresponding result obtained by Jaffe [22,23], in his original notations, is , where Translating this to our notations, we thus find that b 1 = τ 1 , b 2 = σ and Hence, we see that Jaffe's approximation is symmetric in the parameters n and , and corresponds to the sum of our two approximations; q 3 (n, ) is compared to simulation results in Fig. 6 and is clearly more accurate than q 1 (n, ) for lower values of n.In order to better understand why this is the case, we need some context.First, it is important to know that both Cohen [10] and Jaffe [22,23] prove that U (z 1 , z 2 ) has an infinite number of poles, and that the locations of the poles of U (z 1 , z 2 ), considered as a function of z 1 , are independent of z 2 , and vice versa.Thus, the joint pgf U (z 1 , z 2 ) has two countably infinite sets of poles, one set with respect to the variable z 1 and another set with respect to the variable z 2 .Since α = 1 2 , these two sets are identical.Using our notations, it turns out that, in both sets of poles, the dominant pole is τ 1 = τ 2 and the second most dominant pole is σ .The reason why Jaffe obtains a more accurate result is that he -explicitly, from the very beginning -not only uses the most dominant poles, but also the second most dominant poles of U (z 1 , z 2 ) with respect to the dimensions z 1 and z 2 , in his joint tail approximation, whereas we basically only use the most dominant pole in one direction, and then look for the most dominant pole in the other direction, given this earlier choice, which turns out to be the second most dominant pole in that direction.Apparently, as also explicitly noticed by Jaffe, there is no contribution to the joint tail from the two most dominant poles in both directions (i.e., τ 1 = τ 2 ) together.We further note that, if n is much larger than , then, since b 1 = τ 1 = τ 2 < b 2 = σ , q 2 (n, ) becomes negligible with respect to q 1 (n, ); the opposite is true if is much larger than n; in both cases, Jaffe's approximation basically reduces to one of our two approximations.
Finally, we emphasize that our approach does lead to the same, more accurate, approximation as Jaffe's analysis, if we simply add up our two approximations.This adding up is not just a heuristic trick: it corresponds to taking into account the contributions of more singularities (with respect to both z 1 and z 2 ) when inverting the pgf U (z 1 , z 2 ) to find the pmf Prob[u 1 = n, u 2 = ].This approach could be taken for other arrival processes A(x, y) too, if multiple expressions for Prob[u 1 = n, u 2 = ] are obtained.The reason we did not need it in Sect.5.2.1 (global geometric arrivals) is that there, for a given routing probability α, we only got one expression for Prob[u 1 = n, u 2 = ]; in Sect.5.2.2 (global shifted geometric arrivals), the two approximate expressions were identical.

Discussion of results
In the previous section, we have observed that, depending on the precise nature of the arrival pgf A(x, y), different scenarios may occur with respect to the behavior of queue 2 when the number of customers in queue 1 is large (i.e., in the asymptotic system).In case of identical arrivals, queue 2 is always unstable; for global geometric arrivals, queue 2 may be either stable or unstable, depending on the relative values of λ 1 and λ 2 , whereas in case of global shifted geometric or binomial arrivals, queue 2 is always stable.Intuitively, one might expect that the amount and the nature of the intra-slot correlation between the arrivals of type 1 and type 2, should play an important role in the type of behavior that is observed.If this correlation is positive, we expect both queues to exhibit similar behavior, i.e., the system content in queue 2 is large (small) when the system content in queue 1 is large (small); if the correlation is positive and large, we even expect that u 2 → ∞ when u 1 → ∞, and this is what we observe in case of identical arrivals, indeed.On the other hand, if the arrival correlation is negative, queues 1 and 2 are expected to behave differently, in the sense that u 2 typically takes small values when u 1 takes high values, and vice versa, and therefore, queue 2 should remain stable when u 1 → ∞, i.e., in the asymptotic system.This looks like what we observe in case of global shifted geometric or binomial arrivals.But it is not a priori clear what should happen if the arrival correlation is moderately positive; this seems to be what happens in case of global geometric arrivals.
Let us corroborate these intuitive expectations by means of formulas.We start by computing the arrival correlation.Let a 1 and a 2 denote the numbers of type-1 arrivals and type-2 arrivals in one slot, respectively.Then, the correlation coefficient ρ a between a 1 and a 2 is defined as where λ 1 = E[a 1 ] and λ 2 = E[a 2 ] are the quantities defined earlier in Eq. ( 3), and The quantity cov[a 1 , a 2 ] E[a 1 a 2 ] − λ 1 λ 2 is the covariance of a 1 and a 2 , and has the same sign as the correlation coefficient ρ a , but is easier to compute.

Global geometric arrivals
For the case of global geometric arrivals, where A(x, y) is given in Eq. (58), the arrival correlation turns out to be positive, between 0 and +1, i.e., We thus expect queue 2 in the asymptotic system (where u 1 → ∞) to be more heavily loaded than in the unconditional system.This expectation is confirmed by Eq. (63), i.e., Moreover, λ * 2 only remains smaller than one, i.e., the stability of queue 2 is only guaranteed, if λ 2 < λ 1 < 1. Apparently, the correlation coefficient ρ a (geom), which is an increasing function of λ 2 , for given λ 1 , gets too big to ensure stability of queue 2 beyond this point.But, if λ 2 < λ 1 , then E[u 2 |u 1 → ∞] exists and, according to (59) and (65), is greater than E[u 2 ]: Comparing the tail distributions of u 2 in the asymptotic system and the original system, we also see a heavier tail in the asymptotic system, as, according to (60) and (68), for large ,

Global shifted geometric arrivals
The situation is different for global shifted geometric arrivals, where A(x, y) is given by expression (72) and we find which is clearly negative, between 0 and −1, because 1 ≤ μ < 2. The mean asymptotic arrival rate in queue 2 is given by formula (77): where we have used the inequality (106) Equation ( 105) not only implies that queue 2 is always stable in the asymptotic system, but also that it is less heavily loaded than in the original system.This observation is confirmed by a comparison of the expressions (74) and (78) for the mean system content in queue 2: where we have, again, used the inequality (106).A comparison of the tail distributions of u 2 in the asymptotic system and the original system now reveals that both have the same decay rate, because, according to (79), in this case, the dominant poles τ 2 (of U 2 (z)) and σ (of U * 2 (z)) coincide: Nevertheless, we can derive from (81) and (84) that, for large , using, once again, the inequality (106).

Global binomial arrivals
The case of global binomial arrivals of order two can be dealt with similarly.Here, A(x, y) is given by (91) and, again, we find a negative arrival correlation, between 0 and −1, i.e., According to (94), the mean arrival rate in queue 2 in the asymptotic system is given by Here again, queue 2 is always stable in the asymptotic system, but less loaded than in the original system.Expressions (92), (95) for the mean system content in queue 2 confirm this: because the function is a monotonically increasing function of its argument x (since g (x) > 0, for all real x).
Comparing the tail distributions of u 2 in the asymptotic system and the original system, we also see a lighter tail in the asymptotic system, as, according to (99) and (100), for large ,

General global arrivals
The above results for global geometric, shifted geometric and binomial arrivals with probabilistic routing to both queues, suggest that the (sign of) the correlation coefficient (or, the covariance) between a 1 and a 2 gives a good indication of the type of behavior we may expect for queue 2 in the asymptotic system: if the arrival correlation is positive, then queue 2 receives more arrivals in the asymptotic system than in the original system, whereas the opposite implication seems to hold in case of negative arrival correlation.In an attempt to examine this more generally, we have found that the reality is somewhat more subtle, and, in particular, that a second-order characteristic of the joint arrival process such as the covariance or the correlation coefficient cannot completely predict the type of behavior of the asymptotic system.Let us reconsider -more generallythe case where a global arrival stream with pgf C(z) is split into two substreams with routing probabilities α and 1 − α for queues 1 and 2, respectively, as discussed in the introduction of Sect.5.2.Let us denote the mean and the variance of the total number of arrivals per slot in the system as E[c] and var[c], respectively.Then As the joint pgf A(x, y) is given by A(x, y) = C(αx + (1 − α)y) (see ( 55)), we easily get (107) Among other things, this simple result reveals that the arrival correlation is positive or negative depending on var[c] being larger or smaller than E[c].
We have seen in the earlier examples that λ * 2 was larger than λ 2 when the arrival correlation was positive, and vice versa when it was negative.Let us examine whether this also holds more generally.More specifically, we examine whether the quantities 2 λ * 2 − λ 2 and cov[a 1 , a 2 ] always have the same sign.In order to do so, we compute 2 as where we have used (47) and (55).From its definition, it follows that τ 1 satisfies the equation τ 1 = C(1 − α + ατ 1 ), and, therefore where D(z) C (z)−C (1)C(z).We now introduce the power series expansions of the functions C(z) and C (z) about z = 1: where C (n) (z) denotes the nth derivative of the function C(z).Using (109), we get D(z) as Inserting this result in (108), we obtain a useful expression for 2 : where the quantities ω n [a 1 , a 2 ] are, in fact, higher-order covariances of the arrival process, i.e., Expression (110) reveals that, indeed, there is a relationship between the quantities 2 and cov[a 1 , a 2 ], but cov[a 1 , a 2 ] does not completely determine 2 .It turns out that the higher-order covariances ω n [a 1 , a 2 ], n ≥ 2, of the arrival process are also important.This implies that arrival distributions may exist for which 2 and cov[a 1 , a 2 ] do not have the same sign.We have searched for an example, and, indeed, found one.Specifically, consider the case where the total number of arrivals per slot is either zero or two, i.e., the global arrival pgf C(z) takes the form In this case, the dominant pole τ 1 , solution of the equation τ 1 = C(1 − α + ατ 1 ), is given by and the ω n s are According to (110), 2 is equal to We thus have the situation that 2 is clearly negative for all possible values of the system parameters, whereas the arrival covariance can be either positive or negative, depending on μ being smaller or larger than one.So, in this particular example, the arrival covariance only predicts the behavior of the asymptotic system well in case the overall mean arrival rate μ is big enough.Whether this conclusion can be generalized to other global arrival distributions is not clear.The main conclusion from Sect. 6 seems to be that the arrival covariance cov[a 1 , a 2 ] is indeed related to the performance difference of queue 2 in a comparison between the original system and the asymptotic system, but does not simply predict it.The intuition that it should do so is thus not always correct.

Conclusions and further research
In this paper, we have considered a system of two parallel (coupled) discrete-time single-server queues with general independent arrivals per slot, but with possible intra-slot dependence between the arrival streams into both queues.To explicitly derive the joint steady-state pgf U (z 1 , z 2 ) of the system contents in queues 1 and 2 for a general joint arrival pgf A(x, y) is extremely hard for this system.Only for some classes of arrival pgfs, such explicit solutions are known (Bruneel [7]); we have provided two examples in Sect.5, i.e., the cases of identical arrivals (Sect.5.1) and global shifted-geometric arrivals (Sect.5.2.2).We have therefore focused on the derivation of the performance of queue 2 under the condition that the system content u 1 in queue 1 is very large.We have shown that, under this condition, queue 2 behaves as a regular discrete-time single-server queue with independent arrivals from slot to slot, but with a modified pgf (i.e., A * 2 (z) instead of A 2 (z)) for the number of arrivals per slot.Moreover, this is true for any choice of the joint pgf A(x, y), i.e., even in cases where an explicit solution for U (z 1 , z 2 ) is not known, and we have derived a simple formula to deduct A * 2 (z) from the original A(x, y).Using this fundamental result, we were able not only to determine the stability condition but also to derive all the usual performance characteristics of queue 2, such as the pgf U * 2 (z), the mean value E[u 2 |u 1 → ∞] (or any higher-order moment), and a very accurate approximation for the tail distribution of the system content of queue 2, in the so-called "asymptotic system".Moreover, we could also derive closed-form approximative expressions for the joint tail distribution of the system contents (u 1 , u 2 ) in the two queues in the original system.These approximations, too, are very accurate, especially when we combine the two dual formulas (one for the condition u 1 → ∞ and one for the condition u 2 → ∞) that are obtained with our approach in some cases.Whether this concept can be further generalized could be subject of further research.
We have then applied our theoretical results to a number of specific choices of the arrival pgf A(x, y), especially for the setting where one global arrival stream into the two-queue system is split probabilistically into two substreams, each entering one of the two individual queues.This application gave rise to a large number of explicit closed-form formulas for the main performance metrics of queue 2 in the asymptotic system, which could then be compared with their counterparts in the original unconditioned system.A serious attempt was made to intuitively explain the observed differences, based on the intra-slot correlation (or covariance) in the arrival streams, but it turned out that this was not always possible.
For some choices of A(x, y), where the explicit solution for U (z 1 , z 2 ) or approximative closed-form expressions for the joint tail distribution of (u 1 , u 2 ) have been reported earlier in the literature, we have verified that our solution technique gives rise to the same results.However, our technique is applicable for any A(x, y) and can thus generate new -partial -results for arrival processes that were, until now, not amenable to exact solution.In particular, we hope that the dominant-pole approximation formulas that can be derived with our technique for the joint tail distribution of (u 1 , u 2 ) may, in some cases, be helpful in determining more (or even all) singularities of the joint pgf U (z 1 , z 2 ), both with respect to the variable z 1 and the variable z 2 .The results in Sect.5.2.3 with respect to global binomial arrivals point in that direction.This too could be a topic of future work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.From (115) and (116), we conclude that the area of convergence of the joint pgf U (z 1 , z 2 ) contains at least the region (117) Now suppose that (z 1 , z 2 ) ∈ is a zero-tuple of the kernel K (z 1 , z 2 ), i.e., K (z 1 , z 2 ) = 0, then, since |U (z 1 , z 2 )| < ∞, Eq. (112) implies that L(z 1 , z 2 ) = 0, or, equivalently, Further, assume that, for some value of z 2 , both ( ẑ1 , z 2 ) ∈ and ( z1 , z 2 ) ∈ are zero-tuples of K (z 1 , z 2 ).We then refer to ẑ1 and z1 as "connected" z 1 -values.It then follows from (118) that, at the same time, We now show that, in this particular case, such connected z 1 -values can indeed be found.Specifically, we choose two z 1 -values on a circle with center in the origin and radius √ τ 1 in the complex z 1 -plane as follows: ẑ1 = √ τ 1 e ıθ and z1 = √ τ 1 e −ıθ .
Here ı denotes the imaginary unit, i.e., ı 2 = −1.From (113), it is clear that the equation K (z 1 , z 2 ) = 0 is equivalent to If we substitute z 1 in the above equation for either ẑ1 or z1 , we obtain the same equation in both cases: The above equation is a quadratic equation in the variable z 2 .It thus has two solutions for z 2 .
We now show that one of these solutions, say z 2 (θ ), is real and lies in the interval [0, 1].Indeed, for z 2 = 0, we have f (0) = −ατ 1 τ 2 < 0, whereas for z 2 = 1, we get Hence, the polynomial f (z) (a continuous function) goes from negative to positive between z 2 = 0 and z 2 = 1.Therefore, it has a zero, z Identifying equal (positive or negative) powers of z on both sides of the above equation leads to We note that u(0, 0) = U (0, 0) = 0 in this particular model, because the total number of arrivals during a slot into the system is always at least equal to 1 and, hence, the two-queue system can never be completely empty at the beginning of a slot.It thus follows that (121) The remaining unknown u(1, 0) can be computed from the known equality where U 2 (z) and λ 2 were derived from (73) and (56).Choosing z = 1 in (121), we easily get (122)