2.1. The associated dynamical system

We introduce some notation used throughout this section. The state space of a node is given by ${\mathcal X}$ , defined as follows:

(4) \begin{equation} \begin{cases} {\mathcal X} = \{z = (x,y)\in{\mathbb N}^2\,|\, x + y{\leq} C\},\\[5pt] {\mathcal X}_+ = \{z = (x,y)\in{\mathbb N}^2\,|\, x + y < C\}, \text{ and }{\mathcal X}_+^c = {\mathcal X}{\setminus}{\mathcal X}_+ = \{z = (x,y)\in{\mathbb N}^2\,|\, x + y = C\} . \end{cases}\end{equation}

If $z = (x,y)\in{\mathcal X}$ , x (resp. y) is the number of class 1 (resp. class 2) customers.

The set of probability distributions on ${\mathcal X}$ is denoted by ${\mathcal P}({\mathcal X})$ , and ${\mathcal C}({\mathbb R}_+,{\mathcal P}({\mathcal X}))$ is the set of continuous functions with values in ${\mathcal P}({\mathcal X})$ . If $\zeta\in{\mathcal P}({\mathcal X})$ , f is a real-valued function on ${\mathcal X}$ , and $A{\subset}{\mathcal X}$ , we write

\[\left\langle{\zeta,f}\right\rangle = \int f(z)\, \zeta(\mathop{}\mathopen{}\mathrm{d} z)\text{ \ and\ } \zeta(A) = \int_A \zeta(\mathop{}\mathopen{}\mathrm{d} z).\]

The state space of the process describing the whole network is ${\mathcal S}_N\stackrel{{\text {def}.}} = {\mathcal X}^N$ . For $t \ge 0$ and $i\in\{1,\ldots,N\}$ , $(X_i^N(t))$ (resp. $(Y_i^N(t))$ ) denotes the number of class 1 (resp. class 2) customers at node i.

For $\mu\in\{\mu_1,\mu_2\}$ , $({\mathcal N}_{\mu}^{ij})$ is an independent, identically distributed (i.i.d.) sequence of Poisson processes with rate $\mu$ . They are associated to the service times of class 1 and 2 customers at the nodes of the networks. For $i, j \geq 1$ , $\mathcal{N}_{\mu}^{ij}$ is the Poisson process associated to the service times of jth server of the ith node. Similarly, $({\mathcal N}_{\lambda}^i) = ((t_n^i))$ is an i.i.d. sequence of Poisson processes with rate $\lambda$ . For $i\in{\mathbb N}$ , ${\mathcal N}_{\lambda}^i$ is the arrival process at node i. All Poisson processes are assumed to be independent. Additionally, $(U_{nk}^i)$ is an i.i.d. sequence of uniform random variables on [0, 1], also referred to as ‘marks’ of the point process $\overline{\mathcal N}_{\lambda}^i$ defined by

\begin{align*}\overline{\mathcal N}_{\lambda}^i\stackrel{{\text {def}.}} = \sum_{n\in{\mathbb N}} \delta_{(t_n^i,(U_{nj}^i, j\ge 1))}. \end{align*}

$\overline{\mathcal N}_{\lambda}^i$ is a marked Poisson point process. (See Chapter 5 of Kingman [Reference Kingman21].) If ${\mathcal M}\stackrel{{\text {def}.}} = [0,1]^{\mathbb N}$ is the space of marks, then $ {\mathcal N}_{\lambda}^i(\mathop{}\mathopen{}\mathrm{d} t) = \overline{\mathcal N}_{\lambda}^i(\mathop{}\mathopen{}\mathrm{d} t,{\mathcal M})$ is the arrival process at node i.

A mark is associated to an arrival instant; it is used to determine to which queue a customer goes if the node where it arrives is saturated. For node i, this is (formally) done with a functional $T_i^N$ on ${\mathcal S} \times {\mathcal M}$ defined as follows. For $(z,u) = ((x_i,y_i),(u_j))\in{\mathcal S} \times {\mathcal M}$ , define $(v_j(u)) = (\lfloor N u_j\rfloor)$ , with 0 identified to N. For i, $n\in{\mathbb N}$ , $(v_j(U_{nj}^i))$ is an i.i.d. sequence uniformly distributed on $\{1,\ldots,N\}$ . With this notation, we take

(5) \begin{equation} T_i^N(z,u) = \begin{cases} i &\text{ if } x_i + y_i < C,\\[3pt] v_k(u) &\text{ otherwise, if } k = \inf\{j\mid x_{v_j(u)} + y_{v_j(u)} < C\}\le p_0,\\[3pt] + \infty &\text{ otherwise.} \end{cases}\end{equation}

The variable $T_i^N(z,u)$ will be used in what follows only when node i is saturated.

When the system is in a state $z\in{\mathcal S}_N$ such that $x_i + y_i = C$ for all $i\in\{1,\ldots,N\}$ , if there is an arrival with a mark u at node i, then, provided that $T_i^N(z,u)$ is less than $p_0$ , the queue with index $T_i^N(z,u)$ receives the rerouted customer. When $p_0 = + \infty$ , if $x_j + y_j < C$ for some $j\in\{1,\ldots,N\}$ , then it is easy to check that

\[{\mathbb P}(T_i(z,U) = j)={1}{\bigg/}\sum_{k=1}^N \textbf{1}_{{\mathcal X}_+}(x_k,y_k),\]

where ${\mathcal X}_+$ is defined in (4).

If the addition of marks looks somewhat formal, it has the advantage of giving a neat framework to handle the stochastic calculus associated to the evolution equations of the system. The martingale property will be with respect to the filtration $({\mathcal F}_t)$ defined as follows: for $t \ge 0$ ,

\[{\mathcal F}_t=\sigma\left\{\overline{\mathcal N}_{\lambda}^i([0,s] \times B),{\mathcal N}_{\mu_{1}}^{ij}([0,s]), {\mathcal N}_{\mu_{2}}^{ij}([0,s])\bigg|i, j\in{\mathbb N}, s\le t, B\in{\mathcal B}({\mathcal M})\right\}.\]

Evolution equations. We assume from now on that $p_0 = + \infty$ ; that is, a customer is rejected only if all nodes are saturated. We consider briefly the case of a finite $p_0$ in Proposition 2 and in Subsection 2.3.

The state of the system is represented by a process $(Z^N(t)) = ((X_i^N(t),Y_i^N(t))$ which is càdlàg and satisfies the following stochastic differential equations (SDE): for $1\le i\le N$ ,

(6) \begin{align} \mathop{}\mathopen{}\mathrm{d} X_i^N(t) &= {\textbf{1}_{\left\{{(X_i^N(t-),Y_i^N(t-))\in{\mathcal X}_+}\right\}}}{\mathcal N}_{\lambda}^i(\mathop{}\mathopen{}\mathrm{d} t)-\sum_{\ell=1}^{X_i^N(t-)} {\mathcal N}_{\mu_1}^{i,\ell}(\mathop{}\mathopen{}\mathrm{d} t), \end{align}

(7) \begin{align} \mathop{}\mathopen{}\mathrm{d} Y_i^N(t) &= \sum_{j=1}^N {\textbf{1}_{\left\{{\substack{(X_j^N(t-),Y_j^N(t-))\in{\mathcal X}_+^c\\[3pt] T_j(u,Z^N(t-))=i}}\right\}}}\overline{{\mathcal N}}_{\lambda}^{j}(\mathop{}\mathopen{}\mathrm{d} t, \mathop{}\mathopen{}\mathrm{d} u)\\ &\quad -\sum_{\ell=1}^{Y_i^N(t-)} {\mathcal N}_{\mu_2}^{i,\ell}(\mathop{}\mathopen{}\mathrm{d} t),\notag\end{align}

where $U(t-)$ denotes the left limit of U at $t>0$ .

The first term of the right-hand side of (7) corresponds to arrivals that find their arrival node saturated and are therefore allocated to some non-saturated node by means of repeated random sampling of nodes until one is found that can accommodate it. Note that it is not excluded that all nodes are saturated; in this case the customer is rejected since it is not allocated anywhere.

Martingales. We recall classical results on the martingales associated to marked Poisson point processes. (See Sections 4.4 and 4.5 of Jacobsen [Reference Jacobsen19], for example. See also Last and Brandt [Reference Last and Brandt22].) If h is a bounded function on ${\mathbb R}_+ \times \{1,\ldots,N, + \infty\}$ such that $h(\cdot, + \infty)\,{\equiv}\,0$ , and $\Lambda^N$ is the empirical distribution defined by Equation (1), the process

\begin{align*}&\left(\int_{[0,t]\times{\mathcal M}} h\left(Z^N(s-),T_i^N(Z^N(s-),u)\right)\overline{\mathcal N}_{\lambda}^i(\mathop{}\mathopen{}\mathrm{d} s, \mathop{}\mathopen{}\mathrm{d} u)\right. \\ &\quad \left.-\lambda\sum_{j=1}^N \int_0^t \frac{h(Z^N(s),j)}{N\Lambda^N(s)({\mathcal X}_+)}\textbf{1}_{{\mathcal X}_+}\left(X_j^N(s),Y_j^N(s)\right)\mathop{}\mathopen{}\mathrm{d} s \right)\end{align*}

is a martingale whose previsible increasing process (quadratic variation) is given by

\[\left(\lambda\sum_{j=1}^N \int_0^t \frac{h^2(Z^N(s),j)}{N\Lambda^N(s)({\mathcal X}_+)}\textbf{1}_{{\mathcal X}_+}\left(X_j^N(s),Y_j^N(s)\right)\mathop{}\mathopen{}\mathrm{d} s \right).\]

For example, the integration of SDE (7) and the compensation of the Poisson processes give the relation

\begin{align*} Y_i^N(t) &= Y_i^N(0) +M_i^N(t) \\[3pt] &\quad +\lambda \sum_{j=1}^N \int_0^t {\textbf{1}_{\left\{{\substack{(X_j^N(s),Y_j^N(s))\in{\mathcal X}_+^c\\(X_i^N(s),Y_i^N(s))\in{\mathcal X}_+}}\right\}}}\frac{1}{N\Lambda^N(s)({\mathcal X}_+)}\,\mathop{}\mathopen{}\mathrm{d} s -\mu_2\int_0^t Y_i^N(s) \,\mathop{}\mathopen{}\mathrm{d} s,\end{align*}

which can be written in the more compact form

(8) \begin{align} Y_i^N(t) &= Y_i^N(0) +M_i^N(t)\nonumber\\ &\quad + \lambda\int_0^t {\textbf{1}_{\left\{{(X_i^N(s),Y_i^N(s))\in{\mathcal X}_+}\right\}}} \frac{\Lambda^N(s)({\mathcal X}_+^c)}{\Lambda^N(s)({\mathcal X}_+)}\,\mathop{}\mathopen{}\mathrm{d} s -\mu_2\int_0^t Y_i^N(s) \,\mathop{}\mathopen{}\mathrm{d} s,\end{align}

where $(M_i^N(t))$ is a martingale whose previsible increasing process is

(9) \begin{align}&\left(\left\langle{M_i^N}\right\rangle(t)\right)\nonumber\\ &\quad =\left( \lambda\int_0^t {\textbf{1}_{\left\{{(X_i^N(s),Y_i^N(s))\in{\mathcal X}_+}\right\}}} \frac{\Lambda^N(s)({\mathcal X}_+^c)}{\Lambda^N(s)({\mathcal X}_+)}\,\mathop{}\mathopen{}\mathrm{d} s+\mu_2\int_0^t Y_i^N(s) \,\mathop{}\mathopen{}\mathrm{d} s \right).\end{align}

Empirical distributions. The empirical distribution process $(\Lambda^N(t))$ associated to $(Z^N(t))$ is defined as follows: for $t \ge 0$ ,

(10) \begin{equation}\Lambda^N(t)(f) = \int_{{\mathcal X}} f(z)\Lambda^N(t)(\mathop{}\mathopen{}\mathrm{d} z) \stackrel{{\text {def}.}} = \frac{1}{N}\sum_{i=1}^N f\left(X_i^N(t),Y_i^N(t)\right)\end{equation}

for any nonnegative function f on ${\mathcal X}$ . It is a stochastic process with values in the set ${\mathcal P}({\mathcal X})$ of probability distributions on ${\mathcal X}$ .

As in the derivation of (8), the integration of (6) and (7) and the compensation of Poisson processes give the relation

(11) \begin{align} \left\langle{\Lambda^N(t),f}\right\rangle &= \left\langle{\Lambda^N(0),f}\right\rangle+M_f^N(t)\nonumber\\[3pt] &\quad +\lambda\int_0^t\int_{{\mathcal X}_+} \nabla_1^+(f)(z)\,\Lambda^N(s)(\mathop{}\mathopen{}\mathrm{d} z)\,\mathop{}\mathopen{}\mathrm{d} s-\mu_1\int_0^t\int_{z = (x,y)\in{\mathcal X}} x\nabla_1^-(f)(z)\,\Lambda^N(s)(\mathop{}\mathopen{}\mathrm{d} z)\,\mathop{}\mathopen{}\mathrm{d} s\nonumber\\[3pt] &\quad + \lambda\int_0^t \int_{{\mathcal X}_+} \nabla_2^+(f)(z)\frac{\Lambda^N(s)({\mathcal X}_+^c)}{\Lambda^N(s)({\mathcal X}_+)}\,\Lambda^N(s)(\mathop{}\mathopen{}\mathrm{d} z)\,\mathop{}\mathopen{}\mathrm{d} s\nonumber\\[3pt] &\quad -\mu_2\int_0^t \int_{z = (x,y)\in{\mathcal X}} y\nabla_2^-(f)(z)\,\Lambda^N(s)(\mathop{}\mathopen{}\mathrm{d} z)\,\mathop{}\mathopen{}\mathrm{d} s\end{align}

for a real-valued function f on ${\mathcal X}_+$ , where, if $z = (x,y)$ , then $\nabla_1^\pm(f)(z) = f(x{\pm}1,y)-f(x,y)$ and $\nabla_2^\pm(f)(z) = f(x,y{\pm}1)-f(x,y)$ are the gradient operators, and $(M_f^N(t))$ is a martingale.

Proposition 1. (Dynamical system with an unbounded number of retrials) If $p_0 = + \infty$ and $(\Lambda(t))$ is the unique solution of the differential equation

(12) \begin{align} \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t}\left\langle{\Lambda(t),f}\right\rangle &= \lambda\left\langle{\Lambda(t), \nabla_1^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle+\mu_1\left\langle{\Lambda(t),I_1\nabla_1^-(f)}\right\rangle \nonumber\\[3pt] &\quad + \lambda \frac{\Lambda(t)({\mathcal X}_+^c)}{\Lambda(t)({\mathcal X}_+)} \left\langle{\Lambda(t),\nabla_2^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle +\mu_2\left\langle{\Lambda(t),I_2\nabla_2^-(f)}\right\rangle \end{align}

for $t < H_0(\Lambda)$ , where $I_1(x,y)\stackrel{{{\text {def}}.}} = x$ and $I_2(x,y)\stackrel{{{\text {def}}.}} = y$ , and

\[H_0(\zeta)\stackrel{{ {\text {def}}.}} = \inf\{t>0\mid \zeta(t)({\mathcal X}_+) = {0}\} \qquad\text{ for } (\zeta(t))\in{\mathcal C}({\mathbb R}_+,{\mathcal P}({\mathcal X})),\]

then, for the convergence in distribution of processes,

\[\lim_{N\to+\infty} \left(\Lambda^N(t),t < H_0\left(\Lambda^N\right)\right)= (\Lambda(t),t < H_0(\Lambda)).\]

Note that the sequence of processes $(\Lambda^N(t))$ is in fact a sequence of finite-dimensional processes with dimension $\mathop{}\mathopen{}\mathrm{card} ({\mathcal X}_+)$ . The convergence in distribution of the proposition refers to the case when the space of càdlàg functions with values in ${\mathbb R}_+^{\mathop{}\mathopen{}\mathrm{card}({\mathcal X}_+)}$ is endowed with the uniform norm.

The variable $H_0(\Lambda)$ is the blow-up time of the dynamical system (12). If finite, it amounts to the fact that the system is completely saturated at the fluid scale. It will be seen in Subsection (2.4) that the saturated state is a stable equilibrium of the network. Note that, because of its singular aspect, it cannot be really defined through the ODEs associated to (12).

Note that, for $s < t$ ,

\[ \Bigg|\int_s^t \int_{{\mathcal X}_+} \nabla_2^+(f)(z)\frac{\Lambda^N(u)({\mathcal X}_+^c)}{\Lambda^N(u)({\mathcal X}_+)}\,\Lambda^N(u)(\mathop{}\mathopen{}\mathrm{d} z)\Bigg|\leq 2\|f\|_{\infty}(t-s). \]

By using the criterion of the modulus of continuity and Equation (11), we get that the sequence of processes $(\Lambda^N(t))$ is tight in distribution for the topology of uniform convergence on compact sets. (See Theorem 7.3 of Billingsley [Reference Billingsley3], for example.)

The solution of (12) lives in a finite-dimensional state space, the set of probability distributions on ${\mathcal X}$ . The system (12) can be seen as a set of ODEs. It has in particular a unique solution up to time $H_0(\Lambda)$ . If

\[H_{\varepsilon}(\zeta)=\inf\{t>0\mid \zeta(t)({\mathcal X}_+)>{\varepsilon}\}\]

then $H_{\varepsilon}(\Lambda^N)$ converges in distribution to $H_{\varepsilon}(\Lambda)$ . By the continuous mapping theorem used in (11), one sees that on the event $\{H_{\varepsilon}>t\}$ , the relation

\begin{align*}\left\langle{\Lambda(t),f}\right\rangle = &\left\langle{\Lambda(0),f}\right\rangle + \lambda\int_0^t \left\langle{\Lambda(s), \nabla_1^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s + \mu_1\int_0^t\left\langle{\Lambda(s),I_1\nabla_1^-(f)}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s\\[5pt] & + \lambda \int_0^t\frac{\Lambda(s)({\mathcal X}_+^c)}{\Lambda(s)({\mathcal X}_+)} \left\langle{\Lambda(s),\nabla_2^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s +\mu_2\int_0^t\left\langle{\Lambda(s),I_2\nabla_2^-(f)}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s\end{align*}

holds. These equations can be seen as a system of ODEs which clearly has a unique solution up to blow-up time $H_0$ . The proposition is proved.

We state the analogous result when the number of retrials is finite; the notation is the same as in Proposition 1. The proof of the proposition being simpler in this case, it is skipped.

Proposition 2. (Dynamical system with a maximum of $p_0$ retrials) If $(\Lambda(t))$ is the unique solution of the differential equations

(13) \begin{align} \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t}\left\langle{\Lambda(t),f}\right\rangle = & \lambda\left\langle{\Lambda(t), \nabla_1^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle+\mu_1\left\langle{\Lambda(t),I_1\nabla_1^-(f)}\right\rangle \nonumber\\ & + \lambda \Lambda(t)({\mathcal X}_+^c)\left(\frac{1-[\Lambda(t)({\mathcal X}_+^c)]^{p_0}}{1-\Lambda(t)({\mathcal X}_+^c)}\right)\left\langle{\Lambda(t),\nabla_2^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle +\mu_2\left\langle{\Lambda(t),I_2\nabla_2^-(f)}\right\rangle\end{align}

for $t>0$ , then, for the convergence in distribution of processes, the relation

\[\lim_{N\to+\infty} \left(\Lambda^N(t)\right)= (\Lambda(t))\]

holds.

The nonlinear term in the right-hand side of (13) is bounded by $p_0$ , and is therefore without singularity.

2.2. The number of equilibrium points

We investigate the fixed points of the linearized version $(\Lambda_R(t))$ of the dynamical system (12) when the nonlinear term $\Lambda(t)({\mathcal X}_+^c)$ is replaced by a constant $R\in(0,1)$ . This dynamical system satisfies the following equations:

(14) \begin{align}\left\langle{\Lambda_R(t),f}\right\rangle = & \left\langle{\Lambda_R(0),f}\right\rangle + \lambda\int_0^t \left\langle{\Lambda_R(s), \nabla_1^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s\nonumber\\ & + \mu_1\int_0^t\left\langle{\Lambda_R(s),I_1\nabla_1^-(f)}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s + \lambda \frac{R}{1-R}\int_0^t \left\langle{\Lambda_R(s),\nabla_2^+(f)\textbf{1}_{{\mathcal X}_+}}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s\nonumber\\ & +\mu_2\int_0^t\left\langle{\Lambda_R(s),I_2\nabla_2^-(f)}\right\rangle\,\mathop{}\mathopen{}\mathrm{d} s\end{align}

for any real-valued function f on ${\mathcal X}$ , with the notation of Proposition 1. Recall that for any function f and measure $\mu$ on ${\mathcal X}$ ,

\[\left\langle{\mu,f}\right\rangle = \sum_{(x,y){\mathcal X}} f(x,y)\mu((x,y)).\]

The process $(\Lambda_R(t))$ describes the evolution of the law of a classical Erlang model with capacity C where two classes of customers arrive at rates $\lambda$ and $\lambda R/(1-R)$ respectively and are served at rates $\mu_1$ and $\mu_2$ . Its invariant distribution on ${\mathcal X}$ is given by

(15) \begin{equation}\pi_R(x,y)=\frac{1}{Z_R} \cdot \frac{\rho_1^x}{x!}\frac{\rho_2^y}{y!}\left(\frac{R}{1-R}\right)^y,\quad (x,y)\in{\mathcal X},\end{equation}

where $\rho_1\,{=}\,\lambda/\mu_1$ and $\rho_2\,{=}\,\lambda/\mu_2$ , and $Z_R$ is the normalization constant

\[Z_R=\sum_{m=0}^C \sum_{x+y=m} \frac{\rho_1^x}{x!}\frac{\rho_2^y}{y!}\left(\frac{R}{1-R}\right)^y=\sum_{m=0}^C \frac{1}{m!} \left(\rho_1 + \rho_2\frac{R}{1-R}\right)^m.\]

We have

\[\pi_R\left({\mathcal X}_+^c\right)=\pi_R((x,y)\in{\mathcal X}, x + y=C)=\frac{1}{Z_R} \frac{1}{C!} \left(\rho_1 + \rho_2\frac{R}{1-R}\right)^C.\]

Hence $\pi_R$ as defined by (15) is a fixed point of the dynamical system (12) if and only if R satisfies the relation

(16) \begin{equation}R=\pi_R\left({\mathcal X}_+^c\right).\end{equation}

The goal of this section is to characterize completely the solutions of (16).

Note that (16) can be rewritten as $\Phi_{\rho_1,\rho_2}(z_R)\,{=}\,0$ with, for $\rho_1$ , $\rho_2\,{>}\,0$ ,

(17) \begin{equation}\Phi_{\rho_1,\rho_2}(z)\stackrel{{\text {def}.}} = \frac{1}{z}\frac{(\rho_1 + \rho_2 z)^C}{C!}-\sum_{m=0}^{C-1} \frac{ (\rho_1 + \rho_2 z)^m}{m!},\end{equation}

and $z_R \,{=}\,{R}/{(1-R)}$ . The next proposition determines the number of roots of the function $\Phi_{\rho_1,\rho_2}$ . It gives all the nonsingular equilibrium points of the ODEs (12). Subsection 2.4 gives a formal presentation of another regime which yields a stable equilibrium when $\rho_2>C$ ; see the remark below.

We will refer to the different cases of the proposition as, respectively, underloaded if $\rho_2 < C$ , critical if $\rho_2 = C$ , and overloaded if $\rho_2>C$ .

Proof. Once the results to be proved have been properly formulated, the proofs of the statements will be done via real analysis. First, note that

(19) \begin{equation} \lim_{z\searrow 0} z\Phi_{\rho_1,\rho_2}(z)=\frac{\rho_1^C}{C!}>0 \text{ and } \lim_{z\to+\infty} \frac{\Phi_{\rho_1,\rho_2}(z)}{z^{C-1}}=\frac{\rho_2^{C-1}}{(C-1)!}\left(\frac{\rho_2}{C}-1\right).\end{equation}

The underloaded case ${\boldsymbol{\rho}}_{\boldsymbol{2}}\in({\boldsymbol{0}},C)$

Define, for $z \ge 0$ ,

\[f_{\rho_1,\rho_2}(z) = \Phi_{\rho_1,\rho_2}(z)e^{-\rho_2 z}.\]

After some simple calculations with telescoping sums, we get that

(20) \begin{equation}f^{\prime}_{\rho_1,\rho_2}(z)= \frac{1}{z^2}\frac{(\rho_1 + \rho_2 z)^{C-1}}{C!}e^{-\rho_2 z} \left(\rho_2(C-\rho_2) z^2+\rho_2(C-\rho_1-1)z-\rho_1\right).\end{equation}

The last term of the right-hand side of the expression for $f^{\prime}_{\rho_1,\rho_2}(z)$ is a polynomial of degree 2. Its value at $z = 0$ is negative and, since $\rho_2 < C$ , it converges to $ + \infty$ as z gets large; hence there exists a unique $z_0>0$ such that $f^{\prime}_{\rho_1,\rho_2}(z_0) = 0$ . It necessarily corresponds to a unique extremum of $f_{\rho_1,\rho_2}$ on ${\mathbb R}_+$ , which is a minimum, given the variations of $f_{\rho_1,\rho_2}$ . (See (19).) The function $f_{\rho_1,\rho_2}$ is decreasing on $(0,z_0)$ and increasing on $[z_0, + \infty)$ , and since $f_{\rho_1,\rho_2}$ converges to 0 at infinity, the relation $f_{\rho_1,\rho_2}(z) < 0$ holds for $z \ge z_0$ . Thus, there exists a unique root $z(\rho_1,\rho_2)$ of $f_{\rho_1,\rho_2}$ in ${\mathbb R}_+$ , located in $(0,z_0)$ . We have therefore the equivalence of the two relations $z < z(\rho_1,\rho_2)$ and $\Phi_{\rho_1,\rho_2}(z)>0$ . The latter relation can be expressed as

\[z\sum_{m=0}^{C-1}\frac{C!}{m!} \frac{1}{(\rho_1 + \rho_2 z)^{C-m}}<1. \]

Hence, if $\rho^{\prime}_{2}>\rho_2$ , then $\Phi_{\rho_1,\rho_2^{\prime}}(z)>0$ holds if $z < z(\rho_1,\rho_2)$ . The function $z(\rho_1,\cdot)\,:\,{\mathbb R}_+{\to}{\mathbb R}_+$ , $\rho_2{\mapsto}z(\rho_1,\rho_2)$ , is therefore increasing.

The critical case ${\boldsymbol{\rho}}_{\boldsymbol{2}} = C$

We easily get the relation

(21) \begin{equation} \lim_{z\to+\infty} \frac{z^2\Phi_{\rho_1,C}(z)}{(\rho_1 + C z)^C}C!=\frac{\rho_1-C + 1}{C}.\end{equation}

If $\rho_1 < C-1$ , then $\Phi_{\rho_1,C}$ has negative values and, by (19), a root $z_C>0$ . The root is unique; otherwise $f_{\rho_1,C}$ would have two distinct extrema, which is impossible by (20). The limit $z(\rho_1,C)$ is necessarily a root of $\Phi_{\rho_1,C}$ ; hence $z(\rho_1,C) = z_C$ .

If $\rho_1 \ge C-1$ , then (20) gives that $f_{\rho_1,C}$ is strictly decreasing. Since it converges to 0, we conclude that $\Phi_{\rho_1,C}$ does not have a root in this case, and thus, necessarily, $z(\rho_1,C) = + \infty$ .

The overloaded case $\boldsymbol{\rho_2}>C$

By (19), the function $\Phi_{\rho_1,\rho_2}$ converges to $ + \infty$ at 0 and at $ + \infty$ . To determine the number of roots of $\Phi_{\rho_1,\rho_2}$ , one has therefore to obtain the sign of $m(\rho_1,\rho_2)$ as defined in (18).

Clearly, $m(\cdot,\rho_2)\,:\,{\mathbb R}_+{\to}{\mathbb R}_+$ , $\rho_1{\mapsto}m(\rho_1,\rho_2)$ , is a continuous function. For $\rho_1^0 < \rho_1^1$ and $z>\rho_1^1/\rho_2$ , the relation

\[\Phi_{\rho_1^0,\rho_2}(z-\rho_1^0/\rho_2) <\Phi_{\rho_1^1,\rho_2}(z-\rho_1^1/\rho_2) \]

implies that $m_{\rho_1^0}\le m_{\rho_1^1}$ ; thus $m(\cdot,\rho_2)$ is an increasing function. Note that, for $z>0$ ,

\[ \lim_{\rho_1\to 0} \Phi_{\rho_1,\rho_2}(\rho_1 z)=-1; \]

hence if $\rho_1$ is sufficiently small then $m(\rho_1,\rho_2) < 0$ . If $\rho_1>C-1$ , then from (20) we get that the function $f_{\rho_1,\rho_2}$ is decreasing and converges to 0 at infinity. Consequently $\Phi_{\rho_1,\rho_2}$ is positive on ${\mathbb R}_+$ , and therefore $m(\rho_1,\rho_2)>0$ . Hence $ \phi_C$ (as defined in (18)) satisfies $\phi_C(\rho_2)\in(0,C-1)$ . By continuity of $m(\cdot,\rho_2)$ , we have $m_{\phi_C(\rho_2)} = 0$ ; the function $\Phi_{\phi_C(\rho_2),\rho_2}$ has therefore a unique root.

By definition, $\rho_1>\phi_C(\rho_2)$ if and only if $m(\rho_1,\rho_2) < 0$ , or, equivalently, since the minimum of $\Phi_{\rho_1,\rho_2}$ is necessarily achieved in some compact interval of $(0, + \infty)$ , if the relation

\[ \sup_{z>0}\left(z\sum_{m=0}^{C-1}\frac{C!}{m!} \frac{1}{(\rho_1 + \rho_2 z)^{C-m}}\right)<1 \]

holds. We deduce that $\phi_C$ is a non-increasing function on $(C, + \infty)$ and that $\phi_C(\rho_2)$ converges to 0 as $\rho_2$ gets large. If

\[ \delta_C\stackrel{{\text {def}.}} = \lim_{{\rho_2\searrow C}}\uparrow\phi_C(\rho_2) < C-1, \]

one can take $\rho_1\in(\delta_C,C-1 )$ . From the critical case, we know that $\Phi_{\rho_1,C}$ has exactly one root and that there exists $z_1>0$ such that $\Phi_{\rho_1,C}(z_1) < 0$ . One can fix $\rho_2^1>C$ sufficiently close to C so that $\Phi_{\rho_1,\rho_2}(z_1) < 0$ . Hence $\Phi_{\rho_1,\rho_2^1}$ has a root, and consequently $\rho_1\le\phi_C(\rho_2^1)$ , which contradicts the fact that $\rho_1>\delta_C \ge \phi_C(\rho_2^1)$ . The proposition is proved.

Examples.

1. 1. When $C = 2$ and $\rho_1 < 2 < \rho_2$ , one has

   \[ 2\Phi_{\rho_1,\rho_2}(z)= \rho_2\left( \rho_2-2 \right) {z}^{2} + 2\left( \rho_1\rho_2-\rho_1-1\right) z + {\rho_1}^{2}. \]
   The function $\Phi_{\rho_1,\rho_2}$ has two roots on ${\mathbb R}_+$ as long as $\Phi^{\prime}_{\rho_1,\rho_2}(0) < 0$ —that is, $\rho_1 < 1/(\rho_2-1)$ —and the minimum of $\Phi_{\rho_1,\rho_2}$ is negative—that is,
   \[ \rho_1^2+2(1-\rho_2)\rho_1+1>0. \]
   It is then easy to deduce that
   \[ \phi_2(\rho_2)=\rho_2-1-\sqrt{\rho_2(\rho_2-2)}. \]

2. 2. When $C = 3$ and $\rho_1 < 3 < \rho_2$ , one has

   \begin{align*} &6\Phi_{\rho_1,\rho_2}(z) = {\rho_2}^{2} \left(\rho_2-3 \right) {z}^{3} + 3\rho_2 \left( \rho_1\rho_2-2\rho_1-2 \right) {z}^{2}\\ &\quad + 3 \left( {\rho_1}^{2}\rho_2-{\rho_1}^{2}-2\rho_1-2 \right) z + {\rho_1}^{3}.\end{align*}
   The discriminant of this polynomial (in z) is
   \[H_{\rho_2}(u)\stackrel{{\text {def}.}} = 3{\rho_1}^{4} + 2\left(6-5\rho_2\right) {u}^{3} + 3\left(3{\rho_2}^{2}-8(\rho_2-1)\right){u}^{2}-12\left(\rho_2-2\right) u-8\rho_2 + 12.\]
   Since $H_{\rho_2}(0) = -12(\rho_2-2) < 0$ and $H_{\rho_2}(C-1) = 4(9\rho_2-25))(\rho_2-3)>0$ for $\rho_2>3$ , it is then easily seen that
   \[\phi_3(\rho_2)=\inf\{u\,:\,H_{\rho_2}(u) = 0\}.\]

2.3. The RIST algorithm with one retrial

Proposition 3 shows that the dynamical system associated to the RIST algorithm with an infinite number of retrials has at most two equilibrium points. As noted earlier, there is another equilibrium which is not described here; it is investigated in Subsection 2.4.

The purpose of this section is to show that in the case when only one attempt to accommodate a request is allowed, there are also cases with three equilibrium points, but all of them are equilibria of the ‘smooth’ dynamical system.

In view of Proposition 2, and in the same way as in the derivation of (17), an equilibrium point is of the form $\pi_S$ , with

(22) \begin{equation}\pi_S(x,y)=\frac{1}{Z_S} \cdot \frac{\rho_1^x}{x!}\frac{(\rho_2 S)^y}{y!},\quad (x,y)\in{\mathcal X},\end{equation}

where S is a solution $z\in(0,1)$ of

(23) \begin{equation}\frac{(\rho_1 + \rho_2 z)^C}{C!}-z\sum_{m=0}^{C} \frac{ (\rho_1 + \rho_2 z)^m}{m!}=0,\end{equation}

and $Z_S$ is the normalization constant.

With the change of coordinates $z\ {\mapsto}\ (z-\rho_1)/\rho_2$ , this amounts to finding the roots $z\in(\rho_1,\rho_1 + \rho_2)$ of g, where

\[g(z)\stackrel{{\text {def}.}} = \frac{\rho_2}{z-\rho_1}\frac{z^C}{C!}e^{-z}-\sum_{m=0}^{C} \frac{ z^m}{m!}e^{-z}.\]

Note that $g(s){\to} + \infty$ when $s{\searrow}\rho_1$ , and $g(\rho_1 + \rho_2) < 0$ . Simple calculations give the relation

\[g^{\prime}(z)=\frac{z^{C-1}}{C!}e^{-z}f(z),\]

with $f(z) = z^3-(2\rho_1 + \rho_2) z^2 + (\rho_1^2 + \rho_1\rho_2 + \rho_2(C-1)) z-C\rho_1\rho_2$ . This shows in particular that Equation (23) cannot have more than three solutions.

We give a scaling picture of the fixed point equation (23). As with the DAR algorithm, which is investigated in the next section, we study the case when the capacity C is a scaling parameter going to infinity and $\rho_2$ is of the order of C, i.e. $\rho_2 = \nu_2 C$ for some $\nu_2>0$ .

Under some conditions, there is always a solution of Equation (23) close to 0. For that we do the change of variable $z{\mapsto}\rho_1z$ ; the relation becomes

\[\psi_{1,C}(z)\stackrel{{\text {def}.}} = -z+\rho_1^{C-1}\frac{({1+}\nu_2C z)^C}{C!}-z\sum_{m=1}^{C} \rho_1^m\frac{ (1 + \nu_2C z)^m}{m!}=0.\]

Using Stirling’s formula, it is easily seen that, for ${\varepsilon}>0$ ,

\[\limsup_{C\to + \infty} \psi_{1,C}({\varepsilon})\leq -{\varepsilon}+\lim_{C\to + \infty}\rho_1^{C-1}\frac{({1+}\nu_2C {\varepsilon})^C}{C!}=-{\varepsilon},\]

provided that $\nu_2\rho_1\geq e$ . Since $\psi_{1,C}(0)>0$ , we get that, for C sufficiently large, Equation (23) has a solution in the interval $(0,(\rho_1{\varepsilon}){\wedge}1)$ .

Returning to Equation (23), it can be written as

\[\psi_{2,C}(z)\stackrel{{\text {def}.}} = z\sum_{k=0}^{C}\frac{1}{(\rho_1 + \nu_2C z)^k} \frac{C!}{(C-k)!}-1=0.\]

Assuming that $\nu_2 z>1$ , we can check that

\[\lim_{C\to + \infty} \psi_{2,C}(z)=\Psi_{2,\infty}(z)\stackrel{{\text {def}.}} = \frac{\nu_2z^2-\nu_2z + 1}{\nu_2z-1},\]

and the convergence is uniform on any compact interval of $(1/\nu_2,1]$ . If $\nu_2>4$ , the function $\Psi_{2,\infty}$ has two zeroes in the interval $(1/\nu_2,1)$ , given by

\[z^*_{\pm}\stackrel{{\text {def}.}} = \frac{1\pm\sqrt{1-4/\nu_2}}{2}.\]

We now summarize this result in the following proposition.

Proposition 4. Under the assumption that $\rho_2 = \nu_2C$ and $\nu_2>\max(4,e/\rho_1)$ , there exists $C_0>0$ such that if $C \ge C_0$ , the dynamical system (13) associated to the RIST algorithm with one retrial has exactly three equilibrium points, which converge respectively to

\[0,\; \frac{1-\sqrt{1-4/\nu_2}}{2}, \; \frac{1 + \sqrt{1-4/\nu_2}}{2}\]

as C goes to infinity.

The proof is skipped since most of the arguments have been given and the proof of Theorem 2 in the next section is similar and slightly more technical.

2.4. Stability of saturation

In this section it is assumed that $\rho_2>C$ . The main result is that if the initial state is sufficiently congested, so is the state of the network on any finite time interval.

We describe the general strategy of our approach. Note first that the two-dimensional process of the total number of empty places and total number of class 1 customers in the network does not have the Markov property. Indeed, in state (m, n) with $m>0$ , when an external job arrives it is not possible to determine if the transition to $(m-1,n)$ or to $(m-1,n + 1)$ occurs, i.e. if the job is blocked at its arrival queue or not. This process can in fact be compared with an ergodic Markov process in ${\mathbb N}^2$ by using a convenient coupling and a specific order relation in ${\mathbb N}^2$ . The ergodicity property is then used to show that, asymptotically, the total number of empty places and of class 1 customers in the network is negligible with respect to N, so that the entire system is composed of rerouted jobs in the limit.

If $(Z^N(t)) = (X_i^N(t),Y_i^N(t))$ is the solution of the SDE (6) and (7), define

\[\overline{Z}_1^N(t)\stackrel{{\text {def}.}} = \sum_{i=1}^N X_i^N(t),\quad \overline{Z}_2^N(t)\stackrel{{\text {def}.}} = \sum_{i=1}^N Y_i^N(t),\]

and let $\overline{Z}_0^N(t) = CN-\overline{Z}_1^N(t)-\overline{Z}_2^N(t)$ denote the total number of empty places at time t.

A coupling with a two-dimensional Markov process. We introduce a Markov process in ${\mathbb N}^2$ which will be used in the analysis of the asymptotic behavior of the process of the empirical distribution associated to $(Z^N(t))$ .

Definition 1. Let $(U^N(t)\stackrel{{\text {def}.}} = (U_0^N(t),U_1^N(t))$ be the Markov process on ${\mathbb N}^2$ with the initial state $(\overline{Z}_0^N(0),\overline{Z}_1^N(0))$ , and Q-matrix $Q = (q({\cdot},{\cdot}))$ defined on the state space ${\mathcal S}_U = \{u = (u_0,u_1)\in{\mathbb N}^2\,:\,u_0 + u_1{\leq} N\}$ as follows:

\[\begin{cases}q(u,u-e_0 + e_1) = \lambda u_0, \\q(u,u-e_0) = \lambda(N-u_0),\\q(u,u + e_0-e_1) = \mu_1u_1,\\q(u,u + e_0) = \mu_2(CN-u_0-u_1),\end{cases}\]

with $e_0 = (1,0)$ and $e_1 = (0,1)$ , provided that the transitions keep the process in ${\mathcal S}_U$ . Define

\[U_2^N(t)\stackrel{{\text {def}.}} = CN-U_0^N(t)-U_1^N(t).\]

With some abuse of notation, we will also speak of the Z- and U-systems to refer to the associated stochastic processes $(Z^N(t))$ and $(U^N(t))$ , and similarly, we will speak of class 1 and 2 jobs in the U-system with the obvious meaning. Finally,

(24) \begin{equation}S_U^N\stackrel{{\text {def}.}} = \inf\left\{t\geq 0\,:\, U_0^N(t) + U_1^N(t) = N\right\},\end{equation}

and for $a\in{\mathbb N}$ ,

(25) \begin{equation}T_{U}^a\stackrel{{\text {def}.}} = \inf\left\{t\geq 0\,:\, U_0^N(t)\geq a\right\}.\end{equation}

The transition rates of the process $(U^N(t))$ suggest that this process would behave like $(Z^N(t))$ if all nodes had, at most, one empty place. The coupling shows that the process $(Z^N(t))$ can be bounded above in some way by such a process. A key element in the coupling is the explicit use of the fact that the service times of class 1 jobs are ‘smaller’ than the service times of class 2 jobs.

Proposition 5. There exists a coupling of the processes $(U_0^N(t),U_1^N(t))$ and $(Z^N(t))$ such that $(U_0^N(0),U_1^N(0)) = (\overline{Z}_0^N(0),\overline{Z}_1^N(0))$ and the relations

(26) \begin{equation} \begin{cases}U_2^N(t)\leq \overline{Z}_2^N(t),\\[5pt] U_1^N(t) + U_2^N(t)\leq \overline{Z}_1^N(t) + \overline{Z}_2^N(t) \end{cases} \end{equation}

hold for all $t < S_U^N$ , where $S_U^N$ is defined by (24).

Proof. We proceed by induction on the number of jumps. One has to show that if the relation holds initially then it will also hold at the first jump of $(Z^N(t))$ or $(U_0^N(t),U_1^N(t))$ .

For $j\in\{0,1\}$ , define $z_j\stackrel{{\text {def}.}} = \overline{Z}_j^N(0)$ and $u_j\stackrel{{\text {def}.}} = U_j^N(0)$ . By assumption,

\[\begin{cases}u_2\le z_2,\\u_1 + u_2\le z_1 + z_2,\end{cases}\]

and we can assume that $u_0 < N$ since the process is stopped at time $S_U^N$ . We define

\[\widetilde{a}_0=\sum_{i=1}^N {\textbf{1}_{\left\{{X_{i}^N(0)+Y_i^N(0)<C}\right\}}},\]

the number of non-saturated queues. Clearly,

\[\widetilde{a}_0\le z_0 = CN-z_1-z_2\le CN-u_1-u_2=u_0.\]

We will take the convention that $E_{\xi}$ denotes an exponentially distributed random variable with parameter $\xi \ge 0$ and that all exponential random variables constructed are independent. The coupling is done by introducing the following random variables; the minimum of them will define the first jump of the process. For each random variable, the transition is indicated for $(\overline{Z}^N(t)) = (\overline{Z}_0^N(t),\overline{Z}_1^N(t))$ and for $(U^N(t))$ in the case that it has the minimal value.

1. 1. Arrivals

   1. (a) $E_{\lambda \widetilde{a}_0}$ is the minimum of the arrivals of jobs finding a non-congested queue in the Z-system. Transition: $z{\mapsto}z + e_1-e_0$ and $u{\mapsto} u + e_1-e_0$ .

   2. (b) $E_{\lambda ({u}_0-\widetilde{a}_0))}$ is the minimum of the remaining arrivals of jobs finding a non-congested queue in the U-system. Transition: $u{\mapsto} u + e_1-e_0$ and, if $\widetilde{a}_0>0$ , $z{\mapsto}z-e_0$ .

   3. (c) If $u_0>0$ , $E_{\lambda (N-u_0)}$ is the minimum of the arrivals of jobs finding a congested queue in the Z-system and the U-system. Transition: $u{\mapsto}u-e_0$ and, if $\widetilde{a}_0>0$ , $z{\mapsto} z-e_0$ .

2. 2. Services

   1. (a) $E_{\mu_2 u_2}$ is the minimum of the services of $u_2$ class 2 jobs of the U-system and Z-system. Transition: $z{\mapsto}z + e_0$ and $u{\mapsto} u + e_0$ . Recall that $u_2\le z_2$ .

   2. (b) If $u_1\le z_1$ :

      1. i. $E_{\mu_1 u_1}$ is the minimum of the services of $u_1$ class 1 jobs of the U-system and Z-system. Transition: $z{\mapsto}z + e_0-e_1$ and $u{\mapsto} u + e_0-e_1$ .

      2. ii. $E_{\mu_1 (z_1-u_1)}$ is the minimum of the remaining services of the Z-system. In this case, note that one has necessarily $z_1>u_1$ . Transition: $z{\mapsto}z + e_0-e_1$ and $u{\mapsto} u$ .

      3. iii. $E_{\mu_2 (z_2-u_2)}$ is the minimum of the services of the remaining $z_2-u_2$ class 2 jobs of Z-system. In this case $z_2>u_2$ . Transition: $z{\mapsto}z + e_0$ and $u{\mapsto} u$ .

   3. (c) If $u_1>z_1$ , and we fix $F = E_{u_1-z_1}$ :

      1. i. $E_{\mu_1z_1}$ is the minimum of the services of $z_1$ class 1 jobs of the U-system and Z-system. Transition: $z{\mapsto}z + e_0-e_1$ and $u{\mapsto} u + e_0-e_1$ .

      2. ii. $F/\mu_1$ is the minimum of the services of $u_1-z_1$ remaining class 1 jobs of the U-system. Transition: $z{\mapsto}z$ and $u{\mapsto} u + e_0-e_1$ .

      3. iii. $F/\mu_2$ is the minimum of the services of some $u_1-z_1$ class 2 jobs of the Z-system. Transition: This cannot be the next step, since $F/\mu_1 < F/\mu_2$ ; the transition of (2) (c) (ii) therefore occurs earlier.

      4. iv. $E_{\mu_2(z_2-u_2-(u_1-z_1)}$ is the minimum of the services of the remaining $\mu_2(z_2-u_2-(u_1-z_1)$ class 2 jobs of the Z-system. In this case, $u_1 + u_2 < z_1 + z_2$ . Transition: $z{\mapsto}z + e_0$ and $u{\mapsto} u$ .

An easy, but somewhat tedious, check shows that the two processes have the correct time evolution, and furthermore that the order relation is preserved after any of the transitions mentioned above.

An asymptotic analysis of $(U^N(t))$ . It is assumed that the initial state of the process $(U^N(t))$ of Definition (1) satisfies the relation

(27) \begin{equation}\lim_{N\to+\infty} \frac{1}{N}\left(U_0^N(0),U_1^N(0)\right)=(a_0,a_1).\end{equation}

Proposition 6. Under the condition $\rho_1 < C < \rho_2$ , there exists $\eta_0>0$ such that, if the initial condition (27) satisfies the relations $0\le a_0 + a_1\le\eta_0$ , then for any ${\varepsilon}>0$ , there is $t_0>0$ and a constant $K_0$ such that, for any $T>0$ ,

\[\lim_{N\to+\infty}{\mathbb P}\left(\sup_{t_0\le t\le t_0 + T}U_0^N(t){\leq} K_0\log N, \sup_{t_0\le t\le t_0 + T}\frac{U_1^N(t)}{N}\leq {\varepsilon},\quad S_U^N\geq t_0 + T\right)=1.\]

Proof. For $0 < a_0 + a_1 < \eta < 1$ , by Definition 1 of the process $(U^N(t))$ , a simple coupling shows that the process $(U_1^N(t{\wedge}S_U^N{\wedge}T_U^{\eta N}))$ can be stochastically bounded above by $(L_{\eta}^N(t{\wedge}S_U^N{\wedge}T_U^{\eta N}))$ , where $(L_{\eta}^N(t))$ is the process of the number of jobs of an $M/M/\infty$ queue with arrival rate $\lambda\eta N$ , service rate $\mu_1$ , and initial point $U_1^N(0)$ , and $T_U^{\eta N}$ is defined by (25). A classical result (see for example Theorem 6.13 of Robert [Reference Robert29]) gives the following convergence in distribution:

(28) \begin{equation}\lim_{N\to+\infty} \left(\frac{L_{\eta}^N(t)}{N}\right)=\left(\frac{\lambda}{\mu_1}\eta + \left(a_1-\frac{\lambda}{\mu_1}\eta\right)e^{-\mu_1t}\right).\end{equation}

By a similar argument, the process $(U_0^N(t{\wedge}S_U^N{\wedge}T_U^{\eta N}))$ can be bounded above by $(Q((Nt){\wedge}S_U^N{\wedge}T_U^{\eta N}))$ , where (Q(t)) is the process of the number of jobs of an $M/M/1$ queue with respective arrival and service rates $\mu_2C + \eta(\mu_1-\mu_2)$ and $\lambda$ , and initial point $U_0^N(0)$ . Again a classical result (see for example Proposition 5.16 of Robert [Reference Robert29]) gives the following convergence in distribution:

(29) \begin{equation}\lim_{N\to+\infty}\left(\frac{Q(Nt)}{N}\right)=\left(a_0+\left((\mu_1-\mu_2)\eta + \mu_2C-\lambda\right)t\right)^+,\end{equation}

where $a^+ = \max(a,0)$ for $a\in{\mathbb R}$ .

Now, we fix $0 < \eta_0 < 1$ such that

(30) \begin{equation}(\mu_1-\mu_2)\eta_0 + \mu_2C < \lambda \text{ {\rm and } } \eta_0\frac{\lambda}{\mu_1} < 1.\end{equation}

If $(a_0,a_1)$ is such that

\[a_0 < \eta_0{\wedge}(1-\lambda\eta_0/\mu) \text{ {\rm and } } a_1 < {\lambda}\eta_0/{\mu_1},\]

then, by using (28) and (29), we get that

\[\lim_{N\to+\infty} {\mathbb P}\left(\sup_{0\le s\leq t}\frac{L_{\eta_0}^N(s)}{N}<\eta_0, \sup_{0\le s\leq t}\frac{1}{N}(L_{\eta_0}^N(s)+Q(Ns)) < 1\right)=1.\]

This implies, in particular, the relation

(31) \begin{equation}\lim_{N\to+\infty} {\mathbb P}\left(\min\left(S_U^N,T_U^{\eta_0N}\right)>t\right)=1\end{equation}

for all $t \ge 0$ . Therefore, (28) and (29) hold with $(L_{\eta_0}^N(t))$ (resp. (Q(Nt))) replaced by $(U_1^N(t))$ (resp. $(U_0^N(t))$ ). Additionally, (29) shows that

\[\lim_{N\to+\infty} {\mathbb P}\left(T_U^0\leq \frac{\eta_0}{\lambda-(\mu_1-\mu_2)\eta_0 + \mu_2C}\right)=1.\]

With the same coupling as before and the strong Markov property, the process $(U_0^N(T_U^0+t), 0 \le t \le T)$ is bounded above by $(Q(Nt), 0 \le t \le T)$ , where (Q(t)) is the same $M/M/1$ process as before but starting at 0, $Q(0) = 0$ . Let

\[H_b\stackrel{{\text {def}.}} = \inf\{t \ge 0\,:\, Q(t)=b\}.\]

Proposition 5.11 of Robert [Reference Robert29] shows that, if $\rho\stackrel{{\text {def}.}} = (\mu_1-\mu_2)\eta_0 + \mu_2C)/\lambda < 1$ , then as b goes to infinity, the sequence of random variables $(\rho^bH_b)$ converges in distribution to an exponential distribution. Define

\[{\mathcal A}_N\stackrel{{\text {def}.}} = \left\{\sup_{T_U^0\le s\leq T_U^0+T} U_0^N(s) \leq C\log N\right\}\]

by choosing $C>-\log\rho$ ; one then has the relation

\[\liminf_{N\to+\infty}{\mathbb P}\left({\mathcal A}_N\right)\geq \lim_{N\to+\infty}{\mathbb P}\left(\sup_{0\le s\leq NT} Q(s) \leq C\log N\right)=1.\]

The proposition is proved.

Propositions 5 and 6 yield the following proposition.

Proposition 7. (Stability of saturation.) Under the condition $\rho_1 < C < \rho_2$ , there exist some $\eta_0>0$ and $t_0>0$ such that, if the initial condition is such that

(32) \begin{equation}\liminf_{N\to + \infty} \frac{1}{N}\sum_{i=1}^N Y_{i}^N(0) \geq C-\eta_0, \end{equation}

then, for any ${\varepsilon}>0$ and $T \ge 0$ ,

\[ \lim_{N\to+\infty} {\mathbb P}\left(\inf_{t_0\leq s\leq t_0 + T} \frac{1}{N}\sum_{i=1}^N Y_i^N(s) \geq C-{\varepsilon}\right)=1. \]

Proof. With the above notation,

\[\overline{Z}_1^N(t)\stackrel{{\text {def}.}} = \sum_{i=1}^N X_i^N(t),\quad \overline{Z}_2^N(t)\stackrel{{\text {def}.}} = \sum_{i=1}^N Y_i^N(t),\]

and $\overline{Z}_0^N(t) = CN-\overline{Z}_1^N(t)-\overline{Z}_2^N(t)$ . Without loss of generality, by taking a subsequence for example, we can assume that

\[\lim_{N\to + \infty} \frac{\overline{Z}_1^N(0)}{N}=a_1\text{ and }\lim_{N\to + \infty} \frac{\overline{Z}_0^N(0)}{N}=a_0,\]

for some $a_0$ , $a_1\in[0,C]$ . Proposition 5 gives a coupling of the process $(\overline{Z}_0^N(t),\overline{Z}_1^N(t))$ with the process $(U_0^N(t),U_1^N(t))$ , with the same initial conditions, such that the relation $U_2^N(t)\le\overline{Z}_2^N(t)$ holds for $t < S_U^N$ .

Proposition 6 shows that there exist $\eta_0>0$ and $t_0 \ge 0$ such that if (32) holds, then $a_0 + a_1\le\eta_0$ , so that, for $T>0$ ,

\[\lim_{N\to+\infty}{\mathbb P}\left(\sup_{t_0\le t\le t_0 + T}\frac{U_0^N(t)}{N}{\leq} {\varepsilon}, \sup_{t_0\le t\le t_0 + T}\frac{U_1^N(t)}{N}\leq {\varepsilon},\, S_U^N\geq t_0 + T\right)=1.\]

With the coupling, we obtain the relation

\[\lim_{N\to+\infty}{\mathbb P}\left(\inf_{t_0\le t\le t_0 + T}\frac{\overline{Z}_2^N(t)}{N}\geq C-2{\varepsilon}\right)=1.\]

The proposition is proved.

The next result gives a stability property for the limiting points of the sequence of empirical distributions $(\Lambda^N(t))$ .

Corollary 1. Under the assumptions of Proposition 7, any limiting point $(\Lambda(t))$ of the sequence of empirical processes $(\Lambda^N(t))$ of (10) satisfies the following relation for convergence in distribution: for any ${\varepsilon}>0$ there exists $t_0>0$ such that

\[ \Lambda(t)((0,C))>1-{\varepsilon} \qquad \forall t \ge t_0. \]

Proof. This is simply due to the above proposition and the fact that

\[ \frac{1}{N}\sum_{i=1}^N {\textbf{1}_{\left\{{Y_i^N(s){\ne}C}\right\}}} \leq \frac{1}{N}\sum_{i=1}^N (C-Y_i^N(s))= C-\frac{\overline{Z}_2^N(t)}{N}. \]

2.5. A spectral criterion of stability

In this section we investigate the stability properties of the nonlinear dynamical system $(\Lambda(t))$ defined by (12). The corresponding linear system $(\Lambda_R(t))$ is the solution of (14). The goal of this section is to show that if the linear process $(\Lambda_R(t))$ is converging sufficiently fast to equilibrium, the nonlinear process will converge to this invariant distribution provided its initial state is sufficiently close to it.

As before, $\pi_R$ is the probability distribution on ${\mathcal X}$ defined by (15); it is the invariant measure of $(\Lambda_R(t))$ . If R is a solution of Equation (16) (which was studied in Subsection 2.2), $\pi = \pi_R$ is an invariant measure of $(\Lambda(t))$ . Following Aldous and Fill [Reference Aldous and Fill1], the ‘distance’ between $\mu\in{\mathcal P}({\mathcal X})$ and a fixed probability $\pi$ is defined as

\[\|\mu-\pi\|_2^2\stackrel{{\text {def}.}} = \sum_{z\in{\mathcal X}} \left(\frac{\mu(z)}{\pi(z)}-1\right)^2 \pi(z)=\sum_{z\in{\mathcal X}}\frac{(\mu(z)-\pi(z))^2}{\pi(z)}.\]

Lemma 3.26 of Aldous and Fill [Reference Aldous and Fill1] shows that there exists a maximal $\kappa_R>0$ such that, for all $t \ge 0$ ,

(33) \begin{equation}\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} \|\Lambda_R(t)-\pi\|_2^2\leq -2\kappa_R\|\Lambda_R(t)-\pi\|_2^2.\end{equation}

This is the classical exponential convergence to equilibrium for finite Markov processes; the distance $\|\cdot\|_2$ gives the nice inequality (33) describing this phenomenon. The quantity $\kappa_R$ is the spectral gap of the process $(\Lambda_R(t))$ ; see Theorem 3.25 of Aldous and Fill [Reference Aldous and Fill1] for a variational characterization.

Theorem 1. Let R be a solution of Equation (16) and assume that the spectral gap $\kappa_R$ of $(\Lambda_R(t))$ satisfies the condition

(34) \begin{equation}\kappa_R>\frac{\lambda}{1-R}\sqrt{\frac{C}{\rho_2}}.\end{equation}

Then there exist positive constants q and ${\varepsilon}_0$ such that the relation

(35) \begin{equation}\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t}\|\Lambda(t)-\pi\|_2^2\le-q \|\Lambda(t)-\pi\|_2^2 \quad \text{ for all } t\ge 0\end{equation}

holds when $\|\Lambda(0)-\pi\|_2\le{\varepsilon}_0$ , with $\pi = \pi_R$ defined by (24). In particular, $\pi$ is an exponentially stable equilibrium point of $(\Lambda(t))$ .

Strategy of the proof. We first describe the main ideas. Our technique is based on a variation of an argument coming from the theory of attractors. The original approach consists of splitting the half-line ${\mathbb R}_+$ into intervals of large length and decomposing the dynamical system into two parts on each of them. The first one takes the linear part of the equation, for which there is exponential convergence. The second one includes the nonlinearity and is issued from zero on each interval. The convergence of the first part is then used to absorb the second part of the flow originated after the splitting, which allows one to obtain the desired result. We refer the reader to the paper of Zelik [Reference Zelik36], where this powerful approach is used in another context. Unfortunately, the direct application of this technique does not lead to good results in our case because of the strong nonlinearity in (12). To overcome this difficulty, we shall first prove an ‘instantaneous’ absorbing and then use a bootstrap argument to propagate; see Step 2 of the proof below.

Proof. Step1: Splitting of the flow. Let $\Lambda_R(t)$ be the solution of Equation (14) with a fixed initial point $\Lambda(0)$ and

\[ F(t)\stackrel{{\text {def}.}} = \|\Lambda(t)-\pi\|_2^2. \]

For $t \ge 0$ , we have

\[ {F}^{\prime}(t)=\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t}\sum_{z\in{\mathcal X}} \frac{[\Lambda(t)-\pi](z)^2}{\pi{(z)}}=2\sum_{z\in{\mathcal X}} \frac{[\Lambda(t)-\pi](z) \Lambda'(t)(z)}{\pi{(z)}}. \]

We will now state two claims and show how they are used to establish our theorem. The proofs of the claims conclude the proof of the theorem.

Claim 1. The relation

(36) \begin{equation}{\Lambda}^{\prime}(0)(z)={\Lambda}^{\prime}_R(0)(z)+\lambda L_{\Lambda, R}\left(\rule{0mm}{4mm}\Lambda(0)(z-e_2)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}}\Lambda(0)(z)\right)\end{equation}

holds, where

(37) \begin{equation}L_{\Lambda, R} = \frac{\Lambda(0)({\mathcal X}_+^c)}{1-\Lambda(0)({\mathcal X}_+^c)}-\frac{R}{1-R}.\end{equation}

Taking $t = 0$ and using that $\Lambda_R(0) = \Lambda(0)$ together with (36), we obtain

(38) \begin{align}{F}^{\prime}(0) = & 2\sum_{z\in{\mathcal X}} \frac{([\Lambda_R(0)-\pi](z)){\Lambda}^{\prime}_R(0)(z)}{\pi(z)}\nonumber\\ & + 2\lambda L_{\Lambda, R}\sum_{z\in{\mathcal X}} \frac{([\Lambda(0)-\pi](z))(\Lambda(0)(z-e_2)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}}\Lambda(0)(z))}{\pi(z)}={\mathcal I}_1 + {\mathcal I}_2.\end{align}

Note that the term ${\mathcal I}_1$ is the derivative of the function $t{\mapsto}\|\Lambda_R(t)-\pi\|_2^2$ at zero. From Inequality (33), we get therefore that

(39) \begin{equation}\frac{{\mathcal I}_1}{2}\le -\kappa_R\|\Lambda_R(0)-\pi\|_2^2=-\kappa_R\|\Lambda(0)-\pi\|_2^2=-\kappa_R F(0).\end{equation}

Claim 2. There exist ${\varepsilon}>0$ and $\theta < 2\kappa_R$ such that, for $\|\Lambda(0)-\pi\|_2\le {\varepsilon}$ ,

(40) \begin{equation} {\mathcal I}_2\le \theta F(0).\end{equation}

Combining this with (38) and (39), we get

(41) \begin{equation}{F}^{\prime}(0)\le -q F(0)\end{equation}

for $\|\Lambda(0)-\pi\|_2\le {\varepsilon}$ and $q = 2\kappa_R-\theta>0$ .

Step2: Bootstrap argument. First note that we may suppose that $F(0)>0$ . Otherwise $\Lambda(0) = \pi$ and therefore $\Lambda(t) = \pi$ for all $t \ge 0$ ; in this case there is nothing to prove. Let

\[\tau\stackrel{{\text {def}.}} = \inf\left\{s>0\,:\, F(s)>{\varepsilon}^2 \right\},\]

with the convention that $\inf\emptyset = + \infty$ . Since $0 < F(0)\le{\varepsilon}^2$ , we have by (41) that $F^{\prime}(0)\le-qF(0) < 0$ , so that $\tau>0$ .

The constants ${\varepsilon}$ and q in Inequality (41) do not depend on the initial point $\Lambda(0)$ as long as $\|\Lambda(0)-\pi\|_2\le {\varepsilon}$ . Hence we can apply the same argument considering the equation as starting from $\Lambda(s)$ , for $s \ge 0$ , as long as $\|\Lambda(s)-\pi\|_2\le {\varepsilon}$ ; we can thus infer the differential inequality with 0 replaced by s. We have therefore the relation

\[F^{\prime}(s)\le -qF(s);\]

hence $F(s)\le F(0)\exp(-qs)$ for all $0\le s < \tau$ . This implies that $\tau$ is infinite. Inequality (35) is established. It remains to prove our two claims, i.e. Equation (36) and Inequality (40).

Step3: Proof of the identity (36).

From now on, z denotes a generic element $(x,y)\in{\mathcal X}$ of the state space, and $e_1 = (1,0)$ , $e_2 = (0,1)$ are the unit vectors of ${\mathcal X}$ .

From Equation (12), we get

(42) \begin{align} \Lambda^{\prime}(t)(z) = & \lambda\left(\rule{0mm}{4mm}\Lambda(t)(z-e_2)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}}\Lambda(t)(z)\right)\frac{\Lambda(t)({\mathcal X}_+^c)}{1-\Lambda(t)({\mathcal X}_+^c)}\nonumber\\ & +\lambda\left(\rule{0mm}{4mm}\Lambda(t)(z-e_1)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}} \Lambda(t)(z)\right)+\mu_1\left(\rule{0mm}{4mm}(x + 1)\Lambda(t)(z + e_1)-x\Lambda(t)(z)\right)\nonumber\\ & +\mu_2\left(\rule{0mm}{4mm}(y + 1)\Lambda(t)(z + e_2)-y\Lambda(t)(z)\right),\end{align}

with the convention $\Lambda(t)(z^{\prime}) = 0$ if $z^{\prime}{\notin}{\mathcal X}$ .

Similarly, the definition (14) gives

(43) \begin{align} \Lambda^{\prime}_R(t)(z) = & \lambda\left(\rule{0mm}{4mm}\Lambda_R(t)(z-e_2)-{\textbf{1}_{\left\{{{\mathcal X}_+}\right\}}}(z)\Lambda_R(t)(z)\right)\frac{R}{1-R}\nonumber\\ & + \lambda\left(\rule{0mm}{4mm}\Lambda_R(t)(z-e_1)-{\textbf{1}_{\left\{{{\mathcal X}_+}\right\}}}(z)\Lambda_R(t)(z)\right)+\mu_1\left(\rule{0mm}{4mm}(x+1)\Lambda_R(t)(z + e_1)-x\Lambda_R(t)(z)\right)\nonumber\\ & + \mu_2\left(\rule{0mm}{4mm}(y+1)\Lambda_R(t)(z + e_2)-y\Lambda_R(t)(z)\right).\end{align}

Taking $t = 0$ in these two relations and using that $(\Lambda_R(t))$ and $(\Lambda(t))$ have the same initial conditions, we get the identity (36). We now establish the most intricate inequality of our theorem, namely Inequality (40).

Step4: A bound for ${\mathcal I}_2$ . From the expression (37) of $L_{\Lambda, R}$ , we get

\[|L_{\Lambda, R}|\le \frac{1}{(1-R)^2}\frac{|\Lambda(0)({\mathcal X}_+^c)-R|}{1-|\Lambda(0)({\mathcal X}_+^c)-R|/(1-R)},\]

and, by the Cauchy-Schwartz inequality,

\begin{align*}|& \Lambda(0)({\mathcal X}_+^c)-R| = \left|\sum_{z\in{\mathcal X}_+^c}[\Lambda(0)-\pi](z)\right|= \left|\sum_{z\in{\mathcal X}_+^c}\sqrt{\pi(z)}\frac{[\Lambda(0)-\pi](z)}{\sqrt{\pi(z)}}\right|\\[5pt] & \le \left(\sum_{z\in{\mathcal X}_+^c}\pi(z)\right)^{1/2}\left(\sum_{z\in{\mathcal X}_+^c} \frac{([\Lambda(0)-\pi](z))^2}{\pi(z)}\right)^{1/2}=\sqrt{R}\|\Lambda(0)-\pi\|_2=\sqrt{RF(0)}.\end{align*}

Combining these two relations, we get the inequality

(44) \begin{align}|L_{\Lambda, R}|\le \frac{1}{(1-R)^2}\frac{\sqrt{R F(0)}}{1 - \sqrt{R F(0)}/(1-R)}.\end{align}

Since

(45) \begin{align} & \sum_{z\in{\mathcal X}} \frac{[\Lambda(0)-\pi](z)[\Lambda(0)(z-e_2)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}}\Lambda(0)(z)]}{\pi(z)}\nonumber\\[4pt] &\quad =\sum_{z\in{\mathcal X},y \ge 1} \frac{[\Lambda(0)-\pi](z)[\Lambda(0)-\pi](z-e_2)}{\pi(z)} + \sum_{z\in{\mathcal X}} \frac{[\Lambda(0)-\pi](z)[\pi(z-e_2)-{\textbf{1}_{\left\{{z\in{\mathcal X}_+}\right\}}}\pi(z)]}{\pi(z)}\nonumber\\[6pt] & \qquad + \sum_{z\in{\mathcal X}_+} \frac{[\Lambda(0)-\pi](z)[\pi-\Lambda(0)](z)}{\pi(z)}\stackrel{{\text {def}.}} = {\mathcal J}_1 + {\mathcal J}_2 + {\mathcal J}_3,\end{align}

we obtain from (44) and (45) that

(46) \begin{equation}{\mathcal I}_2\le 2\lambda |L_{\Lambda,R}|\left(|{\mathcal J}_1| + |{\mathcal J}_2| + |{\mathcal J}_3|\right).\end{equation}

Step5: Estimates for ${\mathcal J}_i$ and a final bound for ${\mathcal I}_2$ .

Clearly $|{\mathcal J}_3|\le \|\Lambda(0)-\pi\|_2^2$ . By using (15) and the Cauchy-Schwartz inequality, we get

(47) \begin{align} |{\mathcal J}_1|&\le \max_{z\in{\mathcal X}}{\sqrt{\frac{\pi(z-e_2)}{\pi(z)}}}\sum_{z\in{\mathcal X}, y\ge 1} \frac{|[\Lambda(t)-\pi](z)|}{\sqrt{\pi(z)}}\frac{|[\Lambda(0)-\pi](z-e_2)|}{\sqrt{\pi(z-e_2)}}\nonumber\\[4pt] & \le{\sqrt{\frac{C}{R\rho_2}}} \left(\sum_{z\in{\mathcal X}, y \ge 1} \frac{([\Lambda(0)-\pi](z))^2}{\pi(z)}\right)^{1/2}\left(\sum_{z\in{\mathcal X}, y\ge 1} \frac{([\Lambda(0)-\pi](z-e_2))^2}{\pi(z-e_2)}\right)^{1/2}\nonumber\\[4pt] & \le {\sqrt{\frac{C}{R\rho_2}}}\|\Lambda(0)-\pi\|_2^2={\sqrt{\frac{C}{R\rho_2}}}F(0).\end{align}

Another application of the Cauchy-Schwartz inequality gives

\[|{\mathcal J}_2|\le \left(\sum_{z\in{\mathcal X}}\frac{(\pi(z-e_2)-\textbf{1}_{{\mathcal X}_+}(z)\pi(z))^2}{{\pi(z)}}\right)^{1/2}\|\Lambda(0)-\pi\|_2.\]

On the other hand, by using (15) and (16), we obtain

\begin{align*}&\sum_{z\in{\mathcal X}}\frac{(\pi(z-e_2)-\textbf{1}_{{\mathcal X}_+}(z)\pi(z))^2}{{\pi(z)}}\le \sum_{z\in{\mathcal X}}\frac{(\pi(z-e_2))^2}{{\pi(z)}}\\&\quad = \sum_{z\in{\mathcal X}}\pi(z-e_2)\frac{y(1-R)}{R\rho_2}\le \frac{C(1-R)^2}{R\rho_2}, \end{align*}

so that

\[ |{\mathcal J}_2|\le (1-R)\left(\frac{C F(0)}{R\rho_2}\right)^{1/2}.\]

Inequalities (44) and (46) imply

\begin{align*} & {\mathcal I}_2\le 2\lambda \frac{1}{(1-R)^2}\frac{\sqrt{R F(0)}}{1 - \sqrt{R F(0)}/(1-R)} \\[4pt] & \qquad\times\left((1-R)\left(\frac{C F(0)}{R\rho_2}\right)^{1/2} +\left({\sqrt{\frac{C}{R\rho_2}}} + 1\right)F(0)\right); \end{align*}

consequently,

\[ \limsup_{{\varepsilon}{\to}0} \sup_{0 < F(0)\le {\varepsilon}^2} \frac{{\mathcal I}_2}{F(0)} \leq \frac{2\lambda}{1-R}\sqrt{\frac{C}{\rho 2}}<2\kappa_R\]

by the assumption (34). Inequality (40) is thus established. The theorem is proved.

3.1. The basic ODEs

We recall briefly the technical background for this algorithm. (See Gibbens et al. [Reference Gibbens, Hunt and Kelly16].) There is a set V of vertices, and for each couple $(A,B)=(B,A)$ of vertices, referred to as a link or node, there is a Poisson flow of requests with rate $\lambda$ , referred to as calls or jobs, to establish a connection between A and B. The capacity constraint is that there can be at most C jobs at a given time on any node. The state of the process is given by the number of jobs in the nodes of the network. The algorithm works as follows. If a node (A, B) has strictly fewer than C jobs at some instant, then a request arriving at (A, B) is accepted. Otherwise, a vertex $D{\not\in}\{A,B\}$ is chosen at random, and if both links (A, D) and (D, B) have strictly fewer than C jobs, then the job occupies a place in (A, D) and in (D, B) during an exponentially distributed amount of time with parameter 1. If one of the links (A, D), (D, B) is saturated, the request is rejected.

The mean-field result of Graham and Méléard [Reference Graham and Méléard17], described below, shows that from the point of view of the convergence of the empirical distribution process $(\Lambda^N(t))$ defined by Equation (1), the DAR algorithm has the same limiting behavior as the following allocation algorithm. There is a set of $N = |V|(|V|-1)/2$ nodes with finite capacity C; each node receives a Poisson flow of jobs with parameter $\lambda$ to be processed at rate 1. When a request arrives at a saturated node $\ell$ , two other nodes are chosen at random. If neither of them is saturated, a new request is added to each of them. Otherwise the initial request to node $\ell$ is rejected.

For $1{\leq}\ell\le N$ , $L_{\ell}^N(t)$ denotes the number of jobs in node $\ell$ at time $t \ge 0$ . Note that $(L_{\ell}^N(t),1\le\ell\le N)$ is not a Markov process. The mean-field result of Graham and Méléard [Reference Graham and Méléard17] (which was conjectured by Gibbens et al. [Reference Gibbens, Hunt and Kelly16]) is as follows. The initial state is given by i.i.d. random variables with some distribution $\mu$ on $\{0,\ldots,C\}$ for the $L_{\ell}^N(0)$ , $1\le\ell\le N$ , and without any request using two nodes. It has been shown in [Reference Graham and Méléard17] that the convergence in distribution

(48) \begin{equation} \lim_{N\to+\infty} (L_{\ell}^N(t))= (\overline{L}_C(t))\end{equation}

holds for any $\ell \ge 1$ , where, for $t \ge 0$ , the distribution of $\overline{L}_C(t)$ , given by the vector $({\mathbb P}(\overline{L}_C(t) = k), 0\le k\le C) = (x^C_{k}(t),0\le k\le C)$ , is the solution of the following ODEs: for $1{\leq}j < C$ ,

(49) \begin{equation}\begin{cases}\displaystyle \frac{\mathop{}\mathopen{}\mathrm{d} x^C_{0}}{\mathop{}\mathopen{}\mathrm{d} t}(t) \displaystyle = x^C_{1}(t)- \lambda h\left(x^C_{C}(t)\right)x^C_{0}(t),\\[10pt] \displaystyle \frac{\mathop{}\mathopen{}\mathrm{d} x^C_{j}}{\mathop{}\mathopen{}\mathrm{d} t}(t) \displaystyle = \lambda h\left(x^C_{C}(t)\right)x^C_{{j-1}}(t)+ (j + 1) x^C_{{j+1}}(t)-( \lambda h\left(x^C_{C}(t)\right) + j )x^C_{j}(t),\\[10pt] \displaystyle \frac{\mathop{}\mathopen{}\mathrm{d} x^C_{C}}{\mathop{}\mathopen{}\mathrm{d} t}(t)\displaystyle = \lambda h\left(x^C_{C}(t)\right)x^C_{{C-1}}(t)- C x^C_{C}(t),\end{cases}\end{equation}

where $h(x)\stackrel{{\text {def}.}} = (1 + 2 x(1-x))$ , with initial condition $(x^C_k(0)) = (\mu(k))$ .

We give an intuitive explanation of this system, to explain the role of the function h in particular. For $t \ge 0$ , $x_C^C(t)$ is the fraction of saturated nodes, i.e. the number of nodes with C jobs. The mean-field limit expresses an asymptotic independence property: the numbers of jobs at a fixed finite subset of nodes are, in the limit, independent, and $x_C^C(t)$ is the probability that an arbitrary node is saturated. Each non-saturated node accepts external requests arriving at rate $\lambda$ . It may also be occupied with a request which has arrived at a saturated node, if this request has picked (at random) this empty node and another non-saturated node. With the independence approximation, this occurs with probability

\[2x_C^C(t)(1-x_C^C(t)).\]

The equations of the system (49) can then easily be explained.

The equations (49) can be seen as the set of Fokker-Planck equations for a non-homogeneous $M/M/C/C$ queue for which the arrival rate at time t is $\lambda h\left(x^C_{C}(t)\right)$ and the service rate is 1. Equivalently, from a probabilistic point of view, the process $(\overline{L}_C(t))$ has the same distribution as the solution of the following McKean-Vlasov SDE:

(50) \begin{align}\mathop{}\mathopen{}\mathrm{d} \overline{L}_C(t) = & {\textbf{1}_{\left\{{\overline{L}_C(s-) < C}\right\}}}{\mathcal P}_1\left(\left[0,\lambda h\left({\mathbb P}\left(\overline{L}_C(t) = C\right)\right)\right] \times \mathop{}\mathopen{}\mathrm{d} t\right)\nonumber\\[5pt] & -{\mathcal P}_2\left(\left[0,\overline{L}_C(t-)\right] \times \mathop{}\mathopen{}\mathrm{d} t\right),\end{align}

with $\overline{L}_C(0)\stackrel{\text{\rm {dist}.}} = \mu$ , where ${\mathcal P}_i$ , $i = 1$ , 2, are independent Poisson processes on ${\mathbb R}_+^2$ with rate 1.

An equilibrium point $(x_k^*)$ of the dynamical system (49) is given by

\[x_k^*=\frac{1}{Z_C} \frac{\left(\lambda h(x_C^*)\right)^k}{k!}, \quad 0\le k < C,\]

where $Z_C$ is the normalization constant and $z = x_C^*$ is a positive solution of the fixed point equation

(51) \begin{equation}\frac{\left(\lambda h(z)\right)^C}{C!}-z\sum_{k=0}^C \frac{\left(\lambda h(z)\right)^k}{k!}=0.\end{equation}

There always exists a solution to this equation, since the left-hand side of (51) is positive for $z = 0$ and negative for $z = 1$ . The rest of this section is devoted to determining the conditions under which there may exist several solutions for this equation and to investigating their stability properties for the dynamical system (49).

The insightful Gibbens et al. [Reference Gibbens, Hunt and Kelly16] suggests, through approximations and numerical experiments, that this equation may have in fact several solutions:

(The notation has been adapted.) To the best of our knowledge, these statements do not seem to have been established in a more formal way. The rest of this section is devoted to a scaling analysis of the set of ODEs (49). As suggested by these numerical experiments, we will study the case of a large capacity C. Concerning the stability results of this assertion, we have not been able to prove them as such. Subsections 3.3 and 3.4 give only partial results in this domain.

3.2. An asymptotic dynamical system

We denote by $(x^C_k(t))$ the solution of the ODE (49) when $\lambda$ is replaced by $\lambda C$ and $h(x) = (1 + ax(1-x))$ for some $a>1$ and $x\in(0,1)$ . For this algorithm too, there is a kind of analogue of the regime analyzed in Subsection 2.4, in the sense that it has some intuitive explanation. In this regime, in the limit, all jobs are accommodated without rerouting provided that the initial state is not already saturated. As will be seen, for the same parameters, there are nevertheless two other distinct equilibrium regimes for which a job is rerouted with positive probability as C gets large.

1. 1. $\nu\in(0,1)$ , then there exists a solution $x_{C,1}^*\in(0,1)$ of Equation (51) such that $\nu h(x_{C,1}^*) < 1$ , and

   (52) \begin{equation} \lim_{C\to+\infty} x_{C,1}^*=0.\end{equation}

2. 2. If $\nu\in(\nu_a,1)$ , with

   (53) \begin{equation}\nu_a\stackrel{{\text {def}.}} = {3}\left/\left(1 + \frac{2}{9}a + \frac{2}{3}\frac{(a + 3)^{3/2}}{\sqrt{a}}\right),\right.\end{equation}
   then there are three solutions $x_{C,i}^*\in(0,1)$ , $i\in\{1,2,3\}$ , of Equation (51), with the following properties: $\nu h(x_{C,1}^*) < 1$ , and $\nu h(x_{C,i}^*)>1$ for $i = 2$ , 3. The limiting values of $(x_{C,2}^*)$ and $(x_{C,3}^*)$ are the two solutions of the polynomial equation
   (54) \begin{equation}a{z}^{3}-2a{z}^{2} + (a-1)z + 1=\frac{1}{\nu}\end{equation}
   in (0,1), and $(x_{C,1}^*)$ satisfies (52).

3. 3. If $\nu>1$ , there exists a unique solution $x_{C,1}^*\in(0,1)$ of Equation (51), and the sequence $(x_{C,1}^*)$ converges to the unique solution of Equation (54) in the interval (0, 1).

For the Gibbens et al. model, which corresponds to the case $a = 2$ , this gives $\nu_2{\sim}0.937$ ; hence $(0.937,1)$ is the ‘narrow range of the ratio $\lambda/C$ ’ suggested by these authors (see page 41) for which there are three solutions to the fixed point equation.

According to Part (1), when $\nu < 1$ , there is an equilibrium in the light-load regime ( $x_{C,1}^*{\sim}0$ ). We will see a more precise result, Proposition 9, concerning the asymptotic local stability of this equilibrium. When $\nu_a < \nu < 1$ , there are two other equilibrium points, but in a saturated regime: $x_{C,i}^*\geq\eta>0$ for $i = 2$ , 3 and C sufficiently large.

Proof. The fixed point equation (51) with $\lambda$ replaced by $\nu C$ can be rewritten as $\Psi_{C}(z) = 0$ , with, for $z>0$ ,

\[\Psi_C(z)\stackrel{{\text {def}.}} = 1-z\sum_{k=0}^C \frac{C!}{(C-k)!}\frac{1}{\left(\nu C h(z)\right)^k}=1-z\sum_{k=0}^C \prod_{i=0}^{k-1}\left(1-\frac{i}{C}\right) \frac{1}{\left(\nu h(z)\right)^k}.\]

Note that the function $C{\mapsto}\Psi_C(z)$ is decreasing.

If $\nu < 1$ , we can choose ${\varepsilon} < 1/2$ sufficiently small so that $\delta = \nu h({\varepsilon}) < 1/2$ holds. For $z\le{\varepsilon}$ , we have

\[\sum_{k=0}^C \prod_{i=0}^{k-1}\left(1-\frac{i}{C}\right) \frac{1}{\left(\nu h(z)\right))^k}\ge\sum_{k=0}^{\lfloor \delta C\rfloor} \prod_{i=0}^{k-1}\left(1-\frac{i}{C}\right) \frac{1}{\delta^k}\ge\sum_{k=0}^{\lfloor \delta C\rfloor} \left(\frac{1-\delta}{\delta}\right)^k.\]

In particular, for $0 < z < {\varepsilon}$ , $(\Psi_C(z))$ converges to $-\infty$ as C goes to infinity. Since $\Psi_C(0) = 1$ , one can find $C_0>0$ such that if $C \ge C_0$ , then there is a zero of $\Psi_C$ in the interval $(0,{\varepsilon})$ .

For $\delta>1$ , it is easily checked that the convergence

(55) \begin{equation}\lim_{C\to+\infty} \Psi_{C}(z) =\overline{\Psi}(z)\stackrel{{\text {def}.}} = 1-z\frac{1}{1-1/(\nu h(z))}=\frac{\nu(1-z)h(z)-1}{\nu h(z)-1}\end{equation}

holds uniformly for all $z\in(0,1)$ such that $\nu h(z)>\delta$ . Equation (51) becomes, in the limit,

(56) \begin{equation}(1-z)h(z)=\frac{1}{\nu}.\end{equation}

Note that any solution $z < 1$ of such an equation satisfies $\nu h(z)>1$ .

The quantity $(1-z)h(z)$ is the polynomial $a{z}^{3}-2a{z}^{2} + az-z + 1$ . This polynomial is increasing from 1 on the interval $[0,x_0]$ and decreasing on $[x_0,1]$ , where

\[x_0\stackrel{{\text {def}.}} = \frac{1}{3}\left(2-\sqrt{\frac{a + 3}{a}}\right),\quad (1-x_0)h(x_0) = 1{/}\nu_a = \frac{1}{3}\left(1 + \frac{2}{9}a + \frac{2}{3}\frac{(a + 3)^{3/2}}{\sqrt{a}}\right).\]

Equation (56) therefore has two solutions if and only if $\nu\in(\nu_a,1)$ , one solution when $\nu = \nu_a$ or $\nu>1$ , and none if $\nu < \nu_a$ .

If $\nu\in(\nu_a,1)$ , for ${\varepsilon}>0$ sufficiently small, there exist $z_0 < z_1 < z_2$ such that $\nu h(z_i)>1$ for $i\in\{0,1,2\}$ , and $\overline{\Psi}(z_0) < -{\varepsilon}$ , $\overline{\Psi}(z_1)>{\varepsilon}$ and $\overline{\Psi}(z_2) < -{\varepsilon}$ . Consequently, there exists $K_0$ such that, if $C \ge K_0$ , then the last three inequalities hold with $\overline{\Psi}$ replaced by $\Psi_C$ and ${\varepsilon}$ by ${\varepsilon}/2$ . Hence, we get that there are two solutions of the relation $\Psi_C(z) = 0$ such that $\nu h(z)>1$ for $C \ge K_0$ . Assertion (2) is proved.

If $\nu>1$ , the convergence (55) is uniform for $z\in[0,1]$ . Since $\overline{\Psi}(0)>0$ and $\overline{\Psi}(1) < 0$ , by the same argument as before, there exists some $K_1$ such that if $C \ge K_1$ then there is a solution $x_{C,1}^*$ of Equation (51). A simple calculation gives

\[\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} z} \overline{\Psi}(z)=-\frac{\nu}{(\nu h(z)-1)^2}\left(\nu{z}^{2} \left( z-1 \right) ^{2}{a}^{2}-{z}\left( 2\nu (z-1)-3z + 2 \right) a + \nu-1\right),\]

and one has

\begin{align*} &\nu{z}^{2}\left( z-1 \right) ^{2}{a}^{2}-{z}\left( 2\nu(z-1)-3z + 2 \right) a\geq az (3z-2+ \left( {z}^{3}-2{z}^{2}-z + 2 \right) \nu)\\ &\quad \ge az (3z-2+ \left( {z}^{3}-2{z}^{2}-z + 2 \right))=az^2(z^2-2z + 2)\ge 0\end{align*}

by using successively that $a \ge 1$ , that ${z}^{3}-2{z}^{2}-z + 2 = (1-z)(2 + 2z-z^2) \ge 0$ when $z\in[0,1]$ , and finally that $\nu \ge 1$ . Consequently, we get that

(57) \begin{equation}\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} z} \overline{\Psi}(z)\leq-\frac{\nu(\nu-1)}{(\nu h(z)-1)^2}\leq -\frac{\nu(\nu-1)}{(\nu (1 + a/4)-1)^2}<0\end{equation}

holds for all $z\in[0,1]$ .

Note that, as in (55), since $\nu>1$ , the convergence

\[\lim_{C{\to} + \infty} \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} z} {\Psi}_C(z)=\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} z} \overline{\Psi}(z)\]

holds uniformly for $z\in[0,1]$ .

To prove uniqueness, we assume that there is a subsequence $(C_n)$ converging to infinity for which the equation $\Psi_{C_n}(z) = 0$ has two solutions. This implies in particular that we have a sequence $(z_n)$ in (0,1) such that ${\Psi}^{\prime}_{C_n}(z_n) = 0$ . Due to the uniform convergence, this is in contradiction to (57). The proposition is proved.

When $\nu < 1$ , with our method based on the asymptotic behavior of $\Psi_C(z)$ as C gets large, we have not been able to prove that all solutions of Equation (51) are identified, though this is very likely the case.

A scaled version of the dynamical system

For the moment we have given a scaled version of the fixed point equations. It turns out that one can also get some insight from a scaled version of the dynamical system (49) converging to a nontrivial dynamical system whose fixed points are described in Theorem 2.

Let us introduce some notation. The set of bounded sequences is denoted by ${\mathcal B}({\mathbb N})$ , and it is endowed with the following norm: for $z\in{\mathcal B}({\mathbb N})$ ,

\[\|z\|=\sum_{k=0}^{+\infty}\frac{1}{2^k} |z_k|.\]

For $T>0$ and $(z(t)) = (z_k(t))\in{\mathcal C}({\mathbb R}_+,{\mathcal B}({\mathbb N}))$ a continuous function on ${\mathbb R}_+$ , we define

\[\|z\|_{T}\stackrel{{\text {def}.}} = \sup_{0\,{\le}\,t\,{\le}\,T} \|z(t)\|.\]

Additionally, ${\mathcal P}({\mathbb N})$ is the set of probability distributions on ${\mathbb N}$ .

The scaling consists of slowing down the time scale by a factor C and looking at the number of empty places for the McKean-Vlasov process.

Proposition 8. (Asymptotic dynamical system.) If $(x_k^C(t))$ is the solution of the set of ODEs defined by (49) with $\lambda = \nu C$ , with an initial point such that

\[\lim_{C\to+\infty}\|(x^C_{C-k}(0), k\in{\mathbb N})-(q_0(k),k\in{\mathbb N})\|=0\]

for some probability distribution $q_0\in{\mathcal P}({\mathbb N})$ , then, as C goes to infinity, the process

\[ (y_k^C(t),k\in{\mathbb N})\stackrel{{\text {def}.}} = \left(x^C_{C-k}\left({t}/{C}\right), k\in{\mathbb N}\right) \]

converges in distribution for the uniform norm $\|{\cdot}\|_T$ to $(\Gamma(t))\in{\mathcal C}({\mathbb R}_+,{\mathcal P}({\mathbb N}))$ , which is the unique solution of the set of differential equations

(58) \begin{equation}\begin{cases}\displaystyle \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} \Gamma_{0}(t) = \nu h\left(\Gamma_{0}(t)\right)\Gamma_{{1}}(t)- \Gamma_{0}(t),\\[8pt] \displaystyle \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} \Gamma_{k}(t) = \nu h\left(\Gamma_{0}(t)\right)\Gamma_{{k+1}}(t) + \Gamma_{{k-1}}(t)-\left( \nu h\left(\Gamma_{0}(t)\right) + 1\right)\Gamma_{k}(t),\quad k \ge 1,\end{cases} \end{equation}

with $\Gamma(0) = q_0$ .

Proof. Note that the process $(y_k^C(t))$ can be seen as a version of the empirical distribution process of empty places in the nodes of the network with the slowed-down time scale $t{\mapsto}t/C$ .

It is not difficult to check that $(y_k^C(t))$ satisfies the following system of ODEs:

(59) \begin{equation}\begin{cases}\displaystyle \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} y^C_{0}(t) = \nu h\left(y^C_{0}(t)\right)y^C_{{1}}(t)- y^C_{0}(t),\\[12pt] \displaystyle \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} y^C_{k}(t) = \nu h\left(y^C_{0}(t)\right)y^C_{{k+1}}(t) + \left(1-\frac{k + 1}{C}\right) y^C_{{k-1}}(t)\\[12pt] \displaystyle -\left( \nu h\left(y^C_{0}(t)\right) + \left(1-\frac{k}{C}\right)\right)y^C_{k}(t),\\[12pt] \displaystyle \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t} y^C_{C}(t) = \frac{1}{C}y^C_{1}(t)- \nu h\left(y^C_{0}(t)\right)y^C_{C}(t).\end{cases}\end{equation}

The proposition is proved using a classical consequence of the compactness-uniqueness argument: the sequence of functions $(y_k^C(t))$ , $C \ge 1$ , is tight, and any limiting point satisfies the equations (59). If this system of equations has a unique solution, then the convergence is established, since all subsequences have a subsequence converging to this limit. (See Chapter 10 and 11 of Ethier and Kurtz [Reference Ethier and Kurtz14], for example.)

The tightness is due to the Arzelà-Ascoli theorem (see Theorem 7.2 of Billingsley [Reference Billingsley3], for example). The equations (59) show that for $k\in{\mathbb N}$ , the sequence of functions $(y^C_k(t))$ is equicontinuous and therefore is relatively compact for the uniform norm on bounded intervals. The integral form of the ODEs (59) shows that any limiting point satisfies (58).

Let $(x_k(t))$ and $(y_k(t))$ be two solutions of (58) with the same initial condition. Using the fact that the function h is Lipschitz on [0, 1] with parameter 6, the integral form of (58) gives the following inequality: for $(u_k(t)) = (x_k(t)-y_k(t))$ and $t>0$ ,

\begin{align*}|u_k(t)|\leq & 6\nu\int_0^t |u_0(s)|(x_k(s) + x_{k+1}(s))\,\mathop{}\mathopen{}\mathrm{d} s \\ & +6\nu \int_0^t (|u_{k-1}(s)|+|u_{k}(s)| + |u_{k+1}(s)|)\,\mathop{}\mathopen{}\mathrm{d} s.\end{align*}

Since $(x_k(t))$ and $(y_k(t))$ are probability distributions on ${\mathbb N}$ , $(u_k(t))$ is a convergent series and

\[\sum_{k=0}^n |u_k(t)|\leq 30\nu \int_0^t \sum_{k=0}^{n+1} |u_k(t)|\,\mathop{}\mathopen{}\mathrm{d} s,\]

so that

\[U(t)\stackrel{{\text {def}.}} = \sum_{k=0}^{+\infty} |u_k(t)|\leq 30\nu \int_0^t U(s)\,\mathop{}\mathopen{}\mathrm{d} s.\]

Grönwall’s inequality gives the relation $(U(t)) = (0)$ , i.e. the uniqueness of the solution of (58), and, consequently, the desired convergence.

A probabilistic translation of this result can be stated as the fact that if the process $(\overline{L}_C(t))$ defined by (50) satisfies the relation

\[\lim_{C\to+\infty} C-\overline{L}_C(0)=q_0\]

for convergence in distribution, then

\[\lim_{N\to+\infty} (Q_C(t))\stackrel{{\text {def}.}} = (C-\overline{L}_C(t/C))=(\overline{Q}(t)),\]

where $(\overline{Q}(t))$ is the solution of the McKean-Vlasov SDE

(60) \begin{equation} \mathop{}\mathopen{}\mathrm{d} \overline{Q}(t)={\mathcal N}_1\left(\mathop{}\mathopen{}\mathrm{d} t\right)-{\textbf{1}_{\left\{{\overline{Q}(s-)>0}\right\}}}{\mathcal P}_1\left(\left[0,\nu h\left({\mathbb P}\left(\overline{Q}(t) = 0\right)\right)\right] \times \mathop{}\mathopen{}\mathrm{d} t\right),\end{equation}

with $\overline{Q}(0)\stackrel{\text{\rm {dist}.}} = q_0$ , where ${\mathcal N}_1$ and ${\mathcal P}_1$ are independent homogeneous Poisson processes with rate 1 on ${\mathbb R}_+$ and ${\mathbb R}_+^2$ respectively.

The process $(\overline{Q}(t))$ is a nonlinear $M/M/1$ queue with the following jump rates at time t:

(61) \begin{equation}\begin{cases} + 1& 1\\ -1& \nu h\left({\mathbb P}\left(\overline{Q}(t) = 0\right)\right).\end{cases}\end{equation}

It should be noted that this scaling is convenient for studying regimes where the number of empty places is small, i.e. when the system has some level of saturation. Subsection 3.3 studies the case when there is an equilibrium regime with a large number of empty places.

Remark 1. An invariant distribution $\overline{\pi}$ of the nonlinear Markov process $(\overline{Q}(t))$ defined by (60) is the invariant distribution of an $M/M/1$ queue with arrival rate 1 and service rate $\nu h(\overline{\pi}(0))$ . Thus $\overline{\pi}$ is a geometric distribution with parameter $1/(\nu h(\overline{\pi}(0)))$ ; in particular,

\[ \overline{\pi}(0)=1-\frac{1}{\nu h(\overline{\pi}(0))}, \]

which is Equation (56), as can be expected. Indeed, Equation (56) is the limit, as C gets large, of the fixed point equation (51) for $x_{C}$ , which is the equilibrium probability that a node has zero empty places.

3.4 Stability of the underloaded regime

In Theorem 2 we have seen that, under the condition $\nu < 1$ , there is a root of the fixed point equation (51) that is arbitrarily close to 0 as C gets large. This result suggests that the underloaded regime is stable for the dynamical system (49). In this regime most requests are accepted at their arrival node. The following proposition gives a formal characterization of this property.

Proposition 9. (Stability of underloaded regime.) If $\lambda = \nu C$ for some $\nu < 1$ , there exists $\eta\in(0,1)$ such that if the initial state of the dynamical system (49) satisfies the relation

\[ \lim_{C\to + \infty} \sum_{k \ge \eta C}x_k^C(0)=0, \]

then there exists $\eta^*\in(0,1)$ such that

(62) \begin{equation} \lim_{C{\to} + \infty} \sup_{t \ge 0} \left(\sum_{k \ge \eta^* C} x_k^C(t)\right)=0. \end{equation}

Proof. We fix $\delta_0>0$ such that $\eta_0 = \nu(1 + a\delta_0)\in(\eta,1)$ , and we choose $\eta_1$ , $\eta_2\in(\eta_0,1)$ with $\eta_1 < \eta_2$ . For ${\varepsilon} < \delta_0/2$ , we take $C_0$ such that

\[ \sum_{k\geq \eta_0C_0} x_k^C(0)\leq {\varepsilon} \text{ and } \left(\frac{\eta_0}{\eta_1}\right)^{\lfloor (\eta_2-\eta_1)C_0\rfloor}\leq {\varepsilon}. \]

Let (Q(t)) be an $M/M/1$ queue with arrival rate $\eta_0$ and service rate $\eta_1$ with $Q(0) = 0$ , and let $(\overline{L}_C(t))$ be the processes defined by Equation (50). A simple coupling of (Q(t)) with a stationary version of the process (Q(t)) gives the stochastic monotonicity property of the process (Q(t)). (See the proof of Proposition 5.8 of Robert [Reference Robert29], for example.) We thus obtain that, for $C \ge C_0$ ,

\[{\mathbb P}\left(Q(t) \geq (\eta_2-\eta_1)C\right)\\ \leq {\mathbb P}(Q(\infty)\geq (\eta_2-\eta_1)C) \leq \left(\frac{\eta_0}{\eta_1}\right)^{\lfloor (\eta_2-\eta_1)C\rfloor}\leq {\varepsilon},\]

where $Q(\infty)$ is a geometrically distributed random variable with parameter $\eta_0/\eta_1$ . If $(Q_1(t))\stackrel{{\text {def}.}} = \eta_1 C + Q(t)$ , then, for all $t \ge 0$ ,

\[{\mathbb P}(Q_1(t)\geq \eta_2 C)\leq {\varepsilon}\leq \frac{\delta_0}{2}.\]

Note that $(Q_1(C t))$ is a birth-and-death Markov process with birth rate $\eta_0 C$ and death rate $\eta_1 C$ .

We can construct a coupling of the processes $(\overline{L}_C(t))$ and (Q(t)) such that, for all $t \ge 0$ , the relation

\[\overline{L}_C(t)\le Q_1(Ct)\]

holds conditionally on the event $\{\overline{L}_C(0) < \eta_0 C\}$ . This is a simple consequence of the fact that, as long as the relation $x^C_C(t) = {\mathbb P}(\overline{L}_C(t) = C)\le \delta_0$ holds, the input rate of $(\overline{L}_C(t))$ is smaller than $\eta_0 C$ , and when $\overline{L}_C(t)>\eta_1 C$ , the departure rate is at least $\eta_1 C$ . We obtain that, for $C \ge C_0$ and all $t \ge 0$ ,

\begin{align*}\sum_{k \ge \eta_2C} x_k^C(t) = {\mathbb P}\left(\overline{L}_C(t) \ge \eta_2C\right)& \le {\mathbb P}\left(\overline{L}_C(t) \ge \eta_2 C|\overline{L}_C(0)\le \eta_0 C\right) + {\varepsilon}\\ &\leq {\mathbb P}\left(Q_1(Ct) \geq \eta_2 C\right) + {\varepsilon}\leq 2{\varepsilon}.\end{align*}

The proposition is proved.

3.4. Nonlinear $M/M/1$ queues

The equations (58) defining the asymptotic process $(\Gamma(t))=(\Gamma_k(t),k\in{\mathbb N})$ can be written in a vectorial form as follows: for any function $f\,:\,{\mathbb N}{\to}{\mathbb R}_+$ with finite support,

(63) \begin{equation}\frac{\mathop{}\mathopen{}\mathrm{d} }{\mathop{}\mathopen{}\mathrm{d} t} \left\langle{\Gamma(t),f}\right\rangle =\left\langle{\Gamma(t),\nabla^+(f)}\right\rangle + \nu h(\Gamma(t)(0))\left\langle{\Gamma(t),\nabla^-(f)}\right\rangle,\end{equation}

with

\[\nabla^+(f)(x)\stackrel{{\text {def}.}} = f(x + 1)-f(x), \quad \nabla^-(f)(x)\stackrel{{\text {def}.}} = \left(f(x-1)-f(x)\right){\textbf{1}_{\left\{{x>0}\right\}}}.\]

In (58), we had $h(x) = 1 + ax(1-x)$ with $a>1$ .

We assume that h is a continuously differentiable function from [0, 1] to $[1, + \infty)$ , and that the relation $\nu>1$ holds. If there is an equilibrium $\pi\in{\mathcal P}({\mathbb N})$ for the dynamical system (63), it is the equilibrium of the linear Markov process $(\Gamma_S(t))$ , where $S = \pi(0)$ :

(64) \begin{equation}\frac{\mathop{}\mathopen{}\mathrm{d} }{\mathop{}\mathopen{}\mathrm{d} t} \left\langle{\Gamma_S(t),f}\right\rangle =\left\langle{\Gamma_S(t),\nabla^+(f)}\right\rangle + \nu h(S)\left\langle{\Gamma_S(t),\nabla^-(f)}\right\rangle.\end{equation}

As noted before, this is a classical $M/M/1$ queue which is ergodic since the service rate $\nu h(S)$ is greater than 1, the arrival rate, by assumption. We have therefore a representation for the invariant distribution

(65) \begin{equation} \pi_{S}(n)\stackrel{{\text {def}.}} = \left(\frac{1}{\nu h(S)}\right)^n\left(1-\frac{1}{\nu h(S)}\right), \quad n\in{\mathbb N},\end{equation}

and, consequently, the fixed point equation for S,

(66) \begin{equation}S=1-\frac{1}{\nu h(S)},\end{equation}

which we have already seen (see Equation (54)). It is well known that for the standard $M/M/1$ process $(\Gamma_S(t))$ , for any initial condition $\Gamma_S(0)\in{\mathcal P}({\mathbb N})$ , the following inequality holds for all $t \ge 0$ :

(67) \begin{equation} \frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} t}\|\Gamma_S(t)-\pi_S\|_2^2\le -2\kappa_S \|\Gamma_S(t)-\pi_S\|^2, \end{equation}

where $\kappa_S$ , the spectral gap of the process, has the explicit representation

\[ \kappa_S =\left(\sqrt{\nu h(S)}-1\right)^2.\]

See Chen [Reference Chen8] and Liu and Ma [Reference Liu and Ma23], for example.

Theorem 3. If $h\,:\,[0,1]{\to}[1, + \infty)$ is a $C^1$ function and $S\in(0,1)$ is a solution of (66) such that

(68) \begin{equation}|{h}^{\prime}(S)| < \frac{1}{\nu S}\left(\frac{1}{\sqrt{1-S}}-1\right)^2, \end{equation}

then the probability distribution $\pi_S$ defined by (65) is an exponentially stable equilibrium point of the dynamical system defined by (63): there exist positive constants q and ${\varepsilon}$ such that if $\|\Gamma(0)-\pi_S\|_2\le{\varepsilon}$ , the relation

\[\|\Gamma(t)-\pi_S\|_2^2\le \|\Gamma(0)-\pi_S\|_2^2\cdot e^{-q t},\]

holds for all $t \ge 0$ .

Proof. Identity (63) gives

\[{\Gamma}^{\prime}(t)(x) = \left(\rule{0mm}{4mm}\Gamma(t)(x-1)-\Gamma(t)(x)\right){\textbf{1}_{\left\{{x>0}\right\}}} + \nu h(\Gamma(t)(0))\left(\rule{0mm}{4mm}\Gamma(t)(x + 1)-\Gamma(t)(x)\right);\]

in the same way as with Equation (64), $({\Gamma}_S(t))$ is defined as the solution of

\[{\Gamma}^{\prime}_S(t)(x) = \left(\rule{0mm}{4mm}\Gamma_S(t)(x-1)-\Gamma_S(t)(x)\right){\textbf{1}_{\left\{{x>0}\right\}}} + \nu h(S)\left(\rule{0mm}{4mm}\Gamma_S(t)(x + 1)-\Gamma_S(t)(x)\right),\]

with the same initial condition $\Gamma_S(0) = \Gamma(0)$ . We have

\[{\Gamma}^{\prime}(0)(x)={\Gamma}^{\prime}_S(0)(x) + \nu(h(\Gamma(0)(0))-h(S))\left(\rule{0mm}{4mm}\Gamma(0)(x + 1)-\Gamma(0)(x)\right).\]

As in Subsection 2.5, setting

\[F(t)\stackrel{{\text {def}.}} = \|\Gamma(t)-\pi_S\|_2^2,\]

we have

(69) \begin{align}{F}^{\prime}(0) = & 2\sum_{x\in{\mathbb N}} \frac{[\Gamma_S(0)-\pi_S](x){\Gamma}^{\prime}_S(0)(x)}{\pi_S(x)}\nonumber\\ & + 2\nu (h(\Gamma(0)(0))-h(S))\sum_{x\in{\mathbb N}} \frac{[\Gamma(0)-\pi_S](x)(\Gamma(0)(x + 1)-\Gamma(0)(x))}{\pi_S(x)}\nonumber\\ & \stackrel{{\text {def}.}} = {\mathcal I}_1+2\nu (h(\Gamma(0)(0))-h(S))\cdot {\mathcal I}_2.\end{align}

By (67) and the fact that $(\Gamma_S(t))$ and $(\Gamma(t))$ have the same initial condition, we get

(70) \begin{equation}\frac{{\mathcal I}_1}{2}\le -2\kappa_S F(0)=-2\left(\sqrt{\nu h(S)}-1\right)^2F(0).\end{equation}

Furthermore, we have

\begin{align*} {\mathcal I}_2&\stackrel{{\text {def}.}} = \sum_{x\in{\mathbb N}}\frac{[\Gamma(0)-\pi_S](x)(\Gamma(0)(x + 1)-\Gamma(0)(x))}{\pi_S(x)}\\ &=\sum_{x\in{\mathbb N}}\frac{[\Gamma(0)-\pi_S](x)(\pi_S(x + 1)-\pi_S(x))}{\pi_S(x)}\\&\quad + \sum_{x\in{\mathbb N}}\frac{[\Gamma(0)-\pi_S](x)[\Gamma(0)-\pi_S](x + 1)}{\pi_S(x)}-\|\Gamma(0)-\pi_S\|_2^2.\end{align*}

The Cauchy-Schwartz inequality and (65) give

\begin{align*} |{\mathcal I}_2|&\le \left(\sum_{x\in{\mathbb N}} \frac{(\pi_S(x + 1)-\pi_S(x))^2}{\pi_S(x)}\right)^{1/2}\|\Gamma(0)-\pi_S\|_2 \\ &\quad + \sup_{x\in{\mathbb N}}\left(\sqrt \frac{\pi_S(x + 1)}{\pi_S(x)} + 1\right)\|\Gamma(0)-\pi_S\|_2^2\\ &=\left(1{-}\frac{1}{\nu h(S)}\right)\|\Gamma(0){-}\pi_S\|_2+ \left(\frac{1}{\sqrt{\nu h(S)}}{+}1\right)\|\Gamma(0){-}\pi_S\|_2^2. \end{align*}

Combining this with Equation (66) and Inequalities (69) and (70), as well as the fact that $|\Gamma(0)(0)-S|\le\|\Gamma(0)-\pi_S\|_2$ , we obtain that if $F(0)>0$ ,

\[\frac{{F}^{\prime}(0)}{2F(0)}\leq -\kappa_S +\nu \frac{|h(\Gamma(0)(0))-h(S)|}{|\Gamma(0)(0)-S|}\left(S+ \left(\frac{1}{\sqrt{\nu h(S)}} + 1\right)\sqrt{F(0)}\right);\]

hence

\[\limsup_{{\varepsilon}\to 0} \sup_{F(0)\in(0,{\varepsilon}^2)} \frac{{F}^{\prime}(0)}{F(0)}\leq -q=2\left(-\left(\sqrt{\nu h(S)}-1\right)^2+\nu |h^{\prime}(S)|S\right)<0\]

by the assumption (68). To complete the proof, it remains to apply a similar bootstrap argument as in the proof of Theorem 1 to get the above inequality for the ratio ${F}^{\prime}(s)/F(s)$ . The theorem is proved.

We now apply this result to the asymptotic dynamical system of the DAR algorithm.

Proof. The fixed point equation (66) is

(71) \begin{equation}(1-S)(1 + aS (1-S))=\frac{1}{\nu},\end{equation}

and since $h^{\prime}(x) = a\nu(1-2x)$ , the condition (68) is equivalent to

(72) \begin{equation} aS |1-2S| < \left(1-\sqrt{1-S}\right)^2 (1 + aS(1-S)).\end{equation}

To conclude, note that $S = 1/2$ satisfies the condition.

Remarks.

1. 1. It is easy to check that, for $a>4$ , $u_a\in(\nu_a,1)$ , where $\nu_a$ is defined by (53). In this case there are two positive fixed points for the asymptotic dynamical system; the above corollary gives that one of them is locally stable. We have not been able to prove that the other one is unstable, as suggested by Gibbens et al. [Reference Gibbens, Hunt and Kelly16] (see the quotation on page 41).

2. 2. A little more work can give more precise conditions on $\nu$ for the stability of the fixed point. Let $x = \sqrt{1-S}$ . If $S\in[0,1/2]$ , Condition (72) amounts to

   \[P_1(x)\stackrel{{\text {def}.}} = ax^{5}-ax^{4}-3a{x}^3 -ax^{2} + (a-1) x + a + 1>0.\]
   Notice that $P_1(\sqrt{2}/2) = (2-\sqrt{2})(1 + a/4)/2$ and $P_1(1) = -2a$ . If $S\in[1/2,1]$ , the condition is
   \[P_2(x)\stackrel{{\text {def}.}} = ax^5-a{x}^{4} + a{x}^{3} + 3a{x}^{2}-(1 + a)x + 1-a>0,\]
   with $P_2(0) = 1-a < 0$ and $P_2(\sqrt{2}/2) = P_1(\sqrt{2}/2)>0$ . It is not difficult to check that $P_1$ (resp. $P_2$ ) is concave (resp. convex) on [0, 1]; hence there exists a unique root $z_{a,1}$ of $P_1$ in $(\sqrt{2}/2,1)$ (resp. $z_{a,2}$ of $P_2$ in $(0,\sqrt{2}/2)$ ). Hence the condition (68) is satisfied when $S\in(1-z_{a,1}^2,1-z_{a,2}^2)$ , and, by Equation (71), this holds if $\nu\in(Q(z_{a,2}),Q(z_{a,1}))$ , where $Q(z) = 1/[z^2(1 + az^2-az^4)]$ .

3. 3. For the precise case of Gibbens et al. [Reference Gibbens, Hunt and Kelly16], $h(x) = (1 + 2x(1-x))$ , this shows that when $\nu\in(1.2068,1.5978)$ , the unique fixed point is a locally stable equilibrium.

A toy example with an arbitrary number n of stable equilibrium points.

We fix $(u_k)$ , a set of n distinct points in (0,1). Let f be a $C^1$ function such that, for each $k\in\{1, \ldots, n\}$ , the relation $f(u) = 1 + \ln(1-u) - \ln(1-u_k)$ holds in a small neighborhood of $u_k$ . Note that since $f(u_k) = 1>1-u_k$ for $1\le k\le n$ , we can choose f in such a way that $f(u)>1-u$ holds for all $u\in(0,1)$ . If we define $h(u) = f(u)/(1-u)$ for $u\in(0,1)$ , then h maps (0,1) to $(1, + \infty)$ . Taking $\nu = 1$ , we have that each $u_k$ , $1\le k\le n$ , is clearly a fixed point, and the condition (68) is satisfied since $h'(u_k) = 0$ . It is therefore locally stable.