2.1. Stability of the MSJ system under first-come, first-served

In the MSJ setting, the First-Come First-Served (first-come, first-served) policy with head-of-line blocking does not consistently keep all servers occupied, reducing its stability region relative to more optimized policies. Characterizing this stability region is an important line of work.

An important approach for characterizing the stability region of an MSJ system under first-come, first-served is the saturated system approach of [Reference Baccelli and Foss2]. This work provides a general framework under which one can relate the stability region of a queueing system that completes jobs in near- first-come, first-served order to the throughput of a saturated queueing system. In the MSJ system, jobs enter service in first-come, first-served order, causing the completion order to be very close to the first-come, first-served order. In the saturated system, an unlimited amount of jobs are available at all times, replacing the external arrival process.

This approach has been used to characterize the stability region of a variety of MSJ first-come, first-served systems by Grosof et al. [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15], though often requiring the solution of a system of linear equations to obtain the stability region. More explicit symbolic results have been proven for special cases of job duration distribution and server need, in the works of Grosof et al. [Reference Grosof, Harchol-Balter and Scheller-Wolf14], Rumyantsev and Morozov [Reference Rumyantsev and Morozov28], and [Reference Rumyantsev27].

These results have also been generalized to certain asymptotic limits, in the case of exponential duration and 1-or- $n$ server need in the work of [Reference Olliaro, Marsan, Balsamo and Marin24], specifically the limit as the $n$-server jobs complete much faster than the $1$-server jobs, and the limit as the number of servers $n$ diverges.

We study the same 1-and- $n$ exponential setting as [Reference Olliaro, Marsan, Balsamo and Marin24]. The asymptotic limits in that paper are simpler than our asymptotic limit—we allow both $n \to \infty$ and $p_n \to 0$, where $p_n$ is the probability that a job needs all $n$ servers, allowing us to explore a richer variety of asymptotic behavior. Our results cover both the stability region and the mean response time.

2.2. Response time of MSJ system under first-come, first-served

Much less is known about the MSJ first-come, first-served mean response time. The only setting in which mean response time has been exactly characterized is the setting with $n=2$ servers and where 1-server jobs require an exponential service duration [Reference Brill and Green3, Reference Filippopoulos and Karatza9].

A related but distinct model has been considered where a job’s service must start at the same time at each of its servers, but may end at different times on different servers. If completion times are independent across the servers, exact results on mean waiting time and the distribution of waiting time are known, see Green [Reference Green11] and Seila [Reference Seila29]. In contrast, the MSJ model requires that service ends at the same time across all of a job’s servers.

In more general MSJ settings, one must turn either to asymptotic results or numerical algorithms and approximations.

As for numerical results for response time, the standard technique is to employ matrix geometric and matrix analytic methods, and we see no reason why off-the-shelf methods would not be applicable, such as those described by Latouche and Ramaswami [Reference Latouche and Ramaswami19]. These methods were also used to derive some of the stability results discussed above, such as [Reference Rumyantsev and Morozov28]. Such numerical results are far more efficient than simulation, but are inherently restricted to examining a single set of parameters at a time. They cannot provide asymptotic results, which are the focus of this paper.

Mean response time, or equivalently mean waiting time, has been explored in several asymptotic limits. One type of limit that has received significant attention is a scaling limit where the number of servers, the largest server needed in the system, and the load all scale polynomially, mimicking the classical Halfin–Whitt limit, in the work of [Reference Wang, Xie and Harchol-Balter35] and [Reference Hong and Wang18]. However, results in this limit have been limited to the case where the load is low enough to allow stability with the maximal number of wasted servers: If the largest server need of a job is $m_{\max}$, this line of work requires that the mean number of servers occupied is below $n - m_{\max}$. This limit becomes trivial when $m_{\max} = n$, as in this work. Results in this limit focus on the asymptotic decay of queueing probability in [Reference Wang, Xie and Harchol-Balter35] or of mean waiting time in the scaling limit in [Reference Hong and Wang18].

We refer to this limit as the server-focused multilevel scaling (SFMS) limit. The SFMS limit has a much lighter load than the LFMS limit of this paper.

The remaining asymptotic work has focused on the fixed- $n$ limit, as the load approaches capacity. This limit mirrors the classical heavy-traffic limit. Tight asymptotic characterization of mean response time is known for MSJ first-come, first-served in this limit, proven by Grosof et al. [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15]. However, these results do not immediately generalize to any multilevel-scaling limit. This generalization is the contribution of this paper.

2.3. MSJ system under other scheduling policies

Outside of first-come, first-served scheduling, more results are known in the MSJ setting under nontrivial scheduling policies.

Throughput-optimality results have been established both in a preemptive setting, for the MaxWeight policy in the work of [Reference Maguluri, Srikant and Ying22], and in a nonpreemptive setting, for the Randomized Timers policy introduced by Psychas and Ghaderi [Reference Psychas and Ghaderi25], though the Randomized Timers policy has extremely poor response time in certain settings. The recent Markovian Service Rate policy introduced by Chen et al. [Reference Chen, Grosof and Berg6] uses a simple timer-based scheme to achieve throughput optimality, including in a nonpreemptive setting, with predictable and acceptable response time in all settings.

Moving to policies whose analyses primarily focus on mean response time, a significant line of work has focused on settings where all jobs’ server needs are powers of 2 and $n$ is a power of 2, and settings with related divisibility assumptions. First, the preemptive ServerFilling policy has been shown to match the response time of a resource-pooled M/G/1/ first-come, first-served, in the power-of-2 setting with a fixed number of servers $n$, studied by Grosof and Harchol-Balter [Reference Grosof, Harchol-Balter and Scheller-Wolf13], Grosof et al. [Reference Grosof and Harchol-Balter12]. In this setting, the ServerFilling-SRPT policy was shown by Grosof et al. [Reference Grosof, Scully, Harchol-Balter and Scheller-Wolf16] to achieve optimal mean response time in heavy traffic. Finally, in the power-of-2 setting, the Balanced Splitting dispatching policy was exactly analyzed by Anselmi and Doncel [Reference Anselmi and Doncel1], with both fixed- $n$ and scaling analysis results.

In the SFMS limit discussed in Section 2.2, the Smallest Need First policy has also been studied by Hong and Wang [Reference Hong and Wang18], resulting in an expanded set of asymptotic limits in which expected waiting time decays to 0.

Many practically used scheduling policies are backfilling policies, which can be seen as starting with the first-come, first-served scheduling policy, with additional steps taken to use capacity that first-come, first-served would otherwise leave unused, as discussed by Srinivasan et al. [Reference Srinivasan, Kettimuthu, Subramani and Sadayappan30], Carastan-Santos et al. [Reference Carastan-Santos, De Camargo, Trystram and Zrigui5], and Wang and Guo [Reference Wang and Guo34]. These policies have been investigated in simulation and in practice, but little is known about their theoretical properties. We hope to strengthen the theoretical understanding of the first-come, first-served policy, to build toward a theoretical understanding of backfilling policies.

3.1. Multiserver-job system

We consider an MSJ system with $n$ identical servers. Customers arrive according to a Poisson process with rate $\lambda$. We specifically consider a model in which there are two types of jobs: Jobs which need to be served by one server, and jobs that need all $n$ servers. The probability that a job needs one server is $p_1$, and (as defined in Section 2.1) the probability that a job needs all $n$ servers is $p_n = 1-p_1$, i.i.d. The service time for a 1-server job is $Exp(\mu_1)$, and for an $n$-server job is $Exp(\mu_n)$, both independent of all other jobs. Jobs are served in First-Come First-Served (first-come, first-served) order with head-of-line blocking, meaning that jobs are admitted into service in first-come, first-served order until a job is reached that cannot be served in the remaining available servers. No further jobs are served. We define $q$ to be the queue length at a particular point in time, and $Q$ to be the stationary queue length random variable.

Note that the symbol $p_n$ always refers to the fraction of jobs that use all of the servers, even as the number of servers $n$ changes.

Prior work (see, e.g., [Reference Grosof, Harchol-Balter and Scheller-Wolf14, Reference Rumyantsev27, Reference Rumyantsev and Morozov28], [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15], or [Reference Olliaro, Marsan, Balsamo and Marin24]) has demonstrated the existence of a stability threshold $\mu$, where the system is stable if and only if $\lambda \lt \mu$. Note that $\mu$ has been explicitly characterized by several methods, including analysis of the corresponding saturated system, which we define and discuss further in Section 3.2. We build on these prior explicit results to establish the scaling of the stability threshold $\mu$. We define $\rho = \lambda / \mu$ to be the fraction of the stability region that is in use.

A natural asymptotic limit is the $\lambda \to \mu$ scaling limit, known as the “heavy traffic” limit. The limiting value of the mean scaled queue length, $\lim_{\rho \to 1} \mathbb{E}[Q(1-\rho)]$, has been characterized by analyzing the saturated system, which we likewise discuss in Section 3.2.

However, the focus of this paper is the LFMS limit where we allow $\rho$ to grow quickly enough that $\mathbb{E}[Q(1-\rho)]$ reaches its limiting value, and only then scale the number of servers to study the secondary effects of that scaling.

When scaling the number of servers, it is natural to also consider scaling the fraction of single-server jobs $p_1$ toward 1, as the number of single-server jobs that can be served at once grows asymptotically.

Our goal is to characterize asymptotic behavior of the mean scaled queue length for this system under two secondary asymptotic limits: $p_1 \to 1$, and $n \to \infty$. Both of these asymptotic limits are secondary to our load scaling.

In particular, it makes sense to consider taking both limits simultaneously as there is a natural separation, based on whether the system load is dominated by 1-server jobs, by $n$-server jobs, or is balanced between the two.

We refer to the inherent workload of a job as the product of its server need and its duration. We refer to the overall load as the rate at which this workload is arriving to the system. The overall load of 1-server jobs is $\lambda \frac{p_1}{\mu_1}$, and the overall load of $n$-server jobs is $\lambda \frac{n p_n}{\mu_n}$. As a result, 1-server jobs require $\frac{p_1}{p_1 + n p_n}$ fraction of the total load.

We focus on the limit where $n \to \infty$ and $p_n \to 0$, holding $\mu_1$ and $\mu_n$ constant, and adjusting $\lambda$ as necessary to focus on the $\rho \to 1$ limit. As $p_1 \to 1$, the fraction of the total load required by 1 server jobs can be approximated as $\frac{1}{1 + n p_n}$.

In this paper, we define the notation $f(n) = o(g(n))$ and $f(n) = O(g(n))$ as follows:

\begin{align*} f(n) = o(g(n)) \Leftrightarrow \lim_{n \to \infty} \frac{f(n)}{g(n)} = 0, \quad f(n) = O(g(n)) \Leftrightarrow \lim_{n \to \infty} \frac{f(n)}{g(n)} \in (-\infty, \infty). \end{align*}

Note that when $f(n) = o(g(n))$ or $f(n) = O(g(n))$, $f(n)$ may be negative while $g(n)$ is positive.

We define $\omega(\cdot), \Omega(\cdot),$ and $\theta(\cdot)$ similarly. If $f(n) = o(g(n)),$ then $g(n) = \omega(f(n))$. If $f(n) = O(g(n)),$ then $g(n) = \Omega(f(n))$. If $f(n) = O(g(n))$ and $f(n) = \Omega(g(n))$, then $f(n) = \theta(g(n))$.

If $n p_n = \omega(1)$, or equivalently $p_n = \omega(\frac{1}{n})$, the fraction of load occupied by 1-server jobs converges to 0 as $p_n \to 0, n \to \infty$. We call this the “ $n$-server dominated” regime.

If $p_n = c/n$, for an arbitrary positive constant $c$ (which would imply that $n p_n = \theta(1)$, or equivalently $p_n = \theta(\frac{1}{n})$), the fraction of load occupied by 1-server jobs converges to a value strictly between 0 and 1. We call this the “balanced” regime.

Finally, if $n p_n = o(1)$, or equivalently $p_n = o(\frac{1}{n})$, the fraction of load occupied by 1-server jobs converges to 1. We call this the “1-server dominated” regime. We specifically focus on the polynomial-scaling subset of the 1-server dominated regime, where $p_n = 1/n^\alpha, \alpha \gt 1$.

Note that throughout the paper, all limits are taken as the number of servers $n$ tends to infinity.

3.2. Saturated system

A key tool for analyzing the MSJ system is the “saturated system.” The saturated system is a closed system with the same MSJ completion setup as the original open system. In the saturated system, whenever a server is unoccupied, and there is no job blocking the head of the line, jobs enter the system with the same i.i.d. probabilities $p_1$ and $p_n$ of 1 and $n$-server jobs, respectively. There is no separate stochastic arrival process. Jobs enter in response to a completion until either all $n$ servers are occupied, or an $n$-server job is generated that cannot be served immediately.

Saturated system analysis can be used to derive the stability region of the original MSJ system, because the stability region of the original system is bounded by the throughput of the saturated system, as demonstrated by Grosof et al. [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15], Baccelli and Foss [Reference Baccelli and Foss2], and Foss and Konstantopoulos [Reference Foss and Konstantopoulos10].

Saturated system analysis can also be used to derive the mean queue length of open MSJ systems under first-come, first-served service. To do so, [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15] makes use of the relative completions function $\Delta$, which maps states of the saturated system to real values. Note that the concept of relative arrivals and relative completions was also present in work prior to that result [Reference Dimitrov7, Reference Falin and Falin8].

To define $\Delta$, we must first define several quantities in the saturated system. Let $y$ denote a state of the saturated system and define $\mu_y$ as the completion rate in state $y$. Let $\nu_y$ denote the total transition rate out of state $y$, and $\nu_{y,y'}$ denote the specific transition rate from state $y$ to another state $y'$, where $y'$ is an arbitrary state of the saturated system. Finally, define $\mu$ as the time-average completion rate of the saturated system (e.g., the throughput).

Letting $y$ denote a state of the saturated system, $\Delta(y)$ is defined as the solution to the following system of equations:

(1)\begin{equation} \Delta(y) := \frac{\mu_y - \mu}{\nu_y} + \sum_{y'} \frac{\nu_{y,y'}}{\nu_y} \Delta(y'). \end{equation}

This equation (1) can be seen as the Poisson equation for a Markov Reward Process whose instantaneous reward is the completion rate $\mu_y$. This relative completion function impacts mean queue length in the MSJ system similarly to the impact of job size variance in an M/G/1 queue: It effectively measures the variability in completion rate. If the system is in a state $y$ with a large value of $\Delta(y)$, the system is likely to experience a high completion rate in the near future.

Note that the solution to this system of equations is only defined up to an additive constant. Letting $Y$ denote the time-average random variable for the saturated system state, [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15] adopts the convention that $\mathbb{E}[\Delta(Y)] = 0$, which results in a unique solution.

Let $Y_d$ denote the state-average random variable corresponding to the embedded discrete-time Markov chain (DTMC), updating at completion epochs. [Reference Grosof, Hong, Harchol-Balter and Scheller-Wolf15, Theorem 4.2] proves the following result, for an arbitrary MSJ first-come, first-served system, in the heavy traffic limit:

(2)\begin{equation} \lim_{\rho \to 1} \mathbb{E}[Q(1-\rho)] = \mathbb{E}[\Delta(Y_d)] + 1. \end{equation}

Our goal in this paper is to analyze the asymptotic growth rate of $\mathbb{E}[\Delta(Y_d)]$ in each of our three regimes, within the surrounding LFMS context.

5.1. Transition-average stationary distribution

The state of the saturated system in the transition-based embedded DTMC can be captured by a pair $(a, b)$ where $a$ denotes the number of $n$-server jobs present, $a \in \{0, 1\}$, and $b$ denotes the number of $1$-server jobs present, $0 \le b \le n$.

There are two kinds of states in this embedded DTMC: states where jobs are completing, and states where jobs are arriving. If $a = 0$ and $b \lt n$, then the state is an arriving state, because room for additional jobs is available, while if $a = 1$ or $b = n$, the state is a completing state. If $a = 1$ and $b = 0$, then an $n$-server job is in service, while otherwise only 1-server jobs are in service.

The transition probabilities are as follows:

• • From an arriving state $(0, b)$, where $b \lt n$, the arriving job is a 1-server job with probability $p_1$, resulting in a transition to $(0, b+1)$, or an $n$-server job with probability $p_n$, resulting in a transition to $(1, b)$.

• • From a completing state of the form $(1, b)$, where $1 \le b \lt n$, $b$ 1-server jobs are in service, resulting in a guaranteed transition to $(1, b-1)$.

• • From the state $(1, 0)$, the $n$-server job is in service, resulting in a transition to $(0, 0)$.

• • From the state $(0, n)$, $n$ 1-server jobs are in service, resulting in a transition to $(0, n-1)$.

Now, we characterize the transition-average stationary distribution $\pi$:

Lemma 5.1. The transition-average stationary distribution $\pi$ can be written as:

\begin{align*} \forall b \lt n-1, \pi_{0, b} &= \pi_{0, 0} p_1^b, \\ \pi_{0, n-1} &= \pi_{0, 0} \frac{p_1^{n-1}}{p_n}, \\ \pi_{0, n} &= \pi_{0, 0} \frac{p_1^n}{p_n}, \\ \forall b \lt n, \pi_{1, b} &= \pi_{0, 0} p_1^b. \end{align*}

Proof. Note that the arrivals occur in states $(0, b)$ for $0 \le b \lt n$, and job service completions occur in states $(0, n)$ and $(1, b)$ for $0 \le b \lt n$. Furthermore, in state $(1, 0)$, the $n$-server job is in service. In all other job service completion states, only $1$-server jobs are in service. Then we have the following balance equations

\begin{align*} \pi_{0, 0} &= \pi_{1, 0}, \\ \forall b, 0 \lt b \lt n-1, \pi_{0, b} &= \pi_{0, b-1} p_1, \\ \pi_{0, n-1} &= \pi_{0, n-2} p_1 + \pi_{0, n}, \\ \pi_{0, n} &= \pi_{0, n-1} p_1, \\ \forall b, 0 \le b \lt n-1, \pi_{1, b} &= \pi_{1, b+1} + \pi_{0, b} p_n, \\ \pi_{1, n-1} &= \pi_{0, n-1} p_n, \end{align*}

and the result follows.

5.2. Time-average stationary distribution

We now use our results from Section 5.1 to characterize the time-average stationary distribution.

Lemma 5.2. The time-average stationary distribution $P$ is given as

\begin{align*} P_{1, 0} &= C' \frac{1}{\mu_n},\\ \forall b, 1 \lt b \lt n, P_{1, b} &= \frac{1}{b \mu_1} = C' \frac{1}{b \mu_1} p_1^b, \\ P_{0, n} &= C' \frac{1}{n \mu_1} \frac{p_1^n}{p_n}, \end{align*}

where $\mu_1$ and $\mu_n$ is the service completion rate of size 1 and size $n$ jobs, respectively, and

\begin{align*} C' &= \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right)^{-1}. \end{align*}

Proof. The time-average stationary distribution is proportional to the transition-average stationary distribution, rescaled by the time spent in each state. The only states in which time is spent are completion states, and the average time per visit is the inverse of the completion rate. Thus, we define a proportionality constant $C$, such that $P_y = C \frac{\pi_y}{\mu_y}$, where recall that $\mu_y$ is the completion rate in state $y$. We also define $C' = C \pi_{0, 0}$. Now, applying Lemma 5.1, we have

\begin{align*} P_{1, 0} &= C\pi_{1, 0} \frac{1}{\mu_n} = C' \frac{1}{\mu_n},\\ \forall b, 1 \lt b \lt n, P_{1, b} &= C\pi_{1, b} \frac{1}{b \mu_1} = C' \frac{1}{b \mu_1} p_1^b, \\ P_{0, n} &= C\pi_{0, n}\frac{1}{n \mu_1} = C' \frac{1}{n \mu_1} \frac{p_1^n}{p_n}, \end{align*}

and $C'$ can be computed by using the fact that $P$ is a distribution:

\begin{align*} 1 &= \sum_y P_y = C' \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right) \end{align*}

Lemma 5.3. The explicit expression for $\mu$ is given as

\begin{align*} \mu &= \frac{1}{p_n} \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right)^{-1}. \end{align*}

Proof. We apply Lemma 5.2. Note that we have

\begin{align*} \mu &= P_{1, 0} \mu_n + P_{0, n} n \mu_1 + \sum_{b=1}^{n-1} P_{1, b} b \mu_1 \\ &= C' (1 + \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} p_1^b) = C' (1 + \frac{p_1^n}{p_n} + \frac{p_1-p_1^n}{1-p_1}) = C' \frac{1}{p_n}. \end{align*}

Plugging in the closed-form expression of $C'$ gives the desired result.

Note that semantically, $C'$ is the rate of completion of $n$-server jobs, and $\frac{1}{C'}$ is the expected time per $n$-server job completion. We have

\begin{align*} \frac{1}{C'} = \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}. \end{align*}

Intuitively, one can see that there are three phases to the completion of an $n$-server job, namely, time spent in state $(1, 0)$, time spent in state $(0, n)$, and time spent in the states $(1, b)$ for $0 \lt b \lt n$. The expected lengths of these phases are the first, second, and third terms, respectively, in the above expression.

5.3. Completion-average stationary distribution

Now we characterize the completion-average distribution $\pi^d$. Note that there are only $n+1$ possible states in which a completion can occur, with 0 to $n$ 1-server jobs in service.

Lemma 5.4. The completion-average stationary distribution $\pi^d$ is given as

\begin{align*} \forall b \lt n, \, \pi^d_{1, b} &= p_1^b p_n, \\ \pi^d_{0, n} &= p_1^n. \end{align*}

Proof. Note that the distribution $\pi^d$ is equal to the transition-average stationary distribution $\pi$, conditioned on the saturated system being in a completion state. In particular, there exists a normalizing constant $C^{\prime\prime}$ such that for any completing state $s$,

\begin{align*} \pi^d_s = C^{\prime\prime} \pi_s. \end{align*}

Using the expression for $\pi$ given in Lemma 5.1, we can solve for the proportionality constant $C^{\prime\prime}$ as simply $\frac{p_n}{\pi_{0, 0}}$, giving rise to the stated expression for $\pi^d$.

5.4. Relative completions

We now characterize the relative completions in the saturated system, which is key to the characterization of the scaling behavior of mean queue length in the original open system (see Section 3.2 for details).

Lemma 5.5. The explicit formula for $\mathbb{E}[\Delta(Y_d)]$ is

\begin{align*} \mathbb{E}[\Delta(Y_d)] = p_1^n \left(1 - \frac{\mu}{n \mu_1}\right) \left( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}\right) + \sum_{i=1}^{n-1} p_1^i p_n\left(1 - \frac{\mu}{i \mu_1}\right) \left(i - \frac{\mu}{\mu_1} H_i\right). \end{align*}

where $H_i=\sum_{k=1}^i\frac{1}{k}$ is the $i^{\mbox{th}}$ harmonic number.

Proof. For notational convenience, we will now switch to a one-dimensional state representation. More specifically, define the states of the time-average system and the completion-average embedded DTMC to be the integers $\{0, 1, 2, \ldots, n\}$, where the state is the number of 1-server jobs in the system. Note that we ignore the arrival states $(0, 0), (0, 1), \ldots, (0, n-1)$, as the system does not spend time in these states and completions cannot occur in them.

To characterize $\mathbb{E}[\Delta(Y_d)]$, we first characterize $\Delta$ up to an offset. Let $\widetilde{\Delta}(i) := \Delta(i) - \Delta(0)$, for any state $i$. We will use the following identity:

(3)\begin{equation} \mathbb{E}[\widetilde{\Delta}(Y)] = \mathbb{E}[\Delta(Y)] - \mathbb{E}[\Delta(0)] = - \mathbb{E}[\Delta(0)]. \end{equation}

The second equality holds because $\mathbb{E}[\Delta(Y)]$ is defined to be 0, as discussed in Section 3.2. Thus, we can characterize $\mathbb{E}[\Delta(Y_d)]$ as follows:

(4)\begin{align} \mathbb{E}[\widetilde{\Delta}(Y_d)] &= \mathbb{E}[\Delta(Y_d)] - \mathbb{E}[\Delta(0)] = \mathbb{E}[\Delta(Y_d)] + \mathbb{E}[\widetilde{\Delta}(Y)] \nonumber\\ \implies \mathbb{E}[\Delta(Y_d)] &= \mathbb{E}[\widetilde{\Delta}(Y_d)] - \mathbb{E}[\widetilde{\Delta}(Y)].\end{align}

Thus, it suffices to characterize $\widetilde{\Delta}$ and use (4) to characterize $\mathbb{E}[\Delta(Y_d)]$.

We can rewrite (4) more explicitly as

(5)\begin{equation} \mathbb{E}[\Delta(Y_d)] = \sum_{i=0}^n (\pi_i^d - P_i) \widetilde{\Delta}(i), \end{equation}

where note that $P_i = P(Y = i),$ and $\pi_i^d = P(Y_d = i)$. Recall from Lemma 5.4 that the distribution of $\pi^d$ is:

\begin{align*} \forall i \lt n, \pi_i^d &= p_1^i p_n, \\ \pi_n^d &= p_1^n. \end{align*}

Rewriting the formula for $P_i$ from Lemma 5.2, we have

\begin{align*} P_0 &= \mu \frac{p_n}{\mu_n}, \\ \forall 1 \le i \lt n, P_i &= \mu \frac{p_1^i p_n}{i \mu_1}, \\ P_n &= \mu \frac{p_1^n}{n \mu_1}. \end{align*}

Note that $\pi_i^d = \frac{\mu}{\mu_i} P_i$, where $\mu_i$ is the instantaneous completion rate in state $i$. Applying (5),

(6)\begin{align}\mathbb{E}[\Delta(Y_d)] &= (p_1^n - \mu \frac{p_1^n}{n \mu_1}) \widetilde{\Delta}(n) + \sum_{i=1}^{n-1} (p_1^i p_n - \mu \frac{p_1^i p_n}{i \mu_1}) \widetilde{\Delta}(i) \nonumber\\ &= p_1^n (1 - \frac{\mu}{n \mu_1}) \widetilde{\Delta}(n) + \sum_{i=1}^{n-1} p_1^i p_n(1 - \frac{\mu}{i \mu_1}) \widetilde{\Delta}(i).\end{align}

Next, we calculate $\widetilde{\Delta}(i)$. Using the Poisson equation (1), which characterizes $\Delta(i)$ up to an absolute constant, we obtain

\begin{align*} \forall 1 \le i \lt n, \Delta(i) &= 1 - \frac{\mu}{i \mu_1} + \Delta(i-1) \\ \implies \forall 1 \le i \lt n, \Delta(i) - \Delta(i-1) &= 1 - \frac{\mu}{i \mu_1} \\ \implies \forall 1 \le i \lt n, \widetilde{\Delta}(i) &= i - \frac{\mu}{\mu_1} H_i. \end{align*}

Note that state $n$ is a special case:

\begin{align*} \Delta(n) &= 1 - \frac{\mu}{n \mu_1} + p_1 \Delta(n) + p_n \Delta(n-1) \\ \implies p_n \Delta(n) &= 1 - \frac{\mu}{n \mu_1} + p_n \Delta(n-1) \\ \implies \Delta(n) - \Delta(n-1) &= \frac{1}{p_n} - \frac{\mu}{n p_n \mu_1} \\ \implies \widetilde{\Delta}(n) &= n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}. \end{align*}

Plugging the expressions of $\widetilde{\Delta}(i)$ for $i=1,\ldots,n$ into (6), we have

\begin{align*} \mathbb{E}[\Delta(Y_d)] &= p_1^n (1 - \frac{\mu}{n \mu_1}) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}) \\ & + \sum_{i=1}^{n-1} p_1^i p_n(1 - \frac{\mu}{i \mu_1}) (i - \frac{\mu}{\mu_1} H_i). \end{align*}

In the next three sections, we will characterize the asymptotic growth rate of the $\mu$ expression given in Lemma 5.3 and the $\mathbb{E}[\Delta(Y_d)]$ expression given in Lemma 5.5 for the three regimes discussed in Section 4.

8.1. Regime 3, general setting: Any growth rate $p_n = o(1/n)$

This setting contains a wide variety of behaviors, so while we characterize asymptotic behavior of the throughput $\mu$ in significant detail, our result on scaled mean queue length, $\mathbb{E}[\Delta(Y_d)]$, is less tight. Our tighter results are in the more specific setting of Section 8.2.

We start by characterizing the throughput $\mu$ in this regime.

Lemma 8.1. If $p_n = o(1/n \ln n)$, then

\begin{align*} \mu = (1+o(1)) n \mu_1. \end{align*}

On the other hand, if $p_n = \omega(1/(n \ln n))$ and $p_n = o(1/n)$, then

\begin{align*} \mu = (1+o(1)) \frac{\mu_1}{p_n \ln n}. \end{align*}

Proof. We start with Lemma 5.3:

\begin{align*} \mu &= \frac{1}{p_n} \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right)^{-1} \\ &= \frac{1}{p_n} \left( \frac{1}{\mu_n} + (1+o(1))\frac{1}{n \mu_1 p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right)^{-1}. \end{align*}

This step holds because $p_n = o(1/n)$, so $p_1^n = 1+o(1)$. As a result, we can similarly simplify the summation term—note that $p_1^n \le p_1^b \le 1$. Furthermore, $n p_n = o(1)$, so the $1/\mu_n$ term is negligible. Then

\begin{align*} \mu &= \frac{1+o(1)}{p_n} \left(\frac{1}{n \mu_1 p_n} + \sum_{b=1}^{n-1} \frac{1}{b \mu_1}\right)^{-1} = \frac{1+o(1)}{p_n} \left(\frac{1}{n \mu_1 p_n} + \frac{H_{n-1}}{\mu_1}\right)^{-1} \\ &= \frac{1+o(1)}{p_n} \left( \frac{1}{n \mu_1 p_n} + \frac{\ln n}{\mu_1}\right)^{-1} \end{align*}

Either of these two terms can be larger in asymptotic limit, depending on the asymptotic behavior of $p_n$. If $p_n = o(1/(n \ln n))$, then $\frac{1}{n \mu_1 p_n}$ is larger. If $p_n = \omega(1/(n \ln n))$, then $\frac{\ln n}{\mu_1}$ is larger. In the former case,

\begin{align*} \mu &= \frac{1+o(1)}{p_n} \left( \frac{1}{n \mu_1 p_n}\right)^{-1} = (1+o(1)) \frac{n \mu_1 p_n}{p_n} = (1+o(1)) n \mu_1. \end{align*}

In the latter case,

\begin{align*} \mu &= \frac{1+o(1)}{p_n} \left(\frac{\ln n}{\mu_1}\right)^{-1} = (1+o(1)) \frac{\mu_1}{p_n \ln n}. \end{align*}

Finally, if $p_n = \theta(1/(n \ln n))$, then both terms $\frac{1}{n \mu_1 p_n}$ and $\frac{\ln n}{\mu_1}$ are within a constant factor of each other in the asymptotic limit, and both contribute to the asymptotic behavior of $\mu$.

We now switch to $\mathbb{E}[\Delta(Y_d)]$, where we only focus on the case where $p_n = o(1/(n \ln n))$, and only provide an upper bound on the asymptotic growth of $\mathbb{E}[\Delta(Y_d)]$. Our stronger results for $\mathbb{E}[\Delta(Y_d)]$ are in Section 8.2.

Lemma 8.2. If $p_n = o(1/(n \ln n))$, then

\begin{align*} \mathbb{E}[\Delta(Y_d)] = o(1/p_n). \end{align*}

Proof. We start with the expression of $\mathbb{E}[\Delta(Y_d)]$ from Lemma 5.5, which we separate into two terms.

(28)\begin{equation} \mathbb{E}[\Delta(Y_d)] = p_1^n (1 - \frac{\mu}{n \mu_1}) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}) \end{equation}

(29)\begin{equation} \!\!+ \sum_{i=1}^{n-1} p_1^i p_n(1 - \frac{\mu}{i \mu_1}) (i - \frac{\mu}{\mu_1} H_i).\qquad \end{equation}

First consider (28). Note that by Lemma 8.1, in the setting where $p_n = o(1/(n \ln n))$, we have $\mu = (1+o(1)) n \mu_1$, and hence $1-\frac{\mu}{n \mu_1} = o(1)$. Furthermore, $p_1^n \le 1$. Thus,

\begin{align*} &p_1^n (1 - \frac{\mu}{n \mu_1}) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}) \\ &= o(1) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}). \end{align*}

Note that $\frac{1}{p_n} = \omega(n \ln n)$, so the dominant terms are the two terms involving $\frac{1}{p_n}$.

\begin{align*} ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}) &=\frac{1 - \frac{\mu}{n \mu_1}}{p_n} + O(n \ln n) = o(1/p_n). \end{align*}

Note that we only derive an upper bound on the growth rate of (28), namely $o(1/p_n)$. We do not have any lower bound on the growth rate of this term, due to the cancellation of the large quantities $1/p_n$ and $\mu/n p_n \mu_1$. Since these terms are of the same magnitude, we can only upper bound the result.

Now, consider (29). Note that $p_1^i = 1+o(1)$ for any $i \le n$ and $-\frac{n}{i} \le 1-\frac{\mu}{i \mu_1} \le 1$. Specifically, $1-\frac{\mu}{i \mu_1} = O(n/i)$. Furthermore, $-n \ln i \le i - \frac{\mu}{\mu_1} h_i \le i$. This gives the bound $i - \frac{\mu}{\mu_1} H_i = O(n \ln n).$

Thus, substituting these bounds into (29), we find that

\begin{align*} &\sum_{i=1}^{n-1} p_1^i p_n(1 - \frac{\mu}{i \mu_1}) (i - \frac{\mu}{\mu_1} H_i) \le O(p_n) \sum_{i=1}^{n-1} \frac{n^2 \ln n}{i} \\ =& O(p_n) n^2 \ln nH_{n-1} = O(p_n n^2 \ln^2 n) = o(n \ln n) = o(1/p_n). \end{align*}

8.2. Regime 3, specific setting: Polynomial scaling $p_n = 1/n^\alpha, \alpha \gt 1$

In this section, we focus on the specific setting of polynomial $p_n$ scaling where $p_n = 1/n^\alpha, \alpha \gt 1$. This allows us to obtain much stronger results. To prove Theorem 4.2, we start by characterizing the asymptotic behavior of throughput $\mu$ in Lemma 8.3:

Lemma 8.3. If $p_n = 1/n^\alpha$ for $\alpha \gt 1$, then

\begin{align*} \mu = n \mu_1 - \mu_1 n^{2 - \alpha} \ln n - (1 + o(1)) \mu_1 (\frac{\mu_1}{\mu_n} - 1) n^{2 - \alpha}. \end{align*}

Proof. Consider the general expression for $\mu$ from Lemma 5.3:

\begin{align*} \mu &= \frac{1}{p_n} \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right)^{-1}. \end{align*}

From Lemma 8.1, we know that $\mu = (1+o(1)) n \mu_1$. Note that for sufficiently large $n$, $\mu_n \lt n \mu_1$, and as a result $\mu \lt n \mu_1$, because the completion rate $\mu_i$ in state $i$ is at most $n \mu_1$ for all states $i$, for all $n \ge \frac{\mu_n}{\mu_1}$.

We therefore examine $\frac{n \mu_1}{\mu} - 1$:

(30)\begin{align} \frac{n \mu_1}{\mu} - 1 &= \left(n \mu_1 p_n \left( \frac{1}{\mu_n} + \frac{1}{n\mu_1} \frac{p_1^n}{p_n} + \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1}\right) - 1\right) \nonumber\\ &= \frac{n \mu_1 p_n}{\mu_n} + \frac{n \mu_1 p_n}{n\mu_1} \frac{p_1^n}{p_n} + n \mu_1 p_n \sum_{b=1}^{n-1} \frac{p_1^b}{b \mu_1} - 1 \nonumber\\ &= \frac{n^{1-\alpha} \mu_1}{\mu_n} + p_1^n - 1 + n^{1-\alpha} \sum_{b=1}^{n-1} \frac{p_1^b}{b}. \end{align}

We next expand $p_1^n = (1- n^{-\alpha})^n$ as

\begin{align*} (1- n^{-\alpha})^n = 1 - n n^{-\alpha} + \frac{n(n-1)}{2} n^{-2\alpha} \ldots \end{align*}

Because $\alpha \gt 1$, we can summarize this expansion as:

\begin{align*} p_1^n = 1 - n^{1-\alpha} + O(n^{2-2\alpha}). \end{align*}

Substituting the asymptotic behavior of $p_1^n$ into (30), we find that

\begin{align*} \frac{n^{1-\alpha} \mu_1}{\mu_n} + p_1^n - 1 + n^{1-\alpha} \sum_{b=1}^{n-1} \frac{p_1^b}{b} = \frac{n^{1-\alpha} \mu_1}{\mu_n} - n^{1-\alpha} + O(n^{2-2\alpha}) + n^{1-\alpha} \sum_{b=1}^{n-1} \frac{p_1^b}{b}. \end{align*}

Note also that $p_1^b \in [p_1^n, 1]$ for all $b \le n$. In particular, $p_1^b = 1+ O(n^{1-\alpha})$ (uniformly) for all $b \le n$.

We can therefore simplify the summation as

\begin{align*} n^{1-\alpha} \sum_{b=1}^{n-1} \frac{p_1^b}{b} = n^{1-\alpha} (1 + O(n^{1-\alpha})) \sum_{b=1}^{n-1} \frac{1}{b} = n^{1-\alpha}\ln n + O(n^{2-2\alpha} \ln n)). \end{align*}

Putting it together, we find that

(31)\begin{align} \frac{n \mu_1}{\mu} - 1 &= \frac{n^{1-\alpha} \mu_1}{\mu_n} - n^{1-\alpha} + O(n^{2-2\alpha}) + n^{1-\alpha}\ln n + O(n^{2-2\alpha})) \nonumber\\ &= n^{1-\alpha} \left(\ln n + \frac{\mu_1}{\mu_n} - 1\right) + O(n^{2-2\alpha} \ln n).\end{align}

Let $x$ denote the growth rate of the final expression (31), and use it to solve for $\mu$:

\begin{align*} \frac{n \mu_1}{\mu} - 1 &= x \implies n \mu_1 - \mu = x \mu \implies n \mu_1 = (1 + x) \mu \implies \frac{n \mu_1}{1 + x} = \mu. \end{align*}

Note that for any $x'$ near zero, $\frac{1}{1+x'} = 1 - x' + O(x'^2)$. Furthermore, since $\alpha \gt 1$, $x$ is in fact near 0. Thus, we have

\begin{align*} \mu = \frac{n \mu_1}{1 + x} &= n \mu_1 (1 - x + O(x^2)) \\ &=n \mu_1 - n \mu_1 \left(n^{1-\alpha} \left(\ln n + \frac{\mu_1}{\mu_n} - 1\right) + O(n^{2-2\alpha} \ln n) + O(n^{2-2\alpha} \ln^2 n) \right) \\ &= n \mu_1 - \mu_1 n^{2-\alpha} \ln n - \mu_1\left(\frac{\mu_1}{\mu_n} - 1\right) n^{2 - \alpha} + O(n^{3-2\alpha}\ln^2 n). \end{align*}

Since $\alpha \gt 1$, we have $3 - 2\alpha \lt 2 - \alpha$, so our final expression matches our theorem statement with a more specific error bound.

We are now ready to prove the remainder of Theorem 4.2 by characterizing the asymptotic mean queue length $\mathbb{E}[\Delta(Y_d)]$.

Theorem 8.1. If $p_n = 1/n^\alpha$ for $\alpha \gt 1$, then

\begin{align*} \mathbb{E}[\Delta(Y_d)] = (1+o(1))\frac{1}{2} n^{2 - \alpha} \ln^2 n. \end{align*}

Proof. We start by splitting up $\mathbb{E}[\Delta(Y_d)]$ into five terms, as before:

(32)\begin{equation} \mathbb{E}[\Delta(Y_d)] = p_1^n (1 - \frac{\mu}{n \mu_1}) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}) \end{equation}

(33)\begin{equation} +\, p_n \sum_{i=1}^{n-1} p_1^i i\qquad\qquad\qquad\qquad\qquad \end{equation}

(34)\begin{equation} \!\!-\, \frac{\mu p_n}{\mu_1} \sum_{i=1}^{n-1} p_1^i\qquad\qquad\qquad\qquad\quad \end{equation}

(35)\begin{equation} \!\!\!\!-\, \frac{\mu p_n}{\mu_1} \sum_{i=1}^{n-1} p_1^i H_i\qquad\qquad\qquad\qquad \end{equation}

(36)\begin{equation} \!\!\!+\, \frac{\mu^2 p_n}{\mu_1^2} \sum_{i=1}^{n-1} p_1^i \frac{1}{i} H_i\qquad\qquad\qquad\quad \end{equation}

We will demonstrate the following asymptotic growth rates for each of these terms:

1. (1) Non-summation term (32): $(1+o(1)) (\frac{\mu_1}{\mu_n} - \gamma) n^{2-\alpha} \ln n$.

2. (2) $p_1^i i$ term (33): $O(n^{2-\alpha})$.

3. (3) $p_1^i$ term (34): $- (1+o(1)) p_1 n^{2 - \alpha}$.

4. (4) $p_1^i H_i$ term (35): $-n^{2-\alpha} \ln n + O(n^{2-\alpha})$.

5. (5) $p_1^i H_i / i$ term (36): $\frac{1}{2}n^{2-\alpha}\ln^2 n + O(n^{2-\alpha}\ln n)$.

Note that only Item 1, the asymptotic growth rate for (36), is non-negligible. As a result, these growth rates sum to $\frac{1}{2}n^{2-\alpha}\ln^2 + O(n^{2-\alpha} \ln n)$, as desired.

We start with the non-summation term (32):

\begin{align*} p_1^n (1 - \frac{\mu}{n \mu_1}) ( n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1}). \end{align*}

We will consider the three multiplicative terms one by one, left to right. First, note that $p_1^n = 1+o(1)$. Then, we expand $1-\frac{\mu}{n \mu_1}$. Using Lemma 8.3:

\begin{align*} 1 - \frac{\mu}{n \mu_1} = \frac{n \mu_1 - \mu}{n \mu_1} = n^{1-\alpha}\ln n + (1+o(1)) (\frac{\mu_1}{\mu_n} -1) n^{1-\alpha}. \end{align*}

Next expand the final (rightmost) parenthesis of (32). Again using Lemma 8.3, we have:

(37)\begin{align} &n-1 + \frac{1}{p_n} - \frac{\mu}{\mu_1} H_{n-1} - \frac{\mu}{n p_n \mu_1} = n-1 - \frac{\mu}{\mu_1} H_{n-1} + n^\alpha (1-\frac{\mu}{n \mu_1}) \nonumber\\ =&\, n-1 - \frac{\mu}{\mu_1} H_{n-1} + n \ln n + (1+o(1)) n (\frac{\mu_1}{\mu_n} -1) \nonumber\\ =&\, n-1 - (n - \frac{\ln n}{n^{2-\alpha }} - (1 + o(1)) (\frac{\mu_1}{\mu_n} - 1)) n^{2 - \alpha}) H_{n-1} + n \ln n + (1+o(1)) n (\frac{\mu_1}{\mu_n} -1). \end{align}

Note that $2 - \alpha \lt 1$. Let $\epsilon \gt 0$ be a constant such that $2 - \alpha \lt 1-\epsilon \lt 1$. As a result, we can find the asymptotic growth rate of (37) as follows:

\begin{align*} & n-1 - (n - n^{2 - \alpha} \ln n - (1 + o(1)) (\frac{\mu_1}{\mu_n} - 1)) n^{2 - \alpha}) H_{n-1} + n \ln n + (1+o(1)) n (\frac{\mu_1}{\mu_n} -1)\\ =&\, n-1 - (n - O(n^{1-\epsilon})) H_{n-1} + n \ln n + (1+o(1)) n (\frac{\mu_1}{\mu_n} -1)\\ =&\, n-1 - (n - O(n^{1-\epsilon})) (\ln n + \gamma + O(1/n)) + n \ln n + (1+o(1)) n (\frac{\mu_1}{\mu_n} -1))\\ =&\, n-1 + O(n^{1-\epsilon}) (\ln n + \gamma + O(1/n)) - n (\gamma + O(1/n))+ (1+o(1)) n (\frac{\mu_1}{\mu_n} -1)\\ =&\, n(\frac{\mu_1}{\mu_n} - \gamma) +o(n). \end{align*}

Combining these results, we find that the overall growth rate of (32) is

\begin{align*} (1+o(1)) (\frac{\mu_1}{\mu_n} - \gamma) n^{2-\alpha} \ln n, \end{align*}

confirming Item 1.

We next consider (33) which can be quantified exactly as:

\begin{align*} \sum_{i=1}^{n-1} p_1^i p_n i &= p_n \frac{p_1 - n p_1^n + (n-1) p_1^{1 + n}}{(1 - p_1)^2} = \frac{p_1 - n p_1^n +n p_1 p_1^n - p_1 p_1^n}{p_n} \\ &= \frac{p_1 ( 1-p_1^n) - n p_1^n (1-p_1)}{p_n} = \frac{p_1 ( 1-p_1^n) - n p_1^n (1-p_1)}{p_n}. \end{align*}

We now expand $p_1^n = (1-p_n)^n$ using a binomial expansion, and recalling that $n p_n = o(1)$,

\begin{align*} &\frac{p_1 ( 1-p_1^n) - n p_1^n (1-p_1)}{p_n} = \frac{p_1 (n p_n + O(n^2 p_n^2)) - n p_n (1 - n p_n + O(n^2 p_n^2))}{p_n} \\ &= \frac{-(1-p_1) n p_n + O(n^2 p_n^2)}{p_n} = \frac{-n p_n^2 + O(n^2 p_n^2)}{p_n} = \frac{O(n^2 p_n^2)}{p_n} = O(n^2 p_n) = O(n^{2 - \alpha}), \end{align*}

confirming Item 2.

Next consider (34) which again has an exact expression

\begin{align*} -\sum_{i=1}^{n-1} p_1^i p_n \frac{\mu}{\mu_1} &= - \frac{\mu}{\mu_1} p_n \frac{p_1-p_1^n}{1-p_1} = -\frac{\mu}{\mu_1}(p_1 - p_1^n) = -\frac{\mu}{\mu_1} p_1(1-p_1^{n-1}) \\ &= - \frac{\mu}{\mu_1}p_1 (1+o(1)) n^{1-\alpha} = - p_1 (1+o(1)) n^{2 - \alpha}, \end{align*}

confirming Item 3.

We next consider (35):

(38)\begin{align} -\sum_{i=1}^{n-1} p_1^i p_n \frac{\mu}{\mu_1} H_i &= -\sum_{i=1}^{n-1} p_1^i p_n \frac{\mu}{\mu_1} \ln i + O(n^{2-\alpha}) \nonumber\\ &= -p_n \frac{\mu}{\mu_1} \int_{i=1}^n p_1^i \ln i di + O(n^{2-\alpha}) \nonumber\\ &= -p_n \frac{\mu}{\mu_1} \int_{i=1}^n (1-n^{-\alpha})^i \ln i di + O(n^{2-\alpha}). \end{align}

As usual, we convert the $H_i$ to $\ln i$ and the summation to an integral. To bound the change caused by switching from $H_i$ to $\ln i$, note that $|H_i - \ln i| \le 1$, which gives rise to an error bounded by (34), which is negligible compared to the growth rate we will establish for this term.

The change caused by switching from the summation to the integral is similarly bounded using the fact that $\ln i - \ln (i-1) \le 1$ (which again introduces an error bounded by (34)).

Now, we upper and lower bound the integral as follows:

\begin{align*} \int_{i=1}^n p_1^i \ln i di &\le \int_{i=1}^n \ln i di = n \ln n - n +1,\\ \int_{i=1}^n p_1^i \ln i di &\ge \int_{i=1}^n p_1^n \ln i di = (n \ln n - n +1) p_1^n. \end{align*}

Using a binomial expansion of $p_1^n$, we have

\begin{align*} p_1^n = (1 - n^{-\alpha})^n &= 1 - n n^{-\alpha} + \frac{n(n-1)}{2} n^{-2\alpha} - \ldots = 1 - O(n^{1-\alpha}). \end{align*}

We can therefore analyze the bound above as:

\begin{align*} (n \ln n - n +1) p_1^n &= (n \ln n - n +1) (1- O(n^{1-\alpha})) \\ &= n \ln n - n - O(n^{2-\alpha} \ln n) = n \ln n - O(n) \\ \implies \int_{i=1}^n p_1^i \ln i di &= n \ln n - O(n). \end{align*}

Substituting the integral into (38), multiplying by $-p_n \mu/\mu_1 = n^{1-\alpha} + O(n^{2-2\alpha} \ln n)$, and using Lemma 8.3, we get

\begin{align*} -\sum_{i=1}^{n-1} p_1^i p_n \frac{\mu}{\mu_1} H_i = -n^{2-\alpha} \ln n + O(n^{2-\alpha}), \end{align*}

which matches Item 4, the asymptotic growth rate for (35).

Finally, consider (36). As usual, we convert the $H_i$ to $\ln i$ and the summation to an integral:

(39)\begin{equation} \sum_{i=1}^{n-1} p_1^i p_n \frac{\mu^2}{i \mu_1^2} H_i = p_n \frac{\mu^2}{\mu_1^2} \int_{i=1}^n p_1^i \frac{\ln i}{i} di + O(n^{2-\alpha} \ln n). \end{equation}

To justify the error term, note that both conversions introduce an error in the order of

\begin{align*} p_n \frac{\mu^2}{\mu_1} \sum_{i=1}^{n-1} p_1^i \frac{1}{i}. \end{align*}

Making use of our bound from Lemma 8.3, which states that $\mu = n \mu_1 + O(n^{1-\alpha} \ln)$, and the earlier shown result $p_1^i = 1 + O(n^{1-\alpha})$, we have

\begin{align*} p_n \frac{\mu^2}{\mu_1} \sum_{i=1}^{n-1} p_1^i \frac{1}{i} = O(n^{2-\alpha} \ln n), \end{align*}

which is the error term we used in (39).

We upper and lower bound the integral from (39) as follows:

(40)\begin{align} \int_{i=1}^n p_1^i \frac{\ln i}{i} di &\le \int_{i=1}^n \frac{\ln i}{i} di = \frac{\ln^2 n}{2}, \nonumber\\ \int_{i=1}^n p_1^i \frac{\ln i}{i} di &\ge \int_{i=1}^n p_1^n \frac{\ln i}{i} di = p_1^n \frac{\ln^2 n}{2}.\end{align}

Using our prior result $p_1^n = 1 + O(n^{1-\alpha})$, we get

\begin{align*} \int_{i=1}^n p_1^i \frac{\ln i}{i} di &\ge ( 1 + O(n^{1-\alpha})) \frac{\ln^2 n}{2} = \frac{\ln^2 n}{2} + O(n^{1-\alpha}\frac{\ln^2 n}{2}). \end{align*}

Combining this upper bound with our lower bound in (40), we obtain

(41)\begin{equation} \int_{=1}^n p_1^i \frac{\ln i}{i} di = \frac{\ln^2 n}{2} + O(n^{1-\alpha}\frac{\ln^2 n}{2}). \end{equation}

We now need to incorporate the $p_n \frac{\mu^2}{\mu_1^2}$ term. From Lemma 8.3, we have $\frac{\mu}{\mu_1} = n + O(n^{2-\alpha} \ln n)$. Then

(42)\begin{equation} p_n \frac{\mu^2}{\mu_1^2} = n^{-\alpha} (n + O(n^{2-\alpha} \ln n))^2 = n^{-\alpha}(n^2 + O(n^{3-\alpha} \ln n)) = n^{2-\alpha} + O(n^{3 - 2 \alpha} \ln n). \end{equation}

We now combine (41) and (42) to get our result:

\begin{align*} \sum_{i=1}^{n-1} p_1^i p_n \frac{\mu^2}{i \mu_1^2} H_i &= \frac{1}{2} n^{2-\alpha} \ln^2 n + O(n^{2-\alpha} \ln n), \end{align*}

which matches the claimed asymptotic growth rate for (36), and completing the proof.

10.1. Calculations and simulation for 1-and- $n$ system

In Figure 5, we see both insights verified in the 1-and- $n$ system, via direct calculation of throughput $\mu$ and the scaled queue length $\mathbb{E}[\Delta(Y_d)]$ using Lemmas 5.3 and 5.5. Each curve represents a fixed number of servers $n$ and a varying parameter $\alpha$ controlling the fraction of large jobs $p_n = \frac{1}{n^\alpha}$. We scale throughput $\mu$ and scaled queue length $\mathbb{E}[\Delta(Y_d)]$ by a factor of $n$ to make the curves comparable.

Figure 5. Main setting with server needs $1$ and $n$, $\mu_1 = \mu_n = 1$, and $p_n = 1/n^\alpha$. (a) Throughput $\mu/n$, calculated using Lemma 5.3. (b) Scaled mean queue length $\mathbb{E}[\Delta(Y_d)]/n$, calculated using Lemma 5.5.

In Figure 5a, we observe that throughput increases monotonically with $\alpha$, with an inflection point around $\alpha \approx 1.2$, as predicted by Theorems 4.1 and 4.2.

In Figure 5b, we observe that scaled mean queue length increases with $\alpha$ to a peak around $\alpha=1.2$ and then decreases back to 0 as $\alpha$ gets large, also as predicted by Theorems 4.1 and 4.2.

In Figure 6a, we observe via simulation that these insights are not restricted to calculations or to very large $n$. Here, we simulate $n=10$ servers, with each curve representing a different fraction of capacity demanded, using the following formula:

\begin{align*} \text{Fraction of capacity} := \frac{\mu_1 p_1 + \mu_n p_n n}{n} \end{align*}

Figure 6. Simulation results ( $10^8$ jobs per data point) for mean queue length when $p_{10} = 1/n^\alpha$ and $\lambda$ is chosen so that a constant fraction of capacity is utilized. (a) Main setting with server needs $1$ and 10, $n=10$ servers and $\mu_1 = \mu_{10} = 1$. Equal-area setting with server needs $1$ and $10$, $n=10$ servers and $\mu_1 = 1, \mu_{10} = 10$.

The bell-curve behavior of mean queue length is clearly visible, peaking around $\alpha=0.8$, and is especially clear for capacity fractions around the boundary of the stability region. Our bell-curve insight thus applies not only for scaled mean queue length $\mathbb{E}[\Delta(Y_d)]$, but also for mean queue length $\mathbb{E}[Q]$ itself, including in systems with as few as $n=10$ servers.

10.2. Calculations and simulation for duration-scaling system

In this section, we investigate the realistic scenario where low-resource jobs also require much lower service duration than high-resource jobs. This is realistic in datacenters, for instance, where a single-server analytics job might only take a few seconds, while an LLM-training job might take most of the datacenter for a period of months. To model this scenario, we set the service rate $\mu_1$ of the 1-server jobs equal to $n$, the number of servers, while leaving the server rate $\mu_n$ of the $n$-server jobs constant.

In Figure 7, we again observe that both of our insights are verified even more clearly in the duration-scaling system where throughput $\mu$ and the scaled queue length $\mathbb{E}[\Delta(Y_d)]$ are computed using Lemmas 5.3 and 5.5. Note that we use the scaling $\mu/n^2$ and $\mathbb{E}[\Delta(Y_d)]/n^2$ to account for the fact that the 1-server jobs have both $1/n$ the mean duration and $1/n$ the resource demand of the $n$-server jobs.

Figure 7. Scaled service rate setting with server needs $1$ and $n$, $\mu_1 = n, \mu_n = 1$ and $p_n = 1/n^\alpha$. (a) Throughput $\mu/n^2$, calculated using Lemma 5.3. (b) Scaled mean queue length $\mathbb{E}[\Delta(Y_d)]/n^2$, calculated using Lemma 5.5.

Figure 7a shows that throughput $\mu$ has an inflection point at $\alpha=2$, and Figure 7b shows that scaled mean queue length has a particularly symmetrical bell curve centered on $\alpha=2$ as well. We expect the inflection point to be around $\alpha=2$, because this is where the load of $n$-server and $1$-server jobs is equalized. Thus, our insights are verified in this setting as well, despite leaving the domain of applicability of asymptotic results.

Note that Figure 6b shows that our insights again carry over into simulation and into low- $n$ settings, with mean queue length $\mathbb{E}[Q]$ peaking around $\alpha =1.9$, across a wide range of fractions of capacity for $n=10$ servers. Hence, the bell-curve behavior is even smoother and more consistent in this duration-scaling setting than in Section 10.1.

10.3. Calculations and simulation for half-size large jobs

In this section, we investigate the case where the large-resource jobs require $n/2$ servers, rather than $n$ servers as in our asymptotic results. This mimics the definition of leadership-class jobs on the Frontier supercomputer at Oak Ridge National Labs [23], which consists of jobs that require at least 60% of the capacity of the supercomputer, giving top priority to such large-resource jobs and setting policy goals to attain a large fraction of load devoted to leadership-class jobs.

With half-size large jobs, we can no longer use the exact results in Lemmas 5.3 and 5.5 to compute throughput and scaled queue length. Instead, we compute the throughput $\mu$ using the product-form stationary distribution formula for the MSJ 2-class saturated system from [Reference Grosof, Harchol-Balter and Scheller-Wolf14, Section 5], and [Reference Céline and Isaac4, Section 5.1]. Figure 8a shows the result of this calculation, with throughput increasing monotonically with increasing $\alpha$, as the mixture of load moves from $n$-server dominated to 1-server dominated, with an inflection point around $\alpha=1$.

Figure 8. Calculations and simulation for the half-size large jobs setting with server needs $1$ and $n/2$, $\mu_1 = \mu_{n/2} = 1$ and $p_n = 1/n^\alpha$. (a) Throughput $\mu/n$, calculated using [Reference Grosof, Harchol-Balter and Scheller-Wolf14, Section 5], [Reference Céline and Isaac4, Section 5.1]. (b) Simulation results ( $10^8$ jobs per data point) for mean queue length when $\lambda$ is chosen so that a constant fraction of capacity is utilized and $n=10$ servers.

Since we do not have an exact formula for limiting scaled mean queue length $\mathbb{E}[\Delta(Y_d)]$, we instead turn to simulation in Figure 8b, again with $n=10$ servers. Here we see the same bell-curve behavior of mean queue length $\mathbb{E}[Q]$ near the boundary of the stability region, with a peak around $\alpha=0.4$. Note that in this setting, the fraction of capacity in use is significantly higher than in any of the other settings around the boundary of the stability region for $\alpha$ near that peak—around $78\%$ of capacity in this setting, compared to $50-60\%$ of capacity in the other settings. This follows from the fact that this system will never leave $n/2$ or more servers empty, while settings with jobs that demand $n$ servers may leave up to $n-1$ servers empty. Note also that in this setting, the peak of mean queue length $\mathbb{E}[Q]$ does not occur at the same $\alpha$ value as the inflection point of throughput, in contrast to Sections 10.2 and 10.1. We leave the analysis of this behavior to future work.

10.4. Calculations and simulation for three-class system

In this section, we investigate the case where there are three classes of jobs, requiring $1, n/2,$ and $n$ servers, respectively. This is the most realistic setting because in practice, queueing systems where jobs require varying resources are not limited to two possibilities, but instead vary widely.

In this setting, it is more difficult to calculate throughput $\mu$ than in the half-size large jobs setting of Section 10.3. There is an exact product-formula stationary distribution formula for this setting, given by Rumyantsev and Morozov [Reference Rumyantsev and Morozov28], as well as [Reference Grosof, Harchol-Balter and Scheller-Wolf14, Section 3]. However, the number of terms in this formula is exponential in $n$, while the formula in Lemma 5.3 and the formula used in Section 10.3 each contained a number of terms that was linear in $n$. We leave the problem of efficiently calculating $\mu$ in this setting to future work—a formula with a number of terms that is linear in $n$ is likely doable. No exact formula for $\mathbb{E}[\Delta(Y_d)]$ is known in this setting.

We now turn to simulation in Figure 9, again with $n=10$ servers. There is a larger parameter space of probability distributions over the three classes, but we choose arbitrarily to equalize the fraction of jobs in each of the two larger-resource classes. Here we see the same bell-curve behavior of mean queue length $\mathbb{E}[Q]$ near the boundary of the stability region, with a peak around $\alpha=0.7$, between the peaks in the original setting in Section 10.1 and the half-size large jobs setting in Section 10.3. The fraction of capacity in use is likewise between those two settings around the boundary of the stability region near that peak. Even for this more complicated setting, our key bell-curve insight is again confirmed.

Figure 9. Simulation results ( $10^8$ jobs per data point) for three-class setting with server needs $1, 5$, $10$, $n=10$ servers, $\mu_1 = \mu_5= \mu_{10} = 1$, and $p_5 = p_{10} = n^\alpha/2$.