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 from 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 manysources 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 to the QED (quality and efficiency driven) regime for manyserver 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 non-negative 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 considering 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.i.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. Manysources 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, for example, [4,6,13,17,18] for specific applications).

How to adapt classical DPA?
As it turns out, the many-sources regime changes drastically the nature of the DPA. When the queue is pushed into the many-sources regime for 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 non-degenerate 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.
The approach in this paper to deal with these two issues uses the integral representation of the tail probabilities involving the PGF Q(z) of Q along a circle |z| = 1 + ε < Z 0 . By Cauchy's theorem, the tail probabilities can be written as a residue c 0 /(1 − Z 0 )Z N +1 0 at z = Z 0 and an integral along a contour shifted beyond Z 0 . Using the product expansion of Q(z), involving the zeros of a characteristic equation in |z| ≤ 1, the quantity c 0 /(1 − Z 0 ) can be approximated in terms of an integral of the Pollaczek type that has been considered in [11]. In [11], a dedicated saddle point method has been developed to approximate Pollaczek-type integrals for the mass at zero, the mean, and the variance of Q in the many-sources regime. This dedicated saddle point method can also be used to find a many-sources approximation for c 0 /(1 − Z 0 ). The remaining challenge is then to bound the contribution of the contour integral shifted beyond the dominant pole Z 0 . This depends on how far the integration path can be shifted, and this is determined by the positions of the first dominated poles. It turns out that, in the many-sources regime, the dominated poles have approximations of the type (1.4) as well. The integration path is then chosen, roughly, as the circle that passes halfway between Z 0 and the first dominated pole. This, together with the dedicated saddle point method to approximate c 0 /(1 − Z 0 ), then provides a fully rigorous derivation of the asymptotic expression for P(Q > N ) that is of the form (1.6) 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 Sect. 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 means 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 [2,3,14]. 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 Sect. 2 we describe the discrete queue in more detail and present some preliminary results for its stationary queue length distribution. In Sect. 3 we give an overview of the results and the contour integration representation for the tail distribution. In Sect. 4 we give a rigorous proof of the main result for the leading-order term using the dedicated saddle point method (Subsect. 4.1), and we bound the contour integral with integration path shifted beyond the dominant pole (Subsect. 4.2). In Sect. 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]. Thus, we assume that |X (z)| < X (r 1 ), |z| = r 1 , z = r 1 , for any r 1 ∈ (0, r ); see [11], end of Sect. 2 for a discussion of this condition. 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 by 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. The zeros Z k are indexed using integer k; they come in conjugate pairs and we let Z k = Z * −k where Z k is in the (closed) lower half plane for non-negative k. For our analysis, it is crucial that r > 1, and so certain heavy-tailed distributions for X (with r = 1) are excluded.
There is the product form representation [4,17] that gives the analytic extension of Q(z) to all z, |z| < r and z not a zero of z s − A(z), 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 consider the heavy-traffic setting described by (1.2) and (1.3) with n → ∞. In this setting R and Z 0 in (2.7) tend to 1, and thus both terms on the second line of (2.7) need special attention. The second term here is approximated, using the product form expansion in (2.5), in terms of a contour integral of the Pollaczek type to which the dedicated saddle point method of [11], Sect. 3, can be applied. The method has been developed in [11] from Pollaczek's integral representation of the PGF of Q, which holds when |w| < 1 + ε < Z 0 , to get limiting expressions for moments of Q in the generalized heavy-traffic regimes n/s = 1 − γ n −α , α > 0, and n → ∞.
To this end, the contour integrals are deformed so as to pass through the saddle point where B = exp(sg(z sp )) and η = g (z sp ). It is shown in [11], p. 796, using Lagrange's inversion theorem, that in the heavy-traffic regime n/s = 1 − γ / √ n, for a δ > 0 independent of s one has with real coefficients c j , while (2.10) holds with B positive and bounded away from 0 and 1, and η positive and bounded away from 0.

Overview and results
In order to force the discrete queue to operate in the critical many-sources regime, we shall assume throughout the paper the following relation between the number of sources n and the capacity s = s n (γ ): 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 underpins 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 Sect. 4.
We thus get from (3.8) and Lemma 3.3 where 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 (3.14) Combining (3.2), (3.5), (3.8), (3.9), and (3.14) then gives one of our main results: 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 Sect. 4 that this is achieved by taking R such that the curve |z s | = |A(z)|, on which Z 0 and Z ±1 (= Z +1 , 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 In Subsect. 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 Propositions 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 Subsect. 5.2, we extend for a fixed M = 1, 2, . . . the approach in Sect. 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 leads to 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), (3.21) 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 Sect. 3. In Subsect. 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 Lemmas 3.3 and 3.4, and Proposition 3.5.
In Subsect. 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.

Proof of Lemma 3.3
From With the approximation (3.2), written as we get and Hence, by (4.3-4.5) and (5.32), 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, using that d[ln(1 − z/Z )] = −dz/(Z − z) 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 (4.8) and this completes the proof.
We have, see [11], see [11,Subsect. 5.3]. Hence, As for 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), with 0 < B < 1 and η > 0 bounded away from 1 and 0, respectively, and where c j 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 on 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 (4.28) Now s −μ A = γ √ s, and by Lemma 3.4 and (4.18), we have s−1 with I (z) given by (3.12) and admitting an estimate ∈ (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 |z s − A(z)| ≥ C|z| s (4.31) when z is on a contour K as in Fig. 1, 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ê outside the unit disk, according to Thus on K , we have from (4.31) and we estimate 1 2πi . (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 on 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, 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 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 In Fig. 1, we show the curve K (bold), 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 and by (4.42) and monotonicity of |X (Re iϑ )|, ϑ 0 ≤ |ϑ| ≤ ϑ 1 , Therefore, (4.31) holds on K with positive and bounded away from 0 as s gets large.

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 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, Theorem 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,Subsect. 5 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 (5.10) As in Subsect. 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 to 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.

Proposition 5.3 For bounded k ∈ Z,
(5.14) We aim at approximating I (Z k ), showing, in particular, that c k /(1 − Z k ) = 0 is bounded away from 0 for bounded k and large s. To that end, we conduct the dedicated saddle point analysis for I (Z k ). We have for |Z | ≥ Z 0 , Re(Z ) > z sp , and z(v) as in (4.21) and defined implicitly by g(z(v)) = g(z sp ) − 1 2 v 2 g (z sp ). We then find, by using z(v) = z sp + iv + O(v 2 ) and z (v) = i + O(v), that where in the last step the substitution t = v √ sη/2 has been made. Combining in the last integrand in (5.17) the values at t and −t for t ≥ 0, we get the following result.