# On multiclass spatial birth-and-death processes with wireless-type
interactions

Pierre Popineau1, François Baccelli1,2
Affiliation: 1 INRIA-ENS, Paris, France 2 Telecom Paris, France

###### Abstract

This paper studies a multiclass spatial birth-and-death (SBD) processes on a
compact region of the Euclidean plane modeling wireless interactions. In this
model, users arrive at a constant rate and leave at a rate function of the
interference created by other users in the network. The novelty of this work
lies in the addition of service differentiation, inspired by bandwidth
partitioning present in 5G networks: users are allocated a fixed number of
frequency bands and only interfere with transmissions on these bands.

The first result of the paper is the determination of the critical user
arrival rate below which the system is stochastically stable, and above which
it is unstable. The analysis requires symmetry assumptions which are defined
in the paper. The proof for this result uses stochastic monotonicity and fluid
limit models. The monotonicity allows one to bound the dynamics from above and
below by two adequate discrete-state Markov jump processes, for which we
obtain stability and instability results using fluid limits. This leads to a
closed form expression for the critical arrival rate. The second contribution
consists in two heuristics to estimate the steady-state densities of all
classes of users in the network: the first one relies on a Poisson
approximation of the steady-state processes. The second one uses a cavity
approximation leveraging second-order moment measures, which leads to more
accurate estimates of the steady-state user densities. The Poisson heuristic
also gives a good estimate for the critical arrival rate.

##  1 Introduction

A key feature of wireless networks is the presence of spatial interactions
between users. More precisely, the service rate of a wireless link at any
given time is a function of the interference created by other users
transmitting at the same time. This follows from Shannon’s formula ([1]) for
the transmission rate of a wireless link when treating interference as noise.

In [2], the authors introduced a framework to study such interactions using
the theory of spatial birth-and-death processes. The setting of [2] is that of
a device to device wireless network, where users arrive at random locations in
time and space and where their transmission rate is determined by Shannon’s
formula. They derived the stability region and studied the stationary regime
of this class of dynamics.

An important novelty introduced in 5G networks is the partitioning of the
available bandwidth in frequency bands of equal width and the simultaneous use
of multiple frequency bands ([3, 4]) by transmissions depending on their
nature: transmissions requiring lower data rates (phone calls or text
messages) are allocated less bands than those requiring higher data rates
(video streaming for instance). This more flexible and adaptive frequency
bandwidth allocation is meant to increase the network capacity.

Service differentiation based on users needs has been extensively studied,
e.g., in queueing network theory. For instance, BCMP type queueing networks
[5] extend the Jackson framework to the multi-class setting. For the class of
device-to-device wireless networks introduced in [2], the multiclass paradigm
remains an open problem, which is of central importance in the bandwidth part
setting. As mentioned above, in this setting, transmissions can use several
frequency bands simultaneously. Those requiring a higher data rate use a
larger number of bands. They hence transmit with a stronger signal, which
potentially increases their transmission rate. The main novelty and difficulty
is that they then also generate a higher level of interference to other users,
hence slowing down other transmissions. Note that several new phenomena arise
in this setting. For instance, users which are allocated frequency bands that
do not overlap do not directly interact. In contrast, users which are
allocated the same large set of overlapping bands and are nearby interfere a
lot and slowdown each other quite significantly, which in turn slows down all
other nearby transmissions for a longer time.

In view of all these complex interactions, several important questions arise.
The first one is that of the stability region of this class of dynamics.
Stochastic stability has always been of central importance in communication
networks, because it ensures that e.g. user latency converges to a stationary
distribution. A second important question is that of the analytical
characterization of this stationary distribution. For instance, in queueing
network theory, the stability region of Jackson networks is well known as well
as the stationary regime [6]. Another important in this service
differentiation setting is that of the best allocation of resources to users
in the network. An instance of allocation policy is that which attributes
radio resources proportionally to the needs of users (see [7, 8]). Our
approach is oblivious of this type of assumption. In the present paper,
arriving users in the network have their needs categorized in terms of the
mean size of the file that they have to transmit. All monotonic policies
(where larger file sizes lead to larger sets of bands) can be represented and
analyzed.

The first goal of the present paper is to extend the spatial birth-and-death
wireless network framework of [2] to the multiclass model described above. As
in [2], we use a birth-and-death process of wireless dipoles (see [9, 10]),
where dipoles arrive to a compact according to a Poisson rain of constant
intensity, representing the arrival rate of users in the network, and leave
with a departure rate proportional to the their Shannon rate in the spatial
configuration. This problem is quite complex and has been an open question for
several years. As we will see, there is for instance no way to guess the
stability condition of the multiclass case from that of the monoclass setting
derived in [2].

The present paper, which is a first step in the solution of the general model,
makes the following simplifying assumptions: (i) the traffic satisfies certain
natural symmetry assumptions (ii) the compact is assumed to be a torus, and
(iii) the wireless network operates in the low SINR case. Here is the
rationale for these assumptions. (i) Assuming a compact region of space allows
one to use fluid limit techniques to asses stochastic stability. These
techniques cannot be used when the phase space is the whole Euclidean plane
(see [11, 12] for instance). Assuming a square torus allows one to emulate the
Euclidean plane and is the simplest possible compact state space model since
it is invariant by translations modulo the square side. (ii) The symmetry
assumptions require that bands are in a sense exchangeable and that users with
the same needs use the same number of bands and that bands are chosen
uniformly at random. These assumptions are rather natural and are at the same
time essential in the derivations for the stability region. (iii) the use of
the linearized version of the Shannon formula, which is justified by the low
SINR regime assumption, is probably the least essential of the three
assumptions and we explain in the paper how to relax this.

The first part of the paper is focused on the stability region of such a
multiclass network. This part relies on the use of the definition of certain
interference queueing networks (see [12]) to obtain bounds on the stability
region. Due to the nature of these dynamics, fluid limits ([13]) provide a
natural framework to obtain conditions for stability ([14]).

The second part uses stochastic geometry tools (see [9]) as well as Palm
calculus (see [15]) in order to obtain heuristics for user densities in the
stationary regime. A first heuristic relies on a mean-field approximation for
the interacting point processes. Such approximations are known to provide a
powerful framework to obtain both qualitative and quantitative results. An
instance of Poisson approximation is that arising in replica mean-field
systems (see [16]), where the authors prove that as the number of replicas
goes to infinity, point processes behave like independently distributed
Poisson processes. In another setup, [17] and [18] provide a quantification of
the error of the mean-field approximation.

### Organization of the paper

As mentioned earlier, we model wireless interactions using a space-time
Poisson dipole model (see [2] and [10]). This model and the wireless
interactions we consider are defined in Section 2. _Symmetric_ networks are
introduced in Section 2.2); in such networks, users allocated the same number
of communication channels have similar properties in the network, e.g., have
the same arrival rate and the same file size distribution. This hypothesis
allows to use combinatorics in order to obtain a closed form for the critical
arrival rate in the system (through Lemma 8.5).

The irreducibility and monotonicity properties which allow one to assess the
stability region for the dynamics are gathered in Section 3. Whereas the
initial dynamic is measure-valued, monotonicity allows one to reduce the study
to that of two interference queueing networks (see [12]) through an adequate
tessellation of space (see Section 4.1). These two queueing networks provide
an inner and an outer bound for the stability region, which yield the main
stability result of this paper, stated in Theorem 4.1.

To obtain the upper bound for the critical arrival rate, we rely on a fluid
limit model, presented in Section 4.2 and on results linking fluid limits and
positive recurrence (see [19] and [14]). To obtain instability, we use
monotonicity to tie the instability of the network with the instability of an
appropriate M/M/1 queue.

In the second part of the paper, we give two heuristics allowing one to
quantitatively estimate the first moment of the stationary regime. These
heuristics rely on stochastic geometry and rate conservation arguments.
Section 5.1 discusses a first-order heuristic for the stationary user density.
In Section 5.3, we improve the quality of this heuristic by using a cavity
approximation. We also compare all these heuristics with the results obtained
by discrete event simulations in Section 6.

##  2 Mathematical framework

###  2.1 Network setup

We consider an infrastructureless wireless network, where arrivals of
transmitters follow a Poisson rain of intensity $\lambda>0$ on a compact
subset $\mathcal{D}$ of $\mathbb{R}^{2}$. Throughout the paper, we assume that
$\mathcal{D}$ is a square torus.

To each arriving transmitter, we associate a receiver located at a fixed
distance $r\geq 0$ in a random direction, with $r$ small compared to the side
of the square. This is a Poisson dipole model as introduced in [9]. We
describe the configuration of users present in the network at a given time as
receiver-transmitter pairs, denoted by
$\Phi_{t}=\\{(x_{1},y_{1}),\dots(x_{N_{t}},y_{N_{t}})\\}$, where $N_{t}$ is
the number of pairs present in the network at time $t$, $(x_{i})_{1\leq i\leq
N_{t}}$ denotes the location of receivers and $(y_{i})_{1\leq i\leq N_{t}}$
that of transmitters. Let $\Phi^{T}_{t}$ be the point process describing the
location of transmitters in the network and $\Phi^{R}_{t}$ be that of the
locations of receivers at a given time $t$. We will refer to receiver-
transmitter pairs or dipoles as _users_ in the rest of the paper.

$K$ orthogonal transmission bands/channels of equal width are available to
users arriving in the network. Channels $i\neq j$ do not interfere with one
another. Let $\mathcal{P}(K)$ denote the set of non-empty subsets of
$[1,2,\dots,K]$. An arriving user selects a set of channels
$C\in\mathcal{P}(K)$ on which to transmit according to a given distribution
$\\{p_{C}\\}_{C\in\mathcal{P}(K)}$. This selection is made independently of
the state of the network. In the rest of the paper, we will call $C$ the
_class_ of a user.

A receiver-transmitter pair $(x,y)\in\mathcal{D}$ of class $C$ arrives with a
file of random size $L_{x,y}$ attached. File sizes for users in the same class
$C$ are i.i.d. random variables distributed with an exponential law of mean
$L_{C}$. Once the file has been transmitted to the receiver, the pair or
dipole $(x,y)$ leaves the network.

###  2.2 Symmetric system

We say that the network is _symmetric_ , if all users of classes with the same
cardinality have the same statistical properties in the network. In
particular, whenever $\lvert C\rvert=\lvert D\rvert$, we have $p_{C}=p_{D}$
and $L_{C}=L_{D}$.

An instance of this symmetry property is that of a system where user needs are
split in $K$ categories, depending on their needs (e.g. phone calls, text
messages, web browsing or video streaming). Each arriving user has a
probability $p_{j}$ of having a need of category $j$. Upon arrival, users with
needs of category $j$ are given a set of $j$ bands to transmit, sampled
uniformly at random among the $\binom{K}{j}$ possibilities. Thus, the
probability that a user transmits on a given set $C$ of $j$ bands is equal to
$p_{C}=\frac{p_{j}}{\binom{K}{j}}$. The quantity of information this arriving
user has to transmit to the network is thus sampled from an exponential
distributions with parameter $L_{j}$ with $1\leq j\leq K$. Setting
$L_{C}=L_{j}$ and $p_{C}=\frac{p_{j}}{\binom{K}{j}}$ gives an example of a
symmetric system. The symmetry in the system will be important to simplify the
combinatorics we will encounter in Section 4, when proving Theorem 4.1.

The question of non-symmetric networks, where the bandwidth allocation for
users allows for a differentiated use of communication channels are beyond the
scope of this paper, and will be the subject of future work.

###  2.3 Mathematical definitions and notation

In this paper, vectors will be denoted as bold-faced letters and their
coordinates in regular script (for instance, $\mathbf{x}=(x_{i})_{0\leq i\leq
d-1}\in\mathbb{R}^{d}$), and $\leq$ will denote the coordinate-wise partial
ordering when used for vectors. Fluid limits associated with a stochastic
process will be denoted using the same notation as the stochastic process they
are derived from, in lower script, e.g. $z$ is the fluid limit associated with
the process $\mathbf{Z}$.

We denote by $\mathcal{M}(\mathcal{D})$ the set of counting measures on
$\mathcal{D}$, all defined on a probability space
$(\Omega,\mathcal{F},\mathbb{P})$.

For $N,K\in\mathbb{N}$, we denote by $D([0,\infty),\mathbb{R}^{N\times
2^{K}-1})$ the set of càdlàg functions from $\mathbb{R}_{+}$ to
$\mathbb{R}^{N\times 2^{K}-1}$ and $\mathcal{B}$ is its canonical Borel
$\sigma$-algebra. We consider stochastic processes as measurable maps from
$(\Omega,\mathcal{F})$ to $(D([0,\infty),\mathbb{R}^{N\times
2^{K}-1}),\mathcal{B})$.

Let $d_{\infty}^{0}$ be the infinite norm in $\mathbb{R}^{N\times 2^{K}-1}$,
$\lvert x\rvert$ be the $L^{1}$ norm and $\lVert x\rVert$ denote the $L^{2}$
norm of $x$. For a given compact set $\mathcal{D}\subset\mathbb{R}^{2}$, we
define $\langle\cdot_{\mathcal{D}}\rangle$ as follows:

$\langle f_{\mathcal{D}}\rangle=\begin{cases}\int_{\mathcal{D}}f(\lVert
x\rVert)\mathrm{d}x,\quad\forall f:\mathbb{R}\rightarrow\mathbb{R}\\\
\int_{\mathcal{D}}f(x,0)\mathrm{d}x,\quad\forall
f:\mathbb{R}^{2}\times\mathbb{R}^{2}\rightarrow\mathbb{R},\end{cases}$ (1)

whenever the integrals are defined. We use here the same notation for two
different notions; this choice is motivated by results presented in Section 4.

For a sequence $(x^{n})_{n\geq 0}$ and $x$ in
$D([0,\infty),\mathbb{R}^{N\times 2^{K}-1})$, we write $x^{n}\rightarrow x$
if, for all $T>0$ :

$\lim_{n\rightarrow\infty}\sup_{0\leq t\leq T}\lVert x^{n}(t)-x(t)\rVert=0.$

We denote by $\mathcal{E}(\mathbb{R}^{d})$ the set of probabilities on
$\mathbb{R}^{d}$ (with $d\geq 1$), and we define the partial ordering of
probability measures $\leq_{i}$. We define:

$I=\\{f\in\mathcal{E}(\mathbb{R}^{d}),f\text{ is coordinate-wise non-
decreasing}\\}.$

For each $F,G\in I$, we write $F\leq_{i}G$ if and only if

$\int_{\mathbb{R}^{d}}f(x)F(\mathrm{d}x)\leq\int_{\mathbb{R}^{d}}f(x)G(\mathrm{d}x),$

for all $f\in I$. $\leq_{i}$ is a partial order on the set of random vectors
of $\mathbb{R}^{d}$. Let $\Psi$ and $\Psi^{\prime}$ be two point processes on
$\mathcal{D}$. We write $\Psi\leq\Psi^{\prime}$ if we have
$\Psi(\mathcal{D})\leq_{i}\Psi^{\prime}(\mathcal{D})$. In the rest of the
paper, we will use this partial order when comparing stochastic processes.

Finally, for two random variables $X$ and $Y$, we write $X\overset{d}{\sim}Y$
if $X$ and $Y$ have the same probability distribution. Table 1 summarizes the
notation used throughout the paper.

Notation | Description
---|---
$\Phi^{R}_{t},\Phi^{T}_{t}$ | Receiver (resp. transmitter) locations at time $t$
$\Phi_{C,t}$ | Point process of users of class $C$ at time $t$
$\mathcal{P}(K)$ | Set of subsets of $[1,\dots K]$ with the exception of $\emptyset$
$\ell$ | Path-loss function
$r$ | Receiver-transmitter distance
$\lambda$ | Intensity of the arrival process
$(p_{C})_{C\in\mathcal{P}(K)}$ | Arrival distribution of users of each type
$L_{C}$ | Average file size for users of class $C$
$R$ | Rate-of-transmission function
$\mathcal{N}_{0}$ | Thermal noise density in the network
$\langle\ell_{\mathcal{D}}\rangle$ | $\int_{\mathbf{x}\in\mathcal{D}}\ell(\lVert\mathbf{x}\rVert)\mathrm{d}\mathbf{x}$
$\Phi_{0,C}$ | Stationary point process of users of type $C$
$\mu_{C}$ | Spatial intensity of point process $\Phi_{0,C}$
$\mathbb{E}^{0}_{\Phi_{0,C}}$ | Palm expectation with respect to $\Phi_{0,C}$
Table 1: Table of notations

###  2.4 Wireless interactions and service times

Let $\ell$ be a non-negative, bounded and non-increasing function with
$\langle\ell_{\mathcal{D}}\rangle<\infty$ be a path-loss function. Without
loss of generality, we assume that $\ell(0)=1$.

The interference experienced by a receiver located at $x\in\mathcal{D}$ of
class $C_{x}$, whose transmitter is located at $y\in\Phi_{t}^{T}$, is equal
to:

$I(x,\Phi_{t})=\sum_{z\in\Phi^{T}_{t}\backslash\\{y\\}}\lvert C_{x}\cap
C_{z}\rvert\ell(\rVert x-z\lVert),$ (2)

where $C_{z}$ is the class of the transmitter located at $z$. Note that a
receiver does not interfere with its own transmitter. As mentioned earlier,
the interference is a shot-noise of the transmitter point process, depending
on the number of overlapping channels used by transmitters and receivers. We
thus define the rate-of-transmission function for this receiver as follows:

$R(x,\Phi_{t})=\frac{A\lvert
C_{x}\rvert\ell(r)}{\mathcal{N}_{0}+I(x,\Phi_{t})},$ (3)

where $\mathcal{N}_{0}>0$ is the thermal noise power and $A$ is a
multiplicative constant.

(3) can be seen as a linearization of the Shannon-Hartley formula in the low-
SINR approximation: $\lvert C_{x}\rvert\ell(r)$ is the signal power received
by the receiver located at $x\in\mathcal{D}$ from its transmitter, and
$\mathcal{N}_{0}+I(x,\Phi_{t})$ is the noise and interference power seen by
this receiver.

Under these considerations, in the low-SINR case, i.e., if $\lvert
C_{x}\rvert\ell(r)\ll\mathcal{N}_{0}+I(x,\Phi_{t})$, the Shannon-Hartley (see
[1]) formula gives a channel capacity equal to:

$\mathcal{C}(x)=\log_{2}\left(1+\frac{\lvert
C_{x}\rvert\ell(r)}{\mathcal{N}_{0}+I(x,\Phi_{t})}\right)\approx\frac{1}{\ln(2)}\frac{\lvert
C_{x}\rvert\ell(r)}{\mathcal{N}_{0}+I(x,\Phi_{t})},$

which is the expression from (3) with $A=\frac{1}{\ln(2)}$. In the rest of the
paper, we assume that $A=1$.

The service time for a pair $(x,y)$ with a file of size $L_{x,y}$, denoted by
$t_{s,x}$ is defined by:

$t_{s,x}=\min\left\\{w>t_{x}:\int_{t_{x}}^{w}R(x,\Phi_{u})\mathrm{d}u>L_{x,y}\right\\}.$
(4)

The instantaneous departure rate of a receiver-transmitter pair of class $C$
with receiver located at $x\in\mathcal{D}$, at time $t$, is equal to:

$d(x,t)=\frac{1}{L_{C}}R(x,\Phi_{t}).$ (5)

##  3 Preliminary properties and definitions

We establish a stochastic recurrence for the process $\Phi_{t}$, which will
help us asserting a central property for the network, is _stochastic
monotonicity_.

A random variable $X$ with values in $\mathbb{R}^{d}$ _dominates_ $Y$ if the
cumulative distribution functions of $X$ and $Y$, denoted as $F_{X}$ and
$F_{Y}$, are such that $F_{Y}\leq_{i}F_{X}$. A consequence of stochastic
domination is, for all $\mathbf{x}\in\mathbb{R}^{d}$ and for all probabilities
$\mathbb{P}$ defined on $\mathbb{R}^{d}$, we have
$\mathbb{P}[X>\mathbf{x}]\geq\mathbb{P}[Y>\mathbf{x}]$.

We can define a stochastic recurrence relation for each of the subprocesses
$\Phi_{C,t}$: denote as $(T_{C,n})_{n\leq 0}$ the times at which events
(either an arrival or a departure) happen in $\Phi_{C,t}$. The relation
between $\Phi_{C,T_{C,n}}$ and $\Phi_{C,T_{C,n+1}}$ is given by:

$\begin{cases}\Phi_{C,T_{C,n+1}}&=\Phi_{C,T_{C,n}}+\mathcal{B}_{C,n}-\mathcal{D}_{C,n},\\\
T_{C,n+1}&=T_{C,n}+\delta_{C,n},\end{cases}$ (6)

where $\mathcal{B}_{C,n}$ and $\mathcal{D}_{C,n}$ are
$\mathcal{M}(\mathcal{D})$-valued random variables that are defined as
follows: let $b_{C}$ be a real-valued random variable distributed according to
an exponential distribution with rate $\lambda p_{C}$ and $\mathbf{d}_{C,n}$
be a $\mathbb{R}^{\Phi_{C,n}(\mathcal{D})}$-valued random variable where each
coordinate is distributed according an independent exponential distribution
with mean $\frac{1}{L_{C}}R(x_{i},\Phi_{T_{C,n}})$, for each
$x_{i}\in\Phi_{C,T_{C,n}}$. Then, the values of $\mathcal{B}_{C,n}$ and
$\mathcal{D}_{C,n}$ are given as follows :

  * -

$\mathcal{B}_{C,n}$ is equal to the null measure if
$b_{C}\geq\min\mathbf{d}_{C,n}$, or equal to $\delta_{x}$, where $x$ is
sampled uniformly in $\mathcal{D}$ otherwise;

  * -

$\mathcal{D}_{C,n}$ is equal to the null measure if
$b_{C}\leq\min\mathbf{d}_{C,n}$ or if $\Phi_{C,n}(\mathcal{D})=0$, or equal to
$\delta_{x_{i}}$ where $i$ is the smallest coordinate of $\mathbf{d}_{C,n}$
otherwise.

Finally, the time between events is equal to
$\delta_{C,n}=\min\\{b_{C},\min\mathbf{d}_{C,n}\\}$. More results about
partial ordering of stochastic recurrences can be found in Chapter 4 of [15].

###  3.1 Stochastic monotonicity

Stochastic monotonicity will help us obtaining bounds for the dynamics of the
network by defining queueing systems that dominate or are dominated by our
dynamics. We need the following theorem:

######  Theorem 3.1.

Let $\Delta:(\Phi_{0},\lambda,L,R)\mapsto\Phi_{t}$ define the realization of
the dynamics, where $\Phi_{0}$ is the initial condition of the network,
$\lambda$ is the arrival rate, $\mathbf{L}$ is the vector of average file
sizes and $R$ is the rate-of-transmission function. Let
$\Phi=\Delta(\Phi_{0},\lambda,d)$ and
$\Phi^{\prime}=\Delta(\Phi^{\prime}_{0},\lambda^{\prime},d^{\prime})$. The
following conditions are sufficient for $\Phi^{\prime}$ to dominate $\Phi$
(with all the other parameters taken equal) :

  * $i)$

$\lambda\leq\lambda^{\prime}$;

  * $ii)$

$\mathbf{L}\leq\mathbf{L}^{\prime}$;

  * $iii)$

for all point processes $\Psi\leq\Psi^{\prime}$ on $\mathcal{D}$ and
$x\in\mathcal{D}$, $R^{\prime}(x,\Psi^{\prime})\leq R(x,\Psi)$;

  * $iv)$

$\Phi_{0}\subseteq\Phi^{\prime}_{0}$.

It is to note that condition $iii)$ is met when we have two path-loss
functions $\ell$ and $\ell^{\prime}$ such that for all $r\geq 0$,
$\ell(r)\leq\ell^{\prime}(r)$.

###### Proof.

To obtain domination, we use a coupling argument between two instances of the
dynamics to obtain the domination relation, similar as the ones used in
Appendix B of [12].

Let us take $0<\lambda<\lambda^{\prime}$ with the same fixed initial condition
$\Phi_{0}$, a vector $\mathbf{L}$ and a rate-of-transmission function $R$, and
let $\Phi_{t}=\Delta(\Phi_{0},\lambda,d)$ and
$\Phi^{\prime}_{t}=\Delta(\Phi_{0},\lambda^{\prime},d)$. The quantities
related to $\Phi^{\prime}$ will be denoted with a prime. We want to prove that
for all times $0\leq t$, $\Phi_{t}\subseteq\Phi^{\prime}_{t}$.

We couple both the arrival and the departure processes in the network as
follows: the arrival process of $\Phi$ is a Poisson rain $\mathcal{A}$ with
parameter $\lambda$, and the arrival process for $\Phi^{\prime}$ is
$\mathcal{A}\cup\mathcal{A}^{\prime}$, with $\mathcal{A}^{\prime}$ being a
Poisson rain with intensity $\lambda^{\prime}-\lambda>0$ independent from
$\mathcal{A}$, so that common arrivals in $\Phi_{t}$ and $\Phi^{\prime}_{t}$
happen at the same locations and times. Using this coupling, we will show
that, at all times $t$, $\Phi_{t}\subseteq\Phi^{\prime}_{t}$.

At time $t=0$, we have trivially $\Phi_{0}\subseteq\Phi_{0}$. At time $t>0$,
we know that both $\Phi_{t}(\mathcal{D})$ and $\Phi^{\prime}_{t}(\mathcal{D})$
are almost surely finite. Let us write
$\Phi_{t}=\sum_{C\in\mathcal{P}(K)}\sum_{i=1}^{\Phi_{C,t}(\mathcal{D})}\delta_{x_{C,i}}$
and
$\Phi^{\prime}_{t}=\sum_{C}\sum_{i=1}^{\Phi^{\prime}_{C,t}(\mathcal{D})}\delta_{x^{\prime}_{C,i}}$.

Assume that up to time $t$, $\Phi_{t}\subseteq\Phi^{\prime}_{t}$, the next
event happening at a time $\hat{t}\geq t$ can be one of the following nature:

  * -

The arrival of a user of class $C$ in $\Phi_{t}$;

  * -

The arrival of a user of class $C$ in $\Phi^{\prime}_{t}$;

  * -

The departure of a user of class $C$ in $\Phi_{t}$;

  * -

The departure of a user of class $C$ in $\Phi^{\prime}_{t}$.

Because arrivals are coupled, an arrival in either $\Phi_{t}$ or
$\Phi^{\prime}_{t}$ maintains the inclusion. The same holds if an element of
$\Phi_{t}$ leaves. The last case to look at is the departure of an element $x$
of $\Phi^{\prime}_{t}$.

Let $\hat{t}^{-}$ be such that $\hat{t}^{-}=\hat{t}$ but the departure of $x$
has not happened yet. From our assumptions, we know that
$\Phi_{D,\hat{t}^{-}}\subseteq\Phi^{\prime}_{D,\hat{t}^{-}}$ for all classes
$D\in\mathcal{P}(K)$. Assume the next departure is for a user at $x$ of class
$C\in\mathcal{P}(K)$:

$\displaystyle R(x,\Phi_{\hat{t}^{-}})$ $\displaystyle=\frac{\lvert
C_{x}\rvert\ell(r)}{\mathcal{N}_{0}+\sum_{y\in\Phi_{\hat{t}^{-}}}\lvert
C_{x}\cap C_{y}\rvert\ell(\lVert x-y\rVert)}$
$\displaystyle\overset{(a)}{\geq}\frac{\lvert
C_{x}\rvert\ell(r)}{\mathcal{N}_{0}+\sum_{y\in\Phi^{\prime}_{\hat{t}^{-}}}\lvert
C_{x}\cap C_{y}\ell(\lVert x-y\rVert)}$
$\displaystyle=R(x,\Phi^{\prime}_{\hat{t}^{-}}),$

where $(a)$ uses the fact that
$\lvert\Phi_{t}(\mathcal{D})\rvert\leq\lvert\Phi^{\prime}_{t}(\mathcal{D})\rvert$.
In other words, the departure rates in $\Phi_{\hat{t}^{-}}$ are larger than in
$\Phi^{\prime}_{\hat{t}^{-}}$.

We take the _Poisson imbedding_ of the departure processes (see [20]): let
$\mathcal{D}_{C,t}$ and $\mathcal{D}^{\prime}_{C,t}$ be the point processes of
users of class $C$ that left each system up to time $t$. These processes have
respective stochastic intensities $\frac{1}{L_{C}}R(x,\Phi_{t})$ and
$\frac{1}{L_{C}}R(x,\Phi^{\prime}_{t})$. Using Lemma 3 from [20], we can imbed
them on the same Poisson point process $\mathcal{N}$ of intensity 1 on
$\mathbb{R}^{2}$.

Using this Poisson imbedding, any point $x$ leaving
$\Phi^{\prime}_{C,\hat{t}^{-}}$ has already left $\Phi_{C,\hat{t}^{-}}$, which
proves that the inclusion is maintained if the next event to come is a
departure in $\Phi^{\prime}_{t}$ and this concludes the proof of Condition
$i)$.

To obtain the other conditions, we use the same argument and we compare the
rate-of-transmission functions in each case to get the required inclusion. ∎

Theorem 3.1 is central to our study, because it will allow us to bound from
above and below the dynamics of the network we are studying in order to obtain
bounds for the limits of the stability region in the network.

###  3.2 Irreducibility

The second important property about the process $\Phi_{t}$ is the following:

######  Theorem 3.2.

$\Phi_{t}$ is a $\phi$-irreducible Markov jump process on
$\mathcal{M}(\mathcal{D})$.

It is to note that the similarity in notation between $\phi$-irreducibility
and the point process $\Phi_{t}$ is coincidental, the two notions are
independent.

###### Proof.

Our goal is to define a measure on the set of counting measures
$\mathcal{M}(\mathcal{D})$. We start by letting $\phi$ be the set function
defined as follows:

  * –

$\phi(\\{0_{\mathcal{M}(\mathcal{D})}\\})=1$, where
$0_{\mathcal{M}(\mathcal{D})}$ is the measure associated with the null
counting measure,

  * –

For $k\in\mathbb{N}$ and $B_{1},B_{2},\dots,B_{N}$ disjoint Borel subsets of
$\mathcal{D}$, we define the event
$A_{n}=\\{\Phi\in\mathcal{M}(\mathcal{D}):\Phi(\mathcal{D})=n,\Phi=(y_{1},y_{2},\dots,y_{n})\in
B_{1}\times\dots\times B_{n},\Phi(B_{1})>0,\dots\Phi(B_{n})>0\\}$. We then
set:

$\phi(A_{n})=\frac{1}{2^{n}}H(B_{1})H(B_{2})\dots H(B_{n}),$

where $H$ is a Haar measure on the square torus $\mathcal{D}$.

Here, $\phi$ is a set function on the _semiring_ of sets
$(A_{n})_{n\in\mathbb{N},(B_{i})\in\mathcal{F}^{n}}$. Using Theorem 11.3 from
[21], we can extend $\phi$ in a unique way to a measure on the $\sigma$-field
generated by the set of events
$\\{A_{n}:n\in\mathbb{N},(B_{i})\in\mathcal{F}^{n}\\}$, which is equal to
$\mathcal{F}$. For ease of notation, we will denote this measure by $\phi$ as
well. This way, $\phi$ is a measure defined on the set of counting measures
$\mathcal{M}(\mathcal{D})$.

To obtain $\phi$-irreducibility, we will proceed in two steps: in the first
step, we prove that, from each state $\Phi$ with $N$ points, we can reach the
empty state with positive probability, and in the second step, we prove that,
from the empty state, we can reach any state $A$ with $\Phi(A)>0$ with
positive probability in a finite number of steps, which allows us to conclude.

Let us assume that at a given time $t$, $N$ receiver-transmitter pairs are
present in the system, with locations $(x_{i},y_{i})_{1\leq i\leq N}$,
$x_{i}\in\Phi_{t}^{R}$ and $y_{i}\in\Phi_{t}^{T}$, and let $(C_{i})_{1\leq
i\leq N}$ be the classes of these users. The probability that we have the
departure of $(x_{1},y_{1}),\dots,(x_{n},y_{n})$ in the next $N$ steps is
equal to:

$\mathbb{P}[N\text{ departures in a
row}\left\lvert\Phi_{t}\right.]=\prod_{i=1}^{N}\frac{d_{C_{i}}}{d_{C_{i}}+b_{C_{i}}},$

where $b_{C_{i}}=\lambda p_{C_{i}}\lvert\mathcal{D}\rvert$ and
$d_{C_{i}}=\frac{1}{L_{C_{i}}}R(x_{i},\Phi_{t_{i}})$, with $\Phi_{t_{i}}$
being the network configuration after the departure of the user
$(x_{i},y_{i})$. This probability is non-null, i.e., we can reach the null
measure with nonzero probability from any other state.

Similarly, assume we start from the empty state, let $n$ be an integer,
$B_{1},\dots,B_{n}$ be $n$ disjoint Borel subsets of $\mathcal{D}$, and let us
define $A_{n}$ as previously. By definition, $\Phi(A_{n})>0$. The probability
to reach $A_{n}$ from the empty measure $\emptyset$ can be expressed as:

$\mathbb{P}^{n}\left(\emptyset,A_{n}\right)=\int_{y_{1}\in
B_{1},\dots,y_{n}\in
B_{n}}\int_{t_{1}=0}^{\infty}\dots\int_{t_{n}=t_{n-1}}^{\infty}\mspace{-12.0mu}H(\mathrm{d}y_{1})\dots
H(\mathrm{d}y_{n})e^{-\lambda p_{C_{1}}t_{1}}\dots e^{-\lambda
p_{C_{n}}(t_{n}-t_{n-1})}\prod_{i}\frac{b_{C_{i}}}{b_{C_{i}}+d_{C_{i}}}\lambda
p_{C_{1}}\mathrm{d}t_{1}\dots\lambda p_{C_{n}}\mathrm{d}t_{n},$

where $C_{i}$ denotes the class of the $i$-th user. This probability is
positive, which means that we can reach $A_{n}$ from the empty measure with
positive probability in $n$ steps, which gives us $\phi$-irreducibility for
the process $\Phi_{t}$. ∎

###  3.3 Stability

A first question of interest is that of the stability region. In the rest of
the paper, we use the following definition for stability:

######  Definition 3.1.

An SBD process is said to be _stable_ if the Markov chain $\Phi_{t}$ is
positive recurrent, and _unstable_ if it is null recurrent or transient.

We can now state the following theorem:

######  Theorem 3.3.

Under the foregoing assumptions, there exist two values
$0\leq\lambda_{c}^{-}\leq\lambda_{c}^{+}$ such that:

  * –

$\forall\lambda<\lambda_{c}^{-}$, $\Phi_{t}$ is stable;

  * –

$\forall\lambda>\lambda_{c}^{+}$, $\Phi_{t}$ is unstable.

###### Proof.

This theorem is a consequence of Theorem 3.1: if the system is stable for a
given $\lambda_{0}$, then it is stable for each $\lambda<\lambda_{0}$ and if
the system is unstable for a given $\lambda_{1}$, then it is unstable for all
$\lambda>\lambda_{1}$. The existence of cutoff (or critical) arrival rates is
well-established in monotonic queueing networks, and was adapted to spatial
birth-and-death processes (see Theorems 1 and 2 of [2] for instance). Let
$\lambda_{c}^{-}$ and $\lambda_{c}^{+}$ be defined as:

$\displaystyle\lambda_{c}^{-}$ $\displaystyle=\sup_{\lambda>0}\\{\lambda\text{
such that }\Phi_{t}\text{ is stable}\\},$ $\displaystyle\lambda_{c}^{+}$
$\displaystyle=\inf_{\lambda>0}\\{\lambda\text{ such that }\Phi_{t}\text{ is
unstable}\\}.$

$\Phi_{t}$ is an irreducible Markov chain on $\mathcal{M}(\mathcal{D})$. We
know that the Markov chain is either recurrent or transient (see Theorem 8.2.5
from [22]). Using stochastic monotonicity in the network from Theorem 3.1, we
know that the chain $\Phi_{t}$ is positive recurrent for all
$\lambda\leq\lambda_{c}^{-}$ and transient for all
$\lambda\geq\lambda_{c}^{+}$, which gives
$\lambda_{c}^{-}\leq\lambda_{c}^{+}$. Since the system is trivially stable for
$\lambda=0$ (in which case $\Phi_{t}$ is constantly equal to the null
measure), $\lambda_{c}^{-}$ and $\lambda_{c}^{+}$ are well defined and
$\lambda_{c}^{-}\geq 0$, which concludes the proof.

∎

We can refine this result by obtaining an upper and a lower bound bounds for
$\lambda_{c}$:

######  Lemma 3.4.

$\lambda_{c}^{-}$ and $\lambda_{c}^{+}$ satisfy:

$\frac{\ell(r)}{K\bar{L}\langle\ell_{\mathcal{D}}\rangle}\leq\lambda_{c}^{-}\leq\lambda_{c}^{+}\leq\frac{K\ell(r)}{\underline{L}\langle\ell_{\mathcal{D}}\rangle},$

where $\bar{L}=\max_{C}L_{C}$, $\underline{L}=\min_{C}L_{C}$ and
$\langle\ell_{\mathcal{D}}\rangle=\int_{x\in\mathcal{D}}\ell(\lVert
x\rVert)\mathrm{d}x$.

###### Proof.

The rate-of-transmission is bounded from below by that where all interfering
users transmit on all channels (i.e., are of class $[1,\dots,K]$) and the
transmitting user uses only one channel, i.e., for all $x,t$:

$R(x,\Phi_{t})\geq
R_{u}(x,\Phi_{t})\triangleq\frac{\ell(r)}{\mathcal{N}_{0}+\sum_{y\in\Phi_{t}\backslash\\{x\\}}K\ell(\lVert
x-y\rVert)}.$

Moreover, by definition, $\mathbf{L}<\bar{L}$ component-wise.

Conversely, the rate-of-transmission is bounded from above by that where
interfering users use a single channel and the transmitting users use all $K$
channels, i.e., for all $x,t$:

$R(x,\Phi_{t})\leq
R_{d}(x,\Phi_{t})\triangleq\frac{K\ell(r)}{\mathcal{N}_{0}+\sum_{y\in\Phi_{t}\backslash\\{x\\}}\ell(\lVert
x-y\rVert)},$

and $\underline{L}\leq\mathbf{L}$ component-wise. For any given initial
condition $\Phi_{0}$ and arrival rate $\lambda$, let
$\Phi_{u,t}=\Delta(\Phi_{0},\lambda,\underline{L},R_{u})$ and
$\Phi_{d,t}=\Delta(\Phi_{0},\lambda,\bar{L},R_{d})$. $\Phi_{u}$ and $\Phi_{d}$
two monotype dipolar Poisson networks, as defined in [2]. Using the main
result of [2], we know that the cut-off arrival rate is equal to
$\frac{\ell(r)}{K\bar{L}\langle\ell_{\mathcal{D}}\rangle}$ for $\Phi_{d,t}$
and to $\frac{K\ell(r)}{\underline{L}\langle\ell_{\mathcal{D}}\rangle}$ for
$\Phi_{u,t}$.

We apply Theorem 3.1, which states that $\Phi_{u,t}$ dominates $\Phi_{t}$ and
$\Phi_{d,t}$ is dominated by $\Phi_{t}$ to obtain the intended inequality.
Since $\frac{\ell(r)}{K\bar{L}\langle\ell_{\mathcal{D}}\rangle}>0$ and
$\frac{K\ell(r)}{\underline{L}\langle\ell_{\mathcal{D}}\rangle}<\infty$, we
can conclude that $0<\lambda_{c}^{-}$ and $\lambda_{c}^{+}<\infty$. ∎

##  4 Main results

In this section, we study the stability of the system through stochastic
domination and the use of a discretization of the dynamics of the SBD process.
We state the following theorem:

######  Theorem 4.1.

In the symmetric setup,
$\lambda_{c}^{-}=\lambda_{c}^{+}\triangleq\lambda_{c}$, where $\lambda_{c}$ is
the _critical_ arrival rate, equal to:

$\lambda_{c}=\frac{K\ell(r)}{\langle\ell_{\mathcal{D}}\rangle\mathfrak{L}},$
(7)

where $\mathfrak{L}\triangleq\sum_{C}p_{C}\lvert C\rvert L_{C}$.

The proof for this result happens in three steps. In a first step, we define a
discretization of the dynamics to reduce the study of the jump process
$\Phi_{t}$ to the study of two adequate queueing networks, one that bound the
dynamics from above and another on that bounds the dynamics from below. We
then move on to proving the stability of the former network and the
instability of the latter to obtain two bounds for the critical arrival rate
of the network, which coincide to the value of $\lambda_{c}$.

###  4.1 Discretization of the dynamics

To study our dynamics, we introduce the following discretization: for
$\varepsilon>0$, we tessellate $\mathcal{D}$ in $N_{\varepsilon}$ square cells
of side length $\varepsilon$, such that the origin is at the center of its
cell. We denote by $A_{i}$ the cell centered at $a_{i}\in\mathcal{D}$.
Finally, we introduce the stochastic process
$\mathbf{X}_{\varepsilon}(t)=(X_{i,C}(t))_{0\leq i\leq
N_{\varepsilon}-1,C\in\mathcal{P}(K)}$, where for all $i,C$,
$X_{i,C}(t)=\Phi_{t,C}(A_{i})$ is the number of receiver-transmitter pairs of
type $C$ in cell $i$.

We now use this discretization of the dynamics to define two _interference
queueing networks_ (as defined in [12]) $\bar{\mathbf{X}}_{\varepsilon}$ and
$\underline{\mathbf{X}}_{\varepsilon}$ with state space
$\mathbb{N}^{N_{\varepsilon}\times 2^{K}-1}$ such that:

  * –

$\bar{\mathbf{X}}_{\varepsilon}$ dominates $\mathbf{X}_{\varepsilon}$

  * –

$\underline{\mathbf{X}}_{\varepsilon}$ is dominated by
$\mathbf{X}_{\varepsilon}$.

Let $\bar{X}_{i,C}(t)$ and $\underline{X}_{i,C}(t)$ denote the respective
number of users of type $C$ in cell $i$ in each of the two processes. For a
given configuration of each process, we define the dynamics as follows:

  * –

The birth process of $\bar{\mathbf{X}}_{\varepsilon}$ is a Poisson rain of
intensity $\lambda\lvert\mathcal{D}\rvert$ and the death rate for users in
cell $i$ of class $C$ at time $t$ is
$\frac{\bar{X}_{i,C}(t)}{L_{C}}\bar{R}_{i,C}(t)$

  * –

The birth process of $\underline{\mathbf{X}}_{\varepsilon}$ is a Poisson rain
of intensity $\lambda\lvert\mathcal{D}\rvert$ and the death rate for users in
cell $i$ and class $C$ at time $t$ is
$\frac{\underline{X}_{i,C}(t)}{L_{C}}\underline{R}_{i,C}(t)$.

To define the functions $\bar{R}_{i,C}$ and $\underline{R}_{i,C}$, we start by
defining two path-loss functions $\ell^{\varepsilon}$ and $\ell_{\varepsilon}$
as follows:

$\begin{cases}\ell^{\varepsilon}(x,y)=\ell^{\varepsilon}(a_{i},a_{j})\hfill\text{
if }x\in A_{i},y\in A_{j}\\\
\ell^{\varepsilon}(a_{i},a_{j})=\max\left\\{\ell(\lVert
b_{i}-b_{j}\rVert),b_{j}\in\mathcal{V}_{j},b_{i}\in\mathcal{V}_{i}\right\\},&\end{cases}$

and

$\begin{cases}\ell_{\varepsilon}(x,y)=\ell_{\varepsilon}(a_{i},a_{j})\hfill\text{
if }x\in A_{i},y\in A_{j}\\\
\ell_{\varepsilon}(a_{i},a_{j})=\min\left\\{\ell(\lVert
b_{i}-b_{j}\rVert),b_{j}\in\mathcal{V}_{j},b_{i}\in\mathcal{V}_{i}\right\\},\end{cases}$

where, for $0\leq i\leq N_{\varepsilon}-1$, we define
$\mathcal{V}_{i}=\left\\{b_{i},\lVert
b_{i}-a_{i}\rVert\in\\{0,\varepsilon\\}\right\\}$. Using the dominated
convergence theorem, we get:

$\lim_{\varepsilon\rightarrow
0^{+}}\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle=\lim_{\varepsilon\rightarrow
0^{+}}\langle\ell_{\varepsilon,\mathcal{D}}\rangle=\langle\ell_{\mathcal{D}}\rangle,$
(8)

where $\langle\cdot_{\mathcal{D}}\rangle$ is defined in (1). Because of the
square torus topology of $\mathcal{D}$, we have, for all $i$:

$\sum_{k=0}^{N_{\varepsilon}-1}\ell^{\varepsilon}(a_{i},a_{k})=\sum_{k=0}^{N_{\varepsilon}-1}\ell^{\varepsilon}(0,a_{k})=\frac{1}{\varepsilon^{2}}\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle.$
(9)

The same result holds for $\ell_{\varepsilon}$.

The interferences $\bar{I}_{i,C}$ and $\underline{I}_{i,C}$ experienced by
users in cell $i$ of class $C$ and with path-loss functions
$\ell^{\varepsilon}$ and $\ell_{\varepsilon}$ are respectively equal to:

$\displaystyle\bar{I}_{i,C}(t)$ $\displaystyle=\sum_{k,U}\lvert C\cap
U\rvert\ell^{\varepsilon}(a_{k},a_{i})\left(\bar{X}_{k,U}(t)-\mathbf{1}_{\\{U=C,i=k\\}}\right)$
$\displaystyle\underline{I}_{i,C}(t)$ $\displaystyle=\sum_{k,U}\lvert C\cap
U\rvert\ell^{\varepsilon}(a_{k},a_{i})\left(\underline{X}_{k,U}(t)-\mathbf{1}_{\\{U=C,i=k\\}}\right).$

Finally, the rate-of-transmission functions of users in cell $i$ of class $C$
are defined by:

$\displaystyle\bar{R}_{i,C}(t)$ $\displaystyle=\frac{\lvert
C\rvert}{\mathcal{N}_{0}+\bar{I}_{i,C}(t)}$
$\displaystyle\underline{R}_{i,C}(t)$ $\displaystyle=\frac{\lvert
C\rvert}{\mathcal{N}_{0}+\underline{I}_{i,C}(t)}.$

######  Theorem 4.2.

$\bar{\mathbf{X}}_{\varepsilon}$ and $\underline{\mathbf{X}}_{\varepsilon}$
are two irreducible Markov jump processes with state space
$\mathbb{N}^{N_{\varepsilon}\times 2^{K}-1}$. Moreover,
$\bar{\mathbf{X}}_{\varepsilon}$ stochastically dominates
$\mathbf{X}_{\varepsilon}$, and $\underline{\mathbf{X}}_{\varepsilon}$ is
stochastically dominated by $\mathbf{X}_{\varepsilon}$.

###### Proof.

From the definition of the function $\ell^{\varepsilon}$ and
$\ell_{\varepsilon}$, we know that for all $x,y\in\mathcal{D}$, we have:

$\ell_{\varepsilon}(x,y)\leq\ell(\lVert
x-y\rVert)\leq\ell^{\varepsilon}(x,y).$

From this definition, we can obtain that for all $C\in\mathcal{P}(K)$, $0\leq
i\leq N_{\varepsilon}-1$, $x\in A_{i,C}$ and a given network configuration
$\Phi_{t}$, we have:

$\bar{R}_{i,C}(t)\leq R(x,\Phi_{t})\leq\underline{R}_{i,C}(t).$

To obtain the domination relation, we can now apply Theorem 3.1.
Irreducibility is obtained by developing the same argument as for Theorem 3.2,
which concludes the proof. ∎

Using Theorem 4.2, we can now reduce the study of the SBD process $\Phi_{t}$
to the study of the stability of the two processes
$\underline{\mathbf{X}}_{\varepsilon}$ and $\bar{\mathbf{X}}_{\varepsilon}$.
In the next section, we introduce a framework to obtain a condition on the
arrival rate such that the former is unstable, and the latter is stable.

###  4.2 Fluid limits and fluid model

From the definitions of the dynamics in our system, we can establish an
equation for the dynamics of $\bar{\mathbf{X}}_{\varepsilon}$, which will
allow us to obtain a condition on $\lambda_{c}^{-}$ for the stability of the
system.

Let $(\mathcal{A}_{i,C})$ and $(N_{i,C})$ be two families of independent
Poisson processes with intensity 1. By definition of the dynamics of the
discretized systems, arrivals of users in cell $i$ of class $C$ happen with
rate $\lambda p_{C}\varepsilon^{2}$. If the system is in state
$x=(x_{i,C})\in\mathbb{N}^{N_{\varepsilon}\times 2^{K}-1}$ at time $t$, users
in cell $i$ and class $C$ have a departure rate equal to
$\frac{1}{L_{C}}x_{i,C}(t)\bar{R}_{i,C}(t)$. We thus obtain the following
equation ruling the evolution of the population in the chain
$\bar{\mathbf{X}}_{\varepsilon}$:

$\displaystyle\bar{X}_{i,C}(t)=$
$\displaystyle\bar{X}_{i,C}(0)+\mathcal{A}_{i,C}(\lambda
p_{C}\varepsilon^{2}t)$ $\displaystyle-
N_{i,C}\left(\frac{1}{L_{C}}\int_{0}^{t}\bar{X}_{i,C}(u)\bar{R}_{i,C}(u)\mathrm{d}u\right).$
(10)

The goal of this section is to find a condition on $\lambda$ so that the
Markov chain $\bar{\mathbf{X}}_{\varepsilon}(t)$ is positive recurrent. Using
(10), we get the following theorem for fluid limits in our discretized
network:

######  Theorem 4.3.

The fluid limits for the chain $\bar{\mathbf{X}}_{\varepsilon}$ exist and are
solutions to the following system of equations:

$\begin{cases}\bar{x}^{\prime}_{i,C}(t)=\lambda
p_{C}\varepsilon^{2}-\frac{1}{L_{C}}\frac{\lvert
C\rvert\bar{x}_{i,C}(t)}{\sum_{k,U}\lvert C\cap
U\rvert\ell^{\varepsilon}(a_{k},a_{i})\bar{x}_{k,U}}\hfill\textrm{ if
}\bar{\mathbf{x}}(t)\neq 0\\\
\bar{\mathbf{x}}_{0}=\bar{\mathbf{X}}_{\varepsilon}(0)\end{cases}.$

The proof for this result is detailed in Appendix A. It is to note that these
dynamics are not well-defined in the case $\bar{\mathbf{x}}(t)\neq 0$, but we
do not need such a definition: the results about stability and fluid limits we
use involve fluid limits with a strictly positive initial condition. The
structure of this proof is adapted from [23], where the authors propose a
method to prove the existence and uniqueness of fluid limits for a given range
of dynamics. Fluid limits prove to be powerful tools in order to study the
stability and instability of queueing network.

###  4.3 Stability for symmetric multiclass wireless networks

We start by studying the stability of the chain
$\bar{\mathbf{X}}_{\varepsilon}$, which gives us an lower bound for
$\lambda_{c}^{-}$. We have the following theorem:

######  Theorem 4.4.

Let $\varepsilon>0$ and
$\bar{\lambda}_{\varepsilon}=\frac{K}{\mathfrak{L}\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle}$.
If $\lambda<\bar{\lambda}_{\varepsilon}$, then the chain
$\bar{\mathbf{X}}_{\varepsilon}$ is stable.

The proof for this result, relying on [19], is given in Appendix B.

We know from Theorem 4.2 that the process $\bar{\mathbf{X}}_{\varepsilon}$
dominates the original dynamic. From Theorem 3.1, we obtain that, for all
$\varepsilon>0$,
$\lambda_{c}\geq\frac{K}{\mathfrak{L}\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle}$.
Taking the limit as $\varepsilon$ goes to 0 gives us:

$\lambda_{c}^{-}\geq\frac{K}{\mathfrak{L}\langle\ell_{\mathcal{D}}\rangle}.$

The stability of the fluid scaled model is only a necessary condition for the
stability of the Markov chain and not a sufficient one: in [24], we have an
example of a stable queueing system with an unstable fluid model. An idea to
prove instability in the network would be to rely on the concept of _weak
instability_ , (see Definition 4.1 from [25]). Theorem 4.2 from the same paper
provides a result linking weak instability of a fluid limit model and the
transience of the associated Markov chain.

Unfortunately, to use this result, we have to study fluid limits starting from
0, and the dynamics of our fluid limit model are not defined when
$\bar{\mathbf{x}}=0$. To obtain instability for the system, we will use
stochastic domination and an adequate Markov chain to bound
$\underline{\mathbf{X}}$ from below. We have the following theorem:

######  Theorem 4.5.

Let
$\underline{\lambda}_{\varepsilon}=\frac{K}{\langle\ell_{\varepsilon,\mathcal{D}}\rangle\mathfrak{L}}$.
In the symmetric case, if $\lambda>\underline{\lambda}_{\varepsilon}$, then
$\underline{\mathbf{X}}_{\varepsilon}$ is unstable.

The proof of this theorem is presented in Appendix C. We know that the network
is unstable if
$\lambda>\frac{K}{\langle\ell_{\varepsilon,\mathcal{D}}\rangle\mathfrak{L}}$
for each $\varepsilon>0$. Taking the limit as $\varepsilon$ goes to 0, we
conclude:

$\lambda_{c}^{+}\leq\frac{K}{\langle\ell_{\mathcal{D}}\rangle\mathfrak{L}}.$

When combining the results of Theorem 4.4 and Theorem 4.5, we obtain that:

$\lambda_{c}^{-}=\lambda_{c}^{+}=\frac{K}{\langle\ell_{\mathcal{D}}\rangle\mathfrak{L}},$

which concludes the proof of Theorem 4.1 in the case $r=0$.

###  4.4 Generalization to $r>0$

To generalize the result of Theorem 4.1 to an arbitrary link length $r>0$, we
have to consider both the location of the transmitters and of the receivers in
the system, because a transmitter will not necessarily located in the same
cell as its receiver.

Let $\bar{\mathbf{Z}}_{\varepsilon}(t)=(Z_{i,C}(t))$ be the
$N_{\varepsilon}\times 2^{K}-1$ matrix of transmitter locations in the
dominating system and $\bar{\mathbf{M}}_{\varepsilon}(t)$ be such that
$\bar{M}_{i,C}(t)$ is a $N_{\varepsilon}\times 2^{K}-1$ matrix whose
coordinate $(j,D)$ denotes how many transmitters in cell $j$ of class $D$ have
a receiver in cell $i$ of class $C$ at time $t$ in the network.

The process
$\bar{\mathbf{S}}(t)=(\bar{\mathbf{X}}_{\varepsilon}(t),\bar{\mathbf{Z}}_{\varepsilon}(t),\bar{\mathbf{M}}_{\varepsilon}(t))$
is a Markov chain with countable state space. The evolution of
$\bar{\mathbf{S}}(t)$ is the following: to each receiver of class $C$ arriving
in cell $i$ with rate $\lambda p_{C}\varepsilon^{2}$, we pick uniformly a
point $x$ in the cell $A_{i}$. We then draw a circle with radius $r$ centered
at $x$ and we pick a point $y$ uniformly at random on the circle, which
decides the cell in which the transmitter is located. To obtain the
interference seen in the network, we sum over the locations of the
transmitters, i.e. the interference from (2) involves the vector
$\bar{\mathbf{Z}}_{\varepsilon}(t)$. When a receiver leaves the system, we
delete uniformly at random a transmitter such that its receiver is located in
the cell where the departure happened.

Using the definitions of the process $\bar{\mathbf{S}}(t)$ and the same steps
as for Theorem 8.3, we can obtain fluid limits for the process
$\bar{\mathbf{X}}_{\varepsilon}(t)$ in this setup (which is symmetric):

$\bar{x}^{\prime}_{i,C}(t)=\lambda
p_{C}\varepsilon^{2}-\frac{1}{L_{C}}\frac{\lvert
C\rvert\ell(r)\bar{x}_{i,C}(t)}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert
C\cap U\rvert\bar{z}_{k,U}(t)}.$

For all $i,C$ and $t\geq 0$, let $i^{\star}_{C}(t)$ be such that
$z_{i^{\star}_{C}(t),C}=\max_{i}z_{i,C}(t)$. The process
$(\bar{z}_{i^{\star}_{C}(t),C})$ dominates $\bar{\mathbf{Z}}_{\varepsilon}$,
and we have:

$\bar{z}^{\prime}_{i,C}(t)\leq\lambda
p_{C}\varepsilon^{2}-\frac{1}{L_{C}}\frac{\lvert
C\rvert\ell(r)\bar{z}_{i,C}(t)}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert
C\cap U\rvert\bar{z}_{k^{\star}_{U}(t),U}(t)}.$

We can now use similar arguments as for Theorem 4.4 to obtain the stability of
the Markov process $\bar{\mathbf{S}}(t)$ under the condition:

$\lambda_{c}\leq\frac{K\ell(r)}{\mathfrak{L}\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle}.$

To obtain the reciprocal, we define $i_{C}^{\dagger}(t)$ such that
$\underline{Z}_{i_{C}^{\dagger}(t),C}(t)=\min_{i}\underline{Z}_{i,C}(t)$. We
can then use the same arguments as developed for Theorem 4.5 to obtain that
the system is unstable if
$\lambda>\frac{K\ell(r)}{\mathfrak{L}\langle\ell_{\varepsilon,\mathcal{D}}\rangle}$.
This way, we get the reciprocal inequality and we can conclude that:

$\lambda_{c}=\frac{K\ell(r)}{\mathfrak{L}\langle\ell_{\mathcal{D}}\rangle},$

which concludes the proof of Theorem 4.1 in the case $r>0$.

##  5 Poisson heuristic for the critical arrival rate

In this section, we provide a heuristic the user densities in the stationary
regime in the system that will provide with a heuristic for the critical
arrival rate in the system as the largest possible solution of a given system
of equations.

###  5.1 Poisson heuristic for stationary user densities

Let
$\mu_{C}=\frac{1}{\lvert\mathcal{D}\rvert}\mathbb{E}\left[\Phi_{C,0}(\mathcal{D})\right]$
denote the user density of the point process $\Phi_{C,0}$ of users of class
$C\in\mathcal{P}(K)$ in the stationary regime. The first result we can state
about the user densities is the following :

######  Theorem 5.1.

In the symmetric system, the stationary user densities verify:

$\forall C,D\in\mathcal{P}(K),\quad\lvert C\rvert=\lvert
D\rvert\Rightarrow\mu_{C}=\mu_{D}.$

The result from Theorem 5.1 is a logical consequence of the symmetry
assumption, applied to the stationary regime.

A first heuristic to estimate $\mu_{C}$ is the Poisson heuristic: we assume
that, in the stationary regime, all the point processes $\Phi_{C,0}$ for
$C\in\mathcal{P}(K)$ are independent Poisson point processes with intensity
$\mu_{C}^{f}$. Poisson approximations are tied to replica mean-field methods
(see [16] for instance), where we consider multiple realizations of our
stochastic process and interactions between users are picked uniformly at
random among replicas of the dynamics. We know that (see [16]) as the number
of replicas go to infinity, the processes behave like independent Poisson
point processes, which allow to make this approximation. We define the
following heuristic:

######  Heuristic 1.

We define the _Poisson heuristic_ $\mu_{C}^{f}$ for $\mu_{C}$ as the smallest
solution of the following equation:

$\mu_{C}^{f}\int_{0}^{\infty}e^{-z\mathcal{N}_{0}}e^{-\sum_{U\in\mathcal{P}(K)}\mu_{U}^{f}\mathcal{I}(z,\lvert
C\cap U\rvert)}\mathrm{d}z=\frac{\lambda p_{C}L_{C}}{\lvert C\rvert\ell(r)}.$
(11)

where $\mathcal{I}(z,k)=\int_{\mathcal{D}}(1-e^{-zk\ell(\lVert
x\rVert)})\mathrm{d}x$.

We numerically observe that the values of the intensities $\mu_{C}^{f}$
obtained by solving (11) do not depend on the value of $\mathcal{N}_{0}>0$.
This observation is consistent with the fluid system and the stability
condition in the network being independent of the value of $\mathcal{N}_{0}$.

###### Proof.

Let us apply Miyazawa’s Rate Conservation Principle (cf [26]) to the process
$\Phi_{0,C}$ of users of class $C\in\mathcal{P}(K)$: during a time interval
$dt$, in the stationary regime, a $\lambda p_{C}L_{C}\lvert D\rvert dt$ users
arrive on average in the network, and in a given network configuration
$\Phi_{0,C}$,
$\mathbb{E}\left[\frac{1}{L_{C}}\sum_{x\in\Phi_{C,0}}R(x,\Phi_{0})dt\right]$
users leave, which gives:

$\lambda
p_{C}L_{C}\lvert\mathcal{D}\rvert=\mathbb{E}\left[\sum_{x\in\Phi_{C,0}}R(x,\Phi_{0})\right].$

We use the definition of the Palm probability measure to rewrite this as:

$\lambda
p_{C}L_{C}\lvert\mathcal{D}\rvert=\mathbb{E}^{0}_{\Phi_{C,0}}\left[R(0^{C},\Phi_{0})\right]\mathbb{E}\left[\Phi_{0}^{C}(\mathcal{D})\right].$
(12)

Using this, we get that, in the stationary regime:

$\displaystyle\lambda p_{C}L_{C}$
$\displaystyle=\mu_{C}\mathbb{E}^{0}_{\Phi_{C,0}}\left[R(0^{C},\Phi_{0})\right]$
$\displaystyle=\lvert
C\rvert\mu_{C}\ell(r)\mathbb{E}_{\Phi_{C,0}}^{0}\left[\frac{1}{N_{0}+I_{C}}\right].$

This is equivalent to:

$\mu_{C}=\frac{\lambda p_{C}L_{C}}{\lvert
C\rvert\ell(r)}\left(\mathbb{E}_{\Phi_{C,0}}^{0}\left[\frac{1}{N_{0}+I_{C}}\right]\right)^{-1}.$

We use the result of Lemma 1 from [2]:

######  Lemma 5.2.

Let $Y$ be a positive random variable with finite expectation and $c>0$ be a
real number. Then:

$\mathbb{E}\left[\frac{1}{c+Y}\right]=\int_{0}^{\infty}e^{-zc}\mathbb{E}[e^{-zY}]\mathrm{d}z.$

Using Lemma 5.2 with $Y\equiv I(0^{C},\Phi_{0})$, which has a finite first
moment and with $c\equiv\mathcal{N}_{0}$, we get:

$\mathbb{E}_{\Phi_{C,0}}^{0}\left[\frac{1}{\mathcal{N}_{0}+I_{C}}\right]=\int_{0}^{\infty}e^{-z\mathcal{N}_{0}}\mathbb{E}_{\Phi_{C,0}}^{0}\left[e^{-zI_{C}}\right]\mathrm{d}z.$

Let us assume that the processes $\Phi_{0,U}$ are independent Poisson point
processes denoted by $\Psi^{U}$. Using this, combined with Slivnyak’s theorem,
we get:

$\mathbb{E}_{\Phi_{C,0}}^{0}\left[\frac{1}{N_{0}+I_{C}}\right]\simeq\int_{0}^{\infty}e^{-z\mathcal{N}_{0}}\prod_{U\in\mathcal{P}(K)}\mathbb{E}\left[e^{-z\sum_{x\in\Phi_{0}^{U}\backslash\\{0^{C}\\}}\lvert
C\cap U\rvert\ell(\lVert x\rVert)}\right]\mathrm{d}z.$

Using the formula for the Laplace transform of a Poisson point process, we
get:

$\mathbb{E}\left[e^{-z\sum_{x\in\Phi_{0}^{U}\backslash\\{0^{C}\\}}\lvert C\cap
U\rvert\ell(\lVert
x\rVert)}\right]=\exp\left[-\mu_{U}^{f}\int_{\mathcal{D}}(1-e^{-z\lvert C\cap
U\rvert\ell(\lVert x\rVert)})\mathrm{d}x\right].$

To conclude the proof, we combine these equations:

$\mu_{C}^{f}\int_{0}^{\infty}e^{-z\mathcal{N}_{0}}e^{-\sum_{U\in\mathcal{P}(K)}\mu_{U}^{f}\mathcal{I}(z,\lvert
C\cap U\rvert)}\mathrm{d}z=\frac{\lambda p_{C}L_{C}}{\lvert C\rvert\ell(r)},$

which is the announced result. ∎

We can use the result from Theorem 5.1 to obtain another version of the system
of equations in (11): if we define
$\mu_{j}=\frac{1}{\lvert\mathcal{D}\rvert}\mathbb{E}\left[\sum_{U:\lvert
U\rvert=j}\Phi_{U,0}(\mathcal{D})\right]=\sum_{U:\lvert U\rvert=j}\mu_{U}$ as
the intensity of users communicating on $j$ channels, for $1\leq j\leq K$,
with $L_{j}=L_{C}$ and $p_{j}=\binom{K}{j}p_{C}$ for all $C$ such that $\lvert
C\rvert=j$, we obtain the following system of equations from Heuristic 2:

$\displaystyle\mu_{j}^{f}\int_{0}^{\infty}e^{-z\mathcal{N}_{0}-\sum_{l=1}^{K}\mu_{l}^{f}\sum_{m=1}^{j\wedge
l}\alpha_{m,j,l}\mathcal{I}(z,m)}\mathrm{d}z=\frac{\lambda
p_{j}L_{j}}{j\ell(r)},$ (13)

where $\alpha_{m,j,l}=\frac{\binom{j}{m}\binom{j-l}{l-m}}{\binom{K}{j}}$. This
result if a consequence of Lemma 8.5 applied to the vector of user densities.

This result can be interesting computationally, in the case where we want to
quickly estimate the density of users transmitting on a given class in the
symmetric system: in contrast to the system proposed in Heuristic 1, which
needs to solve $2^{K}-1$ equations, this symmetric system possesses $K$
equations. As $K$ grows larger, the gain in performance becomes non-
negligible.

###  5.2 A heuristic for $\lambda_{c}$

An interesting result arising from Heuristic 1 is an estimate for the value of
$\lambda_{c}$: a necessary condition for the fixed point equation (11) to
admit admit solutions is that the system is in it stationary regime, i.e., if
$\lambda<\lambda_{c}$. We define the following estimate for $\lambda_{C}$:

######  Heuristic 2.

Let $\lambda_{\textrm{P}}$ be the largest value of $\lambda$ such that (11)
admits a solution (or (13) in the symmetric network). Then,
$\lambda_{\textrm{P}}$ is an estimate for $\lambda_{c}$, which we call the
_Poisson heuristic_ estimate.

This estimate is interesting because it comes from a completely different
approach and is only relying on the study of the stationary regime of Poisson
approximates of the network. We will discuss in Section 6 the performance of
this heuristic.

The reason why this heuristic captures the stability region of our dynamics is
not yet understood. We do not know if the value of $\lambda_{\textrm{P}}$ is
equal to the value of $\lambda_{c}$ presented in Theorem 4.1 or if this value
is a precise numerical approximation of the critical arrival rate.

###  5.3 Second order heuristic

In an effort to refine the first order heuristic for user densities from
Heuristic 1, we introduce a _second order_ heuristic, which uses an
approximation for the pair-wise correlation function in the stationary regime.

######  Heuristic 3.

The _second-order heuristic_ for the intensity of $\Phi_{0,C}$ is
$\mu^{s}_{C}$ defined as:

$\mu^{s}_{C}=\frac{\lambda p_{C}L_{C}}{\lvert
C\rvert\ell(R)}\left(\mathcal{N}_{0}+I_{C}\right),$

where
$I_{C}\triangleq\mathbb{E}^{0}_{\Phi_{C,0}}\left[I(0^{C},\Phi_{0})\right]$ is
the interference experienced by the typical user of class $C$ in the
stationary regime. The vector $\mathbf{I}=(I_{C})_{C\in\mathcal{P}(K)}$ is the
solution to the following equation:

$I_{C}=\sum_{U\in\mathcal{P}(K)}\frac{\lvert C\cap
U\rvert}{\mu^{s}_{C}}\int_{x\in\mathcal{D}}\ell(\lVert
x\rVert)\rho^{(2)}_{C,U}(x,0)\mathrm{d}x,$

with the second order moment measure $\rho^{(2)}_{C,U}$ is defined as a
function of $\mathbf{I}$:

$\displaystyle\rho^{(2)}_{C,D}(x,y)=\frac{\lambda(\mu^{s}_{D}p_{C}+\mu^{s}_{C}p_{D})}{d_{C,D}(x,y)},$
(14)

and:

$\displaystyle d_{C,D}(x,y)=$ $\displaystyle\frac{1}{L_{C}}\frac{\lvert
C\rvert\ell(R)}{\mathcal{N}_{0}+\lvert C\cap D\rvert\ell\left(\lVert
x-y\rVert\right)+I_{C}}$ $\displaystyle+\frac{1}{L_{D}}\frac{\lvert
D\rvert\ell(R)}{\mathcal{N}_{0}+\lvert C\cap D\rvert\ell\left(\lVert
x-y\rVert\right)+I_{D}}.$

###### Proof.

To obtain this heuristic, we make the following cavity approximation in the
system: let us consider a pair of points $(x,y)$ in $\mathcal{D}$, where $x$
is of class $C$ and $y$ is of class $D$.

The arrival rate of the pair is $\lambda p_{C}\mu_{D}+\lambda p_{D}\mu_{C}$,
taking into account the contributions of the two processes. To obtain the
departure rate, let us fix the locations of the two points. The interference
experienced by a point located at $x$ of class $C$ in the stationary regime
conditioned on a user of class $D$ being present at $y$ is equal to:

$I(x)=\lvert C\cap D\rvert\ell(\lVert x-y\rVert)+I_{C},$ (15)

where $I_{C}$ is the interference experienced by the typical user of class $C$
in the network, i.e.
$I_{C}=\mathbb{E}^{0}_{\Phi_{0}}\left[I(0^{C},\Phi_{0})\right]$. Using the
definition of the second order moment measure relative to processes
$\Phi_{C,0}$ and $\Phi_{U,0}$, $\rho^{(2)}_{C,U}$, we have:

$\displaystyle I_{C}$
$\displaystyle=\mathbb{E}^{0}_{\Phi_{0}}\left[I(0^{C},\Phi_{0})\right]$
$\displaystyle=\sum_{U}\mathbb{E}^{0}_{\Phi_{0}}\left[I(0^{C},\Phi_{0}^{U})\right]$
$\displaystyle=\frac{1}{\mu_{C}}\sum_{U}\lvert C\cap
U\rvert\int_{x\in\mathcal{D}}\ell(\lVert
x\rVert)\rho^{(2)}_{C,U}(x,0)\mathrm{d}x.$

We obtain the departure rate for the pair $(x,y)$ with $x$ being of type $C$
and $y$ being of type $D$:

$\displaystyle d_{C,D}(x,y)=$ $\displaystyle\frac{1}{L_{C}}\frac{\lvert
C\rvert\ell(R)}{\mathcal{N}_{0}+\lvert C\cap D\rvert\ell\left(\lVert
x-y\rVert\right)+I_{C}}$ $\displaystyle+\frac{1}{L_{D}}\frac{\lvert
D\rvert\ell(R)}{\mathcal{N}_{0}+\lvert C\cap D\rvert\ell\left(\lVert
x-y\rVert\right)+I_{D}}.$

In the stationary regime, on average, the umber of arrivals compensate the
departures, which gives us:

$\lambda(\mu_{D}p_{C}+\mu_{C}p_{D})=\rho^{(2)}_{C,D}(x,y)d_{C,D}(x,y).$

Finally, to obtain the second-order heuristic for $\mu_{C}$, we use Miyazawa’s
Rate Conservation principle for the process $\Phi_{C,0}(\mathcal{D})$,
assuming that all the quantities are distributed according to their stationary
distribution:

$\lambda
p_{C}\lvert\mathcal{D}\rvert=\mathbb{E}\left[\Phi_{C,0}(\mathcal{D})\right]\frac{\lvert
C\rvert\ell(r)}{\mathcal{N}_{0}+I_{C}},$

which gives the intended result.

∎

##  6 Numerical simulations and performance of the heuristics

###  6.1 Stability region

In this section, we discuss the performance of the two heuristics developed in
Heuristics 1 and 3. In 3, a symmetric setup with 2 channels and three classes
denoted as $\\{1\\}$, $\\{2\\}$ and $\\{1,2\\}$ is considered. The average
file sizes are $L_{\\{1\\}}=L_{\\{2\\}}=1$ and $L_{\\{1,2\\}}=2$ and the
probability distribution is $p_{\\{1\\}}=p_{\\{2\\}}=0.4$ and and
$p_{\\{1,2\\}}=0.2$. The path-loss function is $\ell(x)=(1+x)^{-4}$, and the
domain $\mathcal{D}$ is the square torus centered at the origin with side
length 10.

To simulate the dynamics of the system, we simulate the embedded Markov chain
in the process $\Phi_{t}$: at each iteration, we compute the time to the next
birth by drawing a time $t_{b}$ from an exponential distribution with mean
$\lambda\lvert\mathcal{D}\rvert$ and a random vector $\mathbf{t}_{D}$ where
each coordinate corresponds to a user located at $x_{i}$ of class $C_{i}$ and
is drawn from an exponential distribution with mean
$\frac{1}{L_{C_{i}}}R(x_{i},\Phi_{t})$.

We then compare the values of $t_{b}$ and the minimum of $\mathbf{t}_{D}$: if
the former is lower, a user arrives at a location $x$ taken uniformly at
random in $\mathcal{D}$; if the latter is lower, the corresponding user leaves
the system. We then update the value of the interference in the system and
move forward in time to the next event. If $(T_{0}=0,T_{1},\dots,T_{n})$ are
the event times in the network and $X$ is the embedded Markov chain, we obtain
that $\Phi_{T_{n}}=X_{n}$, which gives us the representation of process
$\Phi_{t}$. Figure 1 shows the behavior of the system with $K=2$ for two
values of $\lambda$: one that shows stability, and the other, instability.

Figure 1: Number of active users in the network over time with $K=2$ bands for
two values of $\lambda$. On the left, $\lambda=0.95\lambda_{c}$ and we can see
that the number of users in the network stays bounded. On the right,
$\lambda=1.05\lambda_{c}$. In this case, population grows linearly over time
in the network, showing the instability of the system.

A reliable criterion to study the stability of the system is Little’s law (see
[27]), that links the average staying time $W$ of a user in the system, the
arrival rate $\lambda$ and the number of users $L$ in the system in the
stationary regime:

$L=\lambda W.$ (16)

Using the parameters for our system, (16) becomes:

$\mathbb{E}\left[\Phi_{t,C}(\mathcal{D})\right]=\lambda p_{C}W_{C},$ (17)

where $W_{C}$ is the staying time of a user of class $C$ and
$\mathbb{E}\left[\Phi_{t,C}(\mathcal{D})\right]$ is the average number of
users of class $C$ in the stationary regime. When simulating the dynamics of
the network, we can compute the average time spent by users in the network in
each class. If the value converges to a fixed, finite value, then the network
is stable and the number of users follow (17). If the time lived in the
network grows linearly and diverges to infinity, then the network is not
stable.

Figure 2 shows the average staying time of users in the network for two
different values of $\lambda$. The phase transition in the network happens as
expected: for $\lambda=0.9\lambda_{c}$, the system appears to be stable and
for $\lambda=1.1\lambda_{c}$, the average staying time grows linearly, meaning
that no stationary regime exists for in the network, and that our dynamics are
unstable.

Figure 2: Average staying time of users in the network with $K=2$ channels for
two different values of $\lambda$. In blue and red, users transmitting only on
one channel, in green, users transmitting on the two channels. On the left,
$\lambda=0.9\lambda_{c}$, on the right $\lambda=1.1\lambda_{c}$. The left
network is stable, and the average staying time converges to a finite limit.
The right network is unstable and the staying time diverges.

###  6.2 Poisson heuristics

We can compare the heuristics from Heuristic 1 and 2 to the value we obtain
through simulation. To compute this intensity, we use an ergodic approximation
of the form:

$\mu_{C}\approx\frac{1.28}{t}\int_{u}^{u+t}\frac{\Phi_{w,C}(\mathcal{D})}{\lvert\mathcal{D}\rvert}\mathrm{d}w,$
(18)

where $u$ is taken sufficiently large to ensure we are in the stationary
regime of the system and $t$ is large enough so that the approximation is
correct. The 1.28 factor arising in (18) comes from a Palm bias: in the
simulation, when we compute the average number of users in $\mathcal{D}$, we
introduce an observation bias when estimating the stationary user density,
which we correct by using this multiplicative factor.

Figure 3 displays the results of both the Poisson (in red, dashed) against the
numerical estimation of the user density (in blue, plain) for the process
$\Phi_{\\{1\\}}$ using (18).

Figure 3: Value of $\mu_{\\{1\\}}$ as a function of $\lambda/\lambda_{c}$ in a
symmetric configuration.

As we expected from previous results (see Section IV from [2]), the Poisson
heuristic for spatial wireless networks does not estimate precisely the
stationary user densities. This result can be explained by the existence of
_clustering_ , or attraction, in the stationary system. Clustering implies
that the spatial positions of users in the stationary regime are not
independent, and that connecting users tend to attract other users by slowing
communications of users that stay too close to them. A sensible conjecture to
make is that our spatial multiclass wireless network displays the same
behavior in the stationary regime, which would provide an explanation for this
observation.

The Poisson estimate gives us an interesting result: Figure 4 shows the
respective values of $\lambda_{c}$ and $\lambda_{\textrm{P}}$ in the symmetric
case for different values of $p_{\\{1,2\\}}$, with
$p_{\\{1\\}}=p_{\\{2\\}}=\frac{1-p_{\\{1,2\\}}}{2}$. The two values are very
close, which means that the Poisson hypothesis we make in the stationary
regime is a reasonable approximation to obtain a condition for the stability
of the system. In other words, Heuristic 2 provides an accurate numerical
heuristic for the result of Theorem 4.1.

Figure 4: Comparison between $\lambda_{c}$ and its Poisson heuristic as a
function of $p_{\\{1,2\\}}$ in the network in the symmetric case.

##  7 Discussion and extensions

This paper proposed a first step in the analysis of a multiclass interacting
point process dynamics representing important new features of wireless
networks, and which extends the single-class model of [2] by adding service
differentiation. Service differentiation changes quite significantly the
interactions between users, and new ideas were needed to assess the
monotonicity, the $\phi$-irreducibility, the stochastic stability, and the
steady state properties for the associated dynamics.

Here are, to conclude, a few interpretation comments and natural extensions of
this first step.

###  7.1 Interpretation of the results

Theorem 4.1 gives the critical arrival rate below which the dynamics is stable
and above which it is not. The form of the critical value in question is
somewhat unexpected. The quantity $\mathfrak{L}=\sum_{U}\lvert U\rvert
p_{U}L_{U}$ which shows up in this critical value can be given the following
interpretation: on average, a user has $L_{U}$ bits to transmit over $\lvert
U\rvert$ bands. Hence, a user of type $U$ brings a time-space load equal to
$\lvert U\rvert L_{U}$ to the network. Taking the average over the
distribution of user types, we obtain $\mathfrak{L}$. We can rewrite the
stability condition
$\lambda<\frac{K\ell(R)}{\langle\ell_{\mathcal{D}}\rangle\mathfrak{L}}$ as:

$\lambda\frac{\mathfrak{L}}{K}<\frac{\ell(R)}{\langle\ell_{\mathcal{D}}\rangle},$

We can interpret the left hand side term as the average load arriving per unit
of time and per channel. The right-hand side is the throughput capacity of one
channel in heavy load (the signal power divided by the fluid interference
power).

###  7.2 Possible extensions

In this paper, we made the assumption that the system was in a low-SINR
regime, which justified the linearization of the Shannon formula in (3). This
is not an essential limitation and it is not difficult to show that the
stability results we present here remain valid, up to a multiplicative factor,
when using the non-linearized version of the Shannon formula.

The methodology can also be extended and applied to more complex network
settings.

A first natural extension concerns the relaxation of the symmetry hypothesis
we made. Relaxing this assumption prevents the use of Lemma 8.5 and the
stability condition of Theorem 4.1 does not work anymore. For a non-symmetric
multiclass network, it is not clear whether we will get a closed-form for the
critical arrival rate.

A second question comes from the Poisson heuristic we defined in Heuristic 2.
This heuristic gives a precise estimate for the value of $\lambda_{c}$ in the
symmetric system. So far, we cannot say whether this estimate is a precise
numerical approximation or an exact result for $\lambda_{c}$. A first step to
answer this question would be to study the spatial correlation of points when
reaching the limit of the stability region: if this spatial correlation
vanishes, the Poisson approximation becomes relevant and this could explain
why this heuristic performs well. A result supporting this interpretation
comes from Section V of [2], where the authors showed that clustering becomes
less prominent when $\lambda$ becomes closer to $\lambda_{c}$. The question
whether this extends to the multiclass setting is open. It is also worth
noting that in simulations not presented in this paper, the Poisson heuristic
close to $\lambda_{c}$ still performs well when removing the symmetry
hypothesis, and in different numerical setups as the one presented in Section
6.

Among the other interesting properties we can expect to extend from single to
multiclass, let us quote the existence and quantification of clustering (see
Definition 1 of [2]), the definition of more precise stationary user density
heuristics and the study of higher order moment measures. We already know
that, for single class interference queueing networks, the system has
exponential moments in the stationary regime (see [12]). Extending the result
to our multiclass dynamics willl help quantifying the tail distribution of the
service time, which is critical for 5G networks.

Lastly, another question would be the extension of our results to the whole
Euclidean plane. This will however require completely different proof
techniques as there is no theory for assessing the stability of infinite
dimensional systems using fluid limits.

## Acknowledgements

This work was supported by the ERC NEMO grant, under the European Union’s
Horizon 2020 research and innovation programme, grant agreement number 788851
to INRIA. The authors thank S. Foss and S. Shneer of Heriot Watt University
for their comments and suggestions on this work.

## References

  * [1] C. E. Shannon, “A mathematical theory of communication,” The Bell System Technical Journal, vol. 27, 1948.
  * [2] A. Sankararaman and F. Baccelli, “Spatial birth-and-death wireless networks,” IEEE Transactions on Information Theory, vol. 63, no. 6, pp. 3964 – 3982, 2017.
  * [3] White paper, “Bandwidth part adaptation: 5G NR user experience and power consumption enhancements.” https://newsletter.mediatek.com/hubfs/mwc/download/bandwidth-part-adaptation.pdf, 2015\.
  * [4] F. Baccelli and S. Kalamkar, “Bandwidth allocation and service differentiation in d2d wireless networks,” in Proceedings of the IEEE International Conference on Computer Communications (INFOCOM), Toronto, Canada, July 2020.
  * [5] F. Baskett, K. M. Chandy, R. R. Muntz, and F. G. Palacios, “Open, closed, and mixed networks of queues with different classes of customers,” J. ACM, vol. 22, p. 248–260, apr 1975.
  * [6] F. Kelly, Reversibility and Stochastic Networks.  New-York: Wiley, 1979.
  * [7] F. Kelly, “Charging and rate control for elastic traffic,” European Transactions on Telecommunications, vol. 8, no. 1, pp. 33–37, 1997.
  * [8] L. Massoulie and J. Roberts, “Bandwidth sharing: objectives and algorithms,” in IEEE INFOCOM ’99. Conference on Computer Communications., vol. 3, pp. 1395–1403 vol.3, 1999.
  * [9] F. Baccelli and B. Blaszczyszyn, Stochastic Geometry and Wireless networks.  NoW Publishers, 2009.
  * [10] M. Haenggi, “The meta distribution of the SIR in poisson bipolar and cellular networks,” IEEE Transactions on Wireless Communications, vol. 15, no. 4, pp. 2577–2589, 2016.
  * [11] F. Baccelli and S. Foss, “Poisson hail on a hot ground,” Journal of Applied Probability, vol. 48, no. A, p. 343–366, 2011.
  * [12] A. Sankararaman, F. Baccelli, and S. Foss, “Interference queueing networks on grids,” The Annals of Applied Probabilities, vol. 29, no. 5, pp. 2928–2987, 2019.
  * [13] A. N. Rybko and A. Stolyar, “Ergodicity of stochatsic processes describing the operation of open queueing networks,” Problems of Information Transmission, vol. 28, no. 3, pp. 199–221, 1992.
  * [14] J. G. Dai, “On positive harris recurrence of multiclass queueing networks: a unified approach via fluid limit methods,” The Annals of Applied Probability, vol. 5, no. 1, pp. 49–77, 1995.
  * [15] F. Baccelli and P. Brémaud, Elements of Queueing Theory: Palm Martingale Calculus and Stochastic Recurrences.  Springer, 2003.
  * [16] F. Baccelli, M. Davydov, and T. Taillefumier, “Replica-mean-field limits of fragmentation-interaction-aggregation processes,” Journal of Applied Probability, pp. 1–22, 2022.
  * [17] S. Allmeier and N. Gast, “Mean field and refined mean field approximations for heterogeneous systems: It works!,” Proc. ACM Meas. Anal. Comput. Syst., vol. 6, no. 1, 2022.
  * [18] N. Gast, L. Bortolussi, and M. Tribastone, “Size expansions of mean field approximation: Transient and steady-state analysis,” Performance Evaluation, vol. 129, pp. 60–80, 2019.
  * [19] S. Shneer and A. Stolyar, “Stability conditions for a decentralised medium access algorithm: single- and multi-hop networks,” Queueing Systems, vol. 94, pp. 109–128, 2020.
  * [20] P. Brémaud and L. Massoulié, “Stability of nonlinear Hawkes processes,” The Annals of Probability, vol. 24, no. 3, pp. 1563 – 1588, 1996.
  * [21] P. Billingley, Convergence of Probability Measures.  Wiley, 2nd ed., 1999.
  * [22] S. P. Meyn and R. L. Tweedie, Markov Chains and stochastic stability.  Springer Science and Business Media, 2012.
  * [23] J. Reed and B. Zwart, “Limit theorems for markovian bandwidth sharing networks with rate constraints,” Operations Research, 2014.
  * [24] M. Bramson, “A stable queueing network with unstable fluid model,” The Annals of Applied Probability, vol. 9, no. 3, pp. 818–853, 1999.
  * [25] J. G. Dai, “A fluid limit model criterion for instability of multiclass queueing networks,” The Annals of Applied Probability, vol. 6, no. 3, pp. 751–757, 1996.
  * [26] M. Miyazawa, “Rate conservation laws : A survey,” Queueing Systems, vol. 15, pp. 15–58, 1994.
  * [27] J. D. C. Little, “A Proof for the Queueing Formula: $L=\lambda W$,” Operations Research, vol. 9, no. 3, pp. 383–387, 1961.
  * [28] J. Jacod and A. N. Shiryaev, Limit Theorems for Stochastic Processes, vol. 288.  Springer-Verlag, 1987.  Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].
  * [29] W. Whitt, Stochastic-process limits: an introduction to stochastic-process limits and their application to queues.  Springer-Verlag, 2002.

##  8 Appendix

## A: Proof of Theorem 4.3

Let us take a sequence of initial conditions
$\bar{\mathbf{X}}_{\varepsilon}^{n}(0)=(\bar{x}^{n}_{i,C})$ for our system,
with $\lim_{n\rightarrow\infty}\bar{x}^{n}_{i,C}=\infty$. The goal in this
section is to study the limit of the sequence of fluid-scaled processes
$\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(nt)\right)$ and to obtain a
system of ODEs for which the limit of this sequence is a solution.

Using (10), we get, for all $i,C$:

$\displaystyle\frac{1}{n}\bar{X}_{i,C}(nt)=\frac{\bar{x}^{n}_{i,C}}{n}+\frac{1}{n}\mathcal{A}_{i,C}(\lambda
p_{C}\varepsilon^{2}nt)-\frac{1}{n}N_{i,C}\left(\frac{1}{L_{C}}\int_{0}^{nt}\bar{X}_{i,C}(u)R_{i,C}(\bar{\mathbf{X}}_{\varepsilon}(u))\mathrm{d}u\right).$
(19)

We use the variable change $u=ns$ in the integral to get:

$\displaystyle\int_{0}^{nt}\bar{X}_{i,C}(u)R_{i,C}(\bar{\mathbf{X}}_{\varepsilon}(u))\mathrm{d}u=\int_{0}^{t}\frac{n\lvert
C\rvert\frac{1}{n}\bar{X}_{i,C}(ns)\mathrm{d}s}{\frac{\mathcal{N}_{0}}{n}+\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert
C\cap
U\rvert\left(\frac{1}{n}\bar{X}_{k,U}(ns)-\frac{1}{n}\mathbf{1}_{\\{U=C,i=k\\}}\right)}.$

We define:

$R_{i,C}^{n}(x)=\frac{\lvert
C\rvert}{\frac{\mathcal{N}_{0}}{n}+\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert
C\cap U\rvert\left(x_{k,U}-\frac{1}{n}\mathbf{1}_{\\{U=C,i=k\\}}\right)}.$

Using this, (19) becomes:

$\frac{1}{n}\bar{X}_{i,C}(nt)=\frac{1}{n}\bar{x}^{n}_{i,C}+\frac{1}{n}\mathcal{A}_{i,C}(\lambda
p_{C}\varepsilon^{2}nt)\hfill\\\
-\frac{1}{n}N_{i,C}\left(\frac{n}{L_{C}}\int_{0}^{t}\frac{1}{n}\bar{X}_{i,C}(ns)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(ns)\right)\mathrm{d}s\right),$

Let $M^{n}_{i,C}(z)=\frac{1}{n}N_{i,C}(nz)-z$ for $z\in\mathbb{R}$, and let us
define:

$\bar{Y}^{n}_{i,C}(t)=\frac{1}{n}\bar{x}^{n}_{i,C}+\frac{1}{n}\mathcal{A}_{i,C}(\lambda
p_{C}\varepsilon^{2}nt)-M^{n}_{i,C}\left(\frac{1}{L_{C}}\int_{0}^{t}\frac{1}{n}\bar{X}_{i,C}(ns)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(ns)\right)\mathrm{d}s\right).$

From (19), we have:

$\frac{1}{n}\bar{X}_{i,C}(nt)=\bar{Y}^{n}_{i,C}(t)-\frac{1}{L_{C}}\int_{0}^{t}\frac{1}{n}\bar{X}_{i,C}(ns)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(ns)\right)\mathrm{d}s.$

Let us denote
$\bar{\mathbf{x}}_{0}=\lim_{n\rightarrow\infty}\frac{1}{n}\bar{x}^{n}_{i,C}$
when it exists. We say that $\bar{\mathbf{x}}_{0}$ is finite if and only if
all its coordinates are finite.

######  Lemma 8.1.

If $\bar{\mathbf{x}}_{0}$ exists and is finite, then:

$\forall t\geq 0,\forall
i,C,\quad\lim_{n\rightarrow\infty}\bar{Y}^{n}_{i,C}(t)=\bar{X}_{i,C}(0)+\lambda
p_{C}\varepsilon^{2}t.$ (20)

###### Proof.

Using the strong law of large numbers yields, $\mathbb{P}$-almost surely:

$\lim_{n\rightarrow\infty}\frac{1}{n}\mathcal{A}_{i,C}(\lambda
p_{C}\varepsilon^{2}nt)=\lambda p_{C}\varepsilon^{2}t.$

Moreover, we have:

$\displaystyle\frac{1}{L_{C}}\int_{0}^{t}\frac{1}{n}\bar{X}_{i,C}(ns)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(ns)\right)\mathrm{d}s$
$\displaystyle=\frac{1}{L_{C}}\int_{0}^{t}\frac{\lvert
C\rvert\frac{1}{n}\bar{X}_{i,C}(ns)}{\frac{N_{0}}{n}+\sum_{k,U}\lvert C\cap
U\rvert\ell^{\varepsilon}(a_{k},a_{i})\frac{1}{n}\bar{X}_{k,U}(ns)}\mathrm{d}s$
$\displaystyle\leq\frac{1}{L_{C}}\int_{0}^{t}\frac{\lvert
C\rvert\frac{1}{n}\bar{X}_{i,C}(ns)}{\lvert
C\rvert\ell^{\varepsilon}(a_{i},a_{i})\frac{1}{n}\bar{X}_{i,C}(ns)}\mathrm{d}s$
$\displaystyle=\frac{t}{L_{C}},$

where we use the fact that for all $0\leq i\leq N_{\varepsilon}-1$,
$\ell_{\varepsilon}(a_{i},a_{i})=1$.

We use Prohorov’s Theorem (see Theorems 5.1 and 5.2 of [21]) to obtain that
the sequence of processes $\left\\{\frac{1}{n}\bar{X}_{i,C}(n\cdot)\right.$
$\left.R^{n}_{i,C}(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot))\right\\}$
is tight. Finally, because $M^{n}_{i,C}(z)\rightarrow 0$ as $n$ goes to
infinity for all $z$ (using the strong law of large numbers), we can conclude
that:

$\lim_{n\rightarrow\infty}M^{n}_{i,C}\left(\frac{1}{L_{C}}\int_{0}^{t}\frac{1}{n}\bar{X}_{i,C}(ns)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(ns)\right)\mathrm{d}s\right)=0,$

which leads to the intended result. ∎

For $v:[0,T)\rightarrow\mathbb{R}$ and $0\leq t_{1}\leq t_{2}\leq T$, we
define:

$w(v,[t_{1},t_{2}])=\sup\left\\{\lvert
v(u_{1})-v(u_{2})\rvert,u_{1},u_{2}\in[t_{1},t_{2}]\right\\}.$

Let us define the modulus of continuity
$\omega(v,\delta,T)=\sup\\{w(v,[s,t]),0\leq s,t\leq T,\lvert
s-t\rvert<\delta\\}$. We introduce the concept of _C-tightness_ and a useful
characterization (see Definition VI.3.25 and Proposition VI.3.26 from [28]):

######  Definition 8.1 (C-tightness).

A sequence $\\{V^{n},n\geq 1\\}$ of functions is _C-tight_ if and only if, for
all $t_{0}>0$, $\eta>0$, $T>0$, there exists $K^{0}_{\eta}$, $n^{0}_{\eta}$
and $\delta^{0}_{\eta}$ such that, for all $n\geq n^{0}_{\eta}$:

  * i)

$\mathbb{P}\left[\sup_{0\leq t\leq T}\lvert V^{n}(t)\rvert\geq
K^{0}_{\eta}\right]<\eta$.

  * ii)

$\mathbb{P}\left[\omega(V^{n},\delta^{0}_{\eta},t_{0})>\eta\right]<\eta$.

Condition $ii)$ implies that any limit point of $V^{n}$ has continuous sample
paths, $\mathbb{P}$-almost surely: let $V$ be a limit point of the sequence
$V^{n}$. By definition,
$\mathbb{P}\left[\omega(V,\delta^{0}_{\eta},t_{0})>\eta\right]<\eta$ for all
$\eta>0$, i.e., $\mathbb{P}$-almost surely, $\omega(V,\delta^{0},t_{0})>\eta$
for each $T>0$ and $\delta_{C}>0$. By continuity of the function
$\delta\mapsto\omega(V,\delta,t_{0})$ at 0, we can conclude that $V$ is
$\mathbb{P}$-almost surely continuous.

We can prove the following lemma:

######  Lemma 8.2.

If $\bar{\mathbf{x}}_{0}$ exists and is finite, the sequence of processes
$\\{\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot),n\geq 0\\}$ is C-tight.

###### Proof.

Let $T>0$ and $\eta>0$. As $n$ goes to infinity, we know that
$\bar{Y}^{n}_{i,C}(t)\rightarrow\bar{X}_{i,C}(0)+\lambda
p_{C}\varepsilon^{2}t$ $\mathbb{P}$-almost surely from Lemma 8.1. For each
$\eta>0$, there exists $n^{T,1}_{\eta}$ such that for $n\geq n^{T,1}_{\eta}$:

$\mathbb{P}\left[\sup_{0\leq t\leq
T}\lvert\bar{Y}^{n}_{i,C}(t)\rvert\geq\lambda
p_{C}\varepsilon^{2}T+1\right]\leq\eta.$

We also know that
$\frac{1}{n}\bar{X}_{i,C}(nt)R_{i,C}^{n}(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(nt))\leq
1$ at all times $t$, a fortiori for all $t\leq T$. We deduce that for all
$n\geq n^{T,1}_{\eta}$:

$\mathbb{P}\left[\sup_{0\leq t\leq
T}\left\lvert\frac{1}{n}\bar{X}_{i,C}(nt)\right\rvert\geq\lambda
p_{C}\varepsilon^{2}T+1+1\right]\leq\eta.$

We set $K^{0}_{\eta}=\lambda p_{C}\varepsilon^{2}T+2$ to obtain condition
$i)$.

To obtain condition $ii)$, we remark that proving the continuity of the limits
is equivalent to showing that there exists a $\delta^{0}_{\eta}>0$ such that
for all $i,C$ and at all times $t<T$, we have:

$\mathbb{P}\left[\omega(\frac{1}{n}\bar{X}_{i,C}(nt),\delta^{0}_{\eta},t_{0})>\eta\right]<\eta.$

Let $\delta>0$. We have at all times $t\geq 0$:

$\frac{1}{n}\bar{X}_{i,C}(n(t+\delta))-\frac{1}{n}\bar{X}_{i,C}(nt)=\bar{Y}^{n}_{i,C}(t+\delta)-\bar{Y}^{n}_{i,C}(t)\hfill\\\
-\frac{1}{L_{C}}\int_{t}^{t+\delta}\frac{1}{n}\bar{X}_{i,C}(nu)R_{i,C}^{n}\left(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(nu)\right)\mathrm{d}u.$

Taking the supremum over $t\in[0,T]$, and reminding that
$\frac{1}{n}\bar{X}_{i,C}(nt)R_{i,C}^{n}(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(nt))\leq
1$ yields:

$\sup_{0\leq t\leq
T}\left\lvert\frac{1}{n}\bar{X}_{i,C}(n(t+\delta))-\frac{1}{n}\bar{X}_{i,C}(t)\right\rvert\leq\sup_{0\leq
t\leq
T}\left\lvert\bar{Y}^{n}_{i,C}(t+\delta)-\bar{Y}^{n}_{i,C}(t)\right\rvert+\frac{\delta}{L_{C}}.$

From Lemma 8.1, we know that
$\bar{Y}^{n}_{i,C}(t)\rightarrow\bar{X}_{i,C}(0)+\lambda
p_{C}\varepsilon^{2}t$ $\mathbb{P}$-almost surely. This implies that there
exists $n^{T,2}_{\eta}>0$ and $\kappa_{\eta}^{0}>0$ such that, with
probability at least $1-\eta$,
$\lvert\bar{Y}^{n}_{i,C}(t+\delta)-\bar{Y}^{n}_{i,C}(t)\rvert<\lambda
p_{C}\varepsilon^{2}\delta+\delta\kappa_{\eta}^{0}$.

This implies that for all $i,C$, we have, with probability at least $1-\eta$:

$\omega\left(\frac{1}{n}\bar{X}_{i,C},\delta,t_{0}\right)<\delta\left(\lambda
p_{C}\varepsilon^{2}+\kappa_{\eta}^{0}+\frac{1}{L_{C}}\right).$

Let us set $\delta^{0}_{\eta}=\eta\left(\lambda
p_{C}\varepsilon^{2}+\kappa_{\eta}^{0}+\frac{1}{L_{C}}\right)^{-1}$. Thus, we
get, for all $n\geq n^{T,2}_{\eta}$:

$\mathbb{P}\left[\omega\left(\frac{1}{n}\bar{X}_{i,C},\delta^{0}_{\eta},t_{0}\right)<\eta\right]>1-\eta.$

Setting $n^{0}_{\eta}=\max(n^{1,0}_{\eta},n^{2,0}_{\eta})$ concludes the proof
of the C-tightness of the sequence of processes
$(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(t))$.

∎

We can now establish the equation ruling the evolution of the fluid scaled
model:

######  Theorem 8.3.

If
$\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(0)\rightarrow\bar{\mathbf{x}}_{0}$
as $n$ goes to $\infty$, then the sequence of processes
$\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot)$ converges
$\mathbb{P}$-almost surely to
$\bar{\mathbf{x}}_{\varepsilon}(s)=(\bar{x}_{i,C}(s))$, which is the unique
solution of the following system of differential equations:

$\begin{cases}\bar{x}^{\prime}_{i,C}(t)=\lambda
p_{C}\varepsilon^{2}-\frac{1}{L_{C}}\frac{\lvert
C\rvert\bar{x}_{i,C}(t)}{\sum_{k,U}\lvert C\cap
U\rvert\ell^{\varepsilon}(a_{k},a_{i})\bar{x}_{k,U}}\hfill\textrm{ if
}\bar{\mathbf{x}}(t)\neq 0\\\
\bar{\mathbf{x}}_{0}=\bar{\mathbf{X}}_{\varepsilon}(0)\end{cases}.$

###### Proof.

From Lemma 8.2, the sequences
$\\{\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot),n\geq 1\\}$ and
$\\{\bar{\mathbf{Y}}^{n},n\geq 1\\}$ are both tight. It follows from Theorem
11.6.8 from [29] that the sequence
$\\{(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot),\bar{\mathbf{Y}}^{n}),n\geq
1\\}$ is tight in $D([0,\infty],\mathbb{R}^{N_{\varepsilon}\times 2^{K}-1})$,
and thus, by using Prohorov’s Theorem, relatively compact. Let ${n_{l}}$ be a
subsequence along which
$\left(\frac{1}{n_{l}}\bar{\mathbf{X}}_{\varepsilon}(n_{l}\cdot),\bar{\mathbf{Y}}^{n_{l}}\right)$
converges to a limit point $(\bar{\mathbf{x}}_{\varepsilon},\bar{\mathbf{Y}})$
as $l$ goes to infinity.

We use the Skorokhod representation theorem (see Theorem 6.7 in [21]) to get a
probability space $(\hat{\Omega},\hat{\mathcal{F}},\hat{\mathbb{P}})$ with a
sequence of processes
$\\{(\hat{\mathbf{X}}^{n_{l}},\hat{\mathbf{Y}}^{n_{l}}),l\geq 1\\}$ and two
processes $\hat{\mathbf{x}}_{\varepsilon}$ and $\hat{\mathbf{Y}}$ such that:

  * –

$\left(\hat{\mathbf{X}}^{n_{l}},\hat{\mathbf{Y}}^{n_{l}}\right)\rightarrow(\hat{\mathbf{x}}_{\varepsilon},\hat{\mathbf{Y}})$,
$\hat{\mathbb{P}}$-almost surely;

  * –

$\left(\hat{\mathbf{X}}^{n_{l}},\hat{\mathbf{Y}}^{n_{l}}\right)\overset{d}{\sim}(\frac{1}{n_{l}}\bar{\mathbf{X}}_{\varepsilon}(n_{l},\cdot),\bar{\mathbf{Y}}^{n_{l}})$
for all $l\geq 1$;

  * –

$(\hat{\mathbf{x}}_{\varepsilon},\hat{\mathbf{Y}})\overset{d}{\sim}(\bar{\mathbf{X}}_{\varepsilon},\bar{\mathbf{Y}})$.

The second point gives, for all $l\geq 1$:

$\hat{X}^{n_{l}}_{i,C}(t)=\hat{Y}_{i,C}^{n_{l}}(t)-\frac{1}{L_{C}}\int_{0}^{t}\hat{X}^{n_{l}}_{i,C}(s)R_{i,C}^{n}(\hat{\mathbf{X}}^{n_{l}}(s))\mathrm{d}s.$
(21)

Using the C-tightness of $\bar{\mathbf{X}}_{\varepsilon}$, $\hat{\mathbf{X}}$
and $\hat{\mathbf{Y}}$ are continuous. This implies that:

$\sup_{l\geq 1}\sup_{0\leq t\leq T}\lVert X^{n_{l}}_{i,C}(t)\rVert<\infty.$

Furthermore, we have, for all $x\in\mathbb{N}^{\varepsilon\times 2^{K}-1}$:

$\lim_{l\rightarrow\infty}R_{i,C}^{n_{l}}(x)=\frac{\lvert
C\rvert}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
x_{k,U}}\equiv R_{i,C}(x).$

Combining all these results, we can take the limits as $l$ goes to infinity in
(21) and use the dominated convergence theorem to get the following equation
for $(\hat{\mathbf{x}}_{\varepsilon},\hat{\mathbf{Y}})$:

$\displaystyle\hat{x}_{i,C}(t)=\hat{Y}_{i,C}(t)-\frac{1}{L_{C}}\int_{0}^{t}\frac{\lvert
C\rvert\hat{x}_{i,C}(s)}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap
U\rvert\hat{x}_{k,U}(s)}\mathrm{d}s.$

Let us remind that, from Lemma 8.1, $\bar{Y}_{i,C}(t)=\bar{X}_{i,C}(0)+\lambda
p_{C}\varepsilon^{2}t$, which gives us the following equation for the limit
process:

$\hat{x}_{i,C}(t)=\hat{x}_{i,C}(0)+\lambda
p_{C}\varepsilon^{2}t-\int_{0}^{t}\frac{1}{L_{C}}\frac{\lvert
C\rvert\hat{x}_{i,C}(s)}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap
U\rvert\hat{x}_{k,U}(s)}\mathrm{d}s.$ (22)

To obtain the intended equations, we differentiate (22), which concludes the
proof.

∎

## B: Proof of Theorem 4.4

In this section, we provide a proof for Theorem 4.4 using a lemma from [19]
describing the evolution of dynamics encompassing those of our fluid model.
Let us state the following lemma:

######  Lemma 8.4.

Assume that there exists
$\mathbf{z}\in\mathbb{R}_{+,\star}^{N_{\varepsilon}\times 2^{K}-1}$ such that
for all $i,C$:

$\lambda p_{C}\varepsilon^{2}\leq\frac{1}{L_{C}}\frac{\lvert C\rvert
z_{i,C}}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
z_{k,U}}$

Then, the Markov chain $\bar{\mathbf{X}}_{\varepsilon}$ is positive Harris
recurrent.

###### Proof.

To obtain this result, we use Lemma 2 from [19]: the fluid limit
$\mathbf{x}(t)$ meets the intended requirements. Let us take such a
$\mathbf{z}\in\mathbb{R}_{+,\star}^{N_{\varepsilon}\times 2^{K}-1}$. We know
that for all $0<\delta<M<\infty$, there exists $T>0$ such that whenever
$\lVert\mathbf{x}(0)\rVert=M$, we have $\lvert\mathbf{x}(T)\rVert<\delta$.

Theorem 4.2 from [14] with the fluid limit $\bar{\mathbf{x}}$ allows us to
obtain positive Harris recurrence for the chain $\mathbf{X}$. ∎

We use this result in $\mathbb{R}_{+}^{N_{\varepsilon}\times 2^{K}-1}$ by
setting $\psi_{i,C}(\mathbf{x})=\frac{1}{L_{C}}\frac{\lvert C\rvert
x_{i,C}}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
x_{k,U}}$, which is 0-homogeneous, and non-increasing in $x_{k,U}$ for
$(k,U)\neq(i,C)$.

The goal here is to find an appropriate vector $z$ to get the desired upper
bound on $\lambda$. Let us set $z_{i,C}=p_{C}L_{C}$ for all
$C\in\mathcal{P}(K)$. We want that, for each $C\in\mathcal{P}(K)$:

$\lambda p_{C}\varepsilon^{2}\leq\frac{1}{L_{C}}\frac{\lvert C\rvert
p_{C}L_{C}}{\sum_{k,U}\ell^{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
p_{U}L_{U}}.$ (23)

Using the torus property in our system, (23) becomes:

$\lambda\leq\frac{1}{\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle}\frac{\lvert
C\rvert}{\sum_{U}\lvert C\cap U\lvert p_{U}L_{U}}.$ (24)

We will use the following Lemma to simplify (24):

######  Lemma 8.5.

Let $(y_{C})_{C\in\mathcal{P}(K)}\in\mathbb{R}^{2^{K}-1}$ be such that $\lvert
C\rvert=\lvert D\rvert\Rightarrow y_{C}=y_{D}$. We have, for all
$C\in\mathcal{P}(K)$:

$\sum_{U\in\mathcal{P}(K)}\lvert C\cap U\rvert y_{U}=\frac{\lvert
C\rvert}{K}\sum_{U\in\mathcal{P}(K)}\lvert U\rvert y_{U}.$ (25)

###### Proof.

Let $C\in\mathcal{P}(K)$ with $\lvert C\rvert=j$, and
$(y_{U})_{U\in\mathcal{P}(K)}\in\mathbb{R}^{2^{K}-1}$ be such that
$y_{U}=y_{V}$ whenever $\lvert U\rvert=\lvert V\rvert$. We have:

$\sum_{U\in\mathcal{P}(K)}\lvert C\cap U\rvert
y_{U}=\sum_{l=1}^{K}\sum_{m=1}^{j\wedge l}\sum_{\begin{subarray}{c}U:\lvert
U\rvert=l\\\ \lvert C\cap U\rvert=m\end{subarray}}my_{U}.$ (26)

The number of sets $U\in\mathcal{P}(K)$ such that $\lvert U\rvert=l$ and
$\lvert C\cap U\rvert=m$ is equal to $\binom{j}{m}\binom{K-j}{l-m}$. Using the
symmetry property of $y$, we know that $y_{U}=y_{[1,l]}$. Using the fact that
we have $y_{U}=\frac{1}{\binom{K}{l}}\sum_{\lvert V\rvert=l}y_{V}(t)$. We can
rewrite (26) as:

$\displaystyle\sum_{U\in\mathcal{P}(K)}\lvert C\cap U\rvert y_{U}$
$\displaystyle=\sum_{l=1}^{K}y_{[1,l]}\sum_{m=1}^{j\wedge
l}m\binom{j}{m}\binom{K-j}{l-m}$
$\displaystyle=\sum_{l=1}^{K}\left(\sum_{\lvert
U\rvert=l}y_{U}\right)\sum_{m=1}^{j\wedge
l}m\frac{\binom{j}{m}\binom{K-j}{l-m}}{\binom{K}{l}}.$

We can prove that $\sum_{m=1}^{j\wedge
l}m\frac{\binom{j}{m}\binom{K-j}{l-m}}{\binom{K}{l}}=\frac{lj}{K}$ by using
the formula for the expectation of a hypergeometric variable with parameters
$K$, $j$ and $l$. Replacing $l$ by $\lvert U\rvert$ and $j$ by $\lvert
C\rvert$ leads to:

$\sum_{U\in\mathcal{P}(K)}\lvert C\cap U\rvert y_{U}=\frac{\lvert
C\rvert}{K}\sum_{U\in\mathcal{P}(K)}\lvert U\rvert y_{U},$

which concludes the proof of Lemma 8.5.

∎

The vector $(p_{C})_{C\in\mathcal{P}(K)}$ verifies the conditions of Lemma
8.5, which allows us to obtain, after setting $\mathfrak{L}=\sum_{U}\lvert
U\rvert p_{U}L_{U}$:

$\lambda\leq\frac{1}{\langle\ell^{\varepsilon}_{\mathcal{D}}\rangle}\frac{K}{\mathfrak{L}}.$
(27)

The result of Lemma 8.4 gives us stability for the fluid model from Theorem
8.3 in the sense of Definition 4.1 of [14] with $M=1$. Applying the result
from Theorem 4.2 of the same paper gives us positive Harris recurrence for the
chain $\bar{\mathbf{X}}_{\varepsilon}$. $\bar{\mathbf{X}}_{\varepsilon}$ is an
irreducible Markov jump process that is also positive recurrent. Thus,
$\bar{\mathbf{X}}_{\varepsilon}$ is ergodic, which concludes the proof of
Theorem 4.4.

## C: Proof of Theorem 4.5

To prove the instability of the chain $\underline{\mathbf{X}}_{\varepsilon}$,
we start by introducing, for each $i,C$, the function $r_{i,C}$ defined as:

$r_{i,C}(\mathbf{x})=\frac{1}{p_{C}L_{C}\varepsilon^{2}}\frac{\lvert C\rvert
x_{i,C}}{\sum_{k,U}\ell_{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
x_{k,U}}.$

For all $i,C$, $r_{i,C}$ is continuous and 0-homogenous, and it is not defined
for $\mathbf{x}=0$: let us take the sequences $\mathbf{x}_{n}$ and
$\mathbf{y}_{n}$ such that $x_{0,\\{1\\},n}=y_{0,\\{1\\},n}=\frac{1}{n}$,
$y_{1,\\{1\\},n}=\frac{1}{n}$ and all the other coordinates set to 0.
Immediately,
$\lim_{n\rightarrow\infty}\mathbf{x}_{n}=\lim_{n\rightarrow\infty}\mathbf{y}_{n}=0$.
Using the 0-homogeneity of the functions, we have, for all $n$:

$\displaystyle
r_{0,\\{1\\}}(\mathbf{x}_{n})=\frac{1}{p_{\\{1\\}}L_{\\{1\\}}\varepsilon^{2}},$
$\displaystyle
r_{0,\\{1\\}}(\mathbf{y}_{n})=\frac{1}{p_{\\{1\\}}L_{\\{1\\}}\varepsilon^{2}}\frac{1}{1+\ell_{\varepsilon}(a_{0},a_{1})},$

which have different limits as $n$ goes to infinity.

Let us now consider the function
$\mathbf{x}\in\mathbb{R}_{+,*}^{N_{\varepsilon}\times
2^{K}-1}\mapsto\min_{i,C}r_{i,C}(\mathbf{x})$. It is also 0-homogenous, not
defined at $\mathbf{x}=0$, and continuous on the set
$\mathcal{S}=\\{\mathbf{x}\in\mathbb{R}_{+,*}^{N_{\varepsilon}\times
2^{K}-1}:\lvert\mathbf{x}\rvert=1\\}$, which is compact. Thus, it admits a
maximum on $\mathcal{S}$, and on $\mathbb{R}_{+,*}^{N_{\varepsilon}\times
2^{K}-1}$ as a consequence. We set

$\mathcal{S}=\operatorname*{arg\,max}_{\mathbf{x}\in\mathbb{R}_{+,*}^{N_{\varepsilon}\times
2^{K}-1}}\min_{i,C}\frac{1}{p_{C}L_{C}\varepsilon^{2}}\frac{\lvert C\rvert
x_{i,C}}{\sum_{k,U}\ell_{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
x_{k,U}},$

which is non-empty.

######  Lemma 8.6.

Let $\mathbf{z}\in\mathcal{S}$ and let $i^{\star},C^{\star}$ be the
coordinates of a point where the maximum is attained. Let $\lambda>0$ be such
that:

$\lambda p_{C^{\star}}\varepsilon^{2}>\frac{1}{L_{C^{\star}}}\frac{\lvert
C^{\star}\rvert z_{i^{\star},C^{\star}}}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i}^{\star})z_{k,U}}.$

Then, $\underline{\mathbf{X}}_{\varepsilon}$ is transient.

###### Proof.

Let $\mathbf{z}\in\mathcal{S}$. By definition, there exists
$i^{\star},C^{\star}$ such that the maximum value of $\min_{i,C}r_{i,C}$ is
equal to $r_{i^{\star},C^{\star}}(\mathbf{z})$. Let us take $\lambda>0$ such
that:

$\lambda p_{C^{\star}}\varepsilon^{2}>\frac{1}{L_{C^{\star}}}\frac{\lvert
C^{\star}\rvert z_{i^{\star},C^{\star}}}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i}^{\star})z_{k,U}}.$

We define the process $\underline{\mathbf{Y}}$ such that:

  * –

Arrivals in queue $i,C$ happen with rate $\lambda p_{C}\varepsilon^{2}$, for
all $i,C$;

  * –

Departures in queue $i,C$ happen with rate:

    * –

$\frac{1}{L_{C^{\star}}}\frac{\lvert C^{\star}\rvert
z_{i^{\star},C^{\star}}}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i})z_{k,U}}$ if $i,C=i^{\star},C^{\star}$,

    * –

$\frac{1}{L_{C}}\frac{\lvert C\rvert\underline{Y}_{i,C}(t)}{\sum_{k,U}\lvert
C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i})\underline{Y}_{k,U}(t)+\mathcal{N}_{0}}$
else.

$\underline{\mathbf{Y}}$ is a Markov jump process, with state space
$\mathbb{N}^{N_{\varepsilon}\times 2^{K}-1}$. By definition of $\mathbf{z}$,
we have, at all times $t\geq 0$:

$\frac{1}{L_{C^{\star}}}\frac{\lvert C^{\star}\rvert
z_{i^{\star},C^{\star}}}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i})z_{k,U}}\geq\frac{1}{L_{C^{\star}}}\frac{\lvert
C^{\star}\rvert\underline{X}_{i^{\star},C^{\star}}(t)}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i})\underline{X}_{k,U}(t)+\mathcal{N}_{0}}.$

Using a coupling argument similar to the one used in the proof of Theorem 3.1,
we obtain, for all $i,C$ and at all times $t$:

$\underline{Y}_{i,C}(t)\leq\underline{X}_{i,C}(t),\quad\mathbb{P}\textrm{-a.s.}$

The queue $\underline{Y}_{i^{\star},C^{\star}}$ is an M/M/1 queue with
constant arrival rate $\lambda p_{C^{\star}}\varepsilon^{2}$ and a departure
rate equal to $\frac{1}{L_{C^{\star}}}\frac{\lvert C^{\star}\rvert
z_{i^{\star},C^{\star}}}{\sum_{k,U}\lvert C\cap
U\rvert\ell_{\varepsilon}(a_{k},a_{i})z_{k,U}}$.

It is unstable and
$\mathbb{P}\left[\lim_{t\rightarrow\infty}\underline{Y}_{i^{\star},C^{\star}}(t)=+\infty\right]=1$,
implying the same for $\underline{X}_{i^{\star},C^{\star}}$, which concludes
the proof.

∎

We can note that queue $\underline{X}_{i^{\star},C^{\star}}$ is not the only
queue whose population diverges to infinity. The departure rates in all queues
(except for queue $i^{\star},C^{\star}$) are decreasing functions of
$\underline{X}_{i^{\star},C^{\star}}(t)$. Thus, for a given queue
$\underline{X}_{j,D}$, there exists a time $T_{j,D}$ after which the departure
rate in this queue becomes lower than its arrival rate. We can then bound from
below $\underline{X}_{j,D}$ after time $T_{j,D}$ by an adequate M/M/1 queue
and obtain that the population in queue $j,D$ goes to infinity,
$\mathbb{P}$-almost surely for all $j,D$.

To complete the proof of Theorem 4.5, we characterize the value of
$\mathbf{z}$:

######  Lemma 8.7.

Let $\mathbf{z}\in\mathcal{S}$. Then, $\mathbf{z}$ is a solution to the
following system of equations, for all $0\leq i,j\leq N_{\varepsilon}-1$ and
$C,D\in\mathcal{P}(K)$:

$\frac{1}{p_{C}L_{C}}\frac{\lvert C\rvert
z_{i,C}}{\sum_{k,U}\ell_{\varepsilon}(a_{k},a_{i})\lvert C\cap U\rvert
z_{k,U}}=\frac{1}{p_{D}L_{D}}\frac{\lvert D\rvert
z_{j,D}}{\sum_{k,U}\ell_{\varepsilon}(a_{k},a_{j})\lvert D\cap U\rvert
z_{k,U}}.$ (28)

###### Proof.

Let $\mathbf{z}$ be such that:

$\mathbf{z}\in\operatorname*{arg\,max}_{\mathbf{y}\in\mathbb{R}_{+}^{N_{\varepsilon}\times
2^{K}-1}}\min_{i,C}r_{i,C}(\mathbf{y}).$

To prove this result, we use the maximality of $\mathbf{z}$. Without loss of
generality, let us assume that the minimum over $i,C$ of $r_{i,C}(\mathbf{z})$
is reached for $i=0$ and $C=\\{1\\}$. We know that the function
$r_{0,\\{1\\}}$ is a decreasing function of $z_{i,C}$ for all $1\leq i\leq
N-1$ and $C\neq\\{0\\}$. Hence, these values have to be minimal in order to
maximize $r_{0,\\{1\\}}$.

By definition, for all $1\leq j\leq N_{\varepsilon}-1$, $C\neq\\{1\\}$, we
have $r_{0,\\{1\\}}(\mathbf{z})\leq r_{i,C}(\mathbf{z})$. But, for $i\geq 1$
and $C\neq\\{0\\}$, $r_{i,C}$ is an increasing function of $z_{i,C}$. Thus,
the only possible value for $\mathbf{z}$ is such that:

$r_{0,\\{1\\}}(\mathbf{z})=r_{i,C}(\mathbf{z}),$

which concludes the proof.

∎

Finally, we know that $\mathbf{z}$ such that for all $i,C$,
$z_{i,C}=p_{C}L_{C}$ is a solution to (28). Using Lemma 8.6, the system is
unstable if:

$\displaystyle\lambda p_{C^{\star}}\varepsilon^{2}\geq$
$\displaystyle\frac{1}{L_{C^{\star}}}\frac{\lvert C^{\star}\rvert
p_{C^{\star}}L_{C^{\star}}}{\sum_{k,U}\ell_{\varepsilon}(a_{k},a_{i^{\star}})\lvert
C^{\star}\cap U\rvert p_{U}L_{U}}$ $\displaystyle\overset{(a)}{=}$
$\displaystyle\frac{\lvert C^{\star}\rvert
p_{C^{\star}}}{\frac{1}{\varepsilon^{2}}\langle\ell_{\varepsilon,\mathcal{D}}\rangle\frac{\lvert
C^{\star}\rvert}{K}\sum_{U}\lvert U\rvert p_{U}L_{U}}$ $\displaystyle=$
$\displaystyle\frac{\varepsilon^{2}p_{C^{\star}}K}{\langle\ell_{\varepsilon,\mathcal{D}}\rangle\frac{\lvert
C^{\star}\rvert}{K}\sum_{U}\lvert U\rvert p_{U}L_{U}},$

where $(a)$ uses Lemma 8.5 with the vector $(p_{C}L_{C})_{C\in\mathcal{P}(K)}$
and the square torus property of $\mathcal{D}$. Setting
$\mathfrak{L}=\sum_{U}\lvert U\rvert p_{U}L_{U}$ and rearranging the equation
gives us:

$\lambda\geq\frac{K}{\langle\ell_{\varepsilon,\mathcal{D}}\rangle\mathfrak{L}},$

which is the intended condition for $\lambda$ and concludes the proof.