On the Nonsymmetric Longer Queue Model: Joint Distribution, Asymptotic Properties, and Heavy Traffic Limits

Abstract

We consider two parallel queues, each with independent Poisson arrival rates, that are tended by a single server. The exponential server devotes all of its capacity to the longer of the queues. If both queues are of equal length, the server devotes 𝜈 of its capacity to the first queue and the remaining 1−𝜈 to the second. We obtain exact integral representations for the joint probability distribution of the number of customers in this two-node network. Then we evaluate this distribution in various asymptotic limits, such as large numbers of customers in either/both of the queues, light traffic where arrivals are infrequent, and heavy traffic where the system is nearly unstable.

1. Introduction

We consider a nonsymmetric version of the longer queue model. Here there are two parallel queues, each fed by a Poisson arrival stream. There is but a single server who tends to the longer of the two queues. If the number of customers in each queue is the same, then the server devotes 𝜈 of its capacity to the first queue and 1−𝜈 to the second queue, with 0<𝜈<1. We let 𝑁₁ (𝑁₂) denote the number of customers in the first (second) queue, the two arrival rates are 𝜆₁ and 𝜆₂, and the server works at rate 𝜇. Note that the total number of customers, 𝑁₁+𝑁₂, in the two-node network behaves as the standard 𝑀/𝑀/1 model, so in the steady state we have the geometric distribution Prob[𝑁₁+𝑁₂=𝐽]=(1−𝜌₁−𝜌₂)(𝜌₁+𝜌₂)^𝐽, where 𝜌₁=𝜆₁/𝜇,𝜌₂=𝜆₂/𝜇, assuming the stability condition 𝜌₁+𝜌₂<1. The “symmetric case” corresponds to 𝜆₁=𝜆₂ (thus 𝜌₁=𝜌₂) and 𝜈=1/2, and this was analyzed in detail by Flatto [1].

Such models were proposed by Zheng and Zipkin [2] to study problems in inventory control. In [2] finite capacities were assumed in the two queues, and the authors studied numerically the steady state probabilities 𝑝(𝑚,𝑛)=Prob[𝑁₁=𝑚,𝑁₂=𝑛], in terms of the capacity size and also for different service disciplines, such as the longer queue (LQ) discipline here, and also the first-come-first-served discipline.

In [1] the author used two-dimensional generating functions and analyticity arguments and obtained explicit expressions for 𝑝(𝑚,𝑛), in the symmetric case, as contour integrals. Then asymptotic results were derived for the joint distribution 𝑝(𝑚,𝑛), as 𝑚 and/or 𝑛 becomes large, and also the marginal tails and various conditional limit laws were obtained. In this paper we generalize some of the results of [1] to the nonsymmetric model, and we will show that now many of the asymptotic results become quite different. As in [1] we assume that the model is preemptive, so if 𝑁₁=𝑁₂ (and thus the server works on both queues) and a new arrival occurs to the first queue, then the server switches immediately all its capacity to the first queue. The more difficult nonpreemptive version of the LQ model was studied by Cohen [3], in the case of Poisson arrivals and general service times. The problem is reduced in [3], using generating functions, to a functional equation which is furthermore converted to a Riemann-Hilbert boundary value problem.

The present problem corresponds to a random walk in the quarter plane (as 𝑁₁,𝑁₂⩾0), and other examples of such problems include shortest queue (SQ) problems [4–6], fork-join models [7, 8], two coupled processors with generalized processor sharing [9–11], and two coupled tandem queues [12]. General techniques for solving such problems are discussed in [13–15] and they involve functional equations, analyticity arguments, and singular integral equation methods [16].

For the present model the analyticity arguments are fairly simple, and we focus mostly on the asymptotic properties of the solution. We will show that these asymptotics are quite different from those of the symmetric model in [1]. After obtaining exact integral representations for 𝑝(𝑚,𝑛), and also the marginal probabilities 𝑃(𝑚)=Prob[𝑁₁=𝑚] and 𝒫(𝑛)=Prob[𝑁₂=𝑛], we asymptotically evaluate these integrals for 𝑚 and/or 𝑛 large. We use standard techniques, such as the Laplace method, saddle point method, singular analysis, and the Euler-MacLaurin formula. Good general references on the asymptotic evaluation of integrals and sums are the books [17–21].

In addition to this model being interesting on its own, many variants of shortest queue problems, such as ones with multiple servers and finite capacities, can be asymptotically reduced to LQ models of the type considered here (see [22, 23]). For example, in [23] we showed that the finite capacity version of the standard symmetric SQ model (analyzed in [4, 5]), where 𝑁₁,𝑁₂⩽𝐾 and 𝐾 is the capacity, asymptotically reduces to the symmetric LQ model in [1], if we consider the process (𝐾−𝑁₁,𝐾−𝑁₂), which measures the number of spots available in the two waiting rooms. Then having a thorough understanding of the nonsymmetric LQ model and its asymptotics will allow us to analyze, at least in some asymptotic limits, nonsymmetric variants of SQ models.

The remainder of this paper is organized as follows. In Section 2 we summarize all of the main results, both exact and asymptotic. They are listed in Theorems 1–5, and some discussion/interpretation appears following each theorem. In Section 3 we briefly derive the exact expressions for 𝑝(𝑚,𝑛); in Section 4 we derive asymptotic properties of 𝑝(𝑚,𝑛) for 𝑚,𝑛 large. In Section 5 we derive light traffic (where 𝜌₁,𝜌₂→0) and heavy traffic (where 𝜌₁+𝜌₂↑1) results.

2. Problem Statement and Summary of Results

We let (𝑁₁,𝑁₂) be the numbers of customers in the two parallel queues, and let 𝑝(𝑚,𝑛)=Prob[𝑁₁=𝑚,𝑁₂=𝑛] be the joint queue length distribution in the steady state. The two arrival rates are 𝜆₁ and 𝜆₂, the exponential server works at rate 𝜇, and 𝜌₁=𝜆₁/𝜇,𝜌₂=𝜆₂/𝜇. If 𝑁₁>𝑁₂(𝑁₁<𝑁₂) the server works on the first (second) queue, but if 𝑁₁=𝑁₂ the server works at rate 𝜇𝜈 on the first queue and rate 𝜇(1−𝜈) on the second, with 0<𝜈<1. The symmetric case corresponds to 𝜆₁=𝜆₂ and 𝜈=1/2. We henceforth assume the stability condition 𝜌₁+𝜌₂<1.

In Figure 1 we sketch the transition rates for the random walk (𝑁₁,𝑁₂), which illustrates the discontinuity along the “interface” 𝑁₁=𝑁₂. The main balance equations are 𝜌₁+𝜌₂+1𝑝(𝑚,𝑛)=𝜌₁𝑝(𝑚−1,𝑛)+𝜌₂𝑝()()𝜌𝑚,𝑛−1+𝑝𝑚+1,𝑛,𝑚>𝑛+1,𝑛⩾1,₁+𝜌₂+1𝑝(𝑚,𝑛)=𝜌₁𝑝(𝑚−1,𝑛)+𝜌₂𝑝()()𝑚,𝑛−1+𝑝𝑚,𝑛+1,𝑛>𝑚+1,𝑚⩾1. (1)(2) We also have the following three interface equations: 𝜌₁+𝜌₂+1𝑝(𝑛+1,𝑛)=𝜌₁𝑝(𝑛,𝑛)+𝜌₂𝑝()()+()𝑝()𝜌𝑛+1,𝑛−1+𝑝𝑛+2,𝑛1−𝜈𝑛+1,𝑛+1,𝑛⩾1,₁+𝜌₂+1𝑝(𝑛−1,𝑛)=𝜌₁𝑝(𝑛−2,𝑛)+𝜌₂𝑝()()()𝜌𝑛−1,𝑛−1+𝑝𝑛−1,𝑛+1+𝜈𝑝𝑛,𝑛,𝑛⩾2,₁+𝜌₂+1𝑝(𝑛,𝑛)=𝜌₁𝑝(𝑛−1,𝑛)+𝜌₂𝑝()()()𝑛,𝑛−1+𝑝𝑛+1,𝑛+𝑝𝑛,𝑛+1,𝑛⩾1, (3)(4)(5) the two boundary conditions 𝜌₁+𝜌₂+1𝑝(𝑚,0)=𝜌₁𝑝()()𝜌𝑚−1,0+𝑝𝑚+1,0,𝑚⩾2,₁+𝜌₂+1𝑝(0,𝑛)=𝜌₂𝑝(0,𝑛−1)+𝑝(0,𝑛+1),𝑛⩾2, (6)(7) and the corner conditions 𝜌₁+𝜌₂+1𝑝(1,0)=𝜌₁𝑝()+()𝑝()(),𝜌0,01−𝜈1,1+𝑝2,0₁+𝜌₂+1𝑝(0,1)=𝜌₂𝑝()()(),𝜌0,0+𝜈𝑝1,1+𝑝0,2₁+𝜌₂𝑝(0,0)=𝑝(1,0)+𝑝(0,1). (8) The normalization is ^∞∞_𝑛=0𝑚=0𝑝(𝑚,𝑛)=^∞_𝑛=0𝑝(𝑛,𝑛)+_𝑚>𝑛𝑝()+𝑚,𝑛_𝑚<𝑛𝑝(𝑚,𝑛)=1. (9) Note that the elementary difference equations (6) and (7) may be solved immediately to obtain 𝑝(𝑚,0) and 𝑝(0,𝑛), up to multiplicative constants. In Section 3 we analyze (1)–(9) to obtain 𝑝(𝑚,𝑛) in the forms of contour integrals, which we summarize below.

Figure 1 

A sketch of the transition rates for the random walk.

Theorem 1. For 𝜌₁+𝜌₂<1, the steady state distribution is as follows: (i)𝑚>𝑛: 𝑝(𝑚,𝑛)=1−𝜌₁−𝜌₂𝑓(𝑤)2𝜋𝑖𝑤^𝑛+1𝑍₊(𝑤)^𝑚−𝑛𝑑𝑤, (10) where the integral is over a small loop about 𝑤=0, 𝑍₊(𝑤)=12𝜌₁𝜌₁+𝜌₂+1+𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂,𝑓(𝑤)=𝑤+1num_𝑓(𝑤)(𝑤),dennum_𝑓(𝑤)=−2𝜌₁𝜌₂𝑤+𝜌²₂(1−𝜈)+𝜌₁𝜌₂−𝜌₂(1−𝜈)+𝜌²₁+𝜌₁𝜈+𝜌₂−𝜌₁+𝜌₂𝜈𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁,(𝑤)=𝜌𝑤+1den₁𝜈+𝜌₂(1−𝜈)𝜌₁+𝜌₂+1−𝜌₁𝜌₂𝑤+𝜈𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁𝑤+1−𝜌₂𝜌₁𝑤+1−𝜈𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂𝑤+1. (11)(12)(13)(14) (ii)𝑚<𝑛: 𝑝(𝑚,𝑛)=1−𝜌₁−𝜌₂𝑔(𝑤)2𝜋𝑖𝑤^𝑚+1𝑍^∗₊(𝑤)^𝑛−𝑚𝑔(𝑤)=𝑑𝑤,num_𝑔(𝑤)(𝑤),dennum_𝑔(𝑤)=−2𝜌₁𝜌₂𝑤+𝜌²₂(1−𝜈)+𝜌₁𝜌₂+𝜌₂(1−𝜈)+𝜌²₁−𝜌₁𝜈+𝜌₁+𝜌₂𝜈−𝜌₂×𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂,𝑍𝑤+1_∗+(𝑤)=12𝜌₂𝜌₁+𝜌₂+1+𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁𝑤+1. (15)(16)(17)(18)(iii)𝑚=𝑛: 𝑝(𝑛,𝑛)=1−𝜌₁−𝜌₂ℎ(𝑤)2𝜋𝑖𝑤^𝑛+1ℎ(𝑤)=(𝑤)𝑑𝑤,num(𝑤),(𝑤)dennum=−2𝜌₁𝜌₂𝑤+𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−𝜌₁𝜈𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁𝑤+1−𝜌₂(1−𝜈)𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂𝑤+1. (19)(20)(21)

We next evaluate 𝑝(𝑚,𝑛) in various asymptotic limits, to gain more insight into the structure of the joint distribution. Writing 𝑝(𝑚,𝑛)=𝑝(𝑚,𝑛;𝜌₁,𝜌₂,𝜈) to emphasize the dependence on the model parameters, we clearly have the symmetry relation 𝑝𝑚,𝑛;𝜌₁,𝜌₂,𝜈=𝑝𝑛,𝑚;𝜌₂,𝜌₁,1−𝜈. (22) Thus it is sufficient for the asymptotics to assume that 𝜌₁⩾𝜌₂, (23) and we also note that the expressions in Theorem 1 are consistent with (22).

We will show that the asymptotics are quite different whether (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂, (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ or (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂. In Figure 2 we sketch the curve (𝜌₁+𝜌₂)²=𝜌₁−𝜌₂ in the (𝜌₁,𝜌₂) parameter plane, for 𝜌₁+𝜌₂<1. Note that the curve may also be written as 𝜌₂=12−1−2𝜌₁+1+8𝜌₁ (24) and passes through the points (0,0) and (1,0). We first give results for 𝑝(𝑚,𝑛) when (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, and note that this includes the symmetric case of Flatto [1].

Figure 2 

A sketch of the parameter domain and the transition curve (𝜌1+𝜌2)2=𝜌1−𝜌2.

Theorem 2. For 𝜌₁⩾𝜌₂ and (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ one lets 𝑅=𝑚/𝑛 be the queue length ratio. The following asymptotic expressions then hold:(i)𝑛=𝑂(1), 𝑚→∞𝑝()∼𝑚,𝑛1−𝜌₁−𝜌₂𝑓(0)𝑚^𝑛𝜌𝑛!²𝜌₁+𝜌₂+1²−4𝜌₁^𝑛×2𝜌₁𝜌₁+𝜌₂+1+𝜌₁+𝜌₂+1²−4𝜌₁^𝑚. (25)(ii)𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(𝑅_𝑐,∞), 𝑅_𝑐𝜌=1+¹+𝜌₂²+𝜌₂𝜌₁+𝜌₂²+𝜌₂−𝜌₁𝜌₁𝜌₂𝑝()∼>1,𝑚,𝑛1−𝜌₁−𝜌₂√𝑓2𝜋𝑛⁰(𝑅)𝑀(𝑅)1𝑤_𝑠^𝑛1𝑍₊𝑤_𝑠^𝑚−𝑛, (26)(27) where 𝑤_𝑠=𝑤_𝑠(𝑅)=𝜌₁+𝜌₂+1²𝑅−4(𝑅+1)𝜌₁+𝜌₁+𝜌₂+1𝑅²𝜌₁+𝜌₂+1²−4𝜌₁𝑅²×−12𝜌₁𝜌₂(𝑅+1)²⁻¹,𝑍₊𝑤_𝑠=12𝜌₁()×𝑅𝜌𝑅+1₁+𝜌₂++1𝑅²𝜌₁+𝜌₂+1²−4𝜌₁𝑅²,𝑀(𝑅)=−1𝑅+1×𝑅−1𝑅²𝜌₁+𝜌₂+1²−4𝜌₁𝑅²−1𝑅²𝜌₁+𝜌₂+1²−4𝜌₁𝑅²−1−𝜌₁+𝜌₂+1 (28)(29)(30) and 𝑓₀(𝑅)=𝑓(𝑤_𝑠)=𝑓(𝑤_𝑠(𝑅)) can be computed from (12)–(14).(iii)𝑚,𝑛→∞ with 𝑅=𝑅_𝑐+𝜉/√𝑛 (thus 𝑅−𝑅_𝑐=𝑂(𝑛^−1/2)) 𝑝(𝑚,𝑛)∼𝐶₊𝜌₁+𝜌₂^2𝑛𝜌₁𝜌₁+𝜌₂𝜌₁+𝜌₂²+𝜌₂^𝑚−𝑛×1√2𝜋^∞𝜉₁𝑒^−𝑢2/2𝜉𝑑𝑢,₁1=𝜉𝑅_𝑐𝑅−1𝑅^𝑐_𝑐+−12𝜌₁𝜌₂𝜌₁+𝜌₂²+𝜌₂−𝜌₁^2−1/2,𝐶₊=1−𝜌₁−𝜌₂×𝜌₁+𝜌₂⁴−𝜌₁−𝜌₂²𝜌₁𝜌₁+𝜌₂²×𝜌₁+𝜌₂²(1−𝜈)+𝜌₁𝜌₁+𝜌₂²+𝜌₂𝜈+𝜌₁().1−𝜈 (31)(32)(33)(iv)𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(1,𝑅_𝑐)𝑝(𝑚,𝑛)∼𝐶₊𝜌₁+𝜌₂^2𝑛𝜌₁𝜌₁+𝜌₂𝜌₁+𝜌₂²+𝜌₂^𝑚−𝑛. (34)(v)𝑚=𝑛→∞𝑝(𝑛,𝑛)∼𝐶₀𝜌₁+𝜌₂^2𝑛,𝐶₀=1−𝜌₁−𝜌₂𝜌₁+𝜌₂²𝜌₁+𝜌₂⁴−𝜌₁−𝜌₂²𝜌₁+𝜌₂²+𝜌₂𝜈+𝜌₁(1−𝜈). (35)(36)(vi)𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(𝑅^∗_𝑐,1), 𝑅^∗_𝑐≡𝜌₁𝜌₂𝜌₁𝜌₂+𝜌₁+𝜌₂²+𝜌₁−𝜌₂𝜌₁+𝜌₂²+𝜌₁𝑝()<1,𝑚,𝑛∼𝐶₋𝜌₁+𝜌₂^2𝑚𝜌₂𝜌₁+𝜌₂𝜌₁+𝜌₂²+𝜌₁^𝑛−𝑚,𝐶₋=1−𝜌₁−𝜌₂𝜌₁+𝜌₂⁴−𝜌₁−𝜌₂²𝜌₂𝜌₁+𝜌₂²×𝜌₂+𝜈𝜌₁+𝜌₂²𝜌₁+𝜌₂²+𝜌₂𝜈+𝜌₁().1−𝜈 (37)(38)(39)(vii)𝑚,𝑛→∞ with 𝑅=𝑅^∗_𝑐+𝜂/√𝑛 (thus 𝑅−𝑅^∗_𝑐=𝑂(𝑛^−1/2)) 𝑝(𝑚,𝑛)∼𝐶₋𝜌₁+𝜌₂^2𝑚×𝜌₂𝜌₁+𝜌₂𝜌₁+𝜌₂²+𝜌₁^𝑛−𝑚1√2𝜋^𝜂1−∞𝑒^−𝑢2/2𝜂𝑑𝑢,₁1=𝜂1−𝑅^∗_𝑐𝑅^∗_𝑐1−𝑅^∗_𝑐+2𝜌₁𝜌₂𝑅^∗_𝑐𝜌₁+𝜌₂²+𝜌₁−𝜌₂^2−1/2. (40)(41)(viii)𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(0,𝑅^∗_𝑐)𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂√𝑔2𝜋𝑛⁰(𝑅)𝑀_∗(𝑅)1𝑤^∗_𝑠^𝑚1𝑍^∗₊𝑤^∗_𝑠^𝑛−𝑚,𝑤^∗_𝑠=𝑤^∗_𝑠(𝑅)=𝑅𝜌₁+𝜌₂+1²−4𝜌₂()+𝜌𝑅+1₁+𝜌₂+1𝜌₁+𝜌₂+1²−4𝜌₂1−𝑅²×2𝜌₁𝜌₂(𝑅+1)²⁻¹,𝑍^∗₊𝑤^∗_𝑠=12𝜌₂()×𝑅+11+𝜌₁+𝜌₂+1+𝜌₁+𝜌₂²−4𝜌₂1−𝑅²,𝑀_∗(𝑅)=𝑅1+𝑅×1−𝑅𝜌₁+𝜌₂+1²−4𝜌₂1−𝑅²𝜌₁+𝜌₂+1²−4𝜌₂1−𝑅²−𝑅𝜌₁+𝜌₂+1 (42)(43)(44)(45) and 𝑔₀(𝑅)=𝑔(𝑤^∗_𝑠)=𝑔(𝑤^∗_𝑠(𝑅)) can be computed from (14), (16), and (17).(ix)𝑚=𝑂(1), 𝑛→∞𝑝()∼𝑚,𝑛1−𝜌₁−𝜌₂𝑔(0)𝑛^𝑚𝜌𝑚!¹𝜌₁+𝜌₂+1²−4𝜌₂^𝑚×2𝜌₂1+𝜌₁+𝜌₂+1+𝜌₁+𝜌₂²−4𝜌₂^𝑛. (46)

Note that, in view of (12) and (16), 𝑓(0)=𝜌²₂(1−𝜈)+𝜌₁𝜌₂−𝜌₂(1−𝜈)+𝜌²₁+𝜌₁𝜈+𝜌₂−𝜌₁+𝜌₂𝜈𝜌₁+𝜌₂+1²−4𝜌₂×𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−𝜌₁𝜈𝜌₁+𝜌₂+1²−4𝜌₂−𝜌₂(1−𝜈)𝜌₁+𝜌₂+1²−4𝜌₁⁻¹,𝑔(0)=𝜌²₂(1−𝜈)+𝜌₁𝜌₂+𝜌₂(1−𝜈)+𝜌²₁−𝜌₁𝜈+𝜌₁+𝜌₂𝜈−𝜌₂𝜌₁+𝜌₂+1²−4𝜌₂×𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−𝜌₁𝜈𝜌₁+𝜌₂+1²−4𝜌₂−𝜌₂(1−𝜈)𝜌₁+𝜌₂+1²−4𝜌₁⁻¹. (47)

In Theorem 2 we listed the expansions of 𝑝(𝑚,𝑛) in order of decreasing 𝑅, from 𝑅=∞ in (25) to 𝑅=0 in (46). In the symmetric case 𝜌₁=𝜌₂ and 𝜈=1/2, and then 𝑝(𝑚,𝑛)=𝑝(𝑛,𝑚) and Theorem 2 is consistent with the results of Flatto in [1] (there the cases 𝑚=𝑂(1), 𝑛→∞ and 𝑛=𝑂(1), 𝑚→∞ were not considered, and the limits where 𝑅≈𝑅_𝑐,𝑅^∗_𝑐 were not treated in as much detail). Note that (34) and (38) correspond to “product form” approximations to 𝑝(𝑚,𝑛). The expressions in (25) and (46) are actually exact when 𝑛=0 and 𝑚=0, respectively, in view of our comments below (9). Our analysis in Section 4 will also indicate how to compute higher order correction terms in the various asymptotic series.

Next we take (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂. Now some of the asymptotic results for 𝑝(𝑚,𝑛) will be very different, while in other ranges of 𝑅=𝑚/𝑛 they will be similar to those in Theorem 2. We emphasize below the formulas that are different.

Theorem 3. For 𝜌₁>𝜌₂ and (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂ one has the following asymptotic expansions:(i)𝑚,𝑛→∞ with 𝑚>𝑛 and 𝑚−𝑛=𝑂(1)𝑝(𝑚,𝑛)∼𝐶₊√𝑤_𝑏√𝜋𝑛^−3/2𝑤_𝑏^−𝑛2𝜌₁𝜌₁+𝜌₂+1^𝑚−𝑛×2√𝜌₁𝜌₂𝜌₁+𝜌₂()+1𝑚−𝑛+Δ₊,𝑤_𝑏=𝜌₁+𝜌₂+1²−4𝜌₁4𝜌₁𝜌₂,𝐶₊=1−𝜌₁−𝜌₂21−√𝜌₁−𝜌₂√𝜌₁−𝜌₂−𝜌₁+𝜌₂×√𝜌₁−𝜌₂+1/2−(𝜈−1/2)𝜌₁+𝜌₂(1/2)𝜌₁+𝜌₂+1−(𝜈−1)√𝜌₁−𝜌₂,Δ₊=2√𝜌₁𝜌₂𝜌₂(1−𝜈)−𝜌₁+14𝜌₁+𝜌₂+1²×√𝜌₁−𝜌₂−𝜌₁+𝜌₂×1−√𝜌₁−𝜌₂×12𝜌₁+𝜌₂+1−(𝜈−1)√𝜌₁−𝜌₂⁻¹. (48)(49)(50)(51)(ii)𝑚=𝑛→∞𝑝(𝑛,𝑛)∼𝐶√𝑤_𝑏√𝜋𝑛^−3/2𝑤_𝑏^−𝑛,𝐶=1−𝜌₁−𝜌₂𝑤_𝑏𝜌₁𝜌₂^3/2×12+√𝜌₁−𝜌₂−1𝜈−2𝜌₁+𝜌₂×√𝜌₁−𝜌₂−𝜌₁+𝜌₂²×12𝜌₁+𝜌₂+1−(𝜈−1)√𝜌₁−𝜌₂²⁻¹. (52)(53)(iii)𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(𝑅_𝑑,1) (including the limit 𝑛−𝑚=𝑂(1)), 𝑅_𝑑=𝜌₁+𝜌₂+1²−4𝜌₁𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂²−4𝜌₂𝑝()∼<1,𝑚,𝑛𝐶₋𝑅−𝑅_𝑑^−3/2√𝑤_𝑏√𝜋𝑛^−3/2𝑤_𝑏^−𝑚×2𝜌₂𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂^𝑛−𝑚,𝐶₋=𝜌₁−𝜌₂^3/41−𝜌₁−𝜌₂√𝜌₁𝜌₂×𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂(1/4)𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂²−𝜌₂^3/2×12+√𝜌₁−𝜌₂+𝜌₁+𝜌₂12×1−𝜈4𝜌₁+𝜌₂+1²−𝜌₁(×√1−𝜈)𝜌₁−𝜌₂−𝜌₁+𝜌₂²×12𝜌₁+𝜌₂+1+√𝜌₁−𝜌₂(1−𝜈)²⁻¹. (54)(55)(56)(iv)𝑚,𝑛→∞ with 𝑅=𝑅_𝑑+𝑆/√𝑛 (thus 𝑅−𝑅_𝑑=𝑂(𝑛^−1/2)) 𝑝(𝑚,𝑛)=1−𝜌₁−𝜌₂√𝑤2𝜋𝑛^𝑏−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛×𝑔𝑤_𝑏𝑀_𝑑𝑒^−𝑆21/2−𝐾√𝑤_𝑏𝑀_𝑑^3/4𝑒^−𝑆21/4×𝐷_1/2−𝑆₁𝑛^1/4+𝑂𝑛^−1/2,𝑀_𝑑=𝑅_𝑑1−𝑅_𝑑+𝑅_𝑑𝜌₁𝜌₂𝑤_𝑏2𝜌₁−𝜌₂,𝑆₁=11−𝑅_𝑑1𝑀_𝑑2√𝑆,𝐾=𝜌₁𝜌₂12+√𝜌₁−𝜌₂+𝜌₁+𝜌₂12×1−𝜈4𝜌₁+𝜌₂+1²−𝜌₁(×√1−𝜈)𝜌₁−𝜌₂−𝜌₁+𝜌₂²×12𝜌₁+𝜌₂+1+√𝜌₁−𝜌₂(1−𝜈)²⁻¹, (57)(58)(59)(60) and 𝐷_1/2(⋅) is the parabolic cylinder function of order 1/2.(v) For 𝑛=𝑂(1) and 𝑚→∞, (25) applies. For 𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(1,∞), (27)–(30) hold. For 𝑚,𝑛→∞ with 𝑅=𝑚/𝑛∈(0,𝑅_𝑑), (42)–(45) hold. For 𝑚=𝑂(1) and 𝑛→∞, (46) applies.

Thus when (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂ the asymptotics of 𝑝(𝑚,𝑛) are different in eight ranges of 𝑅=𝑚/𝑛, but the final result in four of these is the same as the case (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂. Note however that if (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂ (27) applies for all 𝑅>1, so the transition for 𝑅≈𝑅_𝑐 is now absent. Also, if (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂, (42) applies for 𝑅∈(0,𝑅_𝑑) rather than 𝑅∈(0,𝑅^∗_𝑐). In (57), 𝑔(𝑤_𝑏) is computed by setting 𝑤=𝑤_𝑏, cf. (49), in (16), and we also note that 2𝜌₂𝑍^∗₊(𝑤_𝑏)=𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂.

The results in (49)–(60) are very different in form from the symmetric case asymptotics in [1]. For example, the diagonal probabilities in (52) contain the additional algebraic factor 𝑛^−3/2. In (57) the parabolic cylinder function can be computed, for example, from 𝐷_𝑝(𝑍)=𝑒^𝑍2/4𝑖√2𝜋^Br𝑒^−𝑍𝑢𝑒^𝑢2/2𝑢^𝑝𝑑𝑢, (61) where Br is a vertical contour in the complex 𝑢-plane with Re(𝑢)>0, and we have 𝑝=1/2 for the correction term in (57). Also, 𝐷₀(𝑍)=𝑒^−𝑍2/4 and for 𝑝=𝑁=0,1,2,3,… we can express 𝐷_𝑁(𝑍) in terms of a Hermite polynomial. In (57) we give a two-term asymptotic approximation in the transition range 𝑅≈𝑅_𝑑. Also, it may be shown that the leading term in (57) is just the limiting form of (42), as 𝑅→𝑅_𝑑. The expression in (55) applies both to 𝑚,𝑛→∞ with 𝑅_𝑑<𝑅<1, and to 𝑚,𝑛→∞ with 𝑛−𝑚=𝑂(1). In the latter case we can approximate the factor (𝑅−𝑅_𝑑)^−3/2 by (1−𝑅_𝑑)^−3/2. The parabolic cylinder function 𝐷_1/2(−𝑆₁) has the asymptotic behaviors 𝐷_1/2−𝑆₁∼𝑒^−𝑆21/4−𝑆₁,𝑆₁𝐷⟶−∞,_1/2−𝑆₁∼𝑒^𝑆21/4𝑆₁^−3/2−1√2,𝑆₁⟶+∞, (62) so as 𝑆₁→+∞ in (57) (corresponding to 𝑅>𝑅_𝑑) the 𝑂(𝑛^−1/4) correction term begins to dominate the leading term (which occurs for 𝑒^−𝑆21/2≪𝑛^−1/4 or |𝑆₁|≫(1/2)log𝑛). Then we can show that the expression in (57) for 𝑆₁→+∞ agrees with (55) for 𝑅↓𝑅_𝑑, and hence the two expansions asymptotically match. Note also that the scaling in (57) has 𝑚=𝑛𝑅_𝑑+𝑂(√𝑛).

Next we study the transition range (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂ in parameter space. This will lead to a new set of asymptotic results which will show, for example, how the formula in (52) for 𝑝(𝑛,𝑛) changes to the purely geometric approximation in (35). To quantify the closeness to the curve (𝜌₁+𝜌₂)²=𝜌₁−𝜌₂ in Figure 2 we write 𝜌₁+𝜌₂²−𝜌₁−𝜌₂=𝛿𝛽, (63) where 𝛿→0⁺ and 𝛽 can have either sign. If 𝛽=0 we are exactly on the transition curve. For small 𝛿, in certain ranges of 𝑅=𝑚/𝑛, the results in Theorems 2 and 3 still apply. For example, if 𝑚=𝑂(1) and 𝑛→∞ or 𝑛=𝑂(1) and 𝑚→∞ Theorems 2 and 3 agree, and then no transition range result is needed. As (𝜌₁+𝜌₂)²↓𝜌₁−𝜌₂, 𝑅_𝑐↓1 in (26), and then (27) will hold for all 𝑅∈(1,∞) in the transition case. Thus for 𝑅=𝑚/𝑛>1 the transition case will require a new asymptotic result only if 𝑅≈1. For 𝑅<1 the asymptotic result in (42) will apply for 𝑅<𝑅^∗_𝑐, where now, since (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂, 𝑅^∗_𝑐≈𝑅_𝑑 (cf. (37) and (54)). Then if 𝑅<1, we will need different asymptotic results only when 𝑅∈[𝑅⁰_𝑑,1], including 𝑅≈𝑅⁰_𝑑 and 𝑅≈1, where 𝑅⁰_𝑑=1−𝑎²1+4𝑎+11𝑎²,𝑎=𝜌₁+𝜌₂. (64) Note that 𝑅⁰_𝑑 is the limit of both 𝑅_𝑑 and 𝑅^∗_𝑐, if we replace 𝜌₁−𝜌₂ by (𝜌₁+𝜌₂)²=𝑎², as (37) leads to (setting 𝑎=𝜌₁+𝜌₂, 𝑏=√𝜌₁−𝜌₂) 𝑅^∗_𝑐=𝜌₁𝜌₂𝜌₁𝜌₂+𝜌₁+𝜌₂²+𝜌₁−𝜌₂𝜌₁+𝜌₂²+𝜌₁=(1/4)𝑎²−𝑏⁴(1/4)𝑎²−𝑏⁴+𝑎²+𝑏²𝑎²+(1/2)𝑎+𝑏²⟶(1/4)𝑎²1−𝑎²(1/4)𝑎²1−𝑎²+𝑎²3𝑎²=+𝑎1−𝑎²1+4𝑎+11𝑎², (65) since 𝛿→0⁺ corresponds to 𝑏→𝑎. Similarly, (54) leads to 𝑅_𝑑→𝑅⁰_𝑑 in this limit.

Since we will now have the state variables 𝑚 and 𝑛 large, and also 𝛿 small, it is necessary to relate these. In Theorem 4 we summarize the transition case results, scaling 𝑚 and 𝑛 in terms of 𝛿.

Theorem 4. For 𝜌₁>𝜌₂ and (𝜌₁+𝜌₂)²−(𝜌₁−𝜌₂)=𝛽𝛿 with 𝛿→0⁺, one has the following, where 𝑎=𝜌₁+𝜌₂ and 𝑏=√𝜌₁−𝜌₂.(i)𝑚=𝑛→∞, 𝑛=𝑦/𝛿²=𝑂(𝛿⁻²)𝑝()∼1𝑛,𝑛√2()√𝜋𝑛1−𝑎1−𝑎²𝑤3𝑎+1−2𝑎𝜈^𝑏−𝑛×1+2𝛽₁√𝑦𝑒^{𝛽21𝑦𝛽1√𝑦}_−∞𝑒^−𝑢2,√𝑑𝑢𝛽=𝑎1−𝑎²𝛽₁. (66)(67)(ii)𝑛=𝑦/𝛿²=𝑂(𝛿⁻²),𝑚−𝑛=𝑧/𝛿=𝑂(𝛿⁻¹),𝑚>𝑛𝑝(𝑚,𝑛)∼2(1−𝑎)1−𝑎11+𝑎√𝑤𝜋𝑛^𝑏−𝑛2𝜌₁𝜌₁+𝜌₂+1^𝑚−𝑛×𝑒^{−∼𝑧2/4}+2𝛽₁√𝑦𝑒^𝛽21𝑦𝑒^{−∼𝑧𝛽1√𝑦𝛽1√𝑦−∼𝑧/2}_−∞𝑒^−𝑢2𝑑𝑢_∼,𝑧=1−𝑎𝑧1+𝑎√𝑦=1−𝑎1+𝑎𝑚−𝑛√𝑛=𝑂(1). (68)(69)(iii)𝑛=𝑦/𝛿²=𝑂(𝛿⁻²),𝑅=𝑚/𝑛∈(𝑅⁰_𝑑,1],𝑚<𝑛𝑝(𝑚,𝑛)∼2(1−𝑎)1−𝑎11+𝑎√𝜋𝑛1−𝑎+2𝜈𝑎3𝑎+1−2𝜈𝑎1−𝑅⁰_𝑑𝑅−𝑅⁰_𝑑×𝑤_𝑏^−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛×1+2𝛽₂√𝑦𝑒^{𝛽22𝛽2√𝑦}_−∞𝑒^−𝑢2,√𝑑𝑢𝛽=𝑎1−𝑎²1−𝑅⁰_𝑑𝑅−𝑅⁰_𝑑𝛽₂. (70)(71)(iv)𝑚,𝑛→∞, 𝑅=𝑚/𝑛=𝑅⁰_𝑑+𝑆⁰/√𝑛 (thus 𝑅−𝑅⁰_𝑑=𝑂(𝑛^−1/2)) 𝑝(𝑚,𝑛)∼21−𝑎²𝜋1−𝑎+2𝜈𝑎𝑤1+3𝑎−2𝜈𝑎^𝑏−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛×𝑅⁰_𝑑11−𝑅⁰_𝑑+1−𝑎²8𝑎^2−1/4𝑛^−1/4𝑒^{−(𝑆01)2/4}×𝐷_−1/2−𝑆⁰₁,𝑆⁰₁=1𝑅⁰_𝑑1−𝑅⁰_𝑑11−𝑅⁰_𝑑+1−𝑎²8𝑎^2−1/2𝑆⁰, (72)(73) and 𝐷_−1/2(⋅) is the parabolic cylinder function of order −1/2  (see (61)).(v)𝑚,𝑛→∞, 𝑛=𝑦₁/𝛿⁴=𝑂(𝛿⁻⁴), 𝑅=𝑅⁰_𝑑+𝑆⁰/√𝑛=𝑅⁰_𝑑+𝑂(𝛿²)𝑝()∼2√𝑚,𝑛1−𝑎²𝑛^1/41−𝑎+2𝜈𝑎𝑅1+3𝑎−2𝜈𝑎^0𝑑11−𝑅⁰_𝑑+1−𝑎²8𝑎^2−1/4×𝑤_𝑏^−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛ℒ𝑆⁰₁;𝑦₁^1/4𝛽₃,𝛽₃=𝛽𝑎√1−𝑎²𝑅⁰_𝑑^−1/411−𝑅⁰_𝑑+1−𝑎²8𝑎^2−1/4 (74)(75) and ℒ is defined by the contour integral ℒ=1𝑍;𝛾2𝜋𝑖_Br+𝑒^𝑍𝑢𝑒^𝑢2/2−𝛾+√𝑢𝑑𝑢, (76) where Br₊ is to the right of all singularities of the integrand, including the pole at 𝑢=𝛾² if 𝛾>0.

If (𝜌₁+𝜌₂)²=𝜌₁−𝜌₂ then 𝛽=𝛽₁=0 and the term with the integral in (66) is absent. In Theorem 4 we wrote most results in terms of 𝑎=𝜌₁+𝜌₂ and 𝛽. As 𝛽(or𝛽₁)→+∞ it can be shown that (66) asymptotically matches to (35), and as 𝛽→−∞, (66) matches to (52). Similarly as 𝛽→∞ (68) matches to (27) while as 𝛽→−∞ (68) matches to (48). Some of these asymptotic matching relations are further discussed in Section 4. In (70), 𝛽₂ depends on 𝛽 and also on the ratio 𝑅=𝑚/𝑛, in view of (71). The results in Theorem 4 show that for 𝛿→0⁺ it is natural to scale 𝑛=𝑂(𝛿⁻²), and then either 𝑚=𝑂(𝛿⁻²) independently, or 𝑚−𝑛=𝑂(𝛿⁻¹). The expression in (74) is similar to, but somewhat more complicated than, that in (72). The function ℒ in (76) reduces to (2𝜋)^−1/2𝐷_−1/2(𝑍)𝑒^−𝑍2/4 in the limit 𝛾→0, and thus (72) may be viewed as a special case of (74)–(76). Note also that the results in Theorem 4 can be rewritten without introducing 𝛽 and 𝛿 in (63), as for example 𝛽₁√1𝑦=𝑎√1−𝑎²×𝜌₁+𝜌₂²−𝜌₁−𝜌₂√𝑛,𝑎=𝜌₁+𝜌₂, (77) and 𝛽₃𝑦₁^1/4 is proportional to 𝛽𝛿𝑛^1/4, and 𝛽𝛿 may be written in terms of the original parameters 𝜌₁,𝜌₂.

Next we consider some different asymptotic limits, those of “light” and “heavy” traffic. Light traffic corresponds to infrequent arrivals, where 𝜌₁,𝜌₂→0. Heavy traffic corresponds to nearly unstable systems, where 𝜌₁+𝜌₂↑1. It turns out that the present model has two possible heavy traffic limits. In the first, which we call HTL1, we have 𝜌₁+𝜌₂↑1 with a fixed 0<𝜌₂<1. Then most of the probability mass will occur in the range where 𝑚 and 𝑛 are large, but with 𝑚−𝑛=𝑂(1). More precisely, if 𝜀=1−𝜌₁−𝜌₂→0⁺ then 𝑚 and 𝑛 must be scaled to be 𝑂(𝜀⁻¹), but with the difference 𝑚−𝑛 fixed. In the second heavy traffic limit (HTL2) we again set 𝜀=1−𝜌₁−𝜌₂→0⁺ but now let 𝜌₂→0⁺, with 𝜌₂=𝑂(𝜀). Now the probability mass will become more spread out, with appreciable mass anywhere in the range 𝑋>𝑌>0, where 𝑚=𝑋/𝜀 and 𝑛=𝑌/𝜀. The light and heavy traffic results are summarized below as Theorem 5.

Theorem 5. (i) Light traffic: For 𝜌₁,𝜌₂→0 one has 𝑝(𝑚,𝑛)∼𝜌^𝑚₁𝜌^𝑛₂𝑛𝑚+𝑛. (78)
(ii) Heavy traffic limit 1:𝜀=1−𝜌₁−𝜌₂→0⁺ with 0<𝜌₂<1, 𝑛=𝑌/𝜀=𝑂(𝜀⁻¹),𝑚−𝑛=ℓ=𝑂(1)𝑝(𝑚,𝑛)∼𝜀4𝜌₂2−𝜈−𝜌₂2−𝜌₂+2𝜌₂−1𝜈1−𝜌₂1+𝜌₂^ℓ𝑒^−2𝑌𝑝(),ℓ⩾1,𝑛,𝑛∼𝜀4𝜌₂1−𝜌₂2−𝜌₂+2𝜌₂𝜈𝑒−1^−2𝑌𝑝()4,ℓ=0,𝑚,𝑛∼𝜀1−𝜌₂𝜌₂+𝜈2−𝜌₂+2𝜌₂𝜈𝜌−12−²𝜌₂^−ℓ𝑒^−2𝑌,ℓ⩽−1. (79)(80)(81)
(iii) Heavy traffic limit 2:𝜀=1−𝜌₁−𝜌₂→0⁺ with 𝜌₂=𝛾𝜀=𝑂(𝜀), 𝑛=𝑌/𝜀, 𝑚=𝑋/𝜀𝑝()∼𝜀𝑚,𝑛₂2𝜋𝑖_Br+𝑓₀(𝜉)𝑒^𝑌𝜉×exp𝑌−𝑋22𝛾+1+2𝛾+1²2𝑓+4𝛾𝜉𝑑𝜉,𝑋>𝑌>0,₀(𝜉)=𝜉𝜉+2𝛾+1−2𝛾+1²,𝑝()∼𝜀+4𝛾𝜉𝑛,𝑛²2𝜋𝑖_Br+ℎ₀(𝜉)𝑒^𝑌𝜉ℎ𝑑𝜉,₀(𝜉)=(2−𝜈)𝜉+2𝛾+1(1−𝜈)−(1−𝜈)2𝛾+1²+4𝛾𝜉(2−𝜈)𝜉+2𝛾+1−2𝛾+1²,𝑝()+4𝛾𝜉𝑚,𝑛∼𝜀𝜀𝛾2^𝑛−𝑚12𝜋𝑖_Br+𝑔(𝜉)𝑒^𝑌𝜉𝑑𝜉,𝑛>𝑚,𝑔(𝜀)=𝜈12−𝜈𝛾2𝛾+1²+4𝛾𝜉−2𝛾+12𝛾+1+𝜉−2𝛾+1²+4𝛾𝜉𝑑𝜉. (82)(83)(84)(85)(86)(87)  Here Br₊ is a vertical contour in the 𝜉-plane, which lies to the right of all singularities.

The expression in (78) applies for 𝑚>𝑛, 𝑚<𝑛, and 𝑚=𝑛, and in the light traffic limit the discontinuity of 𝑝(𝑚,𝑛) along the diagonal will appear only in the higher order terms. In HTL1, (79)–(81) show a piecewise geometric distribution in the ℓ variable, and an exponential density in 𝑌. For HTL2, writing (82) as 𝑝(𝑚,𝑛)∼𝜀²ℱ(𝑋,𝑌) we can easily show that ∫^∞₀∫^∞_𝑌ℱ(𝑋,𝑌)𝑑𝑋𝑑𝑌=1 so that to leading order the probability mass concentrates where 𝑚,𝑛=𝑂(𝜀⁻¹) with 𝑋>𝑌. From (84) we have 𝑝(𝑛,𝑛)=𝑂(𝜀²) but the total mass along the main diagonal is 𝑂(𝜀), which is smaller than the mass in (82). Then also 𝑝(𝑛−1,𝑛)=𝑂(𝜀²) with total mass ∑_𝑛𝑝(𝑛−1,𝑛)=𝑂(𝜀), which is comparable to that along the main diagonal. The diagonals with 𝑛−𝑚⩾2 have mass 𝑂(𝜀^𝑛−𝑚), which is smaller still. The integrands in (82), (84), and (86) have branch points at 𝜉=−(2𝛾+1)²/(4𝛾), are analytic at 𝜉=0, and may have poles at 𝜉=−2, if 𝛾>1/2.

This completes our summarization of the exact and asymptotic results. Despite the seeming complexity and the many separate cases, all the results follow from fairly standard asymptotic evaluations of the integrals in Theorems 1, as we will show in Sections 4 and 5.

3. The Exact Solution

We solve the difference equation(s) in (1)–(9) and thus obtain Theorem 1. We begin by introducing the three generating functions: 𝐹(𝑧,𝑤)=^∞∞_{𝑛=0𝑚=𝑛+1}𝑧^{𝑚−𝑛−1}𝑤^𝑛𝑝(),𝐺()=𝑚,𝑛𝑧,𝑤^∞_𝑛=1^𝑛−1_𝑚=0𝑧^{𝑛−𝑚−1}𝑤^𝑛𝑝(),𝐻(𝑤)=𝑚,𝑛_∞𝑛=0𝑝(𝑛,𝑛)𝑤^𝑛, (88)(89)(90) and we note that 𝐻(0)=𝑝(0,0), and 𝐹(0,𝑤)=^∞_𝑛=0𝑝(𝑛+1,𝑛)𝑤^𝑛,𝐺(0,𝑤)=^∞_𝑛=1𝑝(𝑛−1,𝑛)𝑤^𝑛. (91) Then (9) shows that 𝐹(1,1)+𝐺(1,1)+𝐻(1)=1. (92)

From (1) and (3) we obtain, after some calculation, 𝜌₁+𝜌₂+1−𝜌₁1𝑧−𝑧−𝜌₂𝑤𝑧𝐹(𝑧,𝑤)=𝜌₁𝐻(𝑤)+()𝐻(𝑤)()1−𝜈−𝑝0,0𝑤−1+𝜌₂𝑤𝑧𝐹(0,𝑤), (93) and (2) and (4) lead to 𝜌₁+𝜌₂𝜌+1−¹𝑧−𝜌₂1𝑤𝑧−𝐺()𝑤𝑧𝑧,𝑤=𝜌₂(𝑤)𝐻(𝑤)(−1𝑤𝐻+𝜈−𝑝0,0)𝑧𝜌₁+1𝑤𝐺(0,𝑤). (94) Equation (5) along 𝑚=𝑛 then leads to the following relation between 𝐹(0,𝑤), 𝐺(0,𝑤), and 𝐻(𝑤): 𝜌₁−𝜌₂𝐻(𝑤)()=𝜌+1−𝑝0,0₁+1𝑤𝐺(0,𝑤)+𝜌₂𝑤+1𝐹(0,𝑤). (95)

If we define 𝒮(𝜉)=^∞_{𝑁=0𝑚+𝑛=𝑁}𝑝(𝑚,𝑛)𝜉^𝑁, (96) then from (96) and (88)–(90) we find that 𝒮(𝜉)=𝐻𝜉²+𝜉𝐹𝜉,𝜉²+𝜉⁻¹𝐺𝜉⁻¹,𝜉². (97) Using (93) with (𝑧,𝑤)=(𝜉,𝜉²), (94) with (𝑧,𝑤)=(𝜉⁻¹,𝜉²) and (95) with 𝑤=𝜉², we obtain for 𝒮(𝜉) the simpler equation 𝜌₁+𝜌₂+1−𝜌₁+𝜌₂1𝜉−𝜉𝒮(𝜉)=11−𝜉𝑝(0,0) (98) so that 𝒮(𝜉)=𝑝(0,0)/[1−(𝜌₁+𝜌₂)𝜉] and then 𝒮(1)=1 by (92), and thus 𝑝(0,0)=1−𝜌₁−𝜌₂. Then from (96) we have _{𝑚+𝑛=𝑁}𝑝(𝑚,𝑛)=1−𝜌₁−𝜌₂𝜌₁+𝜌₂^𝑁, (99) so that 𝑁₁+𝑁₂ follows a geometric distribution with parameter 𝜌₁+𝜌₂, and hence 𝑁₁+𝑁₂ behaves as the standard 𝑀/𝑀/1 model.

The coefficient of 𝐹(𝑧,𝑤) in (93) has roots at 1𝑧=2𝜌₁1+𝜌₁+𝜌₂±𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂𝑤+1≡𝑍_±(𝑤), (100) while that of 𝐺(𝑧,𝑤) in (94) has roots at 1𝑧=2𝜌₂𝑤1+𝜌₁+𝜌₂±𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁≡𝑍𝑤+1_∗±(𝑤)𝑤. (101) By using (95) to express 𝐻(𝑤) in terms of 𝐹(0,𝑤) and 𝐺(0,𝑤), eliminating 𝐻(𝑤) in the right sides of (93) and (94), and requiring that 𝐹(𝑧,𝑤) be analytic at 𝑧=𝑍₋(𝑤) and that 𝐺(𝑧,𝑤) be analytic at 𝑧=𝑍^∗₋(𝑤)/𝑤, we obtain two equations for 𝐹(0,𝑤) and 𝐺(0,𝑤): 𝜌₂𝑤+1𝜌₂𝑤+𝜈𝜌₁+𝜌₂𝐹()+𝜌+10,𝑤₁𝑤+1𝜌₂𝑤+𝜈𝑤𝜌₁+𝜌₂−𝜌+1¹𝑤+1𝑍^∗₋(𝑤)𝐺()(𝜌0,𝑤=𝑝0,0)𝜈−²𝑤+𝜈𝜌₁+𝜌₂,𝜌+1₂𝑤+1𝜌₁𝑤+1−𝜈𝜌₁+𝜌₂𝑤−+11+𝜌₂𝑤𝑍₋(𝑤)𝐹()+𝜌0,𝑤₁𝑤+1𝜌₁𝑤+1−𝜈𝑤²𝜌₁+𝜌₂𝐺()(+10,𝑤=𝑝0,0)1−𝜈𝑤−𝜌₁𝑤+1−𝜈𝜌₁+𝜌₂𝑤.+1 (102) Solving the algebraic system in (98) and (99), and then using the result in (95) to compute 𝐻(𝑤), we hence obtain explicitly the right sides of (93) and (94). In particular, 𝐹(𝑧,𝑤)=1−𝜌₁−𝜌₂𝑓(𝑤)𝑍₊(𝑤)=−𝑧1−𝜌₁−𝜌₂^∞_𝐽=0𝑧^𝐽𝑍₊(𝑤)^𝐽+1𝑓(𝑤), (103) where 𝑓(𝑤) is as in (12)–(14). Comparing (100) to (88) we conclude that ^∞_𝑛=0𝑝(𝑛+1+𝐽,𝑛)𝑤^𝑛=1−𝜌₁−𝜌₂𝑓(𝑤)𝑍₊(𝑤)^𝐽+1 (104) for 𝐽⩾0. Setting 𝐽=𝑚−𝑛−1 and inverting the generating function in (101) leads to (10) for 𝑚>𝑛. Similarly, inverting the double generating function for 𝐺(𝑧,𝑤) in (89) leads to (15), and we note that (10) and (15) are consistent with the symmetry (𝑛,𝜌₁,𝜈,𝑍₊,𝑓)↔(𝑚,𝜌₂,1−𝜈,𝑍^∗₊,𝑔). Finally, 𝑝(𝑛,𝑛) in (19) is obtained by inverting the generating function 𝐻(𝑤)=(1−𝜌₁−𝜌₂)ℎ(𝑤). We have thus established Theorem 1.

4. Asymptotics of the Joint Distribution

We derive Theorems 2–4 by expanding asymptotically the integrals in Theorem 1. We will use a combination of the saddle point method and singularity analysis. Good general references on techniques for asymptotically evaluating integrals can be found in [17–21].

We need to understand the singularities of the integrands in (10), (15), and (19). There are clearly branch points where 𝑤=𝑤_𝑏 and 𝑤=𝑤^∗_𝑏, with 𝑤_𝑏=𝜌₁+𝜌₂+1²−4𝜌₁4𝜌₁𝜌₂,𝑤^∗_𝑏=𝜌₁+𝜌₂+1²−4𝜌₂4𝜌₁𝜌₂. (105) Since 𝜌₁⩾𝜌₂ we have 𝑤_𝑏⩽𝑤^∗_𝑏 and for 𝜌₁>𝜌₂ the branch point at 𝑤^∗_𝑏 is farther from the origin than the one at 𝑤_𝑏. In fact, 𝑤^∗_𝑏 will never play a role in the asymptotics. The integrands are also singular at 𝑤=0, where (10) has a pole of order 𝑛+1, and (15) has a pole of order 𝑚+1. The only possible other singular points are at the zeros of den(𝑤). We can easily verify that 𝑤=1 is a simple zero of all four functions num_𝑓(𝑤), num_𝑔(𝑤), num(𝑤), and den(𝑤), so all the integrands are analytic at 𝑤=1. In the appendix we study in detail the algebraic equation den(𝑤)=0, and show that the only possible zero is at 𝑤=𝑤_𝑝𝑤_𝑝=𝜌₁+𝜌₂⁻²,𝜌₁+𝜌₂²>𝜌₁−𝜌₂, (106) but only if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ (if (𝜌₁+𝜌₂)²=𝜌₁−𝜌₂ then den(𝑤) has a branch point at 𝑤=𝑤_𝑝(=𝑤_𝑏)). If (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂ then the two branch points are the only singularities of the functions 𝑓(𝑤), 𝑔(𝑤), and ℎ(𝑤) in (12), (16), and (19). If (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, which is clearly true in the symmetric case, then 𝑤=𝑤_𝑝 is a simple pole of these functions (since den^′(𝑤_𝑝)≠0). In view of (105) and (106), we have 1<𝑤_𝑝<𝑤_𝑏<𝑤^∗_𝑏, if the pole is present and the stability condition 𝜌₁+𝜌₂<1 holds.

Consider first (19). If (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ the pole at 𝑤_𝑝 determines the asymptotic behavior, as it is closer to 𝑤=0 than the branch point at 𝑤_𝑏. Hence as 𝑛→∞𝑝(𝑛,𝑛)∼−1−𝜌₁−𝜌₂lim_{𝑤→𝑤𝑝}𝑤−𝑤_𝑝ℎ(𝑤)𝑤^𝑛+1=−1−𝜌₁−𝜌₂𝑤_𝑝^−𝑛−1num𝑤_𝑝den^′𝑤_𝑝. (107) Straightforward computations, using (14) and (21), show that den^′𝑤_𝑝=2𝜌₁𝜌₂𝜌₁+𝜌₂²𝜌₁+𝜌₂²+𝜌₂𝜈+𝜌₁(1−𝜈)1−𝜌₁−𝜌₂𝜌₁+𝜌₂⁴−𝜌₁−𝜌₂²,num𝑤_𝑝=−2𝜌₁𝜌₂1−𝜌₁−𝜌₂𝜌₁+𝜌₂². (108)(109) Using (107)–(109) leads to (35), with (36).

Next we let 𝑚,𝑛→∞ simultaneously, setting 𝑅=𝑚/𝑛, still taking (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, and first consider 𝑅>1, where (10) applies. The integrand in (10) again has singularities at 𝑤_𝑝, 𝑤_𝑏, and 𝑤^∗_𝑏. But now we write 𝑤^−𝑛𝑍₊(𝑤)^𝑛−𝑚=exp𝑛(1−𝑅)log𝑍₊(𝑤)−log𝑤 (110) and let 𝐹(𝑤)=𝐹(𝑤;𝑅)=(1−𝑅)log(𝑍₊(𝑤))−log𝑤. Then with (110), (10) has saddle points where 𝐹^′(𝑤)=0, or ()𝑍1−𝑅^′+(𝑤)𝑍₊(𝑤)−1𝑤=0. (111) Solving (111) using (11), after some algebra we are led to the saddle 𝑤_𝑠=𝑤_𝑠(𝑅) in (28). There is also a second saddle at 𝑤_𝑠, where 𝑤_𝑠 corresponds to replacing √… in (28) by −√…. By evaluating 𝐹^′′(𝑤_𝑠) and 𝐹^′′(𝑤_𝑠) we can easily show that the directions of steepest descent at the saddle(s) are arg(𝑤−𝑤_𝑠)=±𝜋/2 and arg(𝑤−𝑤_𝑠)=0,𝜋. From (28) we have 𝑤_𝑠(1)=𝑤_𝑏, so as 𝑅=𝑚/𝑛↓1 the saddle approaches a branch point. Also, 𝑤_𝑠→0 as 𝑅→∞. It is possible for the pole and saddle to coalesce. Setting 𝑤_𝑠(𝑅)=𝑤_𝑝=(𝜌₁+𝜌₂)⁻² and solving this equation for 𝑅 we find, again after some algebra, that 𝑤_𝑠(𝑅)=𝑤_𝑝when𝑅=𝑅_𝑐, (112) where 𝑅_𝑐 is given by (26). If 𝑅>𝑅_𝑐(>1), we have 𝑤_𝑠<𝑤_𝑝 and if 𝑅∈(1,𝑅_𝑐), we have 𝑤_𝑠>𝑤_𝑝. From (28) we can show that 𝑤^′_𝑠(𝑅)<0 so that the saddle location decreases with the queue length ratio 𝑚/𝑛=𝑅. If 𝑅>𝑅_𝑐 we deform the contour in (10), which is a small counterclockwise loop about 𝑧=0, into a saddle point contour, on which |𝑤|=𝑤_𝑠(𝑅). Then the new contour traverses the saddle in the direction(s) of steepest descent and the standard Laplace estimate of (10) is 𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂2𝜋2𝜋𝑛𝐹^′′𝑤_𝑠𝑓𝑤_𝑠𝑤_𝑠exp𝑛𝐹𝑤_𝑠. (113) But 𝑒^{𝑛𝐹(𝑤𝑠)}=𝑤_𝑠^−𝑛[𝑍₊(𝑤_𝑠)]^𝑛−𝑚 and, in view of (111), 𝐹^′′(𝑤)=1𝑤²+()𝑍1−𝑅^+′′(𝑤)𝑍₊(𝑤)−𝑍^′₊(𝑤)𝑍₊(𝑤)². (114) Then using (11) we commute 𝑤²_𝑠𝐹^′′(𝑤_𝑠) to find that 𝑤²_𝑠(𝑅)𝐹^′′(𝑤_𝑠(𝑅)) is the same as 𝑀(𝑅) in (30). Thus (113) is the same as the right side of (27), and we have established the asymptotics of 𝑝(𝑚,𝑛) for 𝑅>𝑅_𝑐.

If 1<𝑅<𝑅_𝑐,𝑤_𝑝>𝑤_𝑠 and in deforming the loop in (10) to the saddle point contour |𝑤|=𝑤_𝑠(𝑅) we must take into account the contribution from the residue at the pole 𝑤_𝑝. But we have 𝐹(𝑤_𝑝)⩾𝐹(𝑤_𝑠) with equality only if 𝑅=𝑅_𝑐, when 𝑤_𝑝=𝑤_𝑠. Thus the pole contribution dominates the saddle contribution and we have 𝑝(𝑚,𝑛)∼−1−𝜌₁−𝜌₂𝑤_𝑝^−𝑛−1𝑍₊𝑤_𝑝^𝑛−𝑚num_𝑓𝑤_𝑝den^′𝑤_𝑝. (115) From (13) we find that num_𝑓𝑤_𝑝=−2𝜌₂𝜌₁+𝜌₂²𝜌₁+𝜌₂²(1−𝜈)+𝜌₁×1−𝜌₁−𝜌₂ (116) and from (11) 𝑍₊𝑤_𝑝=12𝜌₁𝜌₁+𝜌₂+1+𝜌₁+𝜌₂+1²−4𝜌₁−4𝜌₁𝜌₂𝜌₁+𝜌₂²=𝜌₁+𝜌₂²+𝜌₂𝜌₁𝜌₁+𝜌₂, (117) when (𝜌₁+𝜌₂)²+𝜌₂−𝜌₁>0. With (115)–(117) and (108) we have established (34) for 𝑅∈(1,𝑅_𝑐), where 𝑒^{𝑛𝐹(𝑤𝑠;𝑅)} gives the product form expression and (1−𝜌₁−𝜌₂)|num_𝑓(𝑤_𝑝)|/[𝑤_𝑝den^′(𝑤_𝑝)] corresponds to the multiplicative constant 𝐶₊.

As 𝑅→∞ we have 𝑤_𝑠→0 and we must then reconsider the asymptotics of (10). From (28) 𝑤_𝑠=𝑂(𝑅⁻¹) as 𝑅→∞ so we scale 𝑤=𝑢/𝑚 in the integral in (10) and consider the limit 𝑚→∞ with 𝑛=𝑂(1). Then 𝑍₊(𝑤)^𝑛−𝑚=𝑍₊(0)+𝑍^′₊(0)𝑢𝑚+𝑂𝑚^{−2𝑛−𝑚}∼𝑍₊(0)^𝑛−𝑚−𝑍exp^′+(0)𝑍₊(0)𝑢 (118) and (10) becomes asymptotically 𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂𝑓(0𝑍2𝜋𝑖)₊(0)^𝑛−𝑚×𝑚^𝑛𝑢^𝑛+1−𝑍exp^′+(0)𝑍₊(0)𝑢=𝑑𝑢1−𝜌₁−𝜌₂𝑓(0)𝑍₊(0)^𝑛−𝑚−𝑍^′₊(0)𝑍₊(0)^𝑛𝑚^𝑛.𝑛! (119) Since 𝑍^′₊(0)=−𝜌₂/(𝜌₁+𝜌₂+1)²−4𝜌₁ we see that (119) is the same as (25).

Now consider 𝑅=𝑚/𝑛≈𝑅_𝑐, where the saddle and pole are close to each other. This is a standard problem that is discussed, for example, in [17, 18]. We now expand the integrand in (10) about 𝑤=𝑤_𝑝, ultimately scaling 𝑤−𝑤_𝑝=𝑂(𝑛^−1/2), and then the integrand will approach a limiting value as 𝑛→∞. We have, by Taylor/Laurent series, 𝑓(𝑤)=num_𝑓𝑤_𝑝den^′𝑤_𝑝1𝑤−𝑤_𝑝+𝑂(1),𝑤⟶𝑤_𝑝,𝑛𝐹(𝑤)=−𝑛log𝑤_𝑝−(𝑚−𝑛)log𝑍₊𝑤_𝑝−𝑛𝑤_𝑝+𝑚−𝑛𝑍₊𝑤_𝑝𝑍^′₊𝑤_𝑝𝑤−𝑤_𝑝+12𝑛𝑤²_𝑝+()𝑍𝑚−𝑛^′+𝑤_𝑝𝑍₊𝑤_𝑝²−()𝑍𝑚−𝑛^+′′𝑤_𝑝𝑍₊𝑤_𝑝𝑤−𝑤_𝑝²+𝑂𝑛𝑤−𝑤_𝑝³. (120)(121) From (11) and (122) we find that 𝑍^′₊𝑤_𝑝𝑍₊𝑤_𝑝=1𝑤_𝑝𝑛𝑛−𝑚_{𝑅=𝑅𝑐}=1𝑤_𝑝11−𝑅_𝑐, (122) since, at 𝑅=𝑅_𝑐, 𝑤_𝑝=𝑤_𝑠. It follows that −𝑛𝑤_𝑝−𝑚−𝑛𝑍₊𝑤_𝑝𝑍^′₊𝑤_𝑝=𝑛𝑤_𝑝𝑅−𝑅_𝑐𝑅_𝑐=√−1𝑛𝜉𝑤_𝑝𝑅_𝑐,−1 (123) where 𝜉 was defined above (31). Also, from (11) we find that 𝑍₊^′′𝑤_𝑝=−2𝜌₁𝜌²₂𝜌₁+𝜌₂³𝜌₁+𝜌₂²+𝜌₂−𝜌₁³. (124) Then setting 𝑤−𝑤_𝑝=1√𝑛𝑤_𝑝𝑅_𝑐𝑅_𝑐+−12𝜌₁𝜌₂𝜌₁+𝜌₂²+𝜌₂−𝜌₁^2−1/2𝑣 (125) and scaling 𝜉 as in (32), (121) becomes 𝑛𝐹(𝑤)=−𝑛log𝑤_𝑝−(𝑚−𝑛)log𝑍₊𝑤_𝑝+𝜉₁1𝑣+2𝑣²+𝑂𝑛^−1/2. (126) Then we use (120) and (126) to get 𝑝(𝑚,𝑛)=𝑤_𝑝^−𝑛𝑍₊𝑤_𝑝^𝑛−𝑚×−𝐶₊12𝜋𝑖_Br−𝑒^𝜉1𝑣𝑒^𝑣2/21𝑣×𝑛𝑑𝑣1+𝑂^−1/2, (127) where Br₋ is a vertical contour with Re(𝑣)<0. To obtain (127) we shifted the original contour in (10) into the circle |𝑤|=𝑤_𝑝−𝜀₁, 𝜀₁>0, and note that |𝑤|<𝑤_𝑝 implies to leading order that Re(𝑣)<0. The integral in (127) can be expressed in terms of a parabolic cylinder function of order 𝑝=−1 (see (61)), which can be expressed in terms of the standard error function, using the identity −12𝜋𝑖_Br−1𝑣𝑒^𝜉1𝑣𝑒^𝑣2/21𝑑𝑣=√2𝜋^∞𝜉₁𝑒^−𝑢2/2𝑑𝑢. (128) With (128) and (117), (127) becomes the same as (31), so we have derived the leading term for the range 𝑅−𝑅_𝑐=𝑂(𝑛^−1/2). We have thus covered all of the necessary ranges that have 𝑅⩾1, and established items (i)–(v) in Theorem 2.

Now we consider ranges with 𝑅<1(𝑚<𝑛) where 𝑝(𝑚,𝑛) is given by the integral in (15). The analysis is completely analogous to the expansion of (10), so we merely sketch the details. The function 𝑔(𝑤) in (16) has a simple pole at 𝑤=𝑤_𝑝=(𝜌₁+𝜌₂)⁻²>1 if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, with residue num_𝑔(𝑤_𝑝)/den^′(𝑤_𝑝) where num_𝑔𝑤_𝑝=−2𝜌₁𝜌₁+𝜌₂²𝜌₁+𝜌₂²𝜈+𝜌₂1−𝜌₁−𝜌₂. (129) We write 𝑤^−𝑚[𝑍^∗₊(𝑤)]^𝑚−𝑛=exp[𝑛𝐹_∗(𝑤;𝑅)] where 𝐹_∗(𝑤;𝑅)=−𝑅log𝑤+(𝑅−1)log𝑍^∗₊(𝑤). (130) Thus for 𝑚,𝑛 simultaneously large 𝐹_∗ will have a saddle where 𝐹^′_∗(𝑤)=0, and this leads to 𝑤=𝑤^∗_𝑠=𝑤^∗_𝑠(𝑅), as in (43). The saddle and pole coalesce where 𝑅=𝑅^∗_𝑐(<1) in (37). For 𝑅∈(𝑅^∗_𝑐,1) we have 𝑤_𝑝<𝑤^∗_𝑠(𝑅) and the pole determines the asymptotic behavior of 𝑝(𝑚,𝑛). Then 𝑝(𝑚,𝑛)∼−1−𝜌₁−𝜌₂num_𝑔𝑤_𝑝den^′𝑤_𝑝𝑤_𝑝^−𝑚−1𝑍^∗₊(𝑤_𝑝)^𝑚−𝑛, (131) and this leads to (38) with (39), as now 𝑍^∗₊(𝑤_𝑝)=[(𝜌₁+𝜌₂)²+𝜌₁]/[𝜌₂(𝜌₁+𝜌₂)]. For 𝑅∈(0,𝑅^∗_𝑐) the saddle determines the asymptotics, and the estimate 𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂√2𝜋𝑛1𝐹_∗^′′𝑤^∗_𝑠𝑔𝑤^∗_𝑠𝑤^∗_𝑠exp𝑛𝐹_∗𝑤^∗_𝑠 (132) leads to (42), with (43)–(45). As 𝑅→0 the saddle 𝑤^∗_𝑠→0 and the result in (46) can be obtained by approximating the integrand in (15) for 𝑤=𝑂(𝑛⁻¹), similarly as in (118) and (119) (with 𝑍₊ replaced by 𝑍^∗₊, 𝑓 by 𝑔, and 𝑚 by 𝑛). Thus we obtain (46). Finally, when 𝑅≈𝑅^∗_𝑐 the pole and saddle are close. Then we scale 𝑅−𝑅^∗_𝑐=𝜂/√𝑛=𝑂(𝑛^−1/2) and expand the integrand in (15) near 𝑤=𝑤_𝑝, with 𝑤−𝑤_𝑝=𝑂(𝑛^−1/2). This leads ultimately to (40) with (41). We have thus sketched the derivation of items (vi)–(ix) in Theorem 2.

Next we take (𝜌₁+𝜌₂)<𝜌₁−𝜌₂(>0) and establish Theorem 3. Now the pole at 𝑤=𝑤_𝑝 is absent. We first consider the diagonal probabilities in (19). The only singularities of ℎ(𝑤) as the branch points at 𝑤_𝑏 and 𝑤^∗_𝑏, and the former determines the large 𝑛 asymptotics, as 𝑤_𝑏<𝑤^∗_𝑏 in view of (105). Expanding ℎ(𝑤) about 𝑤=𝑤_𝑏 will transfer immediately to an asymptotic series for 𝑝(𝑛,𝑛), but here we focus on the leading term. From (14) and (21) we have num(𝑤)=𝐴+𝐵√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏,den(𝑤)=𝐶+𝐷√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏, (133) where 𝐴=−2𝜌₁𝜌₂𝑤_𝑏+𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−2𝜈𝜌₁√𝜌₁−𝜌₂,𝐵=−2√𝜌₁𝜌₂𝜌₂(),𝜌1−𝜈𝐶=₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−2𝜌₁𝜌₂𝑤_𝑏+𝜈√𝜌₁−𝜌₂,𝐷=−2√𝜌₁𝜌₂𝜌₁𝜌₂𝑤_𝑏+𝜌₂(1−𝜈). (134)(135)(136)(137) Note that (𝜌₁+𝜌₂+1)²−4𝜌₂(𝜌₁𝑤_𝑏+1)=2√𝜌₁−𝜌₂. Using (133)–(137) we have ℎ(𝑤)=𝐴𝐶+𝐵𝐶−𝐴𝐷𝐶²√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏. (138)

The expansion of ℎ(𝑤) will be in powers of √𝑤_𝑏−𝑤, but only the odd powers will contribute to the asymptotics. Now, 1√2𝜋𝑖𝑤_𝑏−𝑤𝑤^𝑛+1𝑑𝑤=√𝑤_𝑏𝑤_𝑏^−𝑛1√2𝜋𝑖1−𝑈𝑈^𝑛+1√𝑑𝑈=−𝑤_𝑏𝑤_𝑏^−𝑛1()!2𝑛−12𝑛4^𝑛(𝑛!)²1∼−2𝑤_𝑏𝜋𝑛^−3/2𝑤_𝑏^−𝑛,𝑛⟶∞, (139) where we used the binomial expansion of √1−𝑈 and Stirling's formula. It follows that the leading term for 𝑝(𝑛,𝑛) is, in view of (138) and (139), 𝑝()1𝑛,𝑛∼−2𝑤_𝑏𝜋𝐵𝐶−𝐴𝐷𝐶²𝑛^−3/2𝑤_𝑏^−𝑛, (140) with a correction that is 𝑂(𝑛⁻¹) relative to the leading term, which may be computed from the 𝑂((𝑤_𝑏−𝑤)^3/2) term in (138), and a refined Stirling approximation of the factorials in (139). Some of the algebra in our calculations is simplified by introducing 𝑎=𝜌₁+𝜌₂ and 𝑏=√𝜌₁−𝜌₂. Then 𝐶 in (136) factors as (𝑏1𝐶=𝑎−𝑏)(𝑏−1)𝜈−2−12(√𝑎+𝑏+1)=−𝜌₁−𝜌₂−𝜌₁+𝜌₂×1−√𝜌₁−𝜌₂12𝜌₁+𝜌₂+1+√𝜌₁−𝜌₂(1−𝜈), (141) and 𝐶<0 if (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂. After some calculation we find that 𝐵𝐶−𝐴𝐷=2𝜌₁𝜌₂^3/21−√𝜌₁−𝜌₂²𝑤_𝑏×𝜌₁+𝜌₂1𝜈−2−12−√𝜌₁−𝜌₂ (142) so that 𝐵𝐶−𝐴𝐷<0 and then (𝐴𝐷−𝐵𝐶)/(2𝐶²) is the same as the constant 𝐶 in (53). We have thus established (52).

Next we consider 𝑅=𝑚/𝑛>1. The saddle point calculation that led to (27) does not depend on whether (𝜌₁+𝜌₂)²≷𝜌₁−𝜌₂. If the pole at 𝑤_𝑝 is absent, then (27) will apply for all 𝑅∈(1,∞). For 𝑛=𝑂(1) and 𝑚→∞ (25) again holds, since this calculation is independent of whether 𝑤_𝑝 exists or not. However, now a different analysis is needed for 𝑅≈1, and the appropriate scale is to fix 𝑚−𝑛 and let 𝑛→∞. As 𝑅↓1 we have 𝑤_𝑠(𝑅)→𝑤_𝑠(1)=𝑤_𝑏 and then 𝑍₊𝑤_𝑠⟶𝑍₊𝑤_𝑏=𝜌₁+𝜌₂+12𝜌₁. (143) The expansion in (27) breaks down as 𝑅→1 (if (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂) as (30) shows that 𝑀(𝑅)∼2(𝑅−1)²𝜌₁+𝜌₂+1²𝜌₁+𝜌₂+1²−4𝜌₁,𝑅↓1 (144) and thus the factor 1/𝑀(𝑅) in (27) vanishes linearly as 𝑅↓1, which indicates a problem in the asymptotics.

We thus reexamine (10) for 𝑚−𝑛=𝑂(1) and 𝑛 large. We again employ singularity analysis and expand the integrand about the branch point 𝑤_𝑏. We let num_𝑓(𝑤)=𝐴_𝑓+𝐵_𝑓√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏=𝐴_𝑓+𝑂𝑤−𝑤_𝑏, (145) since 𝐵_𝑓=0 in view of (13). Then 𝑓(𝑤)=𝐴_𝑓𝐶−𝐴_𝑓𝐷𝐶²√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏, (146) where 𝐴_𝑓=−2𝜌₁𝜌₂𝑤_𝑏+𝜌²₁(1−𝜈)+𝜌₁𝜌₂−𝜌₂(1−𝜈)+𝜌²₁+𝜌₁𝜈+2𝜌₂−𝜌₁+𝜌₂𝜈√𝜌₁−𝜌₂=(1−𝑏)²1𝜈𝑎−2()=√𝑎+1−𝑏1−𝜌₁−𝜌₂²×1𝜈−2𝜌₁+𝜌₂−12−√𝜌₁−𝜌₂ (147) and we again used some algebraic factorization in the 𝑎,𝑏 variables. From (11), 𝑍₊(𝑤)=[𝜌₁+𝜌₂+1+2√𝜌₁𝜌₂√𝑤_𝑏−𝑤]/(2𝜌₁) and hence 𝑍₊(𝑤)^𝑚−𝑛=𝜌₁+𝜌₂+12𝜌₁^𝑚−𝑛×2√1+𝜌₁𝜌₂𝜌₁+𝜌₂()√+1𝑚−𝑛𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏. (148) By multiplying (146) by (148) and using (139), the leading term for 𝑝(𝑚,𝑛) becomes 𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂×𝑤_𝑏𝜋𝑛^−3/2𝑤_𝑏^−𝑛2𝜌₁𝜌₁+𝜌₂+1^𝑚−𝑛×𝐴_𝑓2√2𝐶𝜌₁𝜌₂𝜌₁+𝜌₂()+𝐷+1𝑚−𝑛𝐶. (149) But (149) is the same as (48), with 𝐶₊=(1−𝜌₁−𝜌₂)𝐴_𝑓/(2𝐶) and Δ₊=𝐷/𝐶, since 𝐷 in (137) may also be written as 𝐷=−2√𝜌₁𝜌₂𝜌₂(1−𝜈)−𝜌₁+14𝜌₁+𝜌₂+1². (150) We have thus established (48).

Now consider 𝑅∈[0,1) with (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂. The integrand in (15) again has a saddle at 𝑤^∗_𝑠(𝑅), and 𝑤^∗_𝑠→0 as 𝑅→0. For 𝑚=𝑂(1) and 𝑛→∞ we again obtain (46). As 𝑅↑1 the saddle 𝑤^∗_𝑠(𝑅)→𝑤^∗_𝑠(1)=𝑤^∗_𝑏>𝑤_𝑏. Thus there is a critical value of 𝑅, which we call 𝑅_𝑑, such that 𝑤^∗_𝑠𝑅_𝑑=𝑤_𝑏. (151) Using (130), and the facts that 𝐹^′_∗(𝑤_𝑠)=0 and (𝜌₁+𝜌₂+1)²−4𝜌₂(𝜌₁𝑤_𝑏+1)=√𝜌₁−𝜌₂, we obtain 𝑅_𝑑𝑤_𝑏𝜌₁+𝜌₂+1+2√𝜌₁−𝜌₂=𝜌₁𝜌₂√𝜌₁−𝜌₂1−𝑅_𝑑, (152) which is a linear equation for 𝑅_𝑑 whose solution is given by (54). If 0<𝑅<𝑅_𝑑 then 𝑤^∗_𝑠(𝑅)<𝑤_𝑏 and the saddle point approximation in (42) holds. If 𝑅_𝑑<𝑅<1 the branch point at 𝑤=𝑤_𝑏 determines the asymptotics of 𝑝(𝑚,𝑛), and we again use singularity analysis. Now we expand 𝑔(𝑤) about 𝑤=𝑤_𝑏, writing (17) as num_𝑔(𝑤)=𝐴_𝑔+𝐵_𝑔√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏, (153) where 𝐴_𝑔=−2𝜌₁𝜌₂𝑤_𝑏+𝜌²₂(1−𝜈)+𝜌₁𝜌₂+𝜌₂()+𝜌1−𝜈₂₁−𝜌₁𝐵𝜈,_𝑔=𝜌₁+𝜌₂𝜈−𝜌₂2√𝜌₁𝜌₂. (154) If we furthermore scale 𝑤−𝑤_𝑏=𝑈/𝑛=𝑂(𝑛⁻¹), then 𝑍^∗₊(𝑤)^𝑚−𝑛=𝑍^∗₊𝑤_𝑏+𝑈𝑛𝑍₊^∗′𝑤_𝑏+𝑂𝑛^{−2𝑚−𝑛}∼𝑍^∗₊𝑤_𝑏^𝑚−𝑛)𝑍exp(𝑅−1^+∗′𝑤_𝑏𝑍^∗₊𝑤_𝑏𝑈,𝑤^−𝑚−1∼𝑤_𝑏^−𝑚−1exp−𝑅𝑈𝑤_𝑏,𝑔(𝑤)=𝐴_𝑔𝐶+𝐵_𝑔𝐶−𝐴_𝑔𝐷𝐶²√−𝑈√𝑛1+𝑂𝑛. (155)(156)(157) It follows, by using (155)–(157) in (15), that 𝑝()∼𝑚,𝑛1−𝜌₁−𝜌₂𝐵_𝑔𝐶−𝐴_𝑔𝐷𝐶²𝑛^−3/2𝑤_𝑏^−𝑚−1𝑍^∗₊𝑤_𝑏^𝑚−𝑛×12𝜋𝑖_Br−−𝑅exp𝑤_𝑏+()𝑍𝑅−1^+∗′𝑤_𝑏𝑍^∗₊𝑤_𝑏𝑈×√−𝑈𝑑𝑈, (158) where Re(𝑈)<0 on Br₋. After some calculation we find that, again using 𝑎=𝜌₁+𝜌₂ and 𝑏=√𝜌₁−𝜌₂, 𝐴_𝑔=(1𝑏+1)(𝑏−1)2+1𝜈−2𝑎𝐵<0,_𝑔𝐶−𝐴_𝑔1𝐷=−4√𝜌₁𝜌₂(1−𝑏)²()×𝑏−2𝜈𝑎+1+𝑎+2𝑏2𝜈²+𝑎+𝑎²+1−2𝑏²=−2√𝜌₁𝜌₂1−√𝜌₁−𝜌₂²×12+√𝜌₁−𝜌₂+𝜌₁+𝜌₂12×1−𝜈4𝜌₁+𝜌₂+1²−𝜌₁(,−𝑅1−𝜈)𝑤_𝑏+()𝑍𝑅−1^+∗′𝑤_𝑏𝑍^∗₊𝑤_𝑏1=−𝑤_𝑏𝑅−𝑅_𝑑1−𝑅_𝑑. (159)(160)(161) The integral in (158) can then be evaluated as 12𝜋𝑖_Br−√−𝑈exp−1𝑤_𝑏𝑅−𝑅_𝑑1−𝑅_𝑑𝑈𝑑𝑈=𝑤_𝑏^3/21−𝑅_𝑑𝑅−𝑅_𝑑^3/2−12√𝜋. (162) Combining (141) with (158)–(162) we obtain precisely the result in (54)–(56).

Now consider the case 𝑅≈𝑅_𝑑, where the branch point at 𝑤=𝑤_𝑏 is close to the saddle at 𝑤=𝑤^∗_𝑠(𝑅). The standard scaling for such situations with coalescing singular points (see [18]) is 𝑤−𝑤_𝑏=𝑂(𝑛^−1/2) and 𝑅−𝑅_𝑑=𝑂(𝑛^−1/2). By expanding in Taylor series about 𝑤=𝑤_𝑏, setting (𝑅−𝑅_𝑑)√𝑛=𝑆 and scaling 𝑆 as in (59), and letting 𝑤−𝑤_𝑏=𝑤_𝑏√𝑛𝑅_𝑑11−𝑅_𝑑+𝜌₁𝜌₂𝑤_𝑏2𝜌₁−𝜌₂^−1/2𝑉, (163) we find that exp−𝑚log𝑤+(𝑚−𝑛)log𝑍^∗₊(𝑤)=𝑤_𝑏^−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛exp−𝑆₁1𝑉+2𝑉²×1+𝑂𝑛^−1/2. (164) Using (163) and (164) in (15) and again expanding 𝑔(𝑤) as in (157) we obtain 𝑝()=𝑚,𝑛1−𝜌₁−𝜌₂𝑀_𝑑^−1/2𝑤_𝑏^−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛𝑛^−1/2×12𝜋𝑖_Br−𝑒^−𝑆1𝑉𝑒^𝑉2/2×𝐴_𝑔𝐶+𝐵_𝑔𝐶−𝐴_𝑔𝐷𝐶²√𝑤_𝑏𝑀_𝑑^−1/4√−𝑉𝑛^1/4+𝑂𝑛^−1/2𝑑𝑉, (165) where Br₋ is to the left of all singularities. Here we also used the definition of 𝑀_𝑑 in (58). The integral(s) in (165) may be easily evaluated using 12𝜋𝑖_Br−𝑒^−𝑆1𝑉𝑒^𝑉2/2𝑒𝑑𝑉=^−𝑆21/2√,12𝜋2𝜋𝑖_Br−𝑒^−𝑆1𝑉𝑒^𝑉2/2√𝑒−𝑉𝑑𝑉=^−𝑆21/4√𝐷2𝜋^1/2−𝑆₁, (166) where the latter follows from (61), after replacing −𝑉 by 𝑉. But 𝑔(𝑤_𝑏)=𝐴_𝑔/𝐶 by (157), and then (165) with (166) give precisely the approximation in (57)–(60). Note that 𝐾 is the same as (𝐵_𝑔𝐶−𝐴_𝑔𝐷)/𝐶².

We can easily compute higher order terms in the expansion, and our analysis shows that the asymptotic series will now involve powers of 𝑛^−1/4. Actually, the leading term in (57), which has a Gaussian dependence on 𝑆₁ (hence on 𝑆), can be obtained by simply expanding the saddle point approximation in (42), for 𝑅↑𝑅_𝑑. However, the 𝑂(𝑛^−1/4) correction term is necessary to see the transition to the range 𝑅>𝑅_𝑑, where (55) applies. In view of (62) the correction term becomes comparable to the leading term when 𝑒^−𝑆21/2∼𝑛^−1/4𝑆₁^−3/2,𝑆₁>0 (167) or, roughly, when |𝑆₁|∼(1/2)log𝑛. Thus for 𝑆₁=𝑂(1) or 𝑆₁≫1 with |𝑆₁|≪(1/2)log𝑛 the Gaussian dominates while for |𝑆₁|≫(1/2)log𝑛 the 𝑂(𝑛^−1/4) term dominates. Note that 𝐾>0 and 𝐷_1/2(−𝑆₁)<0 for 𝑆₁ sufficiently negative, by (62). Thus for 𝑆₁∼(1/2)log𝑛 both the 𝑂(1) and 𝑂(𝑛^−1/4) terms in (57) are positive. The asymptotic matching of (57) for 𝑆₁→+∞ (with 𝑆₁≫(1/2)log𝑛) with (55) as 𝑅↓𝑅_𝑑 is then easily verified, as the approximation in (55) has an 𝑂[(𝑅−𝑅_𝑑)^−3/2] algebraic singularity as 𝑅↓𝑅_𝑑, which agrees with the second formula in (62). Also, up to a multiplicative constant 𝑛^−3/2(𝑅−𝑅_𝑑)^−3/2 agrees with 𝑁^−3/4𝑆₁^−3/2. This completes the derivation of Theorem 3, where (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂.

Next we analyze how the results in Theorem 2 transition to those in Theorem 3, as (𝜌₁+𝜌₂)² decreases through 𝜌₁−𝜌₂(>0). We could simply assume that (𝜌₁+𝜌₂)²=𝜌₁−𝜌₂ and then obtain the necessary asymptotic results. But to see the transition it is necessary to also analyze cases where (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂. To make this more precise we write 𝜌₁+𝜌₂²−𝜌₁−𝜌₂=𝛿𝛽, (168) as in (63), and assume that 𝛿→0⁺. Then 𝛽≷0 according as (𝜌₁+𝜌₂)²≷𝜌₁−𝜌₂. Since only the product 𝛿𝛽 is important, we can set 𝛽=+1,0,−1, according to the cases (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂,(𝜌₁+𝜌₂)²=𝜌₁−𝜌₂,(𝜌₁+𝜌₂)²<𝜌₁−𝜌₂. Then we must relate the small parameter 𝛿 to the large parameters 𝑚,𝑛, and we show below that a natural scaling is to take 𝑚,𝑛=𝑂(𝛿⁻²) as 𝛿→0⁺.

The asymptotic results for 𝑚→∞ and 𝑛=𝑂(1) are the same in Theorems 2 and 3, and thus no transition is needed here. We can write these results in terms of, say, (𝜌₁,𝛿𝛽) rather than (𝜌₁,𝜌₂) and expand for small 𝛿 to somewhat simplify the expression in (25), but we will not do so here. If (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂ the saddle point approximation in (27) applies for all 𝑅>1, while if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ it applies only for 𝑅>𝑅_𝑐. But 𝑅_𝑐→1 (cf. (26)) as (𝜌₁+𝜌₂)²→𝜌₁−𝜌₂ so it will apply for any fixed 𝑅>1. But 𝑅≳1 will require a separate analysis. We also note that the sector 1<𝑅<𝑅_𝑐, where the product form solution in (34) applies, shrinks to zero. Thus if such an approximation will play a role here, it must be contained near 𝑅=1.

We begin by considering the diagonal probabilities 𝑝(𝑛,𝑛), using the scaling in (168). We will approximate ℎ(𝑤) in (19) in such a way that the integrand approaches a nondegenerate limit. The pole at 𝑤_𝑝=(𝜌₁+𝜌₂)⁻², if it exists, is now close to the branch point 𝑤_𝑏=[(𝜌₁+𝜌₂+1)²−4𝜌₁]/(4𝜌₁𝜌₂), since (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂. Note that 𝑎⁻²=[(𝑎+1)²−2(𝑎+𝑏²)]/(𝑎²−𝑏⁴) if 𝑎(=𝜌₁+𝜌₂)=𝑏(=√𝜌₁−𝜌₂). Introducing the scaling 𝑤=𝑤_𝑏(1+𝑈/𝑛), so that 𝑤−𝑤_𝑏=𝑂(𝑛⁻¹), we have 𝑤^−𝑛−1𝑑𝑤∼𝑛⁻¹𝑤_𝑏^−𝑛𝑒^−𝑈𝑑𝑈. Now we write ℎ(𝑤) near 𝑤_𝑏 as ℎ(𝑤)=𝐴+𝑂√𝑤_𝑏−𝑤𝐶+𝐷√𝑤_𝑏−𝑤+𝑂𝑤−𝑤_𝑏. (169) From (134) we find that 1𝐴⟶−2(1−𝑎)²(1+𝑎)as𝑏⟶𝑎or𝛿⟶0⁺, (170) while (136) or (141) shows that (3𝐶∼𝑏−𝑎)(1−𝑎)21𝑎+2−𝑎𝜈∼−𝛿𝛽(),4𝑎1−𝑎)(3𝑎+1−2𝑎𝜈 (171) since 𝑎²−𝑏²=𝛿𝛽 and thus 𝑎−𝑏=𝛿𝛽/(𝑎+𝑏)∼𝛿𝛽/(2𝑎). From (137) or (151) we also find that as 𝑏→𝑎1𝐷⟶−4𝑎(1−𝑎)√1−𝑎²(3𝑎+1−2𝑎𝜈). (172) Then from (169)–(172) we see that 𝐶 becomes comparable to 𝐷√𝑤_𝑏−𝑤 if 𝛿=𝑂(𝑛^−1/2), or 𝑛=𝑂(𝛿⁻²). Then setting 𝑛=𝑦/𝛿² and using (169)–(172), we can approximate (19) to obtain 𝑝()∼()𝑚,𝑛1−𝑎²2()𝛿𝑛1+𝑎𝑤3𝑎+1−2𝑎𝜈^𝑏−𝑛×12𝜋𝑖_Br−𝑒^−𝑈−𝛽/𝑎+√1−𝑎²√−𝑈/√𝑦𝑑𝑈. (173) Thus it is useful to introduce 𝛽₁=𝛽/(𝑎√1−𝑎²) and then we evaluate the contour integral in (173) as, changing 𝑈→−𝑈, 12𝜋𝑖_Br+𝑒^𝑈√𝑦√𝑈−𝛽₁√𝑦=2√𝑑𝑈𝑦2𝜋𝑖_Br+𝑒^𝑍2𝑍−𝛽₁√𝑦+𝛽₁√𝑦𝑍−𝛽₁√𝑦=√𝑑𝑍𝑦√𝜋+√𝑦𝛽𝜋𝑖¹√𝑦_Br+𝑒^𝑍2𝑍−𝛽₁√𝑦=𝑑𝑍𝑦𝜋𝛽1+¹√𝑦√𝑒𝜋𝑖^𝛽21𝑦_Br+𝑒^𝑉2𝑒^{2𝛽1√𝑦𝑉}𝑉=𝑑𝑉𝑦𝜋1+2𝛽₁√𝑦𝑒^{𝛽21𝑦𝛽1√𝑦}_−∞𝑒^−𝑢2𝑑𝑢. (174) Here we used a conformal map 𝑈=𝑍², some contour deformations, (61) with 𝑝=−1, and an identity that relates 𝐷₋₁(⋅) to the standard error function, which yields the last equality in (174). Using (174) in (173) and noting that √𝑦/(𝛿𝑛)=𝑛^−1/2, we obtain (66).

Now consider 𝑅≳1 and we already discussed the case 𝑅>1 with the transition range scaling in (168). We expand now (10) similarly as we expanded (15). For 𝑤−𝑤_𝑏=𝑂(𝑛⁻¹), the factor [𝑍₊(𝑤)]^𝑛−𝑚 can be approximated by 𝑍₊(𝑤)^𝑛−𝑚∼2𝜌₁𝜌₁+𝜌₂+1^𝑚−𝑛2√exp𝜌₁𝜌₂𝑤_𝑏𝜌₁+𝜌₂√+1−𝑈√𝑛(𝑛−𝑚) (175) for 𝑚−𝑛=𝑂(√𝑛). We thus let 𝑧=𝛿(𝑚−𝑛) so that (𝑚−𝑛)/√𝑛=𝑧/√𝑦. Also, 𝑤_𝑏→𝑎⁻² as 𝛿→0⁺, and then 2√𝜌₁𝜌₂𝑤_𝑏𝜌₁+𝜌₂⟶+11−𝑎1+𝑏,𝑏⟶𝑎. (176) Hence (10) becomes 𝑝(𝑚,𝑛)∼1−𝜌₁−𝜌₂𝑤_𝑏^−𝑛2𝜌₁𝜌₁+𝜌₂+1^𝑚−𝑛1𝑛×12𝜋𝑖_Br−𝐴_𝑓𝐶+𝐷√𝑤_𝑏𝑒−𝑈/𝑛^−𝑈×exp−1−𝑎√1+𝑎𝑍−𝑈√𝑦𝑑𝑈. (177) Now 𝐶 and 𝐷 can be evaluated as in (171) and (172), and, as 𝑏→𝑎, 𝐴_𝑓 in (147) becomes 𝐴_𝑓1∼−2(1−𝑎)²(3𝑎+1−2𝜈𝑎). (178) Then (177) can be evaluated similarly to (173), and we ultimately obtain (68), with 𝑧^∼=(1−𝑎)/(1+𝑎)𝑧√𝑦 as in (69).

For 𝑅<1(𝑛>𝑚) we need only consider the ranges 𝑅∈(𝑅⁰_𝑑,1) and 𝑅≈𝑅⁰_𝑑, where 𝑅⁰_𝑑 is in (64), with 𝑅⁰_𝑑 being the limit of both 𝑅^∗_𝑐 and 𝑅_𝑑, as 𝑏→𝑎 or 𝛿→0⁺. For 𝑅>𝑅⁰_𝑑 the saddle at 𝑤^∗_𝑠(𝑅) exceeds the branch point 𝑤_𝑏 and hence the latter determines the asymptotic behavior of 𝑝(𝑚,𝑛). For a fixed 𝑅, we scale 𝑤=𝑤_𝑏+𝑈/𝑛 and use 𝑤^−𝑚−1∼𝑤_𝑏^−𝑚−1−𝑅exp𝑤_𝑏𝑈,𝑍^∗₊(𝑤)^𝑚−𝑛∼𝑍^∗₊𝑤_𝑏^𝑚−𝑛)𝑍exp(𝑅−1^+∗′𝑤_𝑏𝑍^∗₊𝑤_𝑏𝑈. (179)(180) We again expand 𝑔(𝑤) in the form in (169), with 𝐴 now replaced by 𝐴_𝑔, where 𝐴_𝑔1∼−21−𝑎²(1−𝑎+2𝜈𝑎),𝑏⟶𝑎. (181) Using (161) with 𝑅_𝑑 replaced by 𝑅⁰_𝑑, along with (180) and (181), the integral in (7) becomes (1−𝑎)(1+𝑎)(1−𝑎+2𝜈𝑎)23𝑎+1−2𝜈𝑎𝑤𝑛𝛿^{𝑏−𝑚−1}𝑍^∗₊𝑤_𝑏^𝑚−𝑛×12𝜋𝑖_Br+exp𝑉/𝑤_𝑏𝑅−𝑅⁰_𝑑/1−𝑅⁰_𝑑−𝛽/𝑎+√1−𝑎²𝑎√𝑉/√𝑦𝑑𝑉. (182) Here we let 𝑈=−𝑉 and used (171) and (172) to approximate 𝐶 and 𝐷. Scaling 𝛽 as in (178) and evaluating the integral in (182) similarly to (173) leads to (70). If 𝑛−𝑚=𝑂(1), the same analysis applies, as then we can simply replace 𝑅 by 1 and 𝛽₂ by 𝛽₁ in (70) and (71), but must maintain the factor [𝑍^∗₊(𝑤_𝑏)]^𝑚−𝑛.

When 𝑅≈𝑅⁰_𝑑 we let 𝑆⁰=(𝑅−𝑅⁰_𝑑)√𝑛=𝑂(1). Now the saddle 𝑤^∗_𝑠(𝑅) will be close to the branch point 𝑤_𝑏 (with also 𝑤_𝑏∼𝑎⁻² by (168)). We scale 𝑤−𝑤_𝑏=𝑂(𝑛^−1/2) with 𝑤−𝑤_𝑏=−𝑉√𝑛𝑤_𝑏𝑅_𝑑11−𝑅_𝑑+𝜌₁𝜌₂𝑤_𝑏2𝜌₁−𝜌₂^−1/2 (183) and introduce 𝑆⁰₁ as in (73). For 𝑏∼𝑎 we also have 11−𝑅_𝑑+𝜌₁𝜌₂𝑤_𝑏2𝜌₁−𝜌₂∼1+5𝑎+7𝑎²+19𝑎³4𝑎(1+3𝑎) (184) and then 𝑆⁰₁𝑆⁰=1+4𝑎+11𝑎^23/22𝑎()11+3𝑎√1−𝑎²1√1+5𝑎+7𝑎²+19𝑎³. (185) Since 𝐶=𝑂(𝛿) and 𝐷√𝑤_𝑏−𝑤=𝑂(𝑛^−1/4)=𝑂(√𝛿) we have 𝑔(𝑤)=𝐴_𝑔+𝑂𝑤_𝑏−𝑤𝐶+𝐷√𝑤_𝑏−𝑤+𝑂𝑤_𝑏∼𝐴−𝑤_𝑔𝐷1√𝑤_𝑏∼2−𝑤𝑎1+𝑎1−𝑎1−𝑎+2𝜈𝑎11+3𝑎−2𝜈𝑎√𝑤_𝑏−𝑤 (186) and then 𝑝()∼𝐴𝑚,𝑛^𝑔𝐷()11−𝑎𝑛^1/4√𝑤_𝑏𝑅⁰_𝑑11−𝑅⁰_𝑑+1−𝑎²8𝑎^2−1/4×𝑤_𝑏^−𝑚𝑍^∗₊𝑤_𝑏^𝑚−𝑛12𝜋𝑖_Br+𝑒^𝑆01𝑉𝑒^𝑉2/2√𝑉𝑑𝑉. (187) Expressing the integral as a parabolic cylinder function of order 𝑝=−1/2 (see (61)) we obtain the expression in (72), with (73) or (185).

If we let 𝛿→0⁺ but consider even larger values of 𝑛, with 𝑛=𝑦₁/𝛿⁴=𝑂(𝛿⁻⁴), a slightly different expansion applies. Now both 𝐶 and √𝑤_𝑏−𝑤 in (186) become 𝑂(𝛿). Setting 𝑤−𝑤_𝑏=−𝑉/√𝑛, we use 𝑔(𝑤)∼2(1+𝑎)(1−𝑎+2𝜈𝑎)3𝑎+1−2𝜈𝑎−𝛽𝑦₁^1/4𝑎+𝑎√1−𝑎²𝑉⁻¹ (188) in (15) to ultimately obtain the expression in (74), which involves the contour integral in (76). The function ℒ(𝑍;𝛾) can be expressed as an infinite sum of parabolic cylinder functions, as ℒ𝑍;𝛾=^∞_𝑗=0𝛾^𝑗𝐷_{−(𝑗+1)/2}()𝑒−𝑍^−𝑍2/4√2𝜋. (189)

This completes the analysis of the transition range where (𝜌₁+𝜌₂)²≈𝜌₁−𝜌₂, and we have thus established Theorem 4.

5. Light and Heavy Traffic

We establish Theorem 5. First consider the light traffic limit, where 𝜌₁,𝜌₁→0. We now scale 𝑢𝑤=𝜌₁𝜌₂ (190) and then, from (11), (12), (16), and (20), obtain 𝑍₊(𝑤)∼12𝜌₁1+√,𝑍1−4𝑢_∗+(𝑤)∼12𝜌₂1+√,𝑓(𝑤)(𝑤)(𝑤)∼11−4𝑢∼𝑔∼ℎ√.1−4𝑢 (191) Hence the diagonal probabilities have the limiting form 𝑝(𝑛,𝑛)∼𝜌₁𝜌₂^𝑛1𝑢2𝜋𝑖√^−𝑛−1=𝜌1−4𝑢𝑑𝑢₁𝜌₂^𝑛𝑢^𝑛(1−4𝑢)^−1/2=𝜌^𝑛₁𝜌^𝑛₂𝑛,2𝑛 (192) where we used the binomial theorem, and [𝑢^𝑛] denotes the coefficient of 𝑢^𝑛 in the Taylor expansion of what follows.

For 𝑚>𝑛 we similarly obtain 𝑝(𝑚,𝑛)∼𝜌₁𝜌₂^𝑛2𝜌₁^𝑚−𝑛×112𝜋𝑖𝑢^𝑛+111+√1−4𝑢^𝑚−𝑛𝑑𝑢√1−4𝑢=𝜌^𝑚₁𝜌^𝑛₂𝑛,𝑚+𝑛 (193) and this holds also for 𝑚⩽𝑛. By refining (191) with (190) we can obtain higher order terms, as an expansion involving powers of 𝜌₁ and 𝜌₂. We have thus established (78).

For the first heavy traffic limit (HTL1) we let 𝜀=1−𝜌₁−𝜌₂→0⁺, replace 𝜌₂ by 1−𝜌₁−𝜀 and also set 𝑌𝑛=𝜀,𝑚=𝑛+ℓ,𝑤=1−𝜀𝜉. (194) Then after some calculation we find from (14) that (𝑤)=1den4𝜀²𝜉(2+𝜉)2𝜈𝜌₂−𝜈+2−𝜌₂+𝑂𝜀³, (195) and from (21), (13), and (17) that num(𝑤)=𝜀𝜉𝜌₂1−𝜌₂+𝑂𝜀²,num_𝑓(𝑤)=𝜀𝜉𝜌₂2−𝜈−𝜌₂+𝑂𝜀²,num_𝑔(𝑤)=𝜀𝜉1−𝜌₂𝜈+𝜌₂+𝑂𝜀². (196)(197)(198)

Hence the limiting form of (19) is 𝑝(𝑛,𝑛)∼𝜀4𝜌₂1−𝜌₂2𝜈𝜌₂−𝜈+2𝜌₂12𝜋𝑖_Br+𝑒^𝑌𝜉𝜉+2𝑑𝜉=𝜀4𝜌₂1−𝜌₂2𝜈𝜌₂−𝜈+2−𝜌₂𝑒^−2𝑌, (199) which is the same as (80). With the scaling 𝑤=1−𝜀𝜉 the loop integral in (19) can be approximated by the vertical contour Br₊, on which Re(𝜉)>−2 in this case. Since 𝑤^−𝑛−1∼𝑤^−𝑚−1∼𝑒^𝑌𝜉 and 𝑍₊(𝑤)∼1+𝜌₂1−𝜌₂,𝑍^∗₊(𝑤)∼2−𝜌₂𝜌₂ (200) we use (200), (195), (197), and (198) to approximate the integrands in (10) and (15), and this leads to (79) and (81).

We can also derive the HTL1 limits directly from the asymptotic formulas in (34), (35), and (38), as for 𝜌₁=1−𝜌₂−𝜀𝐶₀∼4𝜌₂1−𝜌₂2−𝜌₂+2𝜌₂𝜈𝐶−1𝜀,₊∼4𝜌₂2−𝜈−𝜌₂2−𝜌₂+2𝜌₂𝜈𝐶−1𝜀,₋∼41−𝜌₂𝜌₂+𝜈2−𝜌₂+2𝜌₂𝜈𝜀−1 (201) and (𝜌₁+𝜌₂)^2𝑛=(1−𝜀)^2𝑌/𝜀∼𝑒^−2𝑌. This shows that it is permissible in this case to first let 𝑚,𝑛→∞ with a fixed 𝜌₁,𝜌₂ and then let 𝜌₁+𝜌₂↑1. Note also that in HTL1 the condition (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂ is certainly satisfied, so Theorem 2 applies.

In the second heavy traffic limit, HTL2, we again let 𝜌₁+𝜌₂=1−𝜀 and 𝑤=1−𝜀𝜉, but now 𝜌₂=𝛾𝜀=𝑂(𝜀) is small. We now set 𝑚=𝑋/𝜀 and 𝑛=𝑌/𝜀, and we have 1+𝜌₁+𝜌₂²−4𝜌₁−4𝜌₁𝜌₂𝑤=𝜀²1+2𝛾²+4𝛾𝜉+𝑂𝜀³,𝑍₊(𝑤)1=1+22𝛾+1+2𝛾+1²+4𝛾𝜉𝜀+𝑂𝜀². (202)(203) Then after some calculation we find that as 𝜀→0ℎ(𝑤)⟶ℎ₀(𝜉)=(2−𝜈)𝜉+2𝛾+1(1−𝜈)−(1−𝜈)2𝛾+1²+4𝛾𝜉(2−𝜈)𝜉+2𝛾+1−2𝛾+1²+4𝛾𝜉 (204) and thus 𝑝(𝑛,𝑛) is 𝑂(𝜀²) and given by the expression in (84) (again we have 𝑤^−𝑛−1∼𝑒^𝑌𝜉). Since ℎ₀(𝜉)→1 as 𝜉→∞, the right side of (84) has probability mass along 𝑌=0, and we can write 12𝜋𝑖_Br+𝑒^𝑌𝜉ℎ₀(𝜉)(𝑌)+1𝑑𝜉=𝛿2𝜋𝑖_Br+𝑒^𝑌𝜉2𝛾+1²+4𝛾𝜉−2𝛾+1(2−𝜈)𝜉+2𝛾+1−2𝛾+1²+4𝛾𝜉𝑑𝜉. (205) The presence of the 𝛿(𝑌) term in (205) indicates that different approximations are needed for 𝑝(𝑛,𝑛) for 𝑌>0 and for small 𝑌=𝜀𝑛, such as 𝑛=𝑂(1). But in HTL2 there is little mass along the diagonal 𝑚=𝑛 as a whole.

In HTL2 we also obtain num_𝑓(𝑤)=𝜀²(2−𝜈)𝛾𝜉+𝑂𝜀³,den(𝑤)=𝜀²𝛾()×2−𝜈𝜉+2𝛾+1−2𝛾+1²+4𝛾𝜉+𝑂𝜀³ (206) so that 𝑓(𝑤)→𝑓₀(𝜉), as in (83), and then 𝑝(𝑚,𝑛) is given by (82). We also note that, in view of (203) and the scaling (𝑋,𝑌)=𝜀(𝑚,𝑛), 𝑍₊(𝑤)^𝑛−𝑚⟶exp𝑌−𝑋22𝛾+1+2𝛾+1²+4𝛾𝜉. (207)

If ℱ(𝑋,𝑌) is the limit of 𝜀⁻²𝑝(𝑚,𝑛), for 𝑚>𝑛, in HTL2, then by contour integration ^∞₀^∞_𝑌ℱ(𝑋,𝑌)𝑑𝑋𝑑𝑌≡^∞₀12𝜋𝑖_Br+2𝑓₀(𝜉)𝑒^𝑌𝜉2𝛾+1+2𝛾+1²=1+4𝛾𝜉𝑑𝜉𝑑𝑌2𝜋𝑖_Br−22𝛾+1+2𝛾+1²−1+4𝛾𝜉𝜉𝑓₀(𝜉)=1𝑑𝜉𝑓2𝛾+1⁰(0)=1 (208) so that ℱ is a proper two-dimensional density function, with support in the range 𝑋>𝑌>0. In (208), on Br₋,Re(𝜉)∈(−(2𝛾+1)²/4𝛾,0).

For 𝑚<𝑛 in HTL2 we have 𝑍^∗₊(𝑤)∼2−𝜌₂𝜌₂∼2𝜀𝛾 (209) so that [𝑍^∗₊(𝑤)]^𝑚−𝑛=𝑂(𝜀^𝑛−𝑚). For 𝑛−𝑚=𝑂(1) we also have 𝑤^−𝑚−1∼𝑒^𝑋𝜉∼𝑒^𝑌𝜉 and num_𝑔(𝑤)=𝜀𝜈2𝛾+1²+4𝛾𝜉−2𝛾+1+𝑂𝜀². (210) Then 𝑔(𝑤)∼𝜀⁻¹𝑔(𝜉) as below (86). Thus for 𝑛>𝑚 we have derived (86), and this completes the analysis of HTL2. Note that in HTL2 both 𝑝(𝑛,𝑛) and 𝑝(𝑛−1,𝑛) are 𝑂(𝜀²).

Appendix

Here we study the roots of den(𝑤)=0 in (14). We have den(1)=0 and den^′(1)=−2𝜌₁𝜌₂1−𝜌₁−𝜌₂1−𝜌₂+𝜌₁1−𝜌₁+𝜌₂×1+𝜌₂𝜈+𝜌₁(1−𝜈), (A.1) so that den^′(1)<0 under the stability condition 𝜌₁+𝜌₂<1. Setting Δ₁=𝜌₁+𝜌₂+1²−4𝜌₁𝜌₂,Δ𝑤+1₂=𝜌₁+𝜌₂+1²−4𝜌₂𝜌₁𝑤+1 (A.2) we have den^′(𝑤)=−𝜌₁𝜌₂Δ₁+Δ₂+2𝜌₁𝜌₂𝜌₁𝜌₂𝑤+𝜈Δ₂+𝜌₂𝜌₁𝑤+1−𝜈Δ₁,den^′′(𝑤)=4𝜌²₁𝜌²₂1Δ₁+1Δ₂+𝜌₁𝜌₂𝑤+𝜈Δ₂^3/2+𝜌₂𝜌₁𝑤+1−𝜈Δ₁^3/2. (A.3)(A.4) Now Δ₁=0 when 𝑤=𝑤_𝑏 and Δ₂=0 when 𝑤=𝑤^∗_𝑏>𝑤_𝑏 (since 𝜌₁>𝜌₂). We thus have den(𝑤) real for 𝑤 real and 𝑤∈(−∞,𝑤_𝑏). From (A.4) we conclude that den(𝑤) is a convex function of 𝑤 for 𝑤∈(−∞,𝑤_𝑏). This is obvious for 𝑤⩾0 and for 𝑤<0 we note that Δ₁+𝜌₂(𝜌₁𝑤+1−𝜈)>0 and Δ₂+𝜌₁(𝜌₂𝑤+𝜈)>0, in view of (A.2). Also, (A.3) shows that den^′(𝑤_𝑏)=+∞ and den^′(−∞)=−∞. In view of the convexity and the fact that den(1)=0 and den^′(1)<0 we conclude that den(𝑤)=0 has at most two real roots. But as we discussed in Section 4, 𝑤=𝑤_𝑝=(𝜌₁+𝜌₂)⁻²>1 is a root of den(𝑤)=0 if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂. Setting 𝑤=0 yields den(0)=𝜌₁𝜈𝜌₁+𝜌₂+1−𝜌₁+𝜌₂+1²−4𝜌₂+𝜌₂(1−𝜈)𝜌₁+𝜌₂−+1𝜌₁+𝜌₂+1²−4𝜌₁>0 (A.5) and setting 𝑤=𝑤_𝑏 leads to den𝑤_𝑏=𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(√1−𝜈)−2𝜌₁−𝜌₂𝜌₁𝜌₂𝑤_𝑏+𝜈. (A.6) To determine the sign of den(𝑤_𝑏) we use (49) and the variables 𝑎=𝜌₁+𝜌₂ and 𝑏=√𝜌₁−𝜌₂. Then den(𝑤_𝑏)>0 is equivalent to (𝑎+1)𝑎+(2𝜈−1)𝑏²>(𝑎+1)²+2(𝜈−1)𝑎+𝑏²𝑏 (A.7) which simplifies to (𝑎−𝑏)(1−𝑏)[1+𝑎+2(1−𝜈)𝑏]>0 (A.8) so that den(𝑤_𝑏)>0 according as 𝑎>𝑏 or (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, that is, when 𝑤_𝑝 is a root. It follows that if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂, 𝑤=1 and 𝑤=𝑤_𝑝 are two roots of den(𝑤)=0, and if (𝜌₁+𝜌₂)²<𝜌₁−𝜌₂, 𝑤=1 is the only root. Since den^′(1)<0 and den^′′(𝑤)>0 there can be no roots for 𝑤⩽0.

Now we consider the possibility of having complex roots of den(𝑤)=0. Writing den(𝑤)=𝜌₁+𝜌₂+1𝜌₁𝜈+𝜌₂(1−𝜈)−𝜌₁𝜌₂𝑤+𝜈Δ₂−𝜌₂𝜌₁𝑤+1−𝜈Δ₁ (A.9) we eliminate the radicals to obtain 𝜌₁+𝜌₂+1²𝜌₁𝜈+𝜌₂(1−𝜈)²−𝜌²₁𝜌₂𝑤+𝜈²Δ₂−𝜌²₂𝜌₁𝑤+1−𝜈²Δ₁²=4𝜌²₁𝜌²₂𝜌₂𝑤+𝜈²𝜌₁𝑤+1−𝜈²Δ₁Δ₂. (A.10) Thus if 𝑤 is a root of den(𝑤) it must be a root of (A.10). After some simplification and factoring using MAPLE, we find that solving (A.10) is equivalent to solving the quartic equation (𝑤−1)𝜌₁+𝜌₂²𝑤−1𝐴𝑤²+𝐵𝑤+𝐶=0, (A.11) where 𝐴=𝜌²₁𝜌²₂(2𝜈−1)²,𝐵=𝜌₁𝜌₂𝜈(1−𝜈)1+𝜌₁+𝜌₂²+𝜌₁𝜈+𝜌₂(1−𝜈)1+(2𝜈−1)²,𝐶=𝜈(1−𝜈)𝜌₁+𝜌₂(1−𝜈)1+𝜌₁+𝜌₂²+𝜌₁𝜈+𝜌₂(1−𝜈)²+𝜌₁−𝜌₂²𝜈²(1−𝜈)². (A.12) Here 𝐴, 𝐵, 𝐶 have different meanings from those in Section 4. Clearly 𝐴, 𝐵, 𝐶 are all positive. If 𝜈=1/2, 𝐴=0 and then 𝑤=−𝐶/𝐵<0 is a root (A.11), but it cannot be a root of den(𝑤)=0.

If 𝜈≠1/2, 𝐴𝑤²+𝐵𝑤+𝐶=0 has no roots where 𝑤⩾0, but we can show that it has two real roots in the range 𝑤<0. Letting 𝐵=𝜌₁𝜌₂𝐵₀ and 𝐴=𝜌²₁𝜌²₂𝐴₀ we must show that the discriminant DISC≡𝐵²₀−4𝐴₀𝐶>0. After some calculations we find that DISC=𝜈²(1−𝜈)²1+𝑎²²+8𝑒1+𝑎²+16𝑒²−4(2𝜈−1)²𝑏⁴, (A.13) where 𝑒=𝜌₁𝜈+𝜌₂(1−𝜈). (A.14) But (1+𝑎²)²>4𝑎²=4(𝜌₁+𝜌₂)²⩾4(𝜌₁−𝜌₂)²⩾4(1−2𝜈)²(𝜌₁−𝜌₂)²=4(1−2𝜈)²𝑏⁴, since 𝜈∈(0,1). Thus DISC>0 for 𝜈∈(0,1) and the last factor in (A.11) has two distinct negative roots. If 𝜈=0 or 𝜈=1, there is a double root in the range 𝑤<0, as then DISC=0. In either case these negative roots cannot be roots of den(𝑤).

We have shown that any root of den(𝑤) must be a root of (A.11), and this quartic has only real roots. Then only one 𝑤=1 and 𝑤=(𝜌₁+𝜌₂)⁻² can be roots of den(𝑤), and we already showed that the former is always a root while the latter is a root if and only if (𝜌₁+𝜌₂)²>𝜌₁−𝜌₂. Note that Δ₁𝑤_𝑝=1+𝜌₁+𝜌₂²−4𝜌₁−4𝜌₁𝜌₂𝜌₁+𝜌₂²=1𝜌₁+𝜌₂²𝜌₁+𝜌₂²−𝜌₁−𝜌₂²,Δ₂𝑤_𝑝=1𝜌₁+𝜌₂²𝜌₁+𝜌₂²−𝜌₂−𝜌₁² (A.15) so that for 𝜌₁>𝜌₂, Δ₂(𝑤_𝑝) is always positive, while Δ₁(𝑤_𝑝) is proportional to |(𝜌₁+𝜌₂)²−(𝜌₁−𝜌₂)|, which vanishes along the transition curve.

Acknowledgments

Chalres Knessl work was partly supported by NSA Grant H 98230-11-1-0184. Haishen Yao was supported by PSC-CUNY research award no. 64349-0042.