Asymptotic analysis by the saddle point method of the Anick-Mitra-Sondhi model

We consider a fluid queue where the input process consists of N identical sources that turn on and off at exponential waiting times. The server works at the constant rate c and an on source generates fluid at unit rate. This model was first formulated and analyzed by Anick, Mitra and Sondhi. We obtain an alternate representation of the joint steady state distribution of the buffer content and the number of on sources. This is given as a contour integral that we then analyze for large N. We give detailed asymptotic results for the joint distribution, as well as the associated marginal and conditional distributions. In particular, simple conditional limits laws are obtained. These shows how the buffer content behaves conditioned on the number of active sources and vice versa. Numerical comparisons show that our asymptotic results are very accurate even for N=20.


Introduction
In traditional queueing theory one usually considers a system with one or more servers, at which customers arrive to receive some kind of service. As all servers may be busy, the customers may have to wait for service in a queue. Since there is uncertainty about the actual arrival times of the customers and/or their service requirements, queues are typically modeled stochastically, as random processes.
In the last twenty years, models have appeared in which a continuous quantity, referred to as a fluid, takes the role of the number of customers or queue length. In these models, fluid flows into a reservoir according to some stochastic process. The server may be viewed as a tap at the bottom of the reservoir, allowing fluid to flow out. The rate at which this happens is often constant, but may also be stochastic. Since the fluid reservoir takes the role of the traditional customer queue, it is often referred to as a fluid queue. A third term which is often encountered is fluid buffer, stressing the fact that the storage of fluid is temporary, to prevent loss of fluid at times when the input rate exceeds the output rate. Such models are used as approximations to discrete queuing models of manufacturing systems, high-speed data networks, transportation systems, in the theory of dams, etc. The distribution of the fluid level or buffer content, the average buffer content, the overflow probability (in case of a finite buffer) and the output process are the main quantities of interest in these fluid models.
In the last decade, the literature on queueing theory has paid considerable attention to Markov-modulated fluid models [30], [41]. In these models, a fluid buffer is either filled or depleted, or both, at rates which are determined by the current state of a background Markov process, also called Markovian random environment.

The Markov-modulated fluid model
In this section we describe a general model of fluid entering and leaving a single buffer system [24]. Let X(t) denote the amount of fluid at time t in the buffer. Furthermore, let Z(t) be a continuous-time Markov process. Z(t) will be said to evolve "in the background". We will assume that Z(t) has a finite state space N with N = {0, 1, . . . , N}.
The buffer content X(t) is regulated (or driven) by Z(t) in such a way that the net input rate into the buffer (i.e., the rate of change of its content) is d [Z(t)]. The function d (·) is called the drift function. When the buffer capacity is infinite, the dynamics of X(t) are given by The condition at X(t) = 0 ensures that the process X(t) does not become negative. When the buffer capacity is finite, say K, the dynamics are given by The condition at X(t) = K prevents the buffer content from exceeding K.
One of the main reasons why Markov-modulated fluid models have attracted so much attention is that they are relevant for modelling telecommunications networks. In the literature, more attention has been paid to models in which the buffer capacity is infinitely large, because this case is easier to analyze. For most practical situations in telecommunications the infinite buffer case is a good approximation for the finite buffer case, since overflow of the buffer is assumed to be extremely rare. For highly loss-sensitive traffic, allowed loss fractions in the order of 10 −6 to 10 −10 are typical [9].
In many applications, the process Z(t) evolves as a finite state birth-death process. Birthdeath processes have wide applications in many practical problems [13]. They can be regarded as continuous time analogs of random walks.
The paper which has become the main reference for birth-death fluid models was written by Anick, Mitra and Sondhi [3]. We shall refer to their model as the AMS model. The AMS model describes an infinitely large fluid buffer which is fed by N identical exponential on-off sources and emptied by an output channel with constant capacity c. An on source turns off at rate 1 and an off source turns on at rate λ. Thus, the net input is regulated by a specific birth-death process Z(t) with state space N , birth rate λ k = λ (N − k) , death rate µ k = k and drift function d(k) = k − c. The rates are conditioned on Z(t) = k, 0 ≤ k ≤ N.
In a way, the AMS model was a generalization of the earlier model of Kosten [19], [20], where the limiting case N → ∞, λ → 0 (with λN fixed) is considered. According to the authors in [3], their model and its variants had also been proposed in other previous papers [7], [15], [40].
In [21], Kosten generalizes the AMS model by considering the same problem for a Multi Group Finite Source (MGFS) system, consisting of m groups of identical sources and a single shared buffer. In [22] Kosten considers a central processor, working at uniform speed, that receives information at uniform rates from a number of sources and a buffer which stores the information that cannot be handled directly. He shows that the stationary probability of overflow G(x) satisfies G(x) ∼ Ce −αx as x → ∞. For fairly complex aggregates of sources, a method is developed (called the decomposition method) to determine the exponent α. The coefficient C is determined by simulation.
Daigle and Landford [10] used the AMS model to study packet voice communication systems. Tucker [50] also used a similar model, but the server capacity was defined as an integer number of information units and the model had finite buffer size. Results on the buffer content distribution, cell loss and delay were given, along with simulation comparisons.
The feature common to the above models is that they all assume a continuous time distribution for the on-off intervals. On the other hand, Li [26] introduced a discrete time model with a finite number of on-off sources and geometric distributions for the on-off intervals. He assumed that in one time unit only one on-off source can change state. The channel capacity was assumed to be an integer number of sources and the buffer size could be either zero (burst switching-clipping case) or infinity (packet switching case).
In [31] Mitra shows how to deal with states i ∈ N for which d(i) = 0. He also generalizes the AMS model in two ways. Again, he considers a buffer which receives input from N i.i.d. exponential sources, called producers. However, now the output is also assumed to be Markov-modulated, by M i.i.d. exponential consumers. A second generalization is that he considers the case of finite buffer capacity.
In the manufacturing literature Wijngaard [55] has analyzed the model treated in [31] In this paper we analyze the joint distribution of the buffer content and number of on sources in the AMS model asymptotically, as the number of sources N → ∞. Our results extend those obtained in [32].
In Section 2 we describe the AMS model in more detail and in Section 3 we derive an integral representation of the solution, from which the spectral representation obtained by the authors in [3] can be easily derived. In Sections 4-9 we obtain asymptotic expansions from the exact solution by using the saddle point method. We have to consider eleven relevant regions of the two-dimensional state space. In Section 10 we summarize our results and numerically compare our approximations with the exact solution. In Section 11 we analyze the conditional and marginal distributions.

The Anick-Mitra-Sondhi (AMS) model
In this model a data-handling switch receives messages from N mutually independent sources, which independently and asynchronously alternate between on and off states. The on and off periods are exponentially distributed for each source. These two distributions, while not necessarily identical, are common to all sources.
The unit of time is selected to be the average on period and the average off period is denoted by 1 λ . Thus, if we denote by On (Off) the random variables measuring the elapsed time a source is in the on (off) state, we have An on source will generate packets of information at a uniform rate that we assume, without loss of generality, to be 1 packet per unit of time. While a source is off, it generates no packets. Thus, when k sources are on simultaneously, the instantaneous receiving rate at the switch is k. The switch stores or buffers the incoming information that is in excess of the maximum transmission rate c of an output channel (thus, c is also the ratio of the output channel capacity to an on source's transmission rate). We assume that c is not an integer. Let Z(t) be the number of on sources at time t, and X(t) be the buffer content. As long as the buffer is not empty, the instantaneous rate of change of the buffer content is Z(t) − c. Once the buffer is empty, it remains so as long as Z(t) < c. We assume that the buffer capacity is infinite, so we need the following stability condition to ensure the existence of a steady state distribution. Note that c < N implies that the buffer may be non-empty and c > λ λ+1 N means that the output rate exceeds the average input rate.
Here we consider only the equilibrium probabilities, defined by The balance equations are [3] (k − c) Moreover, if the number of sources on at any time exceeds c, then the buffer content increases and the buffer cannot stay empty. Hence, Also, since F k (∞) is the probability that k out of N sources are on simultaneously and from (2) this follows a binomial distribution with parameter λ λ+1 . In [3] Anick, Mitra and Sondhi derived the spectral representation of the solution to (4)- (7). They wrote the equation (4) in matrix form and M is an (N + 1) × (N + 1) tridiagonal matrix with entries Using (8), they represented the solution F(x) as where z i is a negative eigenvalue and φ i is the associated right eigenvector of the matrix They showed that (4) and (6) imply conditions on the derivatives of F N (x) at x = 0, From (10) they found the coefficients a i in (9) with the help of Vandermonde determinants. Asymptotic and numerical results for the overflow probability of a predefined buffer backlog were also presented.

Exact solution
Here we shall solve the problem (4)-(7) by a different method than that used in [3]. This will lead to an alternate form of the solution that shall prove useful for asymptotic analysis.
Looking for a solution to (4) of the form we get the following difference equation for h k (θ) To solve (11) we represent h k (θ) as an integral where C is an small loop around the origin in the complex plane and H(θ, w) is an analytic function of w inside C. Using (12) in (11) and integrating by parts gives an ODE for H(θ, w) whose solution is with Using the binomial theorem we can also represent h k (θ) as The polynomials h k (θ) are related to the Krawtchouk polynomials [47]. The functions c k (θ) are entire (polynomial) functions of θ, that may be computed from (17).
Since H is analytic inside C, we have Imposing the boundary condition (5) yields h k (θ) = 0 for k ≤ −1. Thus, and therefore H(θ, w) must be a polynomial in w of degree ≤ N. In view of (14) we must have which gives the equation for the eigenvalues θ j . Since from (7) we need F k (x) → F k (∞) for x → ∞, we consider only the negative eigenvalues. Solving for θ j in (18) we get where and with 0 < ρ < 1 from (3). Combining the above results and using (7), we find the spectral representation of the solution where the constants a j are to be determined. To find a j we employ arguments different from those in [3].
Taking the Laplace transform in (4) and using the boundary conditions (6) for k ≥ ⌊c⌋ + 1 gives Equation (23) has two independent solutions, h k (ϑ) and h * k (ϑ), with h k (ϑ) given by (12) and where C is a small loop around the origin and Using the binomial theorem we get from (24) From (25) we have where the c * k (ϑ) are polynomials that can be identified from (25) using (15) and (16). Therefore, the general solution to (23) is given by We note that h k and h * k are entire functions of ϑ. Since h N +1 (ϑ) = 0 but h * N +1 (ϑ) = 0, we see from (5) that we need A 2 (ϑ) = 0, so that Taking the Laplace transform in (22) yields from which we see that the only singularities of F k (ϑ) are simple poles at ϑ = 0 and ϑ = θ j , 0 ≤ j ≤ N − ⌊c⌋ − 1. Since h k (ϑ) is entire, A 1 (ϑ) must have poles at these points also. We thus rewrite (26) as where B(ϑ) is an entire function. Using the inversion formula, we can represent F k (x) in the form where Br + is a vertical contour on which Re(ϑ) > 0. From (27) the residue of F k (ϑ) at ϑ = 0 is equal to and from (14) we know that Thus, from (28), we get We determine B(ϑ) by examining the limit ϑ → ∞, which correspond to considering the boundary condition at x = 0.
From (15) and (16) we can easily obtain the asymptotic expressions Setting w = u ϑ in (12) and using (31) yields We also have Combining (32) and (33) we get Since (6) implies that lim Setting k = ⌊c⌋ + 1 we see that B(ϑ) is an entire function that is o (ϑ) as ϑ → ∞. By the generalized Liouville theorem [17], B(ϑ) must be a constant Therefore we can write Closing Br + in the half-plane Re(ϑ) ≤ 0 we recover the spectral representation In conclusion, we have proved the following theorem.
Theorem 1 The solution of (4)- (7) is given by where Br + is a vertical contour on which Re(ϑ) > 0, C is a small loop around the origin in the complex w−plane We have thus re-derived, using different arguments, the spectral representation obtained by the authors in [3]. We also derived the integral representation (36) that will prove more useful than (37) for asymptotic analysis. We note that the coefficients in the sum in (37) alternate in sign and this makes it very difficult to obtain asymptotic results for F k (x) from the spectral representation.
We study asymptotic properties of F k (x) as N → ∞ for various ranges of k and x. We shall use (36), approximate the integrand for N large and then evaluate the integral over Br + asymptotically. We will introduce the scaled variables (y, z) = x N , k N and consider the limit N → ∞, with y and z fixed. In Section 4 we will derive asymptotic approximations to

The main approximation
We first evaluate the product in (36) for N → ∞.
Proof. Since we assumed that c is not an integer, we introduce the fractional part of c defined by α = c − ⌊c⌋ , 0 < α < 1 (41) which allows us to write N − ⌊c⌋ as (1 − γ)N + α. Since we can apply the Euler-Maclaurin formula [5] to get Changing variables from k to x = k N we obtain If we denote the indefinite integral in (42) by where Therefore, To find I 1 (z), we use integration by parts Since (18) implies that V [−σ(x)] = x, we can write (48) as or, changing the integrating variable from x to σ, The first integral in (49) is The second integral in (49) is quite complicated and after some calculation we find that Hence, with and taking (47) into account we get Using (45) in (51) we have From (45) and (50) we get Combining these results in (44) yields From (45) we see that as N → ∞ Therefore, we conclude that Exponentiating this yields the formula for P (ϑ; N). We note from (39) and consider the asymptotic approximation to F k (x) inside the domain Writing h k (ϑ) in terms of z we have with We first assume that γ < z < 1 and compute the integral (54) as N → ∞ by the saddle point method [6]. Thus, we locate the saddle points W (ϑ, z) by solving the equation and find two saddle points W + (ϑ, z) and W − (ϑ, z) defined by We next study the motion of the saddle points W ± (ϑ, z) , and also the branch points − 1 R 2 (ϑ) and 1 R 1 (ϑ) of η(w, ϑ, z), as ϑ decreases from ∞ to θ 0 , with z > γ. First we observe from (16) that for all real ϑ, R 1 (ϑ) and R 2 (ϑ) > 0 and thus .
For large ϑ we have When ϑ = 0 the branch point 1 and we get For θ 0 < ϑ < 0, the saddle point W + moves to the right of the branch point 1 Finally, when ϑ = θ 0 with γ < z < γ 2 δ the branch point 1 Note that when ϑ = ϑ * , where the discriminant D(ϑ, z) vanishes and the saddle points coalesce. Thus, Hence, D(ϑ, z) has a double zero at ϑ = ϑ * . From (59)-(62) we conclude that for ϑ > θ 0 with z > γ we must deform the contour C to the steepest descent contour through the saddle W − and therefore Note that the steepest descent directions at w = W − are arg (w − W − ) = ± π 2 and we can shift C into another circle that goes through this saddle point.
Using (55) and (56) we find that To compute the above integral as N → ∞, we shall again use the saddle point method. We and Θ(y, z) is defined to be the solution of the equation where we have used ∂η ∂w [W − (Θ, z) , Θ, z] = 0.
For each value of Θ, (70) defines implicitly a curve in the (y, z)-plane. Since (67) has a pole at ϑ = 0, (68) is only valid if y and z are such that Θ > 0, i.e., the saddle lies to the right of the pole. Then Br + can be deformed to the steepest descent contour through ϑ = Θ. We also note that when Θ = 0, (68) is singular. We shall find an asymptotic approximation valid when the saddle is close to ϑ = 0 in Section 6. From (70) we see that Θ = 0 corresponds to the curve For y > Y 0 (z) the saddle point Θ lies on the negative real axis, to the left of the pole at ϑ = 0. We re-write (36) as where θ 0 < Re(ϑ) < 0 on Br − . Note that the residue at ϑ = 0 was already calculated in (29) and corresponds to F k (∞). Using (38) and (65) in (72) and approximating the integral by the saddle point method, we obtain with Ψ(y, w, ϑ, z) and Θ(y, z) defined by (69) and (70) respectively. From (70) we find that the value Θ = θ 0 corresponds to the curve Our approximation to the integral in (36) assumed that ϑ > θ 0 . However, (73) passes smoothly through the curve y = Y 1 (z) since from (66) we have The discriminant D(Θ, z) vanishes when where ϑ * (z) was defined in (63). But, as we showed in (64), the derivative of the discriminant also vanishes along the curve y = Y * (z) so that (73) can be smoothly continued through this curve by simply replacing W − by W + . We summarize our results below.
and Θ + (y, z) is defined to be the solution of the equation The three regions R j are sketched in Figure 1. Note that as we pass from R 1 to R 2 there is a "phase transition" in the asymptotics. However, the solution is smooth as we pass from R 2 to R 3 .
Theorem 3 ceases to be valid along the curve y = Y 0 (z), γ ≤ z ≤ 1, because the function G 1 (y, z) is singular there. We shall find a corner layer solution near the point (0, γ) in Section 5, a transition layer solution along y = Y 0 (z), γ < z < 1 in Section 6 and a corner layer solution near the point (Y 0 (1), 1) in Section 8.2.

From (57) we have
and from (55) we get Hence, and our approximation develops singularities as z → 0 and z → 1, which corresponds to k ≈ 0 and k ≈ N. We shall find a boundary layer solution near z = 0 in Section 7 and a boundary layer solution near z = 1 in Section 8. In these cases we will need to re-examine the expansion of (54). Finally, from (70) we see that Θ → ∞ as y → 0 if γ < z < 1. We shall find a boundary layer approximation near y = 0 with z > γ in Section 9. To analyze the case of small y we will need an expansion of P (ϑ) in Lemma 2 valid for ϑ large. Note that Theorem 3 does apply for y → 0 as long as 0 < z < γ. In this range we can show that Θ + in (83) remains finite. In particular, the boundary masses F k (0) for k = Nz with 0 < z < γ can be computed by setting y = 0 in (80). We also note that in this range F k (0) ∼ F k (∞) and the exponentially small (as N → ∞) difference can be estimated from (80).
Proof. We write Changing index from j to n = (1 − γ)N + α − 1 − j and defining we have We evaluate individually the three products in (88) beginning with P 1 (S). Setting ε = 1 N and using the Euler-MacLaurin summation formula, we have

Changing variables to
Defining From (20) we have and evaluating this at x = 1 − γ we get After some calculations, we find that From (89) we obtain I 2 (0) = 2γ ln(2) + 1 2 (γ + 1) ln(λ + 1) + γ 2 ln (δ) + γ ln(γ) and I 2 (1 − γ) = 2γ ln(2) + 1 2 (γ + 1) ln(λ + 1) + γ 2 ln (δ) + γ − 1. Thus, and we have From (20) and (87) we get Therefore we conclude that We next consider the second product in (88). Using properties of Gamma functions [1] we have S + 1 and Stirling's formula [54] gives Finally we determine the asymptotic approximation to P 3 (S). The Euler-MacLaurin summation formula yields Changing variables in the integral in (92) from n to We write (93) as the difference of two integrals The first integral in (94) is easily evaluated as ε → 0 to give Introducing the function in (94) and using (15) we have After some work we find that From (95) we have, as ε → 0 We conclude that Finally we compute the remaining part of (92) Hence, Combining (90), (91) and (96) we obtain as in (86). Next, we shall analyze the function h k (ϑ) for N → ∞ and (y, z) close to the corner (0, γ). To do so, we introduce the new variable l, defined by and the new function h (1) (l, S) = h l+c−α (SN). From (12) and (55) we have Using (31) we see that as N → ∞ the branch points are located at and that, for N → ∞ with S and w fixed, Therefore, we replace the contour C in (98) by a Hankel's loop and use (99) to obtain where H begins at w = −∞, encircles w = 0 once in a positive sense and returns to its starting point. We assume that the branch of w φ S +l−α−1 takes its principal value at the point where the contour crosses the positive real axis and is continuous everywhere, except for w real and negative [34].
Introducing the new variable t = β 2γ 1 w , with we get where H ′ is another Hankel contour of the same shape as H. Using an integral representation of the Bessel function J · (·) [52] we obtain Using (86) and (101) we obtain the following theorem.
Br + is a vertical contour in the complex plane with Re(S) > 0, J · (·) is the Bessel function, Γ (·) is the Gamma function , and 6 The transition layer at y = Y 0 (z) (Region II) As we observed previously in Section 4, the approximation (68) is not valid along the curve y = Y 0 (z) given by The condition y ≈ Y 0 (z) corresponds to the saddle Θ being close to the pole at ϑ = 0 in (36). To obtain the expansion of (68) for y ≈ Y 0 (z), we first approximate the integrand in (67) for small ϑ, using Using the above in (67) yields with Φ(z) given by (103). Changing variables to we obtain [14] 1 2πi where Br ′ is another vertical contour, on which Re(v) > 0. Therefore, This is valid for y − Y 0 (z) = O N − 1 2 with γ < z < 1 and gives the transition between regions R 1 and R 2 in Theorem 3.

The boundary z = 0 (Region III)
We shall next analyze the part of ∂D corresponding to z = 0. We start by proving the following lemma.
Proof. From (14) we have for N → ∞, with w = O(N), We replace C by C ′ = C N and change variables from w to u = w N in (12) to get Using (109) we have which evaluates to give (108).

The boundary z = 1
We examine the scale k = N − O(1), which corresponds to 1 − z = O(N −1 ). The analysis is different for three ranges of y.

8.1
The boundary layer at z = 1, 0 < y < Y 0 (1) (Region IV) Proof. From (12) we have Changing variables to u = wN we get where C ′ = C N is a small loop around u = 0. Hence, the result follows by expanding exp 1+(1−γ)ϑ λ u in Taylor series.
Using (38) and (112) in (36) we have We evaluate the integral above for N → ∞ by the saddle point method and obtain 8.2 The corner layer at (Y 0 (1), 1) (Region V) From (114) we see that Θ 1 (y) → 0 as y → Y 0 (1). Therefore, we shall find a solution valid in the neighborhood of the point and from (39) we get Using (38) and the above in (35) we obtain 8. 3 The boundary layer at z = 1, y > Y 0 (1) (Region VI) We shall examine the portion of the boundary z = 1 for y > Y 0 (1). To do so, we write F k (x) as and, repeating the same calculation done in Section 8.1, we conclude that 9 The boundary y = 0 In this Section we analyze the part of the boundary where the boundary condition (6) applies. Theorem 3 is not valid for small y since for γ < z ≤ 1 the saddle Θ(y, z) → ∞ as y → 0 + , but our approximation to the integrand in (36) was only valid for ϑ = O(1). We will scale y = x/N and use Lemma 4, which applies for ϑ = SN = O(N). We note that yϑ = xS. We will need to examine separately the cases γ < z < 1 and 1 − z = O(ε).

9.1
The boundary layer at y = 0, γ < z < 1 (Region VII) Proof. Changing variables in (12) from w to u = wSN, we have Use of the saddle point method gives, for N → ∞ with S and z fixed, where the saddle point is located at Using (86) and (118) in (36), we get where Br ′ is a vertical contour on which Re(S) > 0. Computing the integral in (119) as N → ∞ by the saddle point method, we find that the integrand has a saddle at and hence the leading term for (119) is After some simplification we get Note that (121) is singular as z → γ. We can show that if we expand (121) as (x, z) → (0, γ), with x z−γ fixed, the result asymptotically matches to he corner layer approximation in Theorem 5. The latter must be expanded for χ, l → ∞ with l χ fixed.

The corner layer at (0, 1) (Region VIII)
When z → 1, (121) becomes singular. Therefore, we shall find an approximation to F k (x) close to the point (0, 1). We first observe that when we obtain from (12), by a calculation similar to that in the proof of Lemma 7, Using (86) and (122) in (35) we have Using the saddle point method we find that the integrand of (123) has a saddle at where S * (x, z) was defined in (120). Hence, to leading order, we obtain Note that from (121) and (124) we have and thus we recover the result found in [3] (equation (38) therein), that d n F k dx n (0) = 0, ⌊c⌋ + 1 + n ≤ k ≤ N.

Summary and numerical studies
In most of the strip D = {(y, z) : y ≥ 0, 0 ≤ z ≤ 1} , the asymptotic expansion of F k (x) is given by Theorem 3. Below we summarize our results for the various boundary, corner and transition layer regions, where Theorem 3 does not apply. The paragraph number refers to the corresponding region (see Figure 2).

Theorem 9
As N → ∞, asymptotic approximations of the joint distribution F k (x) are as follows. Br + is a vertical contour in the complex plane with Re(S) > 0, J · (·) is the Bessel function, Γ (·) is the Gamma function , and and Θ 0 (y) is the solution to the equation .
where F (4) j (y) is as in item IV.
From (36) we see that the continuous part of the density (for x > 0) is given by Therefore, we can approximate the density by multiplying the asymptotic approximations by the corresponding saddle point: Here the right sides are the expansions valid in the various regions. We compare our asymptotic approximations to the density with the exact results obtained by evaluating (136). We use N = 20, γ = 0.37987897 and λ = 0.0122448. In Figure 3 we graph f k (x) and the asymptotic approximation ΘG 1 (y, z), for k = 17. We note that the asymptotics predict that the maximum value will be achieved at x = NY 0 k N ≃ 3.052. In Figure 4 we graph f 3 (x) and the asymptotic approximation Θ + G 2 (y, z). We note that the maximum is now achieved at x = 0, which is consistent with the asymptotic analysis.
In Figure 5 we graph f 0 (x) and Θ 0 F . This shows the importance of treating Region III.
In Figure 6 we graph ln [f 20 (x)] and ln Θ 1 F (4) 0 (y) , the latter being the boundary layer approximation to the density at z = 1.
In Figure 7 we graph ln [f k (x)] and ln S * (x, z)NF (7) (x, z) with x = 0.001. We see how the asymptotic approximation breaks down for z ≈ γ and for z ≈ 1 (it does not apply for z < γ).   In 8 we graph f k (x) and the asymptotic approximations ΘG 1 (y, z) and Θ + G 2 (y, z) with x = 5. To show in more detail the different ranges, in Table 1 we compare the exact value of the density with the approximations given by Θ 0 F (3) k (y) and Θ + G 2 (y, z), for x = 1 (y = 0.05) with 0 ≤ k ≤ 9. In Table 2 we compare f k (x) with the approximations given by Θ 1 F (4) j (y) and ΘG 1 (y, z), for x = 1 (y = 0.05) with 10 ≤ k ≤ N. We note that for our choice of parameters c ≃ 7.598 and Y 0 (1) ≃ .4998. Also, when y = 1/N = 0.05, the equation Y * (z) = y has the solution z ≃ .4811, so that the transition between R 2 and R 3 occurs at Nz ≃ 9.622.

Conditional limit laws
In this section we will study the probability of k sources being on, given that the buffer content is equal to x > 0 and the probability of the buffer having some non-zero content x, given that k sources are on, for various ranges of x and k.

The distribution of active sources conditioned on x
We define and the conditional probabilities We first consider x = O(N). In this range we compute M 1 (x) using the approximation in R 3 exp NΨ(y, W + , Θ + , z) dz.
Next we assume a small buffer size, with x = O(1/N). Again the main contribution to the sum (137) comes from k ≈ c, but for x = O(1/N) we must use the Region I approximation f k (x) ∼ NF (1)′ l (χ). Thus, which we can write as a residue series where we have used Using the generating function [52] ∞ j=−∞ J j (x)z j = exp and thus It follows that which is a discrete distribution in l.

The buffer density conditioned on k
We first consider γ < z ≤ 1. In this range F k (0) = 0 and from (141) we have M 2 (k) = F k (∞).
This shows that if the number of on sources is less than the service rate c and the buffer is not empty, the conditional density is approximately exponential in x, with mean 1 −Θ + (0,z) . This exponential shape is also illustrated in Figure 4.
Finally, we consider the case k ≈ c. We shall again use the approximation in Region I and from (142) we get which is a density in χ, consisting of an infinite mixture of exponentials.