Mean-field limits for large-scale random-access networks

We establish mean-field limits for large-scale random-access networks with buffer dynamics and arbitrary interference graphs. While saturated-buffer scenarios have been widely investigated and yield useful throughput estimates for persistent sessions, they fail to capture the fluctuations in buffer contents over time, and provide no insight in the delay performance of flows with intermittent packet arrivals. Motivated by that issue, we explore in the present paper random-access networks with buffer dynamics, where flows with empty buffers refrain from competition for the medium. The occurrence of empty buffers thus results in a complex dynamic interaction between activity states and buffer contents, which severely complicates the performance analysis. Hence we focus on a many-sources regime where the total number of nodes grows large, which not only offers mathematical tractability but is also highly relevant with the densification of wireless networks as the Internet of Things emerges. We exploit time scale separation properties to prove that the properly scaled buffer occupancy process converges to the solution of a deterministic initial-value problem, and establish the existence and uniqueness of the associated fixed point. This approach simplifies the performance analysis of networks with huge numbers of nodes to a low-dimensional fixed-point calculation. For the case of a complete interference graph, we demonstrate asymptotic stability, provide a simple closed-form expression for the fixed point, and prove interchange of the mean-field and steady-state limits. This yields asymptotically exact approximations for key performance metrics, in particular the stationary buffer content and packet delay distributions. The methodological framework that we develop easily extends to various model refinements as will be illustrated by several examples.


Background and related work
Wireless networks are already large and complex today, and being at the heart of the so-called Internet of Things (IoT) [2], are expected to grow even denser in the future [21]. The rapid increase in wireless applications has led to an increasing demand for scarce wireless spectrum, and hence it is necessary to make the most efficient use of the limited available capacity. Obviously, when the number of nodes is large, in the hundreds or even thousands of nodes, a dedicated medium or channel cannot be assigned to each node, and nodes have to share the medium. Multiple simultaneous node transmissions on the same channel will however inevitably give rise to interference and loss of throughput. Medium access control (MAC) mechanisms are crucial to resolve this contention. However, in large networks, a centralized control mechanism is hard to implement and to maintain since it would require constant status updates generating prohibitive communication overhead. For this reason the design of efficient distributed (local) MAC protocols has attracted a lot of attention.
A very popular distributed MAC mechanism is the CSMA (Carrier-Sense Multiple-Access) protocol, which is currently at the core of the IEEE 802.11 and 802.15.4 standards. Its popularity is mostly due to its simplicity and efficiency. The key feature of the CSMA protocol is that each node waits for a random back-off period before initiating a transmission. Interference is avoided since the back-off countdown is interrupted whenever potential interference is sensed, and only resumed once the medium is sensed idle again. This protocol, whilst extremely easy to understand on a local level, generates complex and interesting macroscopic network dynamics.
In the performance analysis of CSMA networks, a common assumption is the existence of an underlying graph that represents interference between the various nodes in the network. An edge between two nodes means that destructive interference is caused by simultaneous transmission. Both empirical and theoretical support for the notion of an interference graph is provided in [26,46].
When the nodes always have packets to transmit, the network is said to be saturated and the macroscopic activity behavior is amenable to analysis under the assumption of an interference graph. In particular, the activity process has an elegant product-form stationary distribution [9,29,40]. The computation of the stationary distribution of the activity process reduces to the identification of all the subsets of nodes which may transmit simultaneously, namely the independent sets of the interference graph.
Real-life scenarios however involve unsaturated networks. Packets arrive at the various nodes according to exogenous processes, and buffers may drain from time to time as packets are transmitted. In particular, in IoT applications, sources are likely to generate packets only sporadically, with fairly tight delay constraints, and often have empty buffers. Since empty nodes temporarily refrain from the medium competition, the activity process is strictly intertwined with the buffer content process. In this situation, the product-form solution no longer holds [14,40] and an exact stationary analysis does not seem tractable. Furthermore, not even stability conditions for the queueing dynamics in the nodes' buffers are known, let alone results for the stationary distribution of the buffer content process.
The analysis of unsaturated CSMA networks simplifies if certain symmetry conditions amongst the various nodes hold. An important instance is when there is a substantial number of nodes with similar traffic and placement in the network, so that the operation of one is equivalent to that of many others. More generally, nodes can be divided into classes with the symmetry conditions now applying to nodes of the same class. The asymptotic regime where the number of such nodes in each class grows to infinity, is often referred to as a mean-field regime. Mean-field theory originated in physics, where it is still largely used in the study of models involving a large number of interacting particles. The aggregated effect of all the other particles on any tagged particle is approximated by a single averaged effect (the mean-field), thus reducing a many-body problem to a tractable one-body problem. Mean-field analysis has since been exploited in various scientific areas, from epidemic models to queueing systems and from biology to game theory. In the context of random-access networks a mean-field regime, not only provides analytical tractability, but is also highly relevant in the context of the envisioned massive numbers of IoT devices.
A thorough survey of mean-field analysis of random-access protocols is presented in [19]. The work of Bianchi [6] is a landmark paper which assumed a propagation of chaos property to hold so as to derive tractable formulae for the key performance measures of the system. The papers surveyed in [19] mostly use mean-field theory to provide either evidence or objection for the propagation of chaos assumption of Bianchi. Among these papers, it is worth mentioning [16], where the authors investigated the existence of a global attractor for the mean-field system and provided sufficient conditions for its existence, deducing the validity of Bianchi's assumption. Further papers which deserve to be mentioned are [10,36], where the authors exploited mean-field theory so as to obtain approximations for key performance measures of large systems. In particular [10] focuses on the characterization of the stability region, while [36] examines the throughput performance of the system. None of the above-mentioned papers considered scenarios with unsaturated buffers, with the exception of [10], which however dealt with systems evolving in discrete time and did not consider performance metrics like packet delays.

Key contributions and paper organization
In the present paper we examine the buffer dynamics in large-scale unsaturated random-access networks. Specifically, we analyze the buffer occupancy processes in a mean-field regime where the number of nodes grows large. Such a mean-field regime not only furnishes analytical tractability, but is also particularly relevant in understanding the collective behavior of massive numbers of interacting IoT sources.
The paper has two main interrelated threads. The first part in Section 4 deals with networks with general interference graphs. We prove that a suitably scaled version of the buffer occupancy processes converges in the mean-field limit to a tractable deterministic initial-value problem. We also establish necessary and sufficient conditions for the existence and uniqueness of a fixed point of the initial-value problem, and provide a characterization the fixed point as the solution of a low-dimensional equation.
The second part focuses on scenarios with a complete interference graph. We demonstrate in Section 5 global asymptotic stability of the initial-value problem, i.e., convergence to the unique fixed point from any initial state with finite mass, and then proceed to show positive recurrence of the pre-limit process and tightness. We combine these properties to prove interchange of the mean-field and steady-state limits. In Section 6 the interchange of limits is then leveraged to establish that the stationary buffer content distributions at the various nodes converge to geometric distributions while the stationary distributions of the scaled waiting time and sojourn time converge to exponential distributions. The parameters of these limiting distributions are directly expressed in terms of the fixed point of the initial-value problem. These results provide asymptotically exact approximations for the stationary waiting-time and sojourn time distributions, which in fact turn out to be remarkably accurate, even for a relatively moderate number of nodes.
Mean-field limits for various model extensions are briefly discussed in Section 7, illustrating the versatility of the methodology. In Section 8 we present some simulation experiments to illustrate the analytical results. Finally, in Section 9 we make a few concluding remarks and offer several suggestions for further research.

Model description
We consider a network of N nodes sharing a wireless medium according to a random medium access protocol. The various nodes may be grouped into a set of classes/clusters C = {1, . . . , C} such that nodes in the same class have the same statistical characteristics. Denote by N . Once a class-c node obtains access to the medium, it transmits one packet, which takes an exponentially distributed time with parameter µ (N ) c= µ c . In between two consecutive transmissions a class-c node must back-off for an exponentially distributed time period with parameter ν Note that a node refrains from the back-off competition whenever its buffer is empty. Also, the back-off period of a class-c node is suspended (frozen) when the medium is occupied by an interfering node, i.e., the class activity state does not belong to Ω −c .
Define the queue length process c,k (t) represents the number of packets in the buffer of the k-th node belonging to class c at time t excluding any possible packet in transmission, and define by Y (N ) (t) the class activity process on the state space Ω. Observe that Q (N ) (t), Y (N ) (t) evolves as a Markov process.
A population process description. To ease the notation we write N c = N (N ) c from now on. Due to the class symmetry, the nodes within a class are statistically exchangeable, and the system state may be described in terms of the numbers of nodes belonging to the same class and with the same number of packets in the buffer. In particular, define the population processX Observe that the population process lies in E 1 ⊆ E where x c,n = 1, ∀c ∈ C}. (2) The possible transitions for the process X (N ) (t), Y (N ) (t) may be described as follows: • A packet arrives at a class-c node having n packets in its buffer: this happens at rate λ (N ) c times the number of class-c nodes in state n, i.e., λ c,n , and generates the transition where e c,n ∈ E 1 has all 0 entries except a 1 in position (c, n).
• A transmission is completed by a class-c node: this happens at rate µ (N ) c = µ c and only if a class-c node is transmitting, i.e., Y (N ) ∈ Ω +c . The transition generated is the following • A back-off is completed by a class-c node having n > 0 packets in its buffer: this happens only if the class activity state allows class-c nodes to back-off, i.e., Y (N ) ∈ Ω −c , and at rate ν (N ) c times the number of class-c nodes in state n, i.e., ν c,n , and generates the transition Preliminary results for saturated scenario. As mentioned earlier, we focus on an unsaturated network where nodes with empty buffers refrain from competition for the medium. For later purposes, it is convenient to also consider a related saturated network where each class behaves as a single node always competing for the medium, and backing-off at rate η c ν c , where η = (η c ) c∈C ∈ R C + . The transmission times of the node associated with class c are exponentially distributed with parameter µ c .
Define now the throughput function θ : represents the fraction of time that the node associated with class c is active in the saturated network. As proved in [27,41], the throughput map θ is globally invertible, i.e., for any achievable throughput vector γ ∈ int(Γ) there exists a unique vector η(γ) ∈ R C + such that θ(η(γ)) = γ.
Intuitively, the c-th coordinate of η(γ) represents by what factor the back-off rate at node c needs to accelerate/decelerate in order for the target throughput vector γ to be achieved in stationarity.

Overview of the main results
In this section we provide an overview of the main results of the paper, along with an interpretation and high-level discussion of their ramifications, before presenting the proofs in the subsequent sections.
In the first part, we analyze networks with general class-based interference graphs. We consider the sequence of population processesX (N ) (t) as the number of nodes in the system increases, and rigorously establish its weak convergence to x(t), the solution of a tractable deterministic initial-value problem as stated in the next theorem.
The major difficulty in the analysis of unsaturated networks arises from the correlation between the population processX (N ) (t) and the activity process Y (N ) (t). A key observation in the proof of Theorem 1 is that these processes evolve on different time scales, i.e., the population process evolves N times slower than the activity process. Hence, in the mean-field regime, the rapidly changing activity process "converges to an instantaneous measure on the activity states" determined by the current fraction of queues which are empty. This is the reason why the initial-value problem (6) does not involve explicitly the activity process and the evolution of the population process at time t depends only on x(t). Specifically, fraction of class-c nodes with n ≥ 1 packets in the buffer "increases" at rate λ c and "decreases" at rate ν c π x0 (Ω −c ), where the measure π x0 (Ω −c ) represents the limiting "instantaneous measure" on the activity states that allows class-c nodes to back-off. This argument based on a "separation of time-scales" is highly flexible and applies to more general and complicated versions of random-access networks as will be illustrated in Section 7.
The initial-value problem (6) is certainly easier to analyze than the population process in a network with a finite number of nodes. Theorem 2 states that, under certain necessary and sufficient conditions, there exists a unique equilibrium point x * for the initial-value problem (6), and provides its characterization in terms of the load of the network.
is the unique equilibrium point with H(x * ) = 0 in E 1 of the initial-value problem (6). Specifically, the condition ρ ∈ int(Γ) is necessary and ensures the load ρ of the system is sustainable.
On the other hand, the condition ξ = η(ρ) < e ensures that the desired throughput vector ρ can be achieved by a feasible back-off vector. Although Theorem 2 only concerns the scaled population process, combined with Theorem 1 it yields that the stationary distribution of the associated class activity process is given by π x * 0 (·) in the limit. In other words, in the limit the stationary distribution of the class activity process is the same as in a scenario where the aggregate back-off rate of class c is a constant fraction ξ c of the nominal back-off rate ν c . This is consistent with the fact that ξ c is the stationary fraction of class-c nodes that have non-empty buffers and compete for the medium, and ξ will therefore in the sequel be referred to as the vector of activity factors. The equation ξ = η(ρ) reflects that the activity factors must be such that each class c is active a fraction ρ c of the time.
The second part of the paper focuses on deriving asymptotically exact results for the key performance measures of the system. For a rigorous treatment, we focus on the scenario where all the classes mutually interfere, i.e., the interference graph G = (C, E) is complete, in which case In this framework, the nodes of every class are allowed to back-off only when the activity process is in the idle state. We will show that the vector of activity factors ξ in Theorem 2 simplifies to so that the condition ξ = η(ρ) < e can be expressed in explicit form as which forces the right-hand side to be positive, i.e., ρ ∈ int(Γ). When condition (8) applies, Theorem 3 holds and asymptotically characterizes the stationary distribution of the population process as the equilibrium point of the initial-value problem (6).
Theorem 3. The sequence of stationary distributions (x (N ) ) N ∈N ⊆ E 1 converges to x * ∈ E 1 , i.e., In order to prove Theorem 3, we show that the interchange of limits displayed in Figure 2 holds. Specifically, our methodology reduces the derivation of the stationary distribution of an intractable N -dimensional Markov process to the computation of the equilibrium point of a low-dimensional fixed-point equation.
Theorem 3 is exploited so as to obtain approximations for the performance measures of the system. In particular, denote by Q 3. In Section 4.1.3 we characterize the limiting process (11). 4. As a last step, in Section 4.1.4, we describe how the limiting process α(t) at time t is uniquely determined by the value x(t). Hence a self-contained deterministic initial-value problem governing the behavior of x(t) is obtained.
In preparation for the analysis, we briefly recall few useful properties of the topological space where the population process evolves. The sample paths of X (N ) (t) lie in D E [0, ∞), i.e., the set of the cadlag functions from [0, ∞) in E = χ C , where the space E ⊆ c∈C R ∞ is defined in (1). Define on E the metric ρ(·, ·) such that where x c = (x c,n ) n∈N and ρ 1 is a metric on R ∞ , i.e., Note that (R ∞ , ρ 1 ) is separable and complete, and χ is compact in (R ∞ , ρ 1 ) as each coordinate lies in [0, 1], see [7, page 219]. From [8, Section M6, pages 240-241], the following lemma is therefore obtained.
Lemma 1. The subset E = χ C is complete, separable, and compact under the product topology induced by the metric ρ in the space c∈C R ∞ .
In the subsequent analysis we will work with D E [0, ∞) using the metric d, an exponentially weighted version of the Skorohod metric, see [20, page 117]. Under this metric D E [0, ∞) is complete and separable as (E, ρ) is itself complete and separable, see [20, Theorem 5.6, page 121].

Unit Poisson process representation
In this subsection we describe the fluid-time marginals of the prelimit Markov population process X (N ) , Y (N ) in terms of infinite sequences of unit Poisson processes. For each c ∈ C, and then for each n ∈ N 0 , the random quantity N A,c,n (t) determines the class-c arrivals to queues with n packets in the buffer during the time interval [0, t). We similarly define N B,c,n (t) for each c ∈ C, n ∈ N, to determine the back-offs and finally N T,c (t) determines the process of transmissions for each c ∈ C. These processes are supposed mutually independent and defined on a common probability space Ω F , F F , P . We also define the centered versions of these processes to beÑ A,c,n (s) = N A,c,n (s) − s, c ∈ C, n ∈ N 0 , with corresponding definitions for the back-offs and the transmissions. These are martingales with respect to their natural filtrations.
To obtain the Poisson process representation, we first fix the initial conditions, which are deterministic. Realizations of the above unit Poisson processes are then drawn and the (unscaled) population random variablesX, Y are obtained as solutions to the following equations, where fluid time is t. That is, and e d,m c,n = e d,m − e c,n with e c,n ∈ E denoting the vector having zero entries everywhere but a 1 in its (c, n) coordinate.
As in the discussion in [32, Section 2.2], the solutions to these equations can be obtained by fixing the sample paths of the Poisson processes N A,c,n , N B,c,n , N T,c above. It can then be seen inductively that the sample paths of the process (X (N ) (t), Y (N ) (t)) lie in D E 1 ×Ω . Note that the total rate of the above Poisson processes is bounded above by c∈C (λ c + ν c + µ c ) < ∞. Hence there can only be finitely many jumps in any finite interval [0, t] w.p. 1 and the sample paths are by construction piecewise constant cadlag functions.
To confirm that the distribution corresponds to that of the multi-dimensional process defined in Section 2, let G t denote the history up to time t, that is It is straightforward to show that the infinitesimal transition rates obtained by conditioning on G t and the current state, yield those of the multi-dimensional process detailed in Section 2, see also [32]. But these conditional rates are the defining property of the process and so the representation follows the same statistical law. Equivalently, as the rates are bounded, the representation could also be verified by considering the uniformized version of the process and then showing that the underlying jump chain matches that of Section 2. We now rewrite these equations in a more compact form, working componentwise. The N superscript has been omitted for conciseness. where for each n ∈ N 0 and c ∈ C. Equation (12) expresses the change in the fraction of class-c nodes with n packets, i.e., arrivals to a queue with n − 1 packets or departures from a queue with n + 1 packets increment this fraction whereas arrivals or departures to component n result in a decrement. Equation (13) describes the dynamics of the class activity process and Y c is incremented whenever there is a back-off from any non-empty class-c queue, and returns to inactivity when the corresponding packet has been transmitted. Equation (14) defines the sequence of arrivals by fluid time t, the (stochastic) intensity is proportional to the fraction of queues for each component, with the convention A c,−1 (t) ≡ 0 for all t. Equation (15) expresses the back-offs as an integral of a previsible process over a Poisson process, again with intensity varying according to the fraction of the corresponding component. Note that D c,0 (t) = 0 as no back-offs can occur from empty queues. Finally the transmission process (16) is again an integral of a previsible process over a scaled version of the original Poisson process.
The expressions (12)- (16) are preferably written in the following martingale form: where The following technical lemma guarantees that the processes M A,n , M B,n , M T , and Ξ are martingales, a property that we will use throughout the following analysis. The proof of the lemma is given in Appendix B.1. For ease of exposition we will consider the single-class case only, but the extension of the argument to the multi-class case is straightforward.
Lemma 2. Given the common probability space Ω F , F F , P , there is a filtration F defined on it, such that M A,n , n ∈ N 0 and M B,n , n ∈ N, are sequences of locally square integrable martingales adapted to F together with M T which is also a locally square integrable F-martingale. In addition Ξ is a L 2 -F martingale.
Observe now equations (17) and (18). Note that on a fluid time scale there are order N transitions in the activity process, which thus evolves much faster than the fluid population process. Hence, in the mean-field regime, the activity process reaches its stationary distribution before the state of the population process changes. Intuitively, the stationary distribution of the activity process is influenced by the current state of the population process, in particular only a fraction (1−x c,0 ) ∈ [0, 1] of class-c nodes has packets and so competes for the medium access. Hence, using the notation introduced in (3), given that the population process is in state x ∈ E, the probability for the activity process to be in state ω ∈ Ω in stationarity is π x0 (ω)= π(ω; e − x 0 ).

Representing Y as a random measure
Since the weak limit of the sequence Y (N ) does not exist in D, we introduce the cumulative time process α (N ) ∈ C R C + ,↑ [0, ∞), which has components α (N ) for ω ∈ Ω, as the cumulative fluid time spent in state ω by the activity process Y (N ) in the interval [0, t]. That is, for each ω ∈ Ω, t ≥ 0, Observe that ω α (N ) is a continuous, increasing and unit Lipschitz function for every ω ∈ Ω. Indeed, α (N ) is unit Lipschitz under the sup norm metric.
Thus Y (N ) induces a sequence of random measures γ (N ) , as determined by α (N ) ∈ C R C + ,↑ [0, ∞), on the space of increasing, nonnegative functions with initial value 0. This holds since the mapping (24) is continuous. For any realization, α (N ) the measure on B([0, ∞) × Ω) is obtained by integrating the derivatives of α (N ) . In fact, α We may therefore rewrite (17) by substituting the left derivatives of α (N ) in the intensity function (compensator), We now proceed to show that a weak limit for the sequence of joint processes exists. This follows from Proposition 1, whose proof is given in Appendix B.2.
Observe that, as a consequence of [20, Theorem 2.4, page 107] (see also [7, Exercise 6, page 41]), in order to prove Proposition 1 it suffices to show the tightness for the marginals, i.e., the tightness of α (N ) and of X (N ) . The proof of the tightness of α (N ) is quite simple, but establishing the tightness of X (N ) requires substantial effort. In particular, the argument in Appendix B.2 relies on Lemma 1 and on the sequence of processes X (N ) being asymptotically Lipschitz in probability, an easy-to-verify property which we detail separately in Appendix A.

Mean-field limit characterization
We just established the existence of (x(t), α(t)), i.e., a weak limit for the sequence X (N ) (t), α (N ) (t) . In this section, as a further step towards Theorem 1 we prove an intermediate mean-field limit, which will be further developed in the next section. Specifically we describe the steps needed to obtain the following result.

Proposition 2.
Given (x, α) a weak limit for the sequence (X (N ) , α (N ) ) then x and α satisfy the differential equation Observe that Equation (25) determines X (N ) (t) as in the following sum, of an initial term plus sample paths lying in D: where I (N ) Furthermore if x, α are weak limits of X (N ) , α (N ) then we make corresponding definitions for B (t) and let I be the corresponding limit. In order to prove Proposition 2, we need to go through the following steps: (a) Show that the weak limits of X (N ) (t) and I (N ) (t) coincide. This is established through Proposition 3 and the following remark.
(b) Derive the weak limit of I (N ) (t) as in Proposition 4 and Lemma 3.
Focus on step (a), we apply the Continuous Mapping Theorem, see [7, Theorem 5.1, page 30] to the weak limits for the terms on the right hand side of (26) as follows. First recall [8, Theorem 3.1, on some common probability space) and taking values in a metric space with metric d then if I (N ) ⇒ I and it also holds that d(I (N ) , X (N ) ) ⇒ 0, then it follows that X (N ) ⇒ I. To apply the above result to (26), we will take d to be the Skorohod metric on D E [0, ∞) as defined in [20,Section 3.5]. In order to show that X (N ) and I (N ) have the same weak limit, it will be enough to establish the following proposition which is proved in Appendix B.3.
We now continue with step (b). Observe that I (N ) (t) is a sequence of continuous paths. We will establish convergence of this latter sequence in C E [0, ∞) under the uniform metric using the Continuous Mapping Theorem. This establishes convergence in D (the Skorohod topology relativized to C coincides with the uniform topology).
Given a subsequence (X (Nm) , α (Nm) ) ⇒ (x, α) (without loss of generality we will take this subsequence to be the whole sequence), we obtain the corresponding weak limit of I (N ) provided that the given mapping is continuous. and The proof of Proposition 4 is given in Appendix B.4. Note that X is Lebesgue integrable over any finite interval and so the first integral is well-defined. Since the α (N ) ω are unit Lipschitz and increasing, the integrals in (28) also exist. Proposition 4 only establishes pointwise convergence to the continuous functions on the right hand side of (27) and of (28). However since these functions are non-negative, non-decreasing as well as unit Lipschitz continuous, then given any fixed T > 0 it follows that convergence is uniform in the interval [0, T ], for each component (c, n) as a consequence of the following lemma which is proved in Appendix B.5.
However uniform convergence of the components implies the same for the sequence so that

Weak limit characterization of α
Proposition 2 determines any weak limit x as the solution to a differential equation for which the corresponding limit occupancy measure α is given. We now proceed to characterize this latter limit in terms of the stationary measure of the activity process.
Recall the properties satisfied by the processes α (N ) ω . In particular, for every t ≥ 0, We begin by substituting α into (18) to obtain that is an F-martingale. Dividing by N and taking weak limits for (x, α), we obtain that for every t ≥ 0. We apply to (30) the Continuous Mapping Theorem as explained in the discussion before Proposition 3, and use the limits established in Proposition 4, to obtain for every t ≥ 0.
Proof. Thanks to the fundamental theorem of calculus applied to (31), it holds that Hence, at every time t, we have to look for (z t ω ) ω∈Ω , solutions of where in view of (29) Observe that, z t ω corresponds to the stationary fraction of time spent by the activity process in state ω ∈ Ω in the saturated network where node c transmits at rate µ c and back-offs at rate ν c (1−x c,0 (t)), i.e.,α which is equivalent to (32). In particular, existence and uniqueness of the stationary distribution Hence,α ω (t) = z t ω is well-defined and, since an increasing unit Lipschitz function starting from 0 is determined by its derivatives, the equality holds for all t, determining α ω . Observe that α ω exists everywhere since it coincides with the integral of a continuous function.
Having characterized the weak limit of α, we obtain Theorem 1 by substituting (32) in Proposition 2.

Analysis of the mean-field limit
In the previous section we proved Theorem 1, establishing the convergence of X (N ) (t) to x(t) ∈ C E [0, ∞), the solution of the initial-value problem (6). In this section we present the proof of Theorem 2. Specifically, we show that a solution to (6) exists and is unique. Moreover, we provide a general condition yielding the existence of a unique equilibrium point x * .
Transient behavior. Observe that by definition X (N ) (t) ∈ E 1 for every N ∈ N and t ≥ 0. Such property is carried over as N → ∞, indeed if x(0) ∈ E 1 , then x(t) ∈ E 1 for every t ≥ 0. Hence, the set E 1 is positive invariant for the initial-value problem (6). It is easy to show that the function H(·) defined in (6) is Lipschitz continuous in the metric space (E 1 , ρ), so that uniqueness and existence of the solution of (6) follows.
Limiting behavior. In order to prove Theorem 2, we observe that given that an equilibrium point x * exists, it must satisfy the relation H(x * ) = 0 component-wise, i.e., for every c ∈ C for every n ∈ N. In other words, any equilibrium point x * must satisfy the above geometric relation.
In addition, such an equilibrium point x * lies in E 1 if and only if ξ c = 1 − x * c,0 < 1 for every c ∈ C, i.e., ξ = e − x * 0 < e. Hence, thanks to (35), the uniqueness of ξ yields the same for x * . According to the product-form solution in (3) and the definition of throughput (4), it holds that In order for x * to be an equilibrium point, ξ must solve i.e., the fraction of time that class c is active in equilibrium must be equal to its load. Noting that ρ ∈ int(Γ), the global invertibility of the throughput map (5) implies that ξ has to satisfy Observe that the existence and uniqueness of ξ follows from the global invertibility property of the map θ(·). We conclude that x * as in (7) is the unique equilibrium point of the initial-value problem (6), and that x * ∈ E 1 if and only if η(ρ) < e. (As a side-remark, we mention that the initial-value problem would not have any equilibrium point, if the assumption ρ ∈ int(Γ) were not satisfied. This makes sense since the latter condition is necessary for the queue lengths to be stable.) As an example, consider the square network displayed in Figure 1 and assume that i.e., ρ ∈ int(Γ). In order to identify the unique equilibrium point of the initial-value problem (6), the following system of equations has to be solved: where The equilibrium point exists if and only if ξ < e. Methods for solving these kinds of equations are discussed more extensively in [15,41].
In this section, we proved the convergence of the population process towards the solution of a deterministic initial-value problem. Necessary and sufficient conditions for the existence and uniqueness of an equilibrium point have been established. In the following section we aim to exploit the equilibrium point so as to obtain an asymptotically exact approximation for the stationary performance measures of the system. The analysis will focus on the case of a complete interference graph where all classes mutually interfere.

The complete interference graph case
The second part of the paper focuses on the analysis of important stationary performance measures such as the queue length, the waiting-time, and the sojourn-time distributions. In particular, we will build upon the equilibrium point x * described in Theorem 2 to obtain approximations which are asymptotically exact in the mean-field regime for the case of a complete interference graph and accurate even in moderately large networks. Specifically, we will prove Theorem 3, i.e., we will show that for any initial state x ∞ ∈ E 1 with finite mass the solution x(t) of the initial value problem satisfies Due to the class symmetry, relation (37) implies that In order to establish (37), we need several steps.
• In Section 5.1 we will show that • Equation (37) is well-posed only if the population process in the finite system possesses a stationary distribution, hence conditions for the positive recurrence of the prelimit processes are established in Section 5.2.
• In Section 5.3 the interchange of limits is then proved so that (37) follows.
Throughout Sections 5 and 6 we assume the interference graph G to be complete and condition (8) to hold. Simulation results suggest that (37) remains valid even in the case of a general interference graph. In the course of the analysis we will point out the difficulties that would need to be resolved in order to obtain a rigorous proof in that case.
For a complete interference graph, it holds that Thus the equilibrium point x * has the following closed-form expression, and (8) is indeed necessary and sufficient so as to guarantee ξ < e.
Corollary 2. The unique equilibrium point in E 1 of the initial-value problem (6) is given by

Global stability of x *
We begin by observing that, when the interference graph is complete, Theorem 1 specializes as follows.
where the function H(·) is defined by x c,n e c,n+1 x c,n e c,n−1 c,n , and π b x0 = 1 1+ c∈C νc µc (1−xc,0) . In this section we will show that for any initial state x ∞ ∈ E 1 with finite mass the solution x(t) of the initial value problem satisfies i.e., the fixed point x * ∈ E 1 as described in Corollary 2 is a stable equilibrium point of the initialvalue problem described in Here we briefly outline the steps necessary in order to prove (40): • In Section 5.1.1 we introduce the concept of mass of a population vector as the average number of packets in the buffer of the various nodes. Observe that the mass of the equilibrium point x * is finite and the mass of a solution x(t) of (39) is a continuous Lipschitz function.
• Then, in Section 5.1.2, we show that stochastic dominance is an invariant property for the initial-value problem (39), see Proposition 5. This property allows us to couple different solutions of (39) and will play a key role in the proof of (40).
• We leverage the previous steps in Section 5.1.3 so as to show that, if x(0) has finite mass, x(t) → x * . This is proved in Proposition 6 and implies that x(0) ∈ B(x * ). We will say that x * is a globally stable equilibrium point. Given a solution x(t) ∈ C E 1 [0, ∞) of (39) such that m c (x(0)) < ∞, we observe that m c (x(t)) is differentiable in t and the following differential equation holdṡ Hence, the evolution of m(x(t)) is governed only by the value of x 0 (t), i.e., by the fraction of nodes with empty buffers. Moreover, observe thatṁ(x(t)) = 0 if and only if , ∀ c ∈ C.
Furthermore, observe thatṁ c (x(t)) is increasing in x c,0 (at least as long as x ∈ E 1 ). In fact Therefore, given two solutions x L (t) and

Stochastic dominance
The set E 1 ⊆ E may be partially ordered via the stochastic dominance relation.
x U c,n , for every m ∈ N 0 , c ∈ C.
We denote this relation by x L ≤ s x U .
This relation allows us to couple different solutions of the initial-value problem governed by (39) as stated in the following proposition, whose proof is provided in Appendix B.6.
Proposition 5. The stochastic dominance is a positive invariant relation for the initial-value problem governed by (39). In particular, given solutions Hence, by exploiting that if x(0) = x * then x(t) = x * for every t ≥ 0, we identify the following invariant sets.

Corollary 4. The sets
are positive invariant for the initial-value problem (39).
Moreover, stochastic dominance provides valuable information about how the mass of a certain solution evolves.
for every c ∈ C.
The first relation is a direct consequence of the mass definition, the second follows from the formula for the mass derivative (41) and from x L c,0 (t) ≥ x U c,0 (t) for c ∈ C, which is a direct implication of the stochastic dominance property. In particular, observe that if x(0) ∈ D (resp. x(0) ∈ D) then m c (x(t)) ≤ m c (x * ) (resp. m c (x(t)) ≥ m c (x * )) andṁ c (x(t)) ≥ 0 (resp.ṁ c (x(t)) ≤ 0) for every c ∈ C and t ≥ 0.

Proposition 6. The basin of attraction
Observe that the stated condition is necessary, indeed m(x * ) is finite and a solution of (39) cannot hit E 1 <∞ in finite time starting from an infinite-mass state. However, apart from this necessary exclusion, all the other solutions are attracted by the equilibrium point x * ∈ E 1 , which is therefore globally stable in E 1 <∞ . The same argument, without the assumption for the interference graph to be complete, would not lead to the proof of Proposition 6. Indeed the argument relies on Proposition 5 which is not true in general, i.e, given solutions x L (t), x U (t) of the initial-value problem governed by (6) for general class-interference graphs However, in the inductive argument used in the proof of Proposition 6, we rely on the invariance of the stochastic dominance only in the basic step, i.e., Lemma 5 in Appendix B.7, where we show that x 0 (t) → x * 0 . Hence the following proposition holds.
Proposition 7. Given a general class-interference graph, and assume the existence of an equilibrium point for the initial-value problem (6). Then, if a solution x(t) of (6) is such that Propagation of chaos and consequences of Proposition 6. In dealing with many-particle interaction problems, a common assumption is the propagation of chaos, i.e., the assumption that the particles become pairwise independent when their number grows and their evolution is influenced only by an aggregated averaged effect of the entire population state. Mean-field theory provides a powerful machinery by which the propagation of chaos assumption can be validated. In our setting, the particles are the nodes and their state is given by the number of packets in their buffer. Consider a chaotic initial condition, i.e., the initial number of packets in the buffers of the various nodes are mutually independent, then the multi-exchangeability property [25, Section 5] and the convergence towards the initial-value problem presented in Theorem 1 suffice to guarantee that the independence propagates over any finite horizon, see [3,39]. However, as described in [3], in order for the independence to hold in the limit (as the time grows to infinity), it is necessary that the solution of the initial-value problem (39) converges to an equilibrium point. We conclude that, in a complete interference graph setting, the queues of the various nodes in the system behave independently and undergo a common averaged effect depending on the equilibrium configuration of the network.
In the case of slotted CSMA networks, Bianchi in [6] assumed chaos propagates provided the network is large enough, and so obtained an analytical expression for the network throughput. Many papers discuss this assumption, including [16], which validated the propagation of chaos under certain assumptions via a mean-field approach. Here we just provided another supporting argument for the case of a complete interference graph.

Positive recurrence
As mentioned earlier, the final objective of this section is to obtain an approximation for the stationary distribution of the queue length processes at the various nodes in large systems. In particular, we will show that x * ∈ E 1 provides an approximation for the stationary distribution x (N ) of the system with N nodes when N is large. First, however, it must be ensured that the population process has a stationary distribution x (N ) under condition (8) for each system with a finite number of N nodes. Hence, we will establish the positive recurrence of the queue length process Q (N ) , from which the positive recurrence of X (N ) follows.
The complete interference graph assumption plays an important role in this analysis. Indeed, since at every moment in time at most one node can be active, the CSMA model can be described as a single-server polling system with random routing. Specifically, the switchover periods of the server correspond to the idle-state of the activity process, i.e., the intervals in which each node is allowed to back-off. When a node completes its back-off period, the server visits that node and serves a packet, when present. Specifically, consider the polling system with the following features: • Packets arrive at the k-th node belonging to class c as a Poisson process of rate λ (N ) c .
• The service policy is 1−limited. In particular, a single packet is served if present, otherwise the server leaves the node immediately; • When a class-c node is visited and its buffer is nonempty, a packet is served for an exponentially distributed time with parameter µ (N ) c ; • The switchover time of the server between two consecutive visits is exponentially distributed with parameter ν : • At the end of a switchover period, the probability for the server to pick the k-th node belonging to class c is equal to p c,k = ν Observe that the above-described polling system behaves exactly as the CSMA model with complete interference and that the latter property is crucial in establishing this connection. Polling systems have been extensively studied in the literature. In particular, the condition for the positive recurrence of the above-described model has already been investigated. In [17,18,22], the authors provide positive recurrence conditions for general classes of polling systems via the fluid limit analysis of the queue length processes [17,18] and via a direct approach [22]. The following result is a special case of Theorem 2.3 in [18], and in Appendix B.8 we adapted our notation to the setting in [18], so as to further illuminate the connection.
(ii) If > 1, then the queue length process Q (N ) is transient for every N ∈ N.
The proof in [18, Theorem 2.3] consists of analyzing the fluid limit of a polling system, and we can conclude the positive recurrence of the CSMA model due to the above-described equivalence.
In Appendix B.9 we provide an alternative and constructive proof of Proposition 8 for the singleclass case. The proof consists in exhibiting explicitly a quadratic Lyapunov function.
When the interference graph is not complete, the connection with polling systems does not apply, and an explicit condition for positive recurrence is an open problem.

Tightness and interchange of limits
In Section 5.2 we showed that the queue length process Q (N ) (t) has a stationary distribution if condition (8) is satisfied. Denote it by the probability measure π (N ) , i.e., lim t→∞ P{Q (N ) (t) = q} = π (N ) (q), q ∈ N N 0 .
For finite N the process Q (N ) (t) completely characterizes X (N ) (t), hence the population process X (N ) (t) also has a stationary distribution. In particular, define x (N ) ∈ E 1 as where e(q) ∈ E 1 is such that Then we know that lim At this point we are close to the target we set at the start of this section, that is proving equation (37). Indeed, in Section 5.1 we characterized x * ∈ E 1 , and thanks to Proposition 8, we proved the existence of the sequence of limiting distributions {x (N ) } N ∈N ∈ E 1 .
Two further steps are needed in order to show (37). First, we need to actually show that the sequence (x (N ) ) N ∈N ⊆ E 1 has at least one limiting point. Then we need to show that each possible limiting point necessarily coincides with x * . Once we settle these issues, we can finally state that the limits can indeed be interchanged and that equation (37) holds.
The proof of the following proposition is given in Appendix B.10.
At this point, via Prohorov's theorem [7, Theorem 6.2, page 37] we deduce that the sequence (x (N ) ) N ∈N ⊆ E 1 possesses converging subsequences. In order to complete the proof of Theorem 3, it suffices to show that these subsequences all converge to x * ∈ E 1 .
Given any of the converging subsequences (x (Nm) ) m>0 , denote by d its limit. The elements of the sequence (x (Nm) ) m>0 may be interpreted as a sequence of initial values for the mean-field limit described in Theorem 1. Consider the sequence of solutions of (39) along this sequence of initial-values. Observe that, as N m → ∞, such a solution has to converge to the solution x(t) of the differential equation (39) with initial condition d. However, d is a stationary distribution for the system with N → ∞ nodes, and therefore d must be a stationary point of the initial value problem (39). Due to Theorem 2, the only possibility is d = x * . We conclude that the sequence (x (N ) ) N ∈N converges to x * .
The above argument is well-known in the literature, and has been detailed in [4]. It has been applied in many contexts, see for instance [23,31,38]. Observe that, apart from the tightness of the sequence (x (N ) ) N ∈N , we relied only on the fact that there exists a unique equilibrium point for (6) (stability of x * is not required). In [28, Section 7.2] the authors provide an example where it is shown that if multiple limiting points exist, then the interchange of limits may lead to drastically different results.

Performance measures
In this section, we leverage the asymptotic results of Section 5 for the population process to obtain various performance measures of interest. In particular, we derive approximations for the stationary distribution and expectation of the queue length of a node, the sojourn time and waiting time of a packet in the system. Denote by Q

Convergence in distribution
Consider a tagged class-c node. In order to obtain an approximation for the stationary performance measures, we rely on the symmetry of the class-c nodes and we leverage the results established in the previous section. In particular, for N sufficiently large, we use the approximation and similar approximations for the waiting time and the sojourn time, and derive exact expressions for the limit on the right-hand side as stated in Theorem 4. In order to prove Theorem 4, we begin with the queue length distribution. Theorem 3 implies that lim We now continue with the waiting-time distribution. Denote by φ X (·) and F X (·) the Laplace-Stieltjes transform and the probability generating function of a random variable X. Define the random variableW see [5]. Moreover, we already proved that Hence, it only remains to be shown that the sequence of random variables {W (N ) c } N converges. We state the following proposition, which is a consequence of Little's law. Then, it follows from Prohorov's theorem, that the sequence {W (N ) c } N is relatively compact. Let W c be a limit random variable for an arbitrary convergent subsequence. It follows from (44) that However it is known that the Laplace-Stieltjes transform of a non-negative random variable is analytic in the positive open half plane (excluding the imaginary axis) [45,Chapter 8]. Moreover, from (45), it coincides with the Laplace-Stieltjes transform of a exponential random variable on a line segment. It is known that two analytic functions agreeing on a line segment must be the same on the entire domain, that is the non-negative half plane, see [1], so Therefore, we have shown that We now conclude with the sojourn time distribution. Define the random variableS c ) ∼ Exp(µ c ) denotes the transmission time of a class-c packet. Therefore Combining (43) and Theorem 4 we obtain an approximation for the stationary distribution of the queue length at a tagged class-c node given by i.e., a geometric distribution with parameter ξ c , for N sufficiently large, and similar approximations are obtained for the waiting and sojourn time distributions.

Convergence in expectation
In the previous subsection, we proved Theorem 4 obtaining the limiting distributions of the stationary queue length, waiting-time, and sojourn time distributions of the various classes. Such convergence results do not allow any immediate conclusions for the corresponding expectations as stated in Theorem 5. Here we show the convergence of the stationary expected waiting time from which the convergence of the stationary expected sojourn time and queue length will follow easily. Thanks to Theorem 4, we know that the rescaled stationary waiting time of a packet at a class-c node converges in distribution toW c ∼ Exp 1−ξc ξc . It is well-known [8] that and we want to show that In order to prove that, we will exploit again the equivalence with a 1-limited polling system with random routing described in Section 5.2. Specifically, the following formula is the pseudo-conservation law proved in [12, Eqn. (5.34)] specialized to a CSMA setting, i.e., p c,k = ν (N ) c /ν and σ c,k = ν, with ν = C c=1 ν c : where ρ = N c=1 ρ c . The above pseudo-conservation law is valid as long as the stationary expectations of the waiting times are finite, which is granted by [17,Theorem 5.5] as long as the stability condition (8) is satisfied.
We will show that there exist a sequence of constants K c,N > 0 for c ∈ C, N > 0 converging to K c > 0, such that Since K c , K c,N > 0 for every c ∈ C, N > 0, (46) and (48) imply It remains to show that (48) holds. As the pseudo-conservation law (47) suggests, we consider and consequently K c = pc µc 1 − λc where the second equality is due to (47). On the other hand, hence (48) holds, and we just showed The result for the sojourn time follows immediately since The argument we just presented implies that the sequence of random variables { λc Nc W (N ) c } N is uniformly integrable. Via the conventional Little's law we can then deduce that, for every c ∈ C, the sequence {Q (N ) c } N is itself uniformly integrable. Hence, it holds that

Model extensions
The theoretical machinery for the derivation of mean-field limits developed in Section 4 is highly flexible and may be applied to a wide range of models. In order to illustrate that, we outline in this section some possible extensions of the model described in Section 2, and state the associated mean-field limits. The proof of these limit results resembles that for the original model, and the arguments are not repeated. In particular, the time-scale separation between the activity process and the population process plays a key role once again. However, the expression for the stationary distribution of the activity process with a fixed population is model-dependent and this is reflected in the corresponding limiting initial-value problems which differ only in the expression for π x (Ω −c ). We recall that π x (Ω −c ) may obtained from the product-form solution for the stationary distribution of the activity process of saturated CSMA models.

Models with multiple transmission rates
In this model, the transmission of the packets may occur via different modes, which are associated with different distributions of the transmission times. Denote the set of modes which a class-c node may use by M c = {m c 1 , . . . , m c nc }. In particular, mode m c l is picked with probabilityp c l and the transmission time with that mode is phase-type distributed with mean (µ l c ) −1 > 0. Observe that a weighted sum of phase-type distributions is itself phase-type. This extension influences our original analysis only in Section 4.1.4, where we characterize the limit of α (N ) ω , yielding a more complicated version of the system of equations (34). However, it is shown in [40] that such a system can be uniquely solved and its solution only depends on the mean of the phase-type distributions involved, due to an insensitivity property for the stationary distribution of the activity process in saturated networks.
Hence, the mean-field limit of this generalized model turns out to be very similar to Theorem 1, and we only need to characterize the counterpart of π x (Ω −c ). Denote by U c the expected value of the transmission time of a class-c node, i.e., where the function H mtr (·) is defined by Note that in the original case we have only one rate per class, yielding U c = 1/µ c for every c ∈ C, in which case Theorem 6 reduces to Theorem 1. Observe that all the results in Sections 4.2 and 5 may be extended to the model with multiple transmission rates, by replacing µ c with 1/U c . Moreover, since phase-type distributions approximate any nonnegative distribution, we deduce that the mean-field limit is insensitive with respect to the transmission time distribution.

Models with queue-based back-off rates
In the basic model, nodes with nonempty buffers back-off at a rate that depends only on their class. However, in [24,35], the authors showed that CSMA networks with suitable queue-based back-off rates are maximally stable without any explicit knowledge of the arrival rates in arbitrary interference graphs.
Consider the model described in Section 2, where the back-off rate of a node is a function of the number of packets in its buffer, i.e., the back-off rate of a class-c node having n packets in the buffer is ν (N ) c (n) = ν c (n)/N c , see [11]. Note that when ν c (n) = ν c for every n we recover the basic model. The system dynamics are very similar to those of the original model. Specifically, we only have to replace the back-off transitions, as follows: A back-off is completed by a class-c node having n > 0 packets in its buffer: this happens only if the class activity state allows class-c nodes to back off, i.e., Y (N ) ∈ Ω −c , and at rate ν In accordance with the basic framework, the following theorem provides a description for the mean-field limit of the CSMA model with queue-based back-off rates.
where the function H(·) is defined by x c,n e c,n+1 x c,n = ν c (1 − x c,0 ) and the original mean-field limit in Theorem 1 may be interpreted as a special case of Theorem 7.
In order to understand the formula for π qb x (ω), consider a fictitious saturated CSMA model, with C nodes and the original class-interference graph and a time-invariant population x ∈ E. The time required to transmit a packet originating from node c is exponentially distributed with parameter µ c , while the length of a back-off period of node c is captured by the population-dependent random and therefore V c (x) ∼ Exp( ∞ n=1 ν c (n)x c,n ). In [40] the authors proved the insensitivity of the stationary distribution of the activity process for saturated CSMA networks with respect to the back-off distribution. Hence, we know that the expression for the stationary probability of the activity process to be in state ω corresponds to π qb x (ω). Intuitively, the cumulative back-off time needed by class-c nodes is given by the minimum of the back-off time required by each node in the class, whose distribution depends on the number of packets in their buffer.
whereν c = ν c η(ρ) c , then there exists a unique equilibrium point x * = (x * c,n ) c∈C,n∈N ∈ E for the initial-value problem (49). In particular We now briefly sketch the main argument needed to prove Theorem 8. The relation H qb (x * ) = 0 must be satisfied component-wise, i.e., for every c ∈ C and n > 0. Hence we only need to show that necessarily ξ c (n) =ν c /ν c (n). Given x ∈ S, define the vector of effective class-back-off ratesν(x) ∈ R C + as follows Observe that , and x * ∈ E is an equilibrium point only ifν(x * ) =ν whereν c = η(ρ) c with σ c = 1/µ c for every c ∈ C. Hence, the uniqueness ofν is ensured by the global invertibility of the map θ(·) over Γ. Therefore, if x * is such that H qb (x * ) = 0, then and consequently ξ c (n) =ν c ν c (n) .
As a simple example, consider the case with a single class and ν(n) = nν. Observe thatν must solve λ µ =ν /µ 1 +ν/µ We now need to verify that L < ∞, Hence, we conclude that

Models with finite capacity
A practical constraint in the implementation of CSMA is that the number of packets stored in the buffer of a node must not exceed a certain capacity. In this section, we assume that the buffer capacity of every node in class c is limited to be no more than a certain value K c ∈ N. In particular, this feature yields immediately X (N ) c,n = 0, ∀ n > K c , c ∈ C.
Then the sequence of Markov processes has a continuous limit x(t) ∈ C E [0, ∞) which is determined by the unique solution of the initial-value problem where the function H f b (·) is defined by x c,n e c,n+1 x c,n e c,n−1 c,n .
Observe that when K c = ∞ for every c ∈ C we recover Theorem 1 as a special case of Theorem 9.
Focus now on the analysis of the equilibrium points for (51). For any x ∈ E 1 , denote .
Consider an equilibrium point x * ∈ E 1 for (51), i.e., H f b (x * ) = 0, then x * c,n = (ξ x * c ) n x * c,0 , ∀ n = 0, . . . , K c , c ∈ C, and x * c,n for every n > K c . Hence, x * ∈ E 1 is an equilibrium point, if and only if We now identify sufficient and necessary conditions on ξ ∈ R C + such that the population vector x ξ ∈ E 1 is an equilibrium point for (51), where Assume the population vector x ξ ∈ E 1 to remain constant over time, i.e., x ξ is an equilibrium point for (51). Then, the fraction of time allocated to the activity of class-c nodes is given bỹ On the other hand, the effective stationary load of class-c nodes is given bỹ i.e., a packet arriving at a class-c node with a full buffer is lost and does not enter the system so that the effective arrival rate at class-c nodes is λ c (1 − x ξ c,Kc ). We conclude that x ξ ∈ E 1 is an equilibrium point for (51) if and only ifθ (ξ) =ρ(ξ). (52) Consider the single-class case, i.e., C = 1, and drop the c subscript. Observe that so thatθ(0) = 0, lim ξ→∞θ (ξ) = (ν/µ)/(1 + ν/µ) > 0, and ∂ ∂ξθ (ξ) > 0. On the other hand, andρ(0) = λ/µ > 0, lim ξ→∞ρ (ξ) = 0, and ∂ ∂ξρ (ξ) < 0. Hence, the uniqueness of ξ satisfying (52) is ensured for every λ, µ, ν, K, and in turn yields the uniqueness of the equilibrium point for the initial-value problem (51).

Numerical experiments
In this section we perform numerical experiments in order to complement the analytical results obtained in the previous sections. In particular, we provide evidence that the results proved in Sections 5-7 for complete interference graphs in fact extend to more general graphs.
Throughout this section we focus on the square network displayed in Figure 1 and considered at the end of Section 4 with the following parameters: Observe that the above system of equations cannot be solved explicitly in general. Hence, we use the numerical algorithm described in [15], which is an adaptation of the algorithm presented in [41], and obtain ξ ≈ (0.4302, 0.2635, 0.6537, 0.3442).
Theorem 2 thus yields that Convergence of population process. In Theorem 1 we proved the convergence of the properly scaled population process {X (N ) (N t)} to the solution x(t) of (6) for arbitrary interference graphs. In Proposition 6 we showed global stability of x * , the unique fixed point of (6), for complete interference graphs. In Figure 3 we display a numerical solution of (6) and observe that x(t) → x * . This supports the hypothesis that the solution x(t) of (6) converges to x * even when the interference graph is not complete.
Convergence of the stationary distribution of the population and the queue length process. In Theorem 3 we proved that lim N →∞ x (N ) (∞) = x * , where x (N ) (∞) is a random variable with the stationary distribution of X (N ) (t) as t → ∞, for complete interference graphs. Also, in Theorem 5, we showed that lim is a random variable with the stationary distribution of the queue-length Q (N ) c at a class-c node as t → ∞. In order to examine these two properties, we conducted 10000 simulation runs and collected 10000 observations (x (N ),h c ) h=1,...,10000 of the population process for large t and values N = 2 l for l = 2, . . . , 9. We use the average observed distribution of the population process for large t in a network with N nodes, i.e., In Figure 4 we display in the bar plot the euclidean distance d c (x * c ,x (N ) c ) between the asymptotic approximation x * c and the empirical estimatex (N ) c for values N = 2 l with l = 2, . . . , 9, where The different colors in the plot correspond to the various classes. Observe that the distance is small even for moderate values of N . In the plot on the right we display the empirical estimate CDF (N ) 1 (n) along with the asymptotic value (ξ * 1 ) n for n = 0, . . . , 12 and N = 2 l with l = 2, . . . , 9. Observe that as N increases, the accuracy of the asymptotic approximation (ξ * 1 ) n for the tail probability P{Q (N ) 1 (∞) ≥ n} (in red) improves, and the difference becomes almost indiscernible. In order to save space, we do not display the tail distribution convergence for the classes c = 2, 3, 4 which show similar behavior. These observations suggest that the convergence of the stationary distribution of the population process and the queue length holds in non-complete interference graphs as well, and that the asymptotic values provide accurate approximations even for moderate values of N .
Convergence of mean waiting time. In Theorem 5 we proved that a properly scaled version of the mean waiting time λ c E[W (N ) c ]/N c converges to ξ c /(1 − ξ c ) as N → ∞ for complete interference graphs. In order to obtain an empirical estimate for the mean waiting time from the simulation runs, we invoke Little's law In Figure 5 we display the empirical estimates E[W (N ) c ]/N c , c = 1, 2, 3, 4, along with the asymptotic approximations ξ * c /λ c (1 − ξ * c ) for values N = 2 l with l = 2, . . . , 9. The asymptotic approximations turn out to be extremely accurate even for N = 2 6 = 64, i.e., N c = 16. The results in Figure 5 support the conjecture that the mean waiting-time approximations in Theorem 5 extend to noncomplete interference graphs, and are accurate even for moderate values of N .

Conclusion
We analyzed distributed random-access networks in a many-sources regime with buffer-dynamics. The difficulty in the analysis is the complex interaction between the medium activity state and the buffer content process. However, in the mean-field regime, these processes are shown to decouple and a tractable initial-value problem describing the buffer content process is obtained. As illustrated through several examples, the methodological framework that we developed is fairly generic, and expected to be applicable to a broader range of problems where a similar decoupling occurs.
In the second part of the paper we focus on the analysis of the key performance measures for the case of a complete interference graph. The stationary buffer content and the packet waiting and sojourn time are shown to converge in distribution and expectation to random variables whose parameters depend on the unique equilibrium point of the initial-value problem and are explicitly provided. These results are obtained via an interchange of limits approach which ensures that the resulting approximations are asymptotically exact. Numerical experiments suggest that these results extend to general interference graphs, but a rigorous proof remains as a challenging problem for further research.
It would be of interest to develop a mean-field analysis for random-access networks in different regimes. In this work we assumed both arrival and back-off rates to scale as the number of nodes increases. Such scaling allows analytical tractability thanks to the decoupling between the medium activity state and the buffer content process. However, the back-off rate scaling entails a degradation in the delay performance which scales with the number of nodes as well. The study of large networks where only the arrival rate scales would be important since we expect this degradation would not happen. However, different analytical tools are required and we will investigate this problem in future work.

A Asymptotically Lipschitz in probability
In the course of obtaining the mean-field limit, we will use the property of a sequence in D E [0, ∞) of being asymptotically Lipschitz in probability and its relation with the relative compactness of a sequence of processes taking values in a metric space E with metric ρ and sample paths in D. Various versions of this property have been described in the literature, e.g., [13,37], and we will use the following definition.
Definition 3. A sequence of random paths X (N ) ∈ D E [0, ∞) is asymptotically Lipschitz in probability for a given T > 0 if there exists K ≥ 0 such that for every t 1 , t 2 ∈ [0, T ] and for every > 0 the following is satisfied: We now reproduce [20,Theorem 3.7.2], and also refer to [8,Theorem 16.4], and the subsequent discussion. (a) For every η ∈ R >0 and b ∈ Q ≥0 , there exists a compact set This leaves us with condition (b) and we therefore proceed to show the next proposition.
Consider condition (b) in Theorem 10 and assume η and T given. Note that it suffices to verify this condition for N ≥ N 0 , as [20, Lemma 6.2 (a)] shows that w is right continuous and therefore we may choose δ < δ N0 to ensure (b) holds for the first n ≤ N 0 cases.
Fix < η/2. Now choose N 0 ( , η) such that Observe that N 0 ( , η) exists due to Definition 3. Now choose a uniform and sufficiently fine partition so that the spacing is δ 1 where δ 1 K < η/2 and denote this partition as Π δ1 . Observe that for N > N 0 ( , η), for every i = 1, . . . , n Π δ 1 . We now choose δ = δ 1 /2 and observe that by definition of w (X (N ) , δ, T ), Indeed, w (X (N ) , δ, T ) is given by the infimum over the suitable partitions and therefore is certainly lower than the outcome obtained by considering the suitable partition Π δ1 , hence For the remaining n < N 0 make δ smaller as described above to obtain the required result as w is non-decreasing in δ.
B Extended proofs B.1 Proof of Lemma 2 As previously stated, we consider the single-class case only. The following arguments are along the same lines as given in [30]. First observe that the unit Poisson processes for the arrivals, the back-offs, and the transmissions are all mutually independent and by assumption are defined on a common probability space which satisfies the usual conditions. Moreover, the arrival rates, back-off rates and the transmission rate are all bounded. In particular, the natural filtrations are defined with the σ-algebra generated by N A,n (u) for all 0 ≤ u ≤ t. We make corresponding definitions for the back-offs H B n (t), n ∈ N, and for the transmission process denoted by H T (t). Given a multi-parameter time s = (s 0 , s 1 , s 2 , . . .) ∈ [0, ∞) ∞ , define H(s) to be the σ-algebra where N C Null is the collection of all null sets in F C . With this construction we thus obtain the multi-parameter filtration H .
= H(s) : s ∈ R ∞ + , see [20, page 85], where the metric space for the index set is R ∞ + with the usual metric. (We may suppose that this filtration is right continuous as otherwise we may work with its right continuous augmentation, see [34,Lemma 67.4] for example.) Also the process N T ,Ñ A,0 ,Ñ A,1 ,Ñ B,1 , . . . is a right continuous multi-parameter martingale adapted to the filtration H.
Since X, Y are the unique solutions to (12)- (16), it follows that we can apply [20, Theorem 2.2, page 313], to show that τ is a multi-parameter H-stopping time, where τ A n (t) = H(τ (t)) where τ is the above H stopping time. Since τ (t) is a right continuous function of t and H(s) is a right continuous filtration, it follows that if η is a F(t) stopping time, then so is τ (η).
It follows from the Optional Sampling Theorem for martingales on a metric lattice [20, Theorem 8.7, page 87] that M A is a right continuous F multi-parameter martingale. This martingale is locally square integrable, because the rates are uniformly bounded over time. In particular if T > 0 the version of these martingales stopped at T is a L 2 F-martingale. The same arguments show that the vector back-off process, is also a locally square integrable F-martingale as of course is the processÑ T N µ t 0 ds . However Y ∈ mF(t) which implies that M B is also a locally square integrable F-martingale as the integrand is previsible and similarly for M T . Finally Ξ is adapted, being a countable sum of F martingales, and it is straightforward to show that it is bounded in L 2 and thus a L 2 martingale.  [44]) that we may extend this result to deduce that L[0, ∞) is compact within C R C + ,↑ [0, ∞). Therefore a sample path of α (N ) lies in a compact set within C R C + ,↑ [0, ∞). This implies that the sequence of measures on α (N ) is tight. We turn to the sequence of population processes X (N ) . First observe that the space E is compact as a consequence of Tychonoff's Theorem as it is a finite product of compact spaces, see Lemma 1 and the preceding discussion. In Appendix A the concept of asymptotically Lipschitz in probability is introduced and its relation with relatively compactness is stated in Proposition 11. Lemma 4, which is stated and proved below implies that the sequence of population processes is asymptotically Lipschitz in probability.
Lemma 4. The sequence of processes (X (N ) ) N is asymptotically Lipschitz in probability.
However since the space (E, ρ) is complete and separable, it follows from the converse half of Prohorov's Theorem that the sequence X (N ) is tight, [7, Theorem 6.2, page 37]. It follows, as already stated, that the sequence of measures (X (N ) , α (N ) ) is therefore tight. Hence, the sequence (X (N ) , α (N ) ) is relatively compact, again as a consequence of Prohorov's Theorem, [7, Theorem 6.2, page 37].
Proof. (of Lemma 4) Recall that E = χ C and C < ∞, see (1). Focus on (25), and observe that its (c, n) component can be written as where the centered Poisson processes are martingales bounded in L 2 . We recall that Doob's submartingale inequality, see [ and it follows thatÑ ·,c,n are non-negative submartingales. Also, yielding that and a similar argument applies for N B,c,n . From (54) we obtain that From (56), we obtain for every t, u ∈ [0, T ]. Observe that, thanks to (57), it holds that for every t, u ∈ [0, T ].

B.3 Proof of Proposition 3
As shown in Lemma 2, the process Ξ is a multi-parameter martingale adapted to F. Since C is finite, there exists K η such that c∈C n>Kη |Ξ (N ) Applying Doob's submartingale inequality with p = 1 (see [33]), we obtain P sup x c,n (s)α ω (s)ds.

B.5 Proof of Lemma 3
Given > 0, we may find a partition of [0, T ], t 0 = 0 < t 1 < · · · < t K = T such that, 0 ≤ u(t k ) − u(t k−1 ) ≤ /2, k = 1, · · · , K, since u is Lipschitz continuous by assumption. Next given any t ∈ (0, T ) the following inequality holds, for some k < K, since if on the one hand, u (N ) (t) ≤ u(t), then the first term on the right is a bound and on the other hand if the inequality is reversed then the second term is a bound. Using the triangle inequality on each respective term, we obtain Since the partition is fixed, we can find N so that the right hand side is smaller than using pointwise convergence and the properties of the partition. Clearly, this N can be taken to cover the case t = 0, T as well and the result is proved.

B.6 Proof of Proposition 5
For notational convenience, the proof is only presented for the single-class case, but the arguments easily extend to the multi-class case.
Consider solutions x L (t) and x U (t) of the initial-value problem (39). These solutions are continuous and differentiable in t, therefore it is sufficient to show that, when the solutions are on the boundary of the stochastic dominance region, the solutions are forced to remain inside such region. Specifically, consider a time t ≥ t 0 such that the stochastic dominance relation x L (t) ≤ s x U (t) holds, i.e., the solutions are inside the stochastic dominance region, and there existsk ∈ N 0 such that i.e., the solutions are on the boundary of the region, then we need to show that i.e., the solutions are forced to remain inside the stochastic dominance region. From (39) we obtain and in fact π x L 0 (t) (0) ≥ π x U 0 (t) (0) since x L 0 (t)) ≥ x U 0 (t). Observe that the stochastic dominance property and (59) yield and therefore Hence, the nonnegativity of (60) follows.

B.7 Proof of Proposition 6
For notational convenience, the proof is only presented for the single-class case, but the arguments easily extend to the multi-class case.
Assume now that x(0) ∈ E 1 <∞ , we will prove that the solution x(t) of the initial-value problem (39) is then such that x n (t) → x * n for every n ∈ N 0 . The proof relies on induction on n, and the base case is proved in the following lemma.
The inductive step is proved in the following lemma.
Lemma 6. Consider a solution x(t) = x n (t) n≥0 of (39) with x(0) ∈ E 1 <∞ and such that x m (t) → x * m for every m ≤ n − 1. It holds that x n (t) → x * n .
The fact that x m (t) → x * m for every m ≤ n − 1 implies that there exists T such that for every t ≥ T , we have that |π x0(t) (0) − π x * 0 (0)| ≤ and |x m (t) − x * m | ≤ for every m = 0, . . . , n − 1. Consider now t ≥ T , then there exist functions (g m ( )) m=2,...,n such that Recall that the fundamental theorem of calculus yields (0) = x * n due to the definition of x * in Corollary 2. Since these inequalities hold for every h > 0, thanks to the Lipschitz continuity of x n (t), we conclude that x n (t) → x * n .
In order to complete the proof of Proposition 6, we have to show that given any solution x(t) of (6) with initial-value x(0) ∈ E 1 \ E 1 <∞ it holds that x(t) x * . Note that a necessary condition for x(t) → x * is the convergence of the masses, i.e., m(x(t)) → m(x * ). Note that m(x(0)) = ∞ and m(x * ) < ∞. Hence, since the function m(x(·)) is Lipschitz continuous, the convergence of the masses cannot happen in finite time, and in turn x(0) / ∈ B(x * ).

B.8 Proof of Proposition 8
Here we relate the notation in our framework with the notation in [18], so as to make the proof of Proposition 8 as transparent as possible. For ease of notation we omit the superscript N .
The N queues in [18] are identified in our framework by (c, k) for c ∈ C and k = 1, . . . , N c . The load of queue (c, k) is ρ c,k = λ c,k µ c,k , where λ c,k = λ c /N c and µ c,k = µ c . In our framework when the server visits a queue with m packets in the buffer S(m) packets are served, where S(m) = 0, if m = 0, 1, if m > 0.
HenceS = lim m→∞ S(m) = 1. Denote by r (c,k)(d,h) the probability that the server visits queue (d, h) after having visited queue (c, k). In our case, the queue visited by a server is not dependent on the past visits, so that where ν = c∈C ν c . Also, the switchover times {δ (c,k)(d,h) (n)} n≥0 , i.e., the time the server needs for switching from queue (c, k) to (d, h) at the n-th occurrence of such a switch, are independent of the queue just left. More explicitly, the random variables δ (c,k)(d,h) (n) are i.i.d. copies of the random variable δ (c,k)(d,h) ∼ Exp ν .
In particular, δ (c,k)(d,h) is distributed as the minimum of N exponential random variables each corresponding to the back-off time of a different node whose parameter depends on the class in which it belongs. The routing probabilities r (c,k)(d,h) define a Markov chain on the state space {(c, k) : c ∈ C, k = 1, . . . , N c }. Let p (c,k) be its unique invariant distribution. In our framework, due to the independence of r (c,k)(d,h) from (c, k), it holds that Lastly, define δ * =

B.9 Alternative proof for Proposition 8 -Single-class case
Fix the number of nodes N . For every node k, define a Poisson process of rate ν/N , which represents the back-off opportunities. By opportunity, we mean that the corresponding event (back-off in this case) happens if the condition of the system allows for it. Denote by b (k,N ) n a realization of the backoff opportunity process at node k, i.e., at time b n ) = 1. The arrivals at node k are governed by a Poisson process of rate λ/N , and in this case we do not need to consider opportunities. Indeed at every realization of the process a new packet arrives, no matter the state of the system. These processes are sufficient for describing the stationary behavior in the single-class case with N nodes, a special case of the model described in Section 2.
Consider now the following random sets of times Observe that different elements of S (N ) may identify the same element in T (N ) . Basically, the backoff opportunities in T (N ) are the only opportunities taken by nodes which are not blocked by the activity process. Observe that the absence of packets in their buffer might still prevent them from initiating a transmission. Drop the superscript (N ), and denote the sequence of elements of T (N ) = T by {t n } n∈N . We want to analyze the embedded queue length Markov chain (Q n ) n∈N on the state space N N where Q n = Q(t n ), and Q k (t) is the length of the k-th queue at time t.
Define the following random variables • T n = t n+1 −t n represents the length of the n-th time slot. Observe that (T n ) n∈N is a sequence of i.i.d. random variables whose distribution is given by the sum of two exponentially distributed random variables of rate µ and ν. Observe that E[T n ] =T = 1 µ + 1 ν .
• A n = (A n k ) k=1,...,N counts the number of packets arriving at the various queues during the n-th time slot. Observe that A n k | Tn=s ∼ Poisson( λ N s) for every k = 1, . . . , N . • D n = (D n k ) k=1,...,N counts the number of packets leaving the various queues during the n-th time slot. Note that D n k ∈ {0, 1} and k D n k ≤ 1.
Define L(Q n ) = N k=1 (Q n k ) 2 . We will show that L(·) is a Lyapunov function for the embedded queue length Markov chain (Q n ) n∈N . In order to prove that, it suffices to show that E[L(Q n+1 )|Q n ]− L(Q n ) < 0 unless Q n ∈ K, where K ∈ R N + is a compact set.
it holds that .
Consider S 1 , and observe its terms individually. As for the first term, it holds that Focus now on E[D n k |Q n ] for a generic node k ∈ {1, . . . , N }. Denote by C n j the random variable indicating the node associated with the j-th back-off opportunity during the n-th time slot, i.e., in [t n , t n+1 ). Denote by J n the random variable counting the number of back-off opportunities during the n-th time slot. Observe that J n ≥ 1 due to the time-slot definition. Finally, define the random variables Thus Consider now S 2 , and observe that its value is bounded by a constant. Indeed, Combining (63) and (64) we obtain Since λ/µ + λ/ν < 1, there exists a compact set of the form of K = {Q : 1 N N i=1 Q i < c} such that if Q / ∈ K then E[L(Q n+1 )|Q n ] − L(Q n ) < 0.