The Join-the-Shortest-Queue System in the Halfin-Whitt Regime: Rates of Convergence to the Diffusion Limit

1. Introduction

Consider a queueing system with n identical servers, each with a finite buffer of length b. Customers arrive according to a Poisson process with rate $n\lambda$, and service times are independent and identically distributed (i.i.d.), exponentially distributed with rate one. Customers cannot change servers after the initial routing decision, and a customer arriving to a system where all servers are busy and all buffers are full is blocked. This is known as a parallel-server system. A load-balancing policy specifies the manner in which arriving customers are assigned to the servers. In this paper, we consider the classical join-the-shortest-queue (JSQ) policy. Under JSQ, an arriving customer enters service immediately if at least one server is idle; if not, they get routed to the server with the smallest number of customers in its buffer. Ties are broken arbitrarily. We refer to this as the JSQ system.

Parallel-server systems have generated immense interest in recent years, and the JSQ policy is fundamental because it minimizes the expected customer delay and maximizes, with respect to stochastic order, the number of customers served in a given time interval; see, for instance, Winston (1977) and Weber (1978). For a sample of recent work on the JSQ policy, we refer readers to Eryilmaz and Srikant (2012), Mukherjee et al. (2016), Eschenfeldt and Gamarnik (2018), Banerjee and Mukherjee (2019), Gupta and Walton (2019), Liu and Ying (2019), Banerjee and Mukherjee (2020), Braverman (2020), Zhou and Shroff (2020a, b), Cao et al. (2021), Hurtado-Lange and Maguluri (2022), and Zhao et al. (2021). Other popular load-balancing policies include the join-the-idle-queue policy (Stolyar 2015, Mukherjee et al. 2016), the idle-one-first policy (Gupta and Walton 2019), and of course the power-of-d policy (Vvedenskaya et al. 1996, Mitzenmacher 2001); but in this paper, we focus on the JSQ policy. We make no attempt to give a comprehensive review of the literature on parallel-server systems, instead referring the reader to van der Boor et al. (2021) for a recent survey.

Understanding the exact performance of the system is known to be difficult, and much attention has been devoted over the past decade to “heavy-traffic” asymptotics. The term heavy-traffic refers to parameter regimes where the system utilization tends to one. “Conventional heavy traffic” assumes that the number of servers n is fixed and $\lambda \uparrow 1$, whereas “many-server heavy traffic” assumes that $n\to \infty$ and $\lambda \uparrow 1$ jointly. For two examples of work in the conventional heavy-traffic setting, see Eryilmaz and Srikant (2012) and Zhou and Shroff (2020b). In this paper, we use the term heavy-traffic to refer to the many-server setting—the setting considered in most of the papers mentioned in the previous paragraph.

There are multiple many-server heavy-traffic regimes, depending on how n and λ jointly converge to their limit. For example, assuming that $\lambda =1-1/n$ yields entirely different asymptotic behavior compared to when $\lambda =1-1/\sqrt{n}$. To capture all the possible heavy-traffic regimes, it is common practice to assume that the per-server load λ is related to the number of servers n through $\lambda =1-\beta /n^{\alpha }\in (0,1)$ for some $\alpha \ge 0$ and $\beta >0$. In this paper, we focus on the case when $\alpha =1/2$; that is, $\lambda =1-\beta /\sqrt{n}$. This regime is known as the Halfin-Whitt regime and is ubiquitous across the queueing theory literature. It derives from the work of Halfin and Whitt (1981) and is also known as the quality-and-efficiency-driven regime because it achieves reasonable customer wait times while maintaining high utilization of servers. The full list of parameter regimes is found in Figure 1.

We now state and discuss our main results. Let $Q_i(t)$ be the number of servers with i or more customers at time $t\ge 0$, noting that $Q_i(t)=0$ for $i>b+1$. The process $\{Q(t)=(Q_1(t),\dots ,Q_{b+1}(t))\}$ is an irreducible continuous-time Markov chain (CTMC) on a finite state space and therefore possesses a unique stationary distribution. We let $Q=(Q_1,\dots ,Q_{b+1})$ be the random vector having the stationary distribution of the CTMC. To describe the asymptotic behavior of Q, we let $\delta =1/\sqrt{n}$ and define the diffusion-scaled random vector $X=(X_1,\dots ,X_{b+1})$ by $X_1=\delta (n-Q_1)$, and $X_i=\delta Q_i$ for $2\le i\le b+1$. The results of Eschenfeldt and Gamarnik (2018) and Braverman (2020) imply that X converges in distribution to some limiting ${\mathbb{R}}_{+}^{b+1}$-valued random vector Y as $n\to \infty$. In this paper, we establish an upper bound of order $1/\sqrt{n}$ on the rate of convergence to Y.

The random variable Y is distributed according to the stationary distribution of the diffusion process $\{Y(t)\in {\mathbb{R}}_{+}^{b+1}\}$, which satisfies

where $\{W(t)\}$ is standard Brownian motion and $\{U(t)\}$ is the unique nondecreasing, nonnegative process in the space of càdlàg functions $D[0,\infty )$ satisfying ${\int }_0^{\infty }1(Y_1(t)>0)dU(t)=0$. The diffusion $\{Y(t)\}$ was shown to be positive recurrent; see Banerjee and Mukherjee (2019) or Braverman (2020). Furthermore, (1) implies that $Y_3=\cdots =Y_{b+1}=0$.

Our main result is that there exists a constant $C(b,\beta )$ such that for all $n\ge 1$, and any function $h:{\mathbb{R}}_{+}^{b+1}\to \mathbb{R}$ whose first-order and second-order partial derivatives are bounded in magnitude by one,

The assumption that b is finite is used frequently in the proof of (2) and, specifically, in the proof of Proposition 1. We deem the finite buffer assumption to be acceptable because it was shown by Braverman (2020) that even with infinite-sized buffers, $\mathbb{E}{Q}_3\le C(\beta )$ for all $n\ge 1$ in the Halfin-Whitt regime, implying that $X_3\Rightarrow 0$ or that the mass concentrates on those states with at most one customer waiting. Moreover, Liu and Ying (2020) showed that assuming finite buffers, $\mathbb{E}{Q}_3\to 0$ as $n\to \infty$ in the even busier super-Halfin-Whitt regime ($1/2<\alpha <1$).

In addition to the novelty of our result, this paper makes a methodological contribution. We prove (2) using Stein’s method, a framework introduced by Stein (1972) that allows one to study the rate of convergence of a sequence of random variables to its limit. Popularized in the area of queueing systems by Gurvich (2014), Ying (2017), Braverman and Dai (2017), and Gast (2017), the generator comparison approach of Stein’s method, attributed to Barbour (1988, 1990) and Götze (1991), is used to study convergence rates of steady-state Markov chain distributions to their diffusion, fluid, or mean-field limits. For a few recent applications of the generator comparison approach in queueing, we refer the reader to Gaunt and Walton (2020), Hurtado-Lange and Maguluri (2022), Lu (2021), and Liu et al. (2022); this list is by no means comprehensive. In this paper, we restrict our attention to the case when the limit is the stationary distribution of a diffusion process, referring the reader to Ying (2017) for a treatment of fluid and mean-field limits.

The generator approach requires bounds on various moments of the prelimit, known as moment bounds, and bounds on the derivatives of the solution to the Poisson equation for the limiting distribution. The latter are called gradient bounds in Braverman and Dai (2017). But in this paper, we stick with the original term “Stein factors” or “Stein factor bounds”; see, for example, Ross (2011). Although moment bounds can be difficult to obtain in some applications, Stein factor bounds are typically the bigger problem. When the limit is one dimensional, Stein factors are bounded using the explicit form of the solution to the Poisson equation—an ordinary differential equation. When the limit is multidimensional, the Poisson equation is a partial differential equation (PDE) that generally does not have an explicit solution, making Stein factor bounds harder to establish. Techniques proposed to obtain multidimensional Stein factor bounds include using a priori Schauder estimates from elliptic PDE theory as in Gurvich (2014), using couplings to analyze and bound the sensitivity of the diffusion to its initial condition as in Barbour (1988) and Mackey and Gorham (2016), and bounding the Stein factors using Malliavin calculus as in Fang et al. (2018) and Jin et al. (2022). A detailed description of these techniques can be found in section 1.1 of Braverman (2022). However, despite progress on multidimensional Stein factor bounds, the JSQ system is not covered by existing results because our limiting diffusion in (1) is constrained to the nonnegative orthant via reflecting boundary conditions.

To deal with the Stein factor bound problem, this paper promotes the use of the prelimit generator comparison approach, which was recently proposed by Braverman (2022) as an alternative to the generator comparison approach. The prelimit approach is the mirror image of the classical generator approach. Whereas the latter requires moment bounds on the prelimit X and Stein factor bounds for limit Y, the former needs moment bounds on Y and Stein factor bounds for the prelimit X. For the moment bounds used in this paper, the result that all moments of Y are finite, proved by Banerjee and Mukherjee (2019), is sufficient because our limit Y does not depend on n. The Stein factor bounds pose a bigger challenge, and we deal with them in Section 3. It was noted in Braverman (2022) that the prelimit and classical generator comparison approaches should be equivalent, in theory, in the sense that any bound on $|\mathbb{E}h(X)-\mathbb{E}h(Y)|$ obtained using one of them should be attainable using the other. However, in practice, one approach could be more tractable, or convenient, to work with; see, for instance, the example in section 4 of Braverman (2022). In the case of the JSQ system, we discuss in Remark 1 of Section 3.2.3 how the discrete state space simplifies the analysis of the couplings we use to establish Stein factor bounds because the initial spacing of the coupled systems is preserved until coupling.

The introduction of the prelimit approach in Braverman (2022) was intended to be gentle, with the only example used there being the $M/M/1$ system. Our application of the approach to the JSQ system exposes all of its moving pieces and can be useful to those who want to apply the prelimit approach to their own setting. For example, some of the technical components of this paper that could be useful in other settings include the regenerative argument used to establish first-order Stein factor bounds in Section 3.1, the approach we use to bound $\mathbb{E}|X|$ in Section 3.2.1, and our treatment of reflecting boundary conditions in Section A.2.2 in Appendix A.

It should be noted that Hurtado-Lange and Maguluri (2022) and Zhou and Shroff (2020a) used the classical generator comparison approach to obtain rates of convergence of the steady-state total customer count to an exponential random variable for $\alpha >2$. The former paper was in the continuous-time setting, whereas the latter considered the discrete-time system; the results in both papers also hold for routing policies other than JSQ, such as the power-of-d policy. Because the limiting random variable in both papers is one dimensional, the Stein factors bounds do not pose a challenge there.

1.1. Literature Review

Let us first review the literature on the analysis of the JSQ system in the various many-server heavy-traffic regimes. Most of the work has been done in the setting with infinite buffer sizes, so, unless otherwise noted, we assume that $b=\infty$. In Eschenfeldt and Gamarnik (2018), the authors established the process-level convergence of ${\{X(t)\}}_{n=1}^{\infty }$ to its diffusion limit in the Halfin-Whitt regime ($\alpha =1/2$). That paper triggered a wave of interest in the many-server heavy-traffic asymptotics of the JSQ system. Convergence of the stationary distributions was later established by Braverman (2020), and the behavior of the stationary distribution of the limiting diffusion was studied by Banerjee and Mukherjee (2019, 2020). Our work fits with this group of papers, elevating the steady-state convergence result to one with rates of convergence.

Outside the Halfin-Whitt regime, Mukherjee et al. (2016) studied the transient and steady-state behavior of the JSQ system’s fluid limit when $\lambda =1-\beta <1$ is a fixed constant (α = 0); Gupta and Walton (2019) established process-level convergence to the diffusion limit when α = 1, known as the nondegenerate slowdown (NDS) regime and introduced by Atar (2012). In the sub-Halfin-Whitt regime when $\alpha \in (0,1/2)$, Liu and Ying (2019) assumed finite buffers and obtained bounds on the steady-state total customer count in the system. A similar result was obtained for Coxian-2 service times by Liu et al. (2022) and by Liu and Ying (2020) for the super-Halfin-Whitt regime $\alpha \in (1/2,1)$. Another recent work in the super-Halfin-Whitt regime was by Zhao et al. (2021), who worked with infinite buffers and established transient and steady-state diffusion limits for the normalized total queue length process. Their analysis exploited the regenerative structure of the JSQ system and contained several hitting-time estimates very close to our own estimates needed for the Stein factor bounds in Section 3. Lastly, both Hurtado-Lange and Maguluri (2022) and Zhou and Shroff (2020a) established rates of convergence to the exponential distribution for the steady-state normalized total customer count. Their results covered the case when $\alpha >2$.

Other works have used Stein’s method in the setting of parallel-server systems beyond Hurtado-Lange and Maguluri (2022) and Zhou and Shroff (2020a). In Liu and Ying (2019, 2020) and Liu et al. (2022), the authors used Stein’s method for mean-field analysis to obtain bounds on steady-state performance metrics of interest, like $\mathbb{E}{Q}_2$ for instance, for the power-of-d system. Another line of work on power-of-d systems was by Gast (2017), Gast and Van Houdt (2017), and Gast et al. (2019), where the authors showed how to derive refined mean-field models for improved steady-state approximations. More recently, Hairi and Ying (2021) provide calculable error bounds for the mean-field approximation of the power-of-two-choices model.

1.2. Notation

We use $\mathbb{Z}$ to denote the set of integers and let $\mathbb{N}=\{0,1,2,\dots \}$. For any $k\in \mathbb{N}$ and $B\subset {\mathbb{R}}^d$, we let $C^k(B)$ be the set of all k-times continuously differentiable functions $f:B\to \mathbb{R}$. We let $e\in {\mathbb{R}}^d$ be the vector whose elements all equal one and let $e^{(i)}$ be the element with one in the ith entry and zeros otherwise. For any $\delta >0$ and integer d > 0, we let $\delta {\mathbb{Z}}^d=\{\delta k:\hspace{0.17em}k\in {\mathbb{Z}}^d\}$ and define $\delta {\mathbb{N}}^d$ similarly. For any function $f:\delta {\mathbb{Z}}^d\to \mathbb{R}$, we define the forward difference operator in the ith direction as

for $j\ge 0$, we define with the convention that ${\mathrm{\Delta }}_i^{0}f(\delta k)=f(\delta k)$. For a vector $a\in {\mathbb{N}}^d$, we also let if $f:{\mathbb{R}}^d\to \mathbb{R}$, then and we adopt the convention that ${\partial }^{0}f(x)/\partial x^0=f(x)$. For any $x\in {\mathbb{R}}^d$, we define ${\Vert x\Vert }_1={\sum }_{i=1}^d|x_i|$ and use $|x|$ to denote the Euclidean norm. For any $f:{\mathbb{R}}^d\to \mathbb{R}$, we let ${\Vert f\Vert }_{\infty }={\sup }_{x\in {\mathbb{R}}^d}|f(x)|$. Throughout the paper, we will often use C to denote a generic positive constant that may change from line to line and that is independent of any parameters not explicitly specified.

2. Main Result

Recall that $Q_i(t)$ is the number of servers with i or more customers at time $t\ge 0$ and that ${\{Q(t)={(Q_i(t))}_{i=1}^{b+1}\}}_{t\ge 0}$ is an irreducible CTMC with state space given by

Figure 2 gives an example of a state $q\in S_Q$.

We assume that $\lambda =1-\beta /\sqrt{n}$ for some fixed $\beta >0$. Let $\delta =1/\sqrt{n}$ and define the diffusion-scaled CTMC $\{X(t)\}$ by

which takes values on the state space

We will often use $x^q\in S$ and $q\in S_Q$ interchangeably. Recalling that ${\mathrm{\Delta }}_{i}f(x^q)=f(x^q+\delta e^{(i)})-f(x^q)$, for any $f:S\to \mathbb{R}$, the infinitesimal generator of $\{X(t)\}$ satisfies

The first line of transitions in (5) corresponds to arrivals. We see that for $j\ge 2$, the jth component of $x_j^q$ only grows provided the preceding j – 1 horizontal levels, as depicted in Figure 2, are full. The transitions in the second line of (5) correspond to service completions. Using Figure 2 again, we interpret $(q_j-q_{j+1})$ as the number of servers (vertical columns) with exactly j customers.

Recall that $X=(X_1,\dots ,X_{b+1})$ and $Y=(Y_1,Y_2,0,\dots ,0)$ are distributed according to the stationary distributions of the scaled CTMC and the diffusion $\{Y(t)\in {\mathbb{R}}_{+}^{b+1}\}$ defined in (1), respectively. Going forward, we note that unless explicitly stated, all expectations are with respect to the stationary distribution at hand, that is, either X or Y. To state our main result, we define

and $d_{{\mathcal{M}}_j}(X,Y)={\sup }_{h^{\ast }\in {\mathcal{M}}_j}|\mathbb{E}{h}^{\ast }(X)-\mathbb{E}{h}^{\ast }(Y)|$. We use an asterisk to emphasize that $h^{\ast }(x)$ is defined on the continuum ${\mathbb{R}}^{b+1}$. Later we will drop the asterisk to refer to functions defined only on the grid $\delta {\mathbb{Z}}^{b+1}$. It was shown in lemma 2.2 of Mackey and Gorham (2016) that ${\mathcal{M}}_3$ is a convergence-determining class; that is, $d_{{\mathcal{M}}_3}(U,V)\to 0$ implies U and V converge in distribution. The following is our main result.

Theorem 1.

For any $0<b<\infty$, there exists a constant $C(b,\beta )$ such that for all $n\ge 1$,

$$d_{{\mathcal{M}}_2}(X,Y)=\underset{h^{\ast }\in {\mathcal{M}}_2}{\sup }|\mathbb{E}{h}^{\ast }(X)-\mathbb{E}{h}^{\ast }(Y)|\le C(b,\beta )/\sqrt{n}.$$(6)

Note that ${\mathcal{M}}_2$ is also a convergence-determining class because ${\mathcal{M}}_3\subset {\mathcal{M}}_2$. We prove Theorem 1 in Section 2.1 using the prelimit generator approach of Stein’s method. Multiple parts of the proof assume that n is large enough, say, $n>N(\beta )$ for some $N(\beta )>0$. We can make this assumption without loss of generality by redefining $C(b,\beta )$ to be larger than ${\max }_{1\le n\le N(\beta )}{d}_{{\mathcal{M}}_2}(X,Y)$.

2.1. Proving Theorem 1

Central to our proof is the ability to extend any grid-valued function to be defined on all of ${\mathbb{R}}_{+}^{b+1}$. Although there are infinitely many such extensions, we use a polynomial spline A that extends grid-valued functions $f:\delta {\mathbb{N}}^{b+1}\to \mathbb{R}$ to functions $Af:{\mathbb{R}}_{+}^{b+1}\to \mathbb{R}$. We leave the detailed construction to Section A.2 in Appendix A because for this section, it suffices to know that A is a linear operator, that $Af\in C^2({\mathbb{R}}_{+}^{b+1})$, and that A applied to a constant equals that constant. Recalling that $\delta =1/\sqrt{n}$, the following auxiliary lemma is needed.

Lemma 1.

Define

$${\mathcal{M}}_{\mathit{disc},j}(c)=\{h:\delta {\mathbb{N}}^{b+1}\to \mathbb{R},\hspace{0.17em}|{\mathrm{\Delta }}^{a}h(\delta k)|\le c{\delta }^{{\Vert a\Vert }_1},\hspace{0.17em}1\le {\Vert a\Vert }_1\le j,\hspace{0.17em}\delta k\in \delta {\mathbb{N}}^{b+1}\}.$$

There exist some $C,C^{\prime }>0$ independent of any JSQ model parameters such that

$$d_{{\mathcal{M}}_2}(X,Y)\le \underset{h\in {\mathcal{M}}_{\mathit{disc},2}(C)}{\sup }|\mathbb{E}h(X)-\mathbb{E}Ah(Y)|+C^{\prime }\delta .$$(7)

Proof of Lemma 1.

The result follows by repeating the arguments used in the proof of lemma 1 in Braverman (2022). □

Going forward, when we write ${\mathcal{M}}_{\mathit{disc},2}(C)$, the constant C is assumed to be the one in Lemma 1. Furthermore, note that if $h(0)\ne 0$, then the linearity of A and the fact that A applied to a constant equals that constant implies that $\tilde{h}(x)=h(x)-h(0)$ satisfies $\mathbb{E}\tilde{h}(X)-\mathbb{E}A\tilde{h}(Y)=\mathbb{E}h(X)-\mathbb{E}Ah(Y)$. We therefore, without loss of generality, consider only those $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$ such that $h(0)=0$.

To prove Theorem 1, we bound the right-hand side of (7) with the help of the following two ingredients. The first ingredient is a rate-conservation law for $\{Y(t)\}$, proved in Appendix A.

Lemma 2.

Given $f\in C^2({\mathbb{R}}_{+}^{b+1})$, define

$$G_{Y}f(x)=(\beta -(x_1+x_2))\frac{\partial }{\partial x_1}f(x)-x_2\frac{\partial }{\partial x_2}f(x)+\frac{{\partial }^2}{\partial x_1^2}f(x),\hspace{1em}x\in {\mathbb{R}}_{+}^{b+1}.$$(8)

If $\mathbb{E}|f(Y)|<\infty$ and $\mathbb{E}|G_{Y}f(Y)|<\infty$, and if Y(0) is initialized according to Y, then

$$\mathbb{E}{G}_{Y}f(Y)+\mathbb{E}\left({\displaystyle {\int }_0^1\left(\frac{\partial }{\partial x_1}f(Y(s))+\frac{\partial }{\partial x_2}f(Y(s))\right)}1(Y_1(s)=0)dU(s)\right)=0.$$(9)

The second ingredient is the Poisson equation. For $h:\delta {\mathbb{N}}^{b+1}\to \mathbb{R}$ and $c\in \mathbb{R}$, let

which is well defined because the CTMC has a finite state space and is therefore exponentially ergodic. Furthermore, lemma 2 of Braverman (2022) (see also lemma 1 of Barbour 1988) implies that

Most applications of Stein’s method have c = 0, but we choose $c=c^{\ast }=-f_h^{(0)}(0)$ and define

Our choice of c yields $f_h(0)=0$, which comes in handy later when we need to bound $|f_h(x^q)|$ in Proposition 1. Going forward, we assume that $c=c^{\ast }$ when referring to (10).

Let us give an informal roadmap for bounding (7), with the formal statement of the bounds left to Proposition 2. We bound (6) by comparing the CTMC and diffusion generators. However, the former is defined only on a subset of ${\mathbb{R}}_{+}^{b+1}$, which requires the following workaround. Suppose that we are given a set $B\subset {\mathbb{R}}_{+}^{b+1}$ such that (a) $\mathbb{E}h(X)-Ah(x)=A{G}_X{f}_h(x)$ for $x\in B$ and (b) the probability that $Y\notin B$ goes to zero rapidly (we will make this precise) as $n\to \infty$. We decompose $\mathbb{E}h(X)-Ah(x)$ as

and take expected values with respect to Y (we will show that these are finite) to get

Now extend $f_h(x^q)$ to $\delta {\mathbb{N}}^{b+1}$ by defining $f_h(x^q)=0$ for $x^q\in \delta {\mathbb{N}}^{b+1}\mathrm{\setminus }S$ and consider $A{f}_h(x)$. Provided that $\mathbb{E}|A{f}_h(Y)|<\infty$ and $\mathbb{E}|G_{Y}A{f}_h(Y)|<\infty$, we can invoke Lemma 2 with $f(x)=A{f}_h(x)$ there to conclude that

where Y(0) in the third line is initialized according to Y. We bound the first line by showing that G_X and G_Y are close to one another. The middle term is small because of our choice of B, and the last term can be bounded because the JSQ system exhibits reflecting behavior similar to $\{Y(t)\}$ at the boundary $\{x\in S:x_1^q=0\}$. As a final remark, our choice of $f_h(x^q)=0$ for $x^q\in \delta {\mathbb{N}}^{b+1}\mathrm{\setminus }S$ is made for convenience and is not essential to the proof because the probability that $Y\notin B$ shrinks rapidly as $n\to \infty$.

To state the following proposition, define $k:{\mathbb{R}}^{b+1}\to {\mathbb{Z}}^{b+1}$ elementwise by $k_j(x)=\lfloor x_j/\delta \rfloor$. For notational convenience, we also define $I=\{i=(i_1,i_2,0,\dots ,0)\in {\mathbb{N}}^{b+1}:\hspace{0.17em}0\le i_1,i_2\le 4\}$. The following proposition is proved in Section A.2 in Appendix A.

Proposition 2.

If $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$, then Ah(Y), $A{f}_h(Y)$, and $G_{Y}A{f}_h(Y)$ are integrable and (12) holds. Furthermore, suppose that n > 16, define

$$B=\{(x_1,x_2,0,\dots ,0)\in {\mathbb{R}}_{+}^{b+1}:x_2+x_1\le \delta (n/2-8)=(n/2-8)/\sqrt{n}\},$$

and let

$$\begin{array}{ll}\hfill & {\epsilon }_1(Y)=(A{G}_X{f}_h(Y)-G_{Y}A{f}_h(Y))1(Y\in B),\hfill \\ \hfill & {\epsilon }_2(Y)=(\mathbb{E}h(X)-Ah(Y)-G_{Y}A{f}_h(Y))1(Y\notin B),\hfill \\ \hfill & {\epsilon }_3(Y)=\left(\frac{\partial }{\partial x_1}A{f}_h(Y)+\frac{\partial }{\partial x_2}A{f}_h(Y)\right)1(Y\in B), \hspace{0.17em}\hspace{0.17em}\mathit{and} \hfill \\ \hfill & {\epsilon }_4(Y)=\left(\frac{\partial }{\partial x_1}A{f}_h(Y)+\frac{\partial }{\partial x_2}A{f}_h(Y)\right)1(Y\notin B).\hfill \end{array}$$

There exist $C(\beta ),C(b,\beta )>0$ independent of h(x) and n such that

$$\begin{array}{ll}|{\epsilon }_1(Y)|\le \hfill & C(\beta )(1+{\delta }^{-1}{Y}_2)\underset{\begin{array}{c}i\in I\\ a_1+a_2=2\end{array}}{\max }|{\mathrm{\Delta }}_1^{a_1}{\mathrm{\Delta }}_2^{a_2}{f}_h(\delta (k(Y)+i))|+C(\beta ){\delta }^{-2}\underset{i\in I}{\max }|{\mathrm{\Delta }}_1^3{f}_h(\delta (k(Y)+i))|\hfill \\ \hfill & +C(\beta ){\delta }^{-2}1(Y_1\le \delta )\underset{\begin{array}{c}i\in I\\ i_1=0\end{array}}{\max }|({\mathrm{\Delta }}_1^2-({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2))f_h(\delta (k(Y)+i))|,\hfill \\ |{\epsilon }_2(Y)|\le \hfill & C(b,\beta )1(Y\notin B){\delta }^{-2}(1+Y_1+Y_2)\underset{i\in I}{\max }|f_h(\delta (k(Y)+i))|,\hfill \\ |{\epsilon }_3(Y)|\le \hfill & C(\beta ){\delta }^{-1}1(Y\in B)(|({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2)f_h(\delta k(Y)|+\underset{\begin{array}{c}i\in I\\ a_1+a_2=2\end{array}}{\max }|{\mathrm{\Delta }}_1^{a_1}{\mathrm{\Delta }}_2^{a_2}{f}_h(\delta (k(Y)+i))|),\hfill \\ |{\epsilon }_4(Y)|\le \hfill & C(\beta ){\delta }^{-1}1(Y\notin B)\underset{i\in I}{\max }|f_h(\delta (k(Y)+i))|.\hfill \end{array}$$

Note that ${\epsilon }_1(Y)$ and ${\epsilon }_2(Y)$ are related to the first and second lines of (12), respectively, whereas ${\epsilon }_3(Y)$ and ${\epsilon }_4(Y)$ are related to the last line there. From the bounds in Proposition 2, we see that the bound on (12) depends on the CTMC through the function $f_h(x^q)$ and its differences and on the diffusion through the distribution of Y. The differences of $f_h(x^q)$ are commonly known as Stein factors; the following proposition, proved in Section 3, exhibits the Stein factor bounds we need to prove Theorem 1.

Proposition 1.

There exists $C(\beta ,b)>0$ such that for any $n\ge 1$ and $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$,

$$|{\mathrm{\Delta }}_1^{a_1}{\mathrm{\Delta }}_2^{a_2}{f}_h(x^q)|\le C(\beta ,b){\delta }^{a_1+a_2}{(1+x_2^q)}^{a_1+a_2},$$

for all $a_1,a_2\ge 0$ with $1\le a_1+a_2\le 2$, and all $x^q\in S$ with $x_1^q\le \delta (n-a_1),\hspace{0.17em}{x}_2^q\le \delta (n-a_2)$, and $x_3^q=0$. Furthermore,

$$\begin{array}{lll}\hfill & |f_h(x^q)|\le C(\beta ,b)(1+x_2^q)(x_1^q+x_2^q)/\delta ,\hfill & \hspace{1em}{x}^q\in S,\hspace{0.17em}{x}_3^q=0,\hfill \\ \hfill & |{\mathrm{\Delta }}_1^3{f}_h(x^q)|\le C(\beta ,b){\delta }^3{(1+x_2^q)}^3,\hfill & \hspace{1em}{x}^q\in S,\hspace{0.17em}{x}_1^q\le \delta (n-3),\hspace{0.17em}{x}_3^q=0,\hfill \end{array}$$

and for all $x^q\in S$ with $x_1^q=0,\hspace{0.17em}0\le x_2^q\le \delta (n-1)$, and $x_3^q=0$,

$$\begin{array}{ll}|({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2)f_h(x^q)|\le \hfill & C(\beta ,b){\delta }^2{(1+x_2^q)}^2\hspace{1em} \mathit{and} \hfill \\ |({\mathrm{\Delta }}_1^2-({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2))f_h(x^q)|\le \hfill & C(\beta ,b){\delta }^3{(1+x_2^q)}^3.\hfill \end{array}$$

The last component needed for the proof of Theorem 1 is the following lemma.

Lemma 3.

All moments of Y₁ and Y₂ are finite. Furthermore, suppose that Y(0) is initialized according to Y. Then for any j > 0,

$$\mathbb{E}{Y}_2^{j+1}=({\displaystyle {\int }_0^1(Y_2(}s){)}^{j}1(Y_1(s)=0)dU(s)).$$(13)

Proof of Lemma 3.

The finiteness of the moments follows from theorem 2.1 of Banerjee and Mukherjee (2019), and (13) is implied by (9) of Lemma 2 with $f(y)=y_2^{j+1}$ there. □

Proof of Theorem 1.

Initialize Y(0) according to Y. Using (12) and the definitions of ${\epsilon }_1(Y),\dots ,{\epsilon }_4(Y)$, it follows that

$$\mathbb{E}h(X)-\mathbb{E}Ah(Y)=\mathbb{E}{\epsilon }_1(Y)+\mathbb{E}{\epsilon }_2(Y)-\mathbb{E}({\displaystyle {\int }_0^1(}{\epsilon }_3(Y(s))+{\epsilon }_4(Y(s)))1(Y_1(s)=0)dU(s)).$$

We argue that $|\mathbb{E}h(X)-\mathbb{E}Ah(Y)|\le C(b,\beta )\delta$ for any $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$, which implies Theorem 1 when combined with Lemma 1. Because $\delta (k_2(Y)+i_2)\le Y_2+4\delta$ for $i\in I$, applying the Stein factor bounds in Proposition 1 with the bounds on ${\epsilon }_1(Y)$ and ${\epsilon }_2(Y)$ in Proposition 2 yields

$$|{\epsilon }_1(Y)|\le C(b,\beta )1(Y\in B)\delta {(1+Y_2)}^3,\hspace{1em}|{\epsilon }_2(Y)|\le 1(Y\notin B)C(\beta ){\delta }^{-3}{(1+Y_1+Y_2)}^3.$$(14)

We point out that

$${\delta }^{-1}\le C(Y_1+Y_2),\hspace{1em} \text{for} \text{any} Y\notin B,$$(15)

which follows from the facts that $Y_1+Y_2\ge \delta (n/2-8)={\delta }^{-1}/2-\delta$ for $Y\notin B$, that $\delta =1/\sqrt{n}$, and that n > 16. Combining (14), (15), and the fact that the moments of Y_i are finite yields

$$\mathbb{E}|{\epsilon }_1(Y)|+\mathbb{E}|{\epsilon }_2(Y)|\le C(b,\beta )\delta \mathbb{E}{(1+Y_1+Y_2)}^7\le C(b,\beta )\delta .$$

Furthermore, applying the Stein factor bounds in Proposition 1 to the bounds on ${\epsilon }_3(Y)$ and ${\epsilon }_4(Y)$ in Proposition 2, and using (15), we get

$$\begin{array}{l}|{\epsilon }_3(Y)+{\epsilon }_4(Y)|\le C(b,\beta )1(Y\in B)\delta {(1+Y_2)}^2+C(b,\beta )1(Y\notin B){\delta }^{-2}{(1+Y_1+Y_2)}^2\hfill \\ \le C(b,\beta )\delta {(1+Y_1+Y_2)}^5.\hfill \end{array}$$

Thus, (13) of Lemma 3 implies that

$$\mathbb{E}({\displaystyle {\int }_0^1|}{\epsilon }_3(Y(s))+{\epsilon }_4(Y(s))|1(Y_1(s)=0)dU(s))\le C(b,\beta )\delta .\hspace{1em}\square$$

3. Stein Factor Bounds

In this section, we prove Proposition 1. We bound the first-order differences in Section 3.1. This requires the most effort. The second-order differences are bounded at the start of Section 3.2, with Section 3.2.1 showing how they can be used to bound $\mathbb{E}|h(X)|$, which may be of independent interest. Section 3.2.2 contains the third-order bounds, and Section 3.2.3 proves two technical lemmas needed for the second-order bounds.

3.1. First-Order Differences

In this section, we bound

by coupling two copies of the JSQ model initialized one customer apart. The coupling is introduced in the following lemma, which is stated in terms of the unscaled CTMC $\{Q(t)\}$.

Lemma 4.

For $1\le i\le b+1$, define ${\mathrm{\Theta }}_i^Q=\{(q,\tilde{q})\in S_Q\times S_Q:q_i<n,\hspace{0.17em}{\tilde{q}}_i=q_i+1\}$. There exists a coupling $\{\tilde{Q}(t)\}$ of $\{Q(t)\}$ whose transient distribution satisfies

$${\{\tilde{Q}(t)|(Q(0),\tilde{Q}(0))\in {\mathrm{\Theta }}_i^Q,\hspace{0.17em}Q(0)=q\}}_{t\ge 0}\stackrel{d}{=}\{Q(t)|Q(0)=(q+e^{(i)})\}.$$(16)

Furthermore, if $(Q(0),\tilde{Q}(0))\in {\displaystyle {\cup }_{i=1}^{b+1}{\mathrm{\Theta }}_i^Q}$, then

1. The equality $\tilde{Q}(t)=Q(t)\mathit{\text{holds for all times}}t\ge {\tau }_C$, where ${\tau }_C=\inf \{t\ge 0:Q(t)=\tilde{Q}(t)\}$.

2. The pair $(Q(t),\tilde{Q}(t))$ belongs to ${\displaystyle {\cup }_{i=1}^{b+1}{\mathrm{\Theta }}_i^Q}$ for all times $t<{\tau }_C$.

3. Let V be a unit-mean exponentially distributed random variable independent of $\{Q(t)\}$. Then

   $${\tau }_C\stackrel{d}{=}\min \left\{\underset{t\ge 0}{\inf }\right\{{\displaystyle {\int }_0^{t}1}((Q(s),\tilde{Q}(s))\in {\mathrm{\Theta }}_1^Q)ds=V\},\hspace{0.17em}\underset{t\ge 0}{\inf }\{Q_{b+1}(t)=n\}\}.$$(17)

Proof of Lemma 4.

Let us construct a joint CTMC $\{(Q(t),\tilde{Q}(t))\}$ by specifying its transitions. For simplicity, we refer to $\{Q(t)\}$ as system 1 and to $\{\tilde{Q}(t)\}$ as system 2. We think of system 2 as a copy of system 1 but with an additional low-priority customer following a preemptive resume rule. That is, service is interrupted, and the extra customer moves to the back of its buffer when a regular customer joins, even if the low-priority customer is currently in service.

Any state in ${\mathrm{\Theta }}_1^Q$ is one where the low-priority customer is in service. The remaining ${\mathrm{\Theta }}_i^Q$ correspond to states where the low-priority customer is assigned to a server with a total of i customers; Figure 3 contains an example of a states in ${\mathrm{\Theta }}_1^Q$ and ${\mathrm{\Theta }}_3^Q$. Assuming $(Q(0),\tilde{Q}(0))=(q,\tilde{q})\in {\mathrm{\Theta }}_i^Q$ for some $1\le i\le b+1$, we now describe the possible transitions of the joint chain.

If i = 1, then the low-priority customer is in service. After a unit-mean exponentially distributed amount of time, the customer leaves system 2 and both systems couple. After coupling, systems 1 and 2 are identical in terms of current and future customers, so they coincide on every sample path. All other transitions of the joint chain are based on the standard transitions of the JSQ model. In other words, a service completion by any of the q₁ servers working in system 1 results in a customer departure from both systems.

Figure 4 illustrates the effect of arrivals when $(q,\tilde{q})\in {\mathrm{\Theta }}_1^Q$. Namely, when $q_1\le n-2$, a new arrival is assigned to the same idle server in both systems. If a customer arrives when $q_1=n-1$, then system 1 has only one idle server and system 2 has none. In system 1, that customer will be assigned to the last remaining idle server. Recall that when defining our JSQ model, we allowed for an arbitrary tie-breaking decision in routing arrivals. Therefore, in system 2, we assign that customer to the server working on the low-priority customer, causing a service preemption and pushing the low-priority customer to the back of the buffer. An arrival when $q_1=n-1$ transitions the joint chain from ${\mathrm{\Theta }}_1^Q$ to ${\mathrm{\Theta }}_2^Q$.

If $2\le i\le b$, then the low-priority customer is in the back of some server’s buffer. A service completion by any of the q₁ servers working in system 1 results in a customer departure from both systems. If, however, the service completion happens at the server containing the low-priority customer, then the chain transitions from ${\mathrm{\Theta }}_i^Q$ to ${\mathrm{\Theta }}_{i-1}^Q$ because the low-priority customer is now assigned to a server with i – 1 customers; see Figure 5 for a depiction of such a transition. All new arrivals get assigned to the same server in each system. Note that if an arrival happens when $q_i=n-1$ and $q_1=\cdots =q_{i-1}=n$, then the system transitions from ${\mathrm{\Theta }}_i^Q$ to ${\mathrm{\Theta }}_{i+1}^Q$.

The final case is when $i=b+1$. All transitions are identical to the $2\le i\le b$ case, except for a customer arrival to a system where $q_{b+1}=n-1$ and $q_1=\cdots =q_b=n$. In that case, system 1 assigns the customer to the last available slot; but system 2 blocks the customer because it is already full. This transition causes the two systems to couple. Note that our construction immediately implies the three claims in Lemma 4. □

Let $\tilde{X}(t)=(\delta (n-{\tilde{Q}}_1(t)),\delta {\tilde{Q}}_2(t),\dots ,\delta {\tilde{Q}}_{b+1}(t))$ be the scaled version of $\tilde{Q}(t)$. For any $x^q\in S$ with $x_1^q>0$, and any $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$,

where ${\mathbb{E}}_{(x,x-\delta e^{(1)})}(·)$ denotes the expectation given $(X(0),\tilde{X}(0))=(x,x-\delta e^{(1)})$. The inequality is true because the gap between $\{X(t)\}$ and $\{\tilde{X}(t)\}$ never increases beyond one customer. The same argument implies that $|{\mathrm{\Delta }}_i{f}_h(x^q)|\le \delta {\mathbb{E}}_{(x,x+\delta e^{(i)})}{\tau }_C$ for $i\ge 2$, and we see that bounding the first-order Stein factors amounts to bounding the expected coupling time τ_C. The following lemma provides the necessary bound. It is worth highlighting that proving this result requires a large amount of effort and JSQ-model-specific insight.

Lemma 5.

For any $(q,\tilde{q})\in {\displaystyle {\cup }_{i=1}^{b+1}{\mathrm{\Theta }}_i^Q}$,

$${\mathbb{E}}_{(q,\tilde{q})}{\tau }_C\le C(b,\beta )(1+\delta q_2).$$

Before proving the lemma, we note that the first-order bounds in Proposition 1 are a consequence of (18) and Lemma 5; that is,

Furthermore, note that for any $x^q\in S$ with $x_3^q=0$,

Recall that $f_h(0)=0$ and that the definition of S implies that $\delta (j_1,j_2,0,\dots ,0)\in S$ for any $0\le j_1\le x_1^q/\delta$ and $0\le j_2\le x_2^q/\delta$. Combining these facts with (19) yields

which proves one of the claims from Proposition 1.

We now describe the main idea and introduce several auxiliary lemmas used to prove Lemma 5. Our discussion communicates the main intuition behind the proof, leaving the technical details to Appendix B. Let $\gamma >0$ be a constant independent of n whose precise value will be specified later and define

Additionally, we define the stopping times

We now describe a sequence of cycles, or attempts, such that in each cycle, the probability of the joint chain coupling is bounded from below by a constant independent of n. Given an initial state $(Q(0),\tilde{Q}(0))=(q,\tilde{q})$ belonging to some ${\mathrm{\Theta }}_i^Q$, we wait until ${\tau }_2({\theta }_2)$, which marks the start of the first cycle. From that point, we wait until $\min ({\tau }_1({\theta }_1),{\tau }_2(2{\theta }_2))$. If ${\tau }_1({\theta }_1)\ge {\tau }_2(2{\theta }_2)$, then we give up trying to couple this cycle and wait until ${\tau }_2({\theta }_2)$ to start a fresh cycle. If ${\tau }_1({\theta }_1)<{\tau }_2(2{\theta }_2)$, then there are $\lfloor \sqrt{n}\beta /2\rfloor$ idle servers and at most $2{\theta }_2$ nonempty buffers. From such a state, we are guaranteed that coupling happens if the joint CTMC enters ${\mathrm{\Theta }}_1^Q$ and spends an exponentially distributed amount of time there before all servers in $\{Q(t)\}$ become busy; that is, ${\tau }_C<{\tau }_1(n)$. If ${\tau }_C\ge {\tau }_1(n)$, we give up trying to couple this cycle and wait until ${\tau }_2({\theta }_2)$ for the next cycle to restart the coupling attempt. Note that this cycle sequence resembles a renewal sequence, but the new cycle times are not renewal times because the values of $Q_3(·),\dots ,Q_{b+1}(·)$ can vary at the start of each new cycle.

From our discussion, it follows that coupling is guaranteed in any given cycle if, starting from a state with $q_2={\theta }_2$, the events $\{{\tau }_1({\theta }_1)<{\tau }_2(2{\theta }_2)\}$ and $\{{\tau }_C<{\tau }_1(n)\}$ occur. In Appendix B, we derive a lower bound, uniform in n, on the probability of coupling in a given cycle, implying that coupling is guaranteed to happen after a geometrically distributed number of cycles. We also derive an upper bound, uniform in n, on the expected time until the start of the first cycle, as well as the expected cycle duration, and then combine these bounds and prove Lemma 5.

3.2. Higher-Order Bounds

To prove the higher-order bounds, we first use the Poisson equation to write ${\mathrm{\Delta }}_1^2{f}_h(x^q)$ in terms of $h(x^q),\hspace{0.17em}\mathbb{E}h(X)$ and first-order differences of $f_h(x^q)$. With the help of this expression, we use the dynamics of the JSQ model to relate all the second-order differences to each other and prove that

for ${\Vert a\Vert }_1=2,\hspace{0.17em}{x}_1^q\le \delta (n-a_1)$ and $x_2^q\le \delta (n-a_2)$, followed by a similar bound for $|({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2)f_h(0,x_2^q,0,\dots ,0)|$. We then bound ${\sum }_{i=1}^{b+1}\mathbb{E}{X}_i$ using the Poisson equation in Section 3.2.1 and bound $|{\mathrm{\Delta }}_1^3{f}_h(x^q)|$ and $|({\mathrm{\Delta }}_1^2-({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2))f_h(x^q)|$ in Section 3.2.2. In Section 3.2.3, we prove two technical lemmas needed to establish (21). We also briefly discuss (see Remark 1) the advantage of using the prelimit generator approach and working with finite differences of $f_h(x^q)$ as opposed to using the classical generator approach and working with the derivatives of the solution to the Poisson equation for the diffusion.

For the following discussion, we assume that $x^q\in S$ with $x_3^q=0$. Recall from (5) that

We rearrange the Poisson equation $G_X{f}_h(x^q)=\mathbb{E}h(X)-h(x^q)$ to see that when $0<q_1<n$, or alternatively $0<x_1^q<\delta n$,

Note that $\mathbb{E}|h(X)|\le C\mathbb{E}(X_1+\cdots +X_{b+1})$ because $h(0)=0$ and $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$. Together with the bound on ${\mathrm{\Delta }}_i{f}_h(x^q)$ from (19), this implies that

Similarly, if $x_1^q=0$,

and therefore

Not all second-order differences can be bounded like this. For example, the equation for ${\mathrm{\Delta }}_2^2{f}_h(x^q)$ would involve the third-order difference ${\mathrm{\Delta }}_2{\mathrm{\Delta }}_1^2{f}_h(x^q)$, which we have not bounded. Instead, the following lemma relates the remaining second-order differences to ${\mathrm{\Delta }}_1^2{f}_h(x^q)$ and $({\mathrm{\Delta }}_2+{\mathrm{\Delta }}_1)f_h(0,x_2^q,0,\dots ,0)$ using the structure of the JSQ system. The proof is postponed to Section 3.2.3.

Lemma 6.

Fix $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$. Then for any $x^q\in S$ with $x_3^q=0$,

$$\begin{array}{l}|{\mathrm{\Delta }}_1^2{f}_h(x^q)|\le C{\delta }^2+\underset{0\le y_2^q\le x_2^q}{\max }|{\mathrm{\Delta }}_1^2{f}_h(0,y_2^q,0,\dots ,0)|,\hspace{1em} \mathit{provided} \hspace{1em}{x}^q+2\delta e^{(1)}\in S,\hfill \\ |{\mathrm{\Delta }}_2{\mathrm{\Delta }}_1{f}_h(x^q)|\le C{\delta }^2+\underset{\begin{array}{c}0\le y_2^q\le x_2^q\\ j=1,2\end{array}}{\max }|{\mathrm{\Delta }}_j^2{f}_h(0,y_2^q,0,\dots ,0)|,\hspace{1em} \mathit{provided} \hspace{1em}{x}^q+\delta e^{(1)}+\delta e^{(2)}\in S,\hfill \\ |{\mathrm{\Delta }}_2^2{f}_h(x^q)|\le C{\delta }^2+\underset{\begin{array}{c}0\le y_2^q\le x_2^q\\ j=1,2\end{array}}{\max }|{\mathrm{\Delta }}_j^2{f}_h(0,y_2^q,0,\dots ,0)|,\hspace{1em} \mathit{provided} \hspace{1em}{x}^q+2\delta e^{(2)}\in S.\hfill \end{array}$$

We see from Lemma 6 that to bound the second-order differences, we only need bounds on $|{\mathrm{\Delta }}_1^2{f}_h(0,x_2^q,0,\dots ,0)|$ and $|{\mathrm{\Delta }}_2^2{f}_h(0,x_2^q,0,\dots ,0)|$. The former is bounded in (24); for the latter term, we note that for any $x^q\in S$ with $x_1^q=x_3^q=0$,

where the inequality follows from (26). The following lemma bounds the last term on the right-hand side, implying that $|{\mathrm{\Delta }}_2^2{f}_h(0,x_2^q,0,\dots ,0)|\le {\delta }^{2}C{\sum }_{i=1}^{b+1}\mathbb{E}{X}_i+C(b,\beta ){\delta }^2{(1+x_2^q)}^2$ and, consequently, (21). It is proved in Section 3.2.3.

Lemma 7.

For all $n\ge 1$,

$$|f_h(0,x_2^q,0,\dots ,0)-f_h(\delta ,x_2^q+\delta ,0,\dots ,0)|\le C(b,\beta ){\delta }^2(1+x_2^q),\hspace{1em}0\le x_2^q<\delta n.$$(28)

3.2.1. Bounding ${\sum }_{i=1}^{b+1}\mathbb{E}{X}_i$.

The bounds in (21) and (26) do not yet look like the stated bounds in Proposition 1 because the term ${\sum }_{i=1}^{b+1}\mathbb{E}{X}_i$ is present. However, we can bound this expectation using the Poisson equation as follows. Recall that $\lambda =1-\beta /\sqrt{n}$; let $x(\infty )=(\delta (n-\lfloor n\lambda \rfloor ),0,\dots ,0)=(\beta +\delta (n\lambda -\lfloor n\lambda \rfloor ),0,\dots ,0)$, and observe that this point is in S. In fact, it is the closest point in S, when rounded up, to the fluid equilibrium of the JSQ system, which happens to be $(\beta ,0,\dots ,0)$; see Braverman (2020). From (22) we have

Choosing $h(x^q)={\sum }_{i=1}^{b+1}{x}_i^q$ and noting that $h(x(\infty ))=\beta +\delta (n\lambda -\lfloor n\lambda \rfloor )$ yields

To bound ${\sum }_{i=1}^{b+1}\mathbb{E}{X}_i$ we need only bound ${\mathrm{\Delta }}_1^2{f}_h(x(\infty )-\delta e^{(1)})$ because $|{\mathrm{\Delta }}_1{f}_h(x(\infty ))|\le \delta C(b,\beta )$ because of (19). Note that we cannot use (24) for the second-order difference bound because ${\sum }_{i=1}^{b+1}\mathbb{E}{X}_i$ is present on the right-hand side there. Instead, we exploit the structure of the JSQ model to bound ${\mathrm{\Delta }}_1^{2}f(x(\infty )-\delta e^{(1)})$ as follows.

Define ${\tau }^{-}(x_1^q)={\inf }_{t\ge 0}\{X(t)=(x_1^q-\delta ,0,\dots ,0)|X(0)=(x_1^q,0,\dots ,0)\}$; let $(X(t),\tilde{X}(t))$ be the scaled version of the coupling defined in Lemma 4, and let V be the unit-rate exponentially distributed random variable defined in the same lemma. Fix $x^q=(x_1^q,0,\dots ,0)$ with $x_1^q\ge 2\delta$, and suppose $X(0)=x^q$ and $\tilde{X}(0)=x^q-\delta e^{(1)}$. Consider the evolution of $(X(t),\tilde{X}(t))$ for $t\in [0,V\wedge {\tau }^{-}(x_1^q)]$. If $V<{\tau }^{-}(x_1^q)$, the two processes couple and become identical. Otherwise, the joint process is in state $(x^q-\delta e^{(1)},x^q-2\delta e^{(1)})$. Using the strong Markov property, we conclude that

Choosing $x_1^q=x_1(\infty )$, we see that

Choosing $h(x^q)={\sum }_{i=1}^{b+1}{x}_i^q$ and using $|{\mathrm{\Delta }}_{1}f(x(\infty ))|\le \delta C(b,\beta )$, we arrive at

The quantities involving ${\tau }^{-}(x_1(\infty ))$ are bounded in the following lemma.

Lemma 8.

There exists a constant $C(\beta )>0$ such that for all $n\ge 1$,

$$\mathbb{E}{\tau }^{-}(x_1^q)\le C(\beta )\delta ,\hspace{1em} \text{and} \hspace{1em}\mathrm{\mathbb{P}}(V\le {\tau }^{-}(x_1^q))\le C(\beta )\delta ,\hspace{1em} \text{for} \hspace{1em}{x}_1^q\in \{x_1(\infty ),\delta ,2\delta \}.$$(31)

Lemma 8 is proved in Section B.5 in Appendix B. It implies that $|{\mathrm{\Delta }}_1^2{f}_h(x(\infty )-\delta e^{(1)})|\le C(b,\beta ){\delta }^2$, and therefore

Combining (32) with (21) proves the second-order bounds in Proposition 1.

Before moving on, let us make a few remarks. The bound in (32) implies that the sequence of steady-state distributions ${\{X\}}_{n=1}^{\infty }$ is tight; when combined with process-level convergence of $\{X(t)\}$ to the diffusion $\{Y(t)\}$, tightness can be used to imply convergence of the steady-state distributions via a limit-interchange argument. For an example of this applied to the JSQ model, see Braverman (2020). Alternatively, (32) can be recast into a result about the convergence rate to the mean-field equilibrium.

Let $h(x)=|x_1+\cdots +x_{b+1}-\beta |-\beta$, noting that $h\in {\mathcal{M}}_{\mathit{disc},1}(1)$ and that $h(0)=0$; suppose for the sake of exposition that $\lfloor n\lambda \rfloor =n\lambda$. One may check that the bound in (30) holds even when $h\in {\mathcal{M}}_{\mathit{disc},1}(1)$, in which case (29) implies that

If we divide both sides by $\sqrt{n}$ to consider the mean-field scaled version of ${\sum }_{i=1}^{b+1}{X}_i$, we get

Thus, we recover the $1/\sqrt{n}$ rate of convergence to the mean field equilibrium that one typically obtains using Stein’s method for the mean-field model, like in Ying (2017). The approach used to show tightness in this section can offer an alternative to the one proposed by Ying (2017), but the difficulty of implementing our approach is directly related to the difficulty of obtaining the relevant Stein factor bounds.

As a final remark, in this section we have shown that establishing tightness, or rates of convergence to the mean-field equilibrium, is equivalent to bounding the first- and second-order differences of $f_h(x^q)$ at a single point near the fluid equilibrium of the CTMC. In contrast, establishing rates of convergence to the diffusion requires bounds on the second- and third-order differences at all points in the support of Y.

3.2.2 Third-Order Bounds.

To bound ${\mathrm{\Delta }}_1^3{f}_h(x^q)$, we recall (23), which says that for $0<x_1^q<\delta n$ with $x_3=0$,

Applying ${\mathrm{\Delta }}_1$ to both sides yields

The bounds on the first- and second-order differences of $f_h(x^q)$, together with the fact that $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$, imply that

which matches the inequality in Proposition 1. The bound on $|({\mathrm{\Delta }}_1^2-({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2))f_h(x^q)|$ when $x_1^q=0$ is proved identically by subtracting $({\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2)f_h(x^q)$ in (25) from ${\mathrm{\Delta }}_1^2{f}_h(x^q)$ in (23). This concludes the proof of Proposition 1. □

3.2.3. Proving Lemmas 6 and 7.

To conclude the section, we prove the auxiliary lemmas from Section 3.2.

Proof of Lemma 6.

Our first task is to bound

$${\mathrm{\Delta }}_1^2{f}_h(x^q)={\displaystyle {\int }_0^{\infty }(}{\mathbb{E}}_{x^q+2\delta e^{(1)}}h(X(t))-2{\mathbb{E}}_{x^q+\delta e^{(1)}}h(X(t))+{\mathbb{E}}_{x^q}h(X(t)))dt.$$

Note that $x^q\in S$ with $x^q+2\delta e^{(1)}\in S$ implies that $q_2\le q_1-2$. Working with the unscaled CTMC, we now construct four processes $\{{\tilde{Q}}^{(1)}(t)\},\dots ,\{{\tilde{Q}}^{(4)}(t)\}$ defined on the time interval $[0,{\tau }_1(n)]$, where

$${\tau }_1(n)=\underset{t\ge 0}{\inf }\{{\tilde{Q}}_1^{(1)}(t)=n\}=\underset{t\ge 0}{\inf }\{{\tilde{X}}_1^{(1)}(t)=0\}.$$(33)

We refer to $\{{\tilde{Q}}^{(i)}(t)\}$ as the ith process. Process four is a copy of $\{Q(t)\}$. Numbers two and three are copies of four but with one extra customer who is assigned to a server with an empty buffer. The extra customer in two is different from the one in three. Lastly, process one is a copy of four but with two extra customers. The extra customers are the same as those in two and three. Figure 6 visualizes the initial condition of the processes.

Let $\{{\tilde{X}}^{(1)}(t)\},\dots ,\{{\tilde{X}}^{(4)}(t)\}$ be the scaled counterparts of these processes. Note that

$$\begin{array}{l}{\mathrm{\Delta }}_1^2{f}_h(x^q)={\displaystyle {\int }_0^{\infty }(}{\mathbb{E}}_{x^q+2\delta e^{(1)}}h(X(t))-2{\mathbb{E}}_{x^q+\delta e^{(1)}}h(X(t))+{\mathbb{E}}_{x^q}h(X(t)))dt\hfill \\ \hspace{1em}\hspace{1em}={\displaystyle {\int }_0^{\infty }{\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x^q}}((h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t))))dt.\hfill \end{array}$$

We refer to the different customers according to their shapes in Figure 6. Define τ_s and τ_d to be the service times of the server with the star and diamond customer, respectively. Both are exponentially distributed with unit mean. Setting ${\tau }_m=\min \{{\tau }_s,{\tau }_d,{\tau }_1(n)\}$, we observe that if ${\tau }_m={\tau }_s$, then

$${\tilde{X}}^{(1)}(t)={\tilde{X}}^{(3)}(t),\hspace{1em}{\tilde{X}}^{(2)}(t)={\tilde{X}}^{(4)}(t),\hspace{1em}t\ge {\tau }_m;$$

if ${\tau }_m={\tau }_d$, then

$${\tilde{X}}^{(1)}(t)={\tilde{X}}^{(2)}(t),\hspace{1em} \text{and} \hspace{1em}{\tilde{X}}^{(3)}(t)={\tilde{X}}^{(4)}(t),\hspace{1em}t\ge {\tau }_m.$$

Therefore,

$$\begin{array}{ll}\hfill & {\displaystyle {\int }_0^{\infty }{\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x}}((h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t))))dt\hfill \\ =\hfill & {\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x}{\displaystyle {\int }_0^{{\tau }_m}(}(h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t))))dt\hfill \\ \hfill & +{\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=x}({\tau }_m={\tau }_1(n)){\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x}[{\mathrm{\Delta }}_1^2{f}_h(0,{\tilde{X}}_2^{(1)}({\tau }_1(n)),0,\dots ,0)|{\tau }_m={\tau }_1(n)].\hfill \end{array}$$(34)

Because ${\tilde{X}}^{(4)}(t)={\tilde{X}}^{(3)}(t)+\delta e^{(1)}={\tilde{X}}^{(2)}(t)+\delta e^{(1)}={\tilde{X}}^{(1)}(t)+2\delta e^{(1)}$ for $0\le t\le {\tau }_m$,

$$|(h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t)))|=|{\mathrm{\Delta }}_1^{2}h({\tilde{X}}^{(1)}(t))|\le C{\delta }^2,$$

where the last inequality follows from $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$. Combining this with the facts that ${\tilde{X}}_2^{(1)}({\tau }_1(n))\le {\tilde{X}}_2^{(1)}(0)$ and ${\mathbb{E}}_x{\tau }_m\le \mathbb{E}{\tau }_s=1$, we conclude that the right-hand side of (34) is bounded by $C{\delta }^2+|{\mathrm{\Delta }}_1^2{f}_h(0,x_2^q,0,\dots ,0)|$, which proves the bound on $|{\mathrm{\Delta }}_1^2{f}_h(x^q)|$.

The remaining bounds are proved similarly, starting with $|{\mathrm{\Delta }}_2{\mathrm{\Delta }}_1{f}_h(x^q)|$. Fix $x^q\in S$ with $x_3^q=0$, and consider

$${\mathrm{\Delta }}_2{\mathrm{\Delta }}_1{f}_h(x^q)=(f_h(x^q+\delta e^{(1)}+\delta e^{(2)})-f_h(x^q+\delta e^{(1)}))-(f_h(x^q+\delta e^{(2)})-f_h(x^q)).$$

We again construct a coupling $\{{\tilde{X}}^{(1)}(t)\},\dots ,\{{\tilde{X}}^{(4)}(t)\}$ corresponding to the four initial states on the right-hand side above. The initial conditions of the unscaled processes are visualized in Figure 7. Our construction yields

$${\mathrm{\Delta }}_2{\mathrm{\Delta }}_1{f}_h(x^q)={\displaystyle {\int }_0^{\infty }{\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x^q}}((h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t))))dt.$$(35)

Let ${\nu }_1={\inf }_{t\ge 0}\{{\tilde{Q}}_1^{(3)}(t)=n\}$. We again let τ_s and τ_d be the remaining service time of the server with the star and diamond customer, respectively, and set ${\tau }_m=\min \{{\tau }_s,{\tau }_d,{\nu }_1\}$. Just like we argued before, if ${\tau }_m={\tau }_s$, then the integrand in (35) is zero after τ_m. If, however, ${\tau }_m={\tau }_d$, then

$$\begin{array}{ll}{\tilde{X}}_2^{(1)}({\tau }_m)=\hfill & {\tilde{X}}_2^{(2)}({\tau }_m)={\tilde{X}}_2^{(3)}({\tau }_m)={\tilde{X}}_2^{(4)}({\tau }_m),\hfill \\ {\tilde{X}}_1^{(2)}({\tau }_m)+2\delta =\hfill & {\tilde{X}}_1^{(1)}({\tau }_m)+\delta ={\tilde{X}}_1^{(4)}({\tau }_m)+\delta ={\tilde{X}}_1^{(3)}({\tau }_m);\hfill \end{array}$$

if ${\tau }_m={\nu }_1$, then ${\tilde{X}}_1^{(i)}({\tau }_m)=0$ for $1\le i\le 4$ and

$${\tilde{X}}_2^{(2)}({\tau }_m)={\tilde{X}}_2^{(1)}({\tau }_m)+\delta ={\tilde{X}}_2^{(4)}({\tau }_m)+\delta ={\tilde{X}}_2^{(3)}({\tau }_m)+2\delta .$$

Therefore,

$$\begin{array}{l}{\mathrm{\Delta }}_2{\mathrm{\Delta }}_1{f}_h(x^q)={\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=x^q}{\displaystyle {\int }_0^{{\tau }_m}(}(h({\tilde{X}}^{(4)}(t))-h({\tilde{X}}^{(3)}(t)))-(h({\tilde{X}}^{(2)}(t))-h({\tilde{X}}^{(1)}(t))))dt\hfill \\ +{\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=x^q}({\tau }_m={\tau }_d){\mathbb{E}}_x[-{\mathrm{\Delta }}_1^2{f}_h({\tilde{X}}^{(2)}({\tau }_d))|{\tau }_m={\tau }_d]\hfill \\ +{\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=x^q}({\tau }_m={\nu }_1){\mathbb{E}}_x[{\mathrm{\Delta }}_2^2{f}_h(0,{\tilde{X}}_2^{(3)}({\nu }_1),0,\dots ,0)|{\tau }_m={\nu }_1]\hfill \\ \le C{\delta }^2+|{\mathrm{\Delta }}_1^2{f}_h(0,x_2^q,0,\dots ,0)|+|{\mathrm{\Delta }}_2^2{f}_h(0,x_2^q,0,\dots ,0)|.\hfill \end{array}$$(36)

Figure 8 illustrates the coupling needed to bound $|{\mathrm{\Delta }}_2^2{f}_h(x^q)|$. The idea of the proof is again to wait until ${\tau }_1(n)$ and analyze what could happen if one of the servers containing the star or diamond customer completes service before ${\tau }_1(n)$. We leave the details to the reader. □

Remark 1.

Let us say a few words on the advantage of using the prelimit generator comparison approach over the classical generator comparison approach. Lemma 6 is proved using a synchronous coupling of four JSQ systems. The four systems are initialized one or two customers apart from one another; because of the discrete state space of the CTMC, all four systems stay one or two customers apart until they couple. Had we used the classical generator comparison approach, we would have needed to carry out a similar analysis by coupling four copies of the diffusion $\{Y(t)\}$. However, unlike the JSQ coupling, the four diffusions would not maintain their initial spacing relative to each other because $\{Y(t)\}$ takes values in a continuous state space. This would further complicate the analysis as we would now need to keep track of the positions of the four diffusions relative to each other.

Proof of Lemma 7.

We want to bound

$$|f_h(0,x_2^q,0,\dots ,0)-f_h(\delta ,x_2^q+\delta ,0,\dots ,0)|=|{\displaystyle {\int }_0^{\infty }(}{\mathbb{E}}_{(0,x_2^q,0,\dots ,0)}h(X(t))-{\mathbb{E}}_{(\delta ,x_2^q+\delta ,0,\dots ,0)}h(X(t)))dt|.$$

As we are accustomed to doing by now, let us construct a coupling $\{{\tilde{Q}}^{(1)}(t),{\tilde{Q}}^{(2)}(t)\}$ with

$${\tilde{Q}}^{(1)}(0)=(n,q_2,0,\dots ,0),\hspace{1em} \text{and} \hspace{1em}{\tilde{Q}}^{(2)}(0)=(n-1,q_2+1,0,\dots ,0).$$

System two has one less idle server and one more customer waiting in a buffer compared to system one, but the total initial customer count is identical across both systems. The initial condition of both systems is visualized in Figure 9. We assume that the diamond and star customers are independent of each other, that the systems see identical arrivals, and that the rest of the customers are identical across both systems.

Now define τ_d and τ_s to be the remaining service times of the server that has the diamond and star customer, respectively; let ${\nu }_1={\inf }_{t\ge 0}\{{\tilde{Q}}_1^{(2)}(t)=n\}$; set ${\tau }_m=\min \{{\tau }_s,{\tau }_d,{\nu }_1\}$. If ${\tau }_m={\nu }_1$ or ${\tau }_m={\tau }_d$, then ${\tilde{Q}}^{(1)}(t)\stackrel{d}{=}{\tilde{Q}}^{(2)}(t)$ for $t\ge {\tau }_m$. Letting $\{{\tilde{X}}^{(i)}(t)\}$ be the scaled version of $\{{\tilde{Q}}^{(i)}(t)\}$, it follows that

$$\begin{array}{ll}\hfill & f_h(0,x_2^q,0,\dots ,0)-f_h(\delta ,x_2^q+\delta ,0,\dots ,0)\hfill \\ =\hfill & {\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=(0,x_2^q,0,\dots ,0)}{\displaystyle {\int }_0^{{\tau }_m}(}h({\tilde{X}}^{(1)}(t))-h({\tilde{X}}^{(2)}(t)))dt\hfill \\ \hfill & +{\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=(0,x_2^q,0,\dots ,0)}({\tau }_m={\tau }_s){\mathbb{E}}_{x^q}[-{\mathrm{\Delta }}_2{f}_h({\tilde{X}}^{(1)}({\tau }_s))|{\tau }_m={\tau }_s].\hfill \end{array}$$

To bound the first term on the right-hand side, note that

$$|{\mathbb{E}}_{{\tilde{X}}^{(1)}(0)=(0,x_2^q,0,\dots ,0)}{\displaystyle {\int }_0^{{\tau }_m}(}h({\tilde{X}}^{(1)}(t))-h({\tilde{X}}^{(2)}(t)))dt|\le C\delta {\mathbb{E}}_{{\tilde{X}}^{(2)}(0)=(\delta ,x_2^q+\delta ,0,\dots ,0)}{\nu }_1\le C(\beta ){\delta }^2.$$

The first inequality is true because $h\in {\mathcal{M}}_{\mathit{disc},2}(C)$, and the last inequality follows from Lemma 8 with $x_1^q=\delta$ there. Furthermore,

$$\begin{array}{ll}\hfill & {\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=(0,x_2^q,0,\dots ,0)}({\tau }_m={\tau }_s)|{\mathbb{E}}_x[-{\mathrm{\Delta }}_2{f}_h({\tilde{X}}^{(1)}({\tau }_s))|{\tau }_m={\tau }_s]|\hfill \\ \hfill & \le {\mathrm{\mathbb{P}}}_{{\tilde{X}}^{(1)}(0)=(0,x_2^q,0,\dots ,0)}({\tau }_s<{\nu }_1)C(b,\beta )\delta (1+x_2^q)\le C(b,\beta ){\delta }^2(1+x_2^q).\hfill \end{array}$$

The first inequality follows from the bound on the first-order difference in (19) together with the fact that ${\tilde{X}}_2^{(1)}(t)\le x_2^q$ for all $t\in [0,{\tau }_m]$. The second inequality follows by noting that τ_s is independent of ν₁ and using Lemma 8 with $x_1^q=\delta ,\hspace{0.17em}{\tau }^{-}(x_1^q)={\nu }_1$, and $V={\tau }_s$ there. □

4. Conclusion

As stated in the introduction, the Stein factor bounds require the bulk of our efforts. Proving the first-order bounds in Section 3.1 amounts to considering two coupled JSQ systems, initialized with a difference of one customer, and bounding the expected coupling time of this joint chain. We bound the coupling time by considering a sequence of coupling attempts where the probability of coupling in a single attempt is bounded away from zero uniformly in n; the expected interattempt times are also bounded from above, uniformly in n. The coupling time can then be bounded by a sum of a geometrically distributed number of random variables representing the interattempt durations. This renewal-like argument applies more generally to settings where (a) there is a region of the state space where the joint chain is guaranteed to couple provided it spends enough time there and (b) one can control the expected time to reach this region and the probability of coupling in the region before leaving it.

With the first-order Stein factor bounds in hand, the higher-order bounds require less effort. Our proofs of the high-order bounds make heavy use of the transition structure of the JSQ system and, in particular, that $Q_2(t),\dots ,Q_{b+1}(t)$ increase only at those times when $Q_1(t)=n$. Readers should not be misled into thinking that high-order Stein factor bounds require less effort than first-order bounds for all models. Indeed, in the classical generator comparison approach, high-order bounds require much more effort; see, for example, Mackey and Gorham (2016), Erdogdu et al. (2018), Jin et al. (2022).

Regarding extending our results, we note that Proposition 2, which compares G_X to G_Y, can be easily adjusted to hold for other parameter regimes and load-balancing policies. The main difficulty would be establishing Stein factor bounds. As mentioned in the introduction, Zhao et al. (2021) considered the super-Halfin-Whitt regime ($1/2<\alpha <1$) and established several hitting-time estimates similar to the ones we use in the proof of Lemma 5 to bound the first-order Stein factors. It may be possible to build on their results and obtain rates of convergence for the super-Halfin-Whitt regime too.

Furthermore, it seems that the sub-Halfin-Whitt regime ($0<\alpha <1/2$) should present less of a challenge than our own setting. Recall from the discussion in Section 3.1 that coupling of the joint CTMC is guaranteed provided it enters ${\mathrm{\Theta }}_1^Q$ and spends an exponentially distributed amount of time there before all servers become busy. Compared to the Halfin-Whitt regime, the rate at which customers arrive in the sub-Halfin-Whitt regime is much smaller, so the event that all servers are busy should happen less frequently. Indeed, Liu and Ying (2020) showed that the steady-state probability that all servers are busy tends to zero in the sub-Halfin-Whitt regime. Consequently, the Stein factor bounds should be simpler to establish.