2. System model and the routeing policy

We consider a system consisting of a large number N of parallel servers. Jobs arrive according to a Poisson process with rate $N\lambda$ and the job lengths are assumed to be i.i.d. from a general distribution $G({\cdot})$ defined on $\mathbb{R}_+$ . A central job dispatcher routes an incoming job to a server according to the SQ(d) policy defined below. We assume that each server has capacity to process up to a number C of jobs simultaneously and each job is processed at unit rate. At any time t, if a server is currently serving i jobs, then we say that the server has occupancy i and vacancy $C-i$ at time t. If an incoming job is routed to a server with occupancy C, then the job is blocked or discarded; otherwise the processing of the job begins immediately and it is processed at unit rate.

In Definition 2.1, we assume sampling with replacement because of notational convenience and it is easy to show that the asymptotic results that are of interest in the paper are not affected whether we sample with or without replacement.

We assume that the service times have finite mean $\frac{1}{\mu}$ and the service time distribution denoted by $G({\cdot})$ on $[0,\infty)$ possesses a continuous density $g({\cdot})$ . We assume that $\overline{G}({\cdot})$ is supported on $[0,\infty)$ , where $\overline{G}(.)$ denotes the complementary distribution. The hazard rate function of $G({\cdot})$ is defined as

\[ \beta(x)=\dfrac{g(x)}{\overline{G}(x)}= \dfrac{g(x)}{1-G(x)}\quad \text{for $x\in[0,\infty)$.} \]

The hazard rate function $\beta$ indicates the instantaneous rate at which the service of a job ends. More precisely, a job with age y (where y denotes the time since its arrival) at time t exits the server in the interval $[t,t+ d t)$ with probability $\beta(y) d t$ .

3. Notation and terminology

We first introduce the notation which is used throughout the paper. Let $\mathbb{Z}$ , $\mathbb{R}$ be the set of integers and real numbers, respectively. Further, let $\mathbb{Z}_+$ , $\mathbb{R}_+$ be the set of non-negative integers and non-negative real numbers, respectively.

3.1. Function and measure spaces

For any given metric space $\mathcal{E}$ , let $\mathcal{K}_{b}(\mathcal{E})$ , $\mathcal{C}_{b}(\mathcal{E})$ , and $\mathcal{C}_{s}(\mathcal{E})$ be the space of bounded measurable real-valued functions, the space of bounded continuous real-valued functions, and the space of continuous real-valued functions with compact support, defined on $\mathcal{E}$ , respectively. Furthermore, let $\mathcal{C}^1(\mathcal{E})$ be the space of once continuously differentiable real-valued functions defined on $\mathcal{E}$ and let the subspace of functions in $\mathcal{C}^1(\mathcal{E})$ which have compact support be denoted by $\mathcal{C}^1_{s}(\mathcal{E})$ . The space of bounded functions in $\mathcal{C}^1(\mathcal{E})$ whose first derivatives are also bounded is denoted by $\mathcal{C}^1_{b}(\mathcal{E})$ . For any function $f\in\mathcal{K}_{b}(\mathcal{E})$ , $h\in\mathcal{C}^1(\mathcal{E})$ , we define

\begin{equation*} \|{{f}}\|=\sup_{x\in \mathcal{E}}{\vert {f(x)} \vert},\quad \|{{h}}\|_1=\|{{h}}\|+\|{{h'}}\|{,} \end{equation*}

where h´ denotes the derivative of h. The space $\mathcal{C}_{b}(\mathcal{E})$ is equipped with the uniform topology, that is, we say that a sequence of functions $(f_n\in\mathcal{C}_b(\mathcal{E}),n\geq 1)$ converges to a function $f\in\mathcal{C}_b(\mathcal{E})$ if $\|{{f_n-f}}\|\to 0$ as $n\to\infty$ . The space $\mathcal{C}^1(\mathcal{E})$ is equipped with the topology induced by the norm $\|{{\cdot}}\|_1$ .

For a given metric space $\mathcal{E}$ , let the Borel $\sigma$ -algebra be denoted by $\mathcal{B}(\mathcal{E})$ . Let the space of finite non-negative measures on $\mathcal{E}$ be denoted by $\mathcal{M}_{F}(\mathcal{E})$ . We use the notation $\nu(B)$ and $\nu(\{y\})$ to denote the measure of a Borel set $B\in\mathcal{B}(\mathcal{E})$ and an element $y\in\mathcal{E}$ with respect to the measure $\nu\in\mathcal{M}_{F}(\mathcal{E})$ , respectively. The space of probability measures is denoted by $\mathcal{M}_{1}(\mathcal{E})$ . Also, let $\mathcal{M}_{1}^N(\mathcal{E})\subset \mathcal{M}_{1}(\mathcal{E})$ be the subspace of probability measures defined as $ \mathcal{M}_{1}^N(\mathcal{E})=\{\nu\in\mathcal{M}_{1}(\mathcal{E}) \colon N\,\nu{(B)}\in\mathbb{Z}_+,\ \text{{for all} } B\in\mathcal{B}(\mathcal{E})\}. $ For any $\phi\in\mathcal{K}_b(\mathcal{E})$ , $\nu\in\mathcal{M}_{F}(\mathcal{E})$ , we define

\begin{equation*} \langle\nu,\phi\rangle=\int_{\mathcal{E}}\phi(y)\nu({\mathrm{d}} y). \end{equation*}

The space of measures $\mathcal{M}_{F}(\mathcal{E})$ is equipped with the weak topology induced by the weak convergence of measures.

The age of an active job is the time elapsed since its arrival. To model the dynamics of an Erlang loss system with capacity C for each server by a Markov process, we define the state of each server as $(n,a_1,a_2,\ldots,a_n)$ , where n denotes the number of jobs that are in progress at the server and $a_i$ denotes the age of the ith job in progress. We now define a space $\mathcal{U}$ that is used in earlier works to study queuing models with general service time distributions by using the classical supplementary variable method [Reference Li and Lin27, Reference Sevastyanov38] such that it contains all the possible server states as elements. The space $\mathcal{U}$ is defined as

\begin{equation*} \mathcal{U}=\bigcup_{n=0}^C \mathcal{U}_n, \end{equation*}

where $\mathcal{U}_0=\{0\}$ and an element in $\mathcal{U}_n$ for $n\geq 1$ is of the form $(n,a_1,\ldots,a_n)$ , where $1\leq n\leq C$ and $a_i\in\mathbb{R}_+$ . The state of a server that has n jobs belongs to the space $\mathcal{U}_n$ . Here, one might omit the variable n and consider just $(a_1,\ldots,a_n)$ to represent a server state, but such representation does not account for idle servers while (0) is the state of idle servers in our representation. Furthermore, the variable n directly gives us information about the number of progressing jobs at a server which changes upon every arrival and departure. Hence, it is convenient to work with the server state representation that has a variable that denotes the number of progressing jobs at a server.

It is also possible to define an element in $\mathcal{U}_n$ by $(n,a_1,\ldots,a_n,0,\ldots,0)$ of size $C+1$ . This allows us to have constant size of $C+1$ for an element in $\mathcal{U}$ . Note that the zeros in the state $(n,a_1,\ldots,a_n,0,\ldots,0)$ act as dummy variables as there are only n jobs. Hence, to make it simple, we consider that an element in $\mathcal{U}_n$ is of the form $(n, a_1,\ldots,a_n)$ with size $n+1$ . Without loss of generality, we refer to an element in the set $\mathcal{U}$ by $\textbf{\textit{u}}$ and an element in the set $\mathcal{U}_n$ by $\textbf{\textit{u}}_n$ . Note that we have $\textbf{\textit{u}}_0=0$ . For $\textbf{\textit{y}}_n=(n,y_1,\ldots,y_{n}),\mathbf{z}_m=(m,z_1,\ldots,z_{m})$ , we define the metric $d_{\mathcal{U}}(\textbf{\textit{y}}_n,\mathbf{z}_m)$ as

\begin{equation*} \label{eq:metric} d_{\mathcal{U}}(\textbf{\textit{y}}_n,\mathbf{z}_m) = \begin{cases} \displaystyle \sum_{i=1}^{n}{ \vert {y_i-z_i} \vert} & \text{if } n=m,\\ \infty & \text{otherwise}. \end{cases} \end{equation*}

For $(n,u_1,\ldots,u_{n})\in\mathcal{U}_n$ and $y\geq 0$ , we use the following notation:

\begin{align*} \textbf{\textit{u}}_n&=(n,u_1,\ldots, u_n),\\ \textbf{\textit{u}}_n^{-j}&=(n-1,u_1,\ldots,u_{j-1},u_{j+1},\ldots,u_n),\\ (\textbf{\textit{u}}_n^{j};\,y)&=(n+1,u_1,\ldots,u_{j-1},y,u_j,\ldots,u_n),\\ (\textbf{\textit{u}}_n^{-j};\,y)&=(n,u_1,\ldots,u_{j-1},y,u_{j+1},\ldots,u_n). \end{align*}

For any Borel set $B\in\mathcal{B}(\mathcal{U})$ , let ${I_{\{{B}\}}}$ be the indicator function of B. Let the function $\textbf{1}$ be defined such that, for all $\textbf{\textit{u}}\in\mathcal{U}$ , we have $ \textbf{1}(\textbf{\textit{u}})=1. $

A function $f \colon \mathcal{U}\mapsto\mathbb{R}$ is said to be differentiable if, for any $n\geq 1$ , the function $\gpfrac{\partial f(\textbf{\textit{u}}_n)}{\partial u_i}$ exists for all $1\leq i\leq n$ at every $\textbf{\textit{u}}_n\in\mathcal{U}_n$ . As a result, the function ${I_{\{{\mathcal{U}_n}\}}},n\geq 1$ is differentiable. For a differentiable function $f \colon \mathcal{U}\mapsto\mathbb{R}$ , we have

\begin{equation*} \|{{f'}}\|=\max_{n\geq 1}\bigg(\sup_{\textbf{\textit{u}}_n\in\mathcal{U}_n}\bigg(\max_{1\leq i\leq n}{ \bigg\vert {\dfrac{\partial f(\textbf{\textit{u}}_n)}{\partial u_i}} \bigg\vert}\bigg)\bigg). \end{equation*}

Further, for a differentiable function $f \colon \mathcal{U}\mapsto\mathbb{R}$ , let the function $\nabla_1f$ be defined as

\begin{equation*} \nabla_1f(n,u_1,\ldots,u_n)=\nabla f\cdot\textbf{1}=\sum_{i=1}^n\dfrac{\partial f(\textbf{\textit{u}}_n)}{\partial u_i}. \end{equation*}

A measure $\nu\in\mathcal{M}_F(\mathcal{U})$ when it is restricted to $\mathcal{U}_0$ is a Dirac measure at $\{0\}$ satisfying $\nu(\mathcal{U}_0)=\nu(\{0\})$ . We say that a measure $\nu$ is absolutely continuous with respect to Lebesgue measure if $\nu(\{\textbf{\textit{x}}_n\})=0$ at every $\textbf{\textit{x}}_n\in \mathcal{U}_n$ for all $n\geq 1$ . For any Borel measurable function f that is defined on $\mathcal{U}$ , we define

\begin{equation*} \langle\nu,f\rangle =f(0)\nu(\{0\})+\sum_{n=1}^C\int_{ \mathcal{U}_n}f(\mathbf{z}_n)\nu({\mathrm{d}} \mathbf{z}_n). \end{equation*}

We now define the function $\mathcal{I} \colon \mathcal{U}\mapsto\mathbb{R}$ as follows:

\begin{equation*} \mathcal{I}(\textbf{\textit{x}}_n) = \begin{cases} \displaystyle \sum_{i=1}^{n}x_i & \text{if } n\geq 1,\\ 0 & \text{otherwise}. \end{cases} \end{equation*}

For $b\geq 0$ , let $\tau_b^+ \colon \mathcal{U}\mapsto\mathcal{U}$ be the transition operator defined as

\begin{equation*} \tau_b^+(\textbf{\textit{x}}_n) = \begin{cases} (n,x_1+b,\ldots,x_n+b) & n\geq 1,\\ 0 & \text{otherwise}. \end{cases} \end{equation*}

Similarly, for any $b\geq 0$ and $f\in\mathcal{K}_b(\mathcal{U})$ , let the mapping $\tau_b \colon \mathcal{K}_b(\mathcal{U})\mapsto \mathcal{K}_b(\mathcal{U})$ be defined as $ \tau_yf(\textbf{\textit{u}})=f(\tau_y^+\textbf{\textit{u}}). $ Also, for $b\geq 0$ , let the measure $\tau_b\nu\in\mathcal{M}_F(\mathcal{U})$ be defined such that, for any Borel set $B\in\mathcal{B}(\mathcal{U})$ , $ \tau_b\nu{(B)}=\nu(\tau_b^+(B)). $ For $\nu\in\mathcal{M}_F(\mathcal{U})$ , the measure $\tau_b\nu\in\mathcal{M}_F(\mathcal{U})$ satisfies $ \label{eq:translation_function} \langle\tau_b\nu,f\rangle=\langle\nu,\tau_bf\rangle $ for all $f\in\mathcal{K}_b(\mathcal{U})$ and the existence of the unique measure $\tau_b\nu$ follows from the Riesz–Markov–Kakutani theorem [37, Theorem 2.14].

4. State descriptor and system dynamics

We index a sequence of systems by N which denotes the total numbers of servers. Incoming jobs arrive according to a Poisson process with rate $N\lambda$ and the job lengths are i.i.d. from a common distribution $G({\cdot})$ defined on $\mathbb{R}_+$ . The state of a server is written as $\mathbf{a}_n=(n,a_1,\ldots,a_n)\in\mathcal{U}$ when there are n progressing jobs and the ith job has age $a_i$ for $1\leq i\leq n$ . A server with state $\mathbf{a}_n$ , say, can be viewed as a particle with the given state. Therefore the system evolution can be considered as the evolution of a system with N particles, where the interactions between particles takes place while routeing an arrival according to the SQ(d) routeing policy.

The age of a job that is in service at a server increases linearly with time at unit rate until its service is completed. We next describe the possible state of a server at time $t+h$ ( $h>0$ ) given that it has state $\mathbf{a}_n$ at time t. We assume that when h is small enough, in the interval $[t,t+h)$ , the probability of having multiple events of arrivals or departures is negligible. In the interval $[t,t+h)$ , if there is no arrival or departure at the given server, then the server state will be equal to $\tau_h^+(\mathbf{a}_n)$ at time $t+h$ . On the other hand, if the ith job expires in the interval $[t,t+h)$ , then the server state will be equal to $\tau_h^+(\mathbf{a}_n^{-i})$ at time $t+h$ . Considering arrivals, if there is an arrival to the server at time $t+r$ ( $0\leq r< h$ ), then the arriving job chooses its position uniformly at random out of $n+1$ possible positions, and if it chooses the jth position, then the server state will be equal to $((\tau_h^+(\mathbf{a}_n))^{ j};\,h-r)$ at time $t+h$ .

Let $S_{(i,t)}^N\in\mathcal{U}$ be the random variable that indicates the state of the server i at time t. One can think of considering $(S^N_{(1,t)},\ldots, S^N_{(N,t)})$ to denote the system state at time t, which is a Markovian representation of the system. The dimension of this state space increases with N as $N\to\infty$ , which is inconvenient to work with since our focus of interest is to study the asymptotic behavior of the system as $N\to\infty$ . Hence, we consider an alternative simple system state representation that can be used to describe the system evolution as the evolution of a Markov process. Note that the system is symmetric with respect to the servers as they are identical and the server identities do not play any role in the evolution with time. Therefore, to model the system evolution by a Markov process, we will show that it is enough to just keep track of the number of servers that lie in each state $\textbf{\textit{u}}\in\mathcal{U}$ in order to establish the mean field limit. Measure-valued Markov processes have also been used to study other interacting particle systems, as in [Reference Gärtner13], [Reference Oelschlager34], and [Reference Kotelenez and Kurtz26], where each particle state $\textbf{\textit{x}}\in\mathbb{R}^n$ , $n>1$ is viewed as a measure-valued Markov process. Following these works, we consider the following system state descriptor.

Definition 4.1. At time t, the state descriptor of the system with index N is a random measure given by

(4.1) \begin{equation} \eta_t^N=\sum_{i=1}^N\delta_{S_{(i,t)}^N}. \label{eqn2} \end{equation}

The interpretation of $\eta_t^N$ is that for any measurable function f defined on $\mathcal{U}$ , we have

\begin{equation*} \langle\eta_t^N,f\rangle=\sum_{i=1}^Nf(S_{(i,t)}^N). \end{equation*}

At time t, conditioned on server states $S_{(i,t)}^{(N)}=s_{(i,t)}$ , say, the system state can be represented by a measure $\nu$ defined as

(4.2) \begin{equation} \nu=\sum_{i=1}^N\delta_{s_{(i,t)}}. \label{eqn3} \end{equation}

For $\eta_t^N=\nu$ , an element $\textbf{\textit{y}}\in\mathcal{U}$ is an atom of $\nu$ if there exists at least one server with the state $\textbf{\textit{y}}$ at time t. The mass of an atom of $\nu$ denotes the number of servers lying at that atom at time t. As a result, since the number of interacting particles in the system is equal to N, the measure $\nu$ defined on $\mathcal{U}$ contains a finite number of atoms which is bounded by N. If all the servers have different states then the number of atoms is equal to N; otherwise the number of atoms is less than N. Let $V_t$ be the number of atoms at time t and let the ith atom be denoted by $\textbf{\textit{v}}^{(i)}_t$ . Further, let the mass of the atom $\textbf{\textit{v}}^{(i)}_t$ be denoted by $a^{(i)}_t$ . Here, $a^{(i)}_t$ denotes the number of servers that lie in the state $\textbf{\textit{v}}^{(i)}_t$ at time t and $a^{(i)}_t\geq 1$ . Hence, for time t, from (4.2), we can also write $\nu$ as

\begin{equation*} \nu=\sum_{i=1}^{V_t}a^{(i)}_t \delta_{\textbf{\textit{v}}^{(i)}_t}. \end{equation*}

For any Borel set $B\in\mathcal{B}(\mathcal{U})$ , the number of servers with ages lying in the set B is equal to $\eta_t^N(B)=\nu(B)=\langle\nu,{I_{\{{B}\}}}\rangle$ . We now define the measure of an element $\textbf{\textit{y}}_n=(n,y_1,\ldots,y_n)$ as below. Let

\[ B^{\epsilon}(\textbf{\textit{y}}_n)=\{(n,r_1,\ldots,r_n) \colon y_i\leq r_i<y_i+\epsilon, 1\leq i\leq n\}. \]

Then, as in [Reference Gromoll, Puha and Williams18], we define

\begin{equation*} \nu(\{\textbf{\textit{y}}_n\})=\lim_{\epsilon\to 0}\nu(B^{\epsilon}(\textbf{\textit{y}}_n)). \end{equation*}

Essentially, $\nu(\{\textbf{\textit{y}}_n\})$ indicates the number of servers with state $\textbf{\textit{y}}_n$ at time t and can be viewed as an occupation count. The notation $d \nu(\textbf{\textit{y}}_n)$ denotes the number of servers with state lying in the interval $[\textbf{\textit{y}}_n,\textbf{\textit{y}}_n+d\textbf{\textit{y}}_n)$ , where $d\textbf{\textit{y}}_n=(d y_1,\ldots,d y_n)$ and $\textbf{\textit{y}}_n+d\textbf{\textit{y}}_n$ is the vector addition of $\textbf{\textit{y}}_n$ and $d\textbf{\textit{y}}_n$ . If there is no server lying in the state $\textbf{\textit{y}}_n$ at time t, then $\nu(\{\textbf{\textit{y}}_n\})=0$ ; otherwise $\textbf{\textit{y}}_n$ is an atom with mass $\nu(\{\textbf{\textit{y}}_n\})$ . The number of servers that have n progressing jobs at time t is given by $\nu(\mathcal{U}_n)=\langle\nu,{I_{\{{\mathcal{U}_n}\}}}\rangle.$

We now obtain the probability that the destination server of an arrival lies in a particular state.

Lemma 4.1. At time t, given that the system state is $\nu$ , that is, $\eta_{t}^N=\nu$ , under the SQ(d) routeing policy, the probability that the destination server of an arrival at time t lies in the state $\mathbf{z}_n=(n,z_1,\ldots,z_n)$ , where $\mathbf{z}_n$ is an atom of $\nu$ , is given by

(4.3) \begin{equation} p_r(\nu \colon \mathbf{z}_n)=\dfrac{\nu(\{\mathbf{z}_n\})}{N}\dfrac{(\overline{R}_n(\frac{\nu}{N})^d-\overline{R}_{n+1}(\frac{\nu}{N})^d)}{(\overline{R}_n(\frac{\nu}{N})-\overline{R}_{n+1}(\frac{\nu}{N}))}, \label{eqn4} \end{equation}

where

\[ \overline{R}_n\bigg(\dfrac{\nu}{N}\bigg)=\sum_{j \colon j=n}^C\dfrac{\nu}{N}(\mathcal{U}_j) \]

represents the fraction of servers with at least n jobs.

Proof. When a potential destination server is chosen uniformly at random from N servers, it will have state $(n,z_1,\ldots,z_n)$ with probability $\gpfrac{\nu(\{(n,z_1,\ldots,z_n)\})}{N}$ . Suppose out of the d potential destination servers, say j servers have occupancy n and the remaining $d-j$ servers have occupancy at least $n+1$ . Further, out of the j $ (j\geq 1)$ potential destination servers with occupancy n, assume r $(r\geq 1)$ servers lie in the state $\mathbf{z}_n$ . Then the probability that the destination server is a server with state $\mathbf{z}_n$ is given by

\begin{equation*}\binom{d }{ j}\binom{j}{ r}\bigg({\dfrac{r}{j}}\bigg)\bigg(\dfrac{\nu(\{\mathbf{z}_n\})}{N}\bigg)^r \bigg(\dfrac{\nu(\{\mathcal{U}_n\})-\nu(\{\mathbf{z}_n\})}{N}\bigg)^{j-r} \bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d-j}. \end{equation*}

Finally, by summing over all the possible values of j ( $j\geq 1$ ) and r ( $r\geq 1$ ), we have

\begin{align*} & \sum_{j=1}^d\sum_{r=1}^j\binom{d }{ j}\binom{j}{ r}\bigg({\dfrac{r}{j}}\bigg)\bigg(\dfrac{\nu(\{\mathbf{z}_n\})}{N}\bigg)^r \bigg(\dfrac{\nu(\{\mathcal{U}_n\})-\nu(\{\mathbf{z}_n\})}{N}\bigg)^{j-r} \bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d-j}\\ & \qquad =\sum_{j=1}^d \binom{d}{ j}\dfrac{1}{j}\bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d-j}\bigg(\dfrac{\nu(\{\mathcal{U}_n\})}{N}\bigg)^{ j}\\ & \qquad\quad\, \times\bigg[\sum_{r=1}^jr\binom{j}{ r}\bigg(\dfrac{\frac{\nu(\{\mathbf{z}_n\})}{N}}{\frac{\nu(\{\mathcal{U}_n\})}{N}}\bigg)^r \bigg(\dfrac{\frac{\nu(\{\mathcal{U}_n\})}{N}-\frac{\nu(\{\mathbf{z}_n\})}{N}}{\frac{\nu(\{\mathcal{U}_n\})}{N}}\bigg)^{j-r}\bigg]. \end{align*}

The term inside the square brackets in the above equation is the average of a binomial random variable, and hence it is equal to

\[ j\bigg(\dfrac{\frac{\nu(\{\mathbf{z}_n\})}{N}}{\frac{\nu(\{\mathcal{U}_n\})}{N}}\bigg). \]

As a result, the above expression simplifies to

\begin{equation*} p_r(\nu \colon \mathbf{z}_n)=\bigg(\dfrac{\frac{\nu(\{\mathbf{z}_n\})}{N}}{\frac{\nu(\{\mathcal{U}_n\})}{N}}\bigg)\sum_{j=1}^d \binom{d}{ j}\bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d-j}\bigg(\dfrac{\nu(\{\mathcal{U}_n\})}{N}\bigg)^{ j}. \end{equation*}

We can further write

\begin{align*} & p_r(\nu \colon \mathbf{z}_n)\\ &\quad =\bigg(\dfrac{\frac{\nu(\{\mathbf{z}_n\})}{N}}{\frac{\nu(\{\mathcal{U}_n\})}{N}}\bigg)\bigg[\bigg(\sum_{j=0}^d \binom{d}{ j}\bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d-j}\bigg(\dfrac{\nu(\{\mathcal{U}_n\})}{N}\bigg)^{ j}\bigg) -\bigg(\sum_{i \colon i=n+1}^{C}\dfrac{\nu(\mathcal{U}_i)}{N}\bigg)^{d}\bigg]. \end{align*}

After simplifications, we get (4.3).

Remark 4.1. We can also interpret the expression of the $p_r(\nu\colon \mathbf{z}_n)$ as follows. The probability that all the potential destination servers have occupancy at least n and there exists at least one potential destination server with occupancy n is equal to $\overline{R}_n(\frac{\nu}{N})^d-\overline{R}_{n+1}(\frac{\nu}{N})^d$ . From the SQ(d) policy, the probability that the destination has occupancy n is equal to $\overline{R}_n(\frac{\nu}{N})^d-\overline{R}_{n+1}(\frac{\nu}{N})^d$ . From the list of the servers with occupancy n, the fraction of the servers with the state $\mathbf{z}_n$ is equal to

\[ \dfrac{(\frac{\nu(\{\mathbf{z}_n\})}{N})}{(\frac{\nu(\{\mathcal{U}_n\})}{N})}. \]

Therefore the probability that the destination server lies in the state $\mathbf{z}_n$ is equal to

\[ \dfrac{(\frac{\nu(\{\mathbf{z}_n\})}{N})}{(\frac{\nu(\{\mathcal{U}_n\})}{N})}\times\bigg(\overline{R}_n\bigg(\dfrac{\nu}{N}\bigg)^d-\overline{R}_{n+1}\bigg(\dfrac{\nu}{N}\bigg)^d\bigg) . \]

For the case of exponential job length distributions, $\mathcal{U}_n=\{n\}$ and $\mathbf{z}_n=n$ . Hence,

\[ p_r(\nu\colon \mathbf{z}_n)=\overline{R}_n\bigg(\dfrac{\nu}{N}\bigg)^d-\overline{R}_{n+1}\bigg(\dfrac{\nu}{N}\bigg)^d {,} \]

coinciding with the analysis for the exponential case in [Reference Mukhopadhyay, Mazumdar and Guillemin33] and [Reference Mitzenmacher28].

As is clear from (4.3), the routeing decision depends only on the number of servers lying in each possible server state. Hence, we get the evolution of the process $(\eta_t^N, t\geq 0)$ by tracking arrival events, routeing decisions, and departure events.

5. Main results

Our aim is to study the limit as $N\to \infty$ of the empirical measure of the server states. For this, we define a sequence of systems such that a system with index N has N servers that process the incoming jobs arriving according to a Poisson process with rate $N\lambda$ , and all other system parameters remain the same for all N as given in Section 2. The system consists of a central job dispatcher that routes an arrival to a server according to the SQ(d) policy. For given N, the process $(\eta_t^N,t\geq 0)$ defined in (4.1) describes the dynamics of the system with index N. The goal is to characterize the limit of the normalized process $(\overline{\eta}_t^N,t\geq 0)$ as $N\to\infty$ , where

\begin{equation*} %%\label{eq:normalized_descriptor} \overline{\eta}_t^N=\dfrac{\eta_t^N}{N}. \end{equation*}

For a Borel set $B\in\mathcal{B}(\mathcal{U})$ , $\overline{\eta}_t^N(B)$ is equal to the fraction of the servers with state lying in the set B at time t.

5.2. Transient regime

In this section, we discuss the results on the transient regime. For given system parameters $\lambda,C,d$ , and the probability density function $g({\cdot})$ of the service time distributions, in Proposition 5.1 we state the mean field equations. The dynamics of a mean field solution $(\overline{\eta}_t,t\geq 0)$ are described by using a set of evolution equations of the real-valued processes $(\langle \overline{\eta}_t,\phi\rangle,t\geq 0)$ for all $\phi\in\mathcal{C}_b^1(\mathcal{U})$ , referred to as the mean field equations.

Proposition 5.1. (Mean field equations.) For given system parameters $(\lambda,C,d,g({\cdot}))$ , the process $(\overline{\eta}_t,t\geq 0)$ satisfies the following.

1. (i) The mapping $t \mapsto \overline{\eta}_t$ is a continuous.

2. (ii) For $\phi\in \mathcal{C}_b^1(\mathcal{U})$ , the process $(\overline{\eta}_t,t\geq 0)$ satisfies

(5.1) \begin{align} \langle\overline{\eta}_t,\phi\rangle & =\langle\overline{\eta}_0,\phi\rangle+\int_{s=0}^t \langle\overline{\eta}_s,\nabla_1\phi\rangle\,{\mathrm{d}} s \notag \\ & \quad\, -\int_{s=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j)(\phi(\textbf{\textit{x}}_n^{-j})-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \overline{\eta}_s(\textbf{\textit{x}}_n)\bigg)\,{\mathrm{d}} s \notag \\ &\quad\, +\int_{s=0}^t\bigg(\overline{\eta}_s(\{0\})\lambda\Phi_0(\overline{\eta}_s)(\phi(1,0)-\phi(0)) +\sum_{n=1}^{C-1}\sum_{i=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)} \notag \\ & \quad\, \times\lambda \Phi_n(\overline{\eta}_s) (\phi(\textbf{\textit{x}}_n^i;\,0)-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \overline{\eta}_s(\textbf{\textit{x}}_n)\bigg) \,{\mathrm{d}} s, \label{eqn5}\end{align}

where the index j is used to denote the position of the departing job when there are n progressing jobs and i denotes the position of the arriving job when there are already n progressing jobs at the server. Further,

\[ \Phi_n(\overline{\eta}_s)=\dfrac{(\overline{R}_n(\overline{\eta}_s)^d-\overline{R}_{n+1}(\overline{\eta}_s)^d)}{(\overline{R}_n(\overline{\eta}_s)-\overline{R}_{n+1}(\overline{\eta}_s))}\quad\text{where $\overline{R}_j(\overline{\eta}_s)=\sum_{n \colon n=j}^C\overline{\eta}_s(\mathcal{U}_n)$.} \]

In (5.1), the second term on the right-hand side is due to the increase of the ages of the progressing jobs linearly with time at unit rate. The third and fourth terms on the right-hand side of (5.1) are due to the departure and arrival of a job, respectively.

Remark 5.1. The t-continuity of the mapping $\overline{\eta}_t$ is equivalent to the continuity of the mapping $t\mapsto\langle \overline{\eta}_t,\phi\rangle$ for all $\phi\in\mathcal{C}_b^1(\mathcal{U})$ since $\mathcal{C}_b^1(\mathcal{U})$ is a separating class of $\mathcal{M}_1(\mathcal{U})$ [12, p. 111].

Although the mean field equation (5.1) is defined for the class of functions $\phi\in\mathcal{C}_b^1(\mathcal{U})$ , it is more useful to obtain an approximation of the process $(\langle\overline{\eta}_t^N,{I_{\{{B}\}}}\rangle, t\geq 0)$ for an open set $B\in\mathcal{B}(\mathcal{U})$ . Therefore we need to obtain the evolution equations of the real-valued process $(\langle\overline{\eta}_t,{I_{\{{B}\}}}\rangle, t\geq 0)$ . In this direction, we first obtain the evolution equations of the real-valued process $(\langle\overline{\eta}_t,\psi\rangle,t\geq 0)$ where $\psi\in\mathcal{C}_b(\mathcal{U})$ . We then proceed to obtain the evolution equations of the process $(\langle\overline{\eta}_t,{I_{\{{B}\}}}\rangle,t\geq 0)$ , where B is an open set, with the help of the monotone convergence theorem, since there exists a sequence of functions in $\mathcal{C}_b(\mathcal{U})$ that increase pointwise to ${I_{\{{B}\}}}$ .

Lemma 5.1. A process $(\nu_t\in\mathcal{M}_1(\mathcal{U}),t\geq 0)$ with continuity of the mapping $t\mapsto\nu_t$ satisfies the mean field equation (5.1) if and only if it satisfies the following equation for all $\phi\in\mathcal{C}_b(\mathcal{U})$ :

(5.2) \begin{align} \langle\nu_t,\phi\rangle & =\langle\nu_0,\tau_t\phi\rangle+\int_{r=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j)(\tau_{t-r}\phi(\textbf{\textit{x}}_n^{-j})-\tau_{t-r}\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \nu_r(\textbf{\textit{x}}_n)\notag \\ &\quad\, +\bigg[\nu_r(\{0\})\lambda\Phi_0(\nu_r)(\tau_{t-r}\phi(1,0)-\tau_{t-r}\phi(0))\notag\\ &\quad\, +\sum_{n=1}^{C-1}\sum_{j=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)}\lambda\Phi_n(\nu_r) (\tau_{t-r}\phi(\textbf{\textit{x}}_n^j;\,0)-\tau_{t-r}\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \nu_r(\textbf{\textit{x}}_n)\bigg]\bigg)\,{\mathrm{d}} r. \label{eqn6}\end{align}

The proof is given in Appendix B.

Using (5.2), we show that starting with an initial measure $\nu_0$ , for $t\geq 0$ , there exists a unique measure $\nu_t\in\mathcal{M}_1(\mathcal{U})$ that satisfies (5.1).

For any finite measure $\nu$ defined on $\mathcal{U}$ , the operator $\langle \nu,\phi\rangle$ is a continuous linear operator on the space of functions $\phi\in\mathcal{C}_b(\mathcal{U})$ , and let

\begin{equation*} \|{{\nu}}\|=\sup_{\phi\in\mathcal{C}_b(\mathcal{U})}\dfrac{{ \vert {\langle \nu,\phi\rangle} \vert}}{\|{{\phi}}\|}. \end{equation*}

Theorem 5.1. There exists a unique solution in $\mathcal{C}_{\mathcal{M}_1(\mathcal{U})}([0,\infty))$ to the mean field equations. In particular, if $(\nu_t^1,t\geq 0)$ and $(\nu_t^2,t\geq 0)$ are two mean field solutions starting at initial measures $\nu_0^1\in\mathcal{M}_1(\mathcal{U})$ , $\nu_0^2\in\mathcal{M}_1{(\mathcal{U})}$ , respectively, then

\begin{equation*} \|{{\nu_t^1-\nu_t^2}}\|\leq \|{{\nu_0^1-\nu_0^2}}\|\,{\mathrm{e}}^{(2C\|{{\beta}}\|+8d^2\lambda)t}. \end{equation*}

The proof is given in Appendix C.

We now show the convergence of the sequence of the processes $(\overline{\eta}_t^N,t\geq 0)$ . For this, we first assume the following.

Assumption 5.2. The sequence of the initial random measures $\{\overline{\eta}_0^N\}$ satisfies

\begin{equation*} (\overline{\eta}_0^N,\langle\overline{\eta}_0^N,\mathcal{I}\rangle)\Rightarrow (\vartheta,\langle\vartheta,\mathcal{I}\rangle), \end{equation*}

where $\vartheta\in\mathcal{M}_1(\mathcal{U})$ is a probability measure that possesses a density (with respect to Lebesgue measure) and $\langle\vartheta,\mathcal{I}\rangle<\infty$ .

Theorem 5.3. If the sequence of random measures $\{\overline{\eta}^N_0\}$ satisfies Assumption 5.2, then

\[ (\overline{\eta}^N_t,t\geq 0)\Rightarrow(\overline{\eta}_t,t\geq 0), \]

where $(\overline{\eta}_t,t\geq 0)$ is the unique solution to (5.1) with the initial point $\vartheta$ . The process $(\overline{\eta}_t,t\geq 0)$ is referred to as the mean field limit.

The proof is given in Section 6.

Remark 5.2. For any time t, a consequence of Theorem 5.3 is that as $N\to\infty$ , any finite set of servers are independent of each other. Furthermore, as $N\to\infty$ , $\overline{\eta}_t$ indicates the probability law of a server’s state at time t and the arrival process to a server is a Poisson process with rate $\lambda\Phi_n(\overline{\eta}_t)$ when there are n ( $n\geq 0$ ) progressing jobs. The proof follows from the same arguments as in the proof of Proposition 2 of [Reference Mukhopadhyay, Mazumdar and Guillemin33].

Lemma 5.2. For any time t, the measure $\overline{\eta}_t$ has a density function with respect to Lebesgue measure for almost all $\textbf{\textit{u}}_n\in\mathcal{U}_n$ , $n\geq 1$ .

The proof is given in Appendix F.

For any subset $B\in \mathcal{B}(\mathcal{U})$ , once $(\overline{\eta}^N_t,t\geq 0)\Rightarrow(\overline{\eta}_t,t\geq 0)$ , since $\overline{\eta}_t$ is absolutely continuous with respect to Lebesgue measure for every $t\geq 0$ , the continuous mapping theorem implies that

\[ (\langle \overline{\eta}_t^N,{I_{\{{B}\}}}\rangle,t\geq 0) \Rightarrow (\langle\overline{\eta}_t,{I_{\{{B}\}}}\rangle,t\geq 0). \]

This shows that for large N, we can approximate $\langle \overline{\eta}^N_t,{I_{\{{B}\}}}\rangle$ by $\langle\overline{\eta}_t,{I_{\{{B}\}}}\rangle$ .

5.3. Stationary regime

We now discuss the stationary behavior of the mean field.

We first show an analogy between the MFEs of the considered multi-server Erlang loss system under the SQ(d) routeing policy and the dynamics of another single-server Erlang loss system with state-dependent arrivals. We then exploit this analogy to prove the uniqueness of the fixed point of the mean field and its insensitivity. We first recall the dynamics of the probability measure of the server state of a single-server Erlang loss system with capacity C, where jobs arrive according to a Poisson process with prespecified state-dependent arrival rates.

Consider a single-server system with capacity C where jobs arrive according to a Poisson process at rate $\alpha_n$ when there are n progressing jobs in the system. The service times are generally distributed as stated in the system model of Section 2. Let $\nu_t^{(\mathrm{single})}$ be the probability measure of the server state at time t defined on $\mathcal{U}$ . For $\phi\in\mathcal{C}_b^1(\mathcal{U})$ , it can be verified that the Kolmogorov equations are given by

(5.3) \begin{align} \langle\nu_t^{(\mathrm{single})},\phi\rangle & =\langle\nu_0^{(\mathrm{single})},\phi\rangle+\int_{s=0}^t \langle\nu^{(\mathrm{single})}_s,\nabla_1\phi\rangle\,{\mathrm{d}} s \notag \\ &\quad\, -\int_{s=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j)(\phi(\textbf{\textit{x}}_n^{-j})-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \nu^{(\mathrm{single})}_s(\textbf{\textit{x}}_n) \notag \\ &\quad\, +\bigg[(\alpha_0\nu^{(\mathrm{single})}_s(\{0\})(\phi(1,0)-\phi(0))) +\sum_{n=1}^{C-1}\sum_{j=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)} \notag \\ &\quad\, \times\alpha_n (\phi(\textbf{\textit{x}}_n^j;\,0)-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \nu^{(\mathrm{single})}_s(\textbf{\textit{x}}_n)\bigg]\bigg) \,{\mathrm{d}} s. \label{eqn7}\end{align}

On comparing the mean field equation (5.1) with the Kolmogorov equation of a single-server system given by (5.3), it is clear that both the dynamics are similar except that $\alpha_n$ in (5.3) is replaced by $\lambda\Phi_n(\overline{\eta}_s)$ when the probability measure of the server state is $\overline{\eta}_s$ at time s. Equation (5.1) only differs from the equation with $\alpha_n$ in that the arrival rates depend on $\overline{\eta}_t$ . This is an example of a non-linear Markov process, which means that the generator of the Markov process at time t depends on the current distribution $\overline{\eta}_t$ of the Markov process [Reference Kolokoltsov25], while in equation (5.3) for fixed $(\alpha_i,0\leq i\leq C)$ denotes a Markov process whose generator does not depend on the current distribution.

We now study the fixed point of the mean field. Let $P_t(0)$ be equal to $\nu_t(\{0\})$ and let $p_t(\textbf{\textit{x}}_n)$ be the probability density of $\nu_t$ with respect to Lebesgue measure at $\textbf{\textit{x}}_n$ . We obtain the differential equations satisfied by the process $(P_t,t\geq 0)$ with $P_t=(P_t(\textbf{\textit{u}}),\textbf{\textit{u}}\in\mathcal{U})$ , where

\begin{equation*} P_t(\textbf{\textit{y}}_n)=\int_{x_1=0}^{y_1}\ldots\int_{x_n=0}^{y_n}p_t(\textbf{\textit{x}}_n)\,{\mathrm{d}} x_1\cdots {\mathrm{d}} x_n. \end{equation*}

Here, from Remark 5.2, since $\overline{\eta}_t$ is the distribution of a server’s state as $N\to\infty$ , it implies that $P_t(\textbf{\textit{y}}_n)$ is the probability that a server has n jobs and the ith job’s age is at most $y_i$ for $1\leq i\leq n$ as $N\to\infty$ . Also, since $\overline{\eta}_t^N({\cdot})\Rightarrow\overline{\eta}({\cdot})$ , for a large value of N, the fraction of servers with n jobs and the ith job’s age being at most $y_i$ for $1\leq i\leq n$ can be approximated by $P_t(\textbf{\textit{y}}_n)$ .

Lemma 5.3. The process $(P_t, t\geq 0)$ satisfies

(5.4) \begin{equation} \dfrac{{\mathrm{d}} P_t(0)}{{\mathrm{d}} t}=\int_{y=0}^{\infty}\beta(y)\bigg(\dfrac{\partial P_t(1,y)}{\partial y}\bigg)\,{\mathrm{d}} y-\lambda\Phi_0(P_t)P_t(0), \label{eqn8} \end{equation}

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

\begin{align*} %%\label{eq:mf_pdes2} \dfrac{{\mathrm{d}} P_t(\textbf{\textit{y}}_n)}{{\mathrm{d}} t} & =-\sum_{i=1}^n\dfrac{\partial P_t(\textbf{\textit{y}}_n)}{\partial y_i} +\sum_{j=1}^{n+1}\int_{x_j=0}^{\infty}\beta(x_j)\bigg(\dfrac{\partial P_t(\textbf{\textit{y}}_n^j;\,x_j)}{\partial x_j}\bigg)\,{\mathrm{d}} x_j \notag \\ &\quad\, -\sum_{j=1}^n\int_{x_j=0}^{y_j}\beta(x_j)\bigg(\dfrac{\partial P_t(\textbf{\textit{y}}_n^{-j};\,x_j)}{\partial x_j}\bigg)\,{\mathrm{d}} x_j \notag \\ &\quad\, +\sum_{j=1}^n\lambda\dfrac{\Phi_{n-1}(P_t)}{n}P_t(\textbf{\textit{y}}_n^{-j}) -\lambda\Phi_n(P_t) P_t(\textbf{\textit{y}}_{n}), \end{align*}

and for $n=C$ ,

(5.5) \begin{equation} \dfrac{{\mathrm{d}} P_t(\textbf{\textit{y}}_n)}{{\mathrm{d}} t} =-\sum_{i=1}^n\dfrac{\partial P_t(\textbf{\textit{y}}_n)}{\partial y_i} -\sum_{j=1}^n\int_{x_j=0}^{y_j}\beta(x_j)\bigg(\dfrac{\partial P_t(\textbf{\textit{y}}_n^{-j};\,x_j)}{\partial x_j}\bigg)\,{\mathrm{d}} x_j +\sum_{j=1}^n\lambda\dfrac{\Phi_{n-1}(P_t)}{n}P_t(\textbf{\textit{y}}_n^{-j}), \label{eqn9} \end{equation}

where

\[ \Phi_n(P_t)=\dfrac{(R^d_n(P_t)-R_{n+1}^d(P_t))}{(R_n(P_t)-R_{n+1}(P_t))} \quad\text{and}\quad R_n(P_t)=\sum_{j \colon j=n}^C\lim_{b\to\infty}P_t(j,b,\ldots,b). \]

The proof is given in Appendix E.

Remark 5.3. Specializing the results to the exponential case with mean $\frac{1}{\mu}$ , $\beta(x)=\mu$ , and denoting $Q_t(n) = \lim_{b\to\infty}P_t(n,b,\ldots,b)$ , it can be verified that the process

\[ (Q_t,t\geq 0)=(Q_t(n),0\leq n\leq C,t\geq 0) \]

is the unique solution of the mean field equations given in [Reference Mukhopadhyay, Mazumdar and Guillemin33] for the case of the exponential distributions with rate $\mu=1$ .

We next state the principal result on the insensitivity of the fixed point of the MFEs. The proof is given in Section 7.

Theorem 5.4. The process $(P_t, t\geq 0)=(P_t(\textbf{\textit{u}}),\textbf{\textit{u}}\in\mathcal{U},t\geq 0)$ has a unique fixed point given by $\boldsymbol{\pi}=(\pi(\textbf{\textit{y}}),\textbf{\textit{y}}\in\mathcal{U})$ , where

\begin{equation*} %%\label{eq:mf_fixed_pt} \pi(\textbf{\textit{y}}_n)=\pi^{(\exp)}_n\mu^n\prod_{i=1}^n\int_{x_i=0}^{y_i}\overline{G}(x_i)\,{\mathrm{d}} x_i \end{equation*}

and $\boldsymbol{\pi}^{(\exp)}=(\pi^{(\exp)}_n,0\leq n\leq C)$ denotes the unique fixed point of the mean field when the service times are exponentially distributed with the mean $\frac{1}{\mu}$ and $\pi^{(\exp)}_n$ is the stationary probability that there are n jobs in the limiting system. Since

\[ \int_{x=0}^{\infty}\overline{G}(x)\,{\mathrm{d}} x=\frac{1}{\mu}, \]

the fixed point of the mean field is insensitive, that is,

\begin{equation*} \lim_{b\to\infty}\pi(n,b,\ldots,b)=\pi^{(\exp)}_n. \end{equation*}

6. Convergence of the normalized processes: proof of Theorem 5.3

By using the results on the construction of a martingale in Appendix D, we now show that the normalized process $(\overline{\eta}_t^N,t\geq 0)$ converges to the mean field limit.

Let $(\overline{\mathcal{F}}_t^N,t\geq 0)$ be the right-continuous filtration associated with the process $(\overline{\eta}_t^N,t\geq 0)$ . Note that we have $ (\overline{\eta}_t^N,t\geq 0)\in\mathcal{D}_{\mathcal{M}_1^N(\mathcal{U})}([0,\infty)). $ We first show that the sequence of processes $(\overline{\eta}_t^N,t\geq 0)$ is relatively compact, and we then prove that every limit point $(\boldsymbol{\chi}_t,t\geq 0)$ almost surely has continuous sample paths with respect to t and they coincide with the unique mean field solution with the initial point $\vartheta$ . For every limit point $(\boldsymbol{\chi}_t,t\geq 0)$ , $\boldsymbol{\chi}_0$ almost surely coincides with the measure $\vartheta$ from Assumption 5.2. Further, we find that the mean field solution is unique for the given initial measure. Hence, we conclude that for all the limit points, almost surely sample paths coincide with the unique mean field solution $(\overline{\eta}_t,t\geq 0)$ with the initial point $\vartheta$ . The process $(\overline{\eta}_t,t\geq 0)$ is referred to as the mean field limit. Therefore $(\overline{\eta}_t^N,t\geq 0)$ converges in distribution to the mean field limit.

For $\phi\in\mathcal{C}_b^1(\mathcal{U})$ , from Proposition D.1, the process $(\overline{M}_t^N(\phi),t\geq 0)$ defined as follows is an RCLL square-integrable $\overline{\mathcal{F}}_t^N$ -martingale:

(6.1) \begin{align} \overline{M}_t^N(\phi) & =\langle\overline{\eta}_t^N,\phi\rangle-\langle\overline{\eta}_0^N,\phi\rangle-\int_{s=0}^t \langle\overline{\eta}_s^N,\nabla_1\phi\rangle\,{\mathrm{d}} s -\int_{s=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j)\notag \\ &\quad\, \times(\phi(\textbf{\textit{x}}_n^{-j})-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \overline{\eta}_s^N(\textbf{\textit{x}}_n)\notag \\ &\quad\, +\bigg[\overline{\eta}_s^N(\{0\})\lambda\Phi_0(\overline{\eta}_s^N)(\phi(1,0)-\phi(0)) +\sum_{n=1}^{C-1}\sum_{j=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)} \notag \\ &\quad\, \times\lambda\Phi_n(\overline{\eta}_s^N) (\phi(\textbf{\textit{x}}_n^j;\,0)-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \overline{\eta}_s^N(\textbf{\textit{x}}_n)\bigg]\bigg) \,{\mathrm{d}} s. \label{eqn10}\end{align}

We further have

(6.2) \begin{align} \langle \overline{M}_\cdot^N(\phi)\rangle_t & = \dfrac{1}{N}\bigg[\int_{s=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j)(\phi(\textbf{\textit{x}}_n^{-j})-\phi(\textbf{\textit{x}}_n))^2 \,{\mathrm{d}} \overline{\eta}_s^N(\textbf{\textit{x}}_n) \notag \\ &\quad\, +\bigg[\overline{\eta}_s^N(\{0\})\lambda\Phi_0(\overline{\eta}_s^N)(\phi(1,0)-\phi(0))^2 \notag \\ &\quad\, +\sum_{n=1}^{C-1}\sum_{j=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)}\lambda\Phi_n(\overline{\eta}_s^N) (\phi(\textbf{\textit{x}}_n^j;\,0)-\phi(\textbf{\textit{x}}_n))^2 \,{\mathrm{d}} \overline{\eta}_s^N(\textbf{\textit{x}}_n)\bigg]\bigg) \,{\mathrm{d}} s\bigg]. \label{eqn11}\end{align}

Since the space $\mathcal{D}_{\mathcal{M}_1(\mathcal{U})}([0,\infty))$ endowed with the Skorokhod topology is complete and separable, by using Prokhorov’s theorem [Reference Billingsley3], establishing the relative compactness of the sequence of the processes $\{(\overline{\eta}_t^N,t\geq 0)\}$ is equivalent to proving the tightness of the processes $\{(\overline{\eta}_t^N,t\geq 0)\}$ . From Theorem 4.6 of [Reference Jakubowski20], Jakubowski’s criteria, recalled below, can be used to establish the relative compactness of the sequence of the processes $\{(\overline{\eta}_t^N,t\geq 0)\}$ .

Jakubowski’s criteria. A sequence of $\{X^N\}$ of $\mathcal{D}_{\mathcal{M}_1(\mathcal{U})}([0,\infty))$ -valued random elements defined on $(\Omega,\mathbb{F},\mathbb{P})$ is tight if and only if the following two conditions are satisfied.

J1 For each $T>0$ and $\gamma>0$ , there exists a compact set $\mathbb{K}_{T,\gamma}\subset\mathcal{M}_1(\mathcal{U})$ such that

\begin{equation*} %%\label{eq:J1} \liminf_{N\to\infty}\mathbb{P}(X^N_t\in\mathbb{K}_{T,\gamma} \quad \text{{for all} } t\in[0,T])>1-\gamma. \end{equation*}

This condition is called the compact-containment condition.

J2 There exists a family $\mathcal{Q}$ of real-valued continuous functions F defined on $\mathcal{M}_1(\mathcal{U})$ that separates points in $\mathcal{M}_1(\mathcal{U})$ and is closed under addition such that, for every $F\in \mathcal{Q}$ , the sequence $\{(F(X_t^N),t\geq 0)\}$ is tight in $\mathcal{D}_{\mathbb{R}}([0,\infty))$ .

To prove condition J2, we define a class of functions $\mathcal{Q}$ as follows:

\begin{equation*} %%\label{def:classQ} \mathcal{Q}\triangleq \{F \colon \text{{there exists} } f\in\mathcal{C}_b^1(\mathcal{U}) \text{ such that } F(\nu)=\langle\nu,f\rangle,\ \text{{for all} } \nu\in\mathcal{M}_1(\mathcal{U}) \}. \end{equation*}

Clearly every function $F\in\mathcal{Q}$ is continuous with respect to the weak topology on $\mathcal{M}_1(\mathcal{U})$ , and further the class of functions $\mathcal{Q}$ separates points in $\mathcal{M}_1(\mathcal{U})$ and is also closed under addition. We next recall the following result (from [Reference Robert35, Theorem C.9]) to prove condition J2.

Tightness in $\mathcal{D}_{\mathbb{R}}([0,T])$ . If $S=\mathcal{D}_{\mathbb{R}}([0,T])$ and $\{\mathbb{P}_n\}$ is a sequence of probability distributions on S, then $\{\mathbb{P}_n\}$ is tight if, for any $\epsilon>0$ ,

C1 There exists $b>0$ such that

\begin{equation*} \mathbb{P}_n({ \vert {X(0)} \vert}>b)\leq\epsilon \end{equation*}

for all $n\in\mathbb{Z}_+$ .

C2 For any $\gamma>0$ , there exists $\rho>0$ such that

\begin{equation*} \mathbb{P}_n(w_X(\rho)>\gamma)\leq\epsilon \end{equation*}

for n sufficiently large, where

\begin{equation*} w_X(\rho)=\sup\{{ \vert {X(t)-X(s)} \vert} \colon s,t\leq T, { \vert {s-t} \vert}\leq\rho\} \end{equation*}

and any limiting point $\mathbb{P}$ satisfies $\mathbb{P}(\mathcal{C}_{\mathbb{R}}([0,T]))=1$ .

We first establish the relative compactness of the sequence $\{(\overline{\eta}^N_t,t\geq 0)\}$ . For this, we next prove conditions C1 and C2 that are sufficient to prove the relative compactness of the sequence $\{(\langle\overline{\eta}^N_t,\phi\rangle,t\geq 0)\}$ for $\phi\in\mathcal{C}_b^1(\mathcal{U})$ in $D_{\mathbb{R}}([0,\infty))$ . For any $T>0$ , $t\in[0,T]$ , we have $ \langle\overline{\eta}_t^N,\phi\rangle\leq \|{{\phi}}\|_1 \langle\overline{\eta}_t^N,\mathbf{1}\rangle $ , and since $\langle\overline{\eta}_t^N,\mathbf{1}\rangle=1$ , condition C1 is trivially satisfied with $b=\|{{\phi}}\|_1$ .

We next prove that condition C2 holds. For $\epsilon>0$ , by using (6.2) and Doob’s inequality [Reference Ethier and Kurtz12, p. 63], we have

\begin{equation*} \mathbb{P}\bigg(\sup_{t\leq T}{ \vert {\overline{M}_t^N(\phi)} \vert}\geq \epsilon\bigg) \leq \dfrac{4}{\epsilon^2}{\mathbb{E} [{\langle \overline{M}^N_{\cdot}(\phi) \rangle_T} ]} \leq 4T\|{{\phi}}\|^2\dfrac{1}{N}(\|{{\beta}}\|+ {\mathrm{d}}\lambda){,} \end{equation*}

and hence $\mathbb{P} (\sup_{t\leq T}{\vert {\overline{M}_t^N(\phi)} \vert}\geq \epsilon)\to 0$ as $N\to\infty$ . Therefore the sequence of processes $\{(\overline{M}_t^N(\phi),t\geq 0)\}$ converges in distribution to the null process from the standard convergence criterion in $\mathcal{D}_{\mathbb{R}}([0,T])$ . Further, the sequence of processes $\{(\overline{M}_t^N(\phi),t\geq 0)\}$ is tight in $\mathcal{D}_{\mathbb{R}}([0,T])$ and hence there exists $\rho'>0$ and $N'>0$ such that, for all $N\geq N'$ , we have

(6.3) \begin{equation} \mathbb{P}\bigg(\sup_{u,v\leq T,{ \vert {u-v} \vert}\leq\rho'}{ \vert {\overline{M}^N_v(\phi)-\overline{M}^N_u(\phi)} \vert}\geq \dfrac{\gamma}{2}\bigg)\leq \dfrac{\epsilon}{2} {.} \label{eqn12} \end{equation}

For any $u<v\leq T$ , from (6.1) we have

\begin{align*} { \vert {\langle\overline{\eta}_v^N,\phi\rangle-\langle\overline{\eta}_u^N,\phi\rangle} \vert}& \leq \int_{s=u}^v{ \vert {\langle\overline{\eta}_s^N,\nabla_1\phi\rangle} \vert} \,{\mathrm{d}} s+2\|{{\beta}}\|\|{{\phi}}\|C{ \vert {u-v} \vert}+2\|{{\phi}}\|\lambda{\vert {u-v} \vert}\\ &\quad\, +{\vert {\overline{M}_v^N(\phi)-\overline{M}_u^N(\phi)} \vert}. \end{align*}

Further, we can write

(6.4) \begin{equation} { \vert {\langle\overline{\eta}_v^N,\phi\rangle-\langle\overline{\eta}_u^N,\phi\rangle} \vert}\leq { \vert {v-u} \vert}C\|{{\phi}}\|_1(1+2\|{{\beta}}\|+2\,{\mathrm{d}}\lambda)+{ \vert {\overline{M}_v^N(\phi)-\overline{M}_u^N(\phi)} \vert}. \label{eqn13} \end{equation}

Therefore, by using equations (6.3) and (6.4), there exist $\rho>0$ and $N_1>0$ such that, for $N\geq N_1$ , we have

\[ %%\label{eq:process_tight} \mathbb{P}\bigg(\sup_{u,v\leq T,{\vert {u-v} \vert}\leq\rho}{ \vert {\langle\overline{\eta}^N_v,\phi\rangle-\langle\overline{\eta}^N_u,\phi\rangle} \vert}\geq \gamma\bigg)\leq \epsilon. \]

This proves condition C2. Since conditions C1 and C2 hold, condition J2 also holds.

We next prove the compact containment condition J1. Let $(n_{i}(t),x_{i1}(t),\ldots,x_{in_i(t)}(t))$ be the state of the ith server at time t where $x_{ij}(t)$ denotes the age of the jth job at the ith server. Clearly, we have

\[ \langle\overline{\eta}_t^N,\mathcal{I}\rangle=\dfrac{1}{N}\sum_{i=1,n_i(t)>0}^N(x_{i1}(t)+\cdots+x_{in_i(t)}(t)). \]

We can classify the progressing jobs into two classes. The jobs that are in service from the beginning ( $t=0$ ) form the first class and the second class of jobs comprises the ones that entered the system in the interval (0,t]. At a server, the number of progressing jobs that belong to each class are upper-bounded by C. Let $Y_t$ be a random variable representing the age of a job belonging to the second class that is in progress at time t, and let Y be a random variable with job length distribution G. Then, for any $b\geq 0$ , we have

(6.5) \begin{equation} \mathbb{P}(Y_t\geq b)\leq \mathbb{P}(Y\geq b). \label{eqn14} \end{equation}

Therefore, using (6.5), since each server has capacity C, for any time $t\geq 0$ we can write

\begin{equation*} \mathbb{P}(\langle\overline{\eta}_t^N,\mathcal{I}\rangle\geq b)\leq \mathbb{P}\bigg(\langle\overline{\eta}_0^N,\mathcal{I}\rangle+Ct+\dfrac{1}{N}\sum_{i=1}^N(Y_{i1}+\cdots+Y_{iC})\geq b\bigg), \end{equation*}

where $(Y_{ij},1\leq i\leq N,1\leq j\leq C)$ are i.i.d. random variables with distribution G. Further, by the weak law of large numbers, we have

\[ \dfrac{1}{N}\sum_{i=1}^N(Y_{i1}+\cdots+Y_{iC})\Rightarrow \dfrac{C}{\mu} \quad\text{as $N\to\infty$.} \]

Therefore, by choosing $Z_T=2\langle\vartheta,\mathcal{I}\rangle+2CT+\frac{2C}{\mu}$ , we have

(6.6) \begin{equation} \mathbb{P}\bigg(\sup_{t\in[0,T]}\langle\overline{\eta}_t^N,\mathcal{I}\rangle \,> Z_T\bigg)\to 0 \label{eqn15} \end{equation}

as $N\to\infty$ . Let us define

\[ \mathcal{L}_{T}\triangleq \{\zeta\in\mathcal{M}_1(\mathcal{U}) \colon \langle\zeta,\mathcal{I}\rangle\leq Z_T\}. \]

Since $\langle\zeta,\mathcal{I}\rangle\leq Z_T$ for $\zeta\in\mathcal{L}_{T}$ , let

\[ B=\mathcal{U}_0\cup(\cup_{n\geq 1}\{(n,y_1,\ldots,y_n) \colon 0\leq y_i\leq r, 1\leq i\leq n\}) \]

and let $\overline{B}$ be the complement of B. Then we have $ \zeta(\overline{B})\leq \frac{Z_T}{ r}. $ Hence, $ \lim_{r\to\infty}\sup_{\zeta\in\mathcal{L}_{T}}\zeta(\overline{B})=0. $ Therefore, from Lemma A7.5 of [Reference Kallenberg21], $\mathcal{L}_{T}$ is relatively compact in $\mathcal{M}_1(\mathcal{U})$ . Further, from (6.6), we have

\[ \liminf_{N\to\infty}\mathbb{P}(\overline{\eta}_t^N\in \mathcal{L}_{T}, \ \text{{for all} } t\in[0,T])> 1-\gamma. \]

Let $\mathbb{K}_{T}$ be the closure of $\mathcal{L}_{T}$ . Then we have a compact set $\mathbb{K}_{T}$ such that

\[ \liminf_{N\to\infty}\mathbb{P}(\overline{\eta}_t^N\in \mathbb{K}_{T}, \ \text{{for all} } t\in[0,T])> 1-\gamma\quad\text{for all $0<\gamma<1$.} \]

This establishes condition J1 and hence the proof of the tightness of the sequence of processes $(\overline{\eta}^N_t,t\geq 0)$ is completed.

Let $(\boldsymbol{\chi}_t,t\geq 0)$ be a limit of a converging subsequence $\{(\overline{\eta}_t^{N_{i_k}},t\geq 0)\}$ . From condition C2, $\chi_t$ is continuous in t, $\mathbf{P}_{\chi}$ -almost surely, where $\mathbf{P}_{\chi}$ is the probability law of $(\chi_t,t\geq 0)$ . Furthermore, from [Reference Jakubowski20, Theorem 1.7], for $f\in\mathcal{C}_b(\mathcal{U}),\nu\in\mathcal{M}_1(\mathcal{U})$ , it follows that, for any $T>0$ , $(\nu_t,0\leq t\leq T)\mapsto(\langle\nu_t,f\rangle,0\leq t\leq T)$ is continuous in the Skorokhod topology. Then, since the martingale $(\overline{M}^{N_{i_k}}_t(\phi), t\geq 0)$ converges to the null process, by the continuous mapping theorem, we conclude that

\begin{align*} %%\label{eq:mf_limit_point} \langle\chi_t,\phi\rangle & =\langle\chi_0,\phi\rangle+\int_{s=0}^t \langle\chi_s,\nabla_1\phi\rangle\,{\mathrm{d}} s \notag \\ &\quad\, -\int_{s=0}^t\bigg(\sum_{n=1}^C\sum_{j=1}^n\idotsint_{\mathcal{U}_n}\beta(x_j) (\phi(\textbf{\textit{x}}_n^{-j})-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \chi_s(\textbf{\textit{x}}_n) \notag \\ &\quad\, +\bigg[(\chi_s(\{0\})\lambda\Phi_0(\chi_s)(\phi(1,0)-\phi(0))) +\sum_{n=1}^{C-1}\sum_{j=1}^{n+1}\idotsint_{\mathcal{U}_n}\dfrac{1}{(n+1)} \notag \\ &\quad\, \times\lambda\Phi_n(\chi_s) (\phi(\textbf{\textit{x}}_n^j;\,0)-\phi(\textbf{\textit{x}}_n))\,{\mathrm{d}} \chi_s(\textbf{\textit{x}}_n)\bigg]\bigg) \,{\mathrm{d}} s. \end{align*}

From Assumption 5.2, $\boldsymbol{\chi}_0= \vartheta$ almost surely and hence the sample paths coincide almost surely with the unique mean field solution with the initial point $\vartheta$ . This argument holds for every limit point, and hence the sample paths of every limit point are almost surely the same as the deterministic mean field solution with the initial point $\vartheta$ . This completes the proof.

7. Insensitivity: proof of Theorem 5.4

We now show that $\boldsymbol{\pi}=(\pi(\textbf{\textit{u}}),\textbf{\textit{u}}\in\mathcal{U})$ is the unique fixed point of the mean field. From [Reference Mukhopadhyay, Mazumdar and Guillemin33], we first recall that under the assumption of exponential service time distributions, there exists a unique probability measure of occupancy $\boldsymbol{\pi}^{(\exp)}=(\pi_n^{(\exp)},0\leq n\leq C)$ on $\{0,1,\ldots,C\}$ to the stationary MFEs given below,

(7.1) \begin{equation} \lambda^{(\exp)}_n(\boldsymbol{\pi}^{(\exp)})\pi_n^{(\exp)}=(n+1)\mu\pi_{n+1}^{(\exp)}, \label{eqn16} \end{equation}

where

(7.2) \begin{equation} \lambda^{(\exp)}_n(\boldsymbol{\pi}^{(\exp)})\!=\!\lambda\bgfrac{\bigg(\!\bigg(\sum_{j=n}^C\pi_j^{(\exp)}\bigg)^d\!-\!\bigg(\sum_{j=n+1}^C\pi_j^{(\exp)}\bigg)^d\bigg)}{\bigg(\!\bigg(\sum_{j=n}^C\pi_j^{(\exp)}\bigg)\!-\!\bigg(\sum_{j=n+1}^C\pi_j^{(\exp)}\bigg)\!\bigg)}. \label{eqn17} \end{equation}

Let $\mathbb{\theta}=(\theta(\textbf{\textit{u}}),\textbf{\textit{u}}\in\mathcal{U})$ be a fixed point of the MFEs of the process $(P_t,t\geq 0)$ under general service time distributions. Using $\boldsymbol{\theta}$ , let the corresponding probability measure of occupancy be $\boldsymbol{\Gamma}=(\Gamma_n,0\leq n\leq C)$ defined such that $\Gamma_n=\lim_{b\to\infty}\theta(n,b,\ldots,b))$ and $\Gamma_0=\theta(0)$ . We now show that

(7.3) \begin{equation} \theta(\textbf{\textit{y}}_n)= \bgfrac{\bigg(\prod_{i=1}^n\dfrac{\lambda^{(\exp)}_{i-1}(\boldsymbol{\Gamma})}{i\mu}\bigg)}{\bigg(1+\sum_{m=1}^C\bigg(\prod_{i=1}^m\dfrac{\lambda^{(\exp)}_{i-1}(\boldsymbol{\Gamma})}{i\mu}\bigg)\bigg)}\mu^n\prod_{i=1}^n\int_{x_i=0}^{y_i}\overline{G}(x_i)\,{\mathrm{d}} x_i \label{eqn18} \end{equation}

and

(7.4) \begin{equation} \theta(0)=\bgfrac{1}{\bigg(1+\sum_{m=1}^C\bigg(\prod_{i=1}^m\dfrac{\lambda^{(\exp)}_{i-1}(\boldsymbol{\Gamma})}{i\mu}\bigg)\bigg)}. \label{eqn19} \end{equation}

Then it implies that $\boldsymbol{\Gamma}$ also satisfies equations (7.1)–(7.2), and hence $\boldsymbol{\Gamma}=\boldsymbol{\pi}^{(\exp)}$ , concluding the insensitivity of the fixed point. Furthermore, we have

\[ \theta(\textbf{\textit{y}}_n)=\pi^{(\exp)}_n\mu^n\prod_{i=1}^n\int_{x_i=0}^{y_i}\overline{G}(x_i)\,{\mathrm{d}} x_i {,} \]

concluding the uniqueness of the fixed point of the mean field under general service time distributions.

To complete the proof, it remains to show the validity of equations (7.3)–(7.4). We now recall the stationary distribution $\boldsymbol{\pi}^{(\mathrm{single})}=(\pi^{(\mathrm{single})}(\textbf{\textit{u}}),\textbf{\textit{u}}\in\mathcal{U})$ of a single-server loss system with state-dependent Poisson arrival process with rate $\alpha_n$ ( $0\leq n\leq C$ ) when there are n progressing jobs and the service time distributions are as in the system model of Section 2. Then from [Reference Brumelle8], the stationary probability that the server has n progressing jobs and the ith job has age at most $y_i$ ( $1\leq i\leq n$ ) is given by

(7.5) \begin{equation} \pi^{(\mathrm{single})}(\textbf{\textit{y}}_n)= \bgfrac{\bigg(\prod_{i=1}^n\dfrac{\alpha_{i-1}}{i\mu}\bigg)}{\bigg(1+\sum_{m=1}^C\bigg(\prod_{i=1}^m\dfrac{\alpha_{i-1}}{i\mu}\bigg)\bigg)}\mu^n\prod_{i=1}^n\int_{x_i=0}^{y_i}\overline{G}(x_i)\,{\mathrm{d}} x_i \label{eqn20} \end{equation}

and

(7.6) \begin{equation} \pi^{(\mathrm{single})}(0)=\bgfrac{1}{\bigg(1+\sum_{m=1}^C\bigg(\prod_{i=1}^m\dfrac{\alpha_{i-1}}{i\mu}\bigg)\bigg)}. \label{eqn21} \end{equation}

For the given fixed point $\boldsymbol{\theta}$ of the mean field and its corresponding occupancy probability measure $\boldsymbol{\Gamma}$ , consider a single-server system under the assumption of a Poisson arrival process with state-dependent rate $\lambda_n^{(\exp)}(\boldsymbol{\Gamma})$ ( $0\leq n\leq C$ ) when there are n progressing jobs. Then the unique stationary distribution is given by equations (7.5)–(7.6) with $\alpha_n$ replaced by $\lambda_n^{(\exp)}(\boldsymbol{\Gamma})$ for all $0\leq n\leq C$ . But from equations (5.1), (5.3), and Lemma 5.3, since $\overline{R}_n(\boldsymbol{\theta})=\sum_{j=n}^C\Gamma_j$ , we find that $\boldsymbol{\theta}$ is also an another stationary distribution for the single-server system with state-dependent Poisson arrival process having rates $\lambda_n^{(\exp)}(\boldsymbol{\Gamma})$ for all $0\leq n\leq C$ . Since the stationary distribution must be unique, equations (7.3)–(7.4) must hold. This completes the proof.

8. Numerical results

Showing that the fixed point of the mean field approximates the stationary distribution of the system with large N remains an open problem. If one can establish that the equilibrium or fixed point of the MFEs is globally asymptotically stable (GAS), then the conclusion of the interchange of limits would follow from Prokhorov’s theorem [Reference Billingsley3]. Proving that the equilibrium point of the MFEs is GAS is a challenging problem because the joint distribution of the occupancy and ages does not possess any monotonicity properties, unlike the case of exponential service time distributions [Reference Mukhopadhyay, Mazumdar and Guillemin33]. In this section, we present numerical results on the validity of the GAS of the mean field for the case in which the service time distributions are mixed-Erlang. In this case, the state of a server is also multi-dimensional and the mean field is also not monotonic, unlike the exponential case. It is numerically easier to solve the MFEs for the case of mixed-Erlang distributions as they are systems of ordinary differential equations, unlike the case of general service time distributions for which the MFEs are partial differential equations, as we have shown. One more reason for using mixed-Erlang distributions is that such distributions are dense in the set of all distributions that have support on $\mathbb{R}_+$ ; see [Reference Asmussen2]. Our numerical results show that the mean field is GAS for the case of mixed-Erlang service time distributions.

We consider the system parameters as follows. The capacity of a server is assumed to be $C=5$ . The average job length is assumed to be equal to one, i.e. $\mu=1$ . The service times have a mixed-Erlang distribution given by sums of independent exponentially distributed random variables (known as an Erlang distribution) where the number of exponential phases (or independent random exponentials) is equal to $i\in \{1,2,\ldots,M\}$ with probability $p_i$ such that $\sum_{i=1}^Mp_i=1$ . Each exponential phase is assumed to have rate $\mu_p$ . Therefore we have

\begin{equation*}\dfrac{1}{\mu}=\dfrac{\sum_{i=1}^Mip_i}{\mu_p}.\label{eq:mean}\end{equation*}

We choose $M=3$ , $p_1=0.3,p_2=0.3,p_3=0.4$ .

Under mixed-Erlang service time distribution assumptions, let S be the set of all possible server states defined as $ S=\cup_{n=0}^CS_n $ where $S_0=\{(0)\}$ and

\[ S_n=\{(n,r_1,\ldots,r_n) \colon 1\leq r_i\leq M, 1\leq i\leq n\}. \]

We refer to an element in the set S as $\textbf{\textit{r}}$ and an element in the set $S_n$ as $\textbf{\textit{r}}_n$ . The system dynamics can be modeled as a Markov process $\textbf{\textit{x}}^N(t)=(x^N_{\textbf{\textit{r}}}(t),\textbf{\textit{r}}\in S)$ , where $x^N_{\textbf{\textit{r}}}(t)$ denotes the fraction of servers with n jobs such that the ith job has $r_i$ remaining phases at time t. Since the Markov process $(\textbf{\textit{x}}^N(t),t\geq 0)$ is defined on a finite-dimensional space, we can establish the mean field limit $\textbf{\textit{x}}(t)=(x_{\textbf{\textit{r}}}(t),\textbf{\textit{r}}\in S)$ by using the same procedure as that of the exponential service times case in [Reference Mukhopadhyay, Mazumdar and Guillemin33]. Hence we recall the following result without proof from [Reference Vasantam, Mukhopadhyay and Mazumdar41].

Proposition 8.1. If $\textbf{\textit{x}}^N(0)$ converges in distribution to a state $\textbf{\textit{u}}$ , then the process $\textbf{\textit{x}}^N({\cdot})$ converges in distribution to a deterministic process $\textbf{\textit{x}}(\cdot,\textbf{\textit{u}})$ as $N\to\infty$ called the mean field. The process $\textbf{\textit{x}}(\cdot,\textbf{\textit{u}})$ is the unique solution of the following system of differential equations,

\begin{align*} %%\label{eq:fmean field} \textbf{\textit{x}}(0,\textbf{\textit{u}})&=\textbf{\textit{u}},\\ {{\dot{x}_{\textbf{\textit{r}}_n}}}(t,\textbf{\textit{u}})&={h_{\textbf{\textit{r}}_n}}(\textbf{\textit{x}}(t,\textbf{\textit{u}})), \end{align*}

and $\textbf{\textit{h}}=(h_{\textbf{\textit{r}}},\textbf{\textit{r}}\in S)$ with the mapping $h_{\textbf{\textit{r}}_n}$ given by

\begin{align*} %%\label{eq:smean field} h_{\textbf{\textit{r}}_n}(\textbf{\textit{x}}) & =\sum_{b=1}^{n}\bigg(\dfrac{ p_{r_b}}{n}\bigg)x_{(\textbf{\textit{r}}_n^{-j})}\lambda^{\mathrm{(ME)}}_{n-1}(\textbf{\textit{x}}) - x_{\textbf{\textit{r}}_n}\lambda^{\mathrm{(ME)}}_{n}(\textbf{\textit{x}}) {I_{\{{{n<C}}\}}} \notag \\ &\quad\, +\sum_{b=1}^{n+1}\mu_p{I_{\{{n<C}\}}}x_{(\textbf{\textit{r}}_n^b;\,1)} +\sum_{b=1}^{n}\mu_px_{(n,r_1,\ldots,r_{b-1},r_b+1,r_{b+1},\ldots,r_{n})}-n\mu_px_{\textbf{\textit{r}}_n}, \end{align*}

where

\begin{equation*} %%\label{eq:arrival_rate} \lambda^{\mathrm{(ME)}}_{n}(\textbf{\textit{u}})=\dfrac{\lambda}{(\sum_{\textbf{\textit{r}}_n}u_{\textbf{\textit{r}}_n})}\bigg[ \bigg(\sum_{i=n}^C\sum_{\textbf{\textit{b}}_i}u_{\textbf{\textit{b}}_i}\bigg)^d-\bigg(\sum_{i=n+1}^C\sum_{\textbf{\textit{b}}_i}u_{\textbf{\textit{b}}_i}\bigg)^d \bigg]. \end{equation*}

In Figure 1, we plot $d_{E}^2(\textbf{\textit{x}}(t,\textbf{\textit{u}}),\mathbf{\pi})$ as a function of t, where $d_{E}$ is the Euclidean distance defined by

\begin{equation*} d_{E}(\textbf{\textit{u}},\textbf{\textit{v}})=\sqrt{\sum_{\mathbf{l}\in S}{ \vert {u_{\mathbf{l}}-v_{\mathbf{l}}} \vert}^2}. \end{equation*}

Fig. 1.

Convergence of mean field to the fixed point.

It is observed that for $d=2$ , $\lambda=1$ , and for four different initial points $\textbf{\textit{u}}_1,\textbf{\textit{u}}_2,\textbf{\textit{u}}_3$ , and $\textbf{\textit{u}}_4$ , the mean field $\textbf{\textit{x}}(t,\textbf{\textit{u}})$ for mixed-Erlang service time distribution converges to its unique fixed point $\boldsymbol{\pi}$ . Note that the computed $\boldsymbol{\pi}$ depends on the chosen value of d. This supports that $\boldsymbol{\pi}$ is globally stable.

We conclude with some numerical results for the blocking probability of the above system showing closeness to the theoretical lower bound. Under asymptotic independence, any finite set of servers is independent and the fixed point of the mean field implies that the fixed point is the stationary distribution of the state of a server. The average blocking probability is then given by $\pi_q(C)^d$ , where $\pi_q(C)=\lim_{b\to\infty}\pi(C,b,\ldots,b)$ . Let us recall the lower bound on the average blocking probability, denoted by $P^{\mathrm{avg}}_{\mathrm{block}}$ for any routeing scheme shown in [Reference Mukhopadhyayay, Karthik, Mazumdar and Guillemin32]. From Little’s law, the average number of customers in the system is equal to $(1-P^{\mathrm{avg}}_{\mathrm{block}})N\lambda$ which is upper-bounded by NC. Hence,

\begin{equation*} P^{\mathrm{avg}}_{\mathrm{block}}\geq \bigg(1-\dfrac{C}{\lambda}\bigg)_+, \end{equation*}

where $(x)_+=\max{(x,0)}$ . In Figure 2, we plot the lower bound $(1-\frac{C}{\lambda})_+$ and the average blocking probability under the SQ(d) routeing, and the state-independent random routeing where a destination server is chosen uniformly at random, as a function of $\lambda$ . It is clear that the resulting average blocking probability under the SQ(d) policy is much lower than the resulting average blocking probability when pure random routeing is employed. Furthermore, the average blocking probability under the SQ(d) routeing approaches the lower bound as d increases.

Fig. 2.

Comparison of the average blocking probability under SQ(d) with lower bound.

9. Concluding remarks

In this paper we have provided a measure-valued process approach to establish the mean field behavior of loss systems with SQ(d) routeing and general service time requirements. The extension of these results to multi-class systems, where servers are classified into different classes based on their capacities and jobs are classified into different classes based on their service requirements, follows in a similar manner, mutatis mutandis, from the approach used here. Establishing the global asymptotic stability of the unique fixed point remains an open problem.