# Study of the dynamics of spatial point processes in wireless communication networks

Pierre Popineau

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

## To cite this version:

Pierre Popineau. Study of the dynamics of spatial point processes in wireless communication networks. Probability [math.PR]. Université Paris sciences et lettres, 2023. English. NNT: 2023UPSLE004. tel-04496165v2

HAL Id: tel-04496165

https://theses.hal.science/tel-04496165v2

Submitted on 8 Mar 2024

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
# Étude de la dynamique de processus ponctuels spatiaux appliqués aux réseaux de communication sans fil

# Study of the dynamics of spatial point processes in wireless communication networks

## Soutenue par

Pierre Popineau

Le 14/06/2023

## École doctorale n°386

Mathématiques Paris Centre

## Spécialité

Mathématiques

## Composition du jury :

<table><tr><td>Sara Alouf<br>INRIA Sophia Antipolis</td><td>Examinatrice</td></tr><tr><td>François Baccelli<br>INRIA Paris</td><td>Directeur de thèse</td></tr><tr><td>Marceau Coupechoux<br>LTCI - Télécom Paris</td><td>Président du jury</td></tr><tr><td>Sergey Foss<br>Heriot-Watt University</td><td>Examineur</td></tr><tr><td>Jean-Marie Gorce<br>INSA Lyon</td><td>Rapporteur</td></tr><tr><td>Derya Malak<br>EURECOM</td><td>Examinatrice</td></tr><tr><td>Pascal Moyal<br>IECL - Université de Lorraine</td><td>Rapporteur</td></tr><tr><td>Luis Guilherme Uzeda Garcia<br>Nokia Bell Labs France</td><td>Examineur</td></tr></table>

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
# RÉSUMÉ

Grâce aux nouveaux paradigmes introduits dans la dernière génération de réseaux sans fil, les attentes concernant le temps de service, la latence et la performance du réseau ont augmenté. Pour modéliser ces réseaux, la théorie des processus ponctuels et la géométrie stochastique se sont avérées utiles car elles fournissent un cadre polyvalent et robuste pour obtenir des résultats lorsque l'on travaille avec ces réseaux sans fil. L'ajout d'une dynamique markovienne pour modéliser les connexions et les temps de service complète le cadre d'analyse de ces réseaux sans fil.

La première contribution du travail présenté dans cette thèse réside dans l'analyse de la différenciation des services : les réseaux 5G NR ont introduit le partitionnement de la bande passante comme outil pour augmenter la performance du réseau. Dans cette configuration de réseau, tous les utilisateurs n'interfèrent pas les uns avec les autres avec la même puissance : les utilisateurs qui émettent avec un spectre de fréquence d'émission plus large auront une plus grande bande passante, mais ils rencontreront également plus d'interférences de la part des autres utilisateurs dans le réseau. En revanche, les utilisateurs dont le spectre est plus étroit subiront moins d'interférences. Nous définissons un cadre markovien pour étudier un tel processus spatial multiclasse de naissance et de mort, et nous décrivons sa région de stabilité. Pour de tels systèmes, les propriétés du régime stationnaire sont analysées, telles que les mesures moment ou l'attraction, ce qui permet de mieux comprendre cette dynamique.

Le deuxième problème que nous examinons est celui de la mobilité, qui est devenue une caractéristique importante des réseaux sans fil en raison de l'utilisation d'antennes hautement directionnelles. En utilisant une architecture simple pour un réseau cellulaire à deux niveaux, nous étudions deux familles de politiques d'association : une première famille qui s'appuie uniquement sur la mobilité de l'utilisateur, et une seconde qui offre un compromis entre la géométrie du réseau et la mobilité de l'utilisateur afin d'augmenter les performances du réseau. Ces politiques sont ensuite comparées à une politique d'association classique de puissance maximale afin d'évaluer leurs performances.
<figure><img src="image_5.png" /></figure>
ABSTRACT

Thanks to the new paradigms introduced in the latest generation of wireless
networks, expectations concerning service time, latency and network perfor-
mance have increased. To model such networks, point process theory and
stochastic geometry have proven to be useful as they provide a versatile and
robust framework to obtain results when working with such wireless net-
works. Adding to this Markov dynamics to model connections and service
times completes this framework to analyze such wireless networks.

The first contribution of the work presented in this thesis lies in the anal-
ysis of service differentiation: 5G NR networks have introduced bandwidth
partitioning as a tool to increase network performance. Under this network
setup, not all users interfere with each other with the same power: users
transmitting with a broader transmitting frequency spectrum will have a
larger bandwidth, but they will also encounter more interference from the
other users in the network. In contrast, users with a narrower spectrum will
experience less interference. We define a Markovian framework to study such
a multiclass spatial birth-and-death process, and we describe its stability re-
gion. For such systems, properties of the stationary regime are analyzed,
such as moment measures or statistical clustering, leading to a better under-
standing of these dynamics.

The second problem we look into is mobility, which has become an impor-
tant feature in wireless networks due to the use of highly directional antennas.
Using a simple architecture for a two-tier cellular network, we study two fam-
ilies of association policies: a first family which only relies on user mobility,
and the second offers a trade-off between network geometry and user mobil-
ity to increase network performance. These policies are then compared to a
classical max-power association policy to assert their performance.
<figure><img src="image_6.png" /></figure>
ACKNOWLEDGEMENTS

The first person I want to thank is my supervisor, François Baccelli, for trust-
ing and accompanying me throughout my PhD. When I started working in
NEMO, I had just graduated from Ecole des Mines de Paris, and had very
little experience in research, especially in stochastic geometry. He decided
to accept me initially as a junior researcher in the team and later as a PhD
student to pursue the work I am presenting in this manuscript. François'
valuable experience in the field and his extremely precise remarks were in-
valuable to me. They guided me throughout my PhD and helped me gain a
better understanding of the concepts I encountered during my research.

My second word of thanks goes to the people of ERC NEMO. Firstly, my
co-PhD students and postdocs: Michel Davydov, Bharath Roy Choudhury,
Sayeh Khaniha, Ke Feng, and Guodong Sun. Even though our research
subjects are relatively different, we maintained unity through our reading
groups, seminars, but most importantly, through random discussions in our
shared office whenever one of us was stuck on a problem (which happened to
be me, most of the time). Secondly, I would like to thank Bartek Blaszczyszyn
for accepting my application as a junior researcher at the beginning of my
journey, and for guiding the operations at the head of the Dyogene team,
Ali Khezeli for animating the scientific life of our group by organizing our
seminar series and Julien Guieu, who has been of great help at the interface
with INRIA. I also wanted to thank Sara Alouf, Marceau Coupechoux, Sergey
Foss, Jean-Marie Gorce, Derya Malak, Pascal Moyal and Luis Uzeda Garcia
for accepting to be a part of the jury for this PhD.

During these years, I divided my time between working at INRIA Paris
and LINCS (Laboratory on Information, Networking, and Communication
Sciences, hosted by Telecom Paris). LINCS holds a special place in these ac-
knowledgments as it is a unique place of scientific exchange, where I had the
chance to meet researchers from various fields and have exciting conversa-
tions. I especially want to thank Anna Frisone for her presence and valuable
help with administrative and non-administrative matters. I also thank all
LINCS researchers whom I met, worked with, talked to, played Concept with,
or had coffee with: Fabien Mathieu, François Durand, Ludovic Noirie, Ana
Bušić, Francesca Bassi, Leonardo Linguaglossa, Marc-Olivier Buob, Rémi
Varloot, Maxime Mouchet, Céline Comte, and many more.

Little did I know, when I registered as a PhD student in February 2020, what the future held for us. Starting a PhD during a pandemic was a peculiar experience, as it is a pivotal moment where I needed guidance the most. But I held on tight, and after the end of the lockdown, I was extremely lucky to have access to an office at LINCS, located then inside Paris, to have a dedicated workplace to escape my small apartment. These times were special, as having a whole laboratory for mostly three persons—Anna, Fabien, and I—was strange, but we made the most out of it.

Next, I would like to thank my co-authors. Sanket Kalamkar was of great help when it came to helping me write my first article and gave me very valuable advice during our discussions. Seva Shneer, whom I had the chance to visit in Edinburgh, had very constructive discussions with me over the past year and a half and provided input on earlier versions of some of the proofs presented here, for which I am grateful.

Finally, I want to express my last word of thanks to my parents, Joëlle and Fabrice, my brother, Maxime, and all of my friends: Irène, Valentine, Guillaume B., Léa, Guillaume dL., Romy, Héloïse, and the members of the West Coast Swing and Rock 4 Temps communities in Paris, and more generally, anyone that kept listening when I randomly started talking about mathematics at some point in a discussion. Your support during these years of PhD was invaluable to me, and I am grateful to have had you by my side, especially during the last year of this PhD.

Many thanks to anyone who embarks on the adventure of discovering this work, and I wish you a pleasant reading!
# INTRODUCTION

Wireless networks have seen an important growth over the past decades, driven by the increasing demand for mobile communication services and the development of wireless technologies. This growth was accompanied by a number of constraints resulting from the intrinsic properties of wireless networks, where network performance depends on the spatial configuration of transmitting users at a given moment. The latest generation of wireless networks have brought another layer of constraints, notably with the introduction of multi-tier networks and service differentiation. These considerations have brought a multitude of questions about network load, communication latency and transmission performance, and the need for robust frameworks to model such networks.

Communication networks have been studied extensively in queuing theory, as information about user latency and service times has become increasingly interesting. Among the properties of interest, stochastic stability ensures to avoid congestion or overload in the network, and provides useful bounds in population of the network. Under the Markov property, the stability condition is also important to ensure the convergence of latency to a stationary distribution which one wants to characterize analytically. Jackson networks ([55]) and their multi-class extension, BCMP networks ([24]), fall within the classes of studied networks, but stability results about networks with service differentiation as intended in this work remains an open question.

Another tool that has risen in popularity to study and model communication networks is stochastic geometry ([12], [48]), as it provides a flexible and reliable framework to encapsulate the features of wireless communications. This has become especially true in the most recent wireless networks, as the standards set in terms of reliability and latency have become very high, for which stochastic geometry has proven to be useful ([53]). A vast literature exists on stochastic geometry ([86]) and its applications to wireless networks ([18]), as the framework it offers allows to obtain bounds and performance guarantees in wireless networks ([42]).

This manuscript will tackle two different aspects of wireless networks combining these two approaches. In the two first parts, we consider a dy-
namical system network modeled on Markov jump process, where immobile
users arrive and leave upon the completion of a transmission, for which we
establish conditions for its stochastic stability and provide results about its
steady-state dynamics through Palm calculus and stochastic geometry. In
the last part of this thesis, we will look to a fixed network with mobile users
and look to obtain results about network performance when taking into ac-
count user velocity. The work presented here is seen as a first step towards
designing a comprehensive mathematical framework to encapsulate all the
features of 5G NR networks and their future releases.

Stochastic stability in communication
networks

Spatial interactions and spatial interacting queuing systems have been stud-
ied for a long time, as they prove to be useful to model a large variety of
physical phenomena, ranging from physics to epidemiology, biology or com-
munication theory. Queuing systems have also been studied extensively, as
many results have been obtained about their behavior, especially under the
Markov assumption, which provides a framework where several additional
properties can be reached (existence of a stationary distribution, converging
times, etc.). Adding spatial interactions to queuing systems defines the the-
ory of *spatial queuing networks* ([82]), which form a class of queuing systems
where the locations of queues are not fixed locations in space, and where
the arrival and departure times depend on the geometry of the system at a
given time. Such systems have been introduced by Preston ([75]) as jump
processes over the set of counting measures on a given region of space.

A particularity of wireless networks comes from the nature of spatial
interactions between users. In a device to device wireless network, users
arrive at random locations in time and space and their service rate at any
time is a function of the interference created by other users transmitting at
the same time in the network. In [80], the authors introduced a framework
to study such wireless spatial birth-and-death networks, derived the stability
region, i.e, the set of arrival rates for which the network is stochastically
stable, which happens to be an interval, and studied the stationary regime
of this class of dynamics. On the other hand, when working with cellular
networks, the study of stability happens to be more complex, as metastable
regimes may emerge ([4]). Both problems have been studied in a setup where
all users use the same communication channels, which ensures that all users
interact with each other.

Various instances of spatial queuing networks have been studied and re-
sults have been obtained using queuing theory tools to assert the recurrence
and stability conditions for such systems in a Markov setup ([3], [16], [17],
[29]). Among the numerous fields in which spatial queuing systems have been
studied, communication systems provide an interesting framework, as many
communication protocols prove to respect the Markov property (ALOHA,
DMAA under natural assumptions - [85]).

Service differentiation in wireless networks

The latest generations of wireless networks, such as 4G LTE and 5G NR ([9]), have introduced service differentiation as a means to increase bandwidth and network performance ([39]), and to decrease user latency and increase communication reliability ([62]). This service differentiation uses bandwidth partitioning algorithms ([67]) in order to assign to users communication channels matching their needs and respecting certain fairness criteria.

In the multiclass wireless networks considered here, users can transmit on several frequency bands. Those using a large number of bands transmit with a strong signal, but generate a lot of interference for other users, hence slowing down their transmission. On the other hand, users transmitting on a small number of bands have a weaker signal, but generate less interference and block less the other transmissions. It also to note that some users do not interfere with each other, as the sets of frequency bands allocated to them do not overlap, while others may interfere with all present users, e.g. when they use the whole available spectrum of bands to transmit.

A first goal of this thesis is to develop frameworks designed at modelling dynamic multiclass wireless networks, and studying some of their properties: stability, stationary regime and establishing results about network performance based on the parameters. To capture the space-time interactions between users, we use a birth-and-death process of wireless dipoles (see [12], [47]). In this model, dipoles arrive according to a Poisson rain of constant intensity, representing the arrival rate of users in the network, and leave with a rate inversely proportional to a shot-noise of the configuration. These dynamics are standard when studying wireless networks and represent a linearized version of the Shannon-Hartley formula in the interference limited case, corresponding to low SINR, modelling the capacity of a transmission channel ([83]). This relatively simple framework captures some of the key features of bandwidth partitioning in a 5G network and has been an open question for several years.

As mentioned, the novelty of this works lies in the study of this network of interacting point processes. Spatial birth-and-death processes in the single-class setup have been studied in [80], but they cannot be immediately generalized to a multiclass network. The main difficulty lies in the presence of multiple available communication channels, which we have to take into account in the proofs. In this work, we will restrict our study to a network in a square torus. Studying spatial queuing networks in a finite region of space allows for an easier study, as most classical methods fail or are harder to generalize when considering the whole Euclidean plane (see [16], [81] for instance). This thesis is a first step towards studying more complex networks: the hypothesis considered in this work are restrictive compared to real life systems, but the network setup we decided to study gathers some key features of wireless communication networks and bandwidth sharing and
partitioning policies used in telecommunication protocols, while allowing us
to establish a methodology for their study.

Association policies in heterogenous networks

Another question of interest when considering dynamics of wireless networks is the study of heterogenous and multi-tier networks. Since the introduction of microcells, picocells until the emergence of femtocells in 4G LTE ([57]), association policies have become a key feature of wireless networks. Such multitier network architectures have emerged as a means to improve signal power and reduce interference ([39]), as they allow the network load to be spread over multiple parallel instances of a network ([40]). Multiple types of association policies have been studied ([65]), but max-power association remains a good performer ([79]).

Among the important features an association policy has to take into ac-
count, user mobility has become a prominent factor with the democratization
of beam-based communications, which rely on narrow beams to concentrate
transmission power. In this network paradigm, when a user moves out of the
beam it uses to communicate with its serving base station, it has to reselect
a new beam, which halts the transmission with the network. These *beam re-
selections* are added to the cell handover, performed when a user moves out
of the coverage of a base station to create an *overhead* in the transmission,
significantly decreasing the network performance over a transmission. Thus,
designing relevant association policies when studying such networks have to
take user mobility into account.

In this line of thoughts, association policies need to be rather simple as
more complex association policies may prove to be computationally heavy
and may need a controller to associate users to different tiers, which be-
comes less interesting with the increase of the number of connected users to
a network. As a result, decentralized, user-centric association scheme have
become prominent in massive MIMO networks ([25]), which have to consider
a large amount of connections simultaneously. We will propose a simple 2-
tier wireless networks with macro and micro base stations, capturing some
essentials features of 5G NR networks (directional beamforming, mobility
and beam misalignment), using a similar analysis as developed in [53].

Organization of the manuscript

This manuscript consists in three parts. The first part of this work is dedi-
cated to establishing a framework to study the stability of a multiclass spatial
queuing system with wireless interactions. There, we extend a framework in-
troduced in [80] to adapt it to the multiclass setup. In a first step, we
introduce a symmetry assumption to simplify the study of the system and
allow us to obtain a closed form fo the critical arrival rate in the system.
This study relies on the use of the definition of adequate interference queue-
ing networks (see [81]) to obtain bounds for the stability region. Due to the
nature of these dynamics, fluid limits ([78]) provide a useful framework to
obtain conditions for stability ([37]). In a second step, we relax this sym-
metry assumption and manage to obtain a general stability condition for
spatial queuing systems for which fluid limits follow a certain evolution, giv-
ing a reciprocal to a result presented in [84]. Finally, we present numerical
simulations for our result and discuss some natural extensions of the model.

The second part focuses around the study of the stationary regime of
the dynamics presented in Part I. We start by presenting two heuristics to
estimate the spatial user densities in the stationary regime with the use of
stochastic geometry tools (see [12]) as well as Palm calculus (see [14]) 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. These approximations are known to provide a powerful framework
to obtain both qualitative and quantitative results in many cases. From [80]
we know that a monotype wireless network displays clustering, i.e. attraction
between users in the stationary regime, which leads us to investigate this
property in the multiclass regime, as the spatial interaction between users
are different. Finally, we prove the existence of exponential moments for the
system. This last property allows us to bound the tails of the steady-state
distribution of the dynamics, which have an important practical implications
for wireless communications.

The third part revolves around the study of association policies for a
two-tier heterogeneous wireless network, based on a study presented in [53].
Because mobility has become a prominent feature to take into account when
designing wireless networks ([10], [23]), we introduce a category of associ-
ation policies based on user velocity to maximize the spectral efficiency in
the network. This class of association policies is interesting as they manage
to outperform the classical max-power association policies in some network
setups. Because these policies only rely on the user velocity and not on net-
work geometry, we introduce *velocity-biased max power* policies to take into
account these two components and we propose a computational framework
to study them under a set of networks parameters.

Included papers

- P. Popineau and F. Baccelli, *On multiclass spatial birth-and-death process with wireless type interactions*, https://arxiv.org/abs/2205.15799, submitted for publication to IEEE Transactions on Information Theory.

- P. Popineau, S. S. Kalamkar and F. Baccelli, *On velocity-biased association policies for multi-tier 5G wireless networks*, 2021 IEEE Global Communications Conference (GLOBECOM), Madrid, Spain, 2021, pp. 1-6, https://ieeexplore.ieee.org/document/9685210.

- P. Popineau, V. Shneer, *An instability condition for queuing systems*
with state-dependent departure rates, https://arxiv.org/abs/2304.08853,
submitted for publication to Queueing Systems.

Acknowledgements

The work presented here is supported by the ERC NEMO grant, under the
European Union's Horizon 2020 research and innovation programme, grant
agreement number 788851 to INRIA.
# CONTENTS

<table><tr><td>Résumé</td><td>i</td></tr><tr><td>Abstract</td><td>iii</td></tr><tr><td>Acknowledgements</td><td>v</td></tr><tr><td>Introduction</td><td>vii</td></tr><tr><td>Stability theory in communication networks</td><td>viii</td></tr><tr><td>Service differentiation</td><td>ix</td></tr><tr><td>Organization of the manuscript</td><td>x</td></tr><tr><td>Contents</td><td>xiii</td></tr><tr><td>List of Figures</td><td>xvii</td></tr><tr><td>A Mathematical Foundations</td><td>1</td></tr><tr><td>A.1 Notation</td><td>1</td></tr><tr><td>A.2 Point process theory</td><td>3</td></tr><tr><td>A.3 Poisson-Voronoi Tessellations</td><td>9</td></tr><tr><td>A.4 Markov Jump Processes</td><td>14</td></tr><tr><td>A.5 Stability theory</td><td>17</td></tr><tr><td>B Wireless Networks</td><td>23</td></tr><tr><td>B.1 Channel capacity and Shannon rate</td><td>23</td></tr><tr><td>B.2 Wireless networks and stochastic geometry</td><td>25</td></tr><tr><td>B.3 Beam-based communications</td><td>26</td></tr><tr><td>B.4 Bandwidth partitioning</td><td>29</td></tr><tr><td>I Multiclass spatial queuing networks</td><td>31</td></tr><tr><td>1 Mathematical framework</td><td>33</td></tr><tr><td>1.1 Network setup: bandwidth partitioning</td><td>33</td></tr><tr><td>1.2 Monotonicity and irreducibility</td><td>35</td></tr></table>
<table><tbody><tr><td>1.3</td><td>Critical arrival rate</td><td>39</td></tr><tr><td>2</td><td>Stability of symmetric systems</td><td>41</td></tr><tr><td>2.1</td><td>Discretization of the dynamics</td><td>42</td></tr><tr><td>2.2</td><td>Fluid model and fluid limits</td><td>43</td></tr><tr><td>2.3</td><td>Main result</td><td>44</td></tr><tr><td>3</td><td>Proofs of Theorems of Chapter 2</td><td>51</td></tr><tr><td>3.1</td><td>Proof of Theorem 2.3</td><td>51</td></tr><tr><td>3.2</td><td>Proof of Theorem 2.4</td><td>56</td></tr><tr><td>3.3</td><td>Proof of Theorem 2.5</td><td>58</td></tr><tr><td>4</td><td>Generalization</td><td>61</td></tr><tr><td>4.1</td><td>General queuing systems</td><td>62</td></tr><tr><td>4.2</td><td>An instability condition</td><td>63</td></tr><tr><td>4.3</td><td>Stability of non-symmetric BWP models</td><td>65</td></tr><tr><td>5</td><td>Extension and conclusion</td><td>69</td></tr><tr><td>5.1</td><td>Interpretation of the results</td><td>69</td></tr><tr><td>5.2</td><td>Restricting band configurations</td><td>70</td></tr><tr><td>5.3</td><td>Extension</td><td>71</td></tr><tr><td>II</td><td>Stationary interacting point processes</td><td>73</td></tr><tr><td>6</td><td>Study of the steady state dynamics</td><td>75</td></tr><tr><td>6.1</td><td>First order heuristic and Poisson estimate for λ<sub>c</sub></td><td>76</td></tr><tr><td>6.2</td><td>Second order heuristic</td><td>78</td></tr><tr><td>6.3</td><td>Numerical simulations and performance of the heuristics</td><td>80</td></tr><tr><td>7</td><td>Clustering</td><td>83</td></tr><tr><td>8</td><td>Higher order moment measures</td><td>89</td></tr><tr><td>III</td><td>Velocity-based association policies</td><td>95</td></tr><tr><td>9</td><td>Two-tier 5G wireless network</td><td>97</td></tr><tr><td>9.1</td><td>Network setup</td><td>97</td></tr><tr><td>9.2</td><td>Shannon rate</td><td>99</td></tr><tr><td>9.3</td><td>Mean Effective Shannon rate under an association policy</td><td>100</td></tr><tr><td>10</td><td>Velocity-based association policies</td><td>101</td></tr><tr><td>10.1</td><td>Coverage probability</td><td>101</td></tr><tr><td>10.2</td><td>Optimal threshold association policy</td><td>103</td></tr><tr><td>10.3</td><td>Load-dependent threshold policies</td><td>104</td></tr><tr><td>10.4</td><td>Discussion and comparison</td><td>106</td></tr><tr><td>11</td><td>VBMP association policies</td><td>109</td></tr></tbody></table>
CONTENTS

xv

<table><tbody><tr><td>11.1 VBMP policies</td><td>109</td></tr><tr><td>11.2 Optimal VBMP policy</td><td>111</td></tr><tr><td>11.3 Load-dependent VBMP policy</td><td>113</td></tr><tr><td>11.4 Discussion and comparison</td><td>116</td></tr><tr><td><b>12 Proof of Theorems of Chapters 10 and 11</b></td><td><b>119</b></td></tr><tr><td>12.1 Proof of Theorem 10.1</td><td>119</td></tr><tr><td>12.2 Proof of Theorem 10.2</td><td>121</td></tr><tr><td>12.3 Proof of Theorem 11.1</td><td>122</td></tr><tr><td>12.4 Proof of Lemma 11.3</td><td>124</td></tr><tr><td><b>13 Conclusion</b></td><td><b>125</b></td></tr><tr><td><b>14 Résumé du manuscrit en français</b></td><td><b>131</b></td></tr><tr><td><b>Bibliography</b></td><td><b>141</b></td></tr></tbody></table>
<figure><img src="image_7.png" /></figure>
LIST OF FIGURES

<table><tr><td>A.1</td><td>Size of the 0-cell against the size of a typical cell</td><td>12</td></tr><tr><td>A.2</td><td>Performance of the heuristic of Lemma A.8</td><td>14</td></tr><tr><td>B.1</td><td>Basic communication system</td><td>23</td></tr><tr><td>B.2</td><td>Sector approximation of the antenna gains</td><td>27</td></tr><tr><td>B.3</td><td>Illustration of the total overhead for a moving MU</td><td>28</td></tr><tr><td>B.4</td><td>Physical layer of 5G NR networks</td><td>30</td></tr><tr><td>2.1</td><td>Population over time in the network with K = 2</td><td>47</td></tr><tr><td>2.2</td><td>Average sojourn time in the system over time with K = 2</td><td>48</td></tr><tr><td>2.3</td><td>Average staying time in the system over time with K = 4</td><td>49</td></tr><tr><td>4.1</td><td>Population over time in the non-symmetric network with K = 2</td><td>68</td></tr><tr><td>6.1</td><td>Relative performance of heuristics compared to simulations</td><td>80</td></tr><tr><td>6.2</td><td>Performance of the heuristic for λ<sub>c</sub></td><td>81</td></tr><tr><td>8.1</td><td>Population distribution in the network</td><td>93</td></tr><tr><td>8.2</td><td>Population tails in the symmetric system, with λ = 0.95λ<sub>c</sub></td><td>93</td></tr><tr><td>10.1</td><td>Coverage probability in each tier</td><td>102</td></tr><tr><td>10.2</td><td>Normalized ESR in the network as a function of v<sub>T</sub></td><td>103</td></tr><tr><td>10.3</td><td>Comparison of the heuristics for the optimal threshold policy</td><td>105</td></tr><tr><td>10.4</td><td>Optimal velocity threshold as a function of MU density</td><td>107</td></tr><tr><td>10.5</td><td>Comparison between MARP and threshold policies</td><td>107</td></tr><tr><td>11.1</td><td>Function K<sup>*</sup></td><td>112</td></tr><tr><td>11.2</td><td>Performance of the optimal VBMP policy</td><td>113</td></tr><tr><td>11.3</td><td>Optimal function K<sup>†</sup></td><td>116</td></tr><tr><td>11.4</td><td>Comparison of the three association policies</td><td>117</td></tr><tr><td>11.5</td><td>Comparison of the three integration procedures</td><td>118</td></tr><tr><td>14.1</td><td></td><td>136</td></tr><tr><td>14.2</td><td></td><td>138</td></tr><tr><td>14.3</td><td></td><td>139</td></tr></table>
<figure><img src="image_8.png" /></figure>
# MATHEMATICAL FOUNDATIONS

In this chapter, we will detail the mathematical framework as well as some
basics properties, results and theorems about the mathematical objects that
will be used in this thesis.

## A.1 Notation

Firstly, let us define some notations. 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. Let $S$ be a second-countable locally compact Hausdorff space. In the
manuscript, $S$ will either be equal to the positive real line $\mathbb{R}_+$, the Euclidean
plane $\mathbb{R}^2$ or a compact subset of it. Let $\mathcal{S}$ be the canonical Borel $\sigma$-algebra
associated with $S$. Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a probability space, and we define $\mathcal{M}(S)$
the set of counting measures on $S$.

For $N \in \mathbb{N}$, we denote by $D([0, \infty), \mathbb{R}^N)$ the set of càdlàg (right-continuous with left limits) functions from $\mathbb{R}_+$ to $\mathbb{R}^N$. We consider stochastic processes as measurable maps from $(\Omega, \mathcal{F})$ to $(D([0, \infty), \mathbb{R}^N), \mathcal{S})$. Let $d_\infty^0$ be the infinite norm in $\mathbb{R}^d$, $|x|$ be the $L^1$ norm and $\|x\|$ denote the $L^2$ norm of $x$.

For a sequence $(x_n)_{n \ge 0}$ and $x$ in $D([0, \infty), \mathbb{R}^N)$, 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. \qquad (A.1) $$

In the rest of this thesis, if no mode of convergence is specified, this
definition of convergence is used.

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 \leq_i G$ if and only if

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

for all $f \in I$. $\leq_i$ is a partial order on the set of random vectors of $\mathbb{R}^d$. The next table resumes all the notations used in this chapter.

Let $\Psi$ and $\Psi'$ be two point processes on $\mathcal{D}$. We write $\Psi \leq \Psi'$ if we have $\Psi(\mathcal{D}) \leq_i \Psi'(\mathcal{D})$. In the rest of this thesis, we will use this partial order when comparing stochastic processes. Finally, we say that a random variable $X$ with values in $\mathbb{R}^d$ dominates $Y$ if the cumulative distribution functions of $X$ and $Y$, denoted as $F_X$ and $F_Y$, are such that $F_Y \leq_i F_X$. A consequence of stochastic domination is, for all $x \in \mathbb{R}^d$ and for all probabilities $P$ defined on $\mathbb{R}^d$, we have $P[X > x] \geq P[Y > x]$.

Finally, for two random variables $X$ and $Y$, we write $X \stackrel{d}{\sim} Y$ if $X$ and $Y$ have the same probability distribution.

## A.2 Point process theory

The main mathematical object we manipulate in this manuscript are point process. The most common way to visualize a point process is to consider a random collection of points of a set, whether it is the real line, the Euclidean plane or any mathematical space. Point processes have been introduced as a means to model and study a vast variety of physical processes, like the position of stars in the sky (which gave the idea to John Mitchell in 1767 of what would become the Poisson point process). Nowadays, point processes are used to study a variety of physical phenomena, ranging from biology ([41], [90]) and epidemiology ([66]) to wireless communication networks ([80]). This section aims at presenting point processes and some basics properties and results of point processes that will be used throughout this manuscript.

The first object we define is the point process:

**Definition A.1.** A point process on $S$ is a map $\Phi : \Omega \times S \mapsto \mathbb{N}$ such that:

- For all $\omega \in \Omega$, $\Phi(\omega, \cdot)$ is a locally finite counting measure on $S$, i.e. for all $x \in S$, there exists an open neighborhood $N_x$ of $x$ such that $\Phi(\omega, N_x) < +\infty$.

- For every $B \in \mathcal{S}$, $\Phi(\cdot, B)$ is a random variable on $\mathbb{N}$.

A point process $\Phi$ is a random variable over $M(S)$ over the set of locally finite counting measures on $S$ in a sense that for all $\omega$, $\Phi(\omega, \cdot)$ is a counting measure on $\mathbb{N}$, selected randomly. In this thesis, all point processes will
<table><caption>Table A.1: Table of notations of Chapter A</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td>S</td><td>A second-countable, separable Hausdorff space (assimilated to ℝ+, ℝ<sup>2</sup> or any compact subset of ℝ<sup>2</sup>)</td></tr><tr><td>S</td><td>Canonical Borel σ-algebra associated with S</td></tr><tr><td>(Ω, F, P)</td><td>A probability space</td></tr><tr><td>D(A, B)</td><td>Set of càdlàg function from A to B</td></tr><tr><td>d<sub>∞</sub><sup>0</sup></td><td>Infinite norm in ℝ<sup>d</sup></td></tr><tr><td>|·|</td><td>L<sup>1</sup> norm</td></tr><tr><td>||·||</td><td>L<sup>2</sup> norm</td></tr><tr><td>M(S)</td><td>Set of counting measures on (S, S)</td></tr><tr><td>Φ</td><td>Point process</td></tr><tr><td>δ</td><td>Dirac measure</td></tr><tr><td>ν</td><td>Lebesgue measure on S</td></tr><tr><td>Λ</td><td>Intensity measure of Φ</td></tr><tr><td>C</td><td>Campbell measure of a marked point process</td></tr><tr><td>P<sub>Φ</sub><sup>0</sup></td><td>Palm probability measure with respect to the point process Φ</td></tr><tr><td>E<sub>Φ</sub><sup>0</sup></td><td>Palm expectation with respect to the point process Φ</td></tr><tr><td>★</td><td>Convolution of probability distributions</td></tr><tr><td>K(·, ·)</td><td>Kernel of a Markov jump process</td></tr><tr><td>λ(x)</td><td>Jump intensity function of a jump process</td></tr><tr><td>Π(·, ·)</td><td>The transition matrix of a jump process</td></tr><tr><td>x̄</td><td>Fluid limit associated with the Markov chain X</td></tr><tr><td>ℓ</td><td>Path loss function</td></tr><tr><td>N<sub>0</sub></td><td>Thermal noise density in the network</td></tr></tbody></table>

be considered as random locally finite counting measures, and for simplicity, we will write $\Phi(\cdot)$ to denote either the random variable $\Phi$ or one of its realizations.

A common way to visualize point processes is as a random set of points in $S$: let $X_1, \cdot, X_n$, where $n$ may be infinite and the $X_i$ may not necessarily be distinct, be the points in $S$ for which $\Phi(\{X_i\}) > 0$. Then, to represent the point process, one can write $\Phi = \{X_1, \dots, X_n\}$ to emphasize this interpretation. Because point processes are also random measures on $S$, we will also use two other representations, which involve measure theory. The first one is:

$$ \Phi = \sum_{i=1}^{n} \delta_{X_i}, $$

where $n$ is an integer-valued random variable denoting the number of points of the realization of $\Phi$ and $\delta$ is the Dirac measure. The second representation for point process we will use is the counting notation. Let $T$ be a compact subset of $S$. Then, we define:
$$ \Phi(T) = \int_{x \in T} \Phi(du). $$

This way, $\Phi(T)$ represents the number of points of $\Phi$ in $T$. This counting notation will be particularly useful when studying stationary point processes in Part II. We will also use the classical set operations on point processes: $x \in \Phi$ means that $x$ belongs to a realization of the point process $\Phi$, and we will use the notations $\Psi \subset \Phi$, $\Phi \cup \Psi$, $\Phi \cap \Psi$ and $\Phi \setminus \Psi$ in their usual sense to respectively denote inclusion, union, intersection and difference between the random sets $\Psi$ and $\Phi$. Finally, the translation operator $x + \Phi$, where $x \in S$ and the dilation operator $\lambda\Phi$ with $\lambda \in \mathbb{R}$ keep their usual meaning.

If, for all $x \in S$, we have $\Phi(\{x\}) \le 1$, then the point process is said to be *simple*. A point process is said to be *stationary* if for all $y \in S$, the process $\tilde{\Phi} = \Phi + y$, defined as the set $\{y + X_i, X_i \in \Phi\}$ has the same distribution as $\Phi$. A point process is said to be *isotropic* if it is invariant by rotation, i.e. if $R(\Phi)$ has the same distribution as $\Phi$, where $R$ is a rotation of $S$. The study of stationary point processes is a vast field of research in itself, as it has a lot of interesting properties (see for instance Chapter 4 from [86] or Chapter 1 from [14]).

Two important quantities for point processes are their *intensity measure* and their *distribution*. The intensity measure of a point process is a measure on $S$ defined, for all $B \in \mathcal{S}$, by:

$$ \Lambda(B) = \mathbb{E}[\Phi(B)]. $$

In other words, the intensity measure of a point process is the mean number of points of the point process is a given Borel set $B$. If the point process is stationary, this intensity measure is translation-invariant. An important result in measure theory is that every isometry invariant Radon measure is absolutely continuous with regard to the Lebesgue measure, which can be applied here. As such, if the point process is stationary, there exists $\lambda \in \mathbb{R}$ such that for all Borel set $B$:

$$ \Lambda(B) = \lambda\nu(B), $$

where $\nu$ is the Lebesgue measure on $S$. $\lambda$ is called the *intensity* of the point process, and it can be interpreted as the average number per unit of space of the process, or its spatial intensity.

The distribution of the Point process is defined as $(\mathbb{P}[Y \in \Phi])$ for all the configuration sets $Y$. This distribution uniquely determines the point process.

## The Poisson point process

Among existing point processes, a common instance is the Poisson point process (PPP). The main idea behind the Poisson point process is to model random events sparsely separated spatially or temporally.

To define a Poisson point process, let us take a locally finite intensity measure $\Lambda$ on $S$, and we define the Point process $\Phi$ such that:
- For all $B \in \mathcal{S}$, $\Phi(B)$ is distributed according a Poisson law, i.e., for all $n \in \mathbb{N}$:

$$ \mathbb{P}[\Phi(B) = n] = \frac{\Lambda(B)^n}{n!} e^{-n}, $$

- The number of points in $n$ disjoint Borel subsets of $S$ are independent random variables, i.e., for all $B_1, \dots, B_n \in \mathcal{S}$ such that $i \neq j \Rightarrow B_i \cap B_j = \emptyset$, we have:

$$ \mathbb{P}[\Phi(B_1) = m_1, \dots, B_n = m_n] = \prod_{i=1}^{n} \mathbb{P}[\Phi(B_i) = m_i]. $$

The second property of Poisson point processes is also referred to as *independent scattering* or *increment independence*, as it captures the lack of interaction between of points in the Poisson process. A consequence of this property is that the points of a Poisson point process are uniformly distributed on $S$.

If there exists $\lambda > 0$ such that $\Lambda = \lambda\nu$, where $\nu$ is the Lebesgue measure on $S$, then we say that the Poisson point process is *homogenous*. An interesting property of homogenous Poisson point processes is that, when used to model arrival times of events the independent scattering tells us that the time between events is exponentially distributed with mean $1/\lambda$.

A second interesting property of the Poisson point process is that it is both stationary and isotropic. This allows for the definition of a *typical* point for the point process, which we usually assume to be located at the origin, useful to obtain results about the process.

A first interesting result about point processes is Campbell's Theorem (named after the work of Campbell, see [33])

**Theorem A.1.** Let $f: S \to \mathbb{R}$ be a non-negative measurable function. Then:

$$ \mathbb{E}\left[\sum_{x \in \Phi} f(x)\right] = \int_{y \in S} f(y) \Lambda(dy). $$

This theorem is an application of Fubini's Theorem to the random sum, after rewriting it as a double integral. Campbell's theorem is a particularly useful too when working with shot-noises over point processes, because it allows for a practical way to express the value of a shot-noise of the point process. It is to note that if $\Phi$ is a homogenous PPP with intensity $\lambda$, we get the following alternative form:

$$ \mathbb{E}\left[\sum_{x \in \Phi} f(x)\right] = \lambda \int_{y \in S} f(y) dy. $$

## Marked point processes

When using point processes, it is common to attach to each point a characteristic (or mark). Let $\mathcal{M}$ be the space of marks and $\mathcal{M}$ its Borel $\sigma$-algebra.
We define a marked point process $\Psi = \{[x_m, y_m]\}$, where the points $\{x_m\}$ form a point process with intensity measure $\Lambda$, called the ground process, and the $y_m$ is the process of marks. We can define an intensity measure for the marked point process for all $B \in \mathcal{S}$ and $M \in \mathcal{M}$ as:

$$ \Upsilon(B \times M) = \mathbb{E}[\Psi(B \times L)]. $$

For all $M \in \mathcal{M}$, the measure $\Upsilon(\cdot \times M)$ is absolutely continuous with regard to the intensity of the ground process. There exists a probability distribution $\mathbb{P}_m$, depending on $x \in S$, on $\mathcal{M}$ such that:

$$ \Upsilon(d(x, m)) = \mathbb{P}_m(x, dm)\Lambda(dx). $$

If the marked point process is stationary, then, $\mathbb{P}_m(x, \cdot) = \mathbb{P}_m(\cdot)$ for all $x \in S$, and we have:

$$ \Upsilon(B \times M) = \lambda \mathbb{P}_m(M)\nu(B). \qquad (A.2) $$

Let $\Psi$ be a stationary marked point process with marks in $(\mathbb{M}, \mathcal{M})$. We can define the following random counting measure on $(S \times \mathbb{M}, \mathcal{S} \otimes \mathcal{M})$:

$$ \mathcal{C}(B \times M) = \int \phi(B) 1_{\{M\}}(B) \mathbb{P}(d\phi). \qquad (A.3) $$

In other words, this measure counts the average number of points in $B$ with marks in $M$. This measure is called the *Campbell measure* of the point process $\Psi$.

## Palm calculus

When studying point processes, an important question is that of the existence of some events, or the value of certain quantities conditioned on a point being present at a given location. For instance, the probability that a point has no neighbor located closer than $r$ is equal to $\mathbb{P}[\Phi(B(x, r)) = 1 \mid x \in \Phi]$. This conditional probability is ambiguous, because the conditioned event has a probability 0. In order to define such conditional events, we introduce the Palm probability distribution. In this section, we will propose two definitions of the Palm probability distribution: a first one relying on the Campbell measure defined above, and a second one based an interpretation from Matthes.

Let us take $\Psi$ a stationary marked point process. From (A.2), we know that the intensity measure of $\Psi$ is a product form involving the Lebesgue measure, the intensity of the ground process of $\Psi$ and a probability distribution $\mathbb{P}_m$. Taking the definition of the Campbell measure from (A.3), we get:

$$ \mathcal{C}(B \times M) = \int \phi(B) 1_{\{M\}}(B) \mathbb{P}(d\phi) = \lambda \int_B \mathbb{P}_m(M)\nu(dx). $$

The probability distribution $\mathbb{P}_m$ is called the Palm probability distribution, and it can be seen as the Radon-Nikodym derivative of the Campbell
measure with regard to the intensity measure of the ground process of $\Psi$. Although this definition is implicit, it underlines its dependency with the Campbell measure of a stationary marked point process. The following definition will give an explicit definition and another interpretation for the Palm probability distribution:

**Definition A.2.** Let $\Phi$ be a stationary point process with a finite, nonzero intensity $\lambda$. The Palm probability distribution (at 0) is the probability distribution defined on $(\mathbb{N}, \mathcal{N})$ as:

$$ \mathbb{P}_{\Phi}^{0}[Y] = \frac{1}{\lambda\nu(B)} \int \sum_{x \in \Phi(B)} 1\{\Phi - x \in Y\} \mathbb{P}(d\Phi), \quad (A.4) $$

where $B$ is an arbitrary Borel set with positive measure. It is to note that this definition does not depend on the choice of set $B$.

The Palm probability distribution introduces a measure which acts as a condition probability given that a point of $\Phi$ is located at 0. To understand why, we can use Matthes' interpretation. Let us associate to each point in $\Phi$ a mark, equal to 0 or 1 whether the shifted point process $\Phi - x$ belongs to $Y$. To come back to the previous example, with $Y = \{\Phi, \Phi(B(0, r) = 0)\}$, a point $x$ will have mark 1 if and only if its nearest neighbor is located at a distance greater than $r$. The process of marks is stationary, and we can extract the subprocess of points with marks 1, which has intensity $\lambda_{\{1\}}$. The mark distribution $M$ is given by:

$$ M(\{1\}) = \frac{\lambda_{\{1\}}}{\lambda}. $$

Let us take a Borel set $B$. The average number of points in $B$ that have mark 1 has mean $\lambda_{\{1\}}\nu(B)$, implying that:

$$ \mathbb{P}_{\Phi}^{0}[Y] = M(\{1\}), $$

and that the choice of $B$ does not play a role in (A.4).

A second definition of the Palm probability distribution can be reached from the Campbell measure (see Section 4.3.4 from [86]). This measure can be proven to be absolutely continuous with respect to the intensity measure of the point process $\Phi$, which allows us to derive the Palm probability distribution.

Using the Palm probability measure, we can define the Palm expectation. This expectation can be interpreted as the average of a variable at the points of the process $\Phi$. This expectation can be expressed in terms of $\Phi$ as:

$$ \mathbb{E}_{\Phi}^{0}[Y] = \frac{1}{\lambda\nu(B)} \mathbb{E} \left[ \sum_{x \in B} Y + x \right], $$

where $Y + x = \{y + x, y \in Y\}$. To illustrate the Palm expectation, let us consider a point process $N$ on the positive real line $\mathbb{R}^+$ with probability $\mathbb{P}$, and let us take a stochastic process (defined in the next chapter) $Z$. If the
point process is ergodic, we can use ergodic theorems (see [27] or Theorem 1.6.4 from [14]) and (A.4), then we can write the Palm expectation and the classical expectation as follows:

$$ \mathbb{E}_N^0[Z(0)] = \lim_{t \to \infty} \frac{1}{N((0, t])} \int_0^t Z(s)N(ds) \qquad \mathbb{E}[Z(0)] = \lim_{t \to \infty} \frac{1}{t} \int_0^t Z(s)ds. $$

The difference between Palm expectation and regular expectation is fundamental in our work, especially when we will study the properties of the stationary regime of the system considered in Part II. The Palm expectation works as a user average, as it only considers the values of Y at locations in the Point process Φ, as the classical expectation can be interpreted as a system average.

An important result linking the Palm probability measure and Poisson point processes is Slivnyak's Theorem (sometimes referred to as Slivnyak-Mecke's theorem):

**Theorem A.2 (Slivnyak).** Let Φ be a point process with intensity measure Λ and distribution P. Then, for all $x \in S$, we have:

$$ \mathbb{P}_x = \mathbb{P} \star \delta_x, $$

where $\star$ is the convolution of distributions and $\delta_x$ is the distribution of the point process consisting of the non-random point x.

In other words, the Palm distribution of Φ with respect to x is the distribution of Φ, with an added point at x, as the convolution of distributions corresponds to the superposition of point processes. This theorem allows to link the Palm probability measure, and subsequently the Palm expectation, with the original distribution of Φ.

A second important result is the inversion formula of Slivnyak and Ryll-Nardzewski:

**Theorem A.3.** Let Y be a non-negative random variable. We have:

$$ \mathbb{E}[Y] = \lambda \mathbb{E}_{\Phi}^{0} \left[ \int_{x \in S} (Y + x) \, dx \right]. \qquad (A.5) $$

This result links the Palm probability measure, the intensity of the point process Φ and the expectation. We can use (A.5) as well as (A.4) to obtain a second equation to express a random sum over a point process. We have, for all $B \subset S$ and for all non-negative function $f: S \to \mathbb{R}$:

$$ \mathbb{E} \left[ \sum_{x \in \Phi(B)} g(x) \right] = \mathbb{E}[\Phi(B)] \mathbb{E}_{\Phi}^{0} [g(0)]. \qquad (A.6) $$

This formula is a mean-value formula and uses the user average interpretation for the Palm probability. To put it in words, to obtain the average
value of the sum of $\Phi$, we multiply the expected value for one user, the typical user, and multiply it by the number of users in the process.

Another result about Palm calculus that we will use in this thesis links the Palm probability distribution of a superposition of point processes and the intensities of these processes:

**Lemma A.4.** Let $\Phi_1$ and $\Phi_2$ be two stationary point processes on $\mathcal{D}$ with respective intensities $\mu_1$ and $\mu_2$, and denote $\Phi = \Phi_1 + \Phi_2$ the superposition of the two processes. Then, we have:

$$ \mathbb{E}_{\Phi}^{0}[F(0, \Phi_0)] = \frac{\mu_1}{\mu_1 + \mu_2} \mathbb{E}_{\Phi_1}^{0}[F(0_1, \Phi_0)] + \frac{\mu_2}{\mu_1 + \mu_2} \mathbb{E}_{\Phi_2}^{0}[F(0_2, \Phi_0)], \quad (A.7) $$

where $0, 0_1$ and $0_2$ respectively denote the typical users of $\Phi$, $\Phi_1$ and $\Phi_2$, and $F$ is a shot-noise of the point processes.

*Proof.* Let $F$ be a shot-noise over $\Phi$. Then, we have, from the definition of the Palm probability distribution:

$$ \mathbb{E} \left[ \sum_{x \in \Phi} F(x, \Phi) \right] = \mathbb{E}_{\Phi}^{0}[F(0, \Phi)] \mathbb{E}[\Phi(\mathcal{D})] = \mathbb{E}_{\Phi}^{0}[F(0, \Phi)] (\mathbb{E}[\Phi_{1}(\mathcal{D})] + \mathbb{E}[\Phi_{2}(\mathcal{D})]). $$

We also know that:

$$ \begin{aligned} \mathbb{E} \left[ \sum_{x \in \Phi} F(x, \Phi) \right] &= \mathbb{E} \left[ \sum_{x \in \Phi_1} F(x, \Phi) \right] + \mathbb{E} \left[ \sum_{x \in \Phi_2} F(x, \Phi) \right] \\ &= \mathbb{E}_{\Phi_1}^{0}[F(0_1, \Phi)] \mathbb{E}[\Phi_1(\mathcal{D})] + \mathbb{E}_{\Phi_2}^{0}[F(0_2, \Phi)] \mathbb{E}[\Phi_2(\mathcal{D})]. \end{aligned} $$

We combine these two equations to conclude the proof. $\square$

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

The last important formula we will mention about Palm calculus in this section is the Papangelou formula, which gives a Radon-Nikodym derivative of the Palm probability measure.

**Theorem A.5.** Let $\Phi$ be a point process and $\mathcal{F}_t$ be a history of $\Phi$ adapted to the flow $\{\theta_t\}$. Assume furthermore that the $\mathcal{F}_t$-predictable structure is adapted to $\{\theta_t\}$. Then, $\Phi$ admits an intensity $\lambda(t)$ if and only if $P^0 \ll P$. In this case, we can chose $\lambda(t)$ of the form $\lambda(t) = (\mu \circ \theta_t)\lambda$, where $\mu$ is the Radon-Nikodym derivative of $P^0$ on $\mathcal{F}_0^-$ with respect to $\mathbb{P}$:

$$ \mu = \left. \frac{d\mathbb{P}^{0}_{\Phi}}{d\mathbb{P}} \right|_{\mathcal{F}_{0}^{-}}. $$

This result allows to reduce the value of an unknown Palm probability measure using possibly known values of the stochastic intensity of other processes. Papangelou's formula, when combined with the inversion formula from (A.5) allows to obtain meaningful results about point processes, as displayed in Part II.
## A.3 Poisson-Voronoi Tessellations

The last notion regarding point process theory we will cover in this chapter is the Poisson-Voronoi tessellation (for more results about random plane tessellations, see [13]). A tessellation is a partition of the plane, i.e. a set $\Theta = \{T_i\}_{i \in I}$ such that:

$$ - \bigcup_{i \in I} T_i = \mathbb{R}^2, $$

$$ - \text{ for all } i \neq j, T_i \cap T_j = \emptyset. $$

Let $\mathcal{T}$ be the set of all possible tessellations of $\mathbb{R}^2$. When considering set of points $S = \bigcup_{i=1}^n \{x_i\}$ with $x_1, \dots, x_n \in S$, one can derive the Voronoi tessellation associated with $S$, $\mathcal{V}(S) = \{c_1, \dots, c_m\}$ where, for all $1 \le i \le m$:

$$ c_i = \{y \in S, \|y - x_i\| < \|y - x_j\|, j \neq i\}. $$

In other words, the Voronoi cell $i$ is the set of points of $S$ that are closer to $x_i$ than to any other point $x_j$, with $j \neq i$. It is to note that the cells $c_i$ may be infinite. Voronoi tessellations (also called Voronoi diagrams or Dirichlet tessellations) have a lot of applications, but the one we will focus on comes from information theory and wireless communications. Let us consider a collection of antennas in a wireless network. If we assume that users in the network always associate with the closest antenna to their location, the result association policy creates a Voronoi tessellation of the plane.

Let $\Phi$ be a homogenous PPP with intensity $\lambda > 0$. As a random collection of points of $S$, we can associate it with a Voronoi tessellation $\mathcal{V}(\Phi)$. This Voronoi diagram is called a *Poisson-Voronoi (PV) tessellation* of $S$.

When studying PV tessellations of the plane, two cells are of interest:

- The typical cell, whose existence comes from the stationarity of the PPP used to generate the diagram, which can be seen as a cell sampled randomly in the tessellation, denoted as $V$,

- The 0-cell, which is the cell containing the origin, denoted as $V_0$.

Multiple results exist about the typical cell of a PV tessellation (see Chapter 9 of [86]). For instance, we know that the average volume of the typical cell is equal to:

$$ \mathbb{E}[|V|] = \frac{1}{\lambda}. \qquad (A.8) $$

An important lemma linking the volumes of $V$ and $V_0$ in PV tessellations is:

**Lemma A.6.** In a PV tessellation, the 0-cell is on average larger than the typical cell:

$$ \mathbb{E}[|V|] \le \mathbb{E}[|V_0|]. $$
*Proof.* Let us define $\alpha_0$ as the point process of vertices of the tessellation, $\alpha_1$, the point process of the edge midpoints and denote as $\lambda_0$ and $\lambda_1$ their respective intensities, and let $\beta = \alpha_0 \oplus \alpha_2 \oplus \Phi$, which has intensity $\lambda_\beta$. The Palm probability of $\beta$ is defined as:

$$ \lambda_{\beta} \mathbb{P}_{\beta}^{0}[Y] = \int_{\theta \in \mathbb{T}} \sum_{\beta([0,1])} 1_{Y}\{\theta - x\} \mathbb{P}(d\theta). \qquad (A.9) $$

Let $f$ be a measurable function on $\mathbb{T}$ and $T = \{\theta \in \mathbb{T}, 0 \in \Phi\}$ be the set of tessellations for which the origin is the center of a cell. Let us define:

$$ \mathbb{E}_2[f(\Theta)] = \frac{\lambda_{\beta}}{\lambda} \int_{T_2} f(\theta) \mathbb{P}(d\theta). $$

We have:

$$ \lambda \mathbb{E}_2 \left[ \int_{V_0} f(\Theta - x) dx \right] = \lambda_{\beta} \int_{\mathbb{T}} \int_{\mathbb{R}^2} 1_{T_2}(\theta) 1_{V_0(\theta)}(-x) f(\theta + x) \mathbb{P}_{\beta}^{0}(d\theta). $$

We use the mass transportation principle. We have, for any measurable function $v$:

$$ \lambda_{\beta} \int_{T} \int_{\mathbb{R}^2} v(x, \theta) dx \mathbb{P}_{\beta}^{0}(d\theta) = \int_{T} \sum_{\beta} v(x, \theta - x) \mathbb{P}(d\theta). $$

With $v(x, \theta) = 1_{T_2}(\theta) 1_{V_0(\theta)}(-x) f(\theta + x)$ we get:

$$ \lambda \mathbb{E}_2 \left[ \int_{V_0} f(\Theta - x) dx \right] = \int_{\mathbb{T}} \sum_{x \in \beta} 1_{T_2}(\theta - x) 1_{V_0(\theta - x)}(-x) f(\theta) \mathbb{P}(d\theta). $$

We know that $1_{T_2}(\theta - x) 1_{V_0(\theta - x)}(-x) = 1$ if and only if $-x$ belongs in the 0-cell of $\theta - x$. In each $\theta$, there is only one such $x$, giving:

$$ \lambda \mathbb{E}_2 \left[ \int_{V_0} f(\Theta - x) dx \right] = \int_{\mathbb{T}} f(\theta) \mathbb{P}(d\theta) = \mathbb{E}[f(\Theta)]. \qquad (A.10) $$

Using this result with $f = 1$ yields:

$$ \lambda E_2 \left[ \int_{V_0} dx \right] = \lambda \mathbb{E}[|V|] = 1, $$

which is the volume of the typical cell of a PV tessellation. Using (A.10) with $f = 1/|V_0|$ gives:

$$ \mathbb{E} \left[ \frac{1}{|V_0|} \right] = \lambda = \frac{1}{\mathbb{E}[|V|]}. $$

Using Jensen's inequality gives:
$$ \mathbb{E}[|V|] \leq \mathbb{E}[|V_0|], $$

which is the intended result and concludes the proof. □

Another interesting result about the link between the 0-cell and the typical cell in a PV tessellation comes from [69]:

**Lemma A.7 (Mecke, 1998).** For any measurable, non-negative function g, we have:

$$ \mathbb{E}[g(|V_0|)] = \frac{1}{\mathbb{E}[|V|]} \mathbb{E}[|V|g(|V|)]. $$

Lemma A.7 gives an interesting link when comparing quantities linked with the volume of the 0-cell and the typical cell, which is usually easier to compute - for instance in the case of homogenous PPP.

Although the PV tessellation is one of the simplest instances of random tessellations of the plane, there are very few analytical results about the size distribution of PV cells, but several heuristics and approximation for this distribution. Among the different existing heuristics, [43] provides a simple, one-parameter yet performant heuristic to approximate the size distribution of PV cells:

$$ p(y) = A y^{5/2} \exp \left( -\frac{7}{2} y \right), \qquad (A.11) $$

where $A = \frac{(\frac{7}{2})^{\frac{7}{2}}}{\Gamma(\frac{7}{2})}$ is a normalization constant. Using this heuristic allows us to derive an approximation for the size of the 0-cell. We have:

$$ \mathbb{E}[|V_0|] \approx \int_{0}^{+\infty} xp(x)dx = \frac{9}{7} \frac{1}{\lambda} \approx 1.28 \mathbb{E}[|V|]. $$

This numerical approximation holds when simulating the average size of the typical cell against the size of the 0-cell of a PV tessellation. Figure A.1 displays the size of the zero-cell (in dashed line) against the size of the typical cell (in plain line) for several values of $\lambda$.

In wireless communications, an interesting problem is to study the number of users connected to the same antenna (for reasons exposed in the next chapter). To do so, if we assume that antennas and mobile users are distributed according to two PPP of respective intensities $\lambda$ and $\nu$, a naive heuristic to obtain the average number of users in the zero-cell of the network is given by:

$$ \mathbb{E}\left[\frac{1}{Z}\right] \approx \frac{1}{1 + 1.28 \frac{\lambda}{\nu}}. \qquad (A.12) $$

Though, as we will see late in this section, this approximation performs poorly to estimate the inverse load in a network. This leads us to defining the following heuristic, using (A.11):
<figure><img src="image_10.png" /><figcaption>Figure A.1: Size of the 0-cell (dashed line) against the size of a typical cell (plain line) for a homogeneous PV tessellation</figcaption></figure>

**Lemma A.8.** Let $\Phi$ and $\Psi$ be two PPPs of respective intensities $\lambda$ and $\nu$. Let $V_0$ be the 0-cell of the PV tessellation associated with the process $\Phi$ and let $Z = 1 + |\Psi \cap V_0|$. Finally, let $L: x \mapsto x \left(1 - \left(\frac{1}{1+\frac{2}{\tau x}}\right)^{7/2}\right)$. The moment of order $-1$ of $Z$ can be approximated by:

$$ \mathbb{E}\left[\frac{1}{Z}\right] \approx L\left(\frac{\lambda}{\nu}\right). $$

*Proof.* The number of users in $\Psi$ in the zero-cell of the Poisson-Voronoi tessellation associated with $\Phi$ is distributed with a Poisson law of parameter $\nu|V_0|$, where $V_0$ is the 0-cell of the tessellation. We have:

$$ \begin{aligned} \mathbb{E}\left[\frac{1}{1 + |\Psi(V_0)|}\right] &= \mathbb{E}\left[\sum_{N=0}^{\infty} \frac{1}{1+N} \mathbb{P}[|\Psi(V_0)| = N]\right] \\ &= \mathbb{E}\left[\sum_{N=0}^{\infty} \frac{1}{1+N} \frac{(\nu|V_0|)^N}{N!} e^{-\nu|V_0|}\right] \\ &= \mathbb{E}\left[\frac{1}{\nu|V_0|} (1 - e^{-\nu|V_0|})\right]. \end{aligned} $$

We can use Lemma A.7, with $g(V) = \frac{1}{V}(1 - e^{-\nu V})$ and $\mathbb{E}[|V|] = 1/\lambda$.

$$ \mathbb{E}\left[\frac{1}{1 + |\Psi(V_0)|}\right] = \frac{\lambda}{\nu} \left(1 - \mathbb{E}\left[e^{-\frac{\nu}{\lambda}|\lambda|V|}\right]\right). \quad (A.13) $$

Equation A.13 is the Laplace transform of the random variable $\lambda|V|$ at $s = \frac{\nu}{\lambda}$. Using (A.11) to evaluate the Laplace transform, we get the following approximation:
<figure><img src="image_11.png" /><figcaption>Figure A.2: Performance of the heuristic of Lemma A.8 (plain line) against the simulation (circles) and the naive heuristic (dashed-dotted line) for different values of $\lambda/\nu$.</figcaption></figure>

$$ \mathbb{E} [e^{-\nu|V|}] = \left(1 + \frac{2\nu}{7\lambda}\right)^{-7/2}. $$

Substituting this result in (A.13) leads to the result. □

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

It is to note that this result only relies on the ratio of intensities $\frac{\lambda}{\nu}$, and not on the value of the intensities themselves. We can use numerical simulation to assert the performance of this heuristic, by simulating the inverse load in a homogenous PPP of intensity $\lambda$ and a homogenous PPP of intensity $\nu$. Figure A.2 shows the comparison between the simulation results (in circles) with the heuristic from Lemma A.8 (in dashed line) and the naive heuristic from (A.12) (in dashed-dotted line).

The reason for this numerical setup will be explained in Part III of this manuscript, as having a performant approximation for the average inverse load in the typical cell of the PV tessellation will allow us to define heuristics to maximize some network metrics.

## A.4 Markov Jump Processes

In this section, we will present some basic results about Markov processes. We will consider stochastic processes as measurable maps from $(\Omega, \mathcal{F})$ to $(S, \mathcal{S})$. The main category of stochastic processes we consider in this work are Markov Jump Processes, which is an extension of Markov chain in continuous time. Intuitively, a Markov jump process is a piece-wise constant stochastic process with exponentially distributed interevent times.
Firstly, we define a Markov kernel (or probability kernel) as any map $P: S \times S \to [0, 1]$ such that:

1. The map $x \mapsto P(x, A)$ is measurable for all $A \in S$.

2. For every $x \in S$, the map $A \mapsto P(x, A)$ is a probability measure on $(S, S)$.

If, instead of having mass 1 on $S$, we have $P(x, S) \le 1$ for all $x \in S$, we say that $P$ is a sub-Markov kernel. We then say that the family of Markov kernels $(P_t)_{t \ge 0}$ is a transition function on $S$ if:

1. For all $x \in S$, $A \in S$ and $s, t \ge 0$, $(P_t)$ satisfies the Chapman-Kolmogorov equations:

$$ \int_{y \in S} P_s(y, A)P_t(x, dy) = P_{s+t}(x, A). $$

2. $P_0(x, A) = \delta_x(A)$ for all $x \in E$ and $A \in S$.

Now, let $(X_t)_{t \ge 0}$ be a stochastic process with state space $S$ and $\mathcal{F}_t = \sigma\{X_s, s \le t\}$ be its canonical filtration. We say that $(X_t)_{t \ge 0}$ satisfies the simple Markov property with respect to $\mathcal{F}$ if for all $t, h \ge 0$, we have:

$$ \mathbb{P}[X_{t+h} \in A | X_s, s \le t] = \mathbb{P}[X_{t+h} \in A | X_t] - \text{a.s.} \quad (A.14) $$

We then say that $(X_t)$ has transition function $(P_t)$ if and only if:

$$ \mathbb{P}[X_{t+h} \in A | X_t] = P_h(X_t \in A). \quad (A.15) $$

Equations (A.14) and (A.15) can be seen as the extension of the well-known Markov property for discrete-time processes to continuous time.

Now, we can move on to defining jump processes. Let $\Pi$ be a Markov Kernel on $(S, S)$, and let $\lambda: S \to [0, \infty)$ be a non-negative, measurable function. We define the *jump kernel* as:

$$ K(x, dy) = \lambda(x)\Pi(x, dy). $$

We can then move to defining the *jump-hold construction* of the Markov jump process: let us take a law $\mu$, random variables $(Z_n)_{n \ge 0}$ with values in $S$ and $(t_k)_{k \ge 1}$ positive random variables (possibly infinite) such that:

- $(Z_n)_{n \ge 0}$ is a Markov chain with transition kernel $\Pi$ and initial law $Z_0 \sim \mu$. $(Z_n)$ is called the *embedded jump chain* or embedded Markov chain.

- The variables $(t_k)_{k \ge 1}$ conditional on $(Z_n)$ are independent exponential variables are exponential random variables with respective parameters $\lambda(Z_{k-1})$. The $(t_k)_{k \ge 1}$ are called the *holding times* of the process.
This jump-hold construction displays the two main properties of a Markov jump process: the existence of an embedded Markov chain, representing the jumps in the process, and a collection of exponentially distributed family of holding times depending at each step only on the previous jump. This exponential distribution of holding times and the fact that jumps are following a Markovian evolution justify the name *Markov jump process*, as each jump and the time to wait until the next jump are entirely determined by the present state of the system.

Another determination of Markov jump process comes from its transition function ($P_t$). Let us take $(x_n)_{n \ge 0}$ be elements of $S$ and take $0 < t_1 < \dots < t_k < t$ be real numbers. Taking a jump process with jump kernel $K(x, dy) = \lambda(x)\Pi(x, dy)$, a path from $x_0$ to $x_k$ in time $t$ in $k$ jumps at times $(t_i)_{1 \le i \le k}$ and locations $x_1, \dots, x_{k-1}$ can be associated with weight:

$$ \lambda(x_0)\Pi(x_0, dx_1)e^{-\lambda(x_1)(t_2-t_1)}\Pi(x_1, dx_2) \dots e^{-\lambda(x_{k-1})(t_k-t_{k-1})}. $$

Let us define the following sequence, for $n \ge 0$, $t > 0$, $x_0 \in S$ and $A \in S$:

$$ \begin{aligned} P_t^n(x_0, A) = & e^{-\lambda(x_0)t} \delta_{x_0}(A) + \sum_{k=1}^{n} \int_{0<t_1<\dots<t_k \le t} dt_1 \dots dt_k \\ & \int_{S^k} \Pi(x_0, dx_1) \dots \Pi(x_{k-1}, dx_k) \mathbb{1}_{x_k \in A} \\ & \times \lambda(x_0) e^{-\lambda(x_0)(t_1)} \lambda(x_1) e^{-\lambda(x_1)(t_2-t_1)} \dots \lambda(x_{k-1}) \\ & e^{-\lambda(x_{k-1})(t_k-t_{k-1})} + e^{-\lambda(x_0)t} \mathbb{1}_{x_0 \in A}. \end{aligned} $$

We can prove that this family of functions is increasing, and that the point-wise limit $P_t^*(x, A) = \lim_{n \to \infty} P_t^n(x, A)$ exists for all $x \in S$ and $A \in S$. We can also prove that $P_t^*$ is a sub-Markov kernel.

The main result about the family $P_t^*$ is the following:

**Theorem A.9.** The $P_t^*$ is the solution to the Kolmogorov backward equation:

$$ \frac{\partial}{\partial t} P_t^*(x, A) = \int_S K(x, dy) (P_t^*(y, A) - P_t^*(x, A)). \quad (A.16) $$

Every other sub-Markov kernel solution to (A.16) is greater than $P_t^*$. Finally, if $(P_t^*)$ is a family of Markov kernels, i.e. if $P_t^*(x, S) = 1$ for all $x \in S$, then it is the unique solution to (A.16).

This Theorem gives us a second characterization of Markov jump processes, through its transition functions. If the transition function of the process is a Markov kernel, then the only solution of (A.16) fully characterizes the jump process, and it uniquely characterizes the transition kernel of the jump process.

## Spatial birth-and-death process

Among the Markov jump processes, birth-and-death processes have been used often to model physical phenomena, especially in biology.
We say that a Markov jump process is a birth-and-death process if the only possible transitions from a state with $n$ elements are towards a state with either $n-1$ or $n+1$ elements. For instance, the continuous-time M/M/1 queue with an arrival parameter $\lambda$ and a departure parameter $\nu$ is a simple example of a birth-and-death process.

Though there exists a large variety of birth-and-death processes, the category we will focus in this work is the spatial birth-and-death (SBD) process. This class of birth-and-death process form a Markov jump process on point process of a region (or the entirety) of $\mathbb{R}^d$, and where arrivals and departures depend on the spatial configuration of the system.

These processes have been introduced for the first time by Preston in [75]. In a spatial birth-and-death process, the law of the jump kernel $K$ depends on the spatial configuration of the system at a given time.

To define such a process, assume that at a given time $t$, points are located at $X_1, \dots, X_n$ in $S$. Then, we need two functions $b$ and $d$ such that:

- $b$ is a finite measure on $(S, \mathcal{S})$ such that the probability that a point arrives in $A \in \mathcal{S}$ between times $t$ and $t+s$ is equal to $b(X_1, \dots, X_n; A)s$.

- $d: S \to \mathbb{R}^+$ is a measurable function such that the probability that a point $x$ leaves the system between times $t$ and $t+s$ is equal to $d(\{X_i, X_i \neq X\}; X)s$.

These two functions, called respectively the *birth* and the *death* functions of the process, define a Markov jump process, with state space $M(S)$, the state of counting measures on $S$. We can consider that at each time, the elements of the Markov jump process form a point process on $S$.

## A.5 Stability theory

The final section we will tackle in this chapter comes from the study of the stability of Markov chains. The results and definitions presented here will be used throughout Part I of this manuscript.

The study of the behavior of Markov chains is an interesting question, as it allows to obtain a lot of asymptotical results about the processes we study. In this section, $\mathbf{X}$ denotes a Markov chain with state space $\mathbb{N}$ (these results can be extended to any Markov chain with countable state space).

The first notion we define is the *irreducibility* of a Markov chain:

**Definition A.3.** The Markov chain $\mathbf{X}$ is $\phi$-irreducible if and only if there exists a measure $\phi$ on $\mathcal{B}$ such that whenever $\phi(A) > 0$ for $A \in \mathcal{B}$, we have $\Pi(x, A) > 0$ for all $x \in S$. A Markov jump process is irreducible if and only if its embedded Markov chain is $\phi$-irreducible for some measure $\phi$.

In other words, a chain is $\phi$-irreducible whenever we can find an irreducibility measure $\phi$ such that any event with positive measure $\phi$ can be reached from any given state. A consequence is that a $\phi$-irreducible Markov
chain has only one communication class, i.e., that any state can be reached from any given starting point with non-zero probability.

We can now move to defining notions about stability of the Markov chain. The time of first return to state $i \in \mathbb{N}$ is the first moment at which the chain hits state $i$ given it started in this state:

$$T_i = \inf\{k \ge 1, X_k = i | X_0 = i\}.$$

A state $i$ is said to be transient if $\mathbb{P}[T_i = \infty] > 0$ and recurrent otherwise. Moreover, we say that a state is positive recurrent if the state $i$ is recurrent and we have:

$$\lim_{n \to \infty} \frac{1}{n} \sum_{0 \le k \le n-1} 1_{\{X_k=i\}} = \pi_i > 0. \qquad (A.17)$$

In other words, a state is recurrent if its return time is almost surely finite, and transient otherwise, i.e., if we have a non-zero probability to never return to it. Finally, a state is positive recurrent if we return to it a significant portion of the time. An equivalent definition to A.17 for positive recurrence is that $\mathbb{E}[T_i|X_0 = i] < \infty$, i.e., that the average return time to state $i$ is always finite. If this mean return time to state $i$ is infinite but it is recurrent, we say that it is null recurrent.

We say that the Markov chain is recurrent if all its states are recurrent, and transient if all its states are transient. Finally, we say that a chain is positive recurrent if all its states are positive recurrent.

Positive recurrence for Markov chain is an important property as it gives an easy characterization for the stationary distribution of the Markov chain: let us define $\pi_i$ as in (A.17). Then, $(\pi_i)_{i \in \mathbb{N}}$ forms a probability distribution on $\mathbb{N}$, which is the stationary distribution for the chain $\mathbf{X}$. It also allows for the use of large deviation inequalities, and gives bounds for the speed of convergence to the stationary distribution of the chain.

**Definition A.4.** We say that a Markov chain is stable if it is positive recurrent. Otherwise, it is unstable.

From the definitions of positive recurrence, transience and null recurrence, this definition of stability ensures that $\mathbb{P}$-almost surely, the value of $X(t)$ is finite.

A first interesting result linking $\phi$-irreducibility and state classification of a Markov chain comes from Theorem 8.2.5 from [70]:

**Theorem A.10 (Meyn, Tweedie).** Every $\phi$-irreducible Markov chain is either transient or recurrent.

This result about the classification of Markov chain will help us in obtaining a description of recurrence and transience of Markov chains. The main result about stability for Markov chains comes from Foster (see [44]):
**Theorem A.11** (Foster, 1953). Let $\mathbf{Y}$ be a Markov process with countable state space $E$. Let $C \subset E$ be a compact subset of the state space, $\alpha, \beta > 0$ be two constants and $V: \mathbb{N}^N \to \mathbb{R}^+$ be a positive definite function. Then, if:

$$ \Delta V(x) = \mathbb{E}[V(\mathbf{Y}(T_1)) - V(\mathbf{Y}(0)) | \mathbf{Y}(0) = x] \le \beta \mathbb{1}_{x \in C} - \alpha \mathbb{1}_{x \notin C}, $$

where $t_1$ is the time of the first event in the system, then, $\mathbf{Y}$ is positive recurrent.

In other words, if we find a positive definite function $V$, i.e., such that $V(x) > 0$ for all $x$, and $V(x) = 0$ if and only if $x = 0$ and a compact subset of the state space such that this subset acts as an attractor for the dynamics, then we obtain the positive recurrence for the Markov chain. Sometimes, this theorem is referred to as the Foster-Lyapunov theorem, as it is reminiscent of arguments used in stability theory to find attractors for dynamical systems, and the condition stated resembles the negative Lie derivative requirement to obtain stable points or intervals for systems of ODEs.

This class of arguments, called *drift arguments*, as we need to bound the average increase in the chain over one (or several) jumps. Unfortunately, the main limitation for the use of this Theorem is the need to find an appropriate Lyapunov for the Markov chain we study. Outside of certain go-to forms for certain classes of dynamics, there are not particular method to find an appropriate Lyapunov function for Markov jump processes. Thus, we will rely on another method to obtain stability for Markov processes: fluid limits.

## Fluid limits

Fluid limits have been introduced by Rybko and Stolyar in [78] to study the stability for queues of networks using a different method than the classical arguments stated above. The main idea behind fluid scaling is the following: let $\mathbf{X}$ be a Markov chain, and we define the *fluid-scaled* sequence of processes $\mathbf{X}^n$ as:

$$ \begin{cases} \mathbf{X}^n(0) = n \\ \mathbf{X}^n(t) = \frac{1}{n}\mathbf{X}(nt). \end{cases} $$

The goal is to study the limit as $n$ goes to infinity of the fluid-scaled sequence $\mathbf{X}^n$. If such a function exists, we denote its value at time $t$ as $\bar{x}(t)$. If the function $\bar{x}$ is continuous, we call it a *fluid limit* for the stochastic process $\mathbf{X}$. Usually, when studying this limit, we obtain function equations defining the function $\bar{x}$, which we call a *fluid model* for the dynamics.

We can generalize fluid scaling to any positive increasing sequence $(a_n)_{n \ge 0}$ with $\lim_{n \to \infty} a_n = \infty$. The sequence of fluid-scaled processes with $(a_n)_{n \ge 0}$ is defined as:

$$ \begin{cases} \mathbf{X}^n(0) = a_n \\ \mathbf{X}^n(t) = \frac{1}{a_n}\mathbf{X}(a_n t). \end{cases} $$
It is to note that multiple fluid limits can exists for a given process, as it may depend on the sequence used for scaling. For the same scaling sequence, fluid limits may not be unique. For instance, [77] show an example of a queuing system where fluid limits are random.

**Definition A.5.** We say that a fluid limit is stable if and only if there exists $\delta > 0$ such that $\bar{x}(\cdot + \delta) = 0$.

The main interest for fluid limits is its ties with positive recurrence in Markov chains. The main result linking the two comes from Theorem 4.4 from [37]:

**Theorem A.12 (Dai, 1995).** If the fluid limit $\bar{x}$ is stable, the Markov chain $X$ is positive recurrent.

An example to display how fluid limits work is the M/M/1 queue with an arrival rate $\lambda > 0$ and a departure rate $\mu > 0$ and an initial condition $X(0)$. The stochastic recurrence for the queue length process $X(t)$ is:

$$X(t) = X(0) + \mathcal{A}(\lambda t) - \mathcal{D}(\mu t),$$

where $\mathcal{A}$ and $\mathcal{D}$ are PPP with intensity 1 (see [31]). Let us take the fluid scaled process with $a_n = n$. We have:

$$X^n(t) = \frac{1}{n}X(0) + \frac{1}{n}\mathcal{A}(\lambda nt) - \frac{1}{n}\mathcal{D}(\mu nt).$$

Using the strong law of large numbers, we get the equation ruling the evolution of the fluid limit $\bar{x}$, called the *fluid model*:

$$\bar{x}(t) = 1 + (\lambda - \mu)t. \qquad (A.18)$$

We can immediately see that the condition for the fluid limit to reach 0 after a certain point is $\lambda < \mu$, which is the well-known stability condition for the M/M/1 queue.

Among the results giving stability for fluid models, an interesting one that we will use later in this work comes from [84]:

**Lemma A.13 (Shneer, Stolyar, 2020).** Let $X = (X_i)_{0 \le i \le N-1}$ be a Markov chain with the following fluid model:

$$\begin{cases} \bar{x}'_i(t) = \lambda_i - \psi_i(\bar{x}(t)) \\ \bar{x}(0) = x_0, \end{cases}$$

where $x_0 \in \mathbb{R}^N$ ($\psi_i$) are function such that:

1. $\psi_i$ is non-increasing in $x_j$, for all $j \neq i$,

2. for all $i$, $\psi_i$ is 0-homogenous, i.e., $\psi_i(st) = \psi_i(t)$ for all $s, t > 0$,

3. $\sum_i x_i(t) > 0$ implies that $x_i(\tau) > 0$ for $\tau$ sufficiently close to $t$.
Let $\mathcal{C} = \{\mathbf{z} \in \mathbb{R}_+^N, z_i \le \psi_i(\mathbf{p}) \text{ for some } \mathbf{p} \in \mathbb{R}_+^N\}$. Assume there exists $\nu \in \mathcal{C}$ such that for all $i$, $\lambda_i \le \nu_i$. Then, for all $0 < \delta < K < \infty$, there exists $T > 0$ such that for any trajectory with $\|\bar{\mathbf{x}}(0)\| = K$, we have:

$$ \|\bar{\mathbf{x}}(T)\| < \delta. $$

To obtain stability for the chain $\mathbf{X}$, we use this results with $K = 1$ and $\delta = 1 - \varepsilon$, for some $\varepsilon > 0$. This result linking positive recurrence and fluid limits will be used in Part I to reach a stability condition for the system.

Lemma A.13 allows to give bounds for the stability region, but in certain cases, these bounds are not tight, and it allows only for a partial description of the stability region. Unfortunately, the reciprocal for this result is not true: Bramson, in [30], showed that there exists transient queuing networks for which the fluid limit match the condition of Theorem A.12.

To obtain a reciprocal, we need to loosen the condition, and define the notion of weak instability, from [38]. We say that a fluid model is weakly unstable if for any fluid limit $\bar{\mathbf{x}}$ such that $\bar{\mathbf{x}}(0) = 0$, there exists $\delta > 0$ such that $\bar{\mathbf{x}}(\delta) > 0$. The reciprocal for Theorem A.12 becomes:

**Theorem A.14.** If the fluid model is weakly unstable, then we have:

$$ \mathbb{P}\left[\lim_{t \to \infty} \mathbf{X}(t) = +\infty\right] = 1. $$

This result implies that the Markov chain $\mathbf{X}$ is transient, which gives a description of regions where the chain $\mathbf{X}$ is positive recurrent or transient. Coming back to the M/M/1 queue, with an empty initial condition, (A.18) becomes:

$$ \bar{\mathbf{x}}(t) = (\lambda - \mu)t. $$

Immediately, we can see that if $\lambda > \mu$, the fluid limit $\bar{\mathbf{x}}(t)$ is strictly increasing and the fluid model is weakly unstable, which gives us instability for the M/M/1 queue. This also shows a limitation of fluid limits, where it is not uncommon to obtain a stability and an instability region that do not cover all cases: here, we cannot answer the case $\lambda = \mu$ (which is null recurrent).

This concludes this chapter about the mathematical foundations of the work presented in this manuscript. The next chapter will cover the basics of wireless communications we modelled in this work.
<figure><img src="image_13.png" /></figure>
# WIRELESS NETWORKS

This chapter focuses on wireless networks as well as some elements of information theory that will be used in this thesis. Results shown here as well as more details can be found in Volume I and II of [12], in [36] and in [8].

In our work, we will consider two types of network architectures: cellular networks and device-to-device networks. These two network architectures have a radically different conception of communications: in cellular network, the network is vertical, i.e, there is one particular user, the base station, to which all users of the network connect, either to transmit or receive data. On the opposite, device-to-device networks are horizontal in the sense that each transmitter has its own dedicated receiver, and all receiver-transmitter pairs in the networks are treated equally in the network. This chapter will cover the basics of the models we use for wireless communications.

## B.1 Channel capacity and Shannon rate

When working with communication channels, a most important quantity to consider is the rate at which information is transmitted in the channel. Let us consider a simple communication system.

Figure B.1 shows a basic communication system, where $W$ is the message to transmit, $X$ is the channel input, $Y$ is the channel output and $\hat{W}$ is the transmitted message, which is an estimate of the original message $W$, and $f$ (resp. $g$) are the encoding (resp. decoding) functions. $p(y|x)$ represents the noisy channel as a conditional probability distribution, i.e. the probability

<figure><img src="image_14.png" /><figcaption>Figure B.1: Basic communication system</figcaption></figure>
that the output is $Y$ given the input is message $X$. This condition probability is a property of the communication channel. In this setup, the distribution of input messages $p_X$ fully determines the joint distribution $p_{(X,Y)}(x) = p(y|x)p_X(x)$

We also define the mutual information of the two random variables $X$ and $Y$ as:

$$I(X; Y) = D_{\mathrm{KL}}(p_{(X,Y)}||p_X \otimes p_Y),$$

where $D_{\mathrm{KL}}$ is the Kullback-Leibler divergence (see [59]). Then, the channel capacity is defined as:

$$C = \sup_{p_X} I(X; Y).$$

In other words, the channel capacity is the maximum amount of information that can be transmitted over the communication channel with arbitrary precision. The Shannon-Hartley theorem (see [83]) gives a useful expression for channel capacity in a particular case:

**Theorem B.1.** The channel capacity $C$ in a communication channel subject to additive white Gaussian noise is equal to:

$$C = B \log_2 \left( 1 + \frac{S}{N} \right), \qquad (\mathrm{B.1})$$

where $B$ is the bandwidth of the channel, $S$ is the average signal power received over the bandwidth and $N$ is the average power of the noise and interference.

In the presence of interference, and if we consider interference as noise, the quantity $S/N$ has to be replaced by $S/(N+I)$, where $N$ is the noise power and $I$ is the interference power. This ratio is called the *Signal-to-Noise and Interference Ratio* (SINR). In the rest of this work, the channel capacity from (B.1) will be noted as $\mathcal{R}$ and referred to as the *Shannon rate*. In a more general case, when working with non-deterministic networks, the result from Theorem B.1 becomes:

$$\mathcal{R} = B\mathbb{E}[\log(1 + \mathrm{SINR})]. \qquad (\mathrm{B.2})$$

## Cover and outage probabilities

When considering SINR and channel capacity, an interesting property to study is coverage. Let us set $T > 0$ a SINR threshold. We say that a point $x$ is in coverage if and only if $\mathrm{SINR}(x) > \beta$. If a point is not in coverage, then we say it is in outage.

The study of SINR cells and spatial coverage is a subject of study in itself (see Chapter 5 of [12]). In our work, we will consider the coverage probability, defined as:
**Definition B.1.** Let $\beta > 0$. The coverage probability at point $x$ is defined as:

$$p_c(x, \beta) = \mathbb{P}[\text{SINR}(x) > \beta].$$

In this definition of the coverage probability, the SINR is a random variable depending on the distribution of the base station process. Using the definition of the channel capacity, we can rewrite (B.1) using the coverage probability:

**Lemma B.2.** Given the coverage probability, the channel capacity can be rewritten as:

$$\mathcal{R}(x) = \int_{0}^{\infty} \frac{p_c(x, \beta)}{\beta + 1} d\beta.$$

This result is obtained by using an integration by parts on (B.2). Lemma B.2 gives a useful link tying the coverage probability and the Shannon rate, which we will exploit later in this manuscript.

## B.2 Wireless networks and stochastic geometry

Stochastic geometry has emerged to become a powerful tool to model wireless networks (see [6], [48]), as it provides a flexible framework to define and compute network metrics and encapsulate the main features of wireless networks.

In this work, we will assume that base stations (BSs) will be spatially distributed according a homogeneous Poisson point process $\Phi$ with intensity $\lambda > 0$. Mobile users (MUs) are distributed according an independent homogeneous PPP $\Phi_u$ with intensity $\lambda_u$. Finally, we will assume that MUs move in a straight line with velocities i.i.d. with distribution $f$. The stationary and isotropy of the MU process is allows us to define the typical MU located at the origin.

When it comes to access policies, we assume that the networks will always be in open access, i.e., that all users can connect without restrictions and that there are no communication abandonments. We also assume that a MU will always connect to its closest antenna. This association policy defines a Voronoi diagram based on the BS process. Because they are distributed according a homogeneous PPP, it forms a PV tessellation and we can use results presented in A.3.

When modeling wireless communications, we have to take into account the attenuation of the signal among the link between the receiver and the transmitter. To take into account signal decay over the distance and electromagnetic and atmospheric phenomena, we use an omnidirectional path-loss (OPL) function to model attenuation over the link. For wireless network, we will use an OPL3 function (see Section 23.1.2 from [12]) of the form
$\ell(x) = Kx^{-\alpha}$, with $K > 0$ being a constant depending on the frequency used for the communications and $\alpha > 2$ denotes the attenuation of the signal. It is to note that this function has a singularity when $x$ becomes close to 0.

When studying systems at a smaller scale, like D2D networks, to avoid singularities, we use the OPL2 function of the form $\ell(x) = \frac{K}{1+x^\alpha}$ for D2D networks, which is defined when $x = 0$. As for the path-loss exponent, for numerical applications, we will use $\alpha = 4$, which gives a satisfactory approximation of the signal attenuation in an environment with a lot of reflection and obstacles, like a city.

A phenomenon we need to take into account is fading, which models the effects of propagation of the signal in the environment. In our models, we will only consider the case of Rayleigh fading, which gives a good approximation of signal fading in urban areas (see [34]). Under Rayleigh fading, the signal over the communication channel is multiplied by a fading coefficient $h_x$ distributed exponentially with mean $\frac{1}{\mu}$, where $\mu$ denotes the fading power in the system. When not specified, we will take a fading $\mu = 1$ in our calculations.

The interference in the network is a shot-noise of the base station process. Let us assume that the base stations in the system form a point process $\Phi$. We can define the interference seen from transmissions at location $x$ as:

$$I(x, \Phi) = \sum_{y \in \Phi \setminus \{x\}} Kh_y \ell(\|x - y\|).$$

Assuming that the transmitter associated with the receiver located at $x$ is at a distance $r_x$ and that the network has a noise power density $N_0$, we have:

$$\text{SINR}(x) = \frac{Kh_x \ell(r_x)}{N_0 + I(x, \Phi)}.$$

As mentioned earlier, this rather simple SINR model encapsulates most of the key phenomena needed to model wireless and radio transmissions, and it proves to be interesting as a lot of closed-form result can be obtained from it, as shown in Part III of this manuscript.

## B.3 Beam-based communications

In the latest generation of wireless communications, antennas use directional beamforming to communicate with MUs (see [49]) as this technology has emerged as a mean to increase the spectral efficiency of millimeter-wave networks ([11]). The main idea behind this technology is for antennas to use narrow and directive beams tracking users in order to increase the bandwidth in the network. The beam in which the user is located is called its reference beam. To ensure that the user is always located in its reference beam, the network sends a synchronization signal at a fixed frequency $\tau$, called a synchronization signal burst (SSB).
<figure><img src="image_15.png" /><figcaption>Figure B.2: Sector approximation of the antenna gains</figcaption></figure>

Directional beam forming and sector approximation

In this setup, each antenna has a fixed number $n$ of beams, usually, we have $N = 2^k$ for some $k > 0$. Under directional beamforming, the MU is tracked by the BS so that communication operate through the reference beam. We model beam-based communications using a sector approximation of the antenna. The beam tracking the user thus has an antenna gain equal to the gain of the main lobe of the antenna. To model the residual power emitted around the antenna, we add side lobes to the antenna, which users outside of the beamwidth of the main lobe see. Each beam has the same angular width $\psi = \frac{2\pi}{n}$. The main lobe is restricted to the beamwidth. The antenna gain in the macro tier is assumed to be:

$$G^M(\theta) = \begin{cases} G_m^M & \text{if } |\theta| \le \psi_M/2 \\ G_s^M & \text{else,} \end{cases}$$

where $G_m^M$ is the main lobe gain and $G_s^M$ is the side lobe gain. The probability that an MU lies within the main lobe of an interfering BS is $p_m = \frac{\psi}{2\pi} = \frac{1}{n}$. The BS antenna gain $g_x$, seen at the origin of the system, of an interfering BS located at $x$ in the macro tier is given by:

$$g_{M,x} = \begin{cases} G_m^M & \text{w.p. } p_{M,m} \\ G_s^M & \text{w.p. } 1 - p_{M,m}. \end{cases} \qquad (B.3)$$

This approximation is a first order approximation of the distribution of antenna gains in space. More complex models exist to obtain more precise calculations, but they fall out of the scope of our work.
Overhead and beam reselections

When an MU crosses the boundary of the cell of its serving BS, it needs to select a new antenna to connect and it performs a *handover* to the next BS. Similarly, when an MU moves outside of its reference beam, it needs to select the next beam of its serving BS, performing a *beam reselection*.

Using stochastic geometry tools, for an MU moving in a straight line at velocity $v$, we can compute the time intensity $\nu_c$ of BS handovers in a wireless network using Theorem 2 of [54]:

$$ \nu_c = \frac{4\sqrt{\lambda}}{\pi}v \qquad \nu_b = \frac{n\sqrt{\lambda}}{\pi}v $$

Let $T_c$ and $T_b$ be the time of a BS handover and a beam reselection, respectively. We define the total overhead per unit of time $T_o^M$ as the fraction of time in which MUs are not transmitting to the network due to cell handovers and beam reselections:

$$ T_o^M = \nu_c T_c + \nu_e T_b. \qquad (B.4) $$

An interesting thing to note in the definition of the overhead is the existence of a velocity $v$ such that $T_o = 1$. In this case, the moving MU spends all of its time reshuffling for new beams or new BS to connect to, and will not be able to transmit to the network. This velocity, denoted as $v_{max}$, is equal to:

$$ v_{max} = \frac{\pi}{\sqrt{\lambda}(4T_c + nT_b)}. \qquad (B.5) $$

With the previous definitions and results, we can define the *effective* Shannon rate (ESR) $\mathcal{R}_{eff}$ as:

$$ \mathcal{R}_{eff}(v) = \mathcal{R}(1 - T_o(v))^+, \qquad (B.6) $$

where $x^+ \triangleq \max(0, x)$. This definition of the ESR comes from an ergodic interpretation: when the MU is reselecting a new beam, or performing a handover at the boundary between two cells in the network, it cannot transmit. Figure B.3 displays the overhead of a moving MU in a single-tier network.

Mobility-induced beam misalignment

Due to intra-cell mobility, an MU has to reselect a beam when it moves from the coverage of one beam to another. Such a beam reselection in a 5G network occurs during a synchronization signal block (SSB) burst with period $\tau$. If the MU moves out of the main lobe of its original connection beam (also called the reference beam) between two consecutive SSB bursts without selecting a new beam, a beam misalignment occurs, namely, the MU receives from the serving BS via a side lobe. We assume, as in [54], that the probability $p_{bm}^M$ that there is a beam misalignment event for an MU moving at velocity $v$ is equal to:
<figure><img src="image_16.png" /><figcaption>Figure B.3: Illustration of the total overhead for a moving MU. Each red (resp. blue) line represents the amount of time needed for a beam reselection (resp. a cell handover).</figcaption></figure>

$$p_{\text{bm}}^{M}(v) = 1 - \exp(-v\nu_{M,b}\tau), \qquad (\text{B.7})$$

where $\nu_{M,b}$ is the time intensity of beam reselections. Hence, taking this beam misalignment into account, the antenna gain $g_{M,0}$ at the serving BS of the typical MU is given by:

$$g_{M,0} = \begin{cases} G_m^M & \text{w.p. } 1 - p_{\text{bm}}^M(v) \\ G_s^M & \text{w.p. } p_{\text{bm}}^M(v). \end{cases} \qquad (\text{B.8})$$

We also have to take beam misalignment into account when computing the overhead of users in the system: if the time intensity of beam reselections is higher than the frequency of SSBs, beam reselections happen during each SSB. If the time intensity of beam reselections is lower, the user may stay in the same reference beam between two consecutive SSB bursts, and thus, there are no beam reselection. Thus the effective time intensity $\nu_e$ of beam reselections in macro tier is:

$$\nu_e = \min \left( \frac{1}{\tau}, \nu_b \right).$$

## B.4 Bandwidth partitioning

The densification of wireless communication networks and the arrival of new generation communication protocols set new standards in terms of network performance, among which reliability and latency take a prominent place. To match these standards, the latest generation of telecommunication networks have introduced service differentiation as a mean to provide more flexibility and broader capacities for wireless networks. An important improvement
introduced in 5G networks is the adaptive use of multiple frequency bands available for communications ([1], [2], [88]). The idea is to increase the network capacity through a more flexible frequency band allocation, adapted to the needs of users. In such a setup, the number of bands allocated to a transmission depends on its nature: users transmitting less information (phone calls or text messages) will be allocated less frequency bands than users transmitting large amounts of data (video streaming for instance).

Service differentiation has been extensively studied in queueing networks. For instance, BCMP type queueing networks ([24]) extend the Jackson framework to the multi-class setting. When it comes to bandwidth sharing, an important problem comes from the allocation of resources to users in the network. A most common allocation is to attribute radio resources proportionally to the needs of users as to maintain fairness in the system (see [56], [67]).

In 4G LTE and 5G wireless networks, this bandwidth partitioning is implemented in a time divisive way: when transmitting, messages are encoded during radio frames, that are each divided in a certain number $N_{sc}$ subframes. Then, each subframe is divided in a resource grid, where frequencies are discretized, of resource blocks. These resource blocks consist in a certain number $N_{sc}^{RB}$ (12 for 5G NR networks for instance) of physical resource blocks (PRBs) that can support in parallel $N_{symb}^{slot} \cdot 2^{\mu}$ transmissions. When a MU transmits on the network, it receives a number of resource blocs, usually contiguous to ease the allocation, of PRBs, on which it will transmit. Figure B.4 gives an illustration of the

<figure><img src="image_17.png" /><figcaption>Figure B.4: Physical layer of 5G NR networks (source: [87])</figcaption></figure>

When considering bandwidth partitioning in our study, we will assume that the frequency space is divided in *K* channels, and that each MU receives a set of bands according to a probability distribution that represents user needs.
In a first time, we will assume that all PRB configuration are possible (al-
though this is not the case in real-world applications), and we may relax this
assumption later in the manuscript.
<figure><img src="image_18.png" /></figure>
# Part I

## Multiclass spatial queuing networks
<figure><img src="image_19.png" /></figure>
# MATHEMATICAL FRAMEWORK: MULTICLASS SBD PROCESSES

## 1.1 Network setup: bandwidth partitioning

Let us 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$. We assume that the region $\mathcal{D}$ is a square torus. This choice is motivated by the difficulty to obtain results on queuing systems over the whole Euclidean plane (see [81]).

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 framework defines a Poisson dipole network, as introduced in [12]. 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*. The process $\Phi_t$ is thus a collection of random segments of fixed length for which one of the ends is a part of a Poisson point process. In this part, we will start by studying the case $r=0$, in which the receivers and the transmitters processes are the same as the user process. In this case, we will denote as $\Phi_t$ this unique process. Later in Part I, we will generalize the result to the case $r > 0$.

$K$ orthogonal transmission channels of equal width are available to users arriving in the network. Channels $i \ne j$ do not interfere with one another. Let $\mathcal{P}(K)$ denote the set of non-empty subsets of $[1, 2, \dots, K]$ (in the dynamics
involved in this thesis, a transmitter using $\emptyset$ as a communication channel cannot transmit any signal to its receiver and will not interfere with other users in the network). An arriving user selects a set of channels $C \in \mathcal{P}(K)$ on which to transmit according 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 this chapter, 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. Finally, let $\mathcal{N}_0$ denote the thermal noise density in the network.

## Wireless interactions and service times

Let $\ell$ be a non-negative, bounded and non-increasing path-loss function with $\langle \ell_{\mathcal{D}} \rangle < \infty$. 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 (1.1)$$

where $C_z$ is the class of the transmitter located at $z$. Note that a receiver does not interfere with its own transmitter. We assume a low-SINR approximation, which allows us to linearize the Shannon-Hartley formula (see Theorem B.1 and Equation (B.2)). Here, $|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, the rate-of-transmission function for the receiver located at $x$ is:

$$R(x, \Phi_t) = \frac{|C_x|\ell(r)}{\mathcal{N}_0 + I(x, \Phi_t)}, \quad (1.2)$$

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 (1.3)$$

Using this setup, we can formulate this problem as follows:

**Definition 1.1.** The stochastic process $\Phi_t$ is a Markov jump process on the set of counting measures on $\mathcal{D}$, $M(\mathcal{D})$, which we call a BWP model. At all times $t > 0$, $\Phi_t$ is a point process on $\mathcal{D}$.

We define $\Delta : (\Phi_0, \lambda, \mathbf{L}, R) \mapsto \Phi_t$ as the realization of the dynamics, where $\Phi_0$ is the initial condition of the network, $\lambda$ is the arrival rate, $\mathbf{L}$ is the vector of average file sizes and $R$ is the rate-of-transmission function.
To understand why $\Phi_t$ is a Markov jump process, we can reformulate the dynamics presented earlier using a stochastic recurrence relation for each of the subprocesses $\Phi_{C,t}$. Let $(T_{C,n})_{n \le 0}$ be 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} + \tau_{C,n}, \end{array} \right. \qquad (1.4) $$

where $\mathcal{B}_{C,n}$ and $\mathcal{D}_{C,n}$ are $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 $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 rate $\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 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 d_{C,n}$ or if $\Phi_{C,n}(\mathcal{D})$ in the null measure, or equal to $\delta_x$, for some $x \in \Phi_{C,n}$.

Finally, the time between events is equal to $\tau_{C,n} = \min\{b_C, \min d_{C,n}\}$. With this stochastic recurrence, we can see that $\Phi_t$ is a jump process, and the Markov property comes from the exponential distribution of the file sizes for each user and of the time between arrivals.

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} \qquad (1.5) $$

whenever the integrals are defined. We use here the same notation for two different notions; this choice is motivated by results presented in this part. Table 1.1 summarizes the notation we use in this part.

## 1.2 Monotonicity and irreducibility

A first result about BWP models we can prove is the following.

**Theorem 1.1.** $\Phi_t$ is a $\phi$-irreducible Markov jump process on $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.* Let $\phi$ be the set function defined as:
<table><caption>Table 1.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, ... 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_{x \in D} \ell(\|x\|)dx$</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>

- $\phi(\{0_{\mathcal{M}(S)}\}) = 1$, where $0_{\mathcal{M}(S)}$ 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 $S$, we define the event $A_n = \{\Phi \in \mathcal{M}(S) : \Phi(S) = 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), \qquad (1.6) $$

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 [26], 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}(S)$.

To obtain $\phi$-irreducibility, we will proceed in two steps: in the first step, we prove that, from each state $X$ 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$, the SBD process is in state $x = \{x_1, \dots, x_n\} \in S^N$ with $N$ being the number of points in the SBD. The probability that we have the departure of the $N$ points 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 $S$, and let us define $A_n$ as previously. By definition, $\phi(A_n) > 0$. The probability to reach $A_n$ from the empty state $\emptyset$ can be expressed as:

$$ \mathbb{P}^n(\emptyset, A_n) = \int_{y_1 \in B_1} \dots \int_{y_n \in B_n} \int_{t_1=0}^{\infty} \int_{t_2=t_1}^{\infty} \dots \int_{t_n=t_{n-1}}^{\infty} H(dy_1) \dots H(dy_n) \\ e^{-\lambda p_{C_1} t_1} e^{-\lambda p_{C_2} (t_2-t_1)} \dots e^{-\lambda p_{C_n} (t_n-t_{n-1})} \lambda p_{C_1} dt_1 \dots \lambda p_{C_n} dt_n \prod_i \frac{b_{C_i}}{b_{C_i} + d_{C_i}}, $$

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$

The second property we need for our study is stochastic monotonicity, which we use to obtain bounds for the dynamics of the network through the definition of queuing systems that dominate or are dominated by the original dynamics. We need the following theorem:

**Theorem 1.2.** 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) $\mathbf{L} \le \mathbf{L}'$;

iii) for all point processes $\Psi \le_i \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 [81].

Let us take $0 < \lambda < \lambda'$ with the same fixed initial condition $\Phi_0$, a vector $\mathbf{L}$ and a rate-of-transmission function $R$, and let $\Phi_t = \Delta(\Phi_0, \lambda, d)$ and $\Phi'_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$ and assume that the next event happens at a time $\hat{t} \ge t$. It 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)$:

$$ 
\begin{aligned}
R(x, \Phi_{\hat{t}^-}) &= \frac{|C|\ell(r)}{\mathcal{N}_0 + \sum_{U \in \mathcal{P}(K)} \sum_{y \in \Phi_{U, \hat{t}^-}} |C \cap U| \ell(\|x - y\|)} \\
&\ge \frac{|C|\ell(r)}{\mathcal{N}_0 + \sum_{U \in \mathcal{P}(K)} |C \cap U| \left( \sum_{y \in \Phi_{U, \hat{t}^-}} \ell(\|x - y\|) + \sum_{y \in \Phi'_{U, \hat{t}^-} \setminus \Phi_{U, \hat{t}^-}} \ell(\|x - y\|) \right)} \\
&= R(x, \Phi'_{\hat{t}^-}).
\end{aligned}
$$

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 [31]): 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 [31], 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.

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

Theorem 1.2 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 of the network.
## 1.3 Critical arrival rate

In the definition of the dynamics we presented here, the value of the arrival parameter $\lambda$ plays an important role on the behavior of the dynamics of the system. Intuitively, if $\lambda$ becomes too large, then no receiver-transmitter pair can transmit and the population of the system explodes, preventing arriving users from transmitting. Conversely, if $\lambda$ becomes very low, transmissions can terminate and we can conjecture that the population of the system can stay bounded. The study of this particular kind of network is also dependent on inter-class interactions, as not all users interfere with each other.

The first result linking the value of $\lambda$ and the stability of the system is the following:

**Theorem 1.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 1.2: 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 [80] 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(\mathcal{D})$. We know that the Markov chain is either recurrent or transient (see Theorem A.10). Using stochastic monotonicity in the network from Theorem 1.2, 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$

Once the two values of $\lambda_c^+$ and $\lambda_c^-$ are defined, a natural question arising is their respective values. In a favorable case, these two values are equal, which would give a clear description of the stability region of the network. We can prove this intermediate result, giving a simple lower bound for $\lambda_c^-$ and an upper bound for $\lambda_c^+$:

**Lemma 1.4.** $\lambda_c^-$ and $\lambda_c^+$ satisfy:

$$ \frac{\ell(r)}{K\bar{L}\langle\ell_{\mathcal{D}}\rangle} \le \lambda_c^- \le \lambda_c^+ \le \frac{K\ell(r)}{\underline{L}\langle\ell_{\mathcal{D}}\rangle}, \quad (1.7) $$

where $\bar{L} = \max_C L_C$, $\underline{L} = \min_C L_C$ and $\langle\ell_{\mathcal{D}}\rangle = \int_{x \in \mathcal{D}} \ell(\|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 [80]. Using the main result of [80], we know that the cut-off arrival rate is equal to $\frac{\ell(r)}{K\bar{\mathbf{L}}\langle\ell_{\mathcal{D}}\rangle}$ for $\Phi_{d,t}$ and to $\frac{K\ell(r)}{\underline{\mathbf{L}}\langle\ell_{\mathcal{D}}\rangle}$ for $\Phi_{u,t}$.

We apply Theorem 1.2, 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_{\mathcal{D}}\rangle} > 0$ and $\frac{K\ell(r)}{\underline{\mathbf{L}}\langle\ell_{\mathcal{D}}\rangle} < \infty$, we can conclude that $0 < \lambda_c^{-}$ and $\lambda_c^{+} < \infty$. $\square$

When taking $K = 1$ in (1.7), we obtain:

$$\lambda_c^{-} = \lambda_c^{+} = \frac{\ell(r)}{\langle\ell_{\mathcal{D}}\rangle L},$$

which is the value of the critical arrival rate obtained in [80]. In the light of this remark, we can hope to obtain that $\lambda_c^{-} = \lambda_c^{+}$ in the multiclass system as well, and in this case, which we will also call critical arrival rate.
STABILITY OF SYMMETRIC
SYSTEMS

In this chapter, we will study the BWP model in a particular case, which
introduces several simplifications in the combinatorics of the system: *sym-
metric* systems.

**Definition 2.1.** A BWP model is symmetric, if all users transmitting on the same number of bands have the same stochastic properties, i.e., whenever $|C| = |D|$, we have $p_C = p_D$ and $L_C = L_D$.

An instance of a symmetric model 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 prob-
ability *p*<sub>j</sub> 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 ran-
dom among the binning possibilities. Thus, the probability that a user transmits
on a given set *C* of *j* bands is equal to *p*<sub>*C*</sub> = <p<sub>j</sub> > / <K>. The quantity of information
this arriving user has to transmit to the network is thus sampled from an ex-
ponential distributions with parameter *L*<sub>*j*</sub> with 1 ≤ *j* ≤ *K*. Setting *L*<sub>*C*</sub> = *L*<sub>*j*</sub>
and *p*<sub>*C*</sub> = <p<sub>j</sub> > / (<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 2.1, when proving Theorem 2.2.

In this chapter, we study the stability of symmetric BWP models as
defined in the previous parts. Using Theorem 1.2, we bound from above
and from below the dynamics in the network in order to obtain meaningful
bounds in the symmetric case to describe the stability region of the system
in terms of the value of λ.
## 2.1 Discretization of the dynamics

To study our dynamics, we introduce the following discretization: for $\varepsilon > 0$, we tessellate $\mathcal{D}$ in $N_\varepsilon$ square cells of side length $\varepsilon$, such that the origin is at the center of its cell. We denote by $A_i$ the cell centered at $a_i \in \mathcal{D}$. Finally, we introduce the stochastic process $\mathbf{X}_\varepsilon(t) = (X_{i,C}(t))_{0 \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 [81]) $\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:

$$ \left\{ \begin{array}{ll} \ell^{\varepsilon}(x, y) = \ell^{\varepsilon}(a_i, a_j) & \forall x \in A_i, y \in A_j, \text{ with} \\ \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{array} \right. \qquad (2.1) $$

and

$$ \left\{ \begin{array}{ll} \ell_{\varepsilon}(x, y) = \ell_{\varepsilon}(a_i, a_j) & \forall x \in A_i, y \in A_j, \text{ with} \\ \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{array} \right. \qquad (2.2) $$

where, for $0 \le i \le N_\varepsilon - 1$, we define $\mathcal{V}_i = \{b_i, \|b_i - a_i\| \in \{0, \varepsilon\}\}$. This way, 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). \qquad (2.3) $$
Using the dominated convergence theorem, we get:

$$ \lim_{\epsilon \to 0^+} \langle \ell_{\mathcal{D}}^{\epsilon} \rangle = \lim_{\epsilon \to 0^+} \langle \ell_{\epsilon, \mathcal{D}} \rangle = \langle \ell_{\mathcal{D}} \rangle, \qquad (2.4) $$

where $\langle \cdot, \mathcal{D} \rangle$ is defined in (1.5). Because of the square torus topology of $\mathcal{D}$, we have:

$$ \forall i, \quad \sum_{k=0}^{N_{\epsilon}-1} \ell^{\epsilon}(a_i, a_k) = \sum_{k=0}^{N_{\epsilon}-1} \ell^{\epsilon}(0, a_k) = \frac{1}{\epsilon^2} \langle \ell_{\mathcal{D}}^{\epsilon} \rangle. \qquad (2.5) $$

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

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^{\epsilon}$ and $\ell_{\epsilon}$ are respectively equal to:

$$ \begin{aligned} \bar{I}_{i,C}(t) &= \sum_{k=0}^{N_{\epsilon}-1} \sum_{U \in \mathcal{P}(K)} |C \cap U| \ell^{\epsilon}(a_k, a_i) (\bar{X}_{k,U}(t) - \mathbf{1}_{\{U=C, i=k\}}) \\ \underline{I}_{i,C}(t) &= \sum_{k=0}^{N_{\epsilon}-1} \sum_{U \in \mathcal{P}(K)} |C \cap U| \ell^{\epsilon}(a_k, a_i) (\underline{X}_{k,U}(t) - \mathbf{1}_{\{U=C, i=k\}}). \end{aligned} $$

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)}, \quad \underline{R}_{i,C}(t) = \frac{|C|}{\mathcal{N}_0 + \underline{I}_{i,C}(t)}. $$

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

*Proof.* From 2.3, we can obtain that for all $C \in \mathcal{P}(K)$, $0 \le i \le N_{\epsilon} - 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). \qquad (2.6) $$

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

Using Theorem 2.1, we can now reduce the study of the SBD process $\Phi_t$ to the study of the stability of the two processes $\underline{X}_{\epsilon}$ and $\bar{X}_{\epsilon}$. 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.
## 2.2 Fluid model and fluid limits

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 (2.7) $$

The goal of this section is to find a condition on $\lambda$ so that the Markov chain $\bar{X}_\epsilon(t)$ is positive recurrent.

## 2.3 Main result

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 2.2.** 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}}, \quad (2.8) $$

where $\mathcal{L} \triangleq \sum_C p_C |C| L_C$.

We start by proving Theorem 2.2 in the case $r = 0$ (we remind that we assumed, without loss of generality, that $\ell(0) = 1$). We then extend the result in the case $r > 0$ to obtain the intended stability condition.

### Stability of the dominating chain

We start by studying the stability of the chain $\bar{X}_\epsilon$. The stability condition for the chain is obtained though the study of its fluid limits. We can state the following result:

**Theorem 2.3.** The fluid limits for the chain $\bar{X}_\epsilon$ exist and are solutions to the following system of equations:
$$ \bar{x}'_{i,C}(t) = \lambda p_C \varepsilon^2 - \frac{1}{L_C \sum_{k=0}^{N_\varepsilon-1} \sum_{U \in \mathcal{P}(K)} |C \cap U|^{\ell_\varepsilon}(a_k, a_i) \bar{x}_{k,U}} \quad \text{if } \bar{x}(t) \neq 0 \\ \bar{x}_0 = \bar{X}_\varepsilon(0). $$

In this theorem, the convergence is taken in the sense of (A.1). The proof of this result, based on a construction presented in [76], is presented in Section 3.1. It is to note that these dynamics are not well-defined in the case $\bar{x}(t) = 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.

Once the evolution of the fluid limit is established, we can obtain the following Theorem:

**Theorem 2.4.** Let $\varepsilon > 0$ and $\bar{\lambda}_\varepsilon = \frac{K}{\mathcal{L}(\ell_\varepsilon)}$. If $\lambda < \bar{\lambda}_\varepsilon$, then the chain $\bar{X}_\varepsilon$ is stable.

*Proof.* To obtain stability for the chain $\bar{X}_\varepsilon$, we will rely fluid limits and Lemma A.13. We start by obtaining the existence of fluid limits for this chain. The proof of this result is presented in Section 3.2 of the next chapter. $\square$

We know from Theorem 2.1 that the process $\bar{X}_\varepsilon$ dominates the original dynamic. From Theorem 1.2, we obtain that, for all $\varepsilon > 0$, $\lambda_c \ge \frac{K}{\mathcal{L}(\ell_\varepsilon)}$. Taking the limit as $\varepsilon$ goes to 0 gives us:

$$ \lambda_c^- \ge \frac{K}{\mathcal{L}(\ell_\varepsilon)}. $$

### Instability of the dominated chain

We know from Theorem A.14 that the conditions to obtain instability for Markov chains are weaker than the ones to obtain stability. 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 use stochastic domination and an adequate Markov chain to bound $X_\varepsilon$ from below. We have the following theorem:

**Theorem 2.5.** Let $\underline{\lambda}_\varepsilon = \frac{K}{(\ell_{\varepsilon, \mathcal{D}})\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 Section 3.3. We know that the network is unstable if $\lambda > \frac{K}{(\ell_{\varepsilon, \mathcal{D}})\mathcal{L}}$ for each $\varepsilon > 0$. Taking the limit as $\varepsilon$ goes to 0, we conclude:

$$ \lambda_c^+ \le \frac{K}{(\ell_\varepsilon)\mathcal{L}}. \qquad (2.9) $$
When combining the results of Theorem 2.4 and Theorem 2.5, we obtain that:

$$ \lambda_c^- = \lambda_c^+ = \frac{K}{\langle \ell_D \rangle \mathcal{L}}, \qquad (2.10) $$

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

## Generalization to $r > 0$

In the previous section, all the calculations were made with the assumption that $r = 0$, i.e., in the case where the receiver and the transmitter are located in the same cell in the discretized system. To generalize the result of Theorem 2.2 to an arbitrary link length $r > 0$, we have to consider 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 us define $\bar{Z}_\varepsilon(t) = (Z_{i,C}(t))$ the $N_\varepsilon \times 2^K - 1$ vector 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$ vector 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.

We define the process:

$$ \bar{S}(t) = (\bar{X}_\varepsilon(t), \bar{Z}_\varepsilon(t), \bar{M}_\varepsilon(t)), $$

which is a Markov chain with countable state space. The dynamics of $\bar{S}(t)$ are as follows: to each receiver of class $C$ arriving in cell $i$, we pick uniformly a point $x \in 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 to decide in which cell the transmitter is located.

To obtain the interference in the network, we have to take into account the positions of the receivers. The interference experienced by a transmitter of class $C$ located in cell $i$ with its receiver located in cell $k_C$ in this new system becomes:

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

The rate-of-transmission function changes to:

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

Using the definitions of the process $\bar{S}(t)$ and the same steps as for Theorem 2.3, we can obtain the existence of fluid limits for the process $\bar{S}_\varepsilon(t)$, denoted as $\bar{s}(t)$:

$$ \begin{aligned} \frac{d}{dt} \bar{x}_{i,C}(t) &= \lambda p_C \varepsilon^2 - \frac{1}{L_C} \frac{|C| \ell(r) \bar{x}_{i,C}(t)}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| \bar{z}_{k,U}(t)} \\ \triangleq \lambda p_C \varepsilon^2 - \psi_{i,C}(\bar{s}(t)). \end{aligned} $$
We can now use similar arguments as for Theorem 2.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}. \qquad (2.11) $$

To obtain the reciprocal, we define similarly a chain $\underline{S}(t)$, with the same interference and rate-of-transmission - using the adequate quantities - as for $\bar{S}(t)$. From this, we can use a similar argument as the ones used to obtain Theorem 2.5 to bound the chain $\underline{S}(t)$ by an appropriate M/M/1 queue and obtain its instability whenever:

$$ \frac{K\ell(r)}{\mathcal{L}\langle\ell_D^\epsilon\rangle} \le \lambda_c^- \qquad (2.12) $$

Letting $\epsilon$ go to 0 gives us the intended result and proves Theorem 2.2.

### Simulations

To obtain a numerical validation of Theorem 2.2, we can simulate the dynamics of the symmetric network. To do so, we consider a simple network setup with 2 channels, meaning there are $2^2 - 1 = 3$ classes denoted as $\{1\}$, $\{2\}$ and $\{1, 2\}$. We take as 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 use the jump-hold construction of the Markov jump process: 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 $\mathbf{t}_D$ where each coordinate corresponds to a user located at $x_i$ of class $C_i$ and is drawn from an exponential distribution with mean $\frac{1}{L_{C_i}}R(x_i, \Phi_t)$.

We then compare the values of $t_b$ and the minimum of $\mathbf{t}_D$: if the former is lower, a user arrives at a location $x$ taken uniformly at random in $\mathcal{D}$; if the latter is lower, the corresponding user leaves the system. We then update the value of the interference in the system and move forward in time to the next event. If $(T_0 = 0, T_1, \dots, T_n)$ are the event times in the network and $X$ is the embedded Markov chain, we obtain that $\Phi_{T_n} = X_n$, which gives us the representation of process $\Phi_t$.

Figure 2.1 shows the behavior of the system for two values of $\lambda$: one that shows stability, and the other, instability.

As expected by the result of Theorem 2.2, when $\lambda$ is lower than the value $\lambda_c$, the population in the network stays bounded, as if $\lambda$ is larger than $\lambda_c$, the population diverges to infinity. A reliable criterion to study the stability of the system is Little's law (see [64]), 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:
<figure><img src="image_21.png" /><figcaption>Figure 2.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>

$$L = \lambda W. \qquad (2.13)$$

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

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

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 (2.14). If the time lived in the network grows linearly and diverges to infinity, then the network is not stable.

Figure 2.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_22.png" /><figcaption>Figure 2.2: Average staying time of users in the network with $K = 2$ channels for two different values of $\lambda$: $\lambda = 0.9\lambda_c$ (left) and $\lambda = 1.1\lambda_c$ (right).</figcaption></figure>

The implementation used for simulation runs the system in linear type with respect to the number of steps desired, so we can further verify this conjecture for other configuration, though, the number of possible configuration
grows exponentially with $K$. To be able to consider a system with a large number of bands, we can consider not allowing all possible configurations in the network, which contradicts our initial assumption that $p_C > 0$ for all $C \in \mathcal{P}(K)$. In Chapter 5, we will discuss how blocking certain configurations of bands may change the stability condition for the system.

Figure 2.3 shows the sojourn time in a system with $K = 4$ bands (the legend has been omitted for readability). In the stable case (on the left), we can observe that users with the same number of bands have sojourn times converging to the same values. This property of the system will be investigated in Part II of this thesis.

<figure><img src="image_23.png" /><figcaption>Figure 2.3: Staying time in the system for $K = 4$ bands, with $\lambda = 0.95\lambda_c$ on the left and $\lambda = 1.05\lambda_c$ on the right.</figcaption></figure>

In this chapter, we have seen how fluid limits and combinatorics allows us to obtain a simple closed-form for the critical arrival rate for the system with the symmetry assumption. As mentioned in the introduction, this symmetry assumption is featured in a variety of resource allocation policies. In the next chapter, we will study the system without the symmetry assumption.
<figure><img src="image_24.png" /></figure>
# PROOFS OF THEOREMS OF CHAPTER 2

This chapter focuses around proving the Theorems about stability of the dominating chain $\bar{X}_{\varepsilon}$ and the instability of the chain $X_{\varepsilon}$ presented in Chapter 2.

We start by establishing the fluid model for the dynamics by proving Theorem 2.3. Then, we move on to proving the stability of the chain $\bar{X}_{\varepsilon}$, establishing Theorem 2.4. The last part of this chapter proposes a proof for the instability of $X_{\varepsilon}$, giving Theorem 2.5.

## 3.1 Proof of Theorem 2.3

Let us take a sequence of initial conditions $\bar{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{X}_{\varepsilon}(nt))$ and to obtain a system of ODEs for which the limit of this sequence is a solution.

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

$$ \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{1}{L_C}\int_0^{nt}\bar{X}_{i,C}(u)R_{i,C}(\bar{X}_{\varepsilon}(u))du\right). \quad (3.1) $$

We use the variable change $u = ns$ in the integral to get:
$$
\begin{align*}
& \int_0^{nt} \bar{X}_{i,C}(u) R_{i,C}(\bar{\mathbf{X}}_{\varepsilon}(u)) du \\
&= \int_0^t \bar{X}_{i,C}(ns) R_{i,C}(\bar{\mathbf{X}}_{\varepsilon}(ns)) nds \\
&= \int_0^t \frac{|C| \bar{X}_{i,C}(ns)}{\mathcal{N}_0 + \sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| (\bar{X}_{k,U}(ns) - \mathbf{1}_{\{U=C, i=k\}})} nds \\
&= n \int_0^t \frac{|C| \frac{1}{n} \bar{X}_{i,C}(ns)}{\frac{\mathcal{N}_0}{n} + \sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| (\frac{1}{n} \bar{X}_{k,U}(ns) - \frac{1}{n} \mathbf{1}_{\{U=C, i=k\}})} ds.
\end{align*}
$$

We define:

$$
R_{i,C}^{n}(x) = \frac{|C|}{\frac{\mathcal{N}_0}{n} + \sum_{k,U} \ell^{\epsilon}(a_k, a_i) |C \cap U| \left(x_{k,U} - \frac{1}{n} \mathbf{1}_{\{U=C, i=k\}}\right)}.
$$

Using this, (3.1) 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 \epsilon^2 nt) \\
\qquad - \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}}_{\epsilon}(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 \epsilon^2 nt) \\
\qquad - 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}}_{\epsilon}(ns) \right) ds \right). \quad (3.2)
$$

We can note that $M_{i,C}^n$ is a martingale and $\bar{Y}_{i,C}^n$ is its compensator, for
which we obtain a limit in order to obtain the fluid limit. By using (3.1), 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}}_{\epsilon}(ns) \right) ds. \quad (3.3)
$$

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 3.1.** If $\bar{\mathbf{x}}_0$ *exists and is finite*, then:

$$
\forall t \ge 0, \forall i, C, \lim_{n \to \infty} \bar{Y}_{i,C}^n(t) = \bar{X}_{i,C}(0) + \lambda p_C \epsilon^2 t. \quad (3.4)
$$
*Proof.* Using the strong law of large numbers yields, $\mathbb{P}$-almost surely:

$$ \lim_{n \to \infty} \frac{1}{n} 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{X}_\varepsilon(ns) \right) \mathrm{d}s \\ &= \frac{1}{L_C} \int_0^t \frac{1}{n} \bar{X}_{i,C}(ns) \frac{|C|}{\frac{N_0}{n} + \sum_{k,U} |C \cap U| \ell^\varepsilon(a_k, a_i) \frac{1}{n} \bar{X}_{k,U}(ns)} \mathrm{d}s \\ &\le \frac{1}{L_C} \int_0^t \frac{1}{n} \bar{X}_{i,C}(ns) \frac{|C|}{|C| \ell^\varepsilon(a_i, a_i) \frac{1}{n} \bar{X}_{i,C}(ns)} \mathrm{d}s \\ &= \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 [26]) to obtain that the sequence of processes $\{\frac{1}{n}\bar{X}_{i,C}(n\cdot) R_{i,C}^n(\frac{1}{n}\bar{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{X}_\varepsilon(ns) \right) \mathrm{d}s \right) = 0, $$

which leads to the intended result. $\square$

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

For $v: [0, T) \to \mathbb{R}$ and $0 \le t_1 \le t_2 \le T$, let us 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 [51]):

**Definition 3.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} \left[ \sup_{0 \le t \le T} |V^n(t)| \ge K_\eta^0 \right] < \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, for all $\eta > 0$, we have:

$$ \mathbb{P}[\omega(V, \delta_\eta^0, t_0) > \eta] < \eta. \qquad (3.5) $$
This means, that, $\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 3.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 3.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. \quad (3.6) $$

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$. Using (3.3), 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. \quad (3.7) $$

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$. (3.3) gives 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 we remind that:

$$ \frac{1}{n}\bar{X}_{i,C}(nt)R_{i,C}^n(\frac{1}{n}\bar{\mathbf{X}}_\varepsilon(nt)) \le 1. \quad (3.8) $$

We have:

$$ \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 3.1, 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. 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). \quad (3.9) $$

Let us set $\delta_{\eta}^{0} = \eta (\lambda p_C \varepsilon^2 + \kappa_{\eta}^0 + \frac{1}{L_C})^{-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. \quad (3.10) $$

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{X}_{\varepsilon}(t))$.

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

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

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

$$ \left\{ \begin{array}{ll} \frac{d}{dt}\bar{x}_{i,C}(t) & = \lambda p_C \varepsilon^2 - \frac{1}{L_C} \frac{|C|\bar{x}_{i,C}(t)}{\sum_{k=0}^{N_{\varepsilon}-1} \sum_{U \in \mathcal{P}(K)} |C \cap U|^{\ell^{\varepsilon}(a_k, a_i)} \bar{x}_{k,U}} & \text{if } \bar{x} \neq 0 \\ \bar{x}(0) & = \bar{x}_0. \end{array} \right. \quad (3.11) $$

*Proof.* From Lemma 3.2, the sequences $\{\frac{1}{n}\bar{X}_{\varepsilon}(n\cdot), n \ge 1\}$ and $\{\bar{Y}^n, n \ge 1\}$ are both tight. It follows from Theorem 11.6.8 from [89] that the sequence $\{(\frac{1}{n}\bar{X}_{\varepsilon}(n\cdot), \bar{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{x}_{\varepsilon}, \bar{Y})$ as $l$ goes to infinity.

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

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

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

- $(\hat{x}_{\varepsilon}, \hat{Y}) \stackrel{d}{\sim} (\bar{X}_{\varepsilon}, \bar{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{X}^{n_l}(s)) ds. \quad (3.12) $$

Using the C-tightness of $\bar{X}_{\varepsilon}$, $\hat{X}$ and $\hat{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). \quad (3.13) $$

Combining all these results, we can take the limits as $l$ goes to infinity in (3.12) 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)} ds. \quad (3.14) $$

Let us remind that, from Lemma 3.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)} ds. \quad (3.15) $$

To obtain the intended equations, we differentiate (3.15), which concludes the proof. $\square$

## 3.2 Proof of Theorem 2.4

In this section, we provide a proof for Theorem 2.4 using Lemma A.13, which is applicable to our model. We have:

**Lemma 3.4.** Assume that there exists $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{X}_{\varepsilon}$ is positive Harris recurrent.

*Proof.* To obtain this result, we use Lemma A.13: the fluid limit $\mathbf{x}(t)$ meets the intended requirements. Let us take such a $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 [37] with the fluid limit $\bar{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| x_{i,C}}{\sum_{k,U} \ell^{\varepsilon}(a_k, a_i) |C \cap U| x_{k,U}}$, which is 0-homogeneous, and non-increasing in $x_{k,U}$ for $(k, U) \ne (i, C)$. $\square$
The goal here is to find an appropriate vector $z$ to get the desired upper bound on $\lambda$. Let us set $z_{i,C} = p_C L_C$ for all $C \in \mathcal{P}(K)$. We want that, for each $C \in \mathcal{P}(K)$:

$$ \lambda p_C \varepsilon^2 \leq \frac{1}{L_C} \frac{|C| p_C L_C}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| p_U L_U}. \quad (3.16) $$

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

$$ \lambda \leq \frac{1}{\langle \ell_D^\varepsilon \rangle} \frac{|C|}{\sum_U |C \cap U| p_U L_U}. \quad (3.17) $$

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

**Lemma 3.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:

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

*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 (3.19) $$

The number of sets $U \in \mathcal{P}(K)$ such that $|U| = l$ and $|C \cap U| = m$ is equal to $\binom{j}{m} \binom{K-j}{l-m}$. Using the symmetry property of $y$, we know that $y_U = y_{[1,l]}$. Using the fact that we have $y_U = \frac{1}{\binom{K}{l}} \sum_{|V|=l} y_V(t)$. We can rewrite (3.19) 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}}. \quad (3.20) $$

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, \quad (3.21) $$

which concludes the proof of Lemma 3.5.

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

$$ \lambda \le \frac{1}{\langle \ell_{\mathcal{D}}^{\varepsilon} \rangle} \frac{K}{\mathcal{L}}. \qquad (3.22) $$

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

## 3.3 Proof of Theorem 2.5

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

$$ r_{i,C} : \mathbf{x} \in \mathbb{R}_{+,*}^{N_{\varepsilon} \times 2^{K}-1} \mapsto \frac{1}{p_C L_C \varepsilon^2 \sum_{k,U} \ell_{\varepsilon}(a_k, a_i) |C \cap U|} \mathbf{x}_{i,C} $$

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:

$$ \forall n, \quad r_{0,\{1\}}(\mathbf{x}_n) = \frac{1}{p_{\{1\}} L_{\{1\}} \varepsilon^2}, \quad 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 is not defined at $\mathbf{x} = 0$, and it is 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

$$ \mathcal{B} = \arg\max_{\mathbf{x} \in \mathbb{R}_{+,*}^{N_{\varepsilon} \times 2^{K}-1}} \left\{ \min_{i,C} \frac{1}{p_C L_C \varepsilon^2 \sum_{k,U} \ell_{\varepsilon}(a_k, a_i) |C \cap U|} \mathbf{x}_{i,C} \right\}, \quad (3.23) $$

which is non-empty.

**Lemma 3.6.** Let $\mathbf{z} \in \mathcal{B}$ 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^*} \sum_{k,U} |C \cap U|} \frac{|C^*| z_{i^*, C^*}}{\ell_{\varepsilon}(a_k, a_i) z_{k,U}}. $$

Then, $\underline{X}_{\varepsilon}$ is transient.
Proof. Let $z \in \mathcal{B}$. By definition, there exists $i^*, C^*$ such that the maximum value of $\min_{i,C} r_{i,C}$ is equal to $r_{i^*,C^*}(z)$. Let us take $\lambda > 0$ such that:

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

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

- Arrivals in queue $i, C$ happen with rate $\lambda p_C \varepsilon^2$, for all $i, C$;

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

$$ - \frac{1}{L_{C^*}} \frac{|C^*| z_{i^*, C^*}}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i) z_{k,U}} \quad \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) + N_0} \quad \text{else.} $$

$\mathbf{Y}$ is a Markov jump process, with state space $\mathbb{N}^{N_{\varepsilon} \times 2^K - 1}$. By definition of $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{1}{L_{C^*}} \frac{|C^*| \underline{X}_{i^*, C^*}(t)}{\sum_{k,U} |C \cap U| \ell_{\varepsilon}(a_k, a_i) \underline{X}_{k,U}(t) + N_0}. $$

Using a coupling argument similar to the one used in the proof of Theorem 1.2, 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$

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

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 $z$:

**Lemma 3.7.** Let $z \in \mathcal{B}$. Then, $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}}. \qquad (3.24) $$

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

$$ \mathbf{z} \in \arg \max_{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 gener-
ality, 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_{\varepsilon}-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_28.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 (3.24). Using Lemma 3.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{|C^*| p_{C^*}}{\frac{1}{\varepsilon^2} \langle \ell_{\varepsilon, \mathcal{D}} \rangle \frac{|C^*|}{K} \sum_U |U| p_U L_U} \\ &= \frac{\varepsilon^2 p_{C^*} K}{\langle \ell_{\varepsilon, \mathcal{D}} \rangle \frac{|C^*|}{K} \sum_U |U| p_U L_U} \\ &= \frac{K}{\langle \ell_{\varepsilon, \mathcal{D}} \rangle \mathcal{L}}, \end{align*} $$

where (a) uses Lemma 3.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$ gives us the intended condition for $\lambda$ and concludes the proof.
# GENERALIZATION TO NON-SYMMETRIC SBD PROCESSES AND OTHER DYNAMICS

In this chapter, we will extend the result presented in Theorem 2.2 when relaxing the symmetry condition presented in Definition 2.1. This work was firstly motivated by the study of BWP models where the symmetry assumption was relaxed. In these systems, we allow a differentiated use of communication channels, with different arrival rates and different file sizes. Relaxing the symmetry hypothesis for BWP models allows to capture a more diverse array of association policies. When investigating this problem, we obtained a more general stability condition for an array of dynamics of queuing networks with state-dependent departure rates.

This chapter thus focuses on obtaining a general stability condition for queuing networks where the departure rate in a given server is a decreasing function of other queue sizes. To do so, we will focus on a queuing network consisting in *N* servers for which the queuing discipline is Markovian. We will consider a constant arrival rate with a given routing policy for jobs to be addressed to each server in the network, and the departure rate for jobs in the system will be depending on the state of the system at a given time. We will then apply the main result of this chapter to obtain a stability condition for the non-symmetric SBD process.
## 4.1 Queuing systems with state dependent departure rates

Let us consider a queueing network of $N \ge 2$ interacting queues, with $(X_i(t))$ denoting the queue length process. In each queue, users carrying jobs arrive at a constant rate and leave the system at a rate depending on the state of the queuing system at a given time. Our model is motivated by wireless networks, so we assume that the departure rate in queue $i$ is a decreasing function of the lengths of queues $X_j$, for $j \ne i$. Finally, we assume that users leave the system once their job is completed, and we assume that queue lengths are not bounded. Finally, we assume that the system is Markovian (such an example can be found in [84] for instance).

Arrivals in the system happen at rate $\lambda > 0$, and upon arrival, a user is routed independently with probability $p_i$ to server $i$. Thus, arrivals in server $i$ happen with rate $\lambda_i = \lambda p_i$. We assume that for all $i$, $p_i \ne 0$.

Departures in server $i$ happen with a state-dependent rate $\Psi_i(\mathbf{X}(t))$, where, for all $0 \le N-1$, $\Psi_i$ is decreasing in $X_j$ for all $j \ne i$. Furthermore, we assume that, for all $i$, $\Psi_i$ is such that the following limit exists for all $\mathbf{x} \in \mathbb{R}_+^N$:

$$ \lim_{n \to \infty} \Psi_i(n\mathbf{x}) = \psi_i(\mathbf{x}), $$

where $\psi_i$ is non-constant, continuous and 0-homogeneous, i.e., for all $s > 0$ and $\mathbf{x} \in \mathbb{R}_+^N$, $\psi_i(s\mathbf{x}) = \psi_i(\mathbf{x})$.

Using similar coupling arguments as the ones presented in Theorem 1.2 in Section 1.2 of Chapter 1, we can prove the following Theorem:

**Theorem 4.1.** Let us take $\mathbf{X}$ and $\mathbf{Y}$ to realizations of the dynamics with respective arrival rates $\lambda_X$ and $\lambda_Y$, initial conditions $\mathbf{X}(0)$ and $\mathbf{Y}(0)$, and departure rates $(\Psi_i)_{0 \le i \le N-1}$ and $(\Phi_i)_{0 \le i \le N-1}$. The following conditions are sufficient for $\mathbf{X}$ to dominate $\mathbf{Y}$ (when the two other parameters are taken equal):

A) $\mathbf{Y}(0) \le \mathbf{X}(0)$ coordinate-wise,

B) $\lambda_Y \le \lambda_X$

C) for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ with $\mathbf{x} \le \mathbf{y}$ coordinate-wise, $\Psi_i(\mathbf{x}) \le \Phi_i(\mathbf{y})$.

## Fluid model

Using a similar proof as the one from Theorem 2.3, we can obtain that the sequence of fluid-scaled processes $(\frac{1}{n}X_i(nt))_{n \ge 0}$ admits a limit $(\bar{x}_i(t))$ as $n$ goes to infinity, solution to the following equation:

$$ \bar{x}_i(t) = \bar{x}_i(0) + \lambda_i t - \int_0^t \psi_i(\bar{x}(s))ds. \quad (4.1) $$
Assumptions made about the functions $(\Psi_i)$ and $(\psi_i)$ give three properties of the fluid system:

i) for all $i$, $\psi_i$ is non-increasing in $x_j$ for all $j \neq i$ and increasing in $x_i$,

ii) for all $i$, $\psi_i$ is 0-homogenous, i.e., for all $t \in \mathbb{R}_+^N$ and $s \ge 0$, $\psi_i(st) = \psi_i(t)$,

iii) the $(\psi_i)$ are such that if $\sum_i \bar{x}_i(\tau) > 0$ for a given $\tau$, there exists $\delta > 0$ such that $\bar{x}_i(t) > 0$ for all $i$ and $t \in [\tau, \tau + \delta]$.

Let us define:

$$ \mathcal{B} = \left\{ \lambda \in \mathbb{R}_+ : \lambda \le \frac{1}{p_i} \psi_i(\mathbf{p}) \text{ for some } \mathbf{p} \in \mathbb{R}_+^N \text{ and all } i \right\}. $$

The main result about the stability of such dynamics is given by Lemma A.13, which allows to conclude about stability for the Markov chain $\mathbf{X}$. The main goal of this chapter is to obtain a converse for this result, and to find a characterization of the maximal value of $\mathcal{B}$ in a particular case.

## 4.2 An instability condition

The main Theorem of this section is a reciprocal to Lemma 3.4, which will give a condition for stability for the non-symmetric BWP model:

**Theorem 4.2.** If $\lambda > \lambda^*$, where $\lambda^* = \sup \mathcal{B}$, then the chain $\mathbf{X}$ is transient.

*Proof.* The first thing to note is that, because the functions $\psi_i$ are 0-homogeneous, finding a maximum for $\frac{1}{p_i}\psi_i$ over $\mathbb{R}_+^N$ can be resumed to looking its maximum over $(0,1]^N$. Because the $\psi_i$ are continuous, $\mathbf{x} \mapsto \min_i \frac{1}{p_i}\psi_i(\mathbf{x})$ is also continuous, and admits a maximum in $(0,1]^N$. Finally, the value for this function at $\mathbf{x} = 0$ is not defined: because $\psi_i$ is not constant, let us take $\mathbf{x}_1, \mathbf{x}_2$ such that $\psi_i(\mathbf{x}_1) \neq \psi_i(\mathbf{x}_2)$.

We define two sequences $\mathbf{x}_n = \frac{1}{n}\mathbf{x}_1$ and $\mathbf{y}_n = \frac{1}{n}\mathbf{x}_2$. We have $\lim_{n \to \infty} \mathbf{x}_n = \lim_{n \to \infty} \mathbf{y}_n = 0$, but $\lim_{n \to \infty} \psi_i(\mathbf{x}_n) = \psi_i(\mathbf{x}_1)$ and $\lim_{n \to \infty} \psi_i(\mathbf{y}_n) = \psi_i(\mathbf{x}_2)$, which are not equal.

Thus, we can take $\mathbf{z} \in \mathbb{R}_+^N$ such that:

$$ \mathbf{z} \in \arg \max_{\mathbf{x} \in \mathbb{R}_+^N} \min_i \frac{1}{p_i} \psi_i(\mathbf{z}). $$

Let $\lambda > \lambda^*$. By definition, there exists $i^*$ such that:

$$ \lambda p_{i^*} > \psi_{i^*}(\mathbf{z}). $$

Let us define the chain $\mathbf{Y}$ as follows:

$$ - \text{ For all } 0 \le i \le N-1, Y_i(0) = X_i(0) $$
- The arrival rate in $Y_i$ is equal to $\lambda p_i$

- The departure rate in $Y_i$ is equal to:

- $\psi_i(\mathbf{Y}(t))$ if $i \neq i^*$

- $\psi_{i^*}(\mathbf{z})$ if $i = i^*$

By definition of $\mathbf{z}$, we know that for all $t \ge 0$, $\psi_{i^*}(\mathbf{z}) \ge \Psi_{i^*}(\mathbf{X}(t))$. Hence, we can apply Lemma 4.1 to obtain, for all $0 \le i \le N-1$:

$$Y_i(t) \le X_i(t) \quad \mathbb{P}\text{-a.s.}$$

To conclude, we can remark that server $i^*$ is a M/M/1 queue with an arrival rate $\lambda p_{i^*}$ and a departure rate $\psi_{i^*}(\mathbf{z})$. We thus have, $\mathbb{P}$-almost surely $\lim_{t \to \infty} Y_{i^*}(t) = \infty$ if $\lambda \ge \lambda^*$, which allows us to conclude that $\mathbf{X}$ is transient.

□

We can remark that the $i^*$-th queue length is not the only one to diverge to almost surely diverge to infinity if $\lambda > \lambda^*$. We know that $X_{i^*}(t)$ goes to infinity. Because all the functions $\Psi_j$ are decreasing functions of $X_{i^*}$, for all $j$, there exists a time $T_j$ after which $\Psi_j(\mathbf{X}(T_j)) < \lambda p_j$. We can then bound from below the $j$-th queue by a M/M/1 queue whose arrival rate is strictly larger than its departure rate. Hence, we get:

$$\forall j, \quad \mathbb{P}\left[\lim_{t \to \infty} X_j(t) = \infty\right] = 1.$$

Combining the results of Lemma 3.4 and Theorem 4.2 give us the following result:

**Theorem 4.3.** Let $\lambda^* = \sup \mathcal{B}$. If $\lambda < \lambda^*$, the Markov chain $\mathbf{X}$ is positive recurrent. If $\lambda > \lambda^*$, it is transient.

The behavior of the Markov chain for $\lambda = \lambda^*$ is not studied in this thesis. We know that in the general case, the study of dynamic systems in the case $\lambda = \lambda^*$, although results can be obtained using fluid limits in some cases (see [35]).

An interesting remark for this Theorem comes from the consequences of Theorem 4.1. If the system meets a separability condition, one can use the saturation rule (see [20]). This rules states that the Markov chain representing our dynamics is either positive recurrent, or transient, and that the stability region is uniquely defined by the value of $\lim_{n \to \infty} \frac{T_n}{n}$, where $(T_n)_{n \ge 0}$ is the sequence of maximal daters in the system. Unfortunately, for our dynamics, defining maximal daters and computing the value of the limit seems to be a hard problem, which leads us to obtain stability through different channels.

The second main result of this Chapter is a better characterization of the value of $\lambda^*$, under a stronger assumption. Let us replace assumption $i$) for the functions $(\psi_i)_{0 \le i \le N-1}$ by assumption $i^*$):

$i^*$) for all $i$, $\psi_i$ is non-increasing in $x_j$ for all $j \neq i$ and increasing in $x_i$.
This stronger assumption is motivated by communication networks, where the departure rates of one user is inversely proportional to a shot-noise of the queue length process and proportional to the queue length of the observed server. To obtain the departure rate of all users in a given server, we have to multiply this departure rate by the queue length of said server, which gives a departure rate in a server increasing with its queue length (see [81], [84] for instance).

**Theorem 4.4.** Under assumptions i\*), ii), iii), $\lambda^*$, the maximal value of $\mathcal{B}$ is obtained for any $\mathbf{z}$ solution to, for all $0 \le i, j \le N-1$:

$$ \frac{1}{p_i} \psi_i(\mathbf{z}) = \frac{1}{p_j} \psi_j(\mathbf{z}). \qquad (4.2) $$

We thus have $\lambda^* = \frac{1}{p_i} \psi_i(\mathbf{z})$, for any $0 \le i \le N-1$.

The proof of this result uses the same arguments as the proof of Lemma 3.7, presented in Section 3.3 for the Appendix. This result gives us a simple condition to characterize the values of $\mathcal{B}$. In the next section, we will use this result to prove stability of non-symmetric BWP models.

## 4.3 Stability of non-symmetric BWP models

We can now move on to obtaining stability for the non-symmetric multiclass SBD process. In this setup, we allow users transmitting on the same number of bands to have different stochastic properties in the system. This way, we can encompass more complex allocation policies as the symmetric system does. We formulate the following result:

**Theorem 4.5.** The non-symmetric multiclass SBD system is stable if $\lambda < \lambda^*$ and unstable if $\lambda > \lambda^*$, where:

$$ \lambda^* = \frac{1}{p_C L_C} \frac{|C| z_C^*}{\langle l_D \rangle \sum_U |C \cap U| z_U^*}, $$

where $z^*$ is a solution to the following system of equations:

$$ \forall C, D \in \mathcal{P}(K) \quad \frac{1}{p_C L_C} \frac{|C| z_C^*}{\sum_U |C \cap U| z_U^*} = \frac{1}{p_D L_D} \frac{|D| z_D^*}{\sum_U |D \cap U| z_U^*}. \qquad (4.3) $$

*Proof.* The goal here is to obtain the stability condition for the two chains $\underline{\mathbf{X}}_\varepsilon$ and $\bar{\mathbf{X}}^\varepsilon$, and then to obtain the wanted stability condition for the original system.

Let us start with obtaining the stability condition for the dominating chain $\bar{\mathbf{X}}^\varepsilon$. We take, for all $i, C$, the function $\Psi_{i,C}(\mathbf{x}) = \frac{1}{L_C} \sum_{k,U} \frac{\mathcal{N}_{0+} |C| x_{i,C}}{\ell^\varepsilon(a_k, a_i) |C \cap U| x_{k,U}}$.

Immediately, we can see that, for all $\mathbf{x} \in \mathbb{R}_+^{N \times 2^K - 1}$, we have:
$$ \lim_{n \to \infty} \Psi_{i,C}(nx) = \frac{1}{L_C} \frac{|C|x_{i,C}}{\sum_{k,U} a_{k,i} |C \cap U| x_{k,U}} \triangleq \psi_{i,C}(\mathbf{x}). $$

Moreover, the functions $(\psi_{i,C})$ meet conditions $i^*) - iii)$, which allow us to use the result from Theorem 4.4 to get that $\mathbf{X}$ is stable whenever $\lambda < \lambda^*$, where $\lambda^*$ is such that:

$$ \lambda^* p_C \varepsilon^2 = \psi_{i,C}(\mathbf{z}), $$

where $\mathbf{z}$ is solution to (4.2).

Let us now prove that for all $i, j$ and for all $C$, $z_{i,C} = z_{j,C}$. Let us take $0 \le j \le N$, and define $\mathbf{z}' = (z_{i+j,C})$, where we denote $i + j$ as remainder modulo $N$ of the sum if $i + j \ge N$. We have, for all $i, C$:

$$ \begin{aligned} \psi_{i,C}(\mathbf{z}') &= \frac{1}{L_C} \frac{|C|z_{i+j,C}}{\sum_{k,U} \ell^\varepsilon(a_i, a_k) |C \cap U| z_{k,U}} \\ &\stackrel{(a)}{=} \frac{1}{L_C} \frac{|C|z_{i+j,C}}{\sum_{k',U} \ell^\varepsilon(a_{i+j}, a_{k'}) |C \cap U| z_{k',U}} = \psi_{i+j}(\mathbf{z}). \end{aligned} $$

To explain (a), we use the fact that $\mathcal{D}$ is a square torus, and that the function $\ell^\varepsilon$ is defined on a grid. This way, we know that, for each $i, j, k$, there exists a unique $k'$ such that $\ell^\varepsilon(a_{i+j}, a_k) = \ell^\varepsilon(a_i, a_{k'})$. We use this reindeixation of the sum to conclude. Finally, we know that $\mathbf{z}$ is solution to (4.2), which gives $\psi_{i+j}(\mathbf{z}) = \psi_i(\mathbf{z})$ and gives:

$$ \psi_{i,C}(\mathbf{z}') = \psi_{i,C}(\mathbf{z}). $$

From this result, we know that for any solution $\mathbf{z}$ to (4.3), the value of $\psi_{i,C}(\mathbf{z})$ is invariant by translation on $i$. From this, we can conclude that for all $i, j, z_{i,C} = z_{j,C} = z_C$. Using this, searching for solutions of the system over $\mathbb{R}^{N_\varepsilon \times 2^K-1}$ amounts to searching for solutions over $\mathbb{R}^{2^K-1}$. The system becomes:

$$ \begin{aligned} \frac{1}{p_C \varepsilon^2} \psi_{i,C}(\mathbf{z}) &= \frac{1}{p_C \varepsilon^2} \frac{1}{L_C} \frac{|C| z_C}{\sum_{k,U} \ell^\varepsilon(a_k, a_i) |C \cap U| z_U} \\ &= \frac{1}{p_C L_C} \frac{|C| z_C}{\langle \ell^\varepsilon_{\mathcal{D}} \rangle \sum_U |C \cap U| z_U}. \end{aligned} $$

After simplification, we get:

$$ \frac{1}{p_C L_C} \frac{|C| z_C}{\sum_U |C \cap U| z_U} = \frac{1}{p_D L_D} \frac{|C| z_C}{\sum_U |C \cap U| z_U}. \quad (4.4) $$

We can use the result from Theorem 4.4 and the domination properties from Theorem 1.2 to get the following stability condition for the system:

$$ \lambda_c^+ \le \frac{1}{p_C L_C} \frac{|C| z_C}{\langle \ell^\varepsilon_{\mathcal{D}} \rangle \sum_U |C \cap U| z_U}. \quad (4.5) $$
We can use the same steps with the dominated chain $\underline{X}_\varepsilon$ to obtain the same system as in (4.4), and the reciprocal condition:

$$ \lambda_c^- \geq \frac{1}{p_C L_C} \frac{|C| z_C}{\langle \ell_{\varepsilon, \mathcal{D}} \rangle \sum_U |C \cap U| z_U}. \quad (4.6) $$

Finally, let $\varepsilon$ tend to 0, which gives the condition:

$$ \lambda_c = \frac{1}{p_C L_C} \frac{|C| z_C}{\langle \ell_{\mathcal{D}} \rangle \sum_U |C \cap U| z_U}. \quad (4.7) $$

where $z$ is solution to (4.4), which concludes the proof. $\square$

We can make two remarks about the result of Theorem 4.4. The first one is that the system of equations we obtain for $z$ does not depend on the path-loss function $\ell$ in the system. This lack of dependency comes from the square torus topology we use in our system, which allows for a simplification of the system presented in Theorem 4.3.

The second remark we can make is that this result can be easily generalized to any link length $r > 0$ by taking the same steps as in Subsection 2.3, to get that:

$$ \lambda_c = \frac{\ell(r)}{p_C L_C} \frac{|C| z_C}{\langle \ell_{\mathcal{D}} \rangle \sum_U |C \cap U| z_U}, $$

To verify this result, we add the symmetry hypothesis to compare it to the result known in [72]: whenever $|C| = |D|$, we have $p_C = p_D$ and $L_C = L_D$, i.e. users with the same number of communication channels have the same stochastic properties in the network. With $K = 2$ bands, we look for a solution to (4.3) in the form $z = (x, y, z)$. The first equation gives:

$$ \frac{1}{p_{\{1\}} L_{\{1\}}} \frac{x}{x+z} = \frac{1}{p_{\{1\}} L_{\{1\}}} \frac{y}{y+z}, $$

which yields $x = y$. Injecting this result in the second equation gives:

$$ \frac{1}{p_{\{1\}} L_{\{1\}}} \frac{x}{x+z} = \frac{1}{p_{\{1,2\}} L_{\{1,2\}}} \frac{2z}{x+y+2z} = \frac{1}{p_{\{1,2\}} L_{\{1,2\}}} \frac{z}{x+z}, $$

which, after solving for $z$, gives $z = \frac{p_{\{1,2\}} L_{\{1,2\}}}{p_{\{1\}} L_{\{1\}}} x$. The solution to the system are of the form $x = (x, x, \frac{p_{\{1,2\}} L_{\{1,2\}}}{p_{\{1\}} L_{\{1\}}} x)$, for $x \in (0, \infty)$. We thus get for $\lambda_c$:

$$ \lambda_c = \frac{2}{\langle \ell_{\mathcal{D}} \rangle (p_{\{1\}} L_{\{1\}} + p_{\{2\}} L_{\{2\}} + 2p_{\{1,2\}} L_{\{1,2\}})}, $$

which is the value of the critical arrival rate for this network. If we generalize this result for $K \ge 2$, we can see that $z = p_C L_C$ is a solution for the system defined in Theorem 4.5. We thus get:

$$ \lambda^* p_C = \frac{1}{L_C} \frac{|C| p_C L_C}{\langle \ell \rangle \sum_U |U| p_U L_U}, $$
which, after some combinatorics, give:

$$ \lambda^* = \frac{K}{\langle l \rangle \sum_{U} |U| p_{U} L_{U}}, $$

which is the critical arrival rate for the symmetric network.

We can also verify this stability region through simulation: Figure 4.1 displays a non-symmetric BWP model with $K = 2$ bands where bands $\{1\}$ and $\{2\}$ are not considered symmetrically by the system: 80% of users transmitting on one band chose band $\{1\}$, and the remaining 20% chose band $\{2\}$. The rest of network parameters are the same as the ones used in Section 2.3.

<figure><img src="image_29.png" /><figcaption>Figure 4.1: Number of active users in the non-symmetric network over time with $K = 2$ bands for two values of $\lambda$. In this setup, $p_{\{1\}} = 0.64$, $p_{\{2\}} = 0.16$ and $p_{\{1,2\}} = 0.2$.</figcaption></figure>
# EXTENSION AND CONCLUSION

In the first part of this manuscript, we studied the stability of a multiclass SBD model with wireless-type interactions. As mentioned multiple times, these dynamics are leveraged by wireless networks, but in Chapter 4, we were able to obtain a general stability condition for systems that have fluid limits meeting certain requirements. In the chapter, we will discuss the results obtained in this part and interpret them.

## 5.1 Interpretation of the results

Theorem 2.2 gives the critical arrival rate $\lambda_c$ below which the dynamics is stable and above which it is not in the symmetric system:

$$ \lambda_c = \frac{K\ell(r)}{\mathcal{L}\langle\ell_D\rangle}. \qquad (5.1) $$

We can note that the form of this critical value in question is somewhat unexpected. The quantity $\mathcal{L} = \sum_U p_U |U| L_U$ appears, while a natural interpretation would have preferred the quantity $\sum_U p_U^{L_U/|U|}$, which denotes the average user load per band in the network. We can give the following interpretation for the value of $\mathcal{L}$: 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}, $$

Under this form and with the previous remark, 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). This way, the stability condition for the symmetric system can be seen as the classical
condition $\rho < 1$ seen in Jackson networks, for which the density of arriving jobs in the system cannot exceed its throughput.

For the non-symmetric system, we can not reach a closed form condition for the value of $\lambda_c$, though we can define it as:

$$ \lambda_c = \frac{1}{p_C L_C \langle \ell_D \rangle \sum_U |C \cap U| z_U} \qquad (5.2) $$

where $z \in \mathbb{R}^{2^K-1}$ is the solution of a system of equations given in Theorem 4.5. This results is a corollary of Theorem 4.2, which states that the vector $z$ is such that:

$$ z \in \arg\max_{x \in \mathbb{R}^{N_{\epsilon}}} \min_{i} \frac{1}{p_C L_C} \frac{|C| \ell(r) x_{i,C}}{\sum_{k,U} \ell^{\epsilon}(a_k, a_i) x_{k,U}}. $$

To interpret this result, we can see $z$ as the maximal capacity of the channel with the lowest throughput in the network. Due to the nature of the interactions in the system, whenever the population of a queue explodes, it propagates to the other queues, leading to the instability of the system. Once again, we obtain a "load-throughput" condition for the stability of the system in the non-symmetric system.

## 5.2 Restricting band configurations

In this work, we considered that all band configurations are allowed for users in the system. On top of allowing an exponential number of configurations for users, this setup is not commonly used for wireless networks. For instance, a usual setup for band allocation is to only allow the use of contiguous bands (this way, with $K = 4$, configurations $\{1, 2\}$ and $\{2, 3, 4\}$ are allowed but not $\{1, 3\}$). We will see in this section how to extend the previous results when restricting the band configurations, i.e., when allowing $p_C = 0$ for certain $C \in \mathcal{P}(K)$. Let us now denote $\mathcal{P}(K)^* = \{C \in \mathcal{P}(K), p_C \neq 0\}$ as the set of used band configurations in the network.

The main result we use to obtain stability in the system is stochastic monotonicity. In simple words, this property states that when we increase the population of one class of users, it slows down transmissions for all the other classes in the system. The second observation about the network we can make is the system of equations obtained in Theorem 4.5: this system needs to have a solution space of dimension 1 (i.e., of the form $x_0t$ for a given $x_0$ and $t > 0$) for the value of $\lambda_c$ to be uniquely defined.

If there exists a subset of classes "disconnected" from the system in the sense that they do not have an impact on the rate of transmissions, i.e. a subset $H \in \mathcal{P}(K)$ of classes $D$ such that for all $H \notin S$, $C \cap D = \emptyset$, it will act as a single-class SBD process independent from the rest of the dynamics, with its own stability condition.

Assuming the Poisson rain of arriving users stays common to all users, this leads us to obtain a condition on the form of the set $\{C \cap U, C, U \in \mathcal{P}(K)\}$.
This condition can be expressed in terms of the matrix $\mathcal{J} = (|C \cap U|)_{C,U \in \mathcal{P}(K)}$: if such classes $D$ exists, it implies that the matrix $\mathcal{J}$ can be rearranged to be a block matrix. Another way to describe it is to consider the graph $\mathcal{G}$ for which the vertices are labelled with the elements of $\mathcal{P}(K)$ and two vertices labelled with $U$ and $V$ are connected if and only if $U \cap V \neq \emptyset$. The conditions then becomes that the graph $\mathcal{G}$ is connected.

Collecting these observations, we can state the following Theorem:

**Theorem 5.1.** If the interaction graph $\mathcal{G}$ is connected (or, if the interaction matrix $\mathcal{J}$ cannot be rearranged in a block matrix), the critical arrival $\lambda_c$ remains unchanged when allowing the $p_C$ to be equal to 0.

If the interaction graph $\mathcal{G}$ is not connected (or, if the interaction matrix $\mathcal{J}$ can be rearranged in a block matrix), let $S_1, \dots, S_n$ be the $n$ sets labelling the connected components of $\mathcal{G}$ (or the blocks of $\mathcal{J}$), and let $\lambda_{c,1}, \dots, \lambda_{c,n}$ be the critical arrival rates obtained when considering separately each of the configuration sets $S_1, \dots, S_n$. The critical arrival rate for the whole system becomes:

$$ \lambda_c = \min_i \lambda_{c,i}. \qquad (5.3) $$

An interesting remark is that if we now assume that the Poisson rains for the arrival of users are independent for each of the possible configuration sets $S_1, \dots, S_n$, instead of taking the minimum critical rate for the system, the stability condition becomes:

$$ \lambda_c = \sum_i \lambda_{c,i}. $$

Using Theorem 5.1 and the previous remark allows to capture even more resource allocation policies for networks, and this result can be used for load balancing purposes, as one can be interested in maximizing the value of $\lambda_c$ by tuning the parameters of each of the multiclass SBD processes to desire.

Finally, the results presented in this part have been proved when assuming that the interfering power between channels was equal to the size of their intersection, but they can be extended to any similarity measure $f(C, U) > 0$ between channels $C$ and $U$ by using it instead of the intersection measure in all the previous equations.

## 5.3 Extension

The work presented in this part is meant to introduce a framework to study multiclass SBD processes. This framework is simple yet it captures some key features leveraged in bandwidth partitioning model. Future work may add more complex features to study a larger variety of networks.

A first direction to look into is the cellular version of the presented network. In the setup we presented in this part, we studied a device-to-device network (D2D). A natural extension would be to study a device-to-network
(D2N) version, where an antenna is located at the origin, and all the transmitter share this antenna as a receiver. In such a setup, if we assume that the antenna has an infinite channel capacity, the rate-of-transmission function becomes:

$$R(x, \Phi_t) = \frac{|C|\ell(\|x\|)}{\mathcal{N}_0 + I(x, \Phi_t)}$$

The behavior of D2N networks is known to be different than that of D2D networks. For instance, [3] shows a D2N version of the SBD process studied in [80], in which meta-stability appears: the system appears to be stable until a moment where the network population explodes. In this newer setup for example, it not sufficient to have $\ell'(x) \le \ell(x)$ anymore to obtain domination, but the other conditions from Theorem 1.2 still hold. Thus, a first extension of the problem would be to study the stability of this new system and see how the stability condition from Theorem 2.2 changes.

Among the other possible extensions, one can think of adding fading to the system, to model blocking or atmospheric phenomena that can alter the quality of transmission, and can add line-of-sight models in order to see how the stability condition changes.

This concludes the first part of this work, focusing on describing the stability region of this array of multiclass SBD processes. As positive recurrent Markov processes, these dynamics have a stationary distribution to which they converge. The next part will focus on the study on this stationary distribution as well as the result point processes and some of their properties.
# Part II

## Stationary interacting point processes
<figure><img src="image_30.png" /></figure>
# STUDY OF THE STEADY STATE DYNAMICS

The main result of Part I, Theorem 2.2, gives the existence of a stationary regime for the multiclass spatial birth-and-death process when $\lambda < \lambda_c$. The existence of a stationary regime of the dynamics raises interesting questions about the nature of this stationary regime and the properties of the resulting point process.

The first question we will tackle in this part is the stationary user densities in the system. We know that the arrival process is a Poisson rain in the system, and we also know the stochastic intensity of the departure process. Given this, we try to estimate the spatial user densities in the stationary regime.

The second question we will answer in this part is revolves around the geometry of the point process in the stationary regime. We know from [80] that the single class system show clustering, i.e. attraction in the stationary regime. In a clustered point process, points "attract" each other, which translates in a condition between the regular and the Palm expectation of shot noises of the point process.

Finally, the last question of interest revolves around the stationary distribution of the queuing process. The discretization of the multiclass SBD process introduced in Part I lead us to define two multiclass interference queuing networks (defined similarly as in [81]). We know from [22] that single-class interference queuing networks have exponential moments, which gives as a consequence bounds on the tails of the distribution. This property is interesting for ultra reliable, low latency (URLLC) networks ([62]), as it gives tight bounds for the tails of the distribution of latency in this system.