2.1. System model

In this section we give details of the system model. We study a large-scale multi-server system with N Erlang loss servers and one central job dispatcher, which routes an incoming request to one of the servers according to a predefined load-balancing policy. A server accepts an incoming request if the occupancy or the number of progressing jobs of the server is less than C, a predefined integer value referred to as the capacity of a server. Otherwise, the request is blocked from service and it is considered to be discarded from the system. Furthermore, an accepted job is processed at a constant unit rate upon its acceptance for service until the job is completed. We assume that the service time distributions are exponential with unit mean. The job dispatcher uses the JSQ(d) load-balancing policy to dispatch the incoming jobs.

The arrival process of jobs is a Poisson process with rate $N\lambda^{(N)}$ , where the parameter $\lambda^{(N)}\in\mathbb{R}_+$ is defined such that, for $\sigma\in\mathbb{R}_+$ and $\beta\in\mathbb{R}$ ,

(3) \begin{equation}\lambda^{(N)}=\sigma-\frac{\beta}{\sqrt{N}\,}.\end{equation}

Clearly, $\lambda^{(N)}$ is a perturbation of $\sigma$ by $\frac{\beta}{\sqrt{N}\,}$ . For our model, we first show that there exists a functional law of large numbers limit referred to as the mean-field limit. Next, we establish an FCLT result, which is exploited to quantify the error between the actual blocking probability of the system with N servers and the asymptotic blocking probability expressed as a function of the unique fixed point of the mean field. The particular form of $\lambda^{(N)}$ in (3) subsumes an important special case of $\sigma=C$ that corresponds to the Halfin–Whitt heavy traffic regime as the resulting traffic intensity $\rho^{(N)}=\frac{\lambda^{(N)}}{C}$ approaches 1 as $N\to\infty$ and $\lim_{N\to\infty}\sqrt{N}(1-\rho^{(N)})=\frac{\beta}{C}$ . In this case, the arrival rate of jobs, $N\lambda^{(N)}$ , and the system capacity, NC, are related as $NC=N\lambda^{(N)}+\beta\sqrt{N}$ , and they converge to $\infty$ as $N\to\infty$ .

2.2. Notation

Since the job dispatcher uses only the knowledge of the server occupancies and does not use their identities to dispatch an arrival, we consider the Markov process $(\mathbf{X}^{(N)}(t), t\geq 0)$ to model the time evolution of the system, where $\mathbf{X}^{(N)}(t)=(\mathbf{X}^{(N)}_i(t),0\leq i\leq C)$ , with $\mathbf{X}^{(N)}_i(t)$ denoting the fraction of the servers with at least i progressing jobs at time t. Let $\mathbb{U}$ be the space defined as

\begin{equation*}\mathbb{U}\triangleq\{(u_0,u_1,\ldots,u_C)\,:\,u_0=1\geq u_1\geq\cdots\geq u_C\geq 0\}.\end{equation*}

It is evident that $\mathbf{X}^{(N)}(t)$ lies in the space $\mathbb{U}$ . Without loss of generality, we write an element of the form $(u_0,\ldots,u_C)$ as $\mathbf{u}$ . The space $\mathbb{U}$ is equipped with the metric generated by the Euclidean norm ${\|{\cdot}\|}_2$ , defined as

(4) \begin{equation}{\|{\mathbf{u}\|}}_2=\sqrt{\sum_{i=0}^C{\left \vert {u_i} \right \vert}^2},\end{equation}

where $\mathbf{u}=(u_0,\ldots,u_C)$ . It can be verified that the space $\mathbb{U}$ is a Polish space.

We study stochastic processes that are defined on $(\Omega,\mathbb{F},\mathbb{P})$ with sample paths belonging to the space of right continuous functions with left limits (such functions are also called càdlàg functions). The space of càdlàg functions is equipped with the Skorohod $J_1$ topology. We write $\mathbf{Y}_n\Rightarrow \mathbf{Y}$ as $n\to\infty$ to indicate that a sequence of random elements $\{\mathbf{Y}_n\}_{n\geq 1}$ converges in distribution to a random element $\mathbf{Y}$ . For two real-valued local martingales $(\mathbf{M}^1(t),t\geq 0)$ and $(\mathbf{M}^2(t),t\geq 0)$ , let the covariation process be denoted by $(\langle\mathbf{M}^1,\mathbf{M}^2\rangle_t,t\geq 0)$ and the quadratic variation process be denoted by $(\langle\mathbf{M}^1\rangle_t,t\geq 0)=(\langle\mathbf{M}^1,\mathbf{M}^1\rangle_t,t\geq0)$ .

2.3 Preliminaries

In this section we present a mathematical model of the main problem and give some preliminary results. We first begin with a discussion on the time evolution of the process $(\mathbf{X}^{(N)}(t), t\geq 0)$ . At an arrival instant t, if the system state is $\mathbf{b}=(b_0,\ldots,b_C)$ , implying that the fraction of the servers with at least i progressing jobs at time t is equal to $b_i$ for $0\leq i\leq C$ , then according to the JSQ(d) policy the destination server of the job will have occupancy n with probability $b_n^d-b_{n+1}^d$ . Since the rate of the arrival process is $N\lambda^{(N)}$ , the total instantaneous rate of arrivals to servers that have n jobs is equal to $N\lambda^{(N)}((\mathbf{X}_{n-1}^{(N)}(t))^d-(\mathbf{X}_{n}^{(N)}(t))^d)$ . Furthermore, since the service times have exponential distributions with unit rate, the total instantaneous departure rate of jobs from servers with n progressing jobs is equal to $nN(\mathbf{X}_{n}^{(N)}(t)-\mathbf{X}_{n+1}^{(N)}(t))$ . As a result, we can model the time evolution of the process $(\mathbf{X}^{(N)}(t), t\geq 0)$ by using random time change of a set of mutually independent unit-rate Poisson processes as in [Reference Pang, Talreja and Whitt17, Section 2.1], which we explain below.

Let $\{(\mathcal{N}_i(t),t\geq 0)\}_{i\geq 1}$ be a set of mutually independent unit-rate Poisson processes where $(\mathcal{N}_i(t),t\geq 0)$ is used to model the arrival process to servers that have $i-1$ progressing jobs. Similarly, let $\{(\mathcal{D}_i(t),t\geq 0)\}_{i\geq 1}$ be the collection of a set of mutually independent unit-rate Poisson processes where $(\mathcal{D}_i(t),t\geq 0)$ is used to model the departure process from servers that have i progressing jobs. Furthermore, the set of processes $\{(\mathcal{D}_i(t),t\geq 0)\}_{i\geq 1}$ is independent of the set of processes $\{(\mathcal{N}_i(t),t\geq 0)\}_{i\geq 1}$ . Also, $\{(\mathcal{N}_i(t),t\geq 0)\}_{i\geq 1}$ and $\{(\mathcal{D}_i(t),t\geq 0)\}_{i\geq 1}$ are independent of $\mathbf{X}^{(N)}(0)$ . Since the arrival process of jobs to the system is a Poisson process with rate $N\lambda^{(N)}$ , and the service time distributions are exponential with unit mean, we can write $\mathbf{X}^{(N)}_0(t)=1$ and, for $n\geq 1$ ,

(5) \begin{align}\mathbf{X}_n^{(N)}(t)= & \mathbf{X}^{(N)}_n(0)+\frac{1}{N}\mathcal{N}_n\bigg(N\lambda^{(N)}\int_{s=0}^t\Big(\Big(\mathbf{X}_{n-1}^{(N)}(s)\Big)^d-\Big(\mathbf{X}_{n}^{(N)}(s)\Big)^d\Big)\,{\text{d}} s\bigg)\\ & -\frac{1}{N}\mathcal{D}_n\bigg(Nn\int_{s=0}^t\Big(\Big(\mathbf{X}_{n}^{(N)}(s)\Big)-\Big(\mathbf{X}_{n+1}^{(N)}(s)\Big)\Big)\,{\text{d}} s\bigg).\nonumber\end{align}

We choose the filtration $(\mathcal{F}^{(N)}(t),t\geq 0)$ where

\begin{align*}\mathcal{F}^{(N)}(t) =& \sigma\left(\mathbf{X}^{(N)}(0), \mathcal{N}_n\bigg(N\lambda^{(N)}\int_{s=0}^r\Big(\Big(\mathbf{X}_{n-1}^{(N)}(s)\Big)^d-\Big(\mathbf{X}_{n}^{(N)}(s)\Big)^d\Big)\,{\text{d}} s\bigg),\right. \\ & \left.\mathcal{D}_n\bigg(Nn\int_{s=0}^r\Big(\Big(\mathbf{X}_{n}^{(N)}(s)\Big)-\Big(\mathbf{X}_{n+1}^{(N)}(s)\Big)\Big)\,{\text{d}} s\bigg),n\geq 1,0\leq r\leq t\right),\nonumber\end{align*}

augmented by all the null sets.

We now present the results on the mean-field analysis of the model without proofs as they directly follow from the case of $\lambda^{(N)}=b$ for $b\in\mathbb{R}_+$ , studied in [Reference Mukhopadhyay, Mazumdar and Guillemin15]. The resulting mean-field equations (MFEs) in our case are the same as for $\lambda^{(N)}=b$ except that $\sigma$ replaces b.

Theorem 1. For $\mathbf{u}\in\mathbb{U}$ , if $\mathbf{X}^{(N)}(0)\Rightarrow \mathbf{u}$ as $N\to\infty$ , then $\big(\mathbf{X}^{(N)}(t),t\geq 0\big)\Rightarrow(\mathbf{x}(t,\mathbf{u}),t\geq 0)$ as $N\to\infty$ , where $(\mathbf{x}(t,\mathbf{u}),t\geq 0)=(x_n(t,\mathbf{u}),t\geq 0,\,0\leq n\leq C)$ is the unique solution to the following: let $\mathbf{h}(\mathbf{x}(t,\mathbf{u}))=(h_n(\mathbf{x}(t,\mathbf{u})),0\leq n\leq C)$ , where

(6) \begin{equation} \mathbf{x}(0,\mathbf{u})=\mathbf{u}, \qquad \frac{{\text{d}} x_n(t,\mathbf{u})}{{\text{d}} t} = h_n(\mathbf{x}(t,\mathbf{u})),\end{equation}

satisfying

(7) \begin{equation}h_0(\mathbf{x}(t,\mathbf{u}))=0 \end{equation}

and, for $n\geq 1$ ,

(8) \begin{equation}h_n(\mathbf{x}(t,\mathbf{u}))=\sigma\big(x_{n-1}^d(t,\mathbf{u})-x_n^d(t,\mathbf{u})\big)-n(x_n(t,\mathbf{u})-x_{n+1}(t,\mathbf{u})) ,\end{equation}

with $x_0(t,\mathbf{u})=1$ and $x_{C+1}(t,\mathbf{u})=0$ . The deterministic process $(\mathbf{x}(t,\mathbf{u}),t\geq 0)$ is referred to as the mean-field limit and equations (6)–(8) are referred to as the mean-field equations with initial point $\mathbf{u}$ .

Without loss of generality, we say that if a process $(\mathbf{y}(t),t\geq 0)$ is a solution to the differential equations (6)–(8) then it is the unique generic solution with initial point $\mathbf{y}(0)$ .

The mean field $(\mathbf{x}(t,\mathbf{u}),t\geq 0)$ has a unique global asymptotically stable fixed point $\boldsymbol{\pi}=(\pi_n,0\leq n\leq C)$ with $\pi_0=1$ . Also, the following exchange of limits holds:

\begin{equation*}\lim_{N\to\infty}\lim_{t\to\infty}\mathbf{X}^{(N)}(t)=\lim_{t\to\infty}\lim_{N\to\infty}\mathbf{X}^{(N)}(t).\end{equation*}

Using this, under the assumption of the exchangeability of initial states of servers, we can show the independence of any finite set of servers as $N\to\infty$ . Also, it can be shown that as $N\to\infty$ a server’s distribution equals $\mathbf{x}(t,\mathbf{u})$ and $\pi$ at time t and in the stationary regime, respectively. As a result, $\pi_C$ denotes the stationary probability that a server is fully occupied as $N\to\infty$ . Since the dispatcher samples d servers upon arrival and selects one of them as the destination server, the stationary average blocking probability of a job as $N\to\infty$ is equal to $\pi_C^d$ , where we use the fact that the chosen d servers are independent of each other. Our objective is to find the gap between the actual blocking probability of the system with N servers and $\pi_C^d$ as a function of the parameters $\beta$ , $\boldsymbol{\pi}$ , $\sigma$ , N, and C.

We can find $\boldsymbol{\pi}$ numerically as follows. The fixed point $\boldsymbol{\pi}$ is the unique solution to

\begin{equation*}\sigma\big(\pi_{n-1}^d-\pi_n^d\big)=n(\pi_n-\pi_{n+1})\end{equation*}

for $n\geq 1$ and $\pi_{C+1}=0$ . Then, from there, we can also write

(9) \begin{equation}\sigma\frac{\big(\pi_{n-1}^d-\pi_n^d\big)}{(\pi_{n-1}-\pi_n)}(\pi_{n-1}-\pi_n)=n(\pi_n-\pi_{n+1})\end{equation}

for $n\geq 1$ and $\pi_{C+1}=0$ . Let us define $\hat{\lambda}_n=\sigma\frac{\big(\pi_n^d-\pi_{n+1}^d\big)}{(\pi_n-\pi_{n+1})}$ . Then, from (8), $\boldsymbol{\pi}$ is the stationary distribution of the single-server loss model with a Poisson arrival process of jobs having rate $\hat{\lambda}_n$ when there are n progressing jobs, and $\pi_n$ is the probability that the server has at least n progressing jobs. Let $\mathbb{M}_1(\{0,1,\ldots,C\})$ be the set of probability measures on $\{0,1,\ldots,C\}$ . Then, from [Reference Mukhopadhyay, Mazumdar and Guillemin15], the fixed point $\boldsymbol{\pi}$ can be computed using the formula for the stationary distribution of a single-server loss system with state-dependent arrival rates. We first define two mappings, $\Theta\,:\,\mathbb{M}_1(\{0,1,\ldots,C\})\mapsto\mathbb{R}_+^{C+1}$ and $\widehat{\Xi}\,:\,\mathbb{R}_+^{C+1}\mapsto \mathbb{M}_1(\{0,1,\ldots,C\})$ , that are used in computing $\boldsymbol{\pi}$ . For every $(p_0,\ldots,p_C)\in \mathbb{M}_1(\{0,1,\ldots,C\})$ , there exists $$({r_0}, \ldots ,{r_C}) \in \mathbb{R}_+^{C+1}$$ such that $\Theta((p_0,\ldots,p_C))=(r_0,\ldots,r_C)$ , where

\begin{equation*}r_n=\sigma\frac{\Big(\Big(\!\sum_{j=n}^Cp_j\Big)^d-\Big(\!\sum_{i=n+1}^Cp_i\Big)^d\Big)}{\Big(\Big(\!\sum_{j=n}^Cp_j\Big)-\Big(\!\sum_{i=n+1}^Cp_i\Big)\Big)}.\end{equation*}

Similarly, for every $$({b_0}, \ldots ,{b_C}) \in \mathbb{R}_+^{C+1}$$ , there exists $(a_0,\ldots,a_C)\in\mathbb{M}_1(\{0,1,\ldots,C\})$ such that $\widehat{\Xi}((b_0,\ldots,b_C))=(a_0,\ldots,a_C)$ , where

\begin{equation*}a_n=\left(\prod_{i=1}^n\bigg(\frac{b_{i-1}}{i}\bigg)\right)a_0\end{equation*}

for $n\geq 1$ and $\sum_{i=0}^Ca_i=1$ . Then, $\boldsymbol{\pi}$ is the unique fixed point of the mapping $\widehat{\Xi}(\Theta)$ , which can be computed numerically.

Our objective is to study the limit of the fluctuation process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ as $N\to\infty$ , where

(10) \begin{equation}\mathbf{Z}^{(N)}(t)=\sqrt{N}(\mathbf{X}^{(N)}(t)-\mathbf{x}(t,\mathbf{u})).\end{equation}

It can be checked that $\mathbf{Z}^{(N)}(t)$ lies in the space $\mathbb{V}$ defined as

\begin{equation*}\mathbb{V}\triangleq\{(r_0,\ldots,r_C)\,:\, r_0=0 \text{ and } r_i\in\mathbb{R}, 1\leq i\leq C \}.\end{equation*}

We equip the space $\mathbb{V}$ with the topology induced by the Euclidean norm (4). Our analysis uses the operator norm ${\|{\cdot}\|}_2$ defined as

\begin{equation*}{\|{K}\|}_2=\sup_{\mathbf{v}\in\mathbb{V}}\frac{{\|{K\mathbf{v}\|}}_2}{{\|{\mathbf{v}\|}}_2},\end{equation*}

where $K\,:\,\mathbb{V}\mapsto\mathbb{V}$ is a linear operator.

Next, we obtain the time evolution of the process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ by using (5), (6), and (10). For this, we first define the following three useful operators $W_1, W_2, W_3\,:\,\mathbb{U}\mapsto\mathbb{V}$ . For $\mathbf{b}\in\mathbb{U}$ with $1\leq i\leq 3$ we define $W_i(\mathbf{b})=((W_i(\mathbf{b}))_n, n\geq 0)$ , where $(W_1(\mathbf{b}))_0=0$ , $(W_2(\mathbf{b}))_0=0$ , $(W_3(\mathbf{b}))_0=0$ , and, for $1\leq n\leq C$ ,

(11) \begin{equation}(W_1(\mathbf{b}))_n=\sigma \big(b_{n-1}^d-b_n^d\big),\,\,(W_2(\mathbf{b}))_n=n (b_{n}-b_{n+1}),\,\,(W_3(\mathbf{b}))_n=\beta \big(b_{n-1}^d-b_n^d\big).\end{equation}

From (8) and (11), we have $h_n(\mathbf{b})=(W_1(\mathbf{b}))_n-(W_2(\mathbf{b}))_n$ . Furthermore, let $W\,:\,\mathbb{U}\mapsto\mathbb{V}$ be the operator defined as $W=W_1-W_2$ . The operator W is Lipschitz continuous, satisfying, for all $\mathbf{a},\mathbf{b}\in\mathbb{U}$ , ${\|{W(\mathbf{a})-W(\mathbf{b})}\|}_2\leq B_W{\|{\mathbf{a}-\mathbf{b}}\|}_2$ , where $B_W=2d\sqrt{\sigma^2+C^2}$ .

Now we define a set of independent square-integrable martingales $(\mathbf{M}^{(N)}(t),t\ge 0)=\{(\mathbf{M}^{(N)}_i(t),t\geq 0)\}_{(i\in{0,1,\ldots,C)}}$ adapted to the filtration $(\mathcal{F}^{(N)}(t),t\geq 0)$ such that $(\mathbf{M}^{(N)}(t),t\ge 0)$ is independent of $\mathbf{Z}^{(N)}(0)$ and, for $i\geq 1$ ,

(12) \begin{equation}\langle\mathbf{M}^{(N)}_i\rangle_t=\int_{s=0}^t\left(\big(W_1\big(\mathbf{X}^{(N)}(s)\big)\big)_i+\big(W_2\big(\mathbf{X}^{(N)}(s)\big)\big)_i-\frac{1}{\sqrt{N}\,}\big(W_3\big(\mathbf{X}^{(N)}(s)\big)\big)_i\right)\,{\text{d}} s.\end{equation}

Then, from (5)–(8) and (10), we get

(13) \begin{align}\mathbf{Z}^{(N)}(t) = &\mathbf{Z}^{(N)}(0)+\int_{s=0}^t\sqrt{N}\big(W\big(\mathbf{X}^{(N)}(s)\big)-W\big(\mathbf{x}(s,\mathbf{u})\big)\big)\,{\text{d}} s\\ & -\int_{s=0}^tW_3\big(\mathbf{X}^{(N)}(s)\big)\,{\text{d}} s+\mathbf{M}^{(N)}(t).\nonumber \end{align}

3.1. Transient regime

In this section we show that the process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ converges to an OU process in the transient regime as $N\to\infty$ . First, we begin with the following result that concludes stochastic boundedness of the process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ when $N\to\infty$ . We use this property in proving the tightness of the sequence $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ .

Lemma 1. For any $T>0$ , if $\limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(0)}\|}_2^2}\right]}<\infty$ , then

(14) \begin{equation} \limsup_{N\to\infty}{\mathbb{E}\left[{\sup_{0\leq t\leq T}{\|{\mathbf{Z}^{(N)}(t)}\|}_2^2}\right]}<\infty. \end{equation}

Proof. See Section 4.1.

As we show later in this section, any limit point of $(\mathbf{Z}^{(N)}(t),t\geq 0)$ is a solution to the stochastic differential equation (SDE) defined in (20). Next, we introduce some notation that is used in (20). The proposed SDE depends on $(\mathbf{s}(t),t\geq 0)$ , which is a solution of (18). Let $(\mathbf{w}(t),t\geq 0)$ be a generic solution of (6)–(8) with an initial point $\mathbf{w}(0)$ , where $\mathbf{w}(t)=(w_n(t), 0\leq n\leq C+1)$ satisfying $w_0(t)=1$ and $w_{C+1}(t)=0$ . Then we have

(15) \begin{align}\frac{{\text{d}} w_n(t)}{{\text{d}} t}&=h_n(\mathbf{w}(t)),\end{align}

where

(16) \begin{equation}h_0(\mathbf{w}(t))=0\end{equation}

and, for $1\leq n\leq C$ ,

(17) \begin{equation}h_n(\mathbf{w}(t))=\sigma\big(w_{n-1}^d(t)-w_n^d(t)\big)-n(w_n(t)-w_{n+1}(t)).\end{equation}

By linearizing (15)–(17) around a solution $(\mathbf{r}(t),t\geq 0)$ of (15)–(17) with an initial point $\mathbf{r}(0)$ , we get

(18) \begin{equation}\frac{{\text{d}}\mathbf{s}(t)}{{\text{d}} t}=H(\mathbf{r}(t))\mathbf{s}(t),\end{equation}

where, for $\mathbf{a}\in\mathbb{U}$ and $\mathbf{b}\in\mathbb{V}$ , the linear operator $H(\mathbf{a})\,:\,\mathbb{V}\mapsto\mathbb{V}$ is defined as

\begin{equation*}(H(\mathbf{a})\mathbf{b})_n=\sigma da_{n-1}^{d-1}b_{n-1}-\big(\sigma d a_n^{d-1}+n\big)b_n+nb_{n+1},\end{equation*}

$n\geq 1$ . Note that any solution $(\mathbf{s}(t), t\geq 0)$ of (18) satisfies $\mathbf{s}(t)=\mathbf{w}(t)-\mathbf{r}(t)$ , where $(\mathbf{w}(t), t\geq0)$ is a solution of (15)–(17) with an initial point $\mathbf{w}(0)$ .

We will show that the limit of the sequence $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ depends on the process $(\mathbf{s}(t),t\geq 0)$ when $(\mathbf{r}(t),t\geq 0)$ in (18) is replaced with the mean field $(\mathbf{x}(t,\mathbf{u}),t\geq 0)$ . The operator $H(\mathbf{a})$ is a matrix in the canonical basis $(0,1,0,\ldots,0)$ , $(0,0,1,0,\ldots,0)$ , $\ldots$ , $(0,0,\ldots,0,1)$ , where the dimension of each vector is $C+1$ . We can write $H(\mathbf{a})$ as the following matrix of size $C\times C$ :

\[H(\mathbf{a})=\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} -\nu_1 & 1 & 0 &0&\cdots&0\\[3pt] \gamma_1 & -\nu_2&2&0&\cdots&0\\[3pt] 0& \gamma_2&-\nu_3&3&\cdots&0\\[3pt] \vdots&\vdots&\ddots&\ddots&\ddots&\vdots\\[3pt] 0& 0&\cdots&\gamma_{C-2}&-\nu_{C-1}&C-1\\[3pt] 0&0&\cdots&0&\gamma_{C-1}&-\nu_C \end{array} } \right],\]

where $\gamma_i=\sigma d a_i^{d-1}$ and $\nu_i=\gamma_i+i$ , $1\leq i\leq C$ .

Let $(\mathbf{M}(t),t\geq 0)=\{(\mathbf{M}_i(t),t\geq 0)\}_{i\in\{0,1,\ldots,C\}}$ be a collection of mutually independent real-valued continuous and centered Gaussian martingales, determined in law by their deterministic quadratic variation process

(19) \begin{equation}\langle\mathbf{M}_n\rangle_t=\int_{s=0}^t((W_1(\mathbf{x}(s,\mathbf{u}))_n+(W_2(\mathbf{x}(s,\mathbf{u}))_n)\,{\text{d}} s \end{equation}

for $n\geq 0$ . Note that both $\mathbf{M}(t)$ and $(\langle\mathbf{M}_i\rangle_t,0\leq i\leq C)$ lie in $\mathbb{V}$ . From (19), the martingale $(\mathbf{M}(t),t\geq 0)$ is square integrable since $((W_1(\mathbf{b}))_n+(W_2(\mathbf{b}))_n)$ for $\mathbf{b}\in\mathbb{U}$ is uniformly bounded in n and $\mathbf{b}$ because $0\leq b_i\leq 1$ for $0\leq i\leq C$ .

Now we introduce an SDE; we later show that the limit of the process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ as $N\to\infty$ in the transient regime is a unique solution of this SDE.

Definition 2. SDE for the transient regime: Let $(\mathbf{Z}(t),t\geq 0)$ be a solution of the SDE

(20) \begin{equation} \mathbf{Z}(t)=\mathbf{Z}(0)+\int_{s=0}^tH(\mathbf{x}(s,\mathbf{u}))\mathbf{Z}(s) \,{\text{d}} s-\int_{s=0}^tW_3(\mathbf{x}(s,\mathbf{u}))\,{\text{d}} s+\mathbf{M}(t). \end{equation}

The solution of (20) is an OU process. Next, we study the SDE (20).

The proof of Theorem 2 follows by arguments similar to the proof of [Reference Vasantam and Mazumdar19, Theorem 2], and hence we omit the proof here.

Now we present the main result on the transient regime.

Proof. See Section 4.2.

Remark 2. For a constant $a\in\mathbb{R}_+$ , if $\lambda^{(N)}=a$ then, for $n\geq 1$ ,

\begin{equation*} (W_1(\mathbf{b}))_n=a (b_{n-1}^d-b_n^d), \qquad (W_2(\mathbf{b}))_n=n (b_{n}-b_{n+1}), \qquad (W_3(\mathbf{b}))_n=0. \end{equation*}

As a result, from the SDE (20), we recover the following SDE obtained in [Reference Vasantam and Mazumdar19] for the $\lambda^{(N)}=a$ case:

\begin{equation*} \mathbf{Z}(t)=\mathbf{Z}(0)+\int_{s=0}^tH(\mathbf{x}(s,\mathbf{u}))\mathbf{Z}(s) \,{\text{d}} s+\mathbf{M}(t). \end{equation*}

3.2. Stationary regime

In this section we present results pertaining to the stationary regime. Here, the mean field is located at $\boldsymbol{\pi}$ , and hence we assume that $\mathbf{u}=\boldsymbol{\pi}$ . We recall that $\boldsymbol{\pi}$ satisfies $W(\boldsymbol{\pi})=W_1(\boldsymbol{\pi})-W_2(\boldsymbol{\pi})=0$ .

Our objective is to show that the sequence of processes $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ converges to a limit in the stationary regime as $N\to\infty$ , where $(\mathbf{Z}^{(N)}(t),t\geq 0)$ is the process defined in (10). First, we will show that the sequence $\{\mathbf{Z}^{(N)}(t)\}_{N\geq 1}$ is relatively compact in the stationary regime by using another process $(\mathbf{Q}^{(N)}(t),t\geq 0)$ , where $\mathbf{Q}^{(N)}(t)=\sqrt{N}(\mathbf{X}^{(N)}(t)-\boldsymbol{\pi})$ . We introduce an SDE, (22), and show that there exists a unique solution to this SDE with a unique invariant law. We then use this result to prove that any limit point of the process $(\mathbf{Z}^{(N)}(t),t\geq 0)$ as $N\to\infty$ in the stationary regime is a stationary OU process with the same invariant law as that of the solution of the proposed SDE, (22). As a result, it would imply that $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ converges as $N\to\infty$ in the stationary regime to the solution of (22).

Next, we state the exponential stability of the mean field in Lemma 2, which we later use in the proof of the result stated in Lemma 3.

Proof. See Section 4.6.

The following result shows the tightness of $\big\{\mathbf{Z}^{(N)}(t)\big\}_{N\geq 1}$ in the stationary regime.

Lemma 3. If $\limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}_2^2}\right]}<\infty$ , then

\begin{equation*} \limsup_{N\to\infty}\sup_{t\geq 0}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t)}\|}_2^2}\right]}<\infty. \end{equation*}

Consequently, in the stationary regime corresponding to $t=\infty$ ,

\begin{equation*} \limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(\infty)}\|}_2^2}\right]}<\infty. \end{equation*}

Proof. See Section 4.3.

Next, we state the SDE that is used to obtain the limit of the sequence $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ in the stationary regime. First, we linearize (8) around $\boldsymbol{\pi}$ to obtain a process $(\mathbf{s}(t), t\geq 0)$ satisfying

(21) \begin{equation}\frac{{\text{d}}\mathbf{s}(t)}{{\text{d}} t}=H(\boldsymbol{\pi})\mathbf{s}(t).\end{equation}

Let $\mathbf{B}(t)=(\mathbf{B}_i(t),0\leq i\leq C)$ with $\mathbf{B}_0(t)=0$ , where $\{(\mathbf{B}_i(t),t\geq 0)\}_{0\leq i\leq C}$ are independent centered Brownian motions and ${\mathbb{E}\left[{\mathbf{B}_i^2(1)}\right]}=V_i=\text{var}(\mathbf{B}_i(1))=2i(\pi_i-\pi_{i+1})$ , $i\geq 1$ . The infinitesimal covariance matrix of $(\mathbf{B}(t),t\geq 0)$ is the diagonal matrix diag( $\mathbf{V})$ , where $\mathbf{V}=(V_n,0\leq n\leq C)$ . From (19), the martingales $(\mathbf{M}(t),t\geq 0)$ with $\mathbf{x}(0,\boldsymbol{\pi})=\boldsymbol{\pi}$ have the same law as $(\mathbf{B}(t),t\geq 0)$ . Now we define the following SDE, which is used to study the process $(\mathbf{Z}(t),t\geq 0)$ in the stationary regime.

Definition 3. An SDE for the stationary regime: Let $(\mathbf{Q}(t),t\geq 0)$ be a solution to the SDE

(22) \begin{equation} \mathbf{Q}(t)=\mathbf{Q}(0)+\int_{s=0}^tH(\boldsymbol{\pi})\mathbf{Q}(s)\,{\text{d}} s-\int_{s=0}^tW_3(\boldsymbol{\pi})\,{\text{d}} s+\mathbf{B}(t). \end{equation}

Then, (22) defines an OU process whose drift and variance depend on $\boldsymbol{\pi}$ .

For an arbitrary $\mathbf{Q}(0)$ in (22), we have the following result. The proof follows by the same arguments as in the proof of Theorem 2, and hence is omitted.

We point out that the transpose $H(\boldsymbol{\pi})^*$ of $H(\boldsymbol{\pi})$ is the generator of a finite-state birth–death process, and the birth, death, and killing rates in state i ( $1\leq i\leq C$ ) are $\gamma_i$ , $i-1$ , and 1, respectively. Let $\mathcal{I}$ be the identity matrix of dimension $C\times C$ . Then, since $H(\boldsymbol{\pi})^*+\mathcal{I}$ is the generator of a birth–death process with zero killing rates, all the eigenvalues of $H(\boldsymbol{\pi})^*+\mathcal{I}$ are negative [Reference Ledermann, Reuter and Mahler11]. Hence, all the eigenvalues of $H(\boldsymbol{\pi})$ are less than $-1$ . As a consequence, we have the following result due to the fact that all the eigenvalues are negative.

Lemma 4. The unique solution to (21) is given by $(\mathbf{s}(t),t\geq 0)$ , where $\mathbf{s}(t)={\text{e}}^{H(\boldsymbol{\pi})t}\mathbf{s}(0)$ . Furthermore, $(\mathbf{s}(t),t\geq 0)$ satisfies, for some $\delta_2>0$ and $D_4<\infty$ ,

\begin{equation*} {\|{\mathbf{s}(t)}\|}_2\leq {\text{e}}^{-\delta_2 t}D_4{\|{\mathbf{s}(0)}\|}_2. \end{equation*}

From Lemma 4 and the unique solution given in Theorem 4, the following result follows immediately. Hence, we omit the proof.

We are now ready to state the main result on the FCLT for the stationary regime.

Proof. See Section 4.4.

Now we use Theorem 6 to provide an approximation to the average blocking probability in the system with N servers.

Theorem 7. Let $P_{\text{block}}^{(N)}$ be the average blocking probability in the system with N servers. Then,

(23) \begin{equation} P_{\text{block}}^{(N)}=\pi_C^d-\frac{1}{\sigma\sqrt{N}\,}\left(\sum_{i=0}^C i(\kappa_i-\kappa_{i+1})\right)-\frac{\beta}{\sigma\sqrt{N}\,}\big(1-\pi_C^d\big)+o\Big(N^{-\frac{1}{2}}\Big), \end{equation}

where the vector $\mathbf{\kappa}=(\kappa_i,0\leq i\leq C)=\int_{0}^{\infty}{\text{e}}^{H(\boldsymbol{\pi})s}W_3(\boldsymbol{\pi})\,{\text{d}} s$ is the mean of the unique solution of (22) in the stationary regime.

Proof. See Section 4.5.

The results for the case $\beta=0$ were presented in [Reference Vasantam and Mazumdar19].

The significance of Theorem 7 is that although the exact blocking formula for $P_{\text{block}}^{(N)}$ is not known, and is also difficult to characterize due to the complex interactions between servers, when N becomes large we can compute approximations to the blocking probability as a function of $\boldsymbol{\pi}$ , $\beta$ , N, $\sigma$ , and C.

4.1. Proof of Lemma 1

We have

(24) \begin{align}\mathbf{Z}^{(N)}(t)= & \mathbf{Z}^{(N)}(0)+\int_{s=0}^t\sqrt{N}(W(\mathbf{X}^{(N)}(s))-W(\mathbf{x}(s,\mathbf{u})))\,{\text{d}} s\\ &- \int_{s=0}^tW_3(\mathbf{X}^{(N)}(s))\,{\text{d}} s+\mathbf{M}^{(N)}(t) \nonumber\end{align}

and, for $n\geq 1$ ,

\begin{equation*}\langle\mathbf{M}^{(N)}_n\rangle_t=\int_{s=0}^t((W_1(\mathbf{X}^{(N)}(s)))_n+(W_2(\mathbf{X}^{(N)}(s)))_n-\frac{1}{\sqrt{N}\,}(W_3(\mathbf{X}^{(N)}(s)))_n)\,{\text{d}} s.\end{equation*}

From (24), we obtain

\begin{equation*}{\|{\mathbf{Z}^{(N)}(t)}\|}_2\leq {\|{\mathbf{Z}^{(N)}(0)}\|}_2+B_W\int_{s=0}^t{\|{\mathbf{Z}^{(N)}(s)}\|}_2\,{\text{d}} s + 2C\beta t+{\|{\mathbf{M}^{(N)}(t)}\|}_2.\end{equation*}

By using Gronwall’s lemma,

\begin{equation*}{\|{\mathbf{Z}^{(N)}(t)}\|}_2\leq \big({\|{\mathbf{Z}^{(N)}(0)}\|}_2+2C\beta t+{\|{\mathbf{M}^{(N)}(t)}\|}_2\big){\text{e}}^{B_Wt}.\end{equation*}

As a result, we get

\begin{equation*}{\|{\mathbf{Z}^{(N)}(t)}\|}_2^2\leq 3\big({\|{\mathbf{Z}^{(N)}(0)}\|}_2^2+4C^2\beta^2 t^2+{\|{\mathbf{M}^{(N)}(t)}\|}_2^2\big){\text{e}}^{2B_Wt}.\end{equation*}

For $T>0$ , we have

\begin{equation*}\sup_{0\leq t\leq T}{\|{\mathbf{Z}^{(N)}(t)}\|}_2^2\leq 3\Big({\|{\mathbf{Z}^{(N)}(0)}\|}_2^2+4C^2\beta^2T^2+\sup_{0\leq t\leq T}{\|{\mathbf{M}^{(N)}(t)}\|}_2^2\Big){\text{e}}^{2B_WT}.\end{equation*}

Finally, Doob’s inequality implies that

\begin{equation*}{\mathbb{E}\left[{\sup_{0\leq t \leq T}{\|{\mathbf{Z}^{(N)}(t)}\|}_2^2}\right]}\leq 3{\text{e}}^{2B_WT}\left({\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(0)}\|}_2^2}\right]}+4C^2\beta^2T^2+4{\mathbb{E}\left[{\sum_{i=1}^C\langle\mathbf{M}^{(N)}_i\rangle_T}\right]}\right).\end{equation*}

Since $\sup_{N\geq 1}{\mathbb{E}\left[{\sum_{i=1}^C\langle\mathbf{M}^{(N)}_i\rangle_T}\right]}<\infty$ , we conclude that

\begin{equation*} \limsup_{N\to\infty}{\mathbb{E}\left[{\sup_{0\leq t \leq T}{\|{\mathbf{Z}^{(N)}(t)}\|}_2^2}\right]}<\infty.\end{equation*}

4.2 Proof of Theorem 3

We recall that since the space $\mathbb{V}$ is a Polish space, the space of càdlàg functions under the Skhorohod topology is a Polish space [Reference Ethier and Kurtz4, Theorem 5.6, p. 121]. Hence, from Prohorov’s theorem [Reference Ethier and Kurtz4], tightness is equivalent to relative compactness. Therefore, it is enough to show the tightness, and then we need to show that every limiting point has the same law as the unique OU process that solves (20) with the initial point $\mathbf{Z}(0)$ .

We use Theorem 4.1 of [Reference Ethier and Kurtz4, p. 354] to show the tightness of $(\mathbf{Z}^{(N)}(t),t\geq 0)$ . First, we establish several useful preliminary results.

Since $\mathbf{Z}^{(N)}(0)\Rightarrow \mathbf{Z}(0)$ , it implies that the sequence $\mathbf{Z}^{(N)}(0)$ is tight. Let $\tilde{B}(r)$ be the closed ball with radius r centered at $\mathbf{0}$ . For every $\varepsilon>0$ , there exists $r_{\varepsilon}<\infty$ such that $\mathbb{P}(\mathbf{Z}^{(N)}(0)\in \tilde{B}(r_{\varepsilon}))>1-\varepsilon$ for all $N\geq 1$ . We now define a random variable $\mathbf{X}^{(N,\varepsilon)}(0)$ such that it coincides with $\mathbf{X}^{(N)}(0)$ on $\{\mathbf{Z}^{(N)}(0)\in \tilde{B}(r_{\varepsilon})\}$ and $\mathbf{Z}^{(N,\varepsilon)}(0)$ is uniformly bounded in N on $\{\mathbf{Z}^{(N)}(0)\notin \tilde{B}(r_{\varepsilon})\}$ . Then, by using coupling arguments, the processes $(\mathbf{Z}^{(N,\varepsilon)}(t),t\geq 0)$ and $(\mathbf{Z}^{(N)}(t),t\geq 0)$ coincide on $\{\mathbf{Z}^{(N)}(0)\in \tilde{B}(r_{\varepsilon})\}$ . Hence, without loss of generality, we assume that $\mathbf{Z}^{(N)}(0)$ is uniformly bounded in N. As a consequence, the result stated in Lemma 1 can be used in the rest of the proof.

Next, we state the following useful result from [Reference Graham9, Lemma 3.3]. For a and h in $\mathbb{R}$ , let $\hat{B}(a,h)=(a+h)^d-a^d-da^{d-1}h$ . Then, if both a and $a+h$ lie in [0, 1], we have

(25) \begin{equation}0\leq \hat{B}(a,h)\leq h^d+(2^d-d-2)ah^2.\end{equation}

We now define a mapping $\hat{G}\,:\,\mathbb{U}\times\mathbb{V}\mapsto\mathbb{V}$ as follows: for $\mathbf{r}\in\mathbb{U}$ and $\mathbf{y}\in\mathbb{V}$ ,

\begin{equation*}(\hat{G}(\mathbf{r},\mathbf{y}))_n=\sigma \hat{B}(r_{n-1},y_{n-1})-\sigma \hat{B}(r_{n},y_{n}).\end{equation*}

Then, if $\mathbf{r}+\mathbf{y}\in\mathbb{U}$ , we have

(26) \begin{equation}W(\mathbf{r}+\mathbf{y})-W(\mathbf{r})=H(\mathbf{r})\mathbf{y}+\hat{G}(\mathbf{r},\mathbf{y}).\end{equation}

Note that since $\mathbf{Z}^{(N)}(t)=\sqrt{N}(\mathbf{X}^{(N)}(t)-\mathbf{x}(t,\mathbf{u}))$ , we have

\begin{equation*}\mathbf{X}^{(N)}(t)=\mathbf{x}(t,\mathbf{u})+\frac{\mathbf{Z}^{(N)}(t)}{\sqrt{N}}.\end{equation*}

Here, $\mathbf{X}^{(N)}(t)$ , $\mathbf{x}(t,\mathbf{u})\in\mathbb{U}$ and $\frac{\mathbf{Z}^{(N)}(t)}{\sqrt{N}}\in\mathbb{V}$ . Hence, we can write the following equation from (26):

(27) \begin{equation}W(\mathbf{X}^{(N)}(t)-W(\mathbf{x}(t,\mathbf{u}))=H(\mathbf{x}(t,\mathbf{u}))\frac{\mathbf{Z}^{(N)}(t)}{\sqrt{N}}+\hat{G}\left(\mathbf{x}(t,\mathbf{u}),\frac{\mathbf{Z}^{(N)}(t)}{\sqrt{N}}\right).\end{equation}

Now we show that the conditions of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354] are satisfied. Let us write $\mathbf{B}_n(t)$ and $\mathbf{A}_n^{i,j}(t)$ of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354] as $\mathbf{D}_{(N)}(t)$ and $\mathbf{A}_{(N)}^{i,j}(t)$ , respectively. Then, from (13) we can write

\begin{equation*}\mathbf{D}_{(N)}(t)=\int_{s=0}^t\sqrt{N}(W(\mathbf{X}^{(N)}(s))-W(\mathbf{x}(s,\mathbf{u})))\,{\text{d}} s-\int_{s=0}^tW_3(\mathbf{X}^{(N)}(s))\,{\text{d}} s.\end{equation*}

From (27), we can also write

\begin{align*}\mathbf{D}_{(N)}(t)= & \int_{s=0}^tH(\mathbf{x}(s,\mathbf{u}))\mathbf{Z}^{(N)}(s)\,{\text{d}} s+\int_{s=0}^t\sqrt{N}\hat{G}\left(\mathbf{x}(s,\mathbf{u}),\frac{\mathbf{Z}^{(N)}(s)}{\sqrt{N}}\right)\,{\text{d}} s \\ & -\int_{s=0}^tW_3(\mathbf{X}^{(N)}(s))\,{\text{d}} s.\end{align*}

Also, from (12) we have

\begin{equation*}\mathbf{A}_{(N)}^{i,i}(t)=\int_{s=0}^t\left((W_1(\mathbf{X}^{(N)}(s)))_i+(W_2(\mathbf{X}^{(N)}(s)))_i-\frac{1}{\sqrt{N}\,}(W_3(\mathbf{X}^{(N)}(s)))_i\right)\,{\text{d}} s,\end{equation*}

and $\mathbf{A}_{(N)}^{i,j}(t)=0$ for $i\neq j$ . Let us define

\begin{align*} \mathbf{D}(t) & = \int_{s=0}^tH(\mathbf{x}(s,\mathbf{u}))\mathbf{Z}(s)\,{\text{d}} s - \int_{s=0}^tW_3(\mathbf{x}(s,\mathbf{u}))\,{\text{d}} s, \\ A^{i,i}(t) & = \int_{s=0}^t\left((W_1(\mathbf{x}(s,\mathbf{u})))_i+(W_2(\mathbf{x}(s,\mathbf{u})))_i\right)\,{\text{d}} s,\end{align*}

with $A^{i,j}(t)=0$ for $i\neq j$ .

Since $(\mathbf{Z}_i^{(N)}(t),t\geq 0)$ has jumps of size $\frac{1}{\sqrt{N}\,}$ , and from the continuity of $\mathbf{D}_{(N)}(t)$ and $\mathbf{A}_{(N)}^{i,j}(t)$ in t, conditions $(4.3)$ – $(4.5)$ of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354] are valid. From condition $(4.6)$ of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354], we need to show that, for $T>0$ , $\mathbf{V}_{(N)}(T)\to 0$ in probability as $N\to\infty$ , where $\mathbf{V}_{(N)}(t)=\mathbf{D}_{(N)}(t)-\mathbf{D}(t)$ is given by

\begin{align*} \mathbf{V}_{(N)}(t) = & \left(\int_{s=0}^tH(\mathbf{x}(t,\mathbf{u}))\mathbf{Z}^{(N)}(t)\,{\text{d}} s + \int_{s=0}^t\sqrt{N}\hat{G}\left(\mathbf{x}(t,\mathbf{u}),\frac{\mathbf{Z}^{(N)}(t)}{\sqrt{N}}\right){\text{d}} s \right. \\ & \left.- \int_{s=0}^tW_3\big(\mathbf{X}^{(N)}(s)\big)\,{\text{d}} s\right) \\ & - \left(\int_{s=0}^tH(\mathbf{x}(s,\mathbf{u}))\mathbf{Z}^{(N)}(s)\,{\text{d}} s - \int_{s=0}^tW_3(\mathbf{x}(s,\mathbf{u}))\,{\text{d}} s\right).\end{align*}

From (14) and (25), we have $\int_{s=0}^T\sqrt{N}\hat{G}\left(\mathbf{x}(s,\mathbf{u}),\frac{\mathbf{Z}^{(N)}(s)}{\sqrt{N}}\right){\text{d}} s\to 0$ in probability. From the existence of the mean-field limit, it follows that $\int_{s=0}^TW_3(\mathbf{X}^{(N)}(s))\,{\text{d}} s - \int_{s=0}^TW_3(\mathbf{x}(s,\mathbf{u}))\,{\text{d}} s\to 0$ in probability. Thus, condition $(4.6)$ of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354] is also true. Condition $(4.7)$ of [Reference Ethier and Kurtz4, Theorem 4.1, p. 354] is also true since $\mathbf{A}_{(N)}^{i,i}(T)-A^{i,i}(T)\to 0$ in probability from the existence of the mean-field limit. Finally, from the fact that the SDE (20) has a unique solution, the sequence $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ converges to the unique solution of the SDE (20) as $N\to\infty$ .

4.3. Proof of Lemma 3

The proof is based on Lemma 2. Although our main objective is to establish the convergence of $\{\mathbf{Z}^{(N)}(t)\}_{N\geq 1}$ as $N\to\infty$ in the stationary regime, we prove the tightness of the stationary sequence $\{\mathbf{Z}^{(N)}(\infty)\}_{N\geq 1}$ by studying an alternative process $(\mathbf{Q}^{(N)}(t),t\geq 0)$ with the help of Lemma 2, where $\mathbf{Q}^{(N)}(t)=\sqrt{N}(\mathbf{X}^{(N)}(t)-\boldsymbol{\pi})$ . Let us write the solution to the mean-field equation (8) at time h with the initial point $\mathbf{v}$ as $\mathbf{y}(h,\mathbf{v})$ . We have

\begin{equation*} \mathbf{y}(h,\mathbf{v})=\mathbf{v}+\int_{s=0}^hW(\mathbf{y}(s,\mathbf{v}))\,{\text{d}} s. \end{equation*}

Also, the process $(\mathbf{X}^{(N)}(t),t\geq 0)$ satisfies

(28) \begin{equation} \mathbf{X}^{(N)}(t)=\mathbf{X}^{(N)}(0)+\int_{s=0}^tW(\mathbf{X}^{(N)}(s))\,{\text{d}} s-\int_{s=0}^t\frac{1}{\sqrt{N}\,}W_3\big(\mathbf{X}^{(N)}(s)\big)\,{\text{d}} s+\frac{\mathbf{M}^{(N)}(t)}{\sqrt{N}}. \end{equation}

For $t_0\geq 0$ , we obtain

\begin{equation*} \mathbf{Q}^{(N)}(t_0+h)=\sqrt{N}\big(\mathbf{X}^{(N)}(t_0+h)-\mathbf{y}\big(h,\mathbf{X}^{(N)}(t_0)\big)\big)\ +\sqrt{N}\big(\mathbf{y}\big(h,\mathbf{X}^{(N)}(t_0)\big)-\boldsymbol{\pi}\big). \end{equation*}

By defining $\mathbf{Z}^{(N)}(t_0, h)=\sqrt{N}\big(\mathbf{X}^{(N)}(t_0+h)-\mathbf{y}\big(h,\mathbf{X}^{(N)}(t_0)\big)\big)$ , we obtain

(29) \begin{equation} \mathbf{Q}^{(N)}(t_0+h)=\mathbf{Z}^{(N)}(t_0, h) +\sqrt{N}\big(\mathbf{y}\big(h,\mathbf{X}^{(N)}(t_0)\big)-\boldsymbol{\pi}\big). \end{equation}

Also, from (4.3),

(30) \begin{equation} \mathbf{Z}^{(N)}(t_0,h)=\sqrt{N}\big(\mathbf{X}^{(N)}(t_0+h)-\mathbf{X}^{(N)}(t_0)\big)-\sqrt{N}\int_{s=0}^t W\big(\mathbf{y}\big(s,\mathbf{X}^{(N)}(t_0)\big)\big)\,{\text{d}} s. \end{equation}

From (28), we have

(31) \begin{align} \mathbf{X}^{(N)}(t_0+h)-\mathbf{X}^{(N)}(t_0)= &\int_{s=t_0}^{t_0+h}W(\mathbf{X}^{(N)}(s))\,{\text{d}} s+\frac{\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)}{\sqrt{N}}\\[3pt] & -\frac{1}{\sqrt{N}\,}\int_{s=t_0}^{t_0+h}W_3(\mathbf{X}^{(N)}(s))\,{\text{d}} s.\nonumber \end{align}

Then, by using (30) and (31), we get

\begin{align*} \mathbf{Z}^{(N)}(t_0,h)= & \sqrt{N}\int_{s=0}^{h}W\big(\mathbf{X}^{(N)}(t_0+s)\big)\,{\text{d}} s - \sqrt{N}\int_{s=0}^{h}W\big(\mathbf{y}\big(s,\mathbf{x}^{(N)}(t_0)\big)\big)\,{\text{d}} s \\ & + \big(\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)\big) -\int_{s=t_0}^{t_0+h}W_3\big(\mathbf{X}^{(N)}(s)\big)\,{\text{d}} s. \end{align*}

After simplifications, we can write

\begin{align*} {\|{\mathbf{Z}^{(N)}(t_0,h)}\|}_2 \leq & \sqrt{N}\int_{s=0}^h{\big\|{W\big(\mathbf{X}^{(N)}(t_0+s)\big)-W\big(\mathbf{y}\big(s,\mathbf{X}^{(N)}(t_0)\big)\big)}\big\|}_2\,{\text{d}} s \\[3pt] & +{\|{\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)}\|}_2+{\left \vert {\beta} \right \vert} d\sqrt{C}h. \end{align*}

Hence, we obtain

\begin{equation*} {\|{\mathbf{Z}^{(N)}(t_0,h)}\|}_2\leq B_W\int_{s=0}^h{\|{\mathbf{Z}^{(N)}(t_0,s)}\|}_2\,{\text{d}} s+{\left \vert {\beta} \right \vert}d\sqrt{C}h+{\|{\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)}\|}_2. \end{equation*}

For any $T\geq 0$ , Gronwall’s inequality implies that there exists a constant $S_T$ such that

(32) \begin{equation} \sup_{0\leq h\leq T}{\|{\mathbf{Z}^{(N)}(t_0,h)}\|}_2\leq S_T\Big(S_T+\sup_{0\leq h\leq T}{\|{\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)}\|}_2\Big). \end{equation}

By using Lemma 2, we can write, from (29),

(33) \begin{equation} {\|{\mathbf{Q}^{(N)}(t_0+h)}\|}_2\leq{\|{\mathbf{Z}^{(N)}(t_0,h)}\|}_2+{\text{e}}^{-\delta_1 h}D_3{\|{\mathbf{Q}^{(N)}(t_0)}\|}_2. \end{equation}

From (32) and (33), we obtain

\begin{equation*} {\|{\mathbf{Q}^{(N)}(t_0+h)}\|}_2\leq S_T\Big(S_T+\sup_{0\leq h\leq T}{\|{\mathbf{M}^{(N)}(t_0+h)-\mathbf{M}^{(N)}(t_0)}\|}_2\Big)+{\text{e}}^{-\delta_1 h}D_3{\|{\mathbf{Q}^{(N)}(t_0)}\|}_2. \end{equation*}

As a result, we can find some constant $L_T$ as a function of T such that, for $0\leq h\leq T$ ,

(34) \begin{equation} {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t_0+h)}\|}_2^2}\right]} \leq L_T +3{\text{e}}^{-2\delta_1h}D_3^2{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t_0)}\|}^2_2}\right]}. \end{equation}

We now select a large value of T such that $3{\text{e}}^{-2\delta_1 T}D_3^2\leq \varepsilon<1$ . Then, for all $N\geq 1$ and an integer m, we get

\begin{equation*} {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}((m+1)T)}\|}_2^2}\right]}\leq L_T+\varepsilon {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(mT)}\|}^2_2}\right]}. \end{equation*}

By using the induction method, we obtain

(35) \begin{align} {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(mT)}\|}_2^2}\right]} & \leq L_T\Bigg(\sum_{j=1}^m\varepsilon^{j-1}\Bigg)+\varepsilon^m {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]} \nonumber \\ & \leq \frac{L_T}{1-\varepsilon}+{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]}. \end{align}

However, from (34),

\begin{equation*} \sup_{0\leq h\leq T}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(mT+h)}\|}_2^2}\right]} \leq L_T+3D_3^2 {\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(mT)}\|}^2_2}\right]}. \end{equation*}

As a consequence, (35) implies

\begin{equation*} \sup_{0\leq h\leq T}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(mT+h)}\|}_2^2}\right]} \leq L_T +3D_3^2 \left(\frac{L_T}{1-\varepsilon}+{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]}\right). \end{equation*}

Since m is arbitrary, we conclude that

\begin{equation*} \sup_{t\geq 0}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t)}\|}_2^2}\right]} \leq L_T +3D_3^2 \left(\frac{L_T}{1-\varepsilon}+{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]}\right). \end{equation*}

From ergodicity and Fatou’s Lemma [Reference Ethier and Kurtz4, p. 492], the stationary random variable $\mathbf{Z}^{(N)}(\infty)$ satisfies

\begin{equation*} {\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(\infty)}\|}_2^2}\right]} \leq \liminf_{t\geq 0}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t)}\|}_2^2}\right]} \leq \sup_{t\geq 0}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(t)}\|}_2^2}\right]}. \end{equation*}

Finally, to show that $\limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(\infty)}\|}_2^2}\right]}<\infty$ , we need to find an $\mathbf{X}^{(N)}(0)$ such that $\limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]}<\infty$ . For $n\geq 1$ , if we select $\mathbf{X}^{(N)}_n(0)=\frac{j}{N}$ so that $\frac{-1}{2N}\leq \pi_n-\frac{j}{N}\leq \frac{1}{2N}$ , then ${\mathbb{E}\left[{{\|{\mathbf{Q}^{(N)}(0)}\|}^2_2}\right]}<\frac{C}{4N}$ . Hence, completing the proof, $$\lim {\sup _{N \to \infty }} \mathbb{E} \left[ {{{\bf{Q}}^{(N)}}(0)_2^2} \right] = 0$$ .

4.4. Proof of Theorem 6

From Lemma 3 and the Markov inequality, the sequence $\{\mathbf{Z}^{(N)}(0)\}_{N\geq 1}$ is tight. As a result, from the Prohorov theorem [Reference Ethier and Kurtz4, p. 104], the sequence $\{\mathbf{Z}^{(N)}(0)\}_{N\geq 1}$ is relatively compact. Consider a converging subsequence and let $\mathbf{Z}^{(\infty)}(0)$ be its limiting point, which is square integrable. Then, from Theorem 3, the considered converging subsequence converges in law to the unique OU process $(\mathbf{Z}^{(\infty)}(t),t\geq 0)$ that solves the SDE (22) with initial point $\mathbf{Z}^{(\infty)}(0)$ . But, we know from [Reference Ethier and Kurtz4, Lemma 7.7 and Theorem 7.8, p. 131] that the limit of a sequence of stationary processes is stationary. Hence, the law of $(\mathbf{Z}^{(\infty)}(t),t\geq 0)$ should be the unique law of the stationary OU process solving the SDE (22). This argument applies for every converging subsequence. Hence, the sequence $\{(\mathbf{Z}^{(N)}(t),t\geq 0)\}_{N\geq 1}$ in the stationary regime converges to the unique stationary OU process solving the SDE (22). This completes the proof.

4.5. Proof of Theorem 7

The proof is based on Little’s law [Reference Asmussen2, Theorem 4.1], and Theorems 5 and 6. Let us consider a random variable ${\tilde{\mathbf{S}}}^{(N)}(\infty)$ which denotes the number of progressing jobs in the system in the stationary regime. From Little’s law, we obtain $(N\lambda^{(N)})\big(1-P_{\text{block}}^{(N)}\big)={\mathbb{E}\left[{{\tilde{\mathbf{S}}}^{(N)}(\infty)}\right]}$ . We now consider a random variable $\mathbf{X}^{(N)}(\infty)=(\mathbf{X}_n^{(N)}(\infty),0\leq n\leq C)$ , where $\mathbf{X}_{n}^{(N)}(\infty)$ is the fraction of servers with at least n progressing jobs in the stationary regime. Then, we can write ${\tilde{\mathbf{S}}}^{(N)}(\infty)=\sum_{n=0}^C Nn( \mathbf{X}^{(N)}_n(\infty)-\mathbf{X}^{(N)}_{n+1}(\infty))$ . Hence, we have

\begin{equation*} \big(N\lambda^{(N)}\big)\big(1-P_{\text{block}}^{(N)}\big)={\mathbb{E}\left[{\sum_{n=0}^C Nn\Big(\mathbf{X}^{(N)}_n(\infty)-\mathbf{X}^{(N)}_{n+1}(\infty)\Big)}\right]}. \end{equation*}

Therefore, we obtain

(36) \begin{equation} \lambda^{(N)}\big(1-P_{\text{block}}^{(N)}\big)=\sum_{n=0}^C n {\mathbb{E}\left[{\mathbf{X}^{(N)}_n(\infty)-\mathbf{X}^{(N)}_{n+1}(\infty)}\right]}. \end{equation}

Also, from Theorem 5, Lemma 3, and Theorem 6, since the diffusion limit in the stationary regime has the mean vector $\kappa=\int_{0}^{\infty}{\text{e}}^{H(\boldsymbol{\pi})s}W_3(\boldsymbol{\pi})\,{\text{d}} s$ and $\limsup_{N\to\infty}{\mathbb{E}\left[{{\|{\mathbf{Z}^{(N)}(\infty)}\|}^2}\right]}<\infty$ , we have $\lim_{N\to\infty}\sqrt{N}\left({\mathbb{E}\left[{\mathbf{X}^{(N)}(\infty)}\right]}-\boldsymbol{\pi}\right)=\kappa$ . Therefore,

(37) \begin{equation} {\mathbb{E}\left[{\mathbf{X}^{(N)}(\infty)}\right]}-\boldsymbol{\pi}=\frac{\kappa}{\sqrt{N}\,}+ o\Big(N^{-\frac{1}{2}}\Big). \end{equation}

Then, from (36) and (37), we get

(38) \begin{equation} \lambda^{(N)}\big(1-P_{\text{block}}^{(N)}\big)=\sum_{n=0}^C n(\pi_n-\pi_{n+1})+\frac{1}{\sqrt{N}\,}\sum_{n=0}^C n(\kappa_n-\kappa_{n+1})+o\Big(N^{-\frac{1}{2}}\Big). \end{equation}

But, from the stationary mean-field equations, the fixed point satisfies

(39) \begin{equation} \pi_0=1, \qquad \sum_{n=0}^C n(\pi_n-\pi_{n+1})=\sigma\big(1-\pi_C^d\big). \end{equation}

Then, from (38) and (39), we can write

\begin{equation*} 1-P_{\text{block}}^{(N)}=\frac{1}{\lambda^{(N)}}\left[\sigma(1-\pi_C^d)+\frac{1}{\sqrt{N}\,}\sum_{n=0}^C n(\kappa_n-\kappa_{n+1})\right]+o\Big(N^{-\frac{1}{2}}\Big). \end{equation*}

However, by using $\lambda^{(N)}=\sigma-\frac{\beta}{\sqrt{N}\,}$ and the fact that $\frac{\beta}{\sigma\sqrt{N}\,}<1$ , we obtain

\begin{equation*} 1-P_{\text{block}}^{(N)}=\frac{1}{\sigma}\left(1+\frac{\beta}{\sigma\sqrt{N}\,}+o\Big(N^{-\frac{1}{2}}\Big)\right)\left[\sigma\big(1-\pi_C^d\big)+\frac{1}{\sqrt{N}\,}\sum_{n=0}^C n(\kappa_n-\kappa_{n+1})\right] +o\Big(N^{-\frac{1}{2}}\Big). \end{equation*}

After some simple calculations, we obtain

\begin{equation*} P_{\text{block}}^{(N)}=\pi_C^d-\frac{1}{\sigma\sqrt{N}\,}\left(\sum_{n=0}^C n(\kappa_n-\kappa_{n+1})\right)-\frac{\beta}{\sigma\sqrt{N}\,}\big(1-\pi_C^d\big)+o\Big(N^{-\frac{1}{2}}\Big). \end{equation*}

This completes the proof.

4.6. Proof of Lemma 2

The proof uses the quasi-monotonicity of the mean field. Let us write the unique solution to the MFEs with the initial point $\mathbf{v}$ as $(\mathbf{y}(t,\mathbf{v}),t\geq 0)$ . From the quasi-monotonicity of the mean field, we have

(40) \begin{equation} \mathbf{y}(t,\min(\mathbf{v},\boldsymbol{\pi}))\leq \mathbf{y}(t,\mathbf{v})\leq \mathbf{y}(t,\max(\mathbf{v},\boldsymbol{\pi})). \end{equation}

For $\mathbf{a},\mathbf{b}\in\mathbb{U}$ , let ${\|{\mathbf{a}-\mathbf{b}}\|}_1=\sum_{i=0}^C{\left \vert {a_i-b_i} \right \vert}$ . From Lemma 4 of [Reference Xie, Dong, Lu and Srikant22], since $\min(\mathbf{v},\boldsymbol{\pi})\leq \boldsymbol{\pi}$ and $\max(\mathbf{v},\boldsymbol{\pi})\geq \boldsymbol{\pi}$ , we have

(41) \begin{equation} \begin{split} {\|{\mathbf{y}(t,\min(\mathbf{v},\boldsymbol{\pi}))-\boldsymbol{\pi}}\|}_1 & \leq {\text{e}}^{-t}{\|{\min(\mathbf{v},\boldsymbol{\pi})-\boldsymbol{\pi}}\|}_1, \\ {\|{\mathbf{y}(t,\max(\mathbf{v},\boldsymbol{\pi}))-\boldsymbol{\pi}}\|}_1 & \leq {\text{e}}^{-t}{\|{\max(\mathbf{v},\boldsymbol{\pi})-\boldsymbol{\pi}}\|}_1. \end{split} \end{equation}

For $t\geq 0$ , let us define two sets $V_+(t)$ and $V_-(t)$ as

\begin{align*} V_+(t) & = \{i\,:\,y_i(t,\mathbf{v})\geq \pi_i\}, \\ V_-(t) & = \{i\,:\,y_i(t,\mathbf{v}) < \pi_i\}. \end{align*}

Then, we can write

\begin{equation*} {\|{\mathbf{y}(t,\mathbf{v})-\boldsymbol{\pi}}\|}_1=\sum_{i\in V_+(t)}(y_i(t,\mathbf{v})-\pi_i)+\sum_{j\in V_-(t)}(\pi_j-y_j(t,\mathbf{v})). \end{equation*}

From (40) and (41), we write

\begin{align*} {\|{\mathbf{y}(t,\mathbf{v})-\boldsymbol{\pi}}\|}_1 & \leq \sum_{i\in V_+(t)}(y_i(t,\max(\mathbf{v},\boldsymbol{\pi}))-\pi_i)+\sum_{j\in V_-(t)}(\pi_j-y_j(t,\min(\mathbf{v},\boldsymbol{\pi}))) \\ &\leq{\|{\mathbf{y}(t,\max(\mathbf{v},\boldsymbol{\pi}))-\boldsymbol{\pi}}\|}_1+{\|{\mathbf{y}(t,\min(\mathbf{v},\boldsymbol{\pi}))-\boldsymbol{\pi}}\|}_1, \\ &\leq {\text{e}}^{-t}{\|{\max(\mathbf{v},\boldsymbol{\pi})-\boldsymbol{\pi}}\|}_1+{\text{e}}^{-t}{\|{\min(\mathbf{v},\boldsymbol{\pi})-\boldsymbol{\pi}}\|}_1 \\ &= {\text{e}}^{-t}{\|{\mathbf{v}-\boldsymbol{\pi}}\|}_1. \end{align*}

Finally, the result follows from the fact that the norms ${\|{\cdot}\|}_1$ and ${\|{\cdot}\|}_2$ are equivalent.