Dominant poles and tail asymptotics in the critical Gaussian many-sources regime

The dominant pole approximation (DPA) is a classical analytic method to obtain from a generating function asymptotic estimates for its underlying coefficients. We apply DPA to a discrete queue in a critical many-sources regime, in order to obtain tail asymptotics for the stationary queue length. As it turns out, this regime leads to a clustering of the poles of the generating function, which renders the classical DPA useless, since the dominant pole is not sufficiently dominant. To resolve this, we design a new DPA method, which might also find application in other areas of mathematics, like combinatorics, particularly when Gaussian scalings related to the central limit theorem are involved.


Introduction
Probability generating functions (PGFs) encode the distributions of discrete random variables. When PGFs are considered analytic objects, their singularities or poles contain crucial information about the underlying distributions. Asymptotic expressions for the tail distributions, related to large-deviations events, can typically be obtained in terms of the so-called dominant singularities, or dominant poles. The dominant pole approximation (DPA) for the tail distribution is then derived from the partial fraction expansion of the PGF and maintaining of this expansion the dominant fraction related to the dominant pole. Dominant pole approximations have been applied in many branches of mathematics, including analytic combinatorics [7] and queueing theory [19]. We apply DPA to a discrete queue that has an explicit expression for the PGF of the stationary queue length. Additionally, this queue is considered in a many-sources regime, a heavy-traffic regime in which both the demand on and the capacity of the systems grow large, while their ratio approaches one. This many-sources regime combines high system utilization and short delays, due to economies of scale. The regime is similar in flavor as the QED (quality and efficiency driven) regime for many-server systems [8], although an important difference is that our discrete queue fed by many sources falls in the class of single-server systems and therefore leads to a manageable closed form expression for the PGF of the stationary queue length Q. Denote this PGF by Q(z) = E(z Q ).
PGFs can be represented as power series around z = 0 with nonnegative coefficients (related to the probabilities). We assume that the radius of convergence of Q(z) is larger than one (in which case all moments of Q exist). This radius of convergence is in fact determined by the dominant singularity Z 0 , the singularity in |z| > 1 closest to the origin. For PGFs, due to Pringsheim's theorem [7], Z 0 is always a positive real number larger than one. Then DPA leads to the approximation with c 0 = lim z→Z 0 (z − Z 0 )Q(z). In many cases the approximation (1.1) can be turned into a more rigorous asymptotic expansion (for N large) for the tail probabilities P(Q > N ). We shall now explain in more detail the many-sources regime, the discrete queue, and when combining both, the mathematical challenges that arise when applying DPA.
Many sources and a discrete queue. Consider a stochastic system in which demand per period is given by some random variable A, with mean µ A and variance σ 2 A . For systems facing large demand one can set the capacity according to the rule s = µ A + βσ A , which consists of a minimally required part µ A and a variability hedge βσ A . Such a rule can lead to economies of scale, as we will now describe in terms of a setting in which the demand per period is generated by many sources. Consider a system serving n independent sources and let X denote the generic random variable that describes the demand per source per period, with mean µ and variance σ 2 . Denote the service capacity by s n , so that the system utilization is given by ρ n = nµ/s n , where the index n expresses the dependence on the scale at which the system operates. The traditional capacity sizing rule would then be with β some positive constant. The standard heavy-traffic paradigm [11,15,16], which builds on the Central Limit Theorem, then prescribes to consider a sequence of systems indexed by n with associated loads ρ n such that (also using that s = s n ∼ nµ) where γ = βσ/ √ µ. We shall apply the many-sources regime given by (1.2) and (1.3) to a discrete queue, in which we divide time into periods of equal length, and model the net input in consecutive periods as i.d.d. samples from the distribution of A, with mean nµ and variance nσ 2 . The capacity per period s n is fixed and integer valued. The scaling rule in (1.3) thus specifies how the mean and variance of the demand per period, and simultaneously s n , will all grow to infinity as functions of n. Many-sources scaling became popular through the Anick-Mitra-Sondhi model [1], as one of the canonical models for modern telecommunications networks, in which a switch may have hundreds of different input flows. But apart from communication networks, the concept of many sources can apply to any service system in which demand can be regarded as coming from many different inputs (see e.g. [4,17,6,13,18] for specific applications. How to adapt classical DPA? As it turns out, the many-sources regime changes drastically the nature of the DPA. While the queue is pushed into the many-sources regime for letting n → ∞, the dominant pole becomes barely dominant, in the sense that all the other poles (the dominated ones) of the PGF are approaching the dominant pole. For the partial fraction expansion of the PGF this means that it becomes hard, or impossible even, to simply discard the contributions of the fractions corresponding to what we call dominated poles: all poles other than the dominant pole. Moreover, the dominant pole itself approaches 1 according to This implies that in (1.1) the factor c 0 /(1 − Z 0 ) potentially explodes, while without imposing further conditions on N , the factor Z −N −1 0 goes to the degenerate value 1. The many-sources regime thus has a fascinating effect on the location of the poles that renders a standard DPA useless for multiple reasons. We shall therefore adapt the DPA in order to make it suitable to deal with the complications that arise in the many-sources regime, with the goal to again obtain an asymptotic expansion for the tail distribution. First observe that the term Z −N −1 0 in (1.1) becomes non-degenerate when we impose that N ∼ K √ nσ 2 , with K some positive constant, in which case The condition N ∼ K √ nσ 2 is natural, because the fluctuations of our stochastic system are of the order √ nσ 2 . Of course, there are many ways in which N and n can be coupled, but due to (1.4), only couplings for which N is proportional to √ n lead to a nondegenerate limit for Z −N −1

0
. Now let us turn to the other two remaining issues: The fact that c 0 /(1 − Z 0 ) potentially explodes and that the dominated poles converge to the dominant pole.
To resolve these two issues we present in this paper an approach that relies on approximations of the type (1.4) for all the poles (which are defined implicitly as the solutions to some equation). The approximations are accurate in the many-sources regime, and can then be substituted into the partial fraction expansion that describes the tail distribution. We replace the partial fraction expansion by a contour integral representation, and subsequently apply a dedicated saddle point method recently introduced in [11], with again a prominent role for the dominant pole (this time in relation to the saddle point). The key challenge is to bound the contributions of the contour integral when shifted beyond the dominant pole, a contribution which is substantial due to the relative large impact of the dominated poles. This saddle point method then provides a fully rigorous derivation of the asymptotic expression for P(Q > N ) and is of the form The function h(β) in this asymptotic expression involves infinite series and Riemann zeta functions that are reminiscent of the reflected Gaussian random walk [5,9,10]. Indeed, it follows from [16, Theorem 3] that our rescaled discrete queue converges under (1.3) to a reflected Gaussian random walk. Hence, the tail distribution of our system in the regime (1.3) should for large n be well approximated by the tail distribution of the reflected Gaussian random walk. We return to this connection in Subsection 5. Our approach thus relies on detailed knowledge about the distribution of all the poles of the PGF of Q, and in particular how this distribution scales with the asymptotic regime (1.2)-(1.3). As it turns out, in contrast with classical DPA, this many-sources regime makes that all poles contribute to the asymptotic characterization of the tail behavior. Our saddle point method leads to an asymptotic expansion for the tail probabilities, of which the limiting form corresponds to the heavy-traffic limit, and pre-limit forms present refined approximations for pre-limit systems (n < ∞) in heavy traffic. Such refinements to heavy-traffic limits are commonly referred to as corrected diffusion approximations [14,3,2]. Compared with the studies that directly analyzed the Gaussian random walk [5,9,10], which is the scaling limit of our queue in the many-sources regime, we start from the pre-limit process description, and establish an asymptotic result which is valuable for a queue with a finite yet large number of sources. Starting this asymptotic analysis from the actual pre-limit process description is mathematically more challenging than directly analyzing the process limit, but in return gives valuable insights into the manner and speed at which the system starts displaying its limiting behavior.
Outline of the paper. In Section 2 we describe the discrete queue in more detail and present some preliminary results for its stationary queue length distribution. In Section 3 we give an overview of the results and the contour integration representation for the tail distribution. In Section 4, we give a rigorous proof of the main result of the leading-order term using the dedicated saddle point method (Subsection 4.1), and of bounding the contour integral with integration paths shifted beyond the dominant pole (Subsection 4.2). In Section 5 we elaborate on the connection between the discrete queue and the Gaussian random walk, and we present an asymptotic series for P(Q > N ) comprising not only the dominant poles but also the dominated poles.

Model description and preliminaries
We consider a discrete stochastic model in which time is divided into periods of equal length. At the beginning of each period k = 1, 2, 3, ... new demand A k arrives to the system. The demands per period A 1 , A 2 , ... are assumed independent and equal in distribution to some non-negative integer-valued random variable A. The system has a service capacity s ∈ N per period, so that the recursion assuming Q 1 = 0, gives rise to a Markov chain (Q k ) k≥1 that describes the congestion in the system over time. The PGF is assumed analytic in a disk |z| < r with r > 1, which implies that all moments of A exist. We assume that A k is in distribution equal to the sum of work generated by n sources, X 1,k + ... + X n,k , where the X i,k are for all i and k i.i.d. copies of a random variable X, of which the PGF X(z) = ∞ j=0 P(X = j)z j has radius of convergence r > 1, and Under the assumption (2.3) the stationary distribution lim k→∞ P(Q k = j) = P(Q = j), j = 0, 1, . . . exists, with the random variable Q defined as having this stationary distribution. We let be the PGF of the stationary distribution.
It is a well-known consequence of Rouché's theorem that under (2.3) z s − A(z) has precisely s zeros in |z| ≤ 1, one of them being z 0 = 1. We proceed in this paper under the same assumptions as in [11]. We assume that |X(z)| < X(r 1 ), |z| = r 1 , z = r 1 , for any r 1 ∈ (0, r). Finally, we assume that the degree of X(z) is larger than s/n. Under these conditions, z 0 is the only zero of z s − A(z) on |z| = 1, and all others in |z| ≤ 1, denoted as z 1 , z 2 , ..., z s−1 , lie in |z| < 1. Furthermore, there are at most countably many zeros Z k of z s − A(z) in 1 < |z| < r, and there is precisely one, denoted by Z 0 , with minimum modulus.
There is the product form representation [4,17] where the right-hand side of (2.5) is analytic in |z| < Z 0 and has a first-order pole at z = Z 0 . We have for the tail probability (using that Q(1) = 1) for N = 0, 1, ...
. By contour integration, Cauchy's theorem and Q(1) = 1, we then get for 0 < ε < Z 0 − 1 where c 0 = Res z=Z 0 [Q(z)] and R is any number between Z 0 and min k =0 |Z k |. When n and s are fixed, we have that the integral on the second line of (2.7) is O(R −N ), and so there is the DPA In this paper we crucially rely on Pollaczek's integral representation for the PGF of Q that holds when |z| < 1 + ε < Z 0 (principal value of ln on |v| = 1 + ε).

Overview and results
In order to force the discrete queue to operate in the critical many-sources regimes, we shall assume throughout the paper the following relation between the number of sources n and the capacity s: with γ > 0 bounded away from 0 and ∞ as s → ∞. In this scaling regime, the zeros z j and Z k of z s − A(z) = 0 start clustering near z = 1, as described in the next lemma (proved in the appendix). Let z * j and Z * k denote the complex conjugates of z j and Z k , respectively.
Due to this clustering phenomenon, the main reasoning that underpin classical DPA cannot be carried over. Starting from the expression (2.8) we need to investigate what becomes of the term c 0 /(1 − Z 0 ), and moreover, the validity of the exponentially-small phrase in (2.8) and the actual N -range both become delicate matters that need detailed information about the distribution of the zeros as in Lemma 3.1.
Let us first present a result that identifies the relevant N -range: when N + 1 = L √ s with L > 0 bounded away from 0 and ∞.
Proof. We have from (3.2) and (3.5) that Z 0 = 1 + 2γµ when L is bounded away from 0 and ∞, and this gives the result.
From (2.5) we obtain the representation The next result will be proved in Section 4.
We thus get from (3.8) and Lemma 3.3 for Z ∈ C, |Z| ≥ 1 (principal logarithm). To handle the product P (z), in Lemma 3.4 below, we evaluate ln P (Z) for |Z| ≥ 1 in terms of the contour integral where ε > 0 is such that 1 < 1 + ε < Z 0 .
Lemma 3.4. Let ε > 0, 1 < 1 + ε < Z 0 and |Z| ≥ 1. Then The dedicated saddle point method, as considered in [11], applied to I(Z), with saddle point z sp = 1 + ε of the function g(z) = − ln z + n s ln(X(Z)), yields The next step consists of bounding the integral on the second line of (2.7), that can be written as by choosing R appropriately. To do this, we consider the product representation (2.5) of Q(z), and we want to choose R such that |z s − A(z)| ≥ C|z| s , |z| = R, for some C > 0 independent of s. It will be shown in Section 4 that this is achieved by taking R such that the curve |z s | = |A(z)|, on which Z 0 and Z ±1 lie, is crossed near a point z (also referred to as Z ±1/2 ), where z s and A(z) have opposite sign. A further analysis, using again the dedicated saddle point method to bound the product s−1 j−1 in (2.5), then yields that the integral in (3.16) decays as R −N . Finally using the asymptotic information in (3.2)-(3.4) for Z 0 and Z ±1 , with Z ±1/2 lying midway between Z 0 and Z ±1 , the integral on the second line of (2.7) can be shown to have relative order exp(−DN/ √ s), for some D > 0 independent of s, compared to the dominant-pole term in (2.8).
To summarize, we have now that of P(Q > N ) thus has a relative error that decays exponentially fast.
In Subsection 5.1 the stationary queue length Q, considered in the many-sources regime, is shown to be connected to the Gaussian random walk. This connection will imply that the front factor of the DPA in (3.17) satisfies where ln H(b 0 ) has a power series in b 0 with coefficients that can be expressed in terms of the Riemann zeta function. Combining this with Proposition 3.2 and (3.17) yields when N + 1 = L √ s with L bounded away from 0 and ∞. The leading term in (3.19) agrees with (1.6) when we identify and H(b 0 ) = h(β). In Subsection 5.2 we extend for a fixed M = 1, 2, . . . the approach in Section 4 by increasing the radius R of the integration contour in (3.16) to R M such that the poles Z 0 , Z ±1 , . . . , Z ±M are inside |z| = R M . this lead to (3.21) The front factors c k /(1 − Z k ) in the series in (3.21) satisfy with H k (b 0 ) some explicitly defined integral. When N + 1 = L √ s with L bounded away from 0 and ∞, we find from (3.22) and Proposition 3.2 that compare with (3.18), and it can be shown that this gives rise to an exp(−DL/ √ k) decay of the right-hand side of (3.23). The results in (3.19) and (3.23) together give precise information as to how the DPA arises, with leading behavior from the dominant pole, and lower order refinements coming from the dominated poles.

DPA through contour integration
In this section we present the details of getting approximations of the tail probabilities using a contour integration approach as outlined in Section 3. In Subsection 4.1, we concentrate on approximation of the front factor c 0 /(1 − Z 0 ) and the dominant pole Z 0 , and combine these to obtain an approximation of the leading-order term in (2.8). This gives Lemma 3.3, Lemma 3.4 and Proposition 3.5.
In Subsection 4.2 we assess and bound the integral on the second line of (2.7) and thereby make precise what exponentially small in (2.8) means in the present setting.

Approximation of the leading-order term
With the approximation (3.2), written as we get and Hence, by (4.3-4.5) and (A.6), where we have used d 0 of (4.3) in the last step. This gives (3.9).

Proof of Lemma 3.4
We have |A(z)| < |z s | when z = 1, 1 < |z| < Z 0 , and so ln(1 − A(z)/z s ) is analytic in z = 1, 1 < |z| < Z 0 . When |Z| > 1 + ε, we have by partial integration and Cauchy's theorem This gives the upper-case formula in (3.13), and the middle case follows in a similar manner by taking the residue at z = Z inside |z| = 1 + ε into account. For the lower case in (3.13), we use the result of the middle case, in which we take 1 < Z < 1 + ε, Z ↓ 1. We have I(Z) → I(1) as Z ↓ 1, and and this completes the proof.

Proof of Proposition 3.5
By (3.10) and Lemma 3.4 we have (4.10) and so ln Next, we consider the integral representation (3.12) of I(Z), where we take ε such that with z sp , see [11,Section 3], the unique point z ∈ (1, Z 0 ) such that d dz − ln z + n s ln(X(z)) = 0. (4.13) Observe that (4.14) and this suggests that I(Z 0 ) ≈ −I(1), a statement made precise below, since the main contribution to I(Z) comes from the z's in (3.12) close to z sp . We have, see [11], see [11,Subsection 5.3]. Hence, As to I(Z 0 ), we observe that, see (4.14), Thus, We next estimate the remaining integral in (4.20). With the substitution z = z(v), − 1 2 δ ≤ v ≤ 1 2 δ, we have A(z(v))/z s (v) = B exp(−sηv 2 ) with 0 < B < 1 and η > 0 bounded away from 1 and 0, respectively, and where c n are real. Then we get with exponentially small error Now we get from (4.14) and (4.21) that Inserting this into the integral on the second line of (4.22), we see that the −v in (4.25) cancels upon integration. Also

Bounding the remaining integral
We have from (2.7) where R ∈ (Z 0 , |Z ±1 |), and we intend to bound the integral at the right-hand side of (4.27). We use in (4.27) the Q(z) as represented by the right-hand side of (2.5) which is defined and analytic in z, |z| < r, z = Z k . We write for |z| < r, z = Z k with I(z) given by (3.12) and admitting an estimate since √ s|z − z sp |, B ∈ (0, 1) and η > 0 are all bounded away from 0. Therefore, there remains to be considered (z s − A(z)) −1 . We show below that there is a C > 0, independent of s, such that when z is on a contour K as in Figure , consisting of a straight line segment and a portion of the circle that are joined at the points (Re[Ẑ(± 1 2 )], ± 1 √ s y 0 ). Herê with a 0 , b 0 > 0 given in (A.10) and independent of s, approximates the solution z = Z(t), for real t small compared to s, of the equation outside the unit disk, according to Thus on K we have from (4.31) and we estimate . (4.38) Here we have used that |z − 1| ≥ Re[Ẑ(± 1 2 )] − 1 ≥ E/ √ s, z ∈ K, for some E > 0 independent of s. Observing thatẐ we see that for someD > 0 independent of s. Hence, by (4.36), we see that the relative error in (4.27) due to ignoring the integral at the right-hand side is of order exp(−DN/ √ s) with some D > 0, independent of s.
We show the inequality (4.31) for z ∈ K using the following property of X: there is a δ > 0 and a ϑ 1 ∈ (0, π/2) such that for any R ∈ [1, 1 + δ] the function |X(Re iϑ )| is decreasing in |ϑ| ∈ [0, This property follows from strict maximality of |X(e iϑ )| in ϑ ∈ [−π, π] at ϑ = 0 and analyticity of X(z) in the disk |z| < r (with r > 1). For the construction of the contour K in (4.31-4.32), we consider the quantity where v is of the form v = 2γµ with x 0 > 0 fixed and varying y ∈ R. We choose x 0 such that the outer curve Z(t) is crossed by z = 1 + v near the points Z(± 1 2 ), where z s − A(z) equals 2z s . Thus, we choose We have, as in the analysis in the appendix, that With x 0 > 0 fixed and independent of s, see (4.45), the leading part of the right-hand side in (4.46) is independent of s and describes, as a function of the real variable y, a parabola in the complex plane with real part bounded from above by its real value at y = 0 and that passes the imaginary axis at the points ±πi. Therefore, this leading part has a positive distance to all points 2πik, integer k. Now take y 0 such that 2γµ In Figure 1 we show the curve K (heavy), the approximationẐ(t) of the outer curve, and the choice y 0 = η 0 √ s for the case that γ = 1, µ/σ 2 = 2, s = 100. It follows from the above analysis, with v as in (4.44), that is bounded away from 0 and has a value 1 − c, 1 − c * at y = ±y 0 , where c is bounded away from 1 and |c| < 1. Now write When s is large enough, we have that R ∈ [1, 1 + δ] and 0 ≤ ϑ 0 ≤ ϑ 1 , where δ and ϑ 1 are as above in (4.42). We have 2 )] and η 0 = y 0 / √ s, and portion of the circle |z| = R with R = (ξ 2 + η 2 0 ) 1/2 . Choice of parameters: γ = 1, µ/σ 2 = 2, s = 100. and by (4.42) and monotonicity of |X(Re iϑ )|, ϑ 0 ≤ |ϑ| ≤ ϑ 1 , |X n (Re iϑ )| ≤ |X n (Re iϑ 0 )| ≤ |c|R s , ϑ 0 ≤ |ϑ| ≤ π.

Correction terms and asymptotic expansion
In this section we give a series expansion for the leading term in (3.15) involving the Riemann zeta function. We also show how to find an asymptotic series for P (Q > N ) as N → ∞ of which the term involving the dominant pole is the leading term. Before we do so, we first discuss how this leading term is related to the Gaussian random walk and a result of Chang and Peres [5].

Connection with Gaussian random walk
We know from [16,Theorem 3] that under the critical many-sources scaling, the rescaled queueing process converges to a reflected Gaussian random walk. The latter is defined as (S β (k)) k≥0 with S β (0) = 0 and copies of a normal random variable with mean −β and variance 1. Assume β > 0 (negative drift), and denote the all-time maximum of this random walk by M β . Denote by Q (s) the stationary congestion level for a fixed s (that arises from taking k → ∞ in (2.1)). Then, using ρ the spatially-scaled stationary queue length reaches the limit Q (s) /(σ √ n) d → M β as s, n → ∞ (see [12,15,16]).
The random variable M β was studied in [5,9,10]. In particular, [9, Thm. 1] yields, for β < 2 √ π, and from [5] we have P( where the second approximation holds for small values of β. We will now show how this second approximation in (5.5) follows from our leading term in the expansion.
Proof. It is shown in [11], Subsection 5.3 that in which we take the drift parameter β according to From [10] we have Then from Proposition 3.5, (5.7) and (5.9), we get the results in (5.6).

Asymptotic series for P (Q > N ) as N → ∞
When inspecting the argument that leads to (2.7), it is obvious that one can increase the radius R of the integration contour to values R M between |Z(±M )| and |Z(±(M + 1))| when M = 1, 2, ... is fixed. Here it must be assumed that s is so large that Z k increases in k = 0, 1, ..., M + 1. Then, the poles of Q(z) at z = Z ±k , k = 0, 1, ..., M , are inside |z| = R M , and we get As in Subsection 4.2, one can argue that the integral on the second line of (5.10) is relatively small compared to |Z M | −N −1 when R M is chosen between but away from |Z M | and |Z M +1 |.
We now need the following result.
Proof. This follows from the appendix with a similar argument as in the proof of Lemma 3.3.
As to the terms in the series in (5.10), we have for bounded k, see Lemma 5.2, Thus, we get the following result.
Theorem 5.6. There is the asymptotic series In the consideration of the terms with k = M − 1, M , it is tacitly accumed that s is so large that |Z k |, k = 0, 1, ..., M is a strictly increasing sequence.
A Proof of Lemma 3.1 We consider the zeros z j , j = 0, 1, ..., s − 1, and Z k , k ∈ Z, of the function z s − A(z) in the unit disk |z| ≤ 1 and in the annulus 1 < |z| < r, respectively, in particular those that are relatively close to 1. These zeros are elements of the set S A,s = {z ∈ C||z| < r, |z s | = |A(z)|}. For z ∈ S A,s , we have that ln(z s X −n (z)) is purely imaginary. We thus consider the equation