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

Pierre Popineau, François Baccelli

<figure><img src="image_1.png" /></figure>

► To cite this version:

Pierre Popineau, François Baccelli. On multiclass spatial birth-and-death processes with wireless-type interactions. IEEE Transactions on Information Theory, 2024, pp.1-1. 10.1109/tit.2024.3392979. hal-04884946

HAL Id: hal-04884946

https://hal.science/hal-04884946v1

Submitted on 3 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On multiclass spatial birth-and-death processes with
wireless-type interactions

Pierre Popineau¹, François Baccelli¹,²

Affiliation: ¹ INRIA-ENS, Paris, France ² 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 \ge 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\le i\le N_t}$ denotes the location of receivers and $(y_i)_{1\le i\le N_t}$ that of transmitters. Let $\Phi_t^T$ be the point process describing the location of transmitters in the network and $\Phi_t^R$ 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 $|C| = |D|$, 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 \le j \le 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 \le i \le d-1} \in \mathbb{R}^d$), and $\le$ 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}$, $|x|$ be the $L^1$ norm and $\|x\|$ 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(\|x\|)dx, & \forall f: \mathbb{R} \to \mathbb{R} \\ \int_{\mathcal{D}} f(x, 0)dx, & \forall f: \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}, \end{cases} \quad (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 \ge 0}$ and $x$ in $D([0, \infty), \mathbb{R}^{N \times 2^K - 1})$, we write $x^n \to x$ if, for all $T > 0$:

$$ \lim_{n \to \infty} \sup_{0 \le t \le T} \|x^n(t) - x(t)\| = 0. $$

We denote by $\mathcal{E}(\mathbb{R}^d)$ the set of probabilities on $\mathbb{R}^d$ (with $d \ge 1$), and we define the partial ordering of probability measures $\le_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 \le_i G$ if and only if

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

for all $f \in I$. $\le_i$ is a partial order on the set of random vectors of $\mathbb{R}^d$. Let $\Psi$ and $\Psi'$ be two point processes on $\mathcal{D}$. We write $\Psi \le \Psi'$ if we have $\Psi(\mathcal{D}) \le_i \Psi'(\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 \stackrel{d}{\sim} Y$ if $X$ and $Y$ have the same probability distribution. Table 2.3 summarizes the notation used throughout the paper.

## 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 \setminus \{y\}} |C_x \cap C_z| \ell(\|x - z\|), \quad (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
<table><caption>Table 1: Table of notations</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td>$\Phi_t^R, \Phi_t^T$</td><td>Receiver (resp. transmitter) locations at time $t$</td></tr><tr><td>$\Phi_{C,t}$</td><td>Point process of users of class $C$ at time $t$</td></tr><tr><td>$\mathcal{P}(K)$</td><td>Set of subsets of $[1, \dots, K]$ with the exception of $\emptyset$</td></tr><tr><td>$\ell$</td><td>Path-loss function</td></tr><tr><td>$r$</td><td>Receiver-transmitter distance</td></tr><tr><td>$\lambda$</td><td>Intensity of the arrival process</td></tr><tr><td>$(p_C)_{C \in \mathcal{P}(K)}$</td><td>Arrival distribution of users of each type</td></tr><tr><td>$L_C$</td><td>Average file size for users of class $C$</td></tr><tr><td>$R$</td><td>Rate-of-transmission function</td></tr><tr><td>$\mathcal{N}_0$</td><td>Thermal noise density in the network</td></tr><tr><td>$\langle \ell_D \rangle$</td><td>$\int_{\mathbf{x} \in \mathcal{D}} \ell(\|\mathbf{x}\|) d\mathbf{x}$</td></tr><tr><td>$\Phi_{0,C}$</td><td>Stationary point process of users of type $C$</td></tr><tr><td>$\mu_C$</td><td>Spatial intensity of point process $\Phi_{0,C}$</td></tr><tr><td>$\mathbb{E}_{\Phi_{0,C}}^0$</td><td>Palm expectation with respect to $\Phi_{0,C}$</td></tr></tbody></table>

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|C_x|\ell(r)}{\mathcal{N}_0 + I(x, \Phi_t)}, \qquad (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: $|C_x|\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 $|C_x|\ell(r) \ll \mathcal{N}_0 + I(x, \Phi_t)$, the Shannon-Hartley (see [1]) formula gives a channel capacity equal to:

$$C(x) = \log_2 \left( 1 + \frac{|C_x|\ell(r)}{\mathcal{N}_0 + I(x, \Phi_t)} \right) \approx \frac{1}{\ln(2)} \frac{|C_x|\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) du > L_{x,y} \right\}. \quad (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). \quad (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 \le_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}] \ge \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 \le 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:

$$\left\{ \begin{array}{l} \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{array} \right. \quad (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 \ge \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 \le \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, $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' = \Delta(\Phi'_0, \lambda', d')$. The following conditions are sufficient for $\Phi'$ to dominate $\Phi$ (with all the other parameters taken equal):

i) $\lambda \le \lambda'$;

ii) $L \le L'$;

iii) for all point processes $\Psi \le \Psi'$ on $\mathcal{D}$ and $x \in \mathcal{D}$, $R'(x, \Psi') \le R(x, \Psi)$;

iv) $\Phi_0 \subseteq \Phi'_0$.

It is to note that condition iii) is met when we have two path-loss functions $\ell$ and $\ell'$ such that for all $r \ge 0$, $\ell(r) \le \ell'(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'$ with the same fixed initial condition $\Phi_0$, a vector $L$ and a rate-of-transmission function $R$, and let $\Phi_t = \Delta(\Phi_0, \lambda, d)$ and $\Phi'_t = \Delta(\Phi_0, \lambda', d)$. The quantities related to $\Phi'$ will be denoted with a prime. We want to prove that for all times $0 \le t$, $\Phi_t \subseteq \Phi'_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'$ is $\mathcal{A} \cup \mathcal{A}'$, with $\mathcal{A}'$ being a Poisson rain with intensity $\lambda' - \lambda > 0$ independent from $\mathcal{A}$, so that common arrivals in $\Phi_t$ and $\Phi'_t$ happen at the same locations and times. Using this coupling, we will show that, at all times $t$, $\Phi_t \subseteq \Phi'_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'_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'_t = \sum_C \sum_{i=1}^{\Phi'_{C,t}(\mathcal{D})} \delta_{x'_{C,i}}$.

Assume that up to time $t$, $\Phi_t \subseteq \Phi'_t$, the next event happening at a time $\hat{t} \ge 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'_t$;
- The departure of a user of class C in $\Phi_t$;

- The departure of a user of class C in $\Phi'_t$.

Because arrivals are coupled, an arrival in either $\Phi_t$ or $\Phi'_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'_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'_{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)$:

$$ R(x, \Phi_{\hat{t}^-}) = \frac{|C_x|\ell(r)}{\mathcal{N}_0 + \sum_{y \in \Phi_{\hat{t}^-}} |C_x \cap C_y|\ell(\|x - y\|)} \\ \stackrel{(a)}{\ge} \frac{|C_x|\ell(r)}{\mathcal{N}_0 + \sum_{y \in \Phi'_{\hat{t}^-}} |C_x \cap C_y|\ell(\|x - y\|)} \\ = R(x, \Phi'_{\hat{t}^-}), $$

where (a) uses the fact that $|\Phi_t(\mathcal{D})| \le |\Phi'_t(\mathcal{D})|$. In other words, the departure rates in $\Phi_{\hat{t}^-}$ are larger than in $\Phi'_{\hat{t}^-}$.

We take the Poisson imbedding of the departure processes (see [20]): let $\mathcal{D}_{C,t}$ and $\mathcal{D}'_{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'_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'_{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'_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. $\square$

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 \le i \le N}$, $x_i \in \Phi_t^R$ and $y_i \in \Phi_t^T$, and let $(C_i)_{1 \le i \le 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 } |\Phi_t] = \prod_{i=1}^{N} \frac{d_{C_i}}{d_{C_i} + b_{C_i}}, $$

where $b_{C_i} = \lambda p_{C_i} |\mathcal{D}|$ 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(\emptyset, A_n) = \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} H(dy_1) \dots H(dy_n) e^{-\lambda p_{C_1} t_1} \dots e^{-\lambda p_{C_n} (t_n - t_{n-1})} $$
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$. $\square$

## 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 \le \lambda_c^- \le \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:

$$ \lambda_c^- = \sup_{\lambda>0}\{\lambda \text{ such that } \Phi_t \text{ is stable}\}, $$

$$ \lambda_c^+ = \inf_{\lambda>0}\{\lambda \text{ such that } \Phi_t \text{ is unstable}\}. $$

$\Phi_t$ is an irreducible Markov chain on $M(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 \le \lambda_c^-$ and transient for all $\lambda \ge \lambda_c^+$, which gives $\lambda_c^- \le \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^- \ge 0$, which concludes the proof. $\square$

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_D\rangle} \le \lambda_c^- \le \lambda_c^+ \le \frac{K\ell(r)}{\underline{L}\langle\ell_D\rangle}, $$

where $\bar{L} = \max_C L_C$, $\underline{L} = \min_C L_C$ and $\langle\ell_D\rangle = \int_{x\in D}\ell(\|x\|)dx$.
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, ..., 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 \setminus \{x\}} K \ell(\|x - y\|)}.$$

Moreover, by definition, $\mathbf{L} < \bar{\mathbf{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 \setminus \{x\}} \ell(\|x - y\|)},$$

and $\underline{\mathbf{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{\mathbf{L}}, R_u)$ and $\Phi_{d,t} = \Delta(\Phi_0, \lambda, \bar{\mathbf{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{\mathbf{L}}\langle\ell_D\rangle}$ for $\Phi_{d,t}$ and to $\frac{K\ell(r)}{\underline{\mathbf{L}}\langle\ell_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{\mathbf{L}}\langle\ell_D\rangle} > 0$ and $\frac{K\ell(r)}{\underline{\mathbf{L}}\langle\ell_D\rangle} < \infty$, we can conclude that $0 < \lambda_c^-$ and $\lambda_c^+ < \infty$. $\square$

# 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_D \rangle \mathcal{L}}, \qquad (7)$$

where $\mathcal{L} \triangleq \sum_C p_C |C| 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 tessel- late $\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\le i\le 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|\mathcal{D}|$ 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|\mathcal{D}|$ 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) & \text{if } x \in A_i, y \in A_j \\ \ell^\varepsilon(a_i, a_j) = \max \{\ell(\|b_i - b_j\|), b_j \in \mathcal{V}_j, b_i \in \mathcal{V}_i\}, \end{cases} $$

and

$$ \begin{cases} \ell_\varepsilon(x, y) = \ell_\varepsilon(a_i, a_j) & \text{if } x \in A_i, y \in A_j \\ \ell_\varepsilon(a_i, a_j) = \min \{\ell(\|b_i - b_j\|), b_j \in \mathcal{V}_j, b_i \in \mathcal{V}_i\}, \end{cases} $$

where, for $0 \le i \le N_\varepsilon - 1$, we define $\mathcal{V}_i = \{b_i, \|b_i - a_i\| \in \{0, \varepsilon\}\}$. Using the dominated convergence theorem, we get:

$$ \lim_{\varepsilon \to 0^+} \langle \ell_{\mathcal{D}}^{\varepsilon} \rangle = \lim_{\varepsilon \to 0^+} \langle \ell_{\varepsilon, \mathcal{D}} \rangle = \langle \ell_{\mathcal{D}} \rangle, \quad (8) $$

where $\langle \cdot, \cdot \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_{\mathcal{D}}^{\varepsilon} \rangle. \qquad (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:

$$ \bar{I}_{i,C}(t) = \sum_{k,U} |C \cap U| \ell^{\varepsilon}(a_k, a_i) (\bar{X}_{k,U}(t) - \mathbf{1}_{\{U=C, i=k\}}) $$

$$ \underline{I}_{i,C}(t) = \sum_{k,U} |C \cap U| \ell^{\varepsilon}(a_k, a_i) (\underline{X}_{k,U}(t) - \mathbf{1}_{\{U=C, i=k\}}). $$

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

$$ \bar{R}_{i,C}(t) = \frac{|C|}{\mathcal{N}_0 + \bar{I}_{i,C}(t)} $$

$$ \underline{R}_{i,C}(t) = \frac{|C|}{\mathcal{N}_0 + \underline{I}_{i,C}(t)}. $$

**Theorem 4.2.** $\bar{X}_{\varepsilon}$ and $\underline{X}_{\varepsilon}$ are two irreducible Markov jump processes with state space $\mathbb{N}^{N_{\varepsilon} \times 2^K - 1}$. Moreover, $\bar{X}_{\varepsilon}$ stochastically dominates $X_{\varepsilon}$, and $\underline{X}_{\varepsilon}$ is stochastically dominated by $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) \le \ell(\|x - y\|) \le \ell^{\varepsilon}(x, y). $$

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

$$ \bar{R}_{i,C}(t) \le R(x, \Phi_t) \le \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 $X_{\varepsilon}$ and $\bar{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{X}_\epsilon$, which will allow us to obtain a condition on $\lambda_c^-$ for the stability of the system.

Let $(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 \epsilon^2$. If the system is in state $x = (x_{i,C}) \in \mathbb{N}^{N_\epsilon \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{X}_\epsilon$:

$$ \bar{X}_{i,C}(t) = \bar{X}_{i,C}(0) + A_{i,C}(\lambda p_C \epsilon^2 t) - N_{i,C} \left( \frac{1}{L_C} \int_0^t \bar{X}_{i,C}(u) \bar{R}_{i,C}(u) du \right). \quad (10) $$

The goal of this section is to find a condition on $\lambda$ so that the Markov chain $\bar{X}_\epsilon(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{X}_\epsilon$ exist and are solutions to the following system of equations:

$$ \begin{cases} \bar{x}'_{i,C}(t) = \lambda p_C \epsilon^2 - \frac{1}{L_C} \frac{|C| \bar{x}_{i,C}(t)}{\sum_{k,U} |C \cap U| \ell^\epsilon(a_k, a_i) \bar{x}_{k,U}} & \text{if } \bar{x}(t) \neq 0 \\ \bar{x}_0 = \bar{X}_\epsilon(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{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{X}_\epsilon$, which gives us a lower bound for $\lambda_c^-$. We have the following theorem:

**Theorem 4.4.** Let $\epsilon > 0$ and $\bar{\lambda}_\epsilon = \frac{K}{\mathcal{L}\langle\ell_D^\epsilon\rangle}$. If $\lambda < \bar{\lambda}_\epsilon$, then the chain $\bar{X}_\epsilon$ 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{X}_\varepsilon$ dominates the original dynamic. From Theorem 3.1, we obtain that, for all $\varepsilon > 0$, $\lambda_c \ge \frac{K}{\mathcal{L}\langle\ell_D\rangle}$. Taking the limit as $\varepsilon$ goes to 0 gives us:

$$ \lambda_c^- \ge \frac{K}{\mathcal{L}\langle\ell_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{x} = 0$. To obtain instability for the system, we will use stochastic domination and an adequate Markov chain to bound $\underline{X}$ from below. We have the following theorem:

**Theorem 4.5.** Let $\underline{\lambda}_\varepsilon = \frac{K}{\langle\ell_{\varepsilon,D}\rangle\mathcal{L}}$. In the symmetric case, if $\lambda > \underline{\lambda}_\varepsilon$, then $\underline{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,D}\rangle\mathcal{L}}$ for each $\varepsilon > 0$. Taking the limit as $\varepsilon$ goes to 0, we conclude:

$$ \lambda_c^+ \le \frac{K}{\langle\ell_D\rangle\mathcal{L}}. $$

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

$$ \lambda_c^- = \lambda_c^+ = \frac{K}{\langle\ell_D\rangle\mathcal{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{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{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{S}(t) = (\bar{X}_\varepsilon(t), \bar{Z}_\varepsilon(t), \bar{M}_\varepsilon(t))$ is a Markov chain with countable state space. The evolution of $\bar{S}(t)$ is the following: to each receiver of class $C$ arriving in cell
i with rate $\lambda p_C \epsilon^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{Z}_\epsilon(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{S}(t)$ and the same steps as for Theorem 8.3, we can obtain fluid limits for the process $\bar{X}_\epsilon(t)$ in this setup (which is symmetric):

$$ \bar{x}'_{i,C}(t) = \lambda p_C \epsilon^2 - \frac{1}{L_C} \frac{|C| \ell(r) \bar{x}_{i,C}(t)}{\sum_{k,U} \ell^\epsilon(a_k, a_i) |C \cap U| \bar{z}_{k,U}(t)}. $$

For all $i, C$ and $t \ge 0$, let $i_C^*(t)$ be such that $z_{i_C^*(t),C} = \max_i z_{i,C}(t)$. The process $(\bar{z}_{i_C^*(t),C})$ dominates $\bar{Z}_\epsilon$, and we have:

$$ \bar{z}'_{i,C}(t) \le \lambda p_C \epsilon^2 - \frac{1}{L_C} \frac{|C| \ell(r) \bar{z}_{i,C}(t)}{\sum_{k,U} \ell^\epsilon(a_k, a_i) |C \cap U| \bar{z}_{k_U^*(t),U}(t)}. $$

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

$$ \lambda_c \le \frac{K \ell(r)}{\mathcal{L} \langle \ell_D^\epsilon \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)}{\mathcal{L} \langle \ell_{\epsilon,D} \rangle}$. This way, we get the reciprocal inequality and we can conclude that:

$$ \lambda_c = \frac{K \ell(r)}{\mathcal{L} \langle \ell_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}{|\mathcal{D}|} \mathbb{E}[\Phi_{C,0}(\mathcal{D})]$ 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 |C| = |D| \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, |C \cap U|)} dz = \frac{\lambda p_C L_C}{|C| \ell(r)}. \quad (11) $$

$$ \text{where } \mathcal{I}(z, k) = \int_{\mathcal{D}} (1 - e^{-zk\ell(\|x\|)}) dx. $$

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 |D| dt$ users arrive on average in the network, and in a given network configuration $\Phi_{0,C}$, $\mathbb{E}[\frac{1}{L_C} \sum_{x \in \Phi_{C,0}} R(x, \Phi_0) dt]$ users leave, which gives:

$$ \lambda p_C L_C |\mathcal{D}| = \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 |\mathcal{D}| = \mathbb{E}_{\Phi_{C,0}}^0 [R(0^C, \Phi_0)] \mathbb{E}[\Phi_0^C(\mathcal{D})]. \quad (12) $$
Using this, we get that, in the stationary regime:

$$ 
\begin{aligned}
\lambda p_C L_C &= \mu_C \mathbb{E}_{\Phi_{C,0}}^0 [R(0^C, \Phi_0)] \\
&= |C| \mu_C \ell(r) \mathbb{E}_{\Phi_{C,0}}^0 \left[ \frac{1}{N_0 + I_C} \right].
\end{aligned}
$$

This is equivalent to:

$$ \mu_C = \frac{\lambda p_C L_C}{|C| \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}] dz. $$

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] dz. $$

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] \approx \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 \setminus \{0^C\}} |C \cap U| \ell(\|x\|)} \right] dz. $$

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 \setminus \{0^C\}} |C \cap U| \ell(\|x\|)} \right] = \exp \left[ -\mu_U^f \int_D (1 - e^{-z|C \cap U| \ell(\|x\|)}) dx \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 I(z, |C \cap U|)} dz = \frac{\lambda p_C L_C}{|C| \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}{|D|} \mathbb{E} \left[ \sum_{U:|U|=j} \Phi_{U,0}(D) \right] = \sum_{U:|U|=j} \mu_U$ as the intensity of users communicating on $j$ channels, for $1 \le j \le K$, with $L_j = L_C$ and $p_j = \binom{K}{j} p_C$ for all $C$ such that $|C| = j$, we obtain the following system of equations from Heuristic 2:
$$ \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)} dz = \frac{\lambda p_j L_j}{j\ell(r)}, \quad (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 its stationary regime, i.e., if $\lambda < \lambda_c$. We define the following estimate for $\lambda_c$:

**Heuristic 2.** Let $\lambda_P$ be the largest value of $\lambda$ such that (11) admits a solution (or (13) in the symmetric network). Then, $\lambda_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_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_C^s$ defined as:

$$ \mu_C^s = \frac{\lambda p_C L_C}{|C|\ell(R)} (\mathcal{N}_0 + I_C), $$
where $I_C \triangleq \mathbb{E}_{\Phi_{C,0}}^0 [I(0^C, \Phi_0)]$ 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{|C \cap U|}{\mu_C^s} \int_{x \in \mathcal{D}} \ell(\|x\|) \rho_{C,U}^{(2)}(x, 0) dx,$$

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

$$\rho_{C,D}^{(2)}(x, y) = \frac{\lambda(\mu_D^s p_C + \mu_C^s p_D)}{d_{C,D}(x, y)}, \qquad (14)$$

and:

$$d_{C,D}(x, y) = \frac{1}{L_C} \frac{|C|\ell(R)}{\mathcal{N}_0 + |C \cap D|\ell(\|x - y\|) + I_C} + \frac{1}{L_D} \frac{|D|\ell(R)}{\mathcal{N}_0 + |C \cap D|\ell(\|x - y\|) + 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) = |C \cap D|\ell(\|x - y\|) + I_C, \qquad (15)$$

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

$$\begin{aligned} I_C &= \mathbb{E}_{\Phi_0}^0 [I(0^C, \Phi_0)] \\ &= \sum_U \mathbb{E}_{\Phi_0}^0 [I(0^C, \Phi_0^U)] \\ &= \frac{1}{\mu_C} \sum_U |C \cap U| \int_{x \in \mathcal{D}} \ell(\|x\|) \rho_{C,U}^{(2)}(x, 0) dx. \end{aligned}$$

We obtain the departure rate for the pair $(x, y)$ with $x$ being of type C and $y$ being of type D:
$$d_{C,D}(x, y) = \frac{1}{L_C} \frac{|C|\ell(R)}{\mathcal{N}_0 + |C \cap D|\ell(\|x - y\|) + I_C} + \frac{1}{L_D} \frac{|D|\ell(R)}{\mathcal{N}_0 + |C \cap D|\ell(\|x - y\|) + I_D}.$$

In the stationary regime, on average, the umber of arrivals compensate the depar- tures, which gives us:

$$\lambda(\mu_D p_C + \mu_C p_D) = \rho_{C,D}^{(2)}(x, y)d_{C,D}(x, y).$$

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

$$\lambda p_C |\mathcal{D}| = \mathbb{E}[\Phi_{C,0}(\mathcal{D})] \frac{|C|\ell(r)}{\mathcal{N}_0 + I_C},$$

which gives the intended result.

<figure><img src="image_3.png" /></figure>

# 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 $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|\mathcal{D}|$ and a random vector $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 $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><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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. \qquad (16)$$

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

$$\mathbb{E}[\Phi_{t,C}(\mathcal{D})] = \lambda p_C W_C, \qquad (17)$$

where $W_C$ is the staying time of a user of class C and $\mathbb{E}[\Phi_{t,C}(\mathcal{D})]$ 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><img src="image_5.png" /><figcaption>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.</figcaption></figure>

## 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})}{|\mathcal{D}|} dw, \qquad (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).

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_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
<figure><img src="image_6.png" /><figcaption>Figure 3: Value of $\mu_{\{1\}}$ as a function of $\lambda/\lambda_c$ in a symmetric configuration.</figcaption></figure>

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.

# 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.
<figure><img src="image_7.png" /><figcaption>Figure 4: Comparison between $\lambda_c$ and its Poisson heuristic as a function of $p_{\{1,2\}}$ in the network in the symmetric case.</figcaption></figure>

## 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 $\mathcal{L} = \sum_{U} |U| 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 $|U|$ bands. Hence, a user of type $U$ brings a time-space load equal to $|U|L_U$ to the network. Taking the average over the distribution of user types, we obtain $\mathcal{L}$. We can rewrite the stability condition $\lambda < \frac{K\ell(R)}{\langle\ell_D\rangle\mathcal{L}}$ as:

$$ \lambda \frac{\mathcal{L}}{K} < \frac{\ell(R)}{\langle\ell_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 will 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 2d 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}_{i,C}^n)$ for our system, with $\lim_{n\to\infty} \bar{x}_{i,C}^n = \infty$. The goal in this section is to study the limit of the sequence of fluid-scaled processes $(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(nt))$ 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$:

$$ \frac{1}{n}\bar{X}_{i,C}(nt) = \frac{\bar{x}_{i,C}^n}{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))du\right). \quad (19) $$

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

$$ \int_0^{nt} \bar{X}_{i,C}(u) R_{i,C}(\bar{\mathbf{X}}_{\varepsilon}(u)) du = \int_0^t \frac{n |C|^{\frac{1}{n}} \bar{X}_{i,C}(ns) ds}{\frac{N_0}{n} + \sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| \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{|C|}{\frac{N_0}{n} + \sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| \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}_{i,C}^n + \frac{1}{n}\mathcal{A}_{i,C}(\lambda p_C \varepsilon^2 nt) - \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)ds\right), $$

Let $M_{i,C}^n(z) = \frac{1}{n}N_{i,C}(nz) - z$ for $z \in \mathbb{R}$, and let us define:
$$ \bar{Y}_{i,C}^n(t) = \frac{1}{n}\bar{x}_{i,C}^n + \frac{1}{n}\mathcal{A}_{i,C}(\lambda p_C \varepsilon^2 nt) - M_{i,C}^n \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) ds \right). $$

From (19), we have:

$$ \frac{1}{n}\bar{X}_{i,C}(nt) = \bar{Y}_{i,C}^n(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) ds. $$

Let us denote $\bar{\mathbf{x}}_0 = \lim_{n \to \infty} \frac{1}{n}\bar{x}_{i,C}^n$ 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 \ge 0, \forall i, C, \quad \lim_{n \to \infty} \bar{Y}_{i,C}^n(t) = \bar{X}_{i,C}(0) + \lambda p_C \varepsilon^2 t. \quad (20) $$

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

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

Moreover, we have:

$$ \begin{aligned} & \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) ds \\ &= \frac{1}{L_C} \int_0^t \frac{|C| \frac{1}{n} \bar{X}_{i,C}(ns)}{\frac{N_0}{n} + \sum_{k,U} |C \cap U| \ell^{\varepsilon}(a_k, a_i) \frac{1}{n} \bar{X}_{k,U}(ns)} ds \\ &\le \frac{1}{L_C} \int_0^t \frac{|C| \frac{1}{n} \bar{X}_{i,C}(ns)}{|C| \ell^{\varepsilon}(a_i, a_i) \frac{1}{n} \bar{X}_{i,C}(ns)} ds \\ &= \frac{t}{L_C}, \end{aligned} $$

where we use the fact that for all $0 \le i \le 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 $\{\frac{1}{n}\bar{X}_{i,C}(n\cdot) R_{i,C}^n(\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot))\}$ is tight. Finally, because $M_{i,C}^n(z) \to 0$ as $n$ goes to infinity for all $z$ (using the strong law of large numbers), we can conclude that:

$$ \lim_{n \to \infty} M_{i,C}^n \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) ds \right) = 0, $$

which leads to the intended result. $\square$
For $v: [0, T) \to \mathbb{R}$ and $0 \le t_1 \le t_2 \le T$, we define:

$$w(v, [t_1, t_2]) = \sup \{|v(u_1) - v(u_2)|, u_1, u_2 \in [t_1, t_2]\}.$$

Let us define the modulus of continuity $\omega(v, \delta, T) = \sup\{w(v, [s, t]), 0 \le s, t \le T, |s - t| < \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 \ge 1\}$ of functions is C-tight if and only if, for all $t_0 > 0$, $\eta > 0$, $T > 0$, there exists $K_\eta^0$, $n_\eta^0$ and $\delta_\eta^0$ such that, for all $n \ge n_\eta^0$:

$$i) \mathbb{P}[\sup_{0 \le t \le T} |V^n(t)| \ge K_\eta^0] < \eta.$$

$$ii) \mathbb{P}[\omega(V^n, \delta_\eta^0, t_0) > \eta] < \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}[\omega(V, \delta_\eta^0, t_0) > \eta] < \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 \ge 0\}$ is C-tight.

*Proof.* Let $T > 0$ and $\eta > 0$. As $n$ goes to infinity, we know that $\bar{Y}_{i,C}^n(t) \to \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_\eta^{T,1}$ such that for $n \ge n_\eta^{T,1}$:

$$\mathbb{P}\left[\sup_{0 \le t \le T} |\bar{Y}_{i,C}^n(t)| \ge \lambda p_C \varepsilon^2 T + 1\right] \le \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)) \le 1$ at all times $t$, a fortiori for all $t \le T$. We deduce that for all $n \ge n_\eta^{T,1}$:

$$\mathbb{P}\left[\sup_{0 \le t \le T} \left|\frac{1}{n}\bar{X}_{i,C}(nt)\right| \ge \lambda p_C \varepsilon^2 T + 1 + 1\right] \le \eta.$$

We set $K_\eta^0 = \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_\eta^0 > 0$ such that for all $i, C$ and at all times $t < T$, we have:

$$\mathbb{P}\left[\omega\left(\frac{1}{n}\bar{X}_{i,C}(nt), \delta_{\eta}^{0}, t_{0}\right) > \eta\right] < \eta.$$
Let $\delta > 0$. We have at all times $t \ge 0$:

$$ \frac{1}{n}\bar{X}_{i,C}(n(t+\delta)) - \frac{1}{n}\bar{X}_{i,C}(nt) = \bar{Y}_{i,C}^n(t+\delta) - \bar{Y}_{i,C}^n(t) - \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) du. $$

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)) \le 1$ yields:

$$ \sup_{0 \le t \le T} \left| \frac{1}{n} \bar{X}_{i,C}(n(t+\delta)) - \frac{1}{n} \bar{X}_{i,C}(t) \right| \le \sup_{0 \le t \le T} \left| \bar{Y}_{i,C}^n(t+\delta) - \bar{Y}_{i,C}^n(t) \right| + \frac{\delta}{L_C}. $$

From Lemma 8.1, we know that $\bar{Y}_{i,C}^n(t) \to \bar{X}_{i,C}(0) + \lambda p_C \varepsilon^2 t$ P-almost surely. This implies that there exists $n_{\eta}^{T,2} > 0$ and $\kappa_{\eta}^0 > 0$ such that, with probability at least $1 - \eta$, $|\bar{Y}_{i,C}^n(t+\delta) - \bar{Y}_{i,C}^n(t)| < \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_{\eta}^0 = \eta\left(\lambda p_C \varepsilon^2 + \kappa_{\eta}^0 + \frac{1}{L_C}\right)^{-1}$. Thus, we get, for all $n \ge n_{\eta}^{T,2}$:

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

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

<figure><img src="image_8.png" /></figure>

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) \to \bar{\mathbf{x}}_0$ as $n$ goes to $\infty$, then the sequence of processes $\frac{1}{n}\bar{\mathbf{X}}_{\varepsilon}(n\cdot)$ converges 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}'_{i,C}(t) = \lambda p_C \varepsilon^2 - \frac{1}{L_C} \frac{|C|\bar{x}_{i,C}(t)}{\sum_{k,U} |C \cap U| \ell^{\varepsilon}(a_k, a_i) \bar{x}_{k,U}} & \text{if } \bar{x}(t) \neq 0 \\ \bar{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 \ge 1\}$ and $\{\bar{\mathbf{Y}}^n, n \ge 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 \ge 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 $(\frac{1}{n_l}\bar{X}_\varepsilon(n_l\cdot), \bar{Y}^{n_l})$ 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 \ge 1\}$ and two processes $\hat{\mathbf{x}}_\varepsilon$ and $\hat{\mathbf{Y}}$ such that:

$$- (\hat{\mathbf{X}}^{n_l}, \hat{\mathbf{Y}}^{n_l}) \to (\hat{\mathbf{x}}_\varepsilon, \hat{\mathbf{Y}}), \hat{\mathbb{P}}\text{-almost surely;}$$

$$- (\hat{\mathbf{X}}^{n_l}, \hat{\mathbf{Y}}^{n_l}) \stackrel{d}{\sim} (\frac{1}{n_l}\bar{\mathbf{X}}_\varepsilon(n_l, \cdot), \bar{\mathbf{Y}}^{n_l}) \text{ for all } l \ge 1;$$

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

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

$$\hat{X}_{i,C}^{n_l}(t) = \hat{Y}_{i,C}^{n_l}(t) - \frac{1}{L_C} \int_0^t \hat{X}_{i,C}^{n_l}(s) R_{i,C}^{n}(\hat{\mathbf{X}}^{n_l}(s)) \mathrm{d}s. \quad (21)$$

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

$$\sup_{l \ge 1} \sup_{0 \le t \le T} \|X_{i,C}^{n_l}(t)\| < \infty.$$

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

$$\lim_{l \to \infty} R_{i,C}^{n_l}(x) = \frac{|C|}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| 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}})$:

$$\hat{x}_{i,C}(t) = \hat{Y}_{i,C}(t) - \frac{1}{L_C} \int_0^t \frac{|C| |\hat{x}_{i,C}(s)|}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| |\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{|C| |\hat{x}_{i,C}(s)|}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| |\hat{x}_{k,U}(s)|} \mathrm{d}s. \quad (22)$$

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

<figure><img src="image_9.png" /></figure>
## 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 \le \frac{1}{L_C} \frac{|C| z_{i,C}}{\sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| 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 $\|\mathbf{x}(0)\| = M$, we have $|\mathbf{x}(T)| < \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}$. $\square$

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{|C| z_{i,C}}{\sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| z_{k,U}}$, which is 0-homogeneous, and non-increasing in $x_{k,U}$ for $(k, U) \ne (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 \le \frac{1}{L_C} \frac{|C| p_C L_C}{\sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| p_U L_U} \quad (23) $$

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

$$ \lambda \le \frac{1}{\langle \ell_D^\varepsilon \rangle} \frac{|C|}{\sum_U |C \cap U| p_U L_U} \quad (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 $|C| = |D| \Rightarrow y_C = y_D$. We have, for all $C \in \mathcal{P}(K)$:

$$ \sum_{U \in \mathcal{P}(K)} |C \cap U| y_U = \frac{|C|}{K} \sum_{U \in \mathcal{P}(K)} |U| y_U \quad (25) $$

*Proof.* Let $C \in \mathcal{P}(K)$ with $|C| = j$, and $(y_U)_{U \in \mathcal{P}(K)} \in \mathbb{R}^{2^K - 1}$ be such that $y_U = y_V$ whenever $|U| = |V|$. We have:
$$ \sum_{U \in \mathcal{P}(K)} |C \cap U| y_U = \sum_{l=1}^{K} \sum_{m=1}^{j \wedge l} \sum_{\substack{U: |U|=l \\ |C \cap U|=m}} m y_U. \quad (26) $$

The number of sets $U \in \mathcal{P}(K)$ such that $|U| = l$ and $|C \cap U| = m$ is equal to $\binom{j}{l} \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_{|V|=l} y_V(t)$. We can rewrite (26) as:

$$ \sum_{U \in \mathcal{P}(K)} |C \cap U| y_U = \sum_{l=1}^{K} y_{[1,l]} \sum_{m=1}^{j \wedge l} m \binom{j}{m} \binom{K-j}{l-m} \\ = \sum_{l=1}^{K} \left( \sum_{|U|=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 $|U|$ and $j$ by $|C|$ leads to:

$$ \sum_{U \in \mathcal{P}(K)} |C \cap U| y_U = \frac{|C|}{K} \sum_{U \in \mathcal{P}(K)} |U| y_U, $$

which concludes the proof of Lemma 8.5. $\square$

The vector $(p_C)_{C \in \mathcal{P}(K)}$ verifies the conditions of Lemma 8.5, which allows us to obtain, after setting $\mathcal{L} = \sum_U |U| p_U L_U$:

$$ \lambda \le \frac{1}{\langle l^{\epsilon} \rangle_D} \frac{K}{\mathcal{L}}. \quad (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{X}_{\epsilon}$. $\bar{X}_{\epsilon}$ is an irreducible Markov jump process that is also positive recurrent. Thus, $\bar{X}_{\epsilon}$ is ergodic, which concludes the proof of Theorem 4.4.

## C: Proof of Theorem 4.5

To prove the instability of the chain $\bar{X}_{\epsilon}$, 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{|C|x_{i,C}}{\sum_{k,U} \ell_{\varepsilon}(a_k, a_i) |C \cap U|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 \to \infty} \mathbf{x}_n = \lim_{n \to \infty} \mathbf{y}_n = 0$. Using the 0-homogeneity of the functions, we have, for all $n$:

$$r_{0,\{1\}}(\mathbf{x}_n) = \frac{1}{p_{\{1\}} L_{\{1\}} \varepsilon^2}, \\ 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 $S = \{\mathbf{x} \in \mathbb{R}_{+,*}^{N_{\varepsilon} \times 2^K - 1} : |\mathbf{x}| = 1\}$, which is compact. Thus, it admits a maximum on $S$, and on $\mathbb{R}_{+,*}^{N_{\varepsilon} \times 2^K - 1}$ as a consequence. We set

$$S = \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{|C|x_{i,C}}{\sum_{k,U} \ell_{\varepsilon}(a_k, a_i) |C \cap U|x_{k,U}},$$

which is non-empty.

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

$$\lambda p_{C^*} \varepsilon^2 > \frac{1}{L_{C^*}} \frac{|C^*| z_{i^*, C^*}}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i^*) z_{k,U}}.$$

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

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

$$\lambda p_{C^*} \varepsilon^2 > \frac{1}{L_{C^*}} \frac{|C^*| z_{i^*, C^*}}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i^*) z_{k,U}}.$$

We define the process $\underline{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^*}} \frac{|C^*| z_{i^*, C^*}}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i) z_{k,U}} \text{ if } i, C = i^*, C^*,$$
$$ -\frac{1}{L_C} \frac{|C|Y_{i,C}(t)}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i) Y_{k,U}(t) + \mathcal{N}_0} \text{ else.} $$

$\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 \ge 0$:

$$ \frac{1}{L_{C^*}} \frac{|C^*| z_{i^*, C^*}}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i) z_{k,U}} \ge \frac{|C^*| \underline{X}_{i^*, C^*}(t)}{L_{C^*} \sum_{k,U} |C \cap U| \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) \le \underline{X}_{i,C}(t), \quad \mathbb{P}\text{-a.s.} $$

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

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

We can note that queue $\underline{X}_{i^*, C^*}$ is not the only queue whose population diverges to infinity. The departure rates in all queues (except for queue $i^*, C^*$) are decreasing functions of $\underline{X}_{i^*, C^*}(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 S$. Then, $\mathbf{z}$ is a solution to the following system of equations, for all $0 \le i, j \le N_{\varepsilon} - 1$ and $C, D \in \mathcal{P}(K)$:

$$ \frac{1}{p_C L_C} \frac{|C| z_{i,C}}{\sum_{k,U} \ell_{\varepsilon}(a_k, a_i) |C \cap U| z_{k,U}} = \frac{1}{p_D L_D} \frac{|D| z_{j,D}}{\sum_{k,U} \ell_{\varepsilon}(a_k, a_j) |D \cap U| z_{k,U}}. \quad (28) $$

*Proof.* Let $\mathbf{z}$ be such that:

$$ \mathbf{z} \in \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 \le i \le N - 1$ and $C \ne \{0\}$. Hence, these values have to be minimal in order to maximize $r_{0,\{1\}}$.
By definition, for all $1 \le j \le N_\epsilon - 1$, $C \ne \{1\}$, we have $r_{0,\{1\}}(\mathbf{z}) \le r_{i,C}(\mathbf{z})$. But, for $i \ge 1$ and $C \ne \{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.

<figure><img src="image_10.png" /></figure>

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:

$$ \begin{align*} \lambda p_{C^*} \varepsilon^2 \ge & \frac{1}{L_{C^*}} \frac{|C^*| p_{C^*} L_{C^*}}{\sum_{k,U} \ell_{\varepsilon}(a_k, a_{i^*}) |C^* \cap U| p_U L_U} \\ = & \frac{\langle a, \ell_{\varepsilon, \mathcal{D}} \rangle \frac{|C^*|}{K} \sum_U |U| p_U L_U}{\varepsilon^2 p_{C^*} K} \\ = & \frac{\langle \ell_{\varepsilon, \mathcal{D}} \rangle \frac{|C^*|}{K} \sum_U |U| p_U L_U}{\langle \ell_{\varepsilon, \mathcal{D}} \rangle \frac{|C^*|}{K} \sum_U |U| p_U L_U}, \end{align*} $$

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 $\mathcal{L} = \sum_U |U| p_U L_U$ and rearranging the equation gives us:

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

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