# Bias and Refinement of Multiscale Mean Field Models

Sebastian Allmeier, Nicolas Gast

► To cite this version:

Sebastian Allmeier, Nicolas Gast. Bias and Refinement of Multiscale Mean Field Models. Proceedings of the ACM on Measurement and Analysis of Computing Systems, 2023, 7 (1), pp.1-29. 10.1145/3579336. hal-03981002

HAL Id: hal-03981002

https://inria.hal.science/hal-03981002v1

Submitted on 8 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Bias and Refinement of Multiscale Mean Field Models

SEBASTIAN ALLMEIER and NICOLAS GAST, Inria, France and Univ. Grenoble Alpes, France

Mean field approximation is a powerful technique which has been used in many settings to study large-scale stochastic systems. In the case of two-timescale systems, the approximation is obtained by a combination of scaling arguments and the use of the averaging principle. This paper analyzes the approximation error of this ‘average’ mean field model for a two-timescale model (X, Y), where the slow component X describes a population of interacting particles which is fully coupled with a rapidly changing environment Y. The model is parametrized by a scaling factor N, e.g. the population size, which as N gets large decreases the jump size of the slow component in contrast to the unchanged dynamics of the fast component. We show that under relatively mild conditions, the ‘average’ mean field approximation has a bias of order $O(1/N)$ compared to $\mathbb{E}[X]$. This holds true under any continuous performance metric in the transient regime, as well as for the steady-state if the model is exponentially stable. To go one step further, we derive a bias correction term for the steady-state, from which we define a new approximation called the refined ‘average’ mean field approximation whose bias is of order $O(1/N^2)$. This refined ‘average’ mean field approximation allows computing an accurate approximation even for small scaling factors, i.e., $N \approx 10-50$. We illustrate the developed framework and accuracy results through an application to a random access CSMA model.

CCS Concepts: • Mathematics of computing → Markov processes; Ordinary differential equations.

Additional Key Words and Phrases: Mean Field; Refined Mean Field; Two Timescale System

## 1 INTRODUCTION

The mean field approximation finds widespread application when interested in analyzing the macroscopic behavior of large-scale stochastic systems composed of interacting particles. Its assets lie in a reduction of the model complexity, simplified analysis of the system due to absence of stochastic components, and reduction of computation time compared to a stochastic simulation. The mean field approximation can even yield closed form solutions for the steady-state, e.g., for the well known JSQ(d) model [25]. The mean field approximation is generally given by a set of ordinary differential equations which arise from the assumption that, for large system sizes, the evolution of the particles are stochastically independent of another. This idea works well if the number of particles is large and if the particles can be clustered into a few groups with statistically identical behavior [2, 14, 21-23]. The framework established by Kurtz [21] to derive (weak) convergence results for the stochastic system justifying the use of the mean field approximations finds sustained attention in the literature.

More recently, the authors of [12, 33] showed that for finite system sizes the bias of the mean field approximation is of order $1/N$ when compared to the mean behavior of the system. Here, N is the scaling parameter, which usually refers to the number of homogeneous particles in the system. Additional works such as [13, 15] introduced corrections, called refinement terms, which effectively increase the rate of accuracy of the approximation and therefore the rate of convergence. The most notable term is the first-order bias refinement, since it offers a convincing trade-off between a significant accuracy gain and additional computation cost.

While these classical mean field results hold for a broad class of models, most of the results can not be transferred to systems with more intricate dynamics such as the two-timescale case, studied for instance in [2, 4]. A two-timescale process consists of two coupled components, one evolving slowly compared to the other. The slowly evolving component is often represented by a system of interacting particles, where the state of each particle evolves as a function of the empirical
distribution of all particles but also as a function of the state of the fast component, e.g., a fast changing environment. These types of processes and their ‘averaged’ mean field adaptation have been of interest since the 1960s and became increasingly relevant in the study of modern and complex systems. We refer to [26] for an extensive literature discussion. An important area of application comes from the field of computer networks. Examples include loss networks [16], large-scale random access networks with interference graphs [8, 10] or storage networks [11]. Another recent field of literature from which we draw inspiration are chemical reaction networks. The works of Kang et al. [19, 20] and Ball et al. [1] establish central limit theorems for large multiscale models motivated by biological and chemical processes. Another application in the field of biology is given by [29] who use the mean field idea to study neural plasticity models.

*Contributions.* The aforementioned papers underline that the mean field idea can be adapted to two-timescale models and prove that the ‘average’ mean field approximation is asymptotically exact as $N$ goes to infinity. They do not, however, provide theoretical results which guarantee the accuracy or performance bounds of the approximation for finite systems sizes. This paper aims at filling this gap. We derive accuracy bounds for the ‘average’ mean field approximation in the transient regime and steady-state which show that its bias is of order $O(1/N)$, $N$ being the scaling parameter of the stochastic system. We further derive bias correction terms for the steady-state, from which we define a new approximation called the refined ‘average’ mean field approximation, whose bias is of order $O(1/N^2)$. To prove these accuracy bounds, we develop a framework for two-timescale stochastic models whose slow component is comparable to the concept of density dependent population processes as introduced by Kurtz [22]. Based on this representation, we utilize a combination of generator comparison techniques, Poisson equations as well as derivative bounds on the solution of the Poisson equation. This allows us to bound the bias with respect to the scaling parameter $N$. To take it a step further, we then prove the existence of correction terms which approximate the bias of the ‘average’ mean field and allow defining the refined ‘average’ mean field. To support the practical application of the refinements, we provide an algorithmic way to compute the correction terms. This includes methods to numerically solve the Poisson equation and obtain its derivatives. We illustrate the computation of the refinement terms and confirm the accuracy of the obtained bounds by considering a random access CSMA model. Using the example we show that even for relatively small $N \approx 10$ the refined ‘average’ mean field approximation almost exactly indicates the steady-state of the model.

*Methodological Advances & Technical Challenges.* To obtain our results, we build on the recent line of work on Stein’s method [30]. This method allows calculating the distance between two random variables by looking at the distance between the generators of two related systems. Recently, the method reemerged in publications within the stochastic network community, in particular the works of Braverman et al. [5–7]. In this paper, we use the Poisson equation idea in two ways. First and foremost, we use a Poisson equation, called the ‘fast’ Poisson equation, to bound the distance between a function – in this case the drift – and its average version given by the steady-state distribution of the fast process as introduced in Section 4.2. This step is integral to our paper and constitutes the building block to obtain the error bounds and closed form expressions as it allows to analyze the distance between the coupled stochastic process and its decoupled counterpart. The analysis however bears many technical intricacies one needs to overcome. This includes solving the Poisson equation, stating its derivative bounds and deducing computable expressions used to calculate the bias term. Second, for our steady-state results, we make use of another Poisson equation to compare the stochastic system to the equilibrium point of the approximation. The latter is related to the methods used in [15, 33] but significantly extends the ideas as the derivation of the bias exhibits novel refinement terms which originate from the coupling and correct the error of
the ‘averaging’ method used for the mean field. The closely linked technical challenges include obtaining numerically feasible formulas for the ‘new’ refinement terms. Here, we utilize the derived solution of the ‘fast’ Poisson equation to further specify the fluctuations of the stochastic system around the equilibrium point. Carefully analyzing the combination of the two Poisson equations enables us to obtain the new refinement terms. The closed-form bias terms that we obtained from the steady-state analysis are significantly more complicated than the ones of [15] as they further correct the error made by the averaging method.

***Applicability & Numerical Difficulties.*** In our paper, we make use of a CSMA model to demonstrate the applicability of our results. There are several other examples captured by our framework, such as the Michaelis-Menten enzyme model of [20] or the storage network investigated in [11]. For the latter, the authors observe that when looking at the right timescale the loss of the network can be characterized by a local equilibrium which is obtained using the averaging method. The biochemical Michaelis-Menten enzyme model describes the dynamics between three time-varying species, the enzyme, a substrate and a product. Following the description of the model as in [20] and by using the right scaling arguments, the model exhibits two timescales, the fast reacting and state changing enzymes and the slower changing concentrations of the substrate and product. Our framework can be used on both examples to guarantee accuracy results and can be used to compute refined approximations.

To compute the ‘average’ mean field and refinement terms, one has to overcome numerical difficulties which arise from the averaging method. First, in order to compute the ‘average’ drift, one needs to compute the steady-state probabilities of the fast system, which might not be available in the closed form as for the CSMA model. For this case we provide computational notes in the appendix which aim to facilitate the computation. Another problem for the refinement term is the need for the first and second derivatives of the ‘average’ drift. For the CSMA model we used symbolic representation of the transition rates from which we define the drift and its average version. This method allows to numerically compute the derivatives using ‘sympy’, a Python library for symbolic computation. This method is relatively easy to implement but has a very large computation time. This computation time could be reduced by implementing a faster computation of the derivatives. In fact, to compute the derivative of the matrix $K^+$ which is closely linked to the solution of the ‘fast’ Poisson, we provide supplementary computational notes which describe how to obtain an efficient implementation.

***Roadmap.*** The paper is organized as follows. In Section 2, we formally introduce the two-timescale model, its corresponding ‘average’ mean field approximation and make regularity assumptions on the system. In Section 3 we state the main results of this paper. Section 3.1 states the results for the transient regime, Section 3.2 for the steady-state and Section 3.3 justifies the existence of bias correction terms. Section 4 holds the proofs of the aforementioned results. In Section 5 we apply our theoretical results to the unsaturated random-access network model of [9, 10]. Some technical lemmas and definitions are postponed to the appendix.

***Reproducibility.*** The code to reproduce the paper along with all figures and the implementation of the unsaturated random-access network model is available at https://gitlab.inria.fr/sallmeie/bias-and-refinement-of-multiscale-mean-field-models.

## 2 STOCHASTIC SYSTEM AND MEAN FIELD APPROXIMATION

We consider a two-timescale, coupled, continuous time Markov chain $(X_s^{(N)}, Y_s^{(N)})_{s \ge 0}$ parametrized by a scaling factor $N$, for which we study the behavior as $N$ tends to infinity. As we will see in
the examples, N typically represents the number of objects that interact together. This section introduces the precise model and fixes notations.

## 2.1 Model

For a fixed scaling parameter N, the stochastic process $(X_s^{(N)}, Y_s^{(N)})_{s \ge 0}$ is a continuous time Markov chain that evolves in a state-space $\mathcal{X}^{(N)} \times \mathcal{Y}$. The set $\mathcal{Y}$ is finite and does not depend on N. We further require that for all N, the sets $\mathcal{X}^{(N)}$ are subsets of a convex and compact set $\mathcal{X} \subset \mathbb{R}^{d_x}$. In what follows, unless it is ambiguous in the context, we drop the dependence on N to lighten the notations.

This model has two-timescales in the sense that the size order of jumps of $Y_s$ is N-times larger than the ones of $X_s$. More precisely, we assume that there exists a finite number of transitions $(\ell, y') \in \mathcal{T}$ with their corresponding transition rate functions $\alpha_{\ell, y'} \ge 0$, both being independent of N, such that for all possible states $(X_s, Y_s) \in \mathcal{X} \times \mathcal{Y}$:

$$ (X_s, Y_s) \text{ jumps to } (X_s + \ell/N, y') \quad \text{at rate } N\alpha_{\ell, y'}(X_s, Y_s). \quad (1) $$

The above defines continuous time Markov chains with discrete state-space whose realizations are Càdlàg, i.e., right continuous with a left limit for every time t. Note that in Equation (1), we assume that the transition rates $\alpha$ are defined for all $x \in \mathcal{X}$ (and not just for $x \in \mathcal{X}^{(N)}$).

The notion of slow-fast system comes from the fact that the jumps of the slow component $\mathcal{X}$ are $O(N)$ times smaller than the jumps of the fast component $\mathcal{Y}$ while transition rates are of the same scale. As we will see later, the different timescales imply that for large N, the slow component $X_s$ will 'see' the fast component $Y_s$ as if it is stationary with distribution $(\pi_y(X_s))_{y \in \mathcal{Y}}$ which we formally define in Section 2.2.

## 2.2 Drift, Average Drift and Mean Field Approximation

The jumps of the stochastic system (1) can affect the fast and/or the slow component. In what follows, we construct an approximation that consists (i) in considering that the slow component is not stochastic but evolves deterministically according to its drift (which is its average change), and (ii) in using a time-averaging method that shows the stochastic process Y as being in some stationary state given $x$. This leads us to two definitions:

(i) We call the drift of the slow system (or more concisely the drift) the average change of $X_s$. It is the sum over all possible transitions of the rate of transition multiplied by the changes that such transitions induce on $X_t$. By the form of the transitions in (1), if the process starts in $(X_t, Y_t) = (x, y)$, the drift is given by:

$$ F(x, y) := \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \ell \in \mathbb{R}^{d_x}. \quad (2) $$

This drift function depends on the state of the fast system $y$.

(ii) For the fast component, we define a transition kernel $K_{y,y'}(x)$ that is the rate at which the process Y jumps from $y$ to $y' \ne y$ (divided by N), with the usual convention that $K_{y,y}(x) = -\sum_{y' \ne y} K_{y,y'}(x)$:

$$ K_{y,y'}(x) = \sum_{\ell} \alpha_{\ell,y'}(x,y). \quad (3) $$

As $\mathcal{Y}$ is finite, for a fixed $x$, $K(x)$ is a matrix that corresponds to the kernel of a continuous time Markov chain. Our assumptions will imply that for all $x$, the process associated with $K(x)$ has a unique stationary distribution, which we denote by the vector $\pi(x) = (\pi_y(x))_{y \in \mathcal{Y}}$.
Based on the drift (2) of the stochastic system, we define its ‘average’ version $\bar{F}$ by averaging over the stationary distribution of the fast component. That is:

$$ \bar{F}(\mathbf{x}) := \sum_{y} \pi_{y}(\mathbf{x})F(\mathbf{x}, \mathbf{y}). $$

For an initial state $\mathbf{x}$ and $t \ge 0$, we call the mean field approximation the solution $\phi_t(\mathbf{x})$ of the initial value problem

$$ \frac{d}{dt}\phi_t(\mathbf{x}) = \bar{F}(\phi_t(\mathbf{x})), \quad \phi_0(\mathbf{x}) = \mathbf{x}. \qquad (4) $$

Such an approximation is also called a fluid approximation.

## 2.3 Main assumptions

As we show later, under mild regularity conditions on the transition rate functions $\alpha$, the mean field approximation captures the dynamics of $X_t$ well and with a decreasing bias of order $1/N$:

$$ \mathbb{E}[h(X_t) | X_0, Y_0 = \mathbf{x}, \mathbf{y}] - h(\phi_t(\mathbf{x})) = \frac{C_h(t)}{N} + o(1/N), $$

for a sufficiently regular $h$. This holds for any finite $t$ under assumption (A$_1$)-(A$_2$) below. It also holds for the steady-state regime $t = +\infty$ under the additional assumption (A$_3$). For the steady-state we also show that $C_h$ can be computed numerically and use it to propose a refined approximation. To obtain these results for finite time, we will assume that:

(A$_1$) The set of transitions $\mathcal{T}$ is finite and for all $\ell, \mathbf{y}' \in \mathcal{T}$, $\alpha_{\ell, \mathbf{y}'} \in \mathcal{D}^2(\mathcal{X} \times \mathcal{Y})$, where $\mathcal{D}^k(\mathcal{X} \times \mathcal{Y})$ the space of functions from $\mathcal{X} \times \mathcal{Y}$ to $\mathbb{R}$ for which the Hölder norm $\|h\|_{k,1}$ is finite$^1$.

(A$_2$) For all $\mathbf{x} \in \mathcal{X}$, the matrix $K(\mathbf{x})$ defined in Equation (3) has a unique irreducible class.

(A$_3$) The ODE (4) has a unique, exponentially stable equilibrium, that we denote by $\phi_\infty$, i.e., there exist $a, b > 0$ s.t. $\|\phi_t(\mathbf{x}) - \phi_\infty\| \le a \exp(-bt)$ for all $\mathbf{x} \in \mathcal{X}$.

Assumptions (A$_1$) and (A$_3$) are classical to ensure mean field convergence results. As stated in [15, 33], requiring that the transition rates are twice differentiable which is necessary to guarantee the existence of derivatives for the drift and the differential equation needed for the proofs of the theorems. Assumption (A$_3$) ensures the existence of a unique equilibrium point to which all trajectories of the differential equation converge. This classical assumption is essentially needed to show that the stationary distribution of the stochastic process converges to a deterministic limit, see [2]. It guarantees the stability of both the ODE and the 'slow' Poisson equation used in the proof.

By Assumption (A$_2$), we mean that for all $\mathbf{x}$ the Markov chain should have a unique subset of states that is irreducible (there can be additional states, but they should all be transient). This assumption is equivalent to assuming the uniqueness of the stationary distribution of the Markov chain induced by the generator matrix $K(\mathbf{x})$ which is essential to define the 'averaged' drift and mean field approximation. For a given $\mathbf{x}$, the stationary distribution $\pi(\mathbf{x})$ will be non-zero for all states that are in the irreducible component ($\pi_{\mathbf{y}}(\mathbf{x}) > 0$ for such $\mathbf{y}$'s) and will be zero for the others ($\pi_{\mathbf{y}}(\mathbf{x}) = 0$ for all states that are transient for $K(\mathbf{x})$). This assumption is slightly more general than assuming that $K(\mathbf{x})$ is irreducible because it allows for transient states.

We will show later that the assumptions (A$_1$) and (A$_2$) imply that $\bar{F}$ is Lipschitz continuous because they imply that $\pi(\mathbf{x})$ is Lipschitz-continuous (see Lemma 4). This implies that the mean field approximation (4) is well defined.

<footnote>$^1$ i.e., functions that are k-times differentiable with Lipschitz continuous derivatives (see a more precise definition of this norm in Definition 9 of Appendix A.1)</footnote>
# 3 MAIN RESULTS

This section includes our main results which are threefold. In 3.1 we obtain accuracy results for the mean field approximation in the transient regime. In 3.2 we obtain comparable results when the stochastic system is in its steady-state. Lastly, in 3.3 we introduce a correction term for the steady-state, give accuracy bounds and display closed form expressions of the corrections.

## 3.1 Transient Regime

Our model is a two timescale model, which makes it amenable to be analyzed by time-averaging methods such as the one used in [2, 8, 10, 29]. Such methods guarantee that the stochastic process $X_t$ converges to the solution of the ODE given by (4). Yet, most of the papers that use averaging methods do not quantify the rate at which this convergence occurs. Our first result, Theorem 1, states that the difference between the 'average' mean field approximation $\phi$ derived from the average drift $\bar{F}$ approximates the average behavior of the slow component of the two-timescale system with a bias asymptotically equivalent to $C_h(t)/N$. This result is similar to the one obtained in [13] for classical density dependent process.

**THEOREM 1.** Consider the two-timescale stochastic system $(X^{(N)}, Y^{(N)})_{t \ge 0}$ as introduced in Section 2 starting at $X_0^{(N)}, Y_0^{(N)} = x, y$, and recall that $\phi(x)$ is the solution of the ODE (4) with initial condition $x$. Further, assume $(A_1)$ and $(A_2)$. For any $h \in D^2(X)$ and $t > 0$, there exists a constant $C_h(t)$ such that:

$$ \lim_{N \to \infty} N \left( \mathbb{E}[h(X_t^{(N)}) | X_0^{(N)}, Y_0^{(N)} = x, y] - h(\phi_t(x)) \right) = C_h(t). $$

**MAIN ELEMENT OF THE PROOFS.** The full proof of this theorem is given in Section 4.3. It is decomposed into two parts that correspond to the two approximations (scaling and averaging):

* The first is to approximate the infinitesimal generator $Lh(X_t, Y_t)$ (defined in Section 4.1) of the stochastic system by a process whose drift is $D_x h(X_t) \cdot F(X_t, Y_t)$. This part uses that the rate functions $\alpha$ and $h$ are differentiable with respect to the slow variable $x$.

* The second is that the actual drift of the mean field approximation is $\bar{F}(x) = \sum_y \pi_y(x) f(x, y)$ and not the $F(x, y)$ obtained from the stochastic system. This leads us to bound a term of the form $\bar{F}(x) - F(x, y)$.

The first error term is bounded by using generator techniques similar to the ones used for classical mean field models as in [12]. It is treated in Section 4.3.1. To bound the second term, we use averaging ideas similar to the ones of [2] that are related to how fast the fast timescale converges to its stationary distribution. This is dealt with in Section 4.3.2. □

## 3.2 Steady-State Results

Theorem 1 guarantees that the mean field is a good approximation for any finite time interval. In order to obtain a similar result for the stationary case, an almost necessary condition is that the ODE (4) has a unique fixed point to which all trajectories converge [2]. To obtain the equivalent of Theorem 1, we assume that this unique fixed point is exponentially stable, as is classically done to obtain steady-state guarantees [15, 33]. This assumption is summarized in $(A_3)$ and leads to the following result.

**THEOREM 2.** Assume $(A_1) - (A_3)$, and assume that for all $N$, the stochastic system has a stationary distribution. Denote by $(X_\infty^{(N)}, Y_\infty^{(N)})$ a pair of variables having this stationary distribution. Then, for
$h \in D^2(X \times Y)$ there exists a constant $C_h$ such that:

$$ \lim_{N \to \infty} N \left( \mathbb{E}[h(X_{\infty}^{(N)}, Y_{\infty}^{(N)})] - \sum_{\mathbf{y}} \pi_{\mathbf{y}}(\phi_{\infty})h(\phi_{\infty}, \mathbf{y}) \right) = C_h. $$

### MAIN ELEMENTS OF THE PROOF.
A detailed proof is provided in Section 4.4. To prove the result, the first essential step is to generate two Poisson equations:

(1) The first is used to study $\mathbb{E}[h(X_{\infty}, Y_{\infty}) - \sum_{\mathbf{y}} \pi_{\mathbf{y}}(X_{\infty})h(X_{\infty}, \mathbf{y})]$, which is the difference between the true expectation and a hypothetical case where $Y_{\infty}$ would be independent from $X_{\infty}$ and distributed according to the stationary distribution $\pi(X_{\infty})$.

(2) The second is used to study $\mathbb{E}[\sum_{\mathbf{y}} \pi_{\mathbf{y}}(X_{\infty})h(X_{\infty}, \mathbf{y}) - \sum_{\mathbf{y}} \pi_{\mathbf{y}}(\phi_{\infty})h(\phi_{\infty}, \mathbf{y})]$, which is the error of the 'slow' process compared to the 'average' mean field.

The rest of the proof treats these Poisson equations to establish the constant $C_h$. This is done by generator comparison techniques and application of derivative bounds for the solutions of the Poisson equation. □

We note that contrary to Theorem 1, Theorem 2 allows for functions $h$ that can depend on both the slow and the fast components: $x$ and $y$. In fact, Theorem 1 would not be true if we allow functions $h$ to depend on $Y_t$ as in Theorem 2, because of the fast transitions of $Y_t$. For the steady-state, we use that $(X_{\infty}, Y_{\infty})$ starts in steady-state.

In Theorem 3 below, we will show that, for the steady-state and for functions that only depend on $x$, we can go further and propose an almost closed-form expression for the term $C_h$. This new bias term provides a refined accuracy of order $O(1/N^2)$. To prove this result, we will use that Theorem 2 allows functions $h$ which depend on $X$ and $Y$ to show that this new approximation is $O(1/N^2)$-accurate. We show below that in fact the constant $C_h$ can be computed by a numerical algorithm, and can therefore be used to define a *refined approximation*, similarly to what is done for classical mean field model in [15].

## 3.3 Steady-State Refinement

To obtain a refined approximation, we utilize ideas introduced in [15] and propose an almost-closed form expression for the term $C_h$ of Theorem 2. As we will see later in the proofs, the bias correction term is composed of two distinct components:

(1) The first one (terms $V$ and $W$) is the analogue of the $V$ and $W$ terms of [15] and corresponds to the difference between the stochastic jumps of the slow system versus having a ODE corresponding to the (non-averaged) drift $f(x, y)$.

(2) The second component (terms $T$, $S$ and $U$) corresponds to approximating $Y$ by its stationary distribution $\pi_y(x)$, and its consequence on the behavior of the slow system $X$.

The proofs of the expression for $V$ and $W$ are essentially those derived in [15]. The second correction component (terms $T$, $S$ and $U$) is related to the difference between drift and its average version. It involves studying the intricate and coupled dynamics of $X$ and $Y$ which, to the best of our knowledge, has not been studied and yields novel results.

### THEOREM 3. Assume $(A_1) - (A_3)$. Then, there exist vectors $V, S, T$ and matrices $W$ and $U$ that are solutions of linear systems of equations such that for $h \in D^2(X)$, we define $C_h$ of Theorem 2 as:

$$ C_h = \sum_i \frac{\partial h}{\partial x_i}(\phi_{\infty}) (V_i + T_i + S_i) + \frac{1}{2} \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j}(\phi_{\infty}) (W_{i,j} + U_{i,j}). $$
PROOF. This theorem is a consequence of Proposition 5 that shows the existence of the above constants, combined with Proposition 6 which provides the linear equations satisfied by V and W, and Proposition 7 which provides the equations satisfied by T, U and V. □

We would like to emphasize that the result of Theorem 3 is only valid for functions $h$ that do not depend on $\mathbf{y}$. This shows that there exists a computable constant $C_h$ such that, in steady-state:

$$ \mathbb{E}[h(X_{\infty})] = \phi_{\infty} + \frac{C_h}{N} + o(1/N). $$

Similar to [15], we call $\phi_{\infty} + \frac{C_h}{N}$ the *refined approximation* of $\mathbb{E}[h(X_{\infty})]$. As we see in our numerical experiment in Section 5, this constant is computable and can provide a more accurate approximation than the classical mean field.

In fact, combined with Theorem 2, we can show that the $o(1/N)$ term is a $O(1/N^2)$ term. This shows that the refined approximation is $O(1/N^2)$-accurate. Note, for functions $h$ that do depend on $\mathbf{y}$, the existence of such a function is guaranteed by Theorem 2 but a closed-form expression is currently out of our reach.

# 4 PROOFS

## 4.1 Stochastic Semi-Groups and Generators

Given the stochastic process of Section 2.1, we define the stochastic semi-group operator which maps a pair of initial values $(\mathbf{x}, \mathbf{y})$ and a function $h$ to the expected value of the system at time $t$. This semi-group $\psi_t$ associates to a function $h: \mathcal{X}^{(N)} \times \mathcal{Y} \to \mathbb{R}$ the function $\psi_t h: \mathcal{X}^{(N)} \times \mathcal{Y} \to \mathbb{R}$ defined as:

$$ \psi_t h(\mathbf{x}, \mathbf{y}) = \mathbb{E}[h(X_t, Y_t) \mid X_0, Y_0 = \mathbf{x}, \mathbf{y}]. \quad (5) $$

Using the right continuity of the slow-fast system, it is easy to verify that $\psi_t$ is indeed a $C_0$-semi-group (see a more precise definition in Definition 8 in Appendix A.1).

The infinitesimal generator of the stochastic process is the operator $L$ that maps a function $h \in \mathcal{D}^k(\mathcal{X} \times \mathcal{Y})$ to $Lh: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ defined by:

$$ Lh(\mathbf{x}, \mathbf{y}) = \sum_{\ell, \mathbf{y}' \in \mathcal{T}} N\alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) \left(h(x + \ell/N, \mathbf{y}') - h(\mathbf{x}, \mathbf{y})\right). \quad (6) $$

Note that $Lh$ is obtained by considering average infinitesimal change of the stochastic system starting in $(\mathbf{x}, \mathbf{y})$, i.e.,

$$ Lh(\mathbf{x}, \mathbf{y}) = \lim_{t \downarrow 0} (\psi_t h(\mathbf{x}, \mathbf{y}) - \psi_0 h(\mathbf{x}, \mathbf{y}))/t = \lim_{t \downarrow 0} (\mathbb{E}[h(X_t, Y_t) \mid X_0, Y_0 = \mathbf{x}, \mathbf{y}] - h(\mathbf{x}, \mathbf{y}))/t. $$

Similarly to the notations of the semi-group and generator of the stochastic process, for a given function $h \in \mathcal{D}^k(\mathcal{X})$, we denote by $\Phi_t h(\mathbf{x}) := h(\phi_t(\mathbf{x}))$ the $C_0$-semi-group corresponding to ODE (4). For differentiable $h$, the infinitesimal generator is given by

$$ \Lambda h(\mathbf{x}) := D_x h(\mathbf{x}) \bar{F}(\mathbf{x}). \quad (7) $$

To strengthen intuition, the time derivative of $h(\phi_t(\mathbf{x}))$ can be expressed as

$$ \frac{d}{dt} h(\phi_t(\mathbf{x})) = D_x h(\phi_t(\mathbf{x})) \bar{F}(\phi_t(\mathbf{x})) = \Phi_t \Lambda h(\mathbf{x}). \quad (8) $$

By commuting $\Lambda$ and $\Phi_t$ we have $\Lambda\Phi_t h(\mathbf{x}) = D_x(h \circ \phi_t)(\mathbf{x})\bar{F}(\mathbf{x})$ which is equal to (8) by the results of Appendix A.1. We will use this property to prove the theorems of the following section.
*Abuse of Notations and Dependence on the Fast Component.* The semi-group and generator of the stochastic system are generally defined for functions $h$ in $D^k(\mathcal{X} \times \mathcal{Y})$, i.e., functions which depend on the slow and fast component whereas the semi-group and generator of the ODE are defined for $h$ in $D^k(\mathcal{X})$ that do not depend on the fast component. What we refer to as abuse of notation is the notation we use for the mapping $\psi_t$ of a function $h \in D^k(\mathcal{X})$. $\psi_t h$ still depends on the fast component even if $h$ does not (since the evolution $X_t$ depends on the state of $Y_t$ and the initial values $(x, y)$):

$$ \psi_t h(x, y) = \mathbb{E}[h(X_t) | X_0, Y_0 = x, y], \qquad (9) $$

and

$$ Lh(x, y) = \sum_{\ell, y'} N\alpha_{\ell, y'}(x, y)(h(x + \ell/N) - h(x)). $$

In essence, this allows to use the notion of the semi-group and generator to functions of the slow process. For consistency, we do the same for the ODE: for an arbitrary $y \in \mathcal{Y}$, $\phi_t(x, y) := \phi_t(x)$ and therefore $\Phi_t h(x, y) = h(\phi_t(x, y)) = h(\phi_t(x))$ which is motivated by the abuse of notation in (9). This will merely be used in the proofs and allows focusing on integral proof ideas instead of complex notations.

## 4.2 Generator of the Fast Process and Regularity of the Poisson Equation

The generator $L$ describes the changes induced by the transitions of the fast and slow process. In our analysis, it will be useful to analyze the changes due to the jumps in $Y$ only. We denote by $L_{\text{fast}}$ the generator of the Markov chain induced by $K(x)$ that takes as input a function $h \in D^k(\mathcal{X} \times \mathcal{Y})$ and associates another function $L_{\text{fast}}h$ defined by:

$$ L_{\text{fast}}h(x, y) = \sum_{y'} K_{y, y'}(x)(h(x, y') - h(x, y)). $$

Compared to (1), $L_{\text{fast}}$ is independent of $N$ because the transition rates are rescaled by $1/N$.

Under Assumption (A$_2$), for all $x$, the matrix $K(x)$ characterizes a Markov chain on the state space $\mathcal{Y}$ that has a unique stationary distribution denoted by $\pi(x)$. $\pi_y(x)$ is the stationary probability of $y \in \mathcal{Y}$ of the Markov chain induced by the kernel $K(x)$. Subsequently, we will study the distance between the 'true' stochastic process $Y_t$ and an averaged system where the distribution of $Y_t$ is replaced by the stationary distribution $\pi(x)$. To quantify the error made when replacing $Y_t$ with the stationary distribution, we consider the following Poisson equation:

$$ h(x, y) - \sum_{y \in \mathcal{Y}} h(x, y)\pi_y(x) = L_{\text{fast}}G_h^{\text{fast}}(x, y). \qquad (10) $$

For a given function $h: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}^n$ ($n$ is arbitrary but finite), a function $G_h^{\text{fast}}: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}^n$ that satisfies the above equation is called a solution to this Poisson equation. We will have particular interest in $h(x, y) = F(x, y) \in \mathbb{R}^{d_x}$ namely when the drift $F(x, y)$ is compared to its 'average' version $\bar{F}(x) = \sum_y F(x, y)\pi(y)$.

The existence of a regular solution to this Poisson equation is guaranteed by the following Lemma 4. Note that the solution of the above Poisson equation is not unique: If $G_h^{\text{fast}}$ is a solution, then for any constant $c \in \mathbb{R}^n$, a function $G_h^{\text{fast}} + c$ is also a solution. Later in the proofs of the theorems, when we talk about 'a solution of the Poisson equation', we refer to the solution given in Lemma 4.

LEMMA 4. Assume (A$_2$). Then for all $x \in \mathcal{X}$:
(1) The Markov chain corresponding to $K(\mathbf{x})$ has a unique stationary distribution that we denote by $\pi(\mathbf{x})$. We denote by $\Pi(\mathbf{x}) = 1$ $\pi^T(\mathbf{x})$ the matrix where each line is equal to $\pi(\mathbf{x})$.

(2) The matrix $(K(\mathbf{x}) + \Pi(\mathbf{x}))$ is invertible and its inverse is a generalized inverse of $K(\mathbf{x})$.

(3) Define $K^+(\mathbf{x}) = (K(\mathbf{x}) + \Pi(\mathbf{x}))^{-1}(I - \Pi(\mathbf{x}))$, then, for all functions $h: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}^n$, $G_h(\mathbf{x}, \mathbf{y}) = \sum_{\mathbf{y}'} K_{\mathbf{y}, \mathbf{y}'}^+(\mathbf{x})h(\mathbf{x}, \mathbf{y}')$ is a solution of the Poisson Equation (10).

If, in addition, Assumption (A₁) holds, then $K^+(\mathbf{x})$ is twice differentiable in $\mathbf{x}$.

The proof is provided in Appendix B.1. Note in particular, this result implies that if $h$ is (twice) differentiable in $\mathbf{x}$ then the same holds true for $G_h^{\text{fast}}$.

## 4.3 Proof of Theorem 1 - Transient State Proof

The proof of Theorem 1 can be decomposed in two main parts. We first use a generator transformation to show that the slow system is well approximated by a system whose drift is $D_x h(X_t) F(X_t, Y_t)$ (4.3.1). This leads us to treat terms of the form $\bar{F}(X_t) - F(X_t, Y_t)$. To study them, we use the solution of the Poisson equation for the fast system (10). The second part is the more technical and novel. It is detailed in 4.3.2. Some technical lemmas are postponed to Appendix B.2.

### 4.3.1 Error due to replacing the stochastic jumps of $X_t$ by the drift.
For $h \in \mathcal{D}^2(\mathcal{X})$ recall that the $C_0$-semi-groups of the stochastic system and of the ODE are defined as

$$ \psi_s h(\mathbf{x}_0, \mathbf{y}_0) = \mathbb{E}[h(\mathbf{X}_s) | \mathbf{X}_0, Y_0 = \mathbf{x}_0, \mathbf{y}_0] \quad \text{and} \quad \Phi_t h(\mathbf{x}_0, \mathbf{y}_0) = h(\phi_t(\mathbf{x}_0)). $$

We define $v_s h(\mathbf{x}, \mathbf{y}) := \psi_s \Phi_{t-s} h(\mathbf{x}, \mathbf{y})$ and rewrite

$$ \mathbb{E}[h(X_t) - h(\phi_t(\mathbf{x})) | X_0, Y_0 = \mathbf{x}_0, \mathbf{y}_0] = v_t h(\mathbf{x}_0, \mathbf{y}_0) - v_0 h(\mathbf{x}_0, \mathbf{y}_0) = \int_0^t \frac{d}{ds} v_s h(\mathbf{x}_0, \mathbf{y}_0) ds. \quad (11) $$

To show that the last equation indeed holds true, observe that

$$ \begin{aligned} \frac{d}{ds} v_s h(\mathbf{x}_0, \mathbf{y}_0) &= \frac{d}{ds} \psi_s \Phi_{t-s} h(\mathbf{x}_0, \mathbf{y}_0) \\ &= L \psi_s \Phi_{t-s} h(\mathbf{x}_0, \mathbf{y}_0) - \psi_s \Lambda \Phi_{t-s} h(\mathbf{x}_0, \mathbf{y}_0). \end{aligned} \quad (12) $$

By regularity assumptions on the transition rates and bounded state space the above equation is finite for all $(\mathbf{x}, \mathbf{y}) \in \mathcal{X} \times \mathcal{Y}$ and time $s \ge 0$. The dominated convergence theorem thus justifies the interchange of derivative and integral and validates the last equality. As pointed out in Appendix A.1, $\psi_s$ and $L$, the stochastic semi-group and its infinitesimal generator, commute, i.e., $L\psi_s = \psi_s L$. Hence,

$$ \begin{aligned} (12) &= \psi_s (L - \Lambda) \Phi_{t-s} h(\mathbf{x}_0, \mathbf{y}_0) \\ &= \mathbb{E}[(L - \Lambda) \Phi_{t-s} h(\mathbf{X}_s, \mathbf{Y}_s) | \mathbf{X}_0, \mathbf{Y}_0 = \mathbf{x}_0, \mathbf{y}_0], \end{aligned} \quad (13) $$

where the last line follows by definition of $\psi_s$. Let $g(\mathbf{x}) = h(\phi_{t-s}(\mathbf{x}))$. $g$ is twice differentiable with respect to the initial condition $\mathbf{x}$ as it lies in $\mathcal{D}^2(\mathbb{R}^d, \mathbb{R})$ by Lemma 10. By the use of Taylor expansion rewrite $Lg$ as

$$ \begin{aligned} Lg(\mathbf{x}, \mathbf{y}) &= \sum_{\ell, \mathbf{y}' \in \mathcal{T}} N \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) \left(g(\mathbf{x} + \frac{\ell}{N}) - g(\mathbf{x})\right) \\ &= \sum_{\ell, \mathbf{y}' \in \mathcal{T}} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) D_x g(\mathbf{x}) \ell + \frac{1}{N} \sum_{\ell, \mathbf{y}' \in \mathcal{T}} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) D_x^2 g(\mathbf{x}) \cdot (\ell, \ell) + o\left(\frac{\|g\|_{2,1} \sup_{\ell} |\ell|^2}{N}\right) \\ &= Dg(\mathbf{x}) F(\mathbf{x}, \mathbf{y}) + \frac{1}{N} \sum_{\ell, \mathbf{y}' \in \mathcal{T}} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) D_x^2 g(\mathbf{x}) \cdot (\ell, \ell) + o\left(\frac{C_R \|g\|_{2,1} \sup_{\ell} |\ell|^2}{N}\right), \end{aligned} $$
with $D^2g(\mathbf{x}) \cdot (\ell, \ell) = \sum_{m,n} \frac{\partial^2 g}{\partial x_m \partial x_n}(\mathbf{x})\ell_m\ell_n$. The convergence depends further on the Hölder norm of $h(\phi(\mathbf{x}))$, bounds on the transition rates $C_R$ and jump sizes $\sup |\ell|$. The generator of the ODE related semi-group, $\Lambda$, introduced in Equation (7), maps $g$ to $\Lambda g(\mathbf{x}, \mathbf{y}) = D_x g(\mathbf{x}) \bar{F}(\mathbf{x})$. Hence, we have

$$ (L - \Lambda)g(\mathbf{x}, \mathbf{y}) = D_x g(\mathbf{x}) (F(\mathbf{x}, \mathbf{y}) - \bar{F}(\mathbf{x})) + \frac{1}{N} \sum_{\ell, y' \in \mathcal{T}} \alpha_{\ell, y'}(\mathbf{x}, \mathbf{y}) D_x^2 g(\mathbf{x}) \cdot (\ell, \ell) + o(1/N). $$

### 4.3.2 Error due to replacing the drift by the average drift.
Next, we have a closer look at the first summand of the right hand side in the above equation. Denote by $G_F^{\text{fast}}$ a solution of the Poisson equation (10) where the function “$h$” of (10) is set to the drift $F$. Since $D_x g(\mathbf{x})$ does not depend on $\mathbf{y}$ we have by definition of the Poisson equation and $L_{\text{fast}}$:

$$ D_x g(\mathbf{x}) (F(\mathbf{x}, \mathbf{y}) - \bar{F}(\mathbf{x})) = D_x g(\mathbf{x}) L_{\text{fast}} G_F^{\text{fast}}(\mathbf{x}, \mathbf{y}) = L_{\text{fast}} D_x g(\mathbf{x}) G_F^{\text{fast}}(\mathbf{x}, \mathbf{y}). \quad (14) $$

Combining the above with equation (13) and plugging everything into the integral of equation (11), we get:

$$ \mathbb{E}[h(\mathbf{X}_t) - h(\phi_t(\mathbf{x}_0))] = \int_0^t \mathbb{E}[L_{\text{fast}} D_x(h \circ \phi_{t-s})(\mathbf{X}_s) \cdot G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s)] \\ + \frac{1}{N} \mathbb{E}[\sum_{\ell, y' \in \mathcal{T}} \alpha_{\ell, y'}(\mathbf{X}_s, \mathbf{Y}_s) D_x^2(h \circ \phi_{t-s})(\mathbf{X}_s) \cdot (\ell, \ell)] ds + o(1/N), \quad (15) $$

where we suppress the conditioning on the initial values $\mathbf{X}_0, \mathbf{Y}_0 = \mathbf{x}_0, \mathbf{y}_0$. Let $H_s(\mathbf{x}, \mathbf{y}) := \mathbb{E}[D_x(h \circ \phi_{t-s})(\mathbf{X}_s) G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s) | \mathbf{X}_0, \mathbf{Y}_0 = \mathbf{x}, \mathbf{y}]$. Applying Lemma 11 with $g_s(\mathbf{X}_s, \mathbf{Y}_s) = D_x(h \circ \phi_{t-s})(\mathbf{X}_s) G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s)$ implies that

$$ 0 = \frac{1}{N} H_t(\mathbf{x}, \mathbf{y}) - \frac{1}{N} H_0(\mathbf{x}, \mathbf{y}) - \int_0^t \mathbb{E} \left[ \frac{1}{N} D_x \left( D_x h(\phi_{t-s}(\mathbf{X}_s)) \cdot \bar{F}(\phi_{t-s}(\mathbf{X}_s)) \right) \cdot G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s) \right] ds \\ + \mathbb{E} \left[ \frac{1}{N} L D_x (h \circ \phi_{t-s})(\mathbf{X}_s) G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s) \right] ds. \quad (16) $$

As $h$ and $\phi_{t-s}$ are twice continuously differentiable, by compactness of $\mathcal{X}$, and as $G_F^{\text{fast}}$ is finite, $D_x(D_x h(\phi_{t-s}(\mathbf{X}_s))\bar{F}(\phi_{t-s}(\mathbf{X}_s)))G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s)$ is bounded. Moreover, by Lemma 12 we have

$$ (L_{\text{fast}} - \frac{1}{N}L)D_x(h \circ \phi_{t-s})(\mathbf{X}_s)G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s) \\ = \frac{1}{N} \sum_{\ell, y'} \alpha_{\ell, y'}(\mathbf{X}_s, \mathbf{Y}_s) \ell D_x(D_x(h \circ \phi_{t-s})(\mathbf{X}_s)G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{y}')) + o(1/N). $$

Define

$$ C_h(t) = \int_0^t \mathbb{E} \left[ \sum_{\ell, y' \in \mathcal{T}} \alpha_{\ell, y'}(\mathbf{X}_s, \mathbf{Y}_s) D_x^2 h(\phi_{t-s}(\mathbf{X}_s)) \cdot (\ell, \ell) \right] ds \\ + \int_0^t \mathbb{E} \left[ D_x(D_x h(\phi_{t-s}(\mathbf{X}_s)) \cdot \bar{F}(\phi_{t-s}(\mathbf{X}_s))) \cdot G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{Y}_s) \right] ds \\ + \int_0^t \mathbb{E} \left[ \sum_{\ell, y'} \alpha_{\ell, y'}(\mathbf{X}_s, \mathbf{Y}_s) \ell D_x(D_x(h \circ \phi_{t-s})(\mathbf{X}_s) G_F^{\text{fast}}(\mathbf{X}_s, \mathbf{y}')) \right] ds + \mathbb{E}[H_t(\mathbf{X}_t, \mathbf{Y}_t)] - H_0(\mathbf{x}, \mathbf{y}), $$
which lets us rewrite

$$N (\mathbb{E}[h(X_t)] - h(\phi_t(\mathbf{x}))) = C_h(t) + o(1).$$

Taking the limit $N \to \infty$ concludes the proof.

## 4.4 Proof of Theorem 2 - Steady-State Proof

Considering $\mathbb{E}[h(X_\infty, Y_\infty) - \sum_y \pi_y(\phi_\infty)h(\phi_\infty, y)]$, first, add and subtract the term $\mathbb{E}[\sum_y \pi_y(X_\infty)h(X_\infty, y)]$ into the equation, yielding

$$\mathbb{E}[h(X_\infty, Y_\infty) - \sum_y \pi_y(\phi_\infty)h(\phi_\infty, y)] = \mathbb{E}[h(X_\infty, Y_\infty) - \sum_y \pi_y(X_\infty)h(X_\infty, y)] \quad (A)$$

$$+ \mathbb{E}[\sum_y \pi_y(X_\infty)h(X_\infty, y) - \sum_y \pi_y(\phi_\infty)h(\phi_\infty, y)]. \quad (B)$$

Treating the two terms (A) and (B) separately, we define $G_h^{\text{fast}}$ and $G_h^{\text{slow}}$ as the solutions to the Poisson equations

$$L_{\text{fast}} G_h^{\text{fast}}(\mathbf{x}, \mathbf{y}) = h(\mathbf{x}, \mathbf{y}) - \sum_{\mathbf{y}} \pi_{\mathbf{y}}(\mathbf{x})h(\mathbf{x}, \mathbf{y}), \quad (\text{fast Poisson Equation})$$

$$\Lambda G_h^{\text{slow}}(\mathbf{x}) = \sum_{\mathbf{y}} \pi_{\mathbf{y}}(\mathbf{x})h(\mathbf{x}, \mathbf{y}) - \sum_{\mathbf{y}} \pi_{\mathbf{y}}(\phi_{\infty})h(\phi_{\infty}, \mathbf{y}). \quad (\text{slow Poisson Equation})$$

Recall that the existence and regularity properties of the function $G_h^{\text{fast}}$ are investigated in Lemma 4. For $G_h^{\text{slow}}$, it is known (e.g., [12]) that the exponential stability of the unique fixed point ($A_3$) along with the smoothness of $F$ guarantees that $G_h^{\text{slow}} \in D^2(X)$.

Next we use that taking the expectation of the generator $L$ applied to an arbitrary function $g$ over the stationary distribution of $(X_\infty, Y_\infty)$ is zero, i.e., $\mathbb{E}[Lg(X_\infty, Y_\infty)] = 0$. By using Lemma 12, it follows

$$\begin{aligned} (A) &= \mathbb{E}[L_{\text{fast}} G_h^{\text{fast}}(X_\infty, Y_\infty) - \frac{1}{N} L G_h^{\text{fast}}(X_\infty, Y_\infty)] \\ &= \mathbb{E}\left[\frac{1}{N} \sum_{\ell, y'} \alpha_{\ell, y'}(X_\infty, Y_\infty) \ell D_x G_h^{\text{fast}}(X_\infty, y') + o(1/N)\right]. \end{aligned}$$

For the second term, we apply the 'slow Poisson Equation' equation and subtract $\mathbb{E}[LG_h^{\text{slow}}(X_\infty, Y_\infty)] = 0$ which yields

$$(B) = \mathbb{E}[\Lambda G_h^{\text{slow}}(X_\infty) - LG_h^{\text{slow}}(X_\infty, Y_\infty)]. \quad (17)$$

Note that even if $G_h^{\text{slow}}$ only depends on $x$, the function $LG_h^{\text{slow}}$ does depend on $x$ and $y$ because the rate functions that appear in the generator $L$ do depend on $x$ and $y$. By using the Taylor expansion, $LG_h^{\text{slow}}$ equals

$$\begin{aligned} LG_h^{\text{slow}}(\mathbf{x}, \mathbf{y}) &= \sum_{\ell, y'} N \alpha_{\ell, y'}(\mathbf{x}, \mathbf{y}) (G_h^{\text{slow}}(\mathbf{x} + \frac{\ell}{N}) - G_h^{\text{slow}}(\mathbf{x})) \\ &= \sum_{\ell, y'} N \alpha_{\ell, y'}(\mathbf{x}, \mathbf{y}) D_x G_h^{\text{slow}}(\mathbf{x}) \frac{\ell}{N} + \frac{1}{2} \sum_{\ell, y'} N \alpha_{\ell, y'}(\mathbf{x}, \mathbf{y}) D_x^2 G_h^{\text{slow}}(\mathbf{x}) \frac{(\ell, \ell)}{N^2} + o\left(\frac{\|\ell\|^2}{N}\right). \end{aligned}$$
By application of the definition of $F(\mathbf{x}, \mathbf{y})$, the first term of the right hand side is equal to $F(\mathbf{x}, \mathbf{y})D_x G_h^{\text{slow}}(\mathbf{x})$. Moreover, by definition, $\Lambda G_h^{\text{slow}}(\mathbf{X}_\infty) = D_x G_h^{\text{slow}}(\mathbf{x})\bar{F}(\mathbf{x})$. This shows that

$$ \mathbb{E}[D_x G_h^{\text{slow}}(\mathbf{X}_\infty)(\bar{F}(\mathbf{X}_\infty) - F(\mathbf{X}_\infty, \mathbf{Y}_\infty)) + \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{X}_\infty, \mathbf{Y}_\infty)(\frac{1}{2}D_x^2 G_h^{\text{slow}}(\mathbf{X}_\infty) \cdot \frac{(\ell, \ell)}{N} + o(\frac{\|\ell\|}{N^2}))]. $$

As $G_h^{\text{slow}}$ has a bounded second derivative, the second term is of order $O(1/N)$. However, to bound the first summand we need to investigate $\bar{F}(\mathbf{x}) - F(\mathbf{x}, \mathbf{y}) = \sum_y \pi_{\mathbf{y}'}(\mathbf{x})F(\mathbf{x}, \mathbf{y}') - F(\mathbf{x}, \mathbf{y})$. Falling back on ('fast' Poisson Equation) allows two write

$$ D_x G_h^{\text{slow}}(\mathbf{x}) (\bar{F}(\mathbf{x}) - F(\mathbf{x}, \mathbf{y})) = -D_x G_h^{\text{slow}}(\mathbf{x}) L_{\text{fast}} G_F^{\text{fast}}(\mathbf{x}, \mathbf{y}) = -L_{\text{fast}} D_x G_h^{\text{slow}}(\mathbf{x}) G_F^{\text{fast}}(\mathbf{x}, \mathbf{y}). $$

The last equality holds due to the definition of $L_{\text{fast}}$, and $D_x G_h^{\text{slow}}(\mathbf{x})$ which depends only on $\mathbf{x}$. By adding the zero term $\mathbb{E}[\frac{1}{N}LD_x G_h^{\text{slow}}(\mathbf{X}_\infty)G_F^{\text{fast}}(\mathbf{X}_\infty, \mathbf{Y}_\infty)]$ and using Lemma 12. We see that

$$ \begin{aligned} & \mathbb{E}[-L_{\text{fast}} D_x G_h^{\text{slow}}(\mathbf{X}_\infty) G_F^{\text{fast}}(\mathbf{X}_\infty, \mathbf{Y}_\infty) + \frac{1}{N} L D_x G_h^{\text{slow}}(\mathbf{X}_\infty) G_F^{\text{fast}}(\mathbf{X}_\infty, \mathbf{Y}_\infty)] \\ & = \mathbb{E}[\frac{1}{N} \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{X}_\infty, \mathbf{Y}_\infty) \ell D_x (D_x G_h^{\text{slow}}(\mathbf{X}_\infty) G_F^{\text{fast}}(\mathbf{X}_\infty, \mathbf{y}')) + o(1/N)]. \end{aligned} $$

Put together, the error term of (B) is

$$ \mathbb{E}[\frac{1}{N} \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{X}_\infty, \mathbf{Y}_\infty) \ell D_x (D_x G_h^{\text{slow}}(\mathbf{X}_\infty) G_F^{\text{fast}}(\mathbf{X}_\infty, \mathbf{y}')) + \frac{1}{2N} \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{X}_\infty, \mathbf{Y}_\infty) D_x^2 G_h^{\text{slow}}(\mathbf{X}_\infty)(\ell, \ell)]. \quad (18) $$

By defining $C_h$ as the sum of the residual term $\mathbb{E}[\sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{X}_\infty, \mathbf{Y}_\infty) \ell D_x G_h^{\text{fast}}(\mathbf{X}_\infty, \mathbf{y}')]$ for (A) and equation (18) by $N$ scaled, the proof concludes.

## 4.5 Proof of Theorem 3 (Refinement Theorem, and Closed form Expressions)

In this section we first decompose the constant $C_h$ in two terms in Proposition 5. We then study these two terms in Proposition 6 and 7 where we obtain the closed form expressions for the correction terms which allow to numerically obtain the corrections. We define two functions $Q(x, y)$ and $J_h(x, y)$ (the latter is defined for a function $h \in \mathcal{D}^2(\mathcal{X})):

$$ Q(\mathbf{x}, \mathbf{y}) = \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) \ell \ell^T, $$

$$ J_h(\mathbf{x}, \mathbf{y}) = - \sum_{\ell, \mathbf{y}'} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) D_x (D_x G_h^{\text{slow}}(\mathbf{x}) G_F^{\text{fast}}(\mathbf{x}, \mathbf{y}')) \ell, $$

and we denote by $\bar{J}_h(x) = \sum_y \pi_y(x)J(x, y)$ and $\bar{Q}(x) = \sum_y \pi_y(x)Q(x, y)$ their “average” version. The following proposition holds.

**PROPOSITION 5.** Assume $(A_1)-(A_3)$. Then, for any $h \in \mathcal{D}^2(\mathcal{X})$, we have:

$$ N(\mathbb{E}[h(\mathbf{X}_\infty)] - h(\phi_\infty)) = \frac{1}{2} D_x^2 G_h^{\text{slow}}(\phi_\infty) \bar{Q}(\phi_\infty) + \bar{J}_h(\phi_\infty) + o(1). \quad (19) $$

**PROOF.** Let $h \in \mathcal{D}^2(\mathcal{X})$ and let us denote $G_h^{\text{slow}}$ the solution of ('slow' Poisson Equation). As $h$ does not depend on $y$, this function is such that for any $x$:

$h(x) - h(\phi_\infty) = \Lambda G_h^{\text{slow}}(x) = D_x G_h^{\text{slow}}(x) \bar{F}(x)$. The following steps are similar to the ones for (B) in the proof of Theorem 2. Hence, applying this for $\mathbf{X}_\infty$ and taking the expectation, we get:

$$ N\mathbb{E}[h(\mathbf{X}_\infty) - h(\phi_\infty)] = N\mathbb{E}[\Lambda G_h^{\text{slow}}(\mathbf{X}_\infty)] = N\mathbb{E}[D_x G_h^{\text{slow}}(\mathbf{X}_\infty)\bar{F}(\mathbf{X}_\infty)]. $$
Recall that for any bounded function $g \in \mathcal{D}(\mathcal{X} \times \mathcal{Y})$, $\mathbb{E}[Lg(\mathcal{X}_{\infty}, \mathcal{Y}_{\infty})] = 0$. Hence,

$$ N\mathbb{E}[D_x G_h^{\text{slow}}(\mathcal{X}_\infty) \bar{F}(\mathcal{X}_\infty) - LG_h^{\text{slow}}(\mathcal{X}_\infty, \mathcal{Y}_\infty)] $$

Similarly to what we do to prove Theorem 2, we plug the definition of $L$ in the above equation and use a Taylor expansion to show that this equals

$$ \begin{aligned} & N\mathbb{E}[D_x G_h^{\text{slow}}(\mathcal{X}_\infty) \bar{F}(\mathcal{X}_\infty) - \sum_{\ell, y'} N\alpha_{\ell, y'}(\mathcal{X}_\infty, \mathcal{Y}_\infty)(G_h^{\text{slow}}(\mathcal{X}_\infty + \frac{\ell}{N}) - G_h^{\text{slow}}(\mathcal{X}_\infty))] \\ &= N\mathbb{E}[D_x G_h^{\text{slow}}(\mathcal{X}_\infty) \bar{F}(\mathcal{X}_\infty) - \sum_{\ell, y'} N\alpha_{\ell, y'}(\mathcal{X}_\infty, \mathcal{Y}_\infty)(D_x G_h^{\text{slow}}(\mathcal{X}_\infty) \frac{\ell}{N} + \frac{1}{2} D_x^2 G_h^{\text{slow}}(\mathcal{X}_\infty) \frac{(\ell, \ell)}{N^2} + o(\|\ell\|))] \\ &= N\mathbb{E}[D_x G_h^{\text{slow}}(\mathcal{X}_\infty)(\bar{F}(\mathcal{X}_\infty) - F(\mathcal{X}_\infty, \mathcal{Y}_\infty)) - \frac{1}{2N} D_x^2 G_h^{\text{slow}}(\mathcal{X}_\infty) Q(\mathcal{X}_\infty, \mathcal{Y}_\infty)] + o(1). \end{aligned} \quad (20) $$

The rest of the proof follows using

* Lemma 13 to show $N\mathbb{E}[D_x G_h^{\text{slow}}(\mathcal{X}_\infty)(\bar{F}(\mathcal{X}_\infty) - F(\mathcal{X}_\infty, \mathcal{Y}_\infty))] = \bar{J}(\phi_\infty) + o(1)$;

* Theorem 2 which states $\mathbb{E}[\frac{1}{2}D_x^2 G_h^{\text{slow}}(\mathcal{X}_\infty)Q(\mathcal{X}_\infty, \mathcal{Y}_\infty)] = \frac{1}{2}D_x^2 G_h^{\text{slow}}(\phi_\infty)\bar{Q}(\phi_\infty) + o(1)$.

Applying the above equations to (20) implies the statement of Proposition 5. $\square$

In the rest of the section, we show that the quantity $\frac{1}{2}D_x^2 G_h^{\text{slow}}(\phi_\infty)\bar{Q}(\phi_\infty) + \bar{J}(\phi_\infty)$ can be easily computed numerically by solving linear systems of equations. As shown in Proposition 6 and 7, we obtain five different correction terms:

* The first two $V \in \mathbb{R}^{d_x}$, $W \in \mathbb{R}^{d_x \times d_x}$ are closely related to the ones obtained in [15][Theorem 3.1] and derived from the term $\frac{1}{2}D_x^2 G_h^{\text{slow}}(\phi_\infty)\bar{Q}(\phi_\infty)$. This term is essentially identical to the refinement term of [15].

* In Proposition 7 we derive three other refinement terms $S, T \in \mathbb{R}^{d_x}$ and $U \in \mathbb{R}^{d_x \times d_x}$ which give closed form descriptions of $\bar{J}(\phi_\infty)$. These terms are novel and take into account the difference of the average drift $\bar{F}$ and the actual drift $F$ obtained from the stochastic system. The proof is based on exact expressions for the Poisson equations $G^{\text{fast}}, G^{\text{slow}}$ and relies on Lemma 4.

### 4.5.1 Correction terms V and W.
To obtain the first result, let A and B be the Jacobian and Hessian matrix of the 'average' drift, i.e.,

$$ A_{i,j} = \frac{\partial \bar{F}_i}{\partial x_j}(\phi_\infty) \quad \text{and} \quad B_{i,k_1,k_2} = \frac{\partial^2 \bar{F}_i}{\partial x_{k_1} \partial x_{k_2}}(\phi_\infty). \quad (21) $$

By the exponential stability of the fixed point (Assumption (A$_3$)), the matrix A is invertible and the Lyapunov equation $AW + WA^T + \bar{Q} = 0$ has a unique solution. We denote by W its solution and define the vector V as

$$ V_i = -\frac{1}{2} \sum_j (A^{-1})_{i,j} \sum_{k_1, k_2} B_{j,k_1, k_2} W_{k_1, k_2}. $$

**PROPOSITION 6.** Assume (A$_1$) - (A$_3$) and $h \in \mathcal{D}^2(\mathcal{X})$. Then:

$$ \frac{1}{2} D_x^2 G_h^{\text{slow}}(\phi_\infty) \bar{Q}(\phi_\infty) = \sum_i \frac{\partial h}{\partial x_i} V_i + \frac{1}{2} \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j} W_{i,j}. $$

**PROOF.** The proposition is a direct consequence of the results of [15][Theorem 3.1] which we apply to the function $\bar{F}(x) = \sum_y \pi_y(x)F(x, y)$. $\square$
### 4.5.2 Correction terms T, S, U.

**PROPOSITION 7.** Assume $(A_1) - (A_3)$ and $h \in D^2(X)$. The closed form solution of $\bar{J}$ at the equilibrium point $\phi_\infty$ is given by

$$ \bar{J}_h(\phi_\infty) = \sum_i \frac{\partial h}{\partial x_i}(\phi_\infty)(T_i + S_i) + \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j}(\phi_\infty)U_{i,j}. $$

U is the unique solution to the Sylvester equation

$$ AX + XA^T = -O, \quad \text{with} \\ O = \sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) K_{y', \ell}^+(\phi_\infty) F(\phi_\infty, \ell) \ell^T, $$

where A, B are the Jacobian, Hessian of the average drift as defined in (21) and $K^+$ as given in Lemma 4. T, S are defined by

$$ T_i := \sum_j A_{i,j}^{-1} \sum_{k_1,k_2} B_{j,k_1,k_2} U_{k_1,k_2}, $$

$$ S_i := \sum_k A_{i,k}^{-1} \sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) \left( \sum_{y'} F_k(\phi_\infty, y') \nabla_x^T K_{y', y'}^+(\phi_\infty) \ell + K_{y', y'}^+(\phi_\infty) \nabla_x^T F_k(\phi_\infty, y') \ell \right). $$

The proof of this proposition is given in Appendix B.5.

## 5 EXAMPLE: CSMA MODEL

To illustrate our results, we consider the unsaturated CSMA random-access networks studied in [10]. In this paper, the authors use a two-scale model to study the performance of a CSMA algorithm with many nodes. The slow process corresponds to the arrival of jobs and the fast process corresponds to the activation and deactivation of nodes. The authors of this paper derive a mean field approximation and show that it is asymptotically exact. With our methods, we go two steps further:

* Theorem 1 and 2 show that not only the mean field approximation is asymptotically exact but also that the error is only of order $O(1/N)$.

* By using Theorem 3, we can compute a refinement term. Our numerical example shows that, similarly to what happens for classical one-scale mean field models [15], this refinement term is extremely accurate. It is much more accurate than the classical mean field approximation when the studied system is not too large.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Two examples of interference graphs.</figcaption></figure>
## 5.1 Model Description

We consider a model with C server types, with N statistically identical servers for each class. All servers communicate through a wireless medium using a random-access protocol and have a finite buffer size B. The C classes form an interference graph G = (C, E) with C = {1, ..., C} the classes and E the network specific edges. (see for instance Figure 1 for examples of graphs with three or five classes). This interference graph indicates that two servers cannot transmit simultaneously if they either are of the same class or belong to an adjacent class. For each class c ∈ C, we will denote by $y_c \in \{0, 1\}$ a variable that equals 1 if a node of class c is transmitting and 0 otherwise. We denote by $\mathcal{Y}$ the set of possible activation vectors $y \in \{0, 1\}^C$. It is equal to the set of the independent sets of the graph, i.e., all activity vectors for which an active node has only inactive neighbors. For instance, for the graph with three classes which are linearly connected as shown in Figure 1a, the set of feasible states is given by

$$ \mathcal{Y} = \{(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 0, 1)\}. $$

Any node can turn active if there are no neighboring active nodes. Once its transmission is finished, it changes back to an idle state. For a given class c, we define $\mathcal{Y}_c^+ := \{y \in \mathcal{Y} : y_c = 1\}$ as the subset of states for which node c is active, and by $\mathcal{Y}_c^- := \{y \in \mathcal{Y} : y_c = y_d = 0 \forall d \text{ s.t. } (c, d) \in E\}$ the subset of states from which node c can turn active, i.e., all neighboring nodes are inactive. For instance, for the linear 3 node graph of Figure 1a and the class c = 1 this yields the following sets

$$ \mathcal{Y}_1^+ = \{(1, 0, 0), (1, 0, 1)\}, $$

$$ \mathcal{Y}_1^- = \{(0, 0, 0), (0, 0, 1)\}. $$

As in [10], in this non-saturated model, we consider that if a node is in class c, new packets arrive to this node at rate $\lambda_c > 0$. If a node has a packet to transmit and no neighboring node is transmitting, then this node becomes active at rate $v_c > 0$. We assume that a transmission from a node of class c takes an exponential time of duration $1/(N\mu_c) > 0$, after which the packet leaves the system.

## 5.2 Two-scale model representation

The model as described above fits in our two-timescale representation. To see why, for each class $c \in C$ and buffer size $b \in \{0...B\}$, we define $X_{t,(c,b)}^{(N)}$ as the fraction of servers of class c that have at least b jobs in their queue at time t. We denote by $X_t$ the vector of all possible $X_{t,(c,b)}$ for all $c \in C$ and $b \in \{0...B\}$. The fast component $Y_t$ is the activation at time t: $Y_{t,(c)} = 1$ if a node of class c is transmitting at time t and 0 otherwise.

Using this representation, we characterize the possible transitions. Given a state pair $(x, y)$, the transitions are represented as a transition vector of the form $(\ell/N, y')$ and a corresponding transition rate $N\alpha_{\ell,y'}(x, y)$ such that the state $(x, y)$ jumps to $(x + \ell/N, y')$ at rate $N\alpha_{\ell,y'}(x, y)$. The transitions can be distinguished into three types:

* Arrival of a packet to a server of class $c \in C$:

$$ \begin{array}{lll} \left( \frac{e_{c,i}}{N}, \mathbf{y} \right) & \text{at rate} & N\lambda_c(1 - x_{c,i}) \quad \text{for } i = 1, \\ \left( \frac{\ell e_{c,i}}{N}, \mathbf{y} \right) & \text{at rate} & N\lambda_c(x_{c,i-1} - x_{c,i}) \quad \text{for } 2 \le i \le B. \end{array} $$

* Back-off of a server of class $c \in C$ with at least one packet if the class activity vector allows the back-off, i.e., $\mathbf{y} \in \mathcal{Y}_c^-$:

$$ \begin{array}{lll} \left(-\frac{e_{c,i}}{N}, \mathbf{y} + e_c\right) & \text{at rate} & N v_c(x_{c,i} - x_{c,i+1}) \quad \text{for } 1 \le i < B, \\ \left(-\frac{\ell e_{c,i}}{N}, \mathbf{y} + e_c\right) & \text{at rate} & N v_c x_{c,i} \quad \text{for } i = B. \end{array} $$
* Transmission completion of an active node of class $c \in C$, i.e., $\mathbf{y} \in \mathcal{Y}_c^+$:

$$ (0, \mathbf{y} - e_c) \text{ at rate } N\mu_c. $$

We denote the set of all possible transition vectors $(\ell/N, \mathbf{y}')$ from a state pair $(\mathbf{x}, \mathbf{y})$ by $\mathcal{T}(\mathbf{x}, \mathbf{y})$.

## 5.3 Steady-State distribution $\pi_y(x)$ and average drift

By using the above transition definitions, the matrix $K$ is given by:

* If $\mathbf{y} \in \mathcal{Y}_c^-$, then $K(\mathbf{y}, \mathbf{y} + e_c) = Nv_c x_{c,1}$,

* If $\mathbf{y} \in \mathcal{Y}_c^+$, then $K(\mathbf{y}, \mathbf{y} - e_c) = N\mu_c$,

all other entries of the matrix being 0.

This representation is used by the authors of [10] to derive the product-form stationary distribution for a fixed server state $\mathbf{x}$. This product form is closely related to the product-form stationary distribution of saturated networks as found in [3, 31, 32]: The quantity $\pi_y(x)$ is calculated as follows:

$$ \pi_y(\mathbf{x}) := \frac{Z(\mathbf{x}, \mathbf{y})}{Z(\mathbf{x})}, \quad \text{with} \quad Z(\mathbf{x}, \mathbf{y}) = \prod_{c \in C} \left( \frac{v_c}{\mu_c} x_{c,1} \right)^{y_c}, \quad Z(\mathbf{x}) = \sum_{\mathbf{y} \in \mathcal{Y}} Z(\mathbf{x}, \mathbf{y}). $$

Following our definition of Section 2.2 the drift and its average version are generically defined by:

$$ F(\mathbf{x}, \mathbf{y}) = \sum_{(t, \mathbf{y}') \in \mathcal{T}(\mathbf{x}, \mathbf{y})} \alpha_{\ell, \mathbf{y}'}(\mathbf{x}, \mathbf{y}) \ell, \quad \text{and} \quad \bar{F}(\mathbf{x}) = \sum_{\mathbf{w} \in \mathcal{Y}} \pi_y(\mathbf{x}) F(\mathbf{x}, \mathbf{y}). $$

For the drift $F(\mathbf{x}, \mathbf{y})$ of the random access model this leads to the closed form expression

$$ F(\mathbf{x}, \mathbf{y}) = \sum_{c \in C} \left( \sum_{1 < i \le B} e_{c,i} \lambda_c (x_{c,i-1} - x_{c,i}) + e_{c,1} \lambda_c (1 - x_{c,1}) - 1_{\{y \in \mathcal{Y}_c\}} \left( \sum_{1 \le i < B} e_{c,i} v_c (x_{c,i} - x_{c,i+1}) + e_{c,B} v_c x_{c,B} \right) \right). $$

It should be clear that assumption (A$_1$) holds in our case because the rates given in Section 5.2 are all continuous in $\mathbf{x}$ (in fact they are all linear). Moreover, the model also satisfies Assumption (A$_2$): For a given $\mathbf{x}$, the set of irreducible states for $K(\mathbf{x})$ contains all the feasible activation vectors $\mathbf{y}$ such that $y_c = 0$ if $x_{c,1} = 0$. The condition $x_{c,1} = 0$ implies that the nodes of class $c$ do not have any packets to transmit. The situation of Assumption (A$_3$) is more complicated. To the best of our knowledge, there has not been a complete stability characterization for the unsaturated random access CSMA model. Cecchi et al. show in [10] that in the case of a complete interference graph stability conditions can be derived which assure global exponential stability. They further conjectured that similar results hold for general interference graphs. In our analysis, we assume that (A$_3$) holds.

## 5.4 Numerical results

To study the accuracy of the mean field approximation and the refined term proposed in Theorem 2, we implemented a Python library<fnref n="2" /> to simulate the system and compute the mean field approximation and the refinement term. This library is generic and can take as an input any instance of the model which we defined above. For instance, in the Code Cell 1, we illustrate how to use this library to construct a model where the interference graph is as in Figure 1b, the rates are $\lambda = [.5, .7, .7, .6, .4]$, $v = [4, 3, 3, 3, 3]$, $\mu = [3, 3, 2, 4, 2]$, and the buffer size is equal to 10. This cell shows how to initialize the 5 node model and obtain the approximation and refinement from

<footnote n="2">²https://gitlab.inria.fr/sallmeie/bias-and-refinement-of-multiscale-mean-field-models</footnote>
our implementation. We also perform the same experiments with the linear 3 node model, for which we provide the results in Appendix D. Note that the results are qualitatively very similar.

Code Cell 1. Initialization and Computation of Mean Field and Refinements.

```R
# Graph structure (this is the five node example)
G = np.array([[0,1,1,1,0],
              [1,0,0,0,0],
              [1,0,0,0,1],
              [1,0,0,0,1],
              [0,0,1,1,0]])

# rates & buffer size
lambda = np.array([.5,.7,.7,.6,.4])
nu = np.array([4,3,3,3,3])
mu = np.array([3,3,2,4,2])
buffer_size = 10

# We define the model, compute a trajectory and the refinement term.
csma = symbolic_CSMA(nu, mu, lambda, G, buffer_size)
T, X = csma.ode(time=200) # mean field ODE
v, s, t, w, u = csma.compute_refinements(x[-1]) # steady-state refinement
```

In order to compute the refinement terms, the library needs to compute various derivatives (of the drift or of the matrix $K^+(x)$). To implement this, we rely on symbolic differentiation provided by the sympy library [24]. As we see later in Table 1, the use of the symbolic differentiation is the performance bottleneck of our implementation. In Appendix C we furthermore show how to obtain the stationary distribution and the derivative of $K^+(x)$ numerically.

*Transient regime and illustration of Theorem 1.* To illustrate the accuracy of the mean field, we use the 5 node model described in the Code Cell 1. We first simulate the CSMA model and compare it with the mean field ODE. The results are reported in Figure 2 where we plot the mean field approximation against a sample mean $E[X_t]$ derived from 1000 simulations. Initially, all servers are idle. The plot shows the share of servers of class $c = 3$ that have at least one job, that is $E[X_{t,(3,1)}]$. We compare the results for a model with $N = 10$ servers per class (top right), $N = 20$ (bottom right), or $N = 50$ (left). We observe that in all cases, the evolution of the stochastic system is very well predicted by the mean field approximation. To quantify this more precisely, each plot contains a zoom on the trajectory between the time $t = 8$ to $t = 13$. These zooms show that for $N = 50$, the quantity $E[X_{t,(3,1)}^{(N)}]$ is almost indistinguishable from the mean field approximation. For $N = 10$ or $N = 20$, the estimated average is slightly above the mean field curve, but the confidence intervals remain almost equal to the error.

*Steady-state and refined accuracy.* While Theorems 1 and 2 provide a guarantee on the accuracy of the mean field approximation, Theorem 3 shows that it is possible to compute an approximation that is more accurate than the original mean field approximation. We illustrate this in Figure 3 where we show the steady-state average queue lengths for the same 5 node graph. The sample mean and confidence interval are computed from 40 steady-state samples which again are obtained from independent time-averages of $7.5 \times 10^6$ events of the Markov chain after a warm-up of $2.5 \times 10^6$ events. For a class $c$ and a buffer size $b$, the quantity $\sum_{b=1}^B E[X_{(c,b)}]$ is equal to the steady-state average queue length of each server of class $c$. In Figure 3, we consider different values of $N$, and calculate the average queue length by the following three methods:

* By using a stochastic simulator of the original CSMA model.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of the transient behavior of the CSMA model. We compare the ‘average’ mean field and stochastic simulations for three different scaling parameters: $N = 10$, $N = 20$ and $N = 50$.</figcaption></figure>

* By using the fixed point of the mean field approximation: $\sum_{b=1}^{B}(\phi_{\infty})_{(c,b)}$.

* By computing the refinement term, $C$, of Theorem 3 and using $\sum_{b=1}^{B}(\phi_{\infty})_{(c,b)} + (C)_{(c,b)} / N$.

When looking at the scale of the y-axis, we see that in all cases, the accuracy of the mean field approximation is already quite good. More importantly, we also observe that in all cases, the refined approximation seems almost exact: For all considered cases, the refined approximation lies within the 95 percent confidence interval of the simulations and seems to work well even for a small number of servers, $N \approx 10, 20$. This result is similar to the one observed for one-timescale mean field models in [15].

<figure><img src="image_4.png" /><figcaption>Fig. 3. Queue Length Distribution for the 5 Node Graph of Figure (1b).</figcaption></figure>

*Computation Time.* While the previous figure shows that the refined approximation provides an increase in accuracy for small values of $N$, it comes at the cost of an increase in computation time because one needs to compute the various derivatives of the rate functions and to solve a new linear systems of equations. In order to quantify the additional computation time, we measure the time taken by our implementation to compute the refinement terms which are reported in Table 1. We compare the 5 node model studied before and a 3-node model whose interference graph is as in Figure 1a. We observe that the time taken to compute the refinement term is significant, in particular for the 5-node model. Yet, when looking more carefully at what takes time, we realize that most of the computation time is taken by the symbolic differentiation. Indeed, to simplify our implementation, we used the automated differentiation method of `sympy`. While this yields simplifications for the implementation, we encountered that it massively slows down the refinement computation times. Through code profiling it showed that around 95 percent of the computing time is taken by `sympy` methods such as differentiation and evaluation of symbolic expressions. For smaller interference graphs, e.g., linear 2/3 node graphs, this effect is not limiting.
For larger graphs, the differentiation turns out to be the restricting factor. In Table 1 we state the computation times for a linear 3 node model and for the setup described before.

We would like to emphasize that the goal of our implementation is to illustrate the theoretical statements, and thus we did not focus on efficiency. The table shows that if one wants to adapt our implementation to work with larger graphs, it would be sufficient to implement a more efficient differentiation method. For instance, this could be done by using closed form expression of the derivatives, or by using automatic differentiation methods, or by using finite difference methods. We believe that such methods would probably be much faster.

<table><caption>Table 1. Refinement computation times for Random Access CSMA model for a 3 node linear graph and the 5 node graph as in Figures 1a, 1b. Times are given in seconds. In parentheses, we indicate the time taken by differentiation and subsequent sympy methods. These numbers show that the Sympy code takes more than 95% of the computation time.</caption><thead><tr><th></th><th>Jakobian (A)</th><th>Hessian (B)</th><th>v + w</th><th>s + t + u</th><th>Total</th><th>Sympy</th></tr></thead><tbody><tr><td>3 Node</td><td>6.08 (6.08)</td><td>30.71 (30.71)</td><td>0.19 (0)</td><td>5.64 (3.81)</td><td>42.63 (40.6)</td><td>95.22%</td></tr><tr><td>5 Node</td><td>97.09 (97.09)</td><td>897.15 (897.15)</td><td>1. (0)</td><td>122.36 (86.09)</td><td>1117.6 (1080.33)</td><td>96.66%</td></tr></tbody></table>

# 6 CONCLUSION

In this paper we investigate the accuracy of the classical *averaging method* that is used to study two timescale models. We study a generic two timescale model and show that under mild regularity conditions, the bias of this ‘average’ mean field approximation is of order $O(1/N)$. This result holds for any finite time-horizon, and extends to the steady-state regime under the classical assumption that the system has a unique and stable fixed point. Our results show the existence of a bias term $C_h(t)$ for any regular function $h$:

$$ \mathbb{E}[h(X_t)] = \underbrace{h(\Phi_t(x))}_{\text{classical 'average' mean field}} + \underbrace{\frac{1}{N}C_h(t)}_{\substack{O(1/N) \text{ expansion of the bias} \\ \text{refined 'average' mean field}}} + o(1/N). $$

For the steady-state regime $t = +\infty$, we propose an algorithmic method to calculate this term $C_h$. This correction term can be computed by solving linear systems and is therefore easily numerically computable. We show on an example that, similarly to what was done for classical one timescale models [15], the bias term leads to an approximation that is almost exact for small values of $N$ like $N = 10, 20$.

An interesting open question would be to obtain a characterization of $C_h(t)$ for the transient regime. Yet, it is not clear to us if those expressions would be usable as their size grows quickly with the system size. From an application point of view, our examples show that the new approximation leads to very accurate estimates for CSMA models. We believe that the same should hold for other multiscale models.

# ACKNOWLEDGMENTS

This work is supported by the French National Research Agency (ANR) through REFINO Project under Grant ANR-19-CE23-0015.
REFERENCES

[1] Karen Ball, Thomas G. Kurtz, Lea Popovic, and Greg Rempala. 2006. Asymptotic analysis of multiscale approximations to reaction networks. *The Annals of Applied Probability* 16, 4 (Nov. 2006), 1925–1961.
https://doi.org/10.1214/105051606000000420

[2] M. Benaim and J. Y. Boudec. 2008. A class of mean field interaction models for computer and communication systems. *Performance Evaluation* 65 (2008). https://doi.org/10.1016/j.peva.2008.03.005

[3] R. Boorstyn, A. Kershenbaum, B. Maglaris, and V. Sahin. 1987. Throughput Analysis in Multihop CSMA Packet Radio Networks. *IEEE Transactions on Communications* 35, 3 (1987), 267–274. https://doi.org/10.1109/TCOM.1987.1096769

[4] Charles Bordenave, David R. McDonald, and Alexandre Proutière. 2010. A particle system in interaction with a rapidly varying environment: Mean field limits and applications. *Networks and Heterogeneous Media* 5, 1 (2010), 31–62.

[5] Anton Braverman. 2022. The Prelimit Generator Comparison Approach of Stein's Method. *Stochastic Systems* 12, 2 (2022), 181–204. https://doi.org/10.1287/stsy.2021.0085

[6] Anton Braverman and J. G. Dai. 2017. Stein's method for steady-state diffusion approximations of $M/Ph/n + M$ systems. *The Annals of Applied Probability* 27, 1 (Feb 2017). https://doi.org/10.1214/16-AAP1211

[7] Anton Braverman, J. G. Dai, and Jiekun Feng. 2016. Stein's method for steady-state diffusion approximations: An introduction through the Erlang-A and Erlang-C models. *Stochastic Systems* 6, 2 (2016), 301–366.
https://doi.org/10.1214/15-SSY212

[8] Eyal Castiel, Sem Borst, Laurent Miclo, Florian Simatos, and Phil Whiting. 2021. Induced idleness leads to deterministic heavy traffic limits for queue-based random-access algorithms. *The Annals of Applied Probability* 31, 2 (2021), 941–971. https://doi.org/10.1214/20-AAP1609

[9] F. Cecchi. 2018. *Mean-field limits for ultra-dense random-access networks*. Technische Universiteit Eindhoven. Proefschrift.

[10] Fabio Cecchi, Sem C. Borst, Johan S. H. van Leeuwaarden, and Philip A. Whiting. 2021. Mean-Field Limits for Large-Scale Random-Access Networks. *Stochastic Systems* 11, 3 (2021), 193–217. https://doi.org/10.1287/stsy.2021.0068
arXiv:https://doi.org/10.1287/stsy.2021.0068

[11] Mathieu Feuillet and Philippe Robert. 2014. A Scaling Analysis of a Transient Stochastic Network. *Advances in Applied Probability* 46, 2 (Jun 2014), 516–535. https://doi.org/10.1239/aap/1401369705

[12] Nicolas Gast. 2017. Expected Values Estimated via Mean-Field Approximation are 1/N-Accurate. *Proceedings of the ACM on Measurement and Analysis of Computing Systems* 1, 1 (Jun 2017), 1–26. https://doi.org/10.1145/3084454

[13] Nicolas Gast, Luca Bortolussi, and Mirco Tribastone. 2019. Size expansions of mean field approximation: Transient and steady-state analysis. *Performance Evaluation* 129 (Feb. 2019), 60–80. https://doi.org/10.1016/j.peva.2018.09.005

[14] Nicolas Gast and Bruno Gaujal. 2012. Markov chains with discontinuous drifts have differential inclusion limits. *Performance Evaluation* 69, 12 (Dec 2012), 623–642. https://doi.org/10.1177/075910639203700105

[15] Nicolas Gast and Benny Van Houdt. 2017. A Refined Mean Field Approximation. *Proceedings of the ACM on Measurement and Analysis of Computing Systems* 1, 2 (Dec 2017), 33:1–33:28. https://doi.org/10.1145/3154491

[16] P.J. Hunt and T.G. Kurtz. 1994. Large loss networks. *Stochastic Processes and their Applications* 53, 2 (Oct 1994), 363–378.
https://doi.org/10.1016/0304-4149(94)90071-X

[17] Jeffrey J. Hunter. 1982. Generalized inverses and their application to applied probability problems. *Linear Algebra Appl.* 45 (1982), 157–198. https://doi.org/10.1016/0024-3795(82)90218-X

[18] Ilse C. F. Ipsen and Carl D. Meyer. 1994. Uniform Stability of Markov Chains. 15, 4 (1994), 1061–1074.
https://doi.org/10.1137/S0895479892237562

[19] Hye-Won Kang and Thomas G. Kurtz. 2013. Separation of time-scales and model reduction for stochastic reaction networks. *The Annals of Applied Probability* 23, 2 (Apr 2013). https://doi.org/10.1214/12-AAP841

[20] Hye-Won Kang, Thomas G. Kurtz, and Lea Popovic. 2014. Central limit theorems and diffusion approximations for multiscale Markov chain models. *The Annals of Applied Probability* 24, 2 (Apr 2014).
https://doi.org/10.1214/13-AAP934

[21] Thomas G. Kurtz. 1970. Solutions of ordinary differential equations as limits of pure jump markov processes. *Journal of Applied Probability* 7, 1 (April 1970), 49–58. https://doi.org/10.2307/3212147

[22] Thomas G. Kurtz. 1978. Strong approximation theorems for density dependent Markov chains. *Stochastic Processes and their Applications* 6, 3 (Feb 1978), 223–240. https://doi.org/10.1016/0304-4149(78)90020-0

[23] Jean-Yves Le Boudec, David McDonald, and Jochen Mundinger. 2007. A Generic Mean Field Convergence Result for Systems of Interacting Objects. In *Fourth International Conference on the Quantitative Evaluation of Systems (QEST 2007)*. IEEE, Edinburgh, Scotland, UK, 3–18. https://doi.org/10.1109/QEST.2007.8

[24] Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. 2017. SymPy: symbolic computing in Python. *PeerJ Computer Science* 3 (2017), e103.
[25] M. Mitzenmacher. 2001. The power of two choices in randomized load balancing. *IEEE Transactions on Parallel and Distributed Systems* 12, 10 (oct 2001), 1094–1104. https://doi.org/10.1109/71.963420

[26] Grigorios A. Pavliotis and Andrew M. Stuart. 2008. *Multiscale Methods: Averaging and Homogenization*. Texts Applied in Mathematics, Vol. 53. Springer New York, New York, NY.

[27] A. Pazy. 1983. *Semigroups of Linear Operators and Applications to Partial Differential Equations*. Applied Mathematical Sciences, Vol. 44. Springer New York, New York, NY. https://doi.org/10.1007/978-1-4612-5561-1

[28] Lawrence Perko. 2001. *Differential Equations and Dynamical Systems* (3rd ed.). Number 7 in Texts in Applied Mathematics. Springer-Verlag, Berlin, Heidelberg.

[29] Philippe Robert and Gaëtan Vignoud. 2021. Stochastic Models of Neural Synaptic Plasticity: A Scaling Approach. *SIAM J. Appl. Math.* 81, 6 (2021), 2362–2386. https://doi.org/10.1137/20M1382891

[30] Charles Stein. 1986. *Approximate Computation of Expectations*. Vol. 7. Institute of Mathematical Statistics. i–164 pages. http://www.jstor.org/stable/4355512

[31] P. M. van de Ven, S. C. Borst, J. S. H. van Leeuwaarden, and A. Proutière. 2010. Insensitivity and stability of random-access networks. *Performance Evaluation* 67, 11 (2010), 1230–1242. https://doi.org/10.1016/j.peva.2010.08.011

[32] Xin Wang and K. Kar. 2005. Throughput modelling and fairness issues in CSMA/CA based ad-hoc networks. In *Proceedings IEEE 24th Annual Joint Conference of the IEEE Computer and Communications Societies*. (2005-03), Vol. 1. 23–34 vol. 1. https://doi.org/10.1109/INFCOM.2005.1497875 ISSN: 0743-166X.

[33] Lei Ying. 2016. On the Approximation Error of Mean-Field Models. In *Proceedings of the 2016 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Science - SIGMETRICS '16*. ACM Press, Antibes Juan-les-Pins, France, 285–297. https://doi.org/10.1145/2896377.2901463

# A DEFINITIONS

In this section we revise some essential definitions and properties used in the paper. As these definitions are well established, we only briefly recall them to provide a self-contained paper.

## A.1 $C_0$-Semi-Group, Hölder Norm, ODE differentiability

**DEFINITION 8** ($C_0$ SEMI-GROUP [27] DEFINITION 2.1). $T_s$ is called strongly continuous semi-group (or $C_0$-semi-group) if

$$T_0 = \mathrm{Id}, \qquad T_{s+t} = T_s T_t \quad \text{for all } s, t \ge 0, \qquad \lim_{t \downarrow 0} T_t z = z \quad \text{for all } z. \qquad (22)$$

Since we will only work with $C_0$-semi-groups we will simply refer to them as semi-groups. **Semi-group & Generator commutation** The generator of a $C_0$-semi-group is defined by

$$Az = \lim_{t \to 0} \frac{1}{t}(T_t z - T_0 z).$$

A direct consequence of the definition of the generator and a standard property is that it commutes with its defining $C_0$-semi-group, i.e.,

$$AT_s h(z) = T_s Ah(z) = \frac{d}{ds} T_s h(z).$$

Note, this therefore holds true for the semi-groups given by the stochastic system $\Psi_s h(x, y) = \mathbb{E}[h(X_s, Y_s) | X_0, Y_0 = x, y]$ with generator $Lh(x, y)$ as in Equation (6) and the semi-group of the ODE $\Phi_s h(x) = h(\phi_s(x))$ with generator $\Lambda h(x) = D_x h(x) \bar{F}(x)$.

**DEFINITION 9** (HÖLDER NORM AND SPACE). For $U \subset \mathbb{R}^n$ and $u \in C^k(U)$

$$\|u\|_{k,\gamma} := \sum_k \sup_{x \in U} \left\| D^k u(x) \right\| + \sum_k \sup \left\{ \frac{\|D^k u(x) - D^k u(y)\|}{\|x-y\|^\gamma} \mid x, y \in U, x \neq y \right\}$$

is called Hölder norm. The space of functions for which the norm is finite is called Hölder space and denoted by $D^k_U$. For the case $\gamma = 1$ the Hölder space encloses all functions who are k-times continuously differentiable with bounded derivatives and who's k-th derivatives are Lipschitz continuous. The latter we simply denote by $D^k_U$.
An important implication is that all Hölder continuous functions are uniformly continuous.

LEMMA 10 (DRIFT INDUCED DIFFERENTIABILITY [28] THEOREM 1 p.80). Let E be an open subset of $\mathbb{R}^n$ containing $x_0$ and assume that $f \in C^k(E)$. Then there exists an $a > 0$ and $\delta > 0$ such that for all $y \in N_\delta(x_0)^3$ the initial value problem

$$ \dot{x} = F(x), \quad x(0) = y, $$

has a unique solution $u$ which is $k$-times continuously differentiable with respect to the initial condition for $t \in [-a, a]$.

PROOF. Theorem 1 p.80-83 and Remark 1 p.83 of [28]. □

# B TECHNICAL LEMMAS AND PROOFS

## B.1 Proof of Lemma 4

By assumption (A$_2$), the transition matrix $K(x)$ has a unique irreducible class. As pointed out in Section 2.3, the corresponding Markov chain has a unique stationary distribution that we denote by $\pi(x)$. Let $\Pi(x) := \mathbb{1}\pi^T(x)$ be the matrix whose lines are all equal to $\pi(x)$. By [17](Theorem 3.5, p.17), $(K(x) + \Pi(x))$ is non-singular and its inverse $(K(x) + \Pi(x))^{-1}$ is a generalized inverse to $K(x)$, which means that it satisfies $K(x)(K(x) + \Pi(x))^{-1}K(x) = K(x)$.

To obtain the solution to the Poisson equation (10) we only consider the case where $h$ takes values in $\mathbb{R}$. The extension to a function $h$ that takes values in $\mathbb{R}^n$ is straightforward as it corresponds to $n$ independent Poisson equations.

Let us suppress the dependence on $x$ for clarity. Recall that $K^+ = (K + \Pi)^{-1}(I - \Pi)$ and let us study the product $KK^+$:

$$
\begin{align*}
KK^+ &= K(K + \Pi)^{-1}(I - \Pi) \tag{by definition} \\
&= K(K + \Pi)^{-1}(I + K - K - \Pi) \\
&= K(K + \Pi)^{-1} + K(K + \Pi)^{-1}K - K(K + \Pi)^{-1}(K + \Pi) \tag{expanding the product} \\
&= K(K + \Pi)^{-1} \tag{the last two terms equal $\pm K$} \\
&= (K + \Pi)(K + \Pi)^{-1} - \Pi(K + \Pi)^{-1} \tag{Adding and subtracting $\Pi(K + \Pi)^{-1}$} \\
&= I - \Pi,
\end{align*}
$$

where the last equality holds because $\pi^T\mathbb{1} = 1$ and therefore $\Pi\Pi = \mathbb{1}\pi^T\mathbb{1}\pi^T = \mathbb{1}\pi^T = \Pi$.

Combined with $\Pi K = 0$, this shows that $\Pi = \Pi(K + \Pi)$ and therefore $\Pi(K + \Pi)^{-1} = \Pi$.

The above computations show that if $G_h^{\text{fast}}(x, y) = \sum_{y'} K_{y,y'}^{+}(x)h(x, y')$, then:

$$ K(x)G_h^{\text{fast}}(x, y) = h(x, y) - \sum_{y'} \pi_{y'}(x)h(x, y'), $$

which shows that $G_h^{\text{fast}}$ is the solution of the Poisson equation.

The differentiability of $G_h^{\text{fast}}$ follows from the differentiability of $h$ and $K^+$: Under Assumption (A$_1$) $K(x)$ is continuously differentiable. By Assumption (A$_2$), $K(x)$ has a unique irreducible class, this implies $\Pi(x)$ is continuously differentiable, for which we refer to [18], and therefore further implies that $(I + \Pi(x))^{-1}$ and $K^+(x)$ are continuously differentiable. This proves that $G_h^{\text{fast}}$ is continuously differentiable in $x$.

<footnote>$^3N_\delta(x_0) := \{x \in \mathbb{R}^n : \|x - x_0\| \le \delta\}$</footnote>
## B.2 Technical Lemmas used to prove Theorem 1 (Transient Regime)

LEMMA 11. For arbitrary but fixed $t > 0$, let $g: (s, x, y) \in [0, t] \times \mathcal{X} \times \mathcal{Y} \mapsto g_s(x, y) \in \mathbb{R}$ be a continuous function that is continuously differentiable in $s$ and let $H_s(x, y) := \mathbb{E}[g_s(X_s, Y_s) | X_0, Y_0 = x, y]$. Then:

$$H_t(x, y) - H_0(x, y) = \int_0^t \mathbb{E}[\frac{d}{dt}g_{s+\tau}(X_s, Y_s) |_{\tau=0}]ds + \int_0^t \mathbb{E}[Lg_s(X_s, Y_s)]ds.$$

PROOF. By definition of the generator $L$, the quantity $H_t(x, y) = \mathbb{E}[g_t(X_t, Y_t) | X_0, Y_0 = x, y]$ is right sided differentiable, i.e., $\frac{d^+}{ds}f(s) = \lim_{ds \to 0^+} \frac{1}{ds}(f(s+ds) - f(s))$, with respect to time. Using semi-group properties and bounds, its derivative is

$$\frac{d^+}{ds}H_s(x, y) = \mathbb{E}[\frac{d}{d\tau}g_{s+\tau}(X_s, Y_s) |_{\tau=0}] + \mathbb{E}[Lg_s(X_s, Y_s)].$$

The first term corresponds to the derivation of $g_s$ with respect to time and the second term to the changes of the stochastic system in $(X_s, Y_s)$. The lemma therefore follows by using that $H_t(x, y) - H_0(x, y) = \int_0^t \frac{d^+}{ds}H_s(X_s, Y_s) \quad \square$

LEMMA 12 (BOUND FOR $L_{\text{fast}}g(x, y)$). Take the two-timescale stochastic system as introduced in Section 2.1 with generator $L$ and assume $(A_1)$ and $(A_2)$. Let $L_{\text{fast}}$ be as defined in Section 4.2. Then, for $g \in \mathcal{D}^1(\mathcal{X} \times \mathcal{Y})$ and $(x, y) \in \mathcal{X} \times \mathcal{Y}$

$$L_{\text{fast}}g(x, y) - \frac{1}{N}Lg(x, y) = \frac{1}{N} \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \ell D_x g(x, y') + o(1/N).$$

PROOF. By definition of $L_{\text{fast}}$, for a continuous function $g$ the values of $L_{\text{fast}}g(x, y)$ and $\frac{1}{N}Lg(x, y)$ coincide in the limit. For finite $N$, we first look at $\frac{1}{N}Lg(x, y)$ which is given by

$$\frac{1}{N}Lg(x, y) = \frac{1}{N} \sum_{\ell, y'} N \alpha_{\ell, y'}(x, y) (g(x + \frac{\ell}{N}, y') - g(x, y)).$$

Using the continuity of $g$ in $x$, and definition of $L_{\text{fast}}g(x, y)$ we have

$$\begin{aligned} L_{\text{fast}}g(x, y) &= \lim_{N \to \infty} \frac{1}{N}Lg(x, y) = \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \left( \lim_{N \to \infty} g\left(x + \frac{\ell}{N}, y'\right) - g(x, y) \right) \\ &= \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \left( g(x, y') - g(x, y) \right). \end{aligned}$$

Using Taylor's theorem

$$\begin{aligned} Lg(x, y) &= N \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \left(g\left(x + \frac{\ell}{N}, y'\right) - g(x, y)\right) \\ &= N \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \left(g\left(x, y'\right) - g(x, y) + \ell D_x g(x, y') + o(\|\ell\|)\right) \\ &= NL_{\text{fast}}g(x, y) + \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) \ell D_x g(x, y') + o(1). \end{aligned}$$

□
## B.3 Technical Lemma used to prove Theorem 3 (steady-state refinement)

The next Lemma justifies that the first term of (20) is approximated in the limit by the 'average' version of $J_h$ as defined in Section 4.5. Second, Lemma 14 gives stability conditions for the solution of the 'slow' Poisson Equation.

**LEMMA 13.** Assume $(A_1)-(A_3)$ in particular assume $h \in D^1(X)$. Further assume that $G_F^{\text{fast}}$, solution to the 'fast' Poisson Equation, is cont. differentiable in $x$ and $G_h^{\text{slow}}$, solution to the 'slow' Poisson Equation, is twice cont. differentiable in $x$. Let $\bar{J}(\phi_\infty) = \sum_y \pi(\phi_\infty)J(\phi_\infty, y)$ be the 'average' version of $J_h(x, y) = \sum_{\ell, y'} \alpha_{\ell, y'}(x, y) D_x(D_x G_h^{\text{slow}}(x) G_F^{\text{fast}}(x, y'))\ell$ in $\phi_\infty$, then

$$ N\mathbb{E}[D_x G_h^{\text{slow}}(X_\infty)(\bar{F}(X_\infty) - F(X_\infty, Y_\infty))] = \bar{J}(\phi_\infty) + o(1). $$

**PROOF.** To prove the lemma, let $G_F^{\text{fast}}$ be the solution to the Poisson equation as given in ('fast' Poisson Equation). We use $G_F^{\text{fast}}$ to rewrite

$$ \begin{aligned} N\mathbb{E}[D_x G_h^{\text{slow}}(X_\infty)(\bar{F}(X_\infty) - F(X_\infty, Y_\infty))] &= -N\mathbb{E}[D_x G_h^{\text{slow}}(X_\infty)(F(X_\infty, Y_\infty) - \bar{F}(X_\infty))] \\ &= -N\mathbb{E}[D_x G_h^{\text{slow}}(X_\infty)(L_{\text{fast}} G_F^{\text{fast}}(X_\infty, Y_\infty))] = -N\mathbb{E}[L_{\text{fast}} D_x G_h^{\text{slow}}(X_\infty)(G_F^{\text{fast}}(X_\infty, Y_\infty))]. \end{aligned} $$

Adding $\mathbb{E}[LD_x G_h^{\text{slow}}(X_\infty)(G_h^{\text{fast}}(X_\infty, Y_\infty))] = 0$ and applying the steps as in the proof of Lemma 12 we see that

$$ \begin{aligned} & -N\mathbb{E}[L_{\text{fast}} D_x G_h^{\text{slow}}(X_\infty)(G_F^{\text{fast}}(X_\infty, Y_\infty)) - LD_x G_h^{\text{slow}}(X_\infty)(G_F^{\text{fast}}(X_\infty, Y_\infty))] \\ &= -\mathbb{E}\left[\sum_{\ell, y'} \alpha_{\ell, y'}(X_\infty, Y_\infty)\left(D_x(D_x G_h^{\text{slow}}(X_\infty)G_F^{\text{fast}}(X_\infty, y'))\right)\ell\right] + o(1). \end{aligned} \quad (23) $$

The last equality follows directly from the definition of $L_{\text{fast}}$ since $G_h$ only depends on $x$. Using Theorem 2 we see that $\mathbb{E}[J_h(X_\infty, Y_\infty)]$ with

$J_h(x, y) := -\sum_{\ell, y'} \alpha_{\ell, y'}(x, y)\ell D_x(D_x G_h^{\text{slow}}(x)G_F^{\text{fast}}(x, y'))$ of equation (23) is approximated by $\bar{J}(\phi_\infty) = \sum_y \pi(\phi_\infty)J(\phi_\infty, y)$. This concludes the proof as it implies

$$ N\mathbb{E}[D_x G_h^{\text{slow}}(X_\infty)(\bar{F}(X_\infty) - F(X_\infty, Y_\infty))] = \mathbb{E}[J_h(X_\infty, Y_\infty)] + o(1) = \bar{J}_h(\phi_\infty) + o(1). \quad \square $$

## B.4 Stability of $G_h^{\text{slow}}$

**LEMMA 14 (STABILITY).** Assume that $\bar{F}$ and $\phi$ are k-times differentiable with uniformly continuous derivatives and that $\phi$ has a unique exponentially stable attractor $\phi_\infty$. Then the k-th derivative of $G_h^{\text{slow}}: x \mapsto \int_0^\infty h(\phi_s x) - h(\phi_\infty)ds$ is bounded and equal to $\int_0^\infty D_x^k(h \circ \phi_s)(x)ds$.

**PROOF.** This is a consequence of [15][Lemma 3.5].

## B.5 Proof of Proposition 7

**PROOF.** In the first part of the proof we find computable expressions of $D_x G_h^{\text{slow}}$, $D_x^2 G_h^{\text{slow}}$ and $D_x G_F^{\text{fast}}$. These expressions allow us to rewrite $J_h$ and construct the closed form representation of $J_h$ as well as $\bar{J}_h$ in the following steps. By definition

$$ \begin{aligned} J_h(x, y) &= -\sum_{\ell, y'} \alpha_{\ell, y'}(x, y) D_x(D_x G_h^{\text{slow}}(x) G_F^{\text{fast}}(x, y'))\ell \\ &= -\sum_{\ell, y'} \alpha_{\ell, y'}(x, y) (G_F^{\text{fast}}(x, y')^T D_x^2 G_h^{\text{slow}}(x) + D_x G_h^{\text{slow}}(x) D_x G_F^{\text{fast}}(x, y'))\ell. \end{aligned} \quad (24) $$
By Lemma 4 there exists a matrix $K^+(x)$ such that $G_F^{\text{fast}}(x, y')$ has the form

$$G_F^{\text{fast}}(x, y') = \sum_{y' \in \mathcal{Y}} K_{y', y''}^{+}(x) F(x, y'').$$

Assumption (A₁) which assures differentiability of the transition rates $\alpha$ with respect to $x$, also implies differentiability for $G_F^{\text{fast}}$. Therefore,

$$D_x G_F^{\text{fast}}(x, y') = \sum_{y'' \in \mathcal{Y}} F(x, y'') \nabla_x^T K_{y', y''}^{+}(x) + K_{y', y''}^{+}(x) D_x F(x, y'') \text{ with}$$

$\nabla_x^T = [\frac{\partial}{\partial x_1}, \dots, \frac{\partial}{\partial x_n}]$, $D_x f = (\frac{\partial}{\partial x_j} f_i)_{i,j=1\dots n}$. Using the results of [15][Lemma 3.6] with A, B the first and second derivative of $\bar{F}$ as defined in (21), it holds that for the equilibrium point $\phi_\infty$

$$D_x G_h^{\text{slow}}(\phi_\infty)_i = \sum_j \frac{\partial h}{\partial x_j}(\phi_\infty) \int_0^\infty (D_x \phi_s(\phi_\infty))_{j,i} ds = \sum_j \frac{\partial h}{\partial x_j}(\phi_\infty) (-A)^{-1}_{j,i}$$

as well as

$$\begin{align}
D_x^2 G_h^{\text{slow}}(\phi_\infty)_{n,m} &= \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j}(\phi_\infty) \int_0^\infty (D_x \phi_s(\phi_\infty))_{j,n} (D_x \phi_s(\phi_\infty))_{i,m} ds \tag{25} \\
&\quad + \sum_i \frac{\partial h}{\partial x_i}(\phi_\infty) \int_0^\infty (D_x^2 \phi_s(\phi_\infty))_{j,m,n} ds \nonumber \\
&= \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j}(\phi_\infty) \int_0^\infty (e^{As})_{j,n} (e^{As})_{i,m} ds \nonumber \\
&\quad + \sum_i \frac{\partial h}{\partial x_i}(\phi_\infty) \left( \sum_j (-A)^{-1}_{i,j} \sum_{k_1,k_2} B_{j,k_1,k_2} \int_0^\infty (e^{As})_{k_1,n} (e^{As})_{k_2,m} ds \right). \nonumber
\end{align}$$

Next, the above equations are used to rewrite (24). To obtain a closed form expression for the left summand, we start by looking at the sum

$$\sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) \sum_{y'} K_{y', y'}^{+} F_m(\phi_\infty, y') \ell_n \int_0^\infty (e^{As})_{k_1, n} (e^{As})_{k_2, m} \tag{26}$$

A solution to the above equation is given by the following Lyapunov equation. To ease notations, define

$$O_{m,n} := \sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) K_{y', :}^{+}(\phi_\infty) F_m(\phi_\infty, :), \ell_n,$$

or, equivalently in matrix notation, $O = \sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) K_{y', :}^{+}(\phi_\infty) F(\phi_\infty, :), \ell^T$. If a matrix $U$ solves<fnref n="4" /> the Sylvester equation (for $X$)

$$AX + XA^T = -O, \tag{27}$$

it is equal to (26). Applying this identity and (25) to the first summand of $\bar{J}(\phi_\infty)$ which is the 'average' version of (24), lets us rewrite

$$\begin{align*}
& \sum_y \pi_y(\phi_\infty) \sum_{\ell, y'} \alpha_{\ell, y'}(\phi_\infty, y) G_F^{\text{fast}}(\phi_\infty, y')^T D_x^2 G_h^{\text{slow}}(\phi_\infty) \ell \\
& = \sum_{i,j} \frac{\partial^2 h}{\partial x_i \partial x_j}(\phi_\infty) U_{i,j} + \sum_i \frac{\partial h}{\partial x_i}(\phi_\infty) \sum_j (-A)^{-1}_{i,j} \sum_{k_1,k_2} B_{j,k_1,k_2} U_{k_1,k_2}.
\end{align*}$$

<footnote n="4">⁴As A is non-singular, A and $-A^T$ don't share any eigenvalues and therefore equation (27) has a unique solution.</footnote>
As shows in Lemma 4 the solution to the 'fast' Poisson Equation 10 has the form
$G_h^{\text{fast}}(x, y) = \sum_{y'} K_{y,y'}^+(x)h(x, y')$. To compute the bias correction terms, it is necessary to calculate the first derivative of $G_h^{\text{fast}}$ with respect to $x$. Since the computation of the derivative of $(K(x) + \Pi(x))^{-1}$ can be non-trivial and time consuming, the note below elaborates how the derivative can be efficiently obtained.

NOTE 16 (COMPUTATION OF $D_x K^+(x)$). By definition $K^+(x) = (K(x) + \Pi(x))^{-1}(I - \Pi(x))$. Using basic matrix derivation rules one has

$$
\begin{align}
\frac{\partial}{\partial x_i} K^+(x) &= \frac{\partial}{\partial x_i} \left( (K(x) + \Pi(x))^{-1} (I - \Pi(x)) \right) \nonumber \\
&= \frac{\partial}{\partial x_i} \left( K(x) + \Pi(x) \right)^{-1} (I - \Pi(x)) - \left( K(x) + \Pi(x) \right)^{-1} \frac{\partial}{\partial x_i} \Pi(x). \tag{29}
\end{align}
$$

Since the numerical difficulty lies only the computation of $D_x(K(x) + \Pi(x))^{-1}$ we focus solely on it.
Define $E(x) = (K(x) + \Pi(x))$ and let $I$ be the identity matrix. As pointed out in Lemma 4
$(K(x) + \Pi(x))$ is indeed invertible and thus for the partials derivative $\frac{\partial}{\partial x_i}$, $i = 0...d_x$ the following
holds:

$$
\frac{\partial}{\partial x_i} I = \frac{\partial}{\partial x_i} (E(x) E^{-1}(x))
$$

$$
\Leftrightarrow \quad 0 = \left( \frac{\partial}{\partial x_i} E(x) \right) E^{-1}(x) + E(x) \frac{\partial}{\partial x_i} E^{-1}(x)
$$

$$
\Leftrightarrow \quad \frac{\partial}{\partial x_i} E^{-1}(x) = -E^{-1}(x) \left( \frac{\partial}{\partial x_i} E(x) \right) E^{-1}(x)
$$

$$
\Leftrightarrow \quad \frac{\partial}{\partial x_i} (K(x) + \Pi(x))^{-1} = -(K(x) + \Pi(x))^{-1} \left( \frac{\partial}{\partial x_i} (K(x) + \Pi(x)) \right) (K(x) + \Pi(x))^{-1}.
$$

The above can now be used to compute Equation (29).

D NUMERICAL RESULTS FOR THE 3 NODE MODEL

For completeness, we give the numerical results of the linear 3 node model of Graph 1a. To obtain
the results we modify the parameters of the code of Code Cell 1 to match with the 3 node setup.
The new model is defined as in Code Cell 2.

Code Cell 2. Initialization of Mean Field for the 3 Node Model.

# Graph structure (this is the 3 node example)

G = np.array([[0,1,0],
              [1,0,1],
              [0,1,0]])

# rates & buffer size

_lambda = np.array([0.4, 0.2, 0.5])

nu = np.array([1.2, 2., 1.5])

mu = np.array([1.4, 1.3, 1.7])

buffer_size = 10

As seen in Figure 4, we obtain similar results as for the 5 node model when considering
steady-state average queue length values. The sample mean and confidence interval are computed
from 40 steady-state samples in the same manner as for the 5 node example. The refined
approximation almost exactly indicates the stationary value of the stochastic process even for
small N while the mean field approximation gets more accurate as N grows.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Stationary Queue Length Distribution for the 3 Node Graph of Figure (1a).</figcaption></figure>