brainstorming

# Exponential Tail Bounds on Queues: A Confluence of Non-Asymptotic Heavy
Traffic and Large Deviations

Prakirt Raj Jhunjhunwala<EMAIL_ADDRESS>Daniela Hurtado-Lange
<EMAIL_ADDRESS>Siva Theja Maguluri<EMAIL_ADDRESS>
###### Abstract

In general, obtaining the exact steady-state distribution of queue lengths is
not feasible. Therefore, our focus is on establish bounds for the tail
probabilities of queue lengths. Specifically, we examine queueing systems
under Heavy-Traffic (HT) conditions and provide exponentially decaying bounds
for the probability $\mathbb{P}(\epsilon q>x)$, where $\epsilon$ is the HT
parameter denoting how far the load is from the maximum allowed load. Our
bounds are not limited to asymptotic cases and are applicable even for finite
values of $\epsilon$, and they get sharper as $\epsilon\to 0$. Consequently,
we derive non-asymptotic convergence rates for the tail probabilities. Unlike
other approaches such as moment bounds based on drift arguments and bounds on
Wasserstein distance using Stein’s method, our method yields sharper tail
bounds. Furthermore, our results offer bounds on the exponential rate of decay
of the tail, given by $-\frac{1}{x}\log\mathbb{P}(\epsilon q>x)$ for any
finite value of $x$. These can be interpreted as non-asymptotic versions of
Large Deviation (LD) results. To obtain our results, we use an exponential
Lyapunov function to bound the moment generating function of queue lengths and
apply Markov’s inequality.

We demonstrate our approach by presenting tail bounds for: (i) a continuous
time Join-the-shortest queue (JSQ) load balancing system, (ii) a discrete time
single-server queue and (iii) an $M/M/n$ queue. We not only bridge the gap
between classical-HT and LD regimes but also explore the large system HT
regimes for JSQ and $M/M/n$ systems. In these regimes, both the system size
and the system load increase simultaneously. Our results also close a gap in
the existing literature on the limiting distribution of JSQ in the super-NDS
(a.k.a. super slowdown) regime. This contribution is of an independent
interest. Here, a key ingredient is a more refined characterization of state
space collapse for JSQ system, achieved by using an exponential Lyapunov
function designed to approximate the $\ell_{\infty}$ norm.

###### keywords:

Classical heavy traffic, Large deviations, Tail probabilities, Join-the-
shortest queue, Transform Method, exponential Lyapunov function

[1]organization=Georgia Institute of Technology, addressline=North Avenue,
city=Atlanta, postcode=30332 GA, country=USA

[2]organization=William and Mary, city=Williamsburg, postcode=VA 23185,
country=USA

[3]organization=Georgia Institute of Technology, addressline=North Avenue,
postcode=30332 GA, city=Atlanta, country=USA

##  1 Introduction

Queueing models are used to study performance of many systems such as cloud
computing, data centers, ride hailing, call centers etc. In general, obtaining
the complete distribution of queue lengths in these systems is intractable.
Therefore, a common approach is to study asymptotic regimes. There are several
popular regimes such as Heavy Traffic (HT), large scale regime, or Large
Deviations (LD). In the HT regime, the system is loaded close to its maximum
capacity while keeping the number of servers fixed. In the large-systems
regime, the system’s load is fixed, but the number of servers is increased to
infinity. And in the LD limit one studies the probability of rare events, that
is, the tail probability for large thresholds. Recently, the Many-Server
Heavy-Traffic (Many-Server-HT) regime has gained more popularity, where the
system is loaded to maximum capacity while simultaneously increasing the
number of servers. The system’s behavior varies greatly depending on how
quickly the load increases relative to the number of servers. As such one
employs very different analysis techniques to study queueing systems in
different regimes.

In the study of HT asymptotics, one typically scales the queue lengths using a
parameter that represents the system’s load. By denoting the load as
$1-\epsilon$, the HT limit is achieved when $\epsilon$ approaches zero. Most
of the literature focuses on systems that satisfy the so-called Complete
Resource Pooling (CRP) condition and behave like a single-server queue in the
limit. For such systems, it is well-known that the scaled queue length follows
an exponential distribution in the HT limit, which gives the tail
probabilities of the limiting system. However, the rate of convergence of the
tail probabilities (of the pre-limit system) to the corresponding HT value
remains unknown.

Most real world systems involve Service Level Agreements (SLA), where
customers are promised a specific level of service, including the maximum
delay they can expect. Motivated by this, in this paper, we focus on
establishing sharp bounds on the tail probabilities of scaled queue length of
the pre-limit system, i.e., for $\epsilon>0$. In particular, we get non-
asymptotic bounds of the form

$\mathbb{P}(\epsilon q>x)\leq\kappa(\epsilon,x)e^{-\theta(\epsilon)x},$

where $q$ represents the total queue length in steady state. Here,
$\theta(\epsilon)$ gives the decay rate of the tail probability of the pre-
limit system, and $\theta(\epsilon)$ converges to the correct HT value as
$\epsilon\rightarrow 0$. Recent results show the rate of convergence to HT in
terms of the mean, moments, or Wasserstein’s distance (for references on each
of these, see Section 1.3). These methods focus on the entire distribution of
the queue lengths and drown the tail. For example, consider the second moment,
and suppose $\epsilon q$ converges in distribution to the random variable
$\Upsilon$. Then, from existing results, one obtains that
$|\mathbb{E}[\epsilon^{2}q^{2}]-\mathbb{E}[\Upsilon^{2}]|$ is $O(\epsilon)$,
which gives a valid bound. From these results, one can obtain bounds in terms
of tail probability of the form $|\mathbb{P}(\epsilon
q>x)-\mathbb{P}(\Upsilon>x)|\leq O(\epsilon)$. However, these are not very
informative as the tail probability itself can be much smaller than
$O(\epsilon)$. Therefore, the rate of convergence of tail probabilities cannot
be obtained using the existing methodologies. In this work, we correctly
characterize $\theta(\epsilon)$ to obtain the rate of convergence of tail
probability to the corresponding HT value. Our results are non-asymptotic in
the sense that they are valid whenever $\epsilon$ is small, and not just when
$\epsilon\rightarrow 0$. Also, our results are precise when $\epsilon$ gets
closer to 0, recovering the HT results.

Our work bridges the gap between the LD, HT, and Many-Server-HT regimes. When
one studies the LD regime, the goal is to find the exponential rate at which
the tail probability decays, which is precisely given by $\theta(\epsilon)$.
As such, our tail bounds can be used to recover the non-asymptotic LD results.
Thus, our tail bounds are at a confluence of non-asymptotic HT and non-
asymptotic LD. This extends the understanding of tail behavior beyond the
classical HT regime. To the best of our knowledge, such comprehensive LD
results have not been previously reported in the existing literature.

###  1.1 Main contribution

We illustrate our methodology by providing results for three well-studied
systems, viz, a load-balancing system under Join the Shortest Queue (hereafter
referred to as the JSQ system), a discrete-time Single-Server Queue (SSQ), and
a multi-server system with a single queue ($M/M/n$ queue). Our contributions
for each of these systems are mentioned below.

####  1.1.1 JSQ system

We consider a continuous-time system with $n$ servers, each of them with its
own queue. Jobs arrive according to a Poisson process with rate $\lambda_{n}$,
and are routed to the server with the shortest queue (breaking ties at
random). Further, the service times are exponentially distributed with rate
$\mu$. In the context of the JSQ system, we consider the Many-Server-HT
regime, where the system size $n$ grows to infinity while the HT parameter
$\epsilon_{n}$ approaches zero. Specifically, we consider
$\epsilon_{n}=n^{-\alpha}$ with $\alpha>1$ constant, and take the limit as
$n\rightarrow\infty$. For the JSQ system, we have the following two
contributions.

  * •

Tail probability:  In Theorem 2, we show that the tail probability of the
steady-state scaled total queue length satisfies

$\displaystyle\frac{1}{1-\epsilon_{n}}e^{-\theta_{n}x}\leq\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}\leq
2ex\Big{(}1+\kappa_{2}n\epsilon_{n}\log\frac{1}{\epsilon_{n}}\Big{)}e^{-\theta_{n}x},$

where $\theta_{n}:=\frac{1}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}$. The
lower bound is valid for all values of $n$ and $\epsilon_{n}$, while the upper
bound holds when the term $n\epsilon_{n}\log\frac{1}{\epsilon_{n}}$ is
sufficiently small. The upper bound leverages the State Space Collapse (SSC)
property of the JSQ system, where the $n$-dimensional state vector collapses
to a one-dimensional subspace, so the JSQ system behaves like an SSQ. The SSC
property holds only when the load is sufficiently large, and so we need
$n\epsilon_{n}\log\frac{1}{\epsilon_{n}}$ to be sufficiently small. This is
satisfied for sufficiently large $n$, when $\epsilon_{n}=n^{-\alpha}$ with
$\alpha>1$.

Figure 1: Visual representation of the condition
$n\epsilon_{n}\log\frac{1}{\epsilon_{n}}<\delta$, where $\delta$ is a small
constant.

Figure 1 visually illustrates the condition discussed above. Region 1
represents values of $n$ and $\alpha$ where
$n\epsilon_{n}\log\frac{1}{\epsilon_{n}}$ is small, and our tail bounds for
the JSQ system hold within this region. As $\alpha$ approaches 1, satisfying
the condition $n\epsilon_{n}\log\frac{1}{\epsilon_{n}}<\delta$ becomes more
challenging, leading to a regime change consistent with existing literature
[1]. In Region 2, we believe that the HT approximation fails to accurately
describe the system dynamics. Finally, Region 3 corresponds to the range
$\alpha\in[0,1]$, where the results and analysis techniques significantly
differ from those presented in this paper.

We also obtain an LD result for the JSQ system, i.e., we have

$\lim_{x\rightarrow\infty}\frac{1}{x}\log\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}=-\theta_{n}.$

Our LD result concerning the scaled total queue length is asymptotically
precise as $\epsilon_{n}$ approaches zero in Many-Server-HT settings (provided
that $\alpha>1$).

Importantly, our work provides an exact characterization of the decay rate of
the tail probabilities for the pre-limit system. This is possible because we
establish both upper and lower bounds on the tail probability. As a result,
our research represents a significant advancement compared to existing
literature. Thus, for JSQ, our result not only bridges the gap between HT and
LD, but also connects these to (some) many server regimes.

  * •

Limiting distribution:  In Theorem 1, we show that as $n\rightarrow\infty$,
all the (scaled) queue lengths become identical and are exponentially
distributed with mean $1$, i.e., in the Many-Server-HT (i.e.,
$\epsilon_{n}=n^{-\alpha}$ with $\alpha>1$), we prove that

$n\epsilon_{n}{\bf q}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\Upsilon\mathbf{1},$

where ${\bf q}$ is the queue length vector in steady state, $\Upsilon$ is an
exponentially distributed random variable, and $\mathbf{1}$ is a vector of all
ones. Prior work [2] establishes such a result only for $\alpha>2$ and leaves
the regime $1<\alpha\leq 2$ open. This result thus fills the gap in the
literature and in conjunction with all the prior work [2], and completes our
understanding of JSQ in all of the Many-Server-HT regimes. A key step to this
end is in establishing SSC for the JSQ system especially in the regimes
$\alpha\in(1,2]$. Our approach in establishing SSC result for the JSQ system
uses a novel Lyapunov function, and differs from the approach based on the
drift arguments of [3] as used in prior work, [4, 5, 6]. Intuitively, instead
of working with the $\ell_{2}$ distance of the queue length vector from the
one dimensional subspace, our key idea is that one should work with the
$\ell_{\infty}$ distance which gives sharper bounds. However, since the
$\ell_{\infty}$ distance is not easily amenable to drift arguments, we work
with a Lyapunov function that can be interpreted as a smooth approximation of
the $\ell_{\infty}$ distance.

####  1.1.2 Discrete time Single-Server Queue

We consider a discrete-time SSQ with general arrival and service
distributions, characterized by mean arrival rate $\lambda_{\epsilon}$ and
mean service rate $\mu_{\epsilon}$, respectively. We introduce the parameter
$\epsilon$ to represent the HT condition, where
$\lambda_{\epsilon}=\mu_{\epsilon}(1-\epsilon)$. Additionally, we define
$\sigma_{\epsilon}^{2}$ as the sum of the variances of the arrival and service
processes.

In Theorem 5, we obtain the following tail bound and LD result for the pre-
limit system:

$\displaystyle\mathbb{P}(\epsilon q>x)\leq
e\theta_{\epsilon}xe^{-\theta_{\epsilon}(1-\kappa_{1}\epsilon)x},$
$\displaystyle\lim_{x\rightarrow\infty}\frac{1}{x}\log\mathbb{P}(\epsilon
q>x)\leq-\theta_{\epsilon}(1-\kappa_{1}\epsilon).$

Here, $\theta_{\epsilon}=\frac{2\mu_{\epsilon}}{\sigma_{\epsilon}^{2}}$ and
$\kappa_{1}$ are defined in terms of the system parameters and moments of the
arrival and service processes. Also, the constant $\kappa_{1}$ depend on the
thrid moments of arrival and service. This result exhibit the influence of a
general arrival and service distribution on the decay rate of the tail
probability under non-asymptotic HT conditions and also provide a non-
asymptotic LD result for the steady-state queue length.

####  1.1.3 $M/M/n$ system

An $M/M/n$ queue is a multi-server system with a single queue, where the
inter-arrival and service times (of each server) are exponentially
distributed, with rates $\lambda_{n}$ and $\mu$, respectively. To parameterize
the Many-Server-HT regime, we use $\lambda_{n}=n\mu(1-\epsilon_{n})$. In this
case, we consider two separate quantities: (i) the number of waiting
customers, denoted by $w_{n}$, and (ii) the number of idle servers, denoted by
$r_{n}$; and provide tail probabilities for both quantities under proper
scaling.

Our approach provides results in three different Many-Server-HT regimes, viz.
Sub-Halfin-Whitt (Sub-HW) regime when $\alpha\in\left(0,\frac{1}{2}\right)$,
Halfin-Whitt (HW) regime when $\alpha=\frac{1}{2}$, and Super-Halfin-Whitt
(Super-HW) regime when $\alpha\in\left(\frac{1}{2},\infty\right)$. A brief
overview of our results for the $M/M/n$ system is presented in Table 1 below.
More details are provided in Section 4.

| Sub-HW  | Halfin-Whitt  | Super-HW
---|---|---|---
$\epsilon_{n}=n^{-\alpha}$ with  |  $\alpha\in\big{(}0,\frac{1}{2}\big{)}$ |  $\alpha=\frac{1}{2}$ | $\alpha>\frac{1}{2}$
$\mathbb{P}(\epsilon_{n}w_{n}>x|w_{n}>0)$ |  $O\big{(}e^{-\theta_{n}x}\big{)}$
$\mathbb{P}(\eta_{n}\tilde{r}_{n}>x|r_{n}>0)$ | $O\big{(}e^{-\frac{1}{2}x^{2}}\big{)}$
$\mathbb{P}(w_{n}>0)$ |  $O\big{(}n^{\alpha-\frac{1}{2}}e^{-n\epsilon_{n}}\big{)}$ |  $\rightarrow p<1$ | $\rightarrow 1$
$\mathbb{P}(r_{n}>0)$ |  $\rightarrow 1$ |  $\rightarrow 1-p<1$ | $O(n^{\frac{1}{2}-\alpha})$
Table 1: A brief overview of our results for the $M/M/n$ queue. Here,
$\tilde{r}_{n}:=r_{n}-n\epsilon_{n}$. The scaling parameter for $w_{n}$ is the
HT parameter $\epsilon_{n}$, while the scaling parameter for $\tilde{r}_{n}$
is $\eta_{n}:=\frac{1}{\sqrt{n(1-\epsilon_{n})}}$. Also, $p$ is a constant
given in Theorem 8.b.

The limiting distribution of the scaled number of waiting customers and the
scaled number of idle servers were established in [7] so we skip mentioning
that here. However, one can easily verify that our results recover the
limiting distribution as $n\rightarrow\infty$. Note that the probability of
there being any idle server, i.e., $\mathbb{P}(r_{n}>0)$ goes to zero in the
Super-HW regime. As such, the tail bound on $r_{n}$ is useful only in Sub-HW
and HW regime. Similarly, the tail bound on $w_{n}$ is valid only in HW and
Super-HW regimes.

###  1.2 Key aspects of our approach

The key idea of this paper is to leverage the existence of the Moment
Generating Function (MGF) of the scaled queue length process in the vicinity
of zero to obtain the decay rate of the tail probability. Let’s consider a
non-negative random variable $X$ and a constant $a>0$ such that
$\mathbb{E}[e^{aX}]\leq\kappa$. By applying Markov’s inequality, we can
establish that $\mathbb{P}(X>x)\leq\kappa e^{-ax}$. This implies that the tail
of $X$ decays exponentially with rate $a$. Thus, determining the appropriate
decay rate of the tail probability reduces to finding the largest $a$ such
that $\mathbb{E}[e^{aX}]$ is bounded by a constant.

We adopt this approach to derive tail bounds for the scaled steady-state
queue-length process. We establish an upper bound on the MGF of the scaled
queue lengths, denoted as $\mathbb{E}[e^{\theta\epsilon q}]$, for a large
range of $\theta$. To determine the correct range of $\theta$, we approximate
the behavior of the actual system by considering its HT counterpart. By doing
so, we also obtain the rate of convergence of the tail probabilities to their
corresponding values in the HT regime.

The use of exponential Lyapunov functions to study queue-length behavior was
introduced as a transform method to obtain HT results in [6]. To get tight
pre-limit results, we use the same exponential test function, but we bound the
error terms in a more refined manner. As such, our technique is divided into
three steps: (i) derive the MGF (or an upper bound) of the scaled queue
length, $\mathbb{E}[e^{\theta\epsilon q}]$, for a large range of values of
$\theta$ , (ii) approximate the MGF with its HT counterpart, and (iii) use
Markov’s inequality and optimize over the values of $\theta$. In the process,
we encounter three multiplicative terms that are required to characterize the
tail probability.

The first term is called the SSC violation. As mentioned before, JSQ system
satisfies SSC in HT. However, in non-asymptotic HT conditions, SSC is not
fully satisfied. As such, the term SSC violation accounts for the level to
which SSC is violated in a pre-limit system. The second one is the pre-limit
tail. When the HT parameter $\epsilon$ is greater than zero, the decay rate of
the tail probability deviates slightly from its HT value, and we refers to
this correct decay rate as the pre-limit tail. The third term is referred to
as the pre-exponent error. Our approach uses Markov’s inequality to compute
bounds tail probabilities from the MGF. However, Markov’s inequality incurs a
cost while obtaining the tail bound, which is captured by the pre-exponent
error.

As a consequence, we easily obtain a LD result after characterizing a bound on
the tail probability. Mathematically, as $\mathbb{E}[e^{aX}]\leq\kappa$
implies $\mathbb{P}(X>x)\leq\kappa e^{-ax}$, we also get
$\lim_{x\rightarrow\infty}\frac{1}{x}\log\mathbb{P}(X>x)\leq-a$. While
characterizing the LD result, the pre-limit error term plays a key role. The
other two error terms, namely the pre-exponent error and the SSC error, vanish
as $x\rightarrow\infty$. Therefore, the pre-limit error term is essential for
accurately characterizing the LD behavior of the system.

###  1.3 Related work

Most of the literature studying the HT behavior of various queues uses a
methodology frequently called ‘diffusion limits.’ A non-exhaustive sample of
articles using this method is [8, 9, 10, 11, 12]. Under this approach, the
scaled system is shown to converge to a Reflected Brownian Motion (RBM) in HT,
and the steady-state behavior of this RBM is studied. The final step is the
so-called interchange-of-limit proof which is usually remarkably challenging.

More recently, there has been an increased number of papers using alternative
methods that do not require the interchange-of-limits step. These are Stein’s
method [13, 14, 15, 16], the BAR approach [17, 18, 19], and the drift method
[4, 6, 20]. All these methods have something in common with the work we
present in this paper and, at the same time, are substantially different.

In the first approach, i.e., the Stein’s method, one derives bounds for the
Wasserstein distance [21] between the pre-limit system and the limiting
distribution, enabling one to obtain the convergence rates. Similar to this,
we also establish bounds that depict the similarity between the pre-limit
system and the limiting distribution. However, our focus lies in directly
characterizing the tail probability. In the second method, i.e., the BAR
approach, the goal is to establish a certain equation, called Basic Adjoint
Relationship (BAR), in terms of the MGF of the steady state distribution of
the limiting distribution. In contrast, we direct work with the MGF of the
pre-limit system.

In the third method, i.e., the drift method, the main idea is to use steady-
state conditions and carefully chosen test functions to compute bounds that
are tight in heavy traffic. Within this method, the use of exponential test
functions contributed to the development of the transform method [6, 20]. Our
work in this paper is inspired by the transform method in the sense that using
exponential test functions and careful manipulation of the queue-length
dynamics are essential to obtain the results. In contrast to the transform
method, which only focuses on the limiting distribution, we carefully compute
the error terms also to obtain the rate of convergence of the tail
probabilities. To the best of our knowledge, the rate of convergence to heavy
traffic in terms of tail probabilities has been not been known before.

$\alpha$ | Regime  | Reference
---|---|---
$\alpha\downarrow 0$ | Mean Field  | [22, 23, 24, 25, 15]
$\alpha\in\big{(}0,\frac{1}{2}\big{)}$ | Sub-Halfin-Whitt  | [26, 27, 28]
$\alpha=\frac{1}{2}$ | Halfin-Whitt  | [29, 30, 31, 32, 33]
$\alpha\in\big{(}\frac{1}{2},1\big{)}$ | Super-Halfin-Whitt  | [34, 35]
$\alpha=1$ | Non-Degenerate Slowdown  | [1]
$\alpha\in(1,2]$ |  Super-Slowdown | This Work
$\alpha>2$ |  | [2]
$\alpha=\infty$ | Classic-HT  | [4, 36, 37, 10, 6, 38, 9]
Table 2: Literature review of asymptotic regimes for Load balancing system.

The literature on Many-Server-HT asymptotics is sub-divided in multiple
categories depending on how fast the load increases with respect to the number
of servers. Using the parameterization $\epsilon_{n}=n^{-\alpha}$ with
$\alpha>0$ introduced above, different regimes are obtained depending on the
value of $\alpha$. The literature on Many-Server-HT with $\alpha\in(0,1]$ is
vast and uses different analysis techniques. We direct the readers to [35] and
references therein for more details on Many-Server-HT with $\alpha\in(0,1]$.
Closest to our work, in [2], it is proved that for $\alpha>2$ the scaled total
queue length is exponentially distribution in limit as $n\rightarrow\infty$.
In our work, we close the gap between the result in [2] and the Non-Degenerate
Slowdown (NDS) [39] regime $(\alpha=1)$, and we obtain the HT behavior for all
$\alpha>1$.

Large Deviations (LD) is a popular regime when one needs to study the
performance of a control policy for routing or scheduling jobs. A
comprehensive report on the application of large deviation theory to queueing
problems can be found in [40]. In [41], it was proved that the if the sequence
of arrivals and services follow a certain Large Deviation Principle (LDP),
then, one can obtain the decay rate of the tail probability of the associated
queue length. This argument was used in [42, 43] to prove a LDP for the queue
length process of a single server process, with further generalization in
[44]. Recently, there has been quite some work on establishing a LDP for the
JSQ system, see [45, 46, 47] and the references therein. In contrast to the
existing literature, our work provides a LD result as a simple closed form
expression in a non-asymptotic LD regime. Notably, we also provide an upper
bound on the pre-exponent term, which is typically absent in previous works.
Moreover, our work considers the JSQ system in a Many-Server regime, and
incorporates the crucial phenomenon of State Space Collapse. Such a phenomenon
was considered in [46], but only for a JSQ system of size $2$.

In addition to proving a LDP, an intriguing research area involves the use of
Lyapunov functions to develop policies that minimize the probability of queue
overflow [48, 49, 50]. The focus in [48, 49] is on minimizing the decay rate
of tail probability, and non-asymptotic large deviation results are not
provided. Another related research field examines queues with Gaussian input,
which are particularly challenging to analyze in most scenarios [51, 52].
These are related to ours in the sense that all compute tight bounds that
characterize the probability of rare events. However, the methodologies are
different, and they focus on the LD regime only. In this paper, we obtain LD
results that are closely connected to the tail probabilities of the queue
lengths in HT and Many-Server-HT.

##  2 Join-the-shortest queue system

###  2.1 Model

We consider a continuous-time queueing system consisting of $n$ SSQs in
parallel, each serving jobs according to first-come-first-serve. At any time
$t$, let ${\bf q}(t)$ denote the queue length vector, where $q_{i}(t)$ is the
queue length of $i^{th}$ queue. For the ease of notation, we use
$\overline{q}(t)$ to denote the total queue length at time $t$, i.e.,
$\overline{q}(t)=\sum_{i=1}^{n}q_{i}(t)$. Jobs arrive to the system according
to a Poisson process with rate $\lambda_{n}$, and service times are
exponentially distributed with rate $\mu$.

When a job arrives, it is dispatched according to JSQ, that is, the job is
sent to the queue with index

$i^{*}(t)\in\arg\min_{i}q_{i}(t),$

breaking ties uniformly at random. Under the JSQ policy, the queue-length
process $\\{{\bf q}(t)\\}_{t\geq 0}$ is a Continuous-Time Markov Chain (CTMC).
Further, it is well-known that the queue-length process is stable (positive
recurrent) if the arrival rate is strictly smaller than the total service rate
$(\lambda_{n}<n\mu)$. In this work, we assume that the system satisfies the
stability condition. Under this assumption, the steady-state distribution of
the queue-length process $\\{{\bf q}(t)\\}_{t=0}^{\infty}$ exists, and we
denote by $\pi_{n}$. We use ${\bf q}$ to denote the steady-state queue length
vector, that is, ${\bf q}$ follows the distribution $\pi_{n}$, and
$\overline{q}:=\sum_{i=1}^{n}q_{i}$.

The system load is $\rho_{n}:=\frac{\lambda_{n}}{n\mu}$, and we define
$\epsilon_{n}=1-\rho_{n}$. Then, the system approaches HT as
$\epsilon_{n}\rightarrow 0$. We consider the JSQ system in HT and Many-Server-
HT. In (classical) HT, the system size $n$ is a constant, and the HT parameter
$\epsilon_{n}$ does not depend on $n$. Then, we drop the subscript $n$ and
take the limit $\epsilon\rightarrow 0$. In the Many-Server-HT, the load and
the number of servers increase together. Then, we consider
$\epsilon_{n}=n^{-\alpha}$ and take the limit as $n\to\infty$. In this work,
in case of Many-Server-HT we only consider the case when $\alpha>1$, which is
also known as super-slowdown regime. In both regimes, we aim to provide a tail
bound on the total queue length, i.e., a bound on
$\mathbb{P}\big{(}\epsilon_{n}\overline{q}>x\big{)}$.

###  2.2 Results for JSQ system

It is well known that the HT distribution of the scaled steady-state total
queue length $\epsilon\overline{q}$ converges to an exponential random
variable as $\epsilon\rightarrow 0$ [6]. Further, as shown in [2], the result
extends to Many-Server-HT, where $\epsilon_{n}\overline{q}$ converges to an
exponential random variable in distribution if $\alpha>2$.

In Theorem 1, we complete the result by demonstrating that
$\epsilon_{n}\overline{q}$ converges in distribution to an exponential random
variable $\alpha>1$. Our result enhances the understanding of the behavior of
the scaled queue length in Many-Server-HT under a broader range of traffic
conditions, encompassing values of $\alpha\in(1,2]$.

######  Theorem 1.

Consider the JSQ system as presented in Section 2.1. Suppose the system
satisfies the condition $\lambda_{n}=n\mu(1-n^{-\alpha})$, i.e.,
$\epsilon_{n}=n^{-\alpha}$, where $\alpha>1$. Then, for any $\theta<1$, we
have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{E}[e^{\theta
n^{-\alpha}\sum_{i=1}^{n}q_{i}}]=\frac{1}{1-\theta},$ $\displaystyle
n^{1-\alpha}{\bf q}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\Upsilon\mathbf{1},\ \text{ as }\ n\rightarrow\infty,$

where $\Upsilon$ is an exponential random variable with mean $1$.

Theorem 1 provides the limiting distribution of the steady-state scaled queue
length vector as $n\rightarrow\infty$ for all $\alpha>1$. The proof of Theorem
1 relies on the fact that the JSQ system satisfies SSC, where the
$n$-dimensional state vector of the system collapses to a one-dimensional
subspace. More precisely, as $n\rightarrow\infty$, we have that
$n^{1-\alpha}q_{i}\approx n^{-\alpha}\sum_{j=1}^{n}q_{j}$ for all
$i\in\\{1,2,\dots,n\\}$. The correct characterization of SSC in the Many-
Server-HT regime is crucial in proving Theorem 1 and in completing the result
of [2] for all $\alpha>1$. More details on SSC for the JSQ system in Many-
Server-HT are provided in Proposition 3 in Section 2.3. A proof sketch for
Theorem 1 is provided in Section 2.4 and the details are provided in Appendix
B.

Next, we provide the tail bound of the scaled steady-state total queue length
for the JSQ system. In Theorem 2, we provide a bound on the tail probability
$\mathbb{P}(\epsilon_{n}\overline{q}>x)$ when the system size is finite.

######  Theorem 2.

Consider the JSQ system as presented in Section 2.1. Suppose the system
satisfies the condition $\lambda_{n}=n\mu(1-\epsilon_{n})$, where
$\epsilon_{n}$ is small enough such that
$n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\frac{\theta_{\perp}}{4}$,
and let $\kappa_{2}:=\frac{4e\kappa_{\perp}^{2}}{\theta_{\perp}}$, where
$\kappa_{\perp}$ and $\theta_{\perp}$ are constants given in Proposition 3.
Suppose $\theta_{n}:=\frac{1}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}$ Then,
for all $x>1-\epsilon_{n}$ we have

$\displaystyle\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}\leq\Big{[}2ex\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)}\Big{]}e^{-\theta_{n}x}.$
(1)

Further, for any $n\geq 1$ and $\epsilon_{n}\in(0,1)$, we have the lower bound

$\displaystyle\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}\geq\frac{1}{1-\epsilon_{n}}e^{-\theta_{n}x}.$
(2)

As a consequence, we have the following large deviation result.

$\displaystyle\lim_{x\rightarrow\infty}-\frac{1}{x}\log\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}=\theta_{n}=\frac{1}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}\in\Big{(}1,\frac{1}{1-\epsilon_{n}}\Big{)}.$
(3)

Theorem 2 establishes the exponential decay of the tail of the total queue
length for a JSQ system. The bounds presented in Eq. (1) and (2) are valid for
both the HT and Many-Server-HT regimes. Further, the result in Theorem 2 is
consistent with the fact that the distribution of the scaled steady-state
total queue length, i.e., $\epsilon_{n}\overline{q}$, converges to an
exponential random variable in distribution, as $n$ grows to $\infty$.

In Theorem 2, we are able to characterize the exact tail decay rate of the
continuous time JSQ system. Our result implies that, in Many-Server-HT with
$\alpha>1$ and when the term
$n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}$ is small enough, the
decay rate of the JSQ system exactly matches the tail decay rate of an SSQ.
This is a significant advancement compared to existing literature. Previous
work primarily focused on comparing the behavior of the JSQ system with an SSQ
in the limiting traffic condition, specifically as $\epsilon_{n}\rightarrow
0$. In contrast, our work examines the behavior of a pre-limit JSQ system and
directly compares it to the corresponding SSQ.

######  Remark 1.

Tail probabilities are often better characterized by multiplicative rather
than additive errors. For instance, a bound of the form
$\Big{|}\mathbb{P}(\epsilon q>x)-\exp\Big{(}-\frac{2\mu
x}{\sigma^{2}}\Big{)}\Big{|}\leq O(\epsilon)$ can become very loose for large
values of $x$. An interesting example is presented in [53, Section 2.1], where
the Central Limit Theorem provides an approximation of $N$ i.i.d. random
variables by a normal distribution, but the error in the approximation can be
as large as $O(1/\sqrt{N})$, which is much larger than the tail itself. Thus,
even if the tail probability converges to an exponentially decaying tail in
the limit, the error in the approximation may still decay slowly. Therefore,
concentration-inequality-type bounds (as in Theorem 5) better characterize the
tail probability.

Additionally, Eq. (3) represents a significant result in the form of an LD
principle. As $\epsilon_{n}\rightarrow 0$ in HT or Many-Server-HT, we have
$\theta_{n}\rightarrow 1$, and so, our LD result for
$\epsilon_{n}\overline{q}$ is asymptotically precise in both, the HT and Many-
Server-HT settings (provided that $\alpha>1$).

####  2.2.1 Discussion on terms in Theorem 2:

Our bounds on tail probability on JSQ system, presented in Theorem 2, can be
decomposed into terms as discussed below.

  * •

SSC violation: For the JSQ system, the SSC violation term is given by
$\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)}$. In non-asymptotic
HT conditions (i.e., when $\epsilon>0$ in HT, or $n<\infty$ in Many-Server-
HT), the SSC property is not fully satisfied. This introduces an additional
multiplicative term in the tail probability bound, which is captured by
$1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}$, and reflects the extent to which
SSC is violated.

To get the SSC violation term below a certain threshold $\delta$ in HT, we
need $\epsilon\log\big{(}\frac{1}{\epsilon}\big{)}$ to be on the order of
$O\big{(}\frac{\delta}{n}\big{)}$. Similarly, in the Many-Server-HT scenario,
it is required that $n^{1-\alpha}\log n$ is on the order of $O(\delta)$, or
alternatively, $n$ must be at least
$\Omega\big{(}\exp\big{(}\frac{1}{\alpha-1}\log\frac{1}{\delta}\big{)}\big{)}$,
or $n\sim\Omega\big{(}\delta^{-\frac{1}{\alpha-1}}\big{)}$ in magnitude.
Satisfying such conditions becomes increasingly challenging as $\alpha$
approaches $1$. This is also shown in Fig. 1 in Section 1.1. Furthermore, it
is important to acknowledge that such conditions cannot be met for $\alpha\leq
1$. These observations provide an intuitive argument for the failure of the
SSC property (as presented in this paper) when $\alpha\leq 1$. Essentially,
for $\alpha\leq 1$, the notion of SSC deviates significantly from the one
considered in this paper.

It is worth noting that although a large system size, represented by $n$, is
required for the SSC error to be small, it remains bounded. Specifically, in
the case of Many-Server-HT, it can be shown that the term
$n\epsilon_{n}\log\frac{1}{\epsilon_{n}}$ is bounded by
$\frac{\alpha}{e(\alpha-1)}$. Consequently, the SSC error is of order at most
$O\left(\frac{1}{\alpha-1}\right)$. Therefore, it is reasonable to expect that
even for moderate system sizes, the tail probability,
$\mathbb{P}(\epsilon\overline{q}>x)$, exhibits exponential decay.

  * •

Pre-limit tail: The pre-limit tail denotes the actual decay rate of the tail
probability of $\epsilon_{n}\overline{q}$ under non-asymptotic HT condition,
i.e., $n<\infty$. For the continuous-time JSQ system, we exactly characterize
the pre-limit tail, which is given by
$\theta_{n}=\frac{1}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}$. In super-
slowdown regime, as $\epsilon_{n}$ goes to zero, the tail of
$\epsilon_{n}\overline{q}$ matches that of an exponential distribution with
mean $1$, as $\lim_{n\rightarrow\infty}\theta_{n}=1$. Further, note that, the
deviation of the pre-limit tail from the corresponding HT value is given by
$|\theta_{n}-1|$, which is of order $O(\epsilon_{n})$.

In the case of continuous-time JSQ, the arrivals and service times follow an
exponential distribution, which is characterized by a single parameter, i.e.,
the arrival or service rate. The specific properties of the exponential
distribution allow for a more accurate analysis of the tail probability.
Intuitively, this distinction is the reason why we are able to precisely
characterize the pre-limit tail for the JSQ system. In case of general arrival
and service distributions, as presented in the case of the discrete-time SSQ,
we do a second-order approximation to approximate the tail probability of the
scaled queue length with its HT counterpart. Consequently, we obtain an upper
bound on the pre-limit tail. More arguments are provided in Section 3.

  * •

Pre-exponent error: In the context of the JSQ system, the pre-exponent error
is represented by the expression $2ex$. This error term arises from using
Markov’s Inequality to obtain tail-probability bounds using MGF. To clarify
this error term, consider a random variable $X$ that follows an exponential
distribution with rate $\lambda$. In this case, the MGF of $X$ is given by
$\mathbb{E}[\exp(\theta X)]=\frac{1}{1-\theta/\lambda}$ for all
$\theta<\lambda$. As shown in Lemma 7 of Appendix A, by Applying Markov’s
Inequality to the MGF and optimizing over the value of $\theta$, we obtain

$\mathbb{P}(X>x)\leq e\lambda xe^{-\lambda x}.$

The upper bound differs from the actual tail of $X$ by a multiplicative factor
of $e\lambda x$, which arises from using Markov’s Inequality. We acknowledge
that it may be possible to eliminate the Markov-Inequality error by employing
more complex techniques. However, we have chosen to rely solely on Markov’s
Inequality for our analysis to maintain simplicity.

###  2.3 State Space Collapse for JSQ system

Next, we mathematically specify SSC for the JSQ system with
$\lambda_{n}=n\mu(1-n^{-\alpha})$ and $\alpha>1$.

######  Proposition 3.

Consider the JSQ system as presented in Section 2.1. Suppose the system
satisfies the condition $\lambda_{n}=n\mu(1-\epsilon_{n})$, and let
$\kappa_{\perp}=128$ and $\theta_{\perp}=1/96$. Then, for $\epsilon_{n}\leq
1/2$, and $\theta\in(0,\theta_{\perp})$, we have that for all
$i\in\\{1,2,\dots,n\\}$,

$\displaystyle\mathbb{E}_{\pi_{n}}\Big{[}\sum_{i=1}^{n}e^{-\theta q_{\perp
i}}+\sum_{i=1}^{n}e^{\theta q_{\perp i}}\Big{]}\leq\kappa_{\perp}n,$
$\displaystyle\mathbb{E}_{\pi_{n}}\big{[}e^{\theta|q_{\perp
i}|}\big{]}\leq\kappa_{\perp},$ (4)

where $q_{\perp i}=q_{i}-\frac{1}{n}\sum_{j=1}^{n}q_{j}$.

As mentioned before, an essential component in establishing the limiting
distribution (Theorem 1) and the tail probability (Theorem 2) is to
characterize SSC accurately. For the JSQ system, when SSC is achieved, the
queue-length vector ${\bf q}$ collapses to a subspace where all its
coordinates become identical, implying that the asymptotic behavior of the JSQ
system closely resembles that of an SSQ.

In Proposition 3, we establish a specific notion of SSC by demonstrating that
the MGF of the deviations of the individual queue lengths and their average,
referred to as the "perpendicular component" and denoted by ${\bf
q}_{\perp}={\bf q}-\frac{1}{n}\overline{q}\mathbf{1}_{n}$, remains uniformly
bounded within an interval around zero for all values of $\epsilon_{n}\leq
1/2$. In particular, this result implies that all the moments of the
perpendicular component remain uniformly bounded, even in the limit as the
system approaches HT conditions ($\epsilon_{n}\rightarrow 0$). In comparison,
the elements of the queue length vector ${\bf q}$ are on the order of
$\frac{1}{\epsilon_{n}}$. Consequently, as $\epsilon_{n}$ becomes small, we
observe that $\epsilon_{n}{\bf
q}\approx\frac{\epsilon_{n}}{n}\overline{q}\,\mathbf{1}_{n}$ (or equivalently,
$\epsilon_{n}{\bf q}_{\perp}\approx\mathbf{0}_{n}$).

##### Intuition behind the choice of the Lyapunov function

In order to establish Proposition 3, we employ the Lyapunov function
$\sum_{i=1}^{n}e^{-\theta q_{\perp i}}+\sum_{i=1}^{n}e^{\theta q_{\perp i}}$,
demonstrating its negative drift. It is worth noting that this Lyapunov
function differs from the one considered in [2], namely
$\exp\big{(}\theta\|{\bf q}_{\perp}\|_{2}\big{)}$, which yields the SSC result
in the Many-Server-HT regime for $\alpha>2$. However, $\exp\big{(}\theta\|{\bf
q}_{\perp}\|_{2}\big{)}$ proves inadequate for establishing SSC in the case of
$\alpha\in(1,2]$. This limitation arises because, to show SSC in Many-Server-
HT for $\alpha\in(1,2]$, it is necessary to obtain a bound for the MGF of
$\|\cdot\|_{\infty}$-norm of the perpendicular component, denoted as $\|{\bf
q}_{\perp}\|_{\infty}$. Intuitively, for an $n$-dim vector $\mathbf{x}$, if
$\|\mathbf{x}\|_{\infty}$ is $O(1)$, then it implies that $\|\mathbf{x}\|_{2}$
is $O(\sqrt{n})$. However, vice-versa need not hold. As such, a bound on
$\ell_{\infty}$ norm is sharper than a bound on $\ell_{2}$ norm, and it plays
a significant role while proving SSC for $\alpha\in(1,2]$. Furthermore, using
$\exp\big{(}\theta\|{\bf q}_{\perp}\|_{\infty}\big{)}$ as the Lyapunov
function proves challenging due to the non-smooth nature of the
$\|\cdot\|_{\infty}$-norm. Consequently, we use the Lyapunov function
$\sum_{i=1}^{n}e^{-\theta q_{\perp i}}+\sum_{i=1}^{n}e^{\theta q_{\perp i}}$
as a suitable alternative. This choice represents a smooth approximation of
$\exp\big{(}\theta\|{\bf q}_{\perp}\|_{\infty}\big{)}$ and simplifies the
technical aspects of the analysis compared to dealing directly with
$\exp\big{(}\theta\|{\bf q}_{\perp}\|_{\infty}\big{)}$. Similar exponential
function was also used in the context of graphical allocation of balls in bins
in [54, 55].

###  2.4 Proof of Theorem 1 and Theorem 2

In this section, provide the proof sketches for Theorem 1 and Theorem 2. More
details and precise mathematical arguments are provided in Appendix B.

######  Lemma 4.

Consider the JSQ system as presented in Section 2.1. Suppose the system
satisfies the condition $\lambda_{n}=n\mu(1-\epsilon_{n})$, where $\alpha>1$.
Let

$\displaystyle\gamma_{n}(\theta):=n\mu-\lambda_{n}e^{\epsilon_{n}\theta}$
$\displaystyle\beta_{n}({\bf
q};\theta):=\mu\sum_{i=1}^{n}\mathds{1}_{\\{q_{i}=0\\}}e^{\theta\epsilon_{n}\sum_{i=1}^{n}q_{i}}.$

Then, for any
$\theta<\theta_{n}:=-\frac{1}{\epsilon_{n}}\log\big{(}1-\epsilon_{n}\big{)}$,
such that $\theta n\epsilon_{n}<\theta_{\perp}$, we have

$\displaystyle\mathbb{E}_{\pi_{n}}\Big{[}e^{\theta\epsilon_{n}\sum_{i=1}^{n}q_{i}}\Big{]}=\frac{1}{\gamma_{n}(\theta)}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf
q};\theta)].$ (5)

Lemma 4 provides the MGF of the scaled total queue length
$\epsilon_{n}\overline{q}$ and is an intermediate step to prove the results in
Theorems 1 and 2. The proof of Lemma 4 is provided in Appendix B.

##### Technical note

It is important to note that the condition $\theta
n\epsilon_{n}<\theta_{\perp}$ in Lemma 4 holds true for all $\theta<0$, and
when $\epsilon_{n}$ is small for $\theta\in(0,\theta_{n})$. Consequently, Eq.
(5) is valid for all $n$ and $\epsilon_{n}$ if $\theta<0$, but if
$\theta\in(0,\theta_{n})$, it can only be used for $\epsilon_{n}$ small
enough. Essentially, we need to use the SSC result given in Proposition 3 to
obtain the bounds in Theorems 1 and 2, and these are valid only under the
condition $\theta n\epsilon_{n}<\theta_{\perp}$.

###### Proof of Theorem 1.

Recall that we use the notation $\lambda_{n}=n\mu(1-\epsilon_{n})$ and
according to the theorem statement, $\epsilon_{n}=n^{-\alpha}$, where
$\alpha>1$. Also, $\overline{q}=\sum_{i=1}^{n}q_{i}$. In order to prove
Theorem 1, we use Lemma 4. We first prove the following claim.

######  Claim 1.

For $n\geq 2$, we have

$\displaystyle\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{]}=\epsilon_{n},\quad\forall
i\in\\{1,2,\dots,n\\}.$ (6)

And, for any $\theta<\theta_{n}$ and $\epsilon_{n}$ small enough such that
$2|\theta|n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\theta_{\perp}$,
we have

$\displaystyle\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{\theta\epsilon_{n}\overline{q}}-1\big{|}\Big{]}\leq\frac{2e\kappa_{\perp}^{2}}{\theta_{\perp}}|\theta|n^{2}\epsilon_{n}^{2}\log\Big{(}\frac{1}{\epsilon_{n}}\Big{)},$

where $\theta_{\perp}$ and $\kappa_{\perp}$ are as defined in Proposition 3.

Proof of Claim 1 is provide in Appendix B. Now, according to the assumption,
$\epsilon_{n}=n^{-\alpha}$, where $\alpha>1$. Under the $\alpha>1$ condition,
the condition
$2|\theta|n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\theta_{\perp}$
can be satisfied for $n$ large enough, i.e., we need $n$ large enough such
that $2\alpha|\theta|n^{1-\alpha}\log n<\theta_{\perp}$, or $\frac{\log
n}{n^{\alpha-1}}<\frac{\theta_{\perp}}{2\alpha|\theta|}$. From Claim 1, we
have that for any $\theta<1$, (as $\theta_{n}\geq 1$), we have,

$\displaystyle\lim_{n\rightarrow\infty}$
$\displaystyle\frac{1}{n\epsilon_{n}}\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{]}=1,$

and

$\displaystyle\lim_{n\rightarrow\infty}$
$\displaystyle\frac{1}{n\mu\epsilon_{n}}\Big{|}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf
q};\theta)]-\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{]}\Big{|}\leq\lim_{n\rightarrow\infty}\frac{1}{n\mu\epsilon_{n}}\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{\theta\epsilon_{n}\overline{q}}-1\big{|}\Big{]}=0.$

Thus, we have
$\lim_{n\rightarrow\infty}\frac{1}{n\mu\epsilon_{n}}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf
q};\theta)]=1$. Further, for any $\theta$, we have

$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n\mu\epsilon_{n}}\gamma_{n}(\theta)$
$\displaystyle=\frac{1}{n\mu\epsilon_{n}}\big{(}n\mu-\lambda_{n}e^{\epsilon_{n}\theta}\big{)}=\frac{1}{\epsilon_{n}}\big{(}1-e^{\epsilon_{n}\theta}+\epsilon_{n}e^{\epsilon_{n}\theta}\big{)}=1-\theta.$

Thus, by using Lemma 4, for any $\theta<1$,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{E}[e^{\theta\epsilon_{n}\overline{q}}]=\lim_{n\rightarrow\infty}\frac{1}{\gamma_{n}(\theta)}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf
q};\theta)]=\frac{1}{1-\theta}.$

This shows that the Moment Generating Function (MGF) of
$\epsilon_{n}\overline{q}$ converges to that of an exponential random variable
with mean $1$. Now, by using [56, Theorem 25.10], we have that
$\epsilon_{n}\overline{q}\stackrel{{\scriptstyle d}}{{\rightarrow}}\Upsilon$,
where $\Upsilon$ is an exponential random variable with mean $1$.

Next, by using Proposition 3, for $\theta\in(0,\theta_{\perp})$ we have that

$\displaystyle
1\leq\lim_{n\rightarrow\infty}\mathbb{E}_{\pi_{\epsilon}}\big{[}e^{n\epsilon_{n}\theta|q_{\perp_{i}}|}\big{]}\leq\lim_{n\rightarrow\infty}\Big{(}\mathbb{E}_{\pi_{\epsilon}}\big{[}e^{\theta|q_{\perp_{i}}|}\big{]}\Big{)}^{n\epsilon_{n}}\leq\lim_{n\rightarrow\infty}\kappa_{\perp}^{n\epsilon_{n}}=1,$

where the second inequality follows by using Jensen’s Inequality. This implies
that, for all $i\in\\{1,2,\dots,n\\}$, we have
$n\epsilon_{n}q_{\perp_{i}}\stackrel{{\scriptstyle d}}{{\rightarrow}}0$. Now,
by definition of ${\bf q}_{\perp}$, we have, $n\epsilon_{n}{\bf
q}=\epsilon_{n}\overline{q}\mathbf{1}+n\epsilon_{n}{\bf q}_{\perp}$, where
$\epsilon_{n}\overline{q}\stackrel{{\scriptstyle d}}{{\rightarrow}}\Upsilon$
and $n\epsilon_{n}q_{\perp i}\stackrel{{\scriptstyle d}}{{\rightarrow}}0$ as
$n\rightarrow\infty$. Thus, we have $n\epsilon_{n}{\bf
q}\stackrel{{\scriptstyle d}}{{\rightarrow}}\Upsilon\mathbf{1}$. This
completes the proof. ∎

###### Proof of Theorem 2.

Note that for $\epsilon_{n}\leq 1/2$, we have $\theta_{n}\leq 2$. Thus, the
condition
$2|\theta|n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\theta_{\perp}$
is satisfied for any $\theta\in(0,\theta_{n})$ whenever $\epsilon_{n}$ is
small enough such that
$4n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\theta_{\perp}$. Then,
from Claim 1, for any $\theta\in(0,\theta_{n})$, we have

$\displaystyle\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf q};\theta)]$
$\displaystyle\leq\mu
n\epsilon_{n}\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)},$

where, as $\theta\leq\theta_{n}\leq 2$, we can take
$\kappa_{2}:=\frac{4e\kappa_{\perp}^{2}}{\theta_{\perp}}$. Now, by using
Markov’s inequality, we have that for all $\theta\in(0,\theta_{n})$,

$\displaystyle\mathbb{P}\big{(}\epsilon_{n}\overline{q}>x\big{)}$
$\displaystyle\leq e^{-\theta
x}\mathbb{E}\big{[}e^{\theta\epsilon_{n}\sum_{i=1}^{n}q_{i}}\big{]}\allowdisplaybreaks$
$\displaystyle\leq\mu
n\epsilon_{n}\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)}\times\frac{1}{\gamma_{n}(\theta)}e^{-\theta
x}\allowdisplaybreaks$ $\displaystyle=\mu
n\epsilon_{n}\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)}\times\frac{1}{n\mu-\lambda_{n}e^{\epsilon_{n}\theta}}e^{-\theta
x}.$ (7)

Next, we optimize the upper bound over the values of
$\theta\in(0,\theta_{n})$. By differentiation,

$\displaystyle\frac{d}{d\theta}\Big{(}\frac{1}{\gamma_{n}(\theta)}e^{-\theta
x}\Big{)}=\frac{1}{\gamma_{n}^{2}(\theta)}\Big{(}-x\gamma_{n}(\theta)+\lambda_{n}\epsilon_{n}e^{\epsilon_{n}\theta}\Big{)}e^{-\theta
x}$

Then, the derivative is equal to zero at

$\displaystyle\theta=\theta_{x,n}:=\frac{1}{\epsilon_{n}}\log\Big{(}\frac{n\mu
x}{\lambda_{n}(x+\epsilon_{n})}\Big{)}=\frac{1}{\epsilon_{n}}\log\Big{(}\frac{x}{(1-\epsilon_{n})(x+\epsilon_{n})}\Big{)}<\theta_{n},$

where the last inequality is valid for any $x>0$. Also, to ensure that
$\theta_{x,n}>0$, we need that

$\displaystyle\frac{x}{(1-\epsilon_{n})(x+\epsilon_{n})}>1\implies
x>1-\epsilon_{n}.$

Then, by substituting $\theta=\theta_{x,n}$, we have

$\displaystyle\frac{1}{n\mu-\lambda_{n}e^{\epsilon_{n}\theta_{x,n}}}e^{-\theta_{x,n}x}$
$\displaystyle=\frac{x+\epsilon_{n}}{\mu
n\epsilon_{n}}\exp\Big{(}-\frac{x}{\epsilon_{n}}\log\Big{(}\frac{x}{(1-\epsilon_{n})(x+\epsilon_{n})}\Big{)}\Big{)}$
$\displaystyle=\frac{x+\epsilon_{n}}{\mu
n\epsilon_{n}}\exp\Big{(}-\frac{x}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}+\frac{x}{\epsilon_{n}}\log\frac{x+\epsilon_{n}}{x}\Big{)}$
$\displaystyle\leq\frac{e(x+\epsilon_{n})}{\mu
n\epsilon_{n}}\exp\Big{(}-\frac{x}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}\Big{)}.$

Using the above result in Eq. (2.4) and $x\geq\epsilon_{n}$, we have

$\displaystyle\mathbb{P}\big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\big{)}$
$\displaystyle\leq
2ex\big{(}1-\kappa_{2}n\epsilon_{n}\log\epsilon_{n}\big{)}\exp\Big{(}-\frac{x}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}\Big{)}.$

This completes the first part of the proof. For the second, i.e., the lower
bound in Eq. (2), we couple the JSQ system with a Single Server Queue (SSQ) as
follows, i.e., we consider a SSQ process $\\{q_{ssq}(t)\\}_{t\geq 0}$, where
all the servers of the original JSQ system are pooled together to create a
single server with service rate $n\mu$. Thus, $\\{q_{ssq}(t)\\}_{t\geq 0}$ is
an $M/M/1$ queue with arrival rate $\lambda_{n}$ and service rate $n\mu$.
Then, the stationary distribution of $\\{q_{ssq}(t)\\}_{t\geq 0}$ is given by,

$\displaystyle\mathbb{P}(q_{ssq}=i)=\Big{(}1-\frac{\lambda_{n}}{n\mu}\Big{)}\Big{(}\frac{\lambda_{n}}{n\mu}\Big{)}^{i}=\epsilon_{n}(1-\epsilon_{n})^{i}.$

This gives us

$\displaystyle\mathbb{P}\Big{(}\epsilon_{n}\sum_{i=1}^{n}q_{i}>x\Big{)}$
$\displaystyle\geq\mathbb{P}(\epsilon_{n}q_{ssq}>x)=(1-\epsilon_{n})^{\big{\lceil}\frac{x}{\epsilon_{n}}\big{\rceil}}\geq(1-\epsilon_{n})^{\frac{x}{\epsilon_{n}}-1}=\frac{1}{1-\epsilon_{n}}\exp\Big{(}-\frac{x}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}\Big{)}.$

This completes the proof. ∎

##  3 Discrete time Single server queue

In this section, we present the bound on the tail probability for a SSQ. The
aim is to show the affect of a general arrival and service distribution on the
decay rate of the tail probability $\mathbb{P}(\epsilon q>x)$. The model is as
follows.

###  3.1 Model for Single Server Queue

We consider a sequence of discrete-time SSQ indexed by the HT parameter
$\epsilon$, where $q(t)$ denotes the number of customers in the system (queue
length) at the beginning of time slot $t$. Customers arrive to the system as
an i.i.d. process $\\{a(t)\\}_{t\geq 0}$, with
$\mathbb{E}[a(t)]=\lambda_{\epsilon}$, Var$(a(t))=\sigma_{\epsilon,a}^{2}$,
and $a(t)\leq A$ almost surely. Once the customers join the queue, the server
decides to serve up to $s(t)$ jobs waiting in the queues, with
$\mathbb{E}[s(t)]=\mu_{\epsilon}$, Var$(s(t))=\sigma_{\epsilon,s}^{2}$, and
$s(t)\leq A$ almost surely. Even though the server can process $s(t)$
customers, it is possible that there are not enough customers in the system,
and some of the service gets wasted. Hence, we call $s(t)$ as the ‘potential’
service. Similar to the arrival process, the potential services are also
independent and identically distributed across time slots. Moreover, we assume
that these potential services are independent of the queue length vector. The
queue evolution process is given by

$q(t+1)=[q(t)+a(t)-s(t)]^{+}=q(t)+a(t)-s(t)+u(t),$ (8)

where $[x]^{+}:=\max\\{x,0\\}$ is used because the queue length cannot be
negative. The term ${u}(t)$ is the unused service and represents the
difference between the potential and actual service. By definition, the unused
service $u(t)$ is positive only if $q(t+1)=0$, which implies $q(t+1)u(t)=0$
for all $t$. Also, the unused service cannot be higher than the potential, and
so, $0\leq u(t)\leq s(t)\leq A$ almost surely.

It is well known that $\lambda_{\epsilon}<\mu_{\epsilon}$ implies stability,
that is, the Markov chain $\\{q(t)\\}_{t\geq 0}$ is positive recurrent. As a
result, the system reaches steady state. We use $\pi_{\epsilon}$ to denote the
steady state distribution of the queue length process $\\{q(t)\\}_{t\geq 0}$.
We drop the symbol $t$ to denote that variables in steady state, i.e., $q$
follows the steady state distribution $\pi_{\epsilon}$, and $q^{+}$ denote the
state that comes after $q$, i.e., $q^{+}=[q+a-s]^{+}=q+a-s+u$, where $a$ and
${s}$ follow the same distribution as $a(t)$ and ${s}(t)$, respectively.

###  3.2 Results for SSQ

In this section, we establish an upper bound on the tail of the queue length
distribution for a system in non-asymptotic HT. We consider
$\lambda_{\epsilon}=\mu_{\epsilon}(1-\epsilon)$, so that as
$\epsilon\rightarrow 0$, the system approaches HT conditions. In Theorem 5, we
provide a bound on the probability $\mathbb{P}(\epsilon q>x)$ when the HT
parameter is bounded away from zero.

######  Theorem 5.

Consider the SSQ as defined in Section 3.1, with
$\lambda_{\epsilon}=\mu_{\epsilon}(1-\epsilon)$. Let
$\sigma_{\epsilon}^{2}:=\sigma_{\epsilon,a}^{2}+\sigma_{\epsilon,s}^{2}$,

$\displaystyle\theta_{\epsilon}:=\frac{2\mu_{\epsilon}}{\sigma_{\epsilon}^{2}},$
$\displaystyle\kappa_{1}:=\frac{2\mu_{\epsilon}E_{3}}{3\sigma_{\epsilon}^{4}}+\epsilon\frac{9\mu_{\epsilon}^{2}A^{4}}{\sigma_{\epsilon}^{6}},$

where $E_{3}:=\max\\{0,\mathbb{E}[(a-s)^{3}]\\}$. Let
$\mathbb{E}[e^{\epsilon\theta q(0)}]<\infty$ for all
$\theta\leq\theta_{\epsilon}$. Then, for all
$x>\frac{\theta_{\epsilon}}{(1+\kappa_{1}\epsilon)}$, we have

$\displaystyle\mathbb{P}(\epsilon q>x)\leq
e\theta_{\epsilon}xe^{-\theta_{\epsilon}(1-\kappa_{1}\epsilon)x},\ \ \ \
\forall\epsilon\in(0,1).$ (9)

Further, as $x\rightarrow\infty$, we have

$\displaystyle\lim_{x\rightarrow\infty}\frac{1}{x}\log\mathbb{P}(\epsilon
q>x)\leq\ -\theta_{\epsilon}(1-\kappa_{1}\epsilon).$ (10)

The exponential decay of the tail of the queue-length distribution of an SSQ
has been established by the tail bound presented in Eq. (9). This result is
consistent with the well-known fact that the distribution of the scaled
steady-state queue length $\epsilon q$ converges to an exponential
distribution as $\epsilon$ tends to zero. However, our result provides a more
detailed insight by characterizing the rate of convergence of the tail
probability to its corresponding HT value.

Note that, SSQ is a one-dimensional system, and so, there is no concept of SSC
in SSQ. Further, in case of SSQ, the pre-exponent error is given by
$e\theta_{\epsilon}x$. The reasoning for this is same as that in the case of
JSQ system. Therefore, we only discuss the pre-limit tail for the discrete
time SSQ.

Pre-limit tail: As mentioned before, pre-limit tail refers to the actual decay
rate of the tail probability for the pre-limit system. Under general arrival
and service distribution, as in the case of discrete time SSQ, we provide an
upper bound on the pre-limit tail, which given by
$\theta_{\epsilon}(1-\kappa_{1}\epsilon)$. The steady-state distribution of
the pre-limit system need not follow an exponential distribution. As such, to
obtain an exponential decay rate on the tail probability, we employ a second
order approximation to approximate the distribution with an exponential. It is
important to note that as $\epsilon$ approaches zero, the upper bound on the
pre-limit tail converges to the correct HT value, i.e.
$\lim_{\epsilon\rightarrow
0}\theta_{\epsilon}(1-\kappa_{1}\epsilon)=\frac{2\mu}{\sigma^{2}}$, where
$\mu$ and $\sigma$ are the limiting values of $\mu_{\epsilon}$ and
$\sigma_{\epsilon}$ as $\epsilon\rightarrow 0$.

Note that the constant $\kappa_{1}$ in the pre-limit tail depends on the third
moment of $(a-s)$, while the actual decay rate of the tail probability in HT
depends on the corresponding first and the second moment. Using exponential
Lyapunov functions, the HT distribution of a queueing system is generally
obtained using a second-order approximation of the MGF of $(a-s)$, which
involves the first two moments. As such, the ‘strength’ of this second-order
approximation depends on the third moment of $(a-s)$, which is essentially
captured by the dependency of the constant $\kappa_{1}$ on $E_{3}$ in Theorem
5.

###  3.3 Proof of Theorem 5

The first step to prove Theorem 5 is characterizing the steady-state
distribution of $\epsilon q$ in terms of its MGF. For this goal, we consider
the Lyapunov function $V(x;\theta):=e^{\theta\epsilon x}$ and obtain the
following result.

######  Lemma 6.

Consider an SSQ as defined in Section 3.1. Suppose the initial state $q(0)$
satisfies $\mathbb{E}[e^{\theta\epsilon q(0)}]<\infty$ for all
$\theta<\theta_{0}$ and for all $\epsilon\in(0,1)$. Then, for any
$\theta<\theta_{0}$, we have

$\displaystyle\mathbb{E}[V(q(t+1);\theta)]$
$\displaystyle=(1-\gamma_{\epsilon}(\theta))^{t+1}\mathbb{E}[V(q(0);\theta)]+\sum_{i=0}^{t}(1-\gamma_{\epsilon}(\theta))^{i}\mathbb{E}[\beta_{\epsilon}(t-i;\theta)],$
(11)

where
$\gamma_{\epsilon}(\theta):=1-\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}$
and $\beta_{\epsilon}(t;\theta):=1-\mathbb{E}\big{[}e^{-\theta\epsilon
u(t)}\big{|}q(t)\big{]}.$ Let
$\Theta_{\epsilon}:=\\{\theta\in\mathbb{R}:\theta<\theta_{0},\
\mathbb{E}[e^{\theta\epsilon(a-s)}]<1\\}$. Then, for all
$\theta\in\Theta_{\epsilon}\cup\\{\theta:\theta\leq 0\\}$, we have

$\displaystyle\mathbb{E}_{\pi_{\epsilon}}[V(q;\theta)]:=\mathbb{E}_{\pi_{\epsilon}}\big{[}e^{\theta\epsilon
q}\big{]}=\frac{1-\mathbb{E}_{\pi_{\epsilon}}[e^{-\theta\epsilon
u}]}{1-\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}}.$ (12)

The expression in Eq. (12) was first introduced in [6] for $\theta$ in an
interval around zero. In this work, we have expanded the range of validity for
Eq. (12) to a much larger interval given by $\Theta_{\epsilon}$. Although it
may be difficult to fully characterize $\Theta_{\epsilon}$ for a fixed
$\epsilon>0$, we can use Taylor’s approximation of the exponential function to
obtain a close enough subset. Accurately characterizing $\Theta_{\epsilon}$ is
crucial to obtaining the correct tail bounds. Our finite-time characterization
of the queue-length distribution, as shown in Eq. (11), allows us to achieve
this expanded range of validity.

The proof of Lemma 6 follows by showing that, for any
$\theta\in\Theta_{\epsilon}$, the right-hand side (RHS) in Eq. (11) converges
to a finite value. The details of the proof of Lemma 6 are presented in
Appendix C.

######  Lemma 7.

Suppose $X$ is a non-negative random variable, and the MGF of $X$ satisfies
the inequality $\mathbb{E}[e^{\theta
X}]\leq\frac{1}{1-\theta/\lambda},\forall\theta\in(0,\lambda)$, where
$\lambda>0$. Then, for any $x>\frac{1}{\lambda}$, we have

$\displaystyle\mathbb{P}(X>x)\leq e\lambda xe^{-\lambda x}.$

The proof of Lemma 7 is provided in Appendix A.

###### Proof of Theorem 5.

By equating the drift of the queue length $q(t)$ to zero in steady state, we
prove the following claim.

######  Claim 2.

For any $\theta\in\mathbb{R}$,
$1-\mathbb{E}_{\pi_{\epsilon}}[e^{-\theta\epsilon
u}]\leq\epsilon^{2}\theta\mu_{\epsilon}.$

Next, using Taylor’s expansion of the exponential function we obtain the
following claim.

######  Claim 3.

For any
$\theta\in\Big{(}0,\frac{\theta_{\epsilon}}{(1+\kappa_{1}\epsilon)}\Big{)}$,
$1-\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}\geq\theta\epsilon^{2}\mu_{\epsilon}\Big{(}1-\frac{\theta}{\theta_{\epsilon}}\big{(}1+\kappa_{1}\epsilon\big{)}\Big{)}.$

The proofs of Claims 2 and 3 are provided in Appendix C. Next, we use the
Lemma 6. Note that, by the assumption $\mathbb{E}[e^{\theta\epsilon
q(0)}]<\infty$ for $\theta<\theta_{\epsilon}$, we have that
$\Big{(}0,\frac{\theta_{\epsilon}}{(1+\kappa_{1}\epsilon)}\Big{)}\subseteq\Theta_{\epsilon}$.
Then, by Lemma 6 and the above mentioned claims, we have

$\displaystyle\mathbb{E}\big{[}e^{\epsilon\theta
q}\big{]}\leq\Big{(}1-\frac{\theta}{\theta_{\epsilon}}\Big{(}1+\kappa_{1}\epsilon\Big{)}\Big{)}^{-1},\quad\forall\theta\in\Big{(}0,\frac{\theta_{\epsilon}}{(1+\kappa_{1}\epsilon)}\Big{)}.$

Afterwards, the result in Theorem 5 follows simply by using Markov’s
inequality as shown in Lemma 7.∎

##  4 $M/M/n$ system

###  4.1 Model

An $M/M/n$ system (also known as Erlang-C) is a multi-server queue with $n$
servers, a single queue, and follows a continuous-time first-come-first-serve
service discipline. Customer arrivals follow a Poisson process with rate
$\lambda_{n}$, and each server has an exponentially distributed service time
with constant service rate $\mu$. The rate $\mu$ is independent of the system
size; thus, the subscript $n$ is omitted.

To represent the system’s dynamics, we use $q_{n}(t)$ to denote the number of
customers in the system at time $t$. We use $w_{n}(t)=[q_{n}(t)-n]^{+}$ to
denote the number of customers waiting in the queue, and
$r_{n}(t)=[n-q_{n}(t)]^{+}$ for the number of idle servers at time $t$.

The queue-length process $\\{q_{n}(t)\\}_{t\geq 0}$ is a CTMC, and it is well
known that it is stable (positive recurrent) if $\lambda_{n}<n\mu$. To be
consequent with the previous sections, we assume that
$\lambda_{n}=n\mu\left(1-\epsilon_{n}\right)$. Further, a stationary
distribution $\pi_{n}$ exists under this stability condition. Similarly to
previous sections, we drop the index $t$ to denote steady-state variables.

Similarly to the JSQ system, to study the Many-Server-HT asymptotics of this
system, we define the load $\rho_{n}:=\frac{\lambda_{n}}{n\mu}=1-\epsilon_{n}$
and consider $\epsilon_{n}=cn^{-\alpha}$, where $c>0$ is a constant, and
$\alpha>0$ a parameter (as in the JSQ system). In this case, we are interested
in $\alpha\in(0,1]$.

As explained in the introduction, the steady-state dynamics of the system vary
with the value of $\alpha$. We study three regimes: (i) Sub-HW, where
$\alpha\in\Big{(}0,\frac{1}{2}\Big{)}$, (ii) HW, where $\alpha=\frac{1}{2}$,
and (iii) Super-HW where $\alpha\in\Big{(}\frac{1}{2},\infty\Big{)}$.

###  4.2 Results on $w_{n}$ and $r_{n}$

In this section, we provide separate tail bounds on the number of idle servers
and the number of waiting customers because the scaling parameter is
different. In Theorem 8 we present our results for the number of idle servers
$r_{n}$.

######  Theorem 8.

Consider the M/M/n queue as given in Section 4.1 with
$\lambda_{n}=n\mu(1-\epsilon_{n})$. Let $\tilde{r}_{n}:=r_{n}-n\epsilon_{n}$,
and $\eta_{n}:=\frac{1}{\sqrt{n(1-\epsilon_{n})}}$.

  1. (a)

Super-HW regime: Suppose $\epsilon_{n}=cn^{-\alpha}$ with
$\alpha>\frac{1}{2}$. Then, if $n^{2\alpha-1}>4c^{2}$, we have

$\displaystyle\mathbb{P}(r_{n}>0)\leq 4e\pi cn^{-\alpha+\frac{1}{2}},$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(r_{n}>0)=0,$ (13)

and there exist a constant $\kappa_{1}$, independent of $n$, such that,

$\displaystyle\mathbb{P}\big{(}\eta_{n}\tilde{r}_{n}>x\big{)}\leq\frac{\kappa_{1}c}{n^{\alpha-\frac{1}{2}}}e^{-\frac{1}{2}x^{2}}.$

  2. (b)

HW regime: Suppose $\epsilon_{n}=cn^{-\frac{1}{2}}$. Then,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(r_{n}>0)=\frac{\sqrt{2\pi}c\exp\left(\frac{c^{2}}{2}\right)\Phi\left(c\right)}{1+\sqrt{2\pi}c\exp\left(\frac{c^{2}}{2}\right)\Phi\left(c\right)},$
(14)

and there exist constant $\kappa_{2}$, such that

$\displaystyle\mathbb{P}\big{(}\eta_{n}\tilde{r}_{n}>x\big{)}\leq\kappa_{2}e^{-\frac{1}{2}x^{2}}.$

  3. (c)

Sub-HW regime: Suppose $\epsilon_{n}=cn^{-\alpha}$ with
$\alpha\in\big{(}0,\frac{1}{2}\big{)}$. Then, there exists constant
$\kappa_{3}$ such that

$\displaystyle\mathbb{P}(r_{n}>0)\geq
1-\frac{\kappa_{3}}{c}n^{\alpha-\frac{1}{2}}e^{-cn^{\frac{1}{2}-\alpha}},$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(r_{n}>0)=1,$

and we have

$\displaystyle\mathbb{P}(\eta_{n}\tilde{r}_{n}>x)\leq e^{-\frac{1}{2}x^{2}},$
$\displaystyle\mathbb{P}(\eta_{n}\tilde{r}_{n}<-x)\leq
e^{-\frac{1}{2}x^{2}+2e\epsilon_{n}x^{2}}.$

Theorem 8 provides the tail bounds on the number of idle servers in the
$M/M/n$ queue in each of the three regimes considered in the paper.

######  Remark 2.

Note that we do not provide a LD result for $r_{n}$, as we did in the case of
SSQ and JSQ systems, as the number of idle servers $r_{n}$ is always bounded
by $n$. As such, for the pre-limit system (i.e., when the system size $n$ is
finite), we have that

$\displaystyle\mathbb{P}(\eta_{n}\tilde{r}_{n}>x)=0,\ \forall
x>\sqrt{n\rho_{n}}.$

Next, we provide the results for the number of waiting customers $w_{n}$.

######  Corollary 9.

Consider the M/M/n queue as given in Section 4.1 with
$\lambda_{n}=n\mu(1-\epsilon_{n})$. Then, for any $\\{\epsilon_{n}\\}_{n\geq
0}$ such that $\epsilon_{n}\rightarrow 0$ as $n\rightarrow\infty$, and for any
$\theta<1$, we have

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{E}\big{[}e^{\theta\epsilon_{n}w_{n}}\big{|}w_{n}>0\big{]}=\frac{1}{1-\theta},$
$\displaystyle[\epsilon_{n}w_{n}|w_{n}>0]\stackrel{{\scriptstyle
d}}{{\rightarrow}}\Upsilon,$

where $\Upsilon$ is an exponential random variable with mean $1$. Further, for
any $x>\epsilon_{n}$, we have

$\displaystyle\frac{1}{(1-\epsilon_{n})^{2}}e^{-\theta_{n}x}\leq\mathbb{P}(\epsilon_{n}w_{n}\geq
x|w_{n}>0)\leq\frac{1}{1-\epsilon_{n}}e^{-\theta_{n}x},$

where $\theta_{n}=\frac{1}{\epsilon_{n}}\log\frac{1}{1-\epsilon_{n}}$. As a
consequence, we get the LD result

$\displaystyle\lim_{x\rightarrow\infty}\frac{1}{x}\mathbb{P}(\epsilon_{n}w_{n}\geq
x)=\lim_{x\rightarrow\infty}\frac{1}{x}\mathbb{P}(\epsilon_{n}w_{n}\geq
x|w_{n}>0)=\theta_{n}.$ (15)

Corollary 9 provides the limiting distribution of the steady state number of
waiting customers in the system. For any $n$ and $\epsilon_{n}$, the
distribution of $[w_{n}|w_{n}>0]$ can be exactly characterized. As will be
depicted in Lemma 10, it turns out that $[w_{n}|w_{n}>0]$ follows a geometric
distribution with parameter $\epsilon_{n}$. Such a result can also be derived
by simply solving for the steady-state distribution of an $M/M/n$ queue.
Furthermore, as we can exactly characterize the distribution of
$[w_{n}|w_{n}>0]$, there is no need for using Markov’s inequality to get a
tail bound on the number of waiting customers. Hence, we do not obtain the
pre-limit nor the pre-exponent error in the case of the number of waiting
customers.

######  Remark 3.

It is worth noting that although the conditional distribution
$[w_{n}|w_{n}>0]$ follows a geometric distribution, it is possible for the
probability $\mathbb{P}(w_{n}>0)$ to approach zero. This can be observed by
examining the bounds on $\mathbb{P}(r_{n}>0)$ provided in Theorem 8.
Specifically, in the Sub-HW regime, we have

$\displaystyle\mathbb{P}(w_{n}>0)\leq
1-\mathbb{P}(r_{n}>0)\leq\frac{\kappa_{3}}{c}n^{\alpha-\frac{1}{2}}e^{-cn^{\frac{1}{2}-\alpha}}.$

It is well-known that $\mathbb{P}(w_{n}>0)$ decreases to zero in this regime.
However, we additionally characterize the rate at which this convergence
occurs.

######  Remark 4.

The results presented in Theorem 8 and Corollary 9 only consider Many-Server-
HT regimes. There, the HT parameter $\epsilon_{n}$ approaches zero as the
system size $n$ grows. One can also consider a conventional HT regime, where
the system size $n$ is constant and the HT parameter $\epsilon\rightarrow 0$
(using the notation $\lambda_{n}=n\mu(1-\epsilon)$) independent of $n$. In
this case, the tail bound for the relevant quantities closely aligns with the
tail-bound results obtained in the Super-HW regime. By employing similar
calculations as, it can be shown that the number of waiting customers exhibits
the same tail bounds as presented in Corollary 9 after replacing
$\epsilon_{n}$ by $\epsilon$. Furthermore, the probability
$\mathbb{P}(r_{n}>0)$ is on the order of $O(\epsilon)$. Moreover, in
conventional HT, there is no need for a tail bound on $r_{n}$ as $r_{n}\leq n$
and $n$ is a constant.

####  4.2.1 Steady-state distribution of $w_{n}$ and $r_{n}$

In this section, we provide the complete characterization of the steady state
distribution of $\\{q_{n}(t)\\}_{t\geq 0}$. Specifically, we provide the MGF
of the steady-state number of waiting customers ($w_{n}$) and idle servers
($r_{n}$).

######  Lemma 10.

Consider the M/M/n queue as given in Section 4.1 and suppose
$\lambda_{n}=n\mu(1-\epsilon_{n})$. Define

$\displaystyle G_{n}(t)$
$\displaystyle:=\exp\Big{(}-n\epsilon_{n}t-n(1-\epsilon_{n})(e^{-t}+t-1)\Big{)}.$

Then, we have the following results.

  1. (a)

For any $\theta\in\mathbb{R}$, we have

$\displaystyle\mathbb{E}\left[e^{\theta
r_{n}}\Big{|}r_{n}>0\right]=G_{n}^{-1}(\theta)\left(\int_{-\infty}^{0}G_{n}dt\right)^{-1}\int_{-\infty}^{\theta}G_{n}dt.$

  2. (b)

For any $\theta<\log\frac{1}{1-\epsilon_{n}}$, we have

$\displaystyle\mathbb{E}\left[e^{\theta
w_{n}}\Big{|}w_{n}>0\right]=\frac{1}{1-\frac{1}{\epsilon_{n}}\left(1-e^{-\theta}\right)}.$

Further, we have
$\mathbb{P}(q_{n}=n)=\left(\frac{1}{\epsilon_{n}}+n\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}$,
and

$\displaystyle\mathbb{P}(w_{n}>0)=\frac{1-\epsilon_{n}}{\epsilon_{n}}\mathbb{P}(q_{n}=n),$
$\displaystyle\mathbb{P}(r_{n}>0)=n\mathbb{P}(q_{n}=n)\int_{-\infty}^{0}G_{n}(t)dt.$

Lemma 10 provides the steady state distribution of $r_{n}$ and $w_{n}$ for any
value of $n$ and $\epsilon_{n}$. The results in Theorem 8 are derived by using
the result in Lemma 10.a, Markov’s inequality. Lemma 10.a indicates that the
distribution of $[r_{n}|r_{n}>0]$ closely resembles that of a truncated normal
random variable. By replacing $G_{n}(t)$ with
$\tilde{G}_{n}(t):=\exp\big{(}-n\eta_{n}\epsilon_{n}t-\frac{1}{2}t^{2}\big{)}$,
we have that the function

$\displaystyle\tilde{G}_{n}^{-1}(\theta)\left(\int_{-\infty}^{0}\tilde{G}_{n}dt\right)^{-1}\int_{-\infty}^{\theta}\tilde{G}_{n}dt$

defines the MGF of a truncated normal random variable, specifically
$[Y_{n}|Y_{n}>0]$, where $Y_{n}$ follows a normal distribution with mean
$n\eta_{n}\epsilon_{n}$ and variance $1$.

Intuitively, for large values of $n$, we can approximate $G_{n}(\eta_{n}t)$ by
$\tilde{G}_{n}(t)$ because
$\frac{1}{\eta_{n}^{2}}\big{(}e^{-\eta_{n}t}+\eta_{n}t-1\big{)}\approx\frac{t^{2}}{2}$,
where $\eta_{n}=\frac{1}{\sqrt{n\rho_{n}}}$. Consequently, we observe that,
for large $n$, the distribution of $[\eta_{n}r_{n}|r_{n}>0]$ closely matches
the distribution of $[Y_{n}|Y_{n}>0]$, where $Y_{n}$ is defined above.
Consequently, for large values of $n$, $[\eta_{n}\tilde{r}_{n}|r_{n}>0]$
closely matches the distribution of a truncated standard normal random
variable. In the proof of Theorem 8, use this idea to establish tail bounds on
$r_{n}$ in each of the three regimes.

The result in Corollary 9 immediately follows from Lemma 10.b just be
observing that the MGF of $[w_{n}|w_{n}>0]$ matches with that of a geometric
random variable with parameter $\epsilon_{n}$. The mathematical details for
results provided in this section (Section 4.2) is provided in Appendix D.

## References

  * [1] V. Gupta, N. Walton, Load balancing in the Nondegenerate Slowdown Regime, Operations Research 67 (1) (2019) 281–294.
  * [2] D. Hurtado-Lange, S. T. Maguluri, A load balancing system in the many-server heavy-traffic asymptotics, Queueing Systems 101 (3-4) (2022) 353–391.
  * [3] B. Hajek, Hitting-time and occupation-time bounds implied by drift analysis with applications, Advances in Applied Probability (1982) 502–525.
  * [4] A. Eryilmaz, R. Srikant, Asymptotically tight steady-state queue length bounds implied by drift conditions, Queueing Systems 72 (3-4) (2012) 311–359.
  * [5] S. T. Maguluri, R. Srikant, Heavy traffic queue length behavior in a switch under the MaxWeight algorithm, Stochastic Systems 6 (1) (2016) 211–250.
URL http://dx.doi.org/10.1214/15-SSY193

  * [6] D. Hurtado-Lange, S. T. Maguluri, Transform methods for heavy-traffic analysis, Stochastic Systems 10 (4) (2020) 275–309.
  * [7] A. Braverman, J. Dai, J. Feng, Stein’s method for steady-state diffusion approximations: An introduction through the Erlang-A and Erlang-C models, Stochastic Systems 6 (2) (2017) 301–366.
  * [8] J. Kingman, On queues in heavy traffic, Journal of the Royal Statistical Society. Series B (Methodological) (1962) 383–392.
  * [9] J. M. Harrison, Heavy traffic analysis of a system with parallel servers: asymptotic optimality of discrete-review policies, Annals of applied probability (1998) 822–848.
  * [10] R. Williams, On dynamic scheduling of a parallel server system with complete resource pooling, Fields Institute Communications 28 (49-71) (2000) 5–1.
  * [11] J. M. Harrison, A. Zeevi, Dynamic scheduling of a multiclass queue in the halfin-whitt heavy traffic regime, Operations Research 52 (2) (2004) 243–257.
  * [12] D. Gamarnik, A. Zeevi, Validity of heavy traffic steady-state approximations in Generalized Jackson Networks, The Annals of Applied Probability (2006) 56–90.
  * [13] A. Braverman, J. Dai, Stein’s method for steady-state diffusion approximations of M/Ph/n+M systems, The Annals of Applied Probability 27 (1) (Feb. 2017).  doi:10.1214/16-AAP1211.
URL https://doi.org/10.1214/16-AAP1211

  * [14] A. Braverman, J. Dai, Stein’s method for steady-state diffusion approximations of M/Ph/n+ M systems, The Annals of Applied Probability 27 (1) (2017) 550–581.
  * [15] L. Ying, Stein’s method for mean field approximations in light and heavy traffic regimes, Proc. ACM Meas. Anal. Comput. Syst. 1 (1) (2017) 12:1–12:27.  doi:10.1145/3084449.
URL http://doi.acm.org/10.1145/3084449

  * [16] R. Gaunt, N. Walton, Stein’s method for the single server queue in heavy traffic, Statistics & Probability Letters 156 (2020) 108566.
  * [17] A. Braverman, J. Dai, M. Miyazawa, Heavy traffic approximation for the stationary distribution of a Generalized Jackson Network: The BAR approach, Stochastic Systems 7 (1) (2017) 143–196.
  * [18] A. Braverman, J. Dai, M. Miyazawa, The bar-approach for multiclass queueing networks with sbp service policies, arXiv preprint arXiv:2302.05791 (2023).
  * [19] J. Dai, P. Glynn, Y. Xu, Asymptotic product-form steady-state for generalized jackson networks in multi-scale heavy traffic, arXiv preprint arXiv:2304.01499 (2023).
  * [20] P. R. Jhunjhunwala, M. Zubeldia, S. T. Maguluri, Join-the-shortest queue with abandonment: Critically loaded and heavily overloaded regimes, arXiv preprint arXiv:2211.15050 (2022).
  * [21] S. Vallender, Calculation of the wasserstein distance between probability distributions on the line, Theory of Probability & Its Applications 18 (4) (1974) 784–786.
  * [22] D. Mukherjee, S. C. Borst, J. S. Van Leeuwaarden, P. A. Whiting, Universality of power-of-d load balancing in many-server systems, Stochastic Systems 8 (4) (2018) 265–292.
  * [23] M. Mitzenmacher, Load balancing and density dependent jump Markov processes, in: focs, IEEE, 1996, p. 213.
  * [24] M. Mitzenmacher, The power of two choices in randomized load balancing, IEEE Transactions on Parallel and Distributed Systems 12 (10) (2001) 1094–1104.
  * [25] A. Stolyar, Pull-based load distribution among heterogeneous parallel servers: The case of multiple routers, Queueing Systems 85 (1-2) (2017) 31–65.
  * [26] X. Liu, L. Ying, A simple steady-state analysis of load balancing algorithms in the sub-Halfin-Whitt regime, ACM SIGMETRICS Performance Evaluation Review 46 (2) (2019) 15–17.
  * [27] S. M. Varma, F. Castro, S. T. Maguluri, Power-of-$d$ choices load balancing in the sub-halfin whitt regime (2022).  arXiv:2208.07539.
  * [28] S. Bhamidi, A. Budhiraja, M. Dewaskar, Near equilibrium fluctuations for supermarket models with growing choices, The Annals of Applied Probability 32 (3) (2022) 2083–2138.
  * [29] S. Banerjee, D. Mukherjee, Join-the-shortest queue diffusion limit in Halfin–Whitt regime: Tail asymptotics and scaling of extrema, The Annals of Applied Probability 29 (2) (2019) 1262–1309.
  * [30] S. Banerjee, D. Mukherjee, Join-the-shortest queue diffusion limit in Halfin–Whitt regime: Sensitivity on the heavy-traffic parameter, The Annals of Applied Probability 30 (1) (2020) 80–144.
  * [31] A. Braverman, Steady-state analysis of the join-the-shortest-queue model in the Halfin–Whitt regime, Mathematics of Operations Research (2020).
  * [32] P. Eschenfeldt, D. Gamarnik, Join the shortest queue with many servers. the heavy-traffic asymptotics, Mathematics of Operations Research 43 (3) (2018) 867–886.
  * [33] S. Halfin, W. Whitt, Heavy-traffic limits for queues with many exponential servers, Operations research 29 (3) (1981) 567–588.
  * [34] X. Liu, L. Ying, Universal scaling of distributed queues under load balancing in the super-halfin-whitt regime, IEEE/ACM Transactions on Networking 30 (1) (2021) 190–201.
  * [35] Z. Zhao, S. Banerjee, D. Mukherjee, Many-server asymptotics for join-the-shortest queue in the super-halfin-whitt scaling window (2023).  arXiv:2106.00121.
  * [36] G. Foschini, J. Salz, A basic dynamic routing problem and diffusion, IEEE Transactions on Communications 26 (3) (1978) 320–327.
  * [37] R. Williams, Diffusion approximations for open multiclass queueing networks: Sufficient conditions involving state space collapse, Queueing Systems Theory and Applications (1998) 27 – 88.
  * [38] J. M. Harrison, M. J. López, Heavy traffic resource pooling in parallel-server systems, Queueing systems 33 (4) (1999) 339–368.
  * [39] V. Gupta, N. Walton, Load balancing in the nondegenerate slowdown regime, Operations Research 67 (1) (2019) 281–294.
  * [40] A. J. Ganesh, Big queues, Springer Science & Business Media, 2004.
  * [41] P. W. Glynn, W. Whitt, Logarithmic asymptotics for steady-state tail probabilities in a single-server queue, Journal of Applied Probability 31 (A) (1994) 131–156.
  * [42] A. Puhalskii, Large deviation analysis of the single server queue, Queueing Systems 21 (1995) 5–66.
  * [43] N. G. Duffield, N. O’connell, Large deviations and overflow probabilities for the general single-server queue, with applications, in: Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 118, Cambridge University Press, 1995, pp. 363–374.
  * [44] P. Dupuis, R. S. Ellis, The large deviation principle for a general class of queueing systems. i, Transactions of the American Mathematical Society 347 (8) (1995) 2689–2751.
  * [45] A. Ridder, A. Shwartz, Large deviations without principle: Join the shortest queue, Mathematical Methods of Operations Research 62 (2005) 467–483.
  * [46] R. D. Foley, D. R. McDonald, Join the shortest queue: stability and exact asymptotics, Annals of Applied Probability (2001) 569–607.
  * [47] A. A. Puhalskii, A. A. Vladimirov, A large deviation principle for join the shortest queue, Mathematics of Operations Research 32 (3) (2007) 700–710.
  * [48] X. Lin, V. Venkataramanan, On the large-deviations optimality of scheduling policies minimizing the drift of a lyapunov function, in: 2009 47th Annual Allerton Conference on Communication, Control, and Computing (Allerton), IEEE, 2009, pp. 919–926.
  * [49] V. Venkataramanan, X. Lin, On the queue-overflow probability of wireless systems: A new approach combining large deviations with lyapunov functions, IEEE transactions on information theory 59 (10) (2013) 6367–6392.
  * [50] S. Banerjee, Y. Kanoria, P. Qian, Large deviations optimal scheduling of closed queueing networks, arXiv preprint arXiv:1803.04959 (2018).
  * [51] M. Zubeldia, M. Mandjes, Large deviations for acyclic networks of queues with correlated gaussian inputs, Queueing Systems 98 (3-4) (2021) 333–371.
  * [52] M. Mandjes, Large deviations for Gaussian queues: modeling communication networks, John Wiley & Sons, 2007.
  * [53] R. Vershynin, High-dimensional probability: An introduction with applications in data science, Vol. 47, Cambridge university press, 2018.
  * [54] Y. Peres, K. Talwar, U. Wieder, Graphical balanced allocations and the (1+ $\beta$)-choice process, Random Structures & Algorithms 47 (4) (2015) 760–775.
  * [55] N. Bansal, O. N. Feldheim, The power of two choices in graphical allocation, in: Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, 2022, pp. 52–63.
  * [56] P. Billingsley, Probability and Measure, 2nd Edition, John Wiley and Sons, 1986\.

##  Appendix A Essential Lemmas and their proofs

###### Proof of Lemma 7.

The proof simply follows by using the Markov’s inequality and then optimizing
over the value of $\theta$. For any $\theta\in(0,1/\lambda)$, we have

$\displaystyle\mathbb{P}(X>x)\leq e^{-\theta x}\mathbb{E}[e^{\theta
X}]\leq\frac{1}{1-\theta/\lambda}e^{-\theta x}.$ (16)

As the above inequality is valid for any $\theta\in(0,\lambda)$, we can
minimize the RHS in the above inequality over the range $(0,\lambda)$. Now,
the derivative of the RHS equates to zero when
$\theta=\theta^{*}:=\lambda\Big{(}1-\frac{1}{\lambda x}\Big{)}$, where
$\theta^{*}\in(0,\lambda)$ whenever $x>1/\lambda$. Thus, for any
$x>1/\lambda$, by substituting $\theta=\theta^{*}$ in Eq. (16), we get the
result. ∎

######  Lemma 11.

Consider a Discrete Time Markov Chain $\\{\mathbf{x}(t)\\}_{t\geq 0}$, where
$\mathbf{x}(t)\in\mathbb{R}^{d}$ for all $t\geq 0$. Let the Markov chain be
positive recurrent, with stationary distribution $\pi$, and
$\mathbf{x}_{\infty}$ be a random variable that follows the distribution
$\pi$. Suppose $f:\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}_{+}$ is
a Lyapunov function such that following conditions hold.

  * C1.

There exists $\theta_{0}>0$ such that $\forall\theta<\theta_{0}$, we have,

$\displaystyle\mathbb{E}[f(\mathbf{x}(t),\theta)]<\infty,\ \ \forall t\geq 0.$

  * C2.

There exists $\theta_{1}>0$ such that for all $\theta<\theta_{1}$, the
function $f$ satisfies the negative drift condition,

$\displaystyle\mathbb{E}[\Delta
f(\mathbf{x}(t+1),\theta)|\mathbf{x}(t)]=-\gamma(\theta)f(\mathbf{x}(t),\theta)+\beta(\mathbf{x}(t),\theta),$
(17)

where for all $\theta<\theta_{1}$, one of the following conditions is
satisfied

    1. 1.

There exists $\bar{\beta}(\theta)$ such that

$\displaystyle\gamma(\theta)>0,$
$\displaystyle\mathbb{E}[\beta(\mathbf{x}(t),\theta)]<\bar{\beta}(\theta),$
$\displaystyle\lim_{t\rightarrow\infty}\mathbb{E}[\beta(\mathbf{x}(t),\theta)]\
\text{ exists}.$

    2. 2.

There exists a function $b(\theta)<\infty$ such that

$\displaystyle\mathbb{E}[f(\mathbf{x}(t),\theta)]\leq b(\theta),\ \forall
t\geq 0.$

Then, we have that for all $\theta<\min\\{\theta_{0},\theta_{1}\\}$,

$\displaystyle\mathbb{E}_{\pi}[f(\mathbf{x}_{\infty},\theta)]=\frac{1}{\gamma(\theta)}\mathbb{E}_{\pi}[\beta(\mathbf{x}_{\infty},\theta)]<\infty.$

###### Proof of Lemma 11.

For the complete part of the proof, we assume that
$\theta<\min\\{\theta_{0},\theta_{1}\\}$. From, we Condition C1, we have that
$\mathbb{E}[f(\mathbf{x}(t),\theta)]<\infty,\ \forall t\geq 0$. Thus, in Eq.
(17) in Condition C2, we can take expectation on both sides to get that,

$\displaystyle\mathbb{E}[f(\mathbf{x}(t+1),\theta)]=(1-\gamma(\theta))\mathbb{E}[f(\mathbf{x}(t),\theta)]+\mathbb{E}[\beta(\mathbf{x}(t),\theta)].$

Depending on the value of $\gamma(\theta)$, we consider two cases.

  * •

Case 1: $\gamma(\theta)\geq 1$.  In this case, we have that
$1-\gamma(\theta)<0$. Then, as $f(\mathbf{x}(t),\theta)\geq 0$, from the
previous equation, we have

$\displaystyle\mathbb{E}[f(\mathbf{x}(t+1),\theta)]\leq\mathbb{E}[\beta(\mathbf{x}(t),\theta)]\leq\bar{\beta}(\theta).$

As $\mathbb{E}[f(\mathbf{x}(t),\theta)]$ is uniformly bounded by $\bar{\beta}$
for all $t\geq 0$,

$\displaystyle\mathbb{E}[f(\mathbf{x}_{\infty},\theta)]=\lim_{t\rightarrow\infty}\mathbb{E}[f(\mathbf{x}(t+1),\theta)]\leq\bar{\beta}(\theta).$

Thus, under steady state, we can rewrite Eq. (17) in Condition C2 as,

$\displaystyle\mathbb{E}[f(\mathbf{x}_{\infty},\theta)]=(1-\gamma(\theta))\mathbb{E}[f(\mathbf{x}_{\infty},\theta)]+\mathbb{E}[\beta(\mathbf{x}_{\infty},\theta)].$

This gives us the result. Note that the above result also holds when we
directly have that $\mathbb{E}[f(\mathbf{x}(t),\theta)]$ is uniformly bounded,
as mentioned in Condition C2 part (2).

  * •

Case 2: $\gamma(\theta)\in(0,1)$.  In this case, we can use Eq. (17) in
Condition C2 recursively, to get that

$\displaystyle\mathbb{E}[f(\mathbf{x}(t+1),\theta)]$
$\displaystyle=(1-\gamma(\theta))^{t+1}\mathbb{E}[f(\mathbf{x}(0),\theta)]$
$\displaystyle\quad+\sum_{i=0}^{t}(1-\gamma(\theta))^{i}\mathbb{E}[\beta(\mathbf{x}(t-i),\theta)].$
(18)

Now, as $\lim_{t\rightarrow\infty}\mathbb{E}[\beta(\mathbf{x}(t),\theta)]$
exists, and is equal to $\mathbb{E}[\beta(\mathbf{x}_{\infty},\theta)]$, we
have that

$\displaystyle\lim_{t\rightarrow\infty}\sum_{i=0}^{t}(1-\gamma(\theta))^{i}\mathbb{E}[\beta(\mathbf{x}(t-i),\theta)]$
$\displaystyle=\frac{1}{\gamma(\theta)}\mathbb{E}[\beta(\mathbf{x}_{\infty},\theta)],$
$\displaystyle\lim_{t\rightarrow\infty}(1-\gamma(\theta))^{t+1}\mathbb{E}[f(\mathbf{x}(0),\theta)]$
$\displaystyle=0.$

Substituting the above in Eq. (• ‣ A) gives us the result.

∎

##  Appendix B Proof of results in Section 2

###### Proof of Proposition 3.

For any $i\in\\{1,2,\dots,n\\}$, we consider the following exponential
Lyapunov functions

$\displaystyle V_{\perp i}(\mathbf{x};\theta)=e^{\theta\mathbf{x}_{\perp
i}}=e^{\theta(x_{i}-\hat{x})},$ $\displaystyle U_{\perp
i}(\mathbf{x};\theta)=e^{\theta\mathbf{x}_{\perp
i}}=e^{-\theta(x_{i}-\hat{x})},$

where $\hat{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$.

First, note that the total queue length at any time $t$ is less than the total
number of arrivals till time $t$. And as the arrivals follow a Poisson
process, the total number of arrivals till time $t$ is a Poisson random
variable with parameter $\lambda_{n}t$. This gives us that

$\displaystyle\mathbb{E}[V_{\perp i}({\bf
q}(t);\theta)]\leq\exp\Big{(}\lambda_{n}t\big{(}e^{\theta}-1\big{)}\Big{)}<\infty.$

Thus, $V_{\perp i}(\cdot)$ satisfies the Condition C1 of Lemma 11. Similarly,
$U_{\perp i}(\cdot)$ also satisfies the Condition C1 of Lemma 11.

Suppose $G$ denotes the generator matrix of the underlying CTMC of the JSQ
system. Then, the drift of $V_{\perp i}({\bf q}(t);\theta)$, denoted by
$\Delta V_{\perp i}({\bf q}(t);\theta)$, is given by

$\displaystyle\mathbb{E}\big{[}\Delta V_{\perp i}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}:=\sum_{{\bf q}^{\prime}}G({\bf
q}^{\prime},{\bf q}(t))\big{[}V_{\perp i}({\bf q}^{\prime};\theta)-V_{\perp
i}({\bf q}(t);\theta)\big{]}.$

Using this, we have

$\displaystyle\mathbb{E}\big{[}$ $\displaystyle\Delta V_{\perp i}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}$
$\displaystyle=\lambda_{n}\left(e^{\theta q_{\perp
i}(t)+\theta\big{(}1-\frac{1}{n}\big{)}}\mathds{1}_{\\{i=i^{*}\\}}+e^{\theta
q_{\perp i}(t)-\frac{\theta}{n}}\mathds{1}_{\\{i\neq i^{*}\\}}-e^{\theta
q_{\perp i}(t)}\right)\allowdisplaybreaks$ $\displaystyle\ \ \ \
+\sum_{j=1}^{n}\mu\left(e^{\theta q_{\perp
i}(t)-\theta\big{(}1-\frac{1}{n}\big{)}}\mathds{1}_{\\{j=i\\}}+e^{\theta
q_{\perp i}(t)+\frac{\theta}{n}}\mathds{1}_{\\{j\neq i\\}}-e^{\theta q_{\perp
i}(t)}\right)\mathds{1}_{\\{q_{j}>0\\}}\allowdisplaybreaks$
$\displaystyle=\lambda_{n}e^{\theta q_{\perp
i}(t)}\left(e^{\theta\big{(}1-\frac{1}{n}\big{)}}\mathds{1}_{\\{i=i^{*}\\}}+e^{-\frac{\theta}{n}}\mathds{1}_{\\{i\neq
i^{*}\\}}-1\right)\allowdisplaybreaks$ $\displaystyle\ \ \ \ +\mu e^{\theta
q_{\perp
i}(t)}\sum_{j=1}^{n}\left(e^{-\theta\big{(}1-\frac{1}{n}\big{)}}\mathds{1}_{\\{j=i\\}}+e^{\frac{\theta}{n}}\mathds{1}_{\\{j\neq
i\\}}-1\right)\mathds{1}_{\\{q_{j}>0\\}}\allowdisplaybreaks$
$\displaystyle=\lambda_{n}e^{\theta q_{\perp
i}(t)}\left(\big{(}e^{\theta}-1\big{)}e^{-\frac{\theta}{n}}\mathds{1}_{\\{i=i^{*}\\}}+e^{-\frac{\theta}{n}}-1\right)$
$\displaystyle\quad+\mu e^{\theta q_{\perp
i}(t)}\Big{(}\big{(}e^{-\theta}-1\big{)}e^{\frac{\theta}{n}}\mathds{1}_{\\{q_{i}(t)>0\\}}+\big{(}e^{\frac{\theta}{n}}-1\big{)}\sum_{j}\mathds{1}_{\\{q_{j}>0\\}}\Big{)},$

where $i^{*}$ is the index chosen by the JSQ dispatcher.

Next, we consider the function
$V_{\perp}(\mathbf{x};\theta)=\sum_{i=1}^{n}V_{\perp i}(\mathbf{x};\theta)$.
Then, the drift of $V_{\perp}({\bf q}(t);\theta)$ is just the sum of the
drifts of $V_{\perp i}({\bf q}(t);\theta)$, that is,

$\displaystyle\mathbb{E}\big{[}$ $\displaystyle\Delta V_{\perp}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}\allowdisplaybreaks$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\big{[}\Delta V_{\perp i}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}\allowdisplaybreaks$
$\displaystyle=\lambda_{n}\left(e^{\theta}-1\right)e^{-\frac{\theta}{n}}e^{\theta
q_{\perp,\min}(t)}+\lambda_{n}\left(e^{-\frac{\theta}{n}}-1\right)\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \
+\mu\sum_{i=1}^{n}e^{\theta q_{\perp
i}(t)}\Big{(}\big{(}e^{-\theta}-1\big{)}e^{\frac{\theta}{n}}\mathds{1}_{\\{q_{i}(t)>0\\}}+\big{(}e^{\frac{\theta}{n}}-1\big{)}\sum_{j}\mathds{1}_{\\{q_{j}>0\\}}\Big{)}\allowdisplaybreaks$
$\displaystyle=\lambda_{n}\left(e^{\theta}-1\right)e^{-\frac{\theta}{n}}e^{\theta
q_{\perp,\min}(t)}+\lambda_{n}\left(e^{-\frac{\theta}{n}}-1\right)\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \
+\mu\left(e^{-\theta}-1\right)e^{\frac{\theta}{n}}\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}+n\mu\left(e^{\frac{\theta}{n}}-1\right)\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \
-\mu\left(e^{-\theta}-1\right)e^{\frac{\theta}{n}}\sum_{i=1}^{n}e^{\theta
q_{\perp
i}(t)}\mathds{1}_{\\{q_{i}(t)=0\\}}-\mu\left(e^{\frac{\theta}{n}}-1\right)\sum_{j=1}^{n}\mathds{1}_{\\{q_{j}(t)=0\\}}\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}.$

Next, by using $\theta>0$,
$q_{\perp,\min}(t)=\min_{i}\big{\\{}q_{i}(t)-\frac{1}{n}\sum_{j=1}^{n}q_{j}(t)\big{\\}}<0$,
and $e^{\theta q_{\perp
i}(t)}\mathds{1}_{\\{q_{i}(t)=0\\}}=e^{-\frac{\theta}{n}\sum_{j=1}^{n}q_{j}(t)}\mathds{1}_{\\{q_{i}(t)=0\\}}\leq
1$, we have

$\displaystyle\mathbb{E}\big{[}\Delta V_{\perp}({\bf q}(t);\theta)\big{|}{\bf
q}(t)\big{]}$
$\displaystyle\leq\lambda_{n}\left(e^{\theta}-1\right)e^{-\frac{\theta}{n}}+n\mu\left(1-e^{-\theta}\right)e^{\frac{\theta}{n}}$
$\displaystyle\quad+\Big{(}\lambda_{n}\left(e^{-\frac{\theta}{n}}-1\right)+\mu\left(e^{-\theta}-1\right)$
$\displaystyle\quad+n\mu\left(e^{\frac{\theta}{n}}-1\right)\Big{)}\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}$

For any $\theta>0$, we have that

$\displaystyle-\frac{1}{n}e^{\theta}\leq\frac{e^{-\frac{\theta}{n}}-1}{1-e^{-\theta}}\leq-\frac{1}{n},$
$\displaystyle\frac{1}{n}\leq\frac{e^{\frac{\theta}{n}}-1}{1-e^{-\theta}}\leq\frac{1}{n}e^{\theta},$

This gives us that, for any $\theta>0$,

$\displaystyle\frac{1}{1-e^{-\theta}}\mathbb{E}\big{[}\Delta V_{\perp}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}$
$\displaystyle\leq\lambda_{n}e^{\theta}+n\mu
e^{\frac{\theta}{n}}+\left(-\frac{\lambda_{n}}{n}+\mu(e^{\theta}-1)\right)\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\leq 2n\mu
e^{\theta}-\frac{\lambda_{n}}{n}\sum_{i=1}^{n}e^{\theta q_{\perp
i}(t)}+\mu\theta e^{\theta}\sum_{i=1}^{n}e^{\theta q_{\perp i}(t)}$ (19)

By using similar calculation as in Eq. (19), and defining
$U_{\perp}(\mathbf{x};\theta):=\sum_{i=1}^{n}U_{\perp i}(\mathbf{x};\theta)$,
we have that

$\displaystyle\mathbb{E}\big{[}$ $\displaystyle\Delta U_{\perp}({\bf
q}(t);\theta)\big{|}{\bf q}(t)\big{]}\allowdisplaybreaks$
$\displaystyle=\lambda_{n}\left(e^{-\theta}-1\right)e^{\frac{\theta}{n}}e^{-\theta
q_{\perp,\min}(t)}+\lambda_{n}\left(e^{\frac{\theta}{n}}-1\right)\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \
+\mu\left(e^{\theta}-1\right)e^{-\frac{\theta}{n}}\sum_{i=1}^{n}e^{-\theta
q_{\perp
i}(t)}+n\mu\left(e^{-\frac{\theta}{n}}-1\right)\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \
-\mu\left(e^{\theta}-1\right)e^{-\frac{\theta}{n}}\sum_{j=1}^{n}e^{-\theta
q_{\perp i}(t)}\mathds{1}_{\\{q_{i}(t)=0\\}}\allowdisplaybreaks$
$\displaystyle\quad-\mu\left(e^{-\frac{\theta}{n}}-1\right)\sum_{j=1}^{n}\mathds{1}_{\\{q_{j}(t)=0\\}}\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}.$

Thus, for any $\theta>0$, by using $q_{\perp
i}(t)=q_{\perp,\min}(t)=-\frac{1}{n}\sum_{j=1}^{n}q_{i}(t)$ whenever
$q_{i}(t)=0$, we have

$\displaystyle\frac{1}{1-e^{-\theta}}\mathbb{E}\big{[}$ $\displaystyle\Delta
U_{\perp}({\bf q}(t);\theta)\big{|}{\bf q}(t)\big{]}\allowdisplaybreaks$
$\displaystyle\leq-\lambda_{n}e^{-\theta
q_{\perp,\min}(t)}+\left(\frac{\lambda_{n}}{n}e^{\theta}+\mu(e^{\theta}-1)\right)\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}\allowdisplaybreaks$ $\displaystyle\ \ \ \ -\mu e^{-\theta
q_{\perp,\min}(t)}\sum_{j=1}^{n}\mathds{1}_{\\{q_{j}(t)=0\\}}+\frac{\mu}{n}e^{\theta}\sum_{j=1}^{n}\mathds{1}_{\\{q_{j}(t)=0\\}}\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}-\lambda_{n}e^{-\theta
q_{\perp,\min}(t)}+\left(\frac{\lambda_{n}}{n}e^{\theta}+2\mu(e^{\theta}-1)\right)\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}-\lambda_{n}e^{-\theta
q_{\perp,\min}(t)}+\frac{\lambda_{n}}{n}\sum_{i=1}^{n}e^{-\theta q_{\perp
i}(t)}+3\mu\theta e^{\theta}\sum_{i=1}^{n}e^{-\theta q_{\perp i}(t)},$ (20)

where $(a)$ follows by using $e^{-\theta
q_{\perp,\min}(t)}\geq\sum_{i=1}^{n}e^{-\theta q_{\perp i}(t)}$ and
$\sum_{j=1}^{n}\mathds{1}_{\\{q_{j}(t)=0\\}}\leq n$; and (b) follows by using
$\frac{\lambda_{n}}{n}e^{\theta}=\frac{\lambda_{n}}{n}+\frac{\lambda_{n}}{n}(e^{\theta}-1)\leq\frac{\lambda_{n}}{n}+\mu\theta
e^{\theta}$.

Before proceeding with the next step, we provide a bound on the following
function.

$\displaystyle D(\mathbf{x};\theta):=-e^{-\theta
x_{\min}}+\frac{1}{n}\sum_{i=1}^{n}e^{-\theta
x_{i}}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}},$

where $\mathbf{x}\in\mathbb{R}^{n}$ such that $\sum_{i=1}^{n}x_{i}=0$ and
$x_{\min}=\min_{i}x_{i}$. Then, $x_{\min}\leq 0$. Suppose
$\mathcal{I}^{+}=\\{i:x_{i}\geq 0\\}$ and
$\mathcal{I}^{-}=[n]\backslash\mathcal{I}^{+}$. Suppose
$n_{+}=|\mathcal{I}^{+}|$ and $n_{-}=|\mathcal{I}^{-}|$

  * 1.

Case 1: Consider the case

$n_{+}\geq n/4,\text{ and }n_{-}\leq 3n/4.$

In this case, we have

$\displaystyle D(\mathbf{x};\theta)$ $\displaystyle=-e^{-\theta
x_{\min}}+\frac{1}{n}\sum_{i\in\mathcal{I}^{+}}e^{-\theta
x_{i}}+\frac{1}{n}\sum_{i\in\mathcal{I}^{-}}e^{-\theta
x_{i}}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}-e^{-\theta
x_{\min}}+\frac{n_{+}}{n}+\frac{n_{-}}{n}e^{-\theta
x_{\min}}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}}$
$\displaystyle=\frac{n_{+}}{n}\big{(}1-e^{-\theta
x_{\min}}\big{)}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\frac{1}{4}\big{(}1-e^{-\theta
x_{\min}}\big{)}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}}$
$\displaystyle\leq\frac{1}{4}-\frac{1}{4n}\sum_{i=1}^{n}e^{-\theta
x_{i}}-\frac{1}{n}\sum_{i=1}^{n}e^{\theta x_{i}},$

where (a) follows by using $x_{i}\geq 0$ for all $i\in\mathcal{I}^{+}$; and
(b) follows by using $x_{\min}\leq 0$.

  * 2.

Case 2: Consider the case

$\frac{1}{n}\sum_{i=1}^{n}e^{-\theta x_{i}}\leq\frac{1}{2}e^{-\theta
x_{\min}}+\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}.$

In this case, we easily have that,

$\displaystyle D(\mathbf{x};\theta)$ $\displaystyle\leq-\frac{1}{2}e^{-\theta
x_{\min}}-\frac{1}{2n}\sum_{i=1}^{n}e^{\theta
x_{i}}\leq-\frac{1}{2n}\sum_{i=1}^{n}e^{-\theta
x_{i}}-\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}$

  * 3.

Case 3: Consider the case

$\frac{1}{n}\sum_{i=1}^{n}e^{-\theta x_{i}}>\frac{1}{2}e^{-\theta
x_{\min}}+\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}\text{ and }n_{+}\leq
n/4.$

In this case, suppose $\sum_{i\in\mathcal{I}^{+}}x_{i}=B$, then
$\sum_{i\in\mathcal{I}^{-}}-x_{i}=B$ as $\sum_{i=1}^{n}x_{i}=0$. Further, as
$n_{-}\geq 3n/4$, suppose $\mathcal{I}^{-}_{1}$ is the collection of indices
of $\frac{n}{2}$ smallest elements in $\mathcal{I^{-}}$ and
$\mathcal{I}^{-}_{2}=\mathcal{I^{-}}\backslash\mathcal{I}^{-}_{1}$. Then, for
any $i\in\mathcal{I}^{-}_{1}$ and $j\in\mathcal{I}^{-}_{2}$, we have
$x_{i}\leq x_{j}$.

We can argue that for any $j\in\mathcal{I}^{-}_{2}$, we have
$-x_{j}\leq\frac{2B}{n}$. This can be shown by contraction. Suppose $\exists
j\in I^{-}_{2}$ such that $-x_{j}>\frac{2B}{n}$. Then, $-x_{i}>\frac{2B}{n}$
for all $i\in I^{-}_{1}$. As $|I^{-}_{1}|=\frac{n}{2}$, we have $\sum_{i\in
I^{-}}-x_{i}\geq\sum_{i\in I^{-}_{1}}-x_{i}>B$. This creates a contradiction
as $\sum_{i\in\mathcal{I}^{-}}-x_{i}=B$. Thus, $-x_{j}\leq\frac{2B}{n}$ for
all $j\in\mathcal{I}^{-}_{2}$.

This implies that,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}e^{-\theta x_{i}}$
$\displaystyle\leq\frac{1}{n}\sum_{i\in\mathcal{I}^{+}}e^{-\theta
x_{i}}+\frac{1}{n}\sum_{i\in\mathcal{I}^{-}_{2}}e^{-\theta
x_{i}}+\frac{1}{n}\sum_{i\in\mathcal{I}^{-}_{1}}e^{-\theta
x_{i}}\allowdisplaybreaks$
$\displaystyle\leq\frac{n_{+}}{n}+\frac{1}{n}\sum_{i\in\mathcal{I}^{-}_{2}}e^{\frac{2\theta
B}{n}}+\frac{1}{2}e^{-\theta x_{\min}}\allowdisplaybreaks$
$\displaystyle\leq\frac{1}{4}+\frac{1}{2}e^{\frac{2\theta
B}{n}}+\frac{1}{2}e^{-\theta x_{\min}}.$ (21)

Next, we have further two cases

    * –

Case 3.1: $\theta B\leq\frac{n}{2}$. In this case, we have that
$e^{\frac{2\theta B}{n}}\leq e$.

    * –

Case 3.2: $\theta B\geq\frac{n}{2}$. In this case, by using the Jensen’s
inequality,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}e^{\theta
x_{i}}\geq\frac{1}{n}\sum_{i\in\mathcal{I}^{+}}e^{\theta
x_{i}}\geq\frac{n_{+}}{n}e^{\frac{\theta
B}{n_{+}}}\geq\frac{1}{4}e^{\frac{4\theta B}{n}},$

where the last inequality holds as $\theta B\geq\frac{n}{2}$, and so,
$\frac{n_{+}}{n}e^{\frac{\theta B}{n_{+}}}$ is minimized by taking
$n_{+}=\frac{n}{4}$. Thus, by using the assumption
$\frac{1}{n}\sum_{i=1}^{n}e^{-\theta x_{i}}>\frac{1}{2}e^{-\theta
x_{\min}}+\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}$ and Eq. (3.), we have
that

$\displaystyle\frac{1}{8}e^{\frac{4\theta
B}{n}}\leq\frac{1}{4}+\frac{1}{2}e^{\frac{2\theta
B}{n}}\leq\frac{3}{4}e^{\frac{2\theta B}{n}}.$

This gives us that $e^{\frac{2\theta B}{n}}\leq 6$.

Thus, in both the above cases, $e^{\frac{2\theta B}{n}}\leq 6$. This gives us

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}e^{-\theta
x_{i}}\leq\frac{1}{4}+3+\frac{1}{2}e^{-\theta x_{\min}}\leq
4+\frac{1}{2}e^{-\theta x_{\min}}.$

And thus,

$\displaystyle D(\mathbf{x};\theta)$ $\displaystyle\leq
4-\frac{1}{2}e^{-\theta x_{\min}}-\frac{1}{2n}\sum_{i=1}^{n}e^{\theta
x_{i}}\leq 4-\frac{1}{2n}\sum_{i=1}^{n}e^{-\theta
x_{i}}-\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}$

By combining the results from Case 1, Case 2 and Case 3, we get

$\displaystyle D(\mathbf{x};\theta)\leq 4-\frac{1}{4n}\sum_{i=1}^{n}e^{-\theta
x_{i}}-\frac{1}{2n}\sum_{i=1}^{n}e^{\theta x_{i}}$ (22)

Now, from Eq. (19) and (20), for any $\theta>0$, we have

$\displaystyle\frac{1}{1-e^{-\theta}}$ $\displaystyle\big{(}\mathbb{E}[\Delta
V_{\perp}({\bf q}(t);\theta)|{\bf q}(t)]+\mathbb{E}[\Delta U({\bf
q}(t);\theta)|{\bf q}(t)]\big{)}$ $\displaystyle\leq 2n\mu
e^{\theta}+\lambda_{n}D({\bf q}(t);\theta)+\mu\theta
e^{\theta}\sum_{i=1}^{n}e^{\theta q_{\perp i}(t)}+3\mu\theta
e^{\theta}\sum_{i=1}^{n}e^{-\theta q_{\perp i}(t)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}2n\mu
e^{\theta}+4\lambda_{n}+\Big{(}\mu\theta
e^{\theta}-\frac{\lambda_{n}}{2n}\Big{)}\sum_{i=1}^{n}e^{\theta q_{\perp
i}(t)}$ $\displaystyle\quad+\Big{(}3\mu\theta
e^{\theta}-\frac{\lambda_{n}}{4n}\Big{)}\sum_{i=1}^{n}e^{-\theta q_{\perp
i}(t)}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}8n\mu+\Big{(}2\mu\theta-\frac{\lambda_{n}}{2n}\Big{)}\sum_{i=1}^{n}e^{\theta
q_{\perp
i}(t)}+\Big{(}6\mu\theta-\frac{\lambda_{n}}{4n}\Big{)}\sum_{i=1}^{n}e^{-\theta
q_{\perp i}(t)}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}8n\mu-\frac{\lambda_{n}}{8n}\sum_{i=1}^{n}e^{\theta
q_{\perp i}(t)}-\frac{\lambda_{n}}{8n}\sum_{i=1}^{n}e^{-\theta q_{\perp
i}(t)},$

where (a) follows using Eq. (22); (b) follows by taking $\theta<1/2$, and so
$e^{\theta}<2$; and (c) follows by using
$\theta\leq\frac{\lambda_{n}}{48n\mu}$. Thus, the Lyapunov function
$V_{\perp}(\cdot)+U_{\perp}(\cdot)$ satisfies the Condition C2 of Lemma 11.
Thus, by using Lemma 11, the drift is zero in steady state, i.e.,

$\displaystyle\mathbb{E}_{\pi_{n}}[\Delta V_{\perp}({\bf
q};\theta)]+\mathbb{E}_{\pi_{n}}[\Delta U({\bf q};\theta)]=0.$

So, for any $n\geq 2$, by picking
$\theta_{\perp}:=\frac{1}{96}\leq\frac{\lambda_{n}}{48n\mu}$, we have that for
any $0<\theta\leq\theta_{\perp}$,

$\displaystyle\mathbb{E}_{\pi_{n}}\Big{[}\sum_{i=1}^{n}e^{-\theta q_{\perp
i}}+\sum_{i=1}^{n}e^{\theta q_{\perp
i}}\Big{]}\leq\frac{64n^{2}\mu}{\lambda_{n}}=\frac{64n}{1-\epsilon_{n}}\leq
128n,$

where the last inequality holds whenever $\epsilon_{n}\leq 1/2$. Now, as the
servers are homogeneous, by using symmetry of the JSQ system, we have

$\displaystyle\mathbb{E}_{\pi_{n}}\big{[}e^{\theta|q_{\perp
i}|}\big{]}\leq\mathbb{E}_{\pi_{n}}\big{[}e^{\theta q_{\perp
i}}\big{]}+\mathbb{E}_{\pi_{n}}\big{[}e^{-\theta q_{\perp i}}\big{]}\leq 128.$

This completes the proof. ∎

###### Proof of Lemma 4.

Recall that we use the notation $\lambda_{n}=n\mu(1-\epsilon_{n})$ and
$\overline{q}=\sum_{i=1}^{n}q_{i}$. Consider the exponential Lyapunov function

$\displaystyle
V(\mathbf{x};\theta):=\exp\Big{(}\theta\epsilon_{n}\sum_{i=1}^{n}x_{i}\Big{)}.$
(23)

Note that the total queue length at any time $t$ is less than the total number
of arrivals till time $t$. And as the arrivals follow a Poisson process, the
total number of arrivals till time $t$ is a Poisson random variable with
parameter $\lambda_{n}t$. This gives us that

$\displaystyle\mathbb{E}[V({\bf
q}(t);\theta)]\leq\exp\Big{(}\lambda_{n}t\big{(}e^{\theta\epsilon_{n}}-1\big{)}\Big{)}<\infty.$

Thus, $V(\cdot)$ satisfies the Condition C1 of Lemma 11.

Next, the drift of the function $V(\cdot)$ is given by

$\displaystyle\mathbb{E}[\Delta V({\bf q}(t);\theta)|{\bf q}(t)]$
$\displaystyle=\Big{[}\lambda_{n}\Big{(}e^{\epsilon_{n}\theta}-1\Big{)}+\mu\Big{(}e^{-\epsilon_{n}\theta}-1\Big{)}\sum_{i=1}^{n}\mathds{1}_{\\{q_{i}(t)>0\\}}\Big{]}V({\bf
q}(t);\theta)$
$\displaystyle=\Big{(}1-e^{-\epsilon_{n}\theta}\Big{)}\Big{[}\lambda_{n}e^{\epsilon_{n}\theta}-n\mu\Big{]}V({\bf
q}(t);\theta)\allowdisplaybreaks$
$\displaystyle\quad+\mu\Big{(}1-e^{-\epsilon_{n}\theta}\Big{)}\sum_{i=1}^{n}\mathds{1}_{\\{q_{i}(t)=0\\}}V({\bf
q}(t);\theta)$
$\displaystyle=\Big{(}1-e^{-\epsilon_{n}\theta}\Big{)}\Big{(}-\gamma_{n}(\theta)V({\bf
q}(t);\theta)+\beta_{n}({\bf q}(t);\theta)\Big{)},$

where

$\displaystyle\gamma_{n}(\theta):=n\mu-\lambda_{n}e^{\epsilon_{n}\theta},$
$\displaystyle\beta_{n}({\bf
q}(t);\theta):=\mu\sum_{i=1}^{n}\mathds{1}_{\\{q_{i}(t)=0\\}}V({\bf
q}(t);\theta).$

Then, we have that

$\displaystyle\gamma_{n}(\theta)>0\text{ for any
}\theta<\theta_{n}:=\frac{1}{\epsilon_{n}}\log\Big{(}\frac{n\mu}{\lambda_{n}}\Big{)}=-\frac{1}{\epsilon_{n}}\log\big{(}1-\epsilon_{n}\big{)}.$
(24)

Note that as $\theta_{n}\geq 1$, $\gamma_{n}(\theta)>0$ holds for any
$\theta\leq 1$.

Next, by using Eq. (23), we have

$\displaystyle\mathbb{E}[\beta_{n}({\bf q}(t);\theta)]\leq
n\mu\mathbb{E}[V({\bf q}(t);\theta)]<\infty.$ (25)

And, in steady state, for $\theta\in(0,\theta_{n})$

$\displaystyle\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf q};\theta)]$
$\displaystyle=\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}e^{\theta\epsilon_{n}\overline{q}}\Big{]}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}e^{\theta
n\epsilon_{n}(q_{i}-q_{\perp i})}\Big{]}\allowdisplaybreaks$
$\displaystyle=\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}e^{-\theta
n\epsilon_{n}q_{\perp i}}\Big{]}\allowdisplaybreaks$
$\displaystyle\leq\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}e^{\theta
n\epsilon_{n}|q_{\perp i}|}\Big{]}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\kappa_{\perp}n\mu,$ (26)

where (a) follows by using $q_{\perp i}=q_{i}-\frac{1}{n}\overline{q}$ by
definition; and (b) follows by using Proposition 3, whenever $\theta
n\epsilon_{n}<\theta_{\perp}$ for all $\theta\in(0,\theta_{n})$, where
$\theta_{\perp}$ is as given in Proposition 3. Then, for any
$\theta\in(0,\theta_{n})$ and $\epsilon_{n}\leq\frac{\theta_{\perp}}{\theta
n}$, by combining Eq. (24), (25) and (B), we have that $V(\cdot)$ satisfies
the Condition C2 of Lemma 11. Thus, by using Lemma 11, we have that for any
$\theta\in(0,\theta_{n})$

$\displaystyle\mathbb{E}_{\pi_{n}}[V({\bf
q};\theta)]=\frac{1}{\gamma_{n}(\theta)}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf
q};\theta)].$ (27)

Further, when $\theta<0$, we easily have

$\displaystyle\mathbb{E}_{\pi_{n}}[V({\bf
q};\theta)]=\mathbb{E}_{\pi_{n}}\Big{[}e^{\theta\epsilon_{n}\overline{q}}\Big{]}\leq
1,$ $\displaystyle\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf q};\theta)]$
$\displaystyle\leq\mu\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\Big{]}\leq
n\mu.$

Thus, as $\mathbb{E}_{\pi_{n}}[V({\bf q};\theta)]$ and
$\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf q};\theta)]$ are bounded, by Condition C2
of Lemma 11, we have that for any $\theta<0$, the Eq. (27) is still satisfied.
This completes the proof. ∎

###### Proof of Claim 1.

By using Lemma 4 and taking $\theta<0$ and $\theta\rightarrow 0$, we have

$\displaystyle\mu\mathbb{E}_{\pi_{n}}\big{[}\sum_{i=1}^{n}\mathds{1}_{\\{q_{i}=0\\}}\big{]}$
$\displaystyle=\lim_{\theta\rightarrow
0^{-}}\mathbb{E}_{\pi_{n}}[\beta_{n}({\bf q};\theta)]\allowdisplaybreaks$
$\displaystyle=\lim_{\theta\rightarrow
0^{-}}\gamma_{n}(\theta)\mathbb{E}_{\pi_{n}}[V({\bf
q};\theta)]\allowdisplaybreaks$
$\displaystyle=n\mu-\lambda_{n}\allowdisplaybreaks$
$\displaystyle=n\epsilon_{n}\mu.$

Also, by symmetry, we have
$\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{]}=\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{j}=0\\}}\big{]}$
for all $i,j$. Thus, we get,

$\displaystyle\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{]}=\epsilon_{n},\quad\forall
i\in\\{1,2,\dots,n\\}.$

Further, we have

$\displaystyle\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{\theta\epsilon_{n}\overline{q}}-1\big{|}\Big{]}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{\theta
n\epsilon_{n}(q_{i}-q_{\perp i})}-1\big{|}\Big{]}\allowdisplaybreaks$
$\displaystyle=\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{-\theta
n\epsilon_{n}q_{\perp i}}-1\big{|}\Big{]}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}n\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{-\theta
n\epsilon_{n}q_{\perp i}}-1\big{|}\Big{]}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}n\mathbb{E}_{\pi_{n}}\big{[}\mathds{1}_{\\{q_{i}=0\\}}^{p}\big{]}^{\frac{1}{p}}\mathbb{E}_{\pi_{n}}\Big{[}\big{|}e^{-\theta
n\epsilon_{n}q_{\perp
i}}-1\big{|}^{r}\Big{]}^{\frac{1}{r}}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}\theta
n^{2}\epsilon_{n}^{1+\frac{1}{p}}\mathbb{E}_{\pi_{n}}\Big{[}|q_{\perp
i}|^{r}e^{r|\theta|n\epsilon_{n}|q_{\perp
i}|}\Big{]}^{\frac{1}{r}}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{\leq}}\theta
n^{2}\epsilon_{n}^{1+\frac{1}{p}}\mathbb{E}_{\pi_{n}}\big{[}|q_{\perp
i}|^{2r}\big{]}^{\frac{1}{2r}}\mathbb{E}_{\pi_{n}}\Big{[}e^{2r|\theta|n\epsilon_{n}|q_{\perp
i}|}\Big{]}^{\frac{1}{2r}}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(f)}}{{\leq}}\theta
n^{2}\epsilon_{n}^{1+\frac{1}{p}}\times\frac{2\kappa_{\perp}r}{\theta_{\perp}}\times\kappa_{\perp},$
(28)

where (a) follows by using $q_{\perp i}=q_{i}-\overline{q}$; (b) follows by
using the symmetry of the system as the servers are homogeneous; (c) and (e)
follows by using Hölder’s inequality by choosing $r,p>0$ and
$\frac{1}{r}+\frac{1}{p}=1$; (d) follows as $|e^{x}-1|\leq|x|e^{|x|}$ for all
$x\in\mathbb{R}$ and Eq. (6); and (f) follows by using Proposition 3, by
choosing $n$ large enough such that $2r|\theta|n\epsilon_{n}<\theta_{\perp}$.

Now, by choosing $r=-\log\epsilon_{n}$ and
$p=\frac{\log\epsilon_{n}}{\log\epsilon_{n}+1}$, we have

$\displaystyle\sum_{i=1}^{n}\mathbb{E}_{\pi_{n}}\Big{[}\mathds{1}_{\\{q_{i}=0\\}}\big{|}e^{\theta\epsilon_{n}\overline{q}}-1\big{|}\Big{]}$
$\displaystyle\leq\frac{2e\kappa_{\perp}^{2}}{\theta_{\perp}}\alpha\theta
n^{2}\epsilon_{n}^{2}\log\Big{(}\frac{1}{\epsilon_{n}}\Big{)},$

where the inequality holds whenever
$2|\theta|n\epsilon_{n}\log\big{(}\frac{1}{\epsilon_{n}}\big{)}<\theta_{\perp}$.
∎

##  Appendix C Proof of results in Section 3

###### Proof of Lemma 6.

First, we show that the MGF of $q(t)$ exists for any $t\geq 0$ and for any
$\theta<\theta_{0}$. Note that as $a(t)\leq A$ almost surely, we have,
$q(t)\leq q(0)+At,\forall t\geq 0$. Thus, by assumption
$V(q(0);\theta)<\infty$, we have

$\displaystyle\mathbb{E}\big{[}e^{\theta\epsilon
q(t)}\big{]}\leq\mathbb{E}[V(q(0);\theta)]e^{\theta\epsilon At}<\infty,\
\forall t\geq 0.$

This shows that $V(\cdot)$ satisfies the Condition C1 of Lemma 11. Further,
the MGF of $a(t),s(t)$ and $u(t)$ also exists as they are bounded random
variables, as given in Section 3.1. Next, the recursion in Eq. (8) gives us
that, $q(t+1)-u(t)=q(t)+a(t)-s(t)$. We write this in an exponential form and
take expectation to get

$\displaystyle\mathbb{E}\big{[}e^{\theta\epsilon(q(t+1)-u(t))}\big{|}q(t)\big{]}$
$\displaystyle=\mathbb{E}\big{[}e^{\theta\epsilon(q(t)+a(t)-s(t))}\big{|}q(t)\big{]}$
$\displaystyle=\mathbb{E}\big{[}e^{\theta\epsilon(a-s)}\big{]}V(q(t);\theta)$
$\displaystyle=(1-\gamma_{\epsilon}(\theta))V(q(t);\theta),$ (29)

where the second equality follows as the arrival and the potential service are
identically distributed across time, and are independent of the queue length
process. Using the relation $q(t+1)u(t)=0,\forall t\geq 0$, we have

$\displaystyle\mathbb{E}\big{[}e^{\theta\epsilon(q(t+1)-u(t))}\big{|}q(t)\big{]}$
$\displaystyle=\mathbb{E}[V(q(t+1);\theta)|q(t)]+\mathbb{E}\big{[}e^{-\theta\epsilon
u(t)}\big{|}q(t)\big{]}-1$
$\displaystyle=\mathbb{E}[V(q(t+1);\theta)|q(t)]-\beta_{\epsilon}(t;\theta).$

Substituting this in the Eq. (C), we arrive at

$\displaystyle\mathbb{E}[V(q(t+1);\theta)|q(t)]=(1-\gamma_{\epsilon}(\theta))V(q(t);\theta)+\beta_{\epsilon}(t;\theta).$
(30)

By taking expectation on both sides, we get

$\displaystyle\mathbb{E}[V(q(t+1);\theta)]=(1-\gamma_{\epsilon}(\theta))\mathbb{E}[V(q(t);\theta)]+\mathbb{E}[\beta_{\epsilon}(t;\theta)].$
(31)

Now the result in Eq. (11) follows by recursively using the Eq. (31).

Further, we can rewrite Eq. (30) as

$\displaystyle\mathbb{E}[\Delta
V(q(t);\theta)|q(t)]=-\gamma_{\epsilon}(\theta)V(q(t);\theta)+\beta_{\epsilon}(t;\theta).$

As $u(t)\leq A$, for any $\theta\in\mathbb{R}$,
$\beta_{\epsilon}(t;\theta)\leq e^{|\theta|\epsilon A}$. Further,
$\gamma_{\epsilon}(\theta)>0$ for any $\theta<\theta_{\epsilon}$. Thus,
$V(\cdot)$ satisfies the Condition C2 of Lemma 11. Thus, from Lemma 11, we get
that for all $\theta<\min\\{\theta_{0},\theta_{\epsilon}\\}$,

$\displaystyle\mathbb{E}[V(q;\theta)]=\lim_{t\rightarrow\infty}\mathbb{E}[V(q(t);\theta)]=\frac{1}{\gamma_{\epsilon}(\theta)}\lim_{t\rightarrow\infty}\mathbb{E}[\beta_{\epsilon}(t;\theta)]=\frac{1-\mathbb{E}_{\pi_{\epsilon}}[e^{-\theta\epsilon
u}]}{1-\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}}.$

This completes the proof. ∎

###### Proof of Claim 2.

By simply equating the drift of the queue length to zero in steady state, that
is

$\displaystyle\mathbb{E}_{\pi_{\epsilon}}\big{[}\mathbb{E}[q(t+1)-q(t)|q(t)]\big{]}=0,$

we have that

$\displaystyle\mathbb{E}_{\pi_{\epsilon}}[u]=\mathbb{E}[s-a]=\epsilon\mu_{\epsilon}.$

Since $e^{y}\geq 1+y$ for all $y\in\mathbb{R}$, we have

$\displaystyle 1-\mathbb{E}_{\pi_{\epsilon}}[e^{-\theta\epsilon
u}]\leq\theta\epsilon\mathbb{E}_{\pi_{\epsilon}}[u]=\epsilon^{2}\theta\mu_{\epsilon}.$

∎

###### Proof of Claim 3.

We have

$\displaystyle\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}$
$\displaystyle\leq
1+\theta\epsilon\mathbb{E}[a-s]+\frac{1}{2}\epsilon^{2}\theta^{2}\mathbb{E}[(a-s)^{2}]+\frac{1}{6}\epsilon^{3}\theta^{3}\mathbb{E}[(a-s)^{3}]$
(32)
$\displaystyle\quad+\epsilon^{4}\theta^{4}\mathbb{E}[(a-s)^{4}e^{\epsilon\theta(a-s)}]$
$\displaystyle\leq
1-\theta\epsilon^{2}\mu_{\epsilon}+\frac{1}{2}\epsilon^{2}\theta^{2}(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})+\frac{1}{6}\epsilon^{3}\theta^{3}E_{3}+(\epsilon\theta
A)^{4}\mathbb{E}[e^{\epsilon\theta(a-s)}]$ (33)

where the second inequality holds by the properties of the arrival and the
service process and by using $E_{3}:=\max\\{0,\mathbb{E}[(a-s)^{3}]\\}$. This
implies that

$\displaystyle\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}\leq\frac{1}{1-(\epsilon\theta
A)^{4}}\Big{(}1-\theta\epsilon^{2}\mu_{\epsilon}+\frac{1}{2}\epsilon^{2}\theta^{2}(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})+\frac{1}{6}\epsilon^{3}\theta^{3}E_{3}\Big{)}.$

Thus, to ensure $\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}<1$, we only
need

$\displaystyle
1-\theta\epsilon^{2}\mu_{\epsilon}+\frac{1}{2}\epsilon^{2}\theta^{2}(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})+\frac{1}{6}\epsilon^{3}\theta^{3}E_{3}\leq
1-(\epsilon\theta A)^{4}.$

This translates to the condition,

$\displaystyle-\mu_{\epsilon}+\frac{1}{2}\theta(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})+\frac{1}{6}\epsilon\theta^{2}E_{3}\leq-\epsilon^{2}\theta^{3}A^{4}.$

By simple calculations, we have that the above inequality holds true whenever

$\displaystyle\theta\in\Bigg{(}0,\mu_{\epsilon}\Big{(}\frac{1}{2}\sigma_{\epsilon}^{2}+\epsilon\frac{\mu_{\epsilon}E_{3}}{3\sigma_{\epsilon}^{2}}+\epsilon^{2}\mu_{\epsilon}^{2}+\epsilon^{2}\frac{4\mu_{\epsilon}^{2}A^{4}}{\sigma_{\epsilon}^{4}}\Big{)}^{-1}\Bigg{)}.$

Thus, by choosing

$\displaystyle\kappa_{1}:=\frac{2\mu_{\epsilon}E_{3}}{3\sigma_{\epsilon}^{4}}+\epsilon\frac{9\mu_{\epsilon}^{2}A^{4}}{\sigma_{\epsilon}^{6}},$
$\displaystyle\tilde{\theta}_{\epsilon}:=\frac{2\mu_{\epsilon}}{\sigma_{\epsilon}^{2}(1+\kappa_{1}\epsilon)},$

we get that

$\displaystyle\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}<1,\quad\forall\theta\in\big{(}0,\tilde{\theta}_{\epsilon}\big{)}.$
(34)

Further, from Eq. (32), for any
$\theta\in\big{(}0,\tilde{\theta}_{\epsilon}\big{)}$, we also have

$\displaystyle 1-\mathbb{E}\big{[}e^{\epsilon\theta(a-s)}\big{]}$
$\displaystyle\geq\theta\epsilon^{2}\mu_{\epsilon}-\frac{1}{2}\epsilon^{2}\theta^{2}(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})-\frac{1}{6}\epsilon^{3}\theta^{3}E_{3}-(\epsilon\theta
A)^{4}\mathbb{E}[e^{\epsilon\theta(a-s)}]\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}\theta\epsilon^{2}\mu_{\epsilon}-\frac{1}{2}\epsilon^{2}\theta^{2}(\sigma_{\epsilon}^{2}+\epsilon^{2}\mu_{\epsilon}^{2})-\frac{1}{6}\epsilon^{3}\theta^{3}E_{3}-(\epsilon\theta
A)^{4}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\geq}}\theta\epsilon^{2}\mu_{\epsilon}\Big{(}1-\theta\frac{\sigma_{\epsilon}^{2}}{2\mu_{\epsilon}}\big{(}1+\kappa_{1}\epsilon\big{)}\Big{)},$

where (a) follows by using Eq. (34); and (b) follows by using
$\theta<\tilde{\theta}_{\epsilon}\leq\frac{2\mu_{\epsilon}}{\sigma_{\epsilon}^{2}}$.
∎

##  Appendix D Proof of results in Section 4

###  D.1 Proof of Lemma 10

###### Proof.

Consider the function $V_{1}(x;\theta)=e^{\theta[x-n]^{+}}$. Then, from the
definition of $q_{n}(t)$ and $w_{n}(t)$, $V_{1}(q_{n}(t);\theta)=e^{\theta
w_{n}(t)}$, and similarly, $V_{1}(q_{n};\theta)=e^{\theta w_{n}}$. Suppose
$a_{n}(t)$ is the total number of arrivals till time $t$. Then, $a_{n}(t)$
follows Poission distribution with parameter $\lambda_{n}t$, and $w_{n}(t)\leq
q_{n}(t)\leq a_{n}(t)$. This gives us that,

$\displaystyle\mathbb{E}[V_{1}(q_{n}(t);\theta)]=\mathbb{E}[e^{\theta
w_{n}(t)}]\leq\mathbb{E}[e^{\theta
a_{n}(t)}]=e^{\lambda_{n}t(e^{\theta}-1)}<\infty.$

Thus, $V_{1}(\cdot)$ satisfies Condition C1 of Lemma 11.

Suppose $G$ denotes the generator matrix of the underlying CTMC of $M/M/n+M$
queue. When the number of jobs in the system is $q$, the drift of
$V_{1}(q;\theta)$, denoted by $\Delta V_{1}(q;\theta)$ is defined as

$\displaystyle\mathbb{E}[\Delta
V_{1}(q(t);\theta)|q(t)=q]:=\sum_{q^{\prime}\in\mathbb{W}}G(q^{\prime},q)\left(V_{1}(q^{\prime};\theta)-V_{1}(q;\theta)\right).$

Then,

$\displaystyle\mathbb{E}[$ $\displaystyle\Delta V_{1}(q(t);\theta)|q(t)=q]$
$\displaystyle=\lambda_{n}V_{1}(q+1;\theta)+\mu\min\\{n,q\\}V_{1}(q-1;\theta)-(\lambda_{n}+\mu\min\\{n,q\\})V_{1}(q;\theta)$
$\displaystyle=e^{\theta[q-n]^{+}}\left(\lambda_{n}\left(e^{\theta}-1\right)\mathds{1}_{\\{q\geq
n\\}}+n\mu\left(e^{-\theta}-1\right)\mathds{1}_{\\{q>n\\}}\right)$
$\displaystyle=\left(e^{\theta}-1\right)\left[e^{\theta[q-n]^{+}}\left(\lambda_{n}-n\mu
e^{-\theta}\right)\mathds{1}_{\\{q>n\\}}+\lambda_{n}\mathds{1}_{\\{q=n\\}}\right]$
$\displaystyle=\left(e^{\theta}-1\right)\left(\lambda_{n}-n\mu
e^{-\theta}\right)e^{\theta[q-n]^{+}}$ $\displaystyle\ \ \ \
-\left(e^{\theta}-1\right)\left(\lambda_{n}\mathds{1}_{\\{q<n\\}}-n\mu
e^{-\theta}\mathds{1}_{\\{q\leq n\\}}\right).$

We define

$\displaystyle\gamma_{n}(\theta)$ $\displaystyle:=-\left(\lambda_{n}-n\mu
e^{-\theta}\right),$ $\displaystyle\beta_{n}(q;\theta)$
$\displaystyle:=\left(n\mu e^{-\theta}\mathds{1}_{\\{q\leq
n\\}}-\lambda_{n}\mathds{1}_{\\{q<n\\}}\right).$

Then, we have

$\displaystyle\mathbb{E}[\Delta
V_{1}(q_{n}(t);\theta)|q_{n}(t)]=\left(e^{\theta}-1\right)\Big{(}-\gamma_{n}(\theta)V_{1}(q_{n}(t);\theta)+\beta_{n}(q_{n}(t);\theta)\Big{)}.$

Note that for any

$\displaystyle\theta<\theta_{n}:=-\log\Big{(}\frac{\lambda_{n}}{n\mu}\Big{)}=\log\frac{1}{1-\epsilon_{n}},$

we have $\gamma_{n}(\theta)>0$. Further, $\beta_{n}(\cdot)$ is a bounded
function as, $\beta_{n}(q_{n}(t);\theta)\leq n\mu\left(e^{\theta}-1\right)$,
and then, as $q_{n}(t)$ converges to $q_{n}$ in distribution, we have

$\displaystyle\mathbb{E}[\beta_{n}(q_{n};\theta)]=\left(e^{\theta}-1\right)\left(n\mu
e^{-\theta}\mathbb{P}(q_{n}\leq n)-\lambda_{n}\mathbb{P}(q_{n}<n)\right)$

Thus, $V_{1}(\cdot)$ also satisfies Condition C2 of Lemma 11. Then, we get
that for all $\theta<\theta_{n}$, we have

$\displaystyle\mathbb{E}[e^{\theta w_{n}}]=\frac{n\mu
e^{-\theta}\mathbb{P}(q_{n}\leq n)-\lambda_{n}\mathbb{P}(q_{n}<n)}{n\mu
e^{-\theta}-\lambda_{n}}.$ (35)

By putting $\theta=0$ in the above equation, we get

$\displaystyle n\mu\mathbb{P}(q_{n}\leq
n)-\lambda_{n}\mathbb{P}(q_{n}<n)=n\mu-\lambda_{n}=n\mu\epsilon_{n}.$

This implies that,

$\displaystyle-n\mu\epsilon_{n}\mathbb{P}(w_{n}>0)+\lambda_{n}\mathbb{P}(q_{n}=n)=0.$
(36)

Further, $\mathbb{E}[e^{\theta w_{n}}]=\mathbb{E}[e^{\theta
w_{n}}\mathds{1}_{\\{w_{n}>0\\}}]+\mathbb{P}(q_{n}\leq n)$. Combining this
with Eq. (35) and (36), we arrive at

$\displaystyle\mathbb{E}\left[e^{\theta
w_{n}}\mathds{1}_{\\{w_{n}>0\\}}\right]=\frac{n\mu\epsilon_{n}\mathbb{P}(w_{n}>0)}{n\mu
e^{-\theta}-\lambda_{n}}.$

This implies that,

$\displaystyle\mathbb{E}\left[e^{\theta
w_{n}}\Big{|}w_{n}>0\right]=\frac{n\mu\epsilon_{n}}{n\mu
e^{-\theta}-\lambda_{n}}=\frac{1}{1-\frac{1}{\epsilon_{n}}\left(1-e^{-\theta}\right)}.$

Next, we evaluate the MGF of $r_{n}$. Consider the function
$V_{2}(x;\theta)=e^{\theta[n-x]^{+}}$. Then, from the definition of $q_{n}(t)$
and $r_{n}(t)$, $V_{2}(q_{n}(t);\theta)=e^{\theta r_{n}(t)}$, and similarly,
$V_{2}(q_{n};\theta)=e^{\theta r_{n}}$. The drift of $V_{2}(q;\theta)$,
denoted by $\Delta V_{2}(q;\theta)$ is defined as

$\displaystyle\mathbb{E}[\Delta
V_{2}(q_{n}(t);\theta)|q_{n}(t)=q]:=\sum_{q^{\prime}\in\mathbb{W}}G(q^{\prime},q)\left(V_{2}(q^{\prime};\theta)-V_{2}(q;\theta)\right).$

Then,

$\displaystyle\mathbb{E}[$ $\displaystyle\Delta
V_{2}(q_{n}(t);\theta)|q_{n}(t)=q]$
$\displaystyle=\lambda_{n}V_{2}(q+1;\theta)+\mu\min\\{n,q\\}V_{2}(q-1;\theta)-(\lambda_{n}+\mu\min\\{n,q\\})V_{2}(q;\theta)$
$\displaystyle=e^{\theta[n-q]^{+}}\left(\lambda_{n}\left(e^{-\theta}-1\right)\mathds{1}_{\\{q<n\\}}+\mu\min\\{n,q\\}\left(e^{\theta}-1\right)\mathds{1}_{\\{q\leq
n\\}}\right)$
$\displaystyle=e^{\theta[n-q]^{+}}\left(\lambda_{n}\left(e^{-\theta}-1\right)+\mu(n-[n-q]^{+})\left(e^{\theta}-1\right)\right)$
$\displaystyle\ \ \ \
+\lambda_{n}\left(1-e^{-\theta}\right)\mathds{1}_{\\{q\geq
n\\}}-n\mu\left(e^{\theta}-1\right)\mathds{1}_{\\{q>n\\}}$
$\displaystyle=(1-e^{-\theta})\Big{(}-\tilde{\gamma}_{n}(q;\theta)V_{2}(q;\theta)+\tilde{\beta}_{n}(q;\theta)\Big{)},$

where

$\displaystyle\tilde{\gamma}_{n}(q;\theta)$
$\displaystyle:=\left(\lambda_{n}-\mu(n-[n-q]^{+})e^{\theta}\right),$
$\displaystyle\tilde{\beta}_{n}(q;\theta)$
$\displaystyle:=\big{(}\lambda_{n}\mathds{1}_{\\{q\geq n\\}}-n\mu
e^{\theta}\mathds{1}_{\\{q>n\\}}\big{)}.$

Note that $V_{2}(q;\theta)\leq e^{\theta n}$. Thus, as $V_{2}(\cdot)$ is
uniformly bounded, in steady-state, for any $\theta\in\mathbb{R}$, we have

$\displaystyle\mathbb{E}[\Delta V_{2}(q_{n};\theta)]=0.$

and so,

$\displaystyle\mathbb{E}[\tilde{\gamma}_{n}(q_{n};\theta)V_{2}(q_{n};\theta)]=\mathbb{E}[\tilde{\beta}_{n}(q_{n};\theta)].$

This gives us that for any $\theta\in\mathbb{R}$,

$\displaystyle\left(\lambda_{n}e^{-\theta}-n\mu\right)\mathbb{E}\left[e^{\theta
r_{n}}\right]+\mu\mathbb{E}\left[r_{n}e^{\theta
r_{n}}\right]=\lambda_{n}e^{-\theta}\mathbb{P}(q_{n}\geq
n)-n\mu\mathbb{P}(q_{n}>n).$ (37)

Further, $\mathbb{E}[e^{\theta r_{n}}]=\mathbb{E}[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}]+\mathbb{P}(q_{n}\geq n)$, and

$\mathbb{E}[r_{n}e^{\theta r_{n}}]=\mathbb{E}[r_{n}e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}]=\frac{d}{d\theta}\mathbb{E}[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}].$

Substituting this in Eq. (37), we have

$\displaystyle\left(n(1-\epsilon_{n})e^{-\theta}-n\right)\mathbb{E}[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}]+\frac{d}{d\theta}\mathbb{E}[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}]=n\mathbb{P}(q_{n}=n).$

Solving the above differential equation gives us that

$\displaystyle\mathbb{E}\left[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}\right]=n\mathbb{P}(q_{n}=n)\tilde{G}_{n}^{-1}(\theta)\int_{-\infty}^{\theta}G_{n}(t)dt,$

where

$G_{n}(t)=\exp\left(-nt-(1-\epsilon_{n})(e^{-t}-1)\right)=\exp\Big{(}-n\epsilon_{n}t-(1-\epsilon_{n})(e^{-t}+t-1)\Big{)}.$

By substituting $\theta=0$, we get

$\displaystyle\mathbb{P}(r_{n}>0)=n\mathbb{P}(q_{n}=n)\int_{-\infty}^{0}G_{n}(t)dt.$
(38)

Thus, we finally get

$\displaystyle\mathbb{E}\left[e^{\theta
r_{n}}\mathds{1}_{\\{r_{n}>0\\}}\right]=\mathbb{P}(r_{n}>0)\tilde{G}_{n}^{-1}(\theta)\left(\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}\int_{-\infty}^{\theta}G_{n}(t)dt.$

This implies that, for any $\theta\in\mathbb{R}$,

$\displaystyle\mathbb{E}\left[e^{\theta
r_{n}}|r_{n}>0\right]=\tilde{G}_{n}^{-1}(\theta)\left(\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}\int_{-\infty}^{\theta}G_{n}(t)dt.$

Next, we calculate the steady state values of $\mathbb{P}(r_{n}>0)$,
$\mathbb{P}(w_{n}>0)$ and $\mathbb{P}(q_{n}=n)$. Note that

$\displaystyle\mathbb{P}(r_{n}>0)+\mathbb{P}(q_{n}=n)+\mathbb{P}(w_{n}>0)=\mathbb{P}(q_{n}<n)+\mathbb{P}(q_{n}=n)+\mathbb{P}(q_{n}>n)=1.$

We use this to solve for $\mathbb{P}(q_{n}=n)$. By using Eq. (36) and (38) to
get that

$\displaystyle\left(\frac{\lambda_{n}}{n\mu\epsilon_{n}}+1+n\int_{-\infty}^{0}G_{n}(t)dt\right)\mathbb{P}(q_{n}=n)=1.$

This gives us that

$\displaystyle\mathbb{P}(q_{n}=n)=\left(\frac{1}{\epsilon_{n}}+n\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}.$

Now, the result follows using Eq. (36) and (38). ∎

###### Proof of Theorem 8.

First note that for any $a,b$ and $s_{1},s_{2}$, we have

$\displaystyle\int_{s_{1}}^{s_{2}}$
$\displaystyle\exp\Big{(}at-\frac{1}{2}bt^{2}\Big{)}$
$\displaystyle=\sqrt{\frac{2\pi}{b}}\exp\Big{(}\frac{a^{2}}{2b}\Big{)}\Big{[}\Phi\Big{(}\sqrt{b}s_{2}-\frac{a}{\sqrt{b}}\Big{)}-\Phi\Big{(}\sqrt{b}s_{1}-\frac{a}{\sqrt{b}}\Big{)}\Big{]},$
(39)

where $\Phi(\cdot)$ denotes the CDF of a standard normal random variable.
First, we prove the result for the number of idle servers when the system is
in Super-HW regime. In this regime, $\epsilon_{n}=cn^{-\alpha}$ where
$\alpha>\frac{1}{2}$. For convenience, we use $\rho_{n}=1-\epsilon_{n}$ and
$\zeta_{n}:=\sqrt{\frac{n\epsilon_{n}^{2}}{\rho_{n}}}$. For system in Super-HW
regime,

$\displaystyle\zeta_{n}=\Big{(}\frac{c^{2}n^{1-2\alpha}}{1-cn^{-\alpha}}\Big{)}^{\frac{1}{2}}\leq
2cn^{\frac{1}{2}-\alpha},$

where the last inequality follows by assuming $n\geq 4c^{2}$. Further, if
$n^{\alpha-\frac{1}{2}}>2c$, we have $\zeta_{n}\leq 1$ Now,

$\displaystyle\int_{-\infty}^{\theta}$ $\displaystyle G_{n}(t)dt$
$\displaystyle=\int_{-\infty}^{\theta}\exp\Big{(}-n\epsilon_{n}t-n\rho_{n}(e^{-t}+t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle=\int_{-\theta}^{\infty}\exp\Big{(}n\epsilon_{n}t-n\rho_{n}(e^{t}-t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle=\int_{-\theta}^{\infty}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}\exp\Big{(}-n\rho_{n}\big{(}e^{t}-\frac{1}{2}t^{2}-t-1\big{)}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\int_{0}^{\infty}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\quad+\exp\Big{(}-\frac{1}{6}n\rho_{n}\theta^{3}\Big{)}\int_{-\theta}^{0}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\sqrt{\frac{2\pi}{n\rho_{n}}}\exp\Big{(}\frac{1}{2}\zeta_{n}^{2}\Big{)}\Phi\Big{(}\zeta_{n}\Big{)}\allowdisplaybreaks$
$\displaystyle\quad+\sqrt{\frac{2\pi}{n\rho_{n}}}\exp\Big{(}\frac{1}{2}\zeta_{n}^{2}-\frac{1}{6}n\rho_{n}\theta^{3}\Big{)}\Big{[}\Phi\Big{(}\sqrt{n\rho_{n}}\theta+\zeta_{n}\Big{)}-\Phi\Big{(}\zeta_{n}\Big{)}\Big{]},$

where (a) follows by using $e^{t}-\frac{1}{2}t^{2}-t-1\geq\frac{1}{6}t^{3}$
for $t<0$; and (b) follows by using Eq. (D.1). This also gives us that

$\displaystyle n\epsilon_{n}\int_{-\infty}^{0}$ $\displaystyle
G_{n}(t)dt\leq\sqrt{2\pi}\zeta_{n}\exp\Big{(}\frac{1}{2}\zeta_{n}^{2}\Big{)}\Phi\Big{(}\zeta_{n}\Big{)}\leq\sqrt{2\pi
e}\zeta_{n},$ (40)

where the last inequality follows by using $\zeta_{n}\leq 1$ whenever
$n^{\alpha-\frac{1}{2}}>2c$. This implies that, in Super-HW regime,

$\displaystyle\mathbb{P}(r_{n}>0)\stackrel{{\scriptstyle(a)}}{{=}}\frac{n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt}{1+n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt}\leq