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

Dominant poles and tail asymptotics in the critical Gaussian many-sources regime
A. J. E. M. Janssen 0
J. S. H. van Leeuwaarden 0
0 Department of Mathematics and Computer Science, Eindhoven University of Technology , P.O. Box 513, 5600 MB Eindhoven , The Netherlands
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.
Heavy traffic; Many sources; Asymptotics; Dominant pole approximation; Saddle point method; QED regime
1 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
B J. S. H. van Leeuwaarden
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 Z0, the singularity in |z| > 1 closest
to the origin. For PGFs, due to Pringsheim’s theorem [7], Z0 is always a positive real
number larger than one. Then DPA leads to the approximation
for large N ,
with c0 = limz→Z0 (z − Z0)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.
1.1 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 σA2 . 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 sn, so that the system
utilization is given by ρn = nμ/sn, where the index n expresses the dependence on
the scale at which the system operates. The traditional capacity sizing rule would then
be
√n
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 = sn ∼
nμ)
βσ γ
∼ 1 − μ√n = 1 − √sn + O(n−1), as 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 sn 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 sn, 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).
1.2 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
Z0 ∼ 1 + √
, as n → ∞.
∼ 1 + √
→ e−2β K ∈ (0, 1) as n → ∞.
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. Now let us turn to the other two remaining issues:
0
the fact that c0/(1 − Z0) 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 + ε < Z0. By Cauchy’s theorem, the tail probabilities can be written as
a residue c0/(1 − Z0)Z0N +1 at z = Z0 and an integral along a contour shifted beyond
Z0. Using the product expansion of Q(z), involving the zeros of a characteristic
equation in |z| ≤ 1, the quantity c0/(1 − Z0) 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 c0/(1 − Z0). The remaining challenge is then to bound the contribution of the
contour integral shifted beyond the dominant pole Z0. 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 Z0 and the first dominated pole.
This, together with the dedicated saddle point method to approximate c0/(1 − Z0),
then provides a fully rigorous derivation of the asymptotic expression for P(Q > N )
that is of the form
nσ 2) ∼ h(β) · e−2β K , as n → ∞.
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.
1.3 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.
2 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 Ak arrives to the
system. The demands per period A1, A2, ... 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
Qk+1 = max{Qk + Ak − s, 0},
k = 1, 2, ...,
assuming Q1 = 0, gives rise to a Markov chain (Qk )k≥1 that describes the congestion
in the system over time. The PGF
A(z) =
j=0
P( A = j )z j
0 < μA = nμ = n X (1) < s.
Under the assumption (2.3), the stationary distribution limk→∞ P(Qk = 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) zs − A(z)
has precisely s zeros in |z| ≤ 1, one of them being z0 = 1. We proceed in this
paper under the same assumptions as in [11]. Thus, we assume that |X (z)| < X (r1),
|z| = r1, z = r1, for any r1 ∈ (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, z0 is the only zero of zs − A(z) on |z| = 1, and all others in |z| ≤ 1,
denoted by z1, z2, ..., zs−1, lie in |z| < 1. Furthermore, there are at most countably
many zeros Zk of zs − A(z) in 1 < |z| < r , and there is precisely one, denoted by
Z0, with minimum modulus. The zeros Zk are indexed using integer k; they come in
conjugate pairs and we let Zk = Z −∗k where Zk 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]
Q(z) =
(s − μA)(z − 1) s−1 z − z j
zs − A(z)
j=1
1 − z j
that gives the analytic extension of Q(z) to all z, |z| < r and z not a zero of zs − A(z),
where the right-hand side of (2.5) is analytic in |z| < Z0 and has a first-order pole at
z = Z0. We have, for the tail probability (using that Q(1) = 1), for N = 0, 1, ...
i=N +1
P(Q > N ) =
P(Q = i ) = CzN
where CzN [ f (z)] denotes the coefficient of z N of the function f (z). By contour
integration, Cauchy’s theorem and Q(1) = 1, we then get for 0 < ε < Z0 − 1
|z|=R
where c0 = Resz=Z0 [Q(z)] and R is any number between Z0 and mink=0 |Zk |. 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
P(Q > N ) =
(1 + exponentially small),
N → ∞.
In this paper we consider the heavy-traffic setting described by (1.2) and (1.3) with
n → ∞. In this setting R and Z0 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,
Q(w) = exp
exp(sg(z(v))) = B exp(− 21 sηv2),
where B = exp(sg(zsp)) and η = g (zsp). 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.
3 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 = sn(γ ):
z(v) = zsp + i v +
j=2
with γ > 0 bounded away from 0 and ∞ as s → ∞. In this scaling regime, the zeros
z j and Zk of zs − A(z) = 0 start clustering near z = 1, as described in the next lemma
(proved in the appendix). Let z∗j and Zk∗ denote the complex conjugates of z j and Zk ,
respectively.
z0 = 1,
Z0 = 1 +
+ O(s−1),
a0 =
b0 =
and principal roots in (3.3)–(3.5).
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 c0/(1 − Z0), 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:
Proposition 3.2
2γ μ
Proof We have from (3.2) and (3.5) that Z0 = 1 + σ 2√s + O 1s . Hence
1 2γ μ
Z N +1 = exp −L√s ln 1 + σ 2√s + O(s−1)
0
From (2.5), we obtain the representation
s−1 Z0 − z j
The next result will be proved in Sect. 4.
1 + O(s−1/2) ,
s−1
j=1
s−1
j=1
P(Z ) =
(1 − z j /Z ) = exp
ln(1 − z j /Z )
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
ln(1 − z−s A(z))
Z − z
We thus get from (3.8) and Lemma 3.3
The dedicated saddle point method, as considered in [11], applied to I (Z ), with
saddle point zsp = 1 + ε of the function g(z) = − ln z + ns ln(X (Z )), yields
Combining (3.2), (3.5), (3.8), (3.9), and (3.14) then gives one of our main results:
Proposition 3.5
− ln(1 − Z −1) + I (Z ), 1 + ε < |Z | < r,
ln(γ √s) + I (1), Z = 1.
1 − Z0
= − ln(4b02) + 2 ln[Q(0)] + O(s−1/2).
|z|=R
z N +1(1 − z)
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 |zs − 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 |zs | = | A(z)|, on which Z0 and Z±1(= Z+1, Z−1) lie, is crossed
near a point z (also referred to as Z±1/2) where zs and A(z) have opposite sign. A
further analysis, using again the dedicated saddle point method to bound the product
sj−=11 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 Z0 and Z±1, with Z±1/2 lying midway
between Z0 and Z±1, the integral on the second line of (2.7) can be shown to have
relative order exp(−D N /√s), for some D > 0 independent of s, compared to the
dominant pole term in (2.8).
To summarize, we have now that
1 + O(e−DN/√s ) ,
1 − Z0
= H (b0) 1 + O(s−1/2) ,
where ln H (b0) has a power series in b0 with coefficients that can be expressed in
terms of the Riemann zeta function. Combining this with Propositions 3.2 and 3.17
yields
P(Q > N ) = H (b0) exp
1 + O(s−1/2)
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
1 − Zk
P(Q > N ) = Re
+ O |Z M+1|−N .
1 − Zk
with Hk (b0) 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
= H (bk ) exp
2
− La0( b0 − 2π i k + b0)
1 + O(s−1/2) ;
compare with (3.18), and it can be shown that this gives rise to an exp(−D L/√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.
4 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 c0/(1 − Z0) and the dominant pole
Z0, 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.
4.1 Approximation of the leading-order term
4.1.1 Proof of Lemma 3.3 From
Z 0s = A(Z0) = X n (Z0), μA = nμ = s 1 − √γ
s
, A (Z0) = n X (Z0)X n−1(Z0),
With the approximation (3.2), written as
s Z 0s−1 − A (Z0)
X (Z0)Z0
X (1)X (Z0)
d0 =
X (Z0) = X (1) + X (1)(Z0 − 1) + O(s−1) = μ +
+ O(s−1) (4.4)
Hence, by (4.3–4.5) and (5.32),
1 −
γ X (Z0)Z0
1 − √s X (1)X (Z0)
μd0
X (Z0) = X (1) + X (1)(Z0 − 1) + O(s−1) = 1 + √s + O(s−1).
= 1 −
= 1 −
= 1 −
+ O(s−1)
+ O(s−1)
where we have used d0 of (4.3) in the last step. This gives (3.9).
4.1.2 Proof of Lemma 3.4
We have | A(z)| < |zs | when z = 1, 1 < |z| < Z0, and so ln(1 − A(z)/zs ) is analytic
in z = 1, 1 < |z| < Z0. When |Z | > 1 + ε, we have by partial integration, using that
d[ln(1 − z/Z )] = −d z/(Z − z) and Cauchy’s theorem,
z (1 − A(z)/zs )
ln 1 − Z 1 − A(z)/zs d z
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
+ ln 1 −
= ln[(Z s − X n(Z )) |Z=1] = ln(s − nμ) = ln(γ √s),
and this completes the proof.
4.1.3 Proof of Proposition 3.5 By (3.10) and Lemma 3.4, we have ln
1 − Z0
From (3.2), it follows that
= ln P(Z0) − ln P(1) + O(s−1/2)
− ln(γ √s) + I (Z0) − I (1) + O(s−1/2). (4.9)
+ O(s−1/2) = − ln(4b02) + O(s−1/2),
= I (Z0) − I (1) − ln(4b02) + O(s−1/2).
We have g(1) = 0 = g(Z0), and so z = zsp for which g (z) = 0 lies about midway
between 1 and Z0. More specifically, by expanding g (z) = g (1) + (z − 1)g (1) +
O(s−1), 1 ≤ z ≤ Z0, and expressing g (1), g (1) in terms of γ , μ, σ, s, while using
(4.3), it follows that
zsp = 21 (1 + Z0) + O(s−1),
Z0 − zsp = zsp − 1 + O(s−1).
ln 1 − z−s A(z)
z(z − 1)
|z|=zsp
ln 1 − z−s A(z) |d z| = O(s−1/2);
|z|=zsp
see [11, Subsect. 5.3]. Hence,
−1
|z|=zsp
ln 1 − z−s A(z)
z − 1
As for I (Z0), we observe that, see (4.14),
d z = ln[ P(Q = 0)] + O(s−1/2).
(z(v) − zsp)z (v) = −v + O(v2).
I (Z0) = −I (1) + 2
|z|=zsp
z − zsp
(Z0 − z)(z − 1)
ln 1 − z−s A(z) d z + O(s−1/2). (4.20)
j=2
We next estimate the remaining integral in (4.20). With the substitution z = z(v),
− 2 δ ≤ v ≤ 21 δ, as explained at the end of Sect. 3, we have A(z(v))/zs (v) =
1
B exp(− 21 sηv2) with 0 < B < 1 and η > 0 bounded away from 1 and 0, respectively,
and
z(v) = zsp + i v +
where c j are real. Then we get with exponentially small error
|z|=zsp
z − zsp
(Z0 − z)(z − 1)
ln 1 − z−s A(z) d z
(z(v) − zsp)z (v)
(Z0 − z(v))(z(v) − 1)
ln(1 − Be− 21 sηv2 )dv.
Now we get from (4.14) and (4.21) that
(z(v) − zsp)z (v)
(Z0 − z(v))(z(v) − 1) =
−v + O(v2)
|zsp − 1|2 + v2 + O
s−1 + v2
Inserting this into the integral on the second line of (4.22), we see that the −v in (4.25)
cancels upon integration. Also
v2
|zsp − 1|2 + v2 ln(1 − Be− 21 sηv2 )dv = O(s−1/2),
and this finally shows that the integral in (4.22) is O(s−1/2). Then combining (4.11),
(4.18), and (4.20), we get the result.
4.2 Bounding the remaining integral We have from (2.7)
P(Q > N ) =
|z|=R
where R ∈ (Z0, |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 = Zk . We write for |z| < r , z = Zk
Q(z) −1
(1 − z)z N +1 = z N +2−s
sj−=11(1 − z j ) zs − A(z) j=1
s−1
Now s − μA = γ √s, and by Lemma 3.4 and (4.18), we have
C γ √s for some C > 0 independent of s. Hence (s − μA)/
in s. Next, for |z| ≥ Z0, we have by Lemma 3.4
(1 − z j /z).
sj−=11(1 − z j ) = P(1) ≥
sj−=11(1 − z j ) is bounded
with I (z) given by (3.12) and admitting an estimate
s−1
j=1
|I (z)| = O |z − zsp|−1
= O(1)
since √s|z − zsp|, B ∈ (0, 1) and η > 0 are all bounded away from 0. Therefore,
there remains to be considered (zs − A(z))−1. We show below that there is a C > 0,
independent of s, such that
|zs − A(z)| ≥ C |z|s
when z is on a contour K as in Fig. 1, consisting of a straight line segment
and a portion of the circle
|z| = R =
Zˆ (t ) = 1 + √a0s ((b02 − 2π i t )1/2 + b0),
with a0, b0 > 0 given in (5.36) 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)
ln X (z) − ln z =
Z (t ) = Zˆ (t ) + O
and we estimate
= O
Re[Zˆ (± 21 )]
Zˆ (0) = 1 +
we see that
Re[Zˆ (± 21 )]
1 −
( 21 b02 + 21 (b04 + π 2)1/2)1/2 − b0
1 + √a0s [( 21 b02 + 21 (b04 + π 2)1/2)1/2 + b0]
= O(exp(−Dˆ N /√s))
for some Dˆ > 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(−D N /√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 (Reiϑ )|
is decreasing in |ϑ | ∈ [0, ϑ1] while
ϑ1 ≤ |ϑ | ≤ π ⇒ |X (Reiϑ )| ≤ |X (Reiϑ1 )|.
This property follows from strict maximality of |X (eiϑ )| 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
n ln[X (1 + v)] − ln(1 + v),
= Zˆ (0) − 1 +
with x0 > 0 fixed and varying y ∈ R. We choose x0 such that the outer curve Z (t ) is
crossed by z = 1 + v near the points Z (± 21 ), where zs − A(z) equals 2zs . Thus, we
choose
We have, as in the analysis in the appendix, that
n ln[X (1 + v)] − s ln(1 + v) =
γ μ
+ 2i σ 2 − x0 y + O(sv3) + O(v2√s).
With x0 > 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π i k, integer k. Now take y0 such that
In Fig. 1, we show the curve K (bold), the approximation Zˆ (t ) of the outer curve, and
the choice y0 = η0√s for the case that γ = 1, μ/σ 2 = 2, s = 100.
It follows from the above analysis, with v as in (4.44), that
1 −
− y0 ≤ y ≤ y0,
Reiϑ0 = Zˆ (0) +
= 1 + v0.
is bounded away from 0 and has a value 1 − c, 1 − c∗ at y = ±y0, where c is bounded
away from 1 and |c| < 1. Now write
C = min 1 − |c|, |ym|≤iny0 1 −
positive and bounded away from 0 as s gets large.
5 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].
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
|X n(1 + v0)| ≤ |c| |1 + v0|s ,
|X n(Reiϑ )| ≤ |X n(Reiϑ0 )| ≤ |c|Rs ,
Therefore, (4.31) holds on K with
|z| = R
Fig. 1 Integration curve K (bold) consisting of line segment z = ξ + i η, −η0 ≤ η ≤ η0, where ξ =
Re[Zˆ (± 21 )] and η0 = y0/√s, and portion of the circle |z| = R with R = (ξ 2 + η02)1/2. Choice of
parameters: γ = 1, μ/σ 2 = 2, s = 100
5.1 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
with Y1, Y2, . . . i.i.d. 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 ρs = 1 − γ /√s, with
r=0
exponentially fast as K → ∞.
Hence, there are the approximations
nσ 2) ≈ P(Mβ > K ) ≈ h(β) · e−2β K , as n → ∞,
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.
+ O(s−1/2),
Proposition 5.1
1 − Z0
r=0
Proof It is shown in [11, Subsect. 5.3] that
ln[Q(0)] = ln[ P(Mβ = 0)] + O(s−1/2),
From [10], we have
b0
ln[ P(Mβ = 0)] = ln(2b0) + √
r=0
Then from Proposition 3.5, (5.7) and (5.9), we get the results in (5.6).
5.2 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 RM between |Z (±M )| and |Z (±(M +1))|
when M = 1, 2, ... is fixed. Here it must be assumed that s is so large that |Zk | increases
in k = 0, 1, ..., M + 1. Then, the poles of Q(z) at z = Z±k , k = 0, 1, ..., M , are inside
|z| = RM , and we get
(1 − Zk )ZkN +1
|z|=RM
when k = o(s).
− s Zks−1 − A (Zk )
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,
s−1 Zk − z j
s−1 1 − z j /Zk
1 − z j
· 1 + O(s−1/2) .
Furthermore, according to Lemma 3.4,
s−1 1 − z j /Zk
j=1
1 − z j
= I (Zk ) − I (1) − ln 2b0 b0 +
+ O(s−1/2).
Thus, we get the following result.
Proposition 5.3 For bounded k ∈ Z,
1 − Zk =
exp(I (Zk ) − I (1))
1 + O(s−1/2) .
We aim at approximating I (Zk ), showing, in particular, that ck /(1 − Zk ) = 0 is
bounded away from 0 for bounded k and large s. To that end, we conduct the dedicated
saddle point analysis for I (Zk ). We have for |Z | ≥ Z0, Re(Z ) > zsp,
ln 1 − z−s A(z)
Z − z
with exponentially small error in the last identity as s → ∞. With g(z) = − ln z +
ns ln X (z), we let
and z(v) as in (4.21) and defined implicitly by g(z(v)) = g(zsp) − 21 v2g (zsp). We
then find, by using z(v) = zsp + i v + O(v2) and z (v) = i + O(v), that
−1
v + i (Z − zsp)
dv + O(s−1/2)
−1
ln(1 − Be−t2 )
dt + O(s−1/2),
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.
Proposition 5.4 For |Z | ≥ Z0, Re(Z ) > zsp,
I (Z ) = J (d) + O(s−1/2),
where d = (Z − zsp)√sη/2, and
In the context of Proposition 5.3, we consider
d = dk = (Zk − zsp) sη/2,
We have that
dk = dˆk + O(s−1/2) ;
|dˆk | ≥ b0,
k ∈ Z.
Since for t ≥ 0 and arg(d) ∈ (− 41 π, 41 π ) we have
ln(1 − Be−t2 ) < 0,
we see that we have complete control on the quantities J (dˆk ) (also note (5.16) for this
purpose). Using that −I (1) = I (Z0) + O(s−1/2), see (4.20), we get the following
result.
Proposition 5.5 For bounded k ∈ Z,
1 − Zk =
exp( J (dˆk ) + J (dˆ0))
2dˆk (dˆk + dˆ0)
+ O(s−1/2),
where the leading quantity in (5.25) = 0, ∞, dˆk is given in (5.22) with dˆ0 = b0, and
J is given in (5.19).
Theorem 5.6 There is the asymptotic series
P(Q > N ) ∼ Re
Proof This follows from (5.10), in which the integral is o(|Z M |−N ) and the term
with k = M is O (|Z M |−N ), while the reciprocal of the term with k = M − 1
is O (|Z M−1|−N ) by Proposition 5.5. In the consideration of the terms with k =
M − 1, M , it is tacitly assumed that s is so large that |Zk |, k = 0, 1, ..., M is a strictly
increasing sequence.
Acknowledgments This work was financially supported by The Netherlands Organization for Scientific
Research (NWO) and by an ERC Starting Grant.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons license, and indicate if changes were made.
Appendix: Proof of Lemma 3.1
We consider the zeros z j , j = 0, 1, ..., s − 1, and Zk , k ∈ Z, of the function zs − 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 SA,s = {z ∈ C||z| <
r, |zs | = | A(z)|}. For z ∈ SA,s , we have that ln(zs X −n (z)) is purely imaginary. We
thus consider the equation
with z near 1 and t small compared to s. Writing
z = 1 + v,
√
Dividing by s and using that n X (1) = 1 − γ / s yields
s
v − 21 v2 + O (v3) =
2 X (1)
where we have used that
= O
σ 2 = X (1) − ( X (1))2 + X (1) > 0.
Taking square roots on both sides of (5.33) and using that the leading term on the
right-hand side of (5.33) has order (1 + |u|)/s, we get
v =
1 −
Irrespective of the ±-sign, the leading part of the right-hand side of (5.34) has order
((1 + |u|)/s)1/2 (and even (|u|/s)1/2 in the case of the −-sign), and the O -term has
order (1 + |u|)/s which is o(((1 + |u|)/s)1/2) as long as |u|/s = o(1). Thus, in that
regime of u, we have
z = 1 + v = 1 + σ 2√s
1 −
1 + |u|
,
s
where we have inserted
a0 =
b0 =
In the case of the minus sign in (5.35), the O -term may be replaced by O (|u|/s).
Choosing u = 2π j with j = 0, 1, ..., and j = o(s), we get from (5.35) with the
minus sign, (3.3). Choosing u = 2π k with k ∈ Z and k = o(s), we get from (5.35)
with the plus sign, (3.4).
v −
1 −
1. Anick , D. , Mitra , D. , Sondhi , M.M. : Stochastic theory of a data-handling system with multiple sources . Bell Syst. Tech. J . 61 ( 8 ), 1871 - 1894 ( 1982 )
2. Asmussen , S. : Applied Probability and Queues , 2nd edn. Springer, New York ( 2003 )
3. Blanchet , J. , Glynn , P. : Complete corrected diffusion approximations for the maximum of a random walk . Ann. Appl. Probab . 16 ( 2 ), 951 - 983 ( 2006 )
4. Bruneel , H. , Kim , B.G. : Discrete-Time Models for Communication Systems Including ATM . Kluwer Academic Publishers, Boston ( 1993 )
5. Chang , J.T. , Peres , Y. : Ladder heights, Gaussian random walks and the Riemann zeta function . Ann. Probab . 25 ( 2 ), 787 - 802 ( 1997 )
6. Dai , J.G. , Shi , P. : A two-time-scale approach to time-varying queues for hospital inpatient flow management . Soc. Sci. Res. Network ( 2015 ). doi:10.2139/ssrn.2489533
7. Flajolet , P. , Sedgewick , R.: Analytic Combinatorics, 1st edn. Cambridge University Press, New York ( 2009 )
8. Halfin , S. , Whitt , W.: Heavy-traffic limits for queues with many exponential servers . Oper. Res . 29 , 567 - 588 ( 1981 )
9. Janssen , A.J.E.M. , van Leeuwaarden , J.S.H.: On Lerch's transcendent and the Gaussian random walk. Ann. Appl. Probab . 17 , 421 - 439 ( 2006 )
10. Janssen , A.J.E.M. , van Leeuwaarden , J.S.H.: Cumulants of the maximum of the Gaussian random walk . Stoch. Process. Appl. 117(12) , 1928 - 1959 ( 2007 )
11. Janssen , A.J.E.M. , van Leeuwaarden , J.S.H. , Mathijsen , B.W.J. : Novel heavy-traffic regimes for largescale service systems . SIAM J. Appl. Math. 75 , 212 - 787 ( 2015 )
12. Jelenkovic´ , P. , Mandelbaum , A. , Momcˇilovic , P. : Heavy traffic limits for queues with many deterministic servers . Queueing Syst . 47 , 53 - 69 ( 2004 )
13. Newell , G.F. : Queues for a fixed-cycle traffic light . Ann. Math. Stat . 31 , 589 - 597 ( 1960 )
14. Siegmund , D. : Sequential Analysis . Springer Series in Statistics. Springer, New York ( 1985 )
15. Sigman , K. , Whitt , W.: Heavy-traffic limits for nearly deterministic queues . J. Appl. Probab . 48 ( 3 ), 657 - 678 ( 2011 )
16. Sigman , K. , Whitt , W.: Heavy-traffic limits for nearly deterministic queues: stationary distributions . Queueing Syst. Theory Appl . 69 ( 2 ), 145 - 173 ( 2011 )
17. van Leeuwaarden , J.S.H.: Queueing models for cable access networks . Ph.D. thesis , Eindhoven University of Technology ( 2005 )
18. van Leeuwaarden , J.S.H.: Delay analysis for the fixed-cycle traffic light queue . Transp. Sci . 40 , 189 - 199 ( 2006 )
19. Van Mieghem , P. : The asymptotic behavior of queueing systems: large deviations theory and dominant pole approximation . Queueing Syst. Theory Appl . 23 ( 1-4 ), 27 - 55 ( 1996 )