Nicolas Gast
email: nicolas.gast@inria.fr

Luca Bortolussi
email: lbortolussi@units.it

Mirco Tribastone
email: mirco.tribastone@imtlucca.it

Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis

.

In this paper, we improve this result in two directions. First, we show how to obtain the same result for the transient regime. Second, we provide a further refinement by expanding the term in 1/N 2 (both for transient and steady-state regime). Our derivations are inspired by moment-closure approximation, a popular technique in theoretical biochemistry. We provide a number of examples that show: (1) that this new approximation is usable in practice for systems with up to a few tens of dimensions, and (2) that it accurately captures the transient and steady state behavior of such systems.

INTRODUCTION

Mean field approximation is a widely used technique in the performance evaluation community. The focus of this approximation is to study the performance of systems composed of a large number of interacting objects. Applications range from biological models [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF] to epidemic spreading [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF] and computer-based systems [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF]. In the performance evaluation community, this approximation has successfully been used to characterize the performance of CSMA protocols, [START_REF] Cecchi | Mean-field analysis of ultra-dense csma networks[END_REF], information spreading algorithms and peer-to-peer networks [START_REF] Chaintreau | The Age of Gossip: Spatial Mean Field Regime[END_REF][START_REF] Massoulié | Coupon Replication Systems[END_REF], caching [START_REF] Chen | A projection method and Kronecker product preconditioner for solving Sylvester tensor equations[END_REF][START_REF] Fagin | Asymptotic miss ratios over independent references[END_REF][START_REF] Gast | Transient and steady-state regime of a family of list-based cache replacement algorithms[END_REF] or a quite popular subject such as load balancing strategies [18, 32, 34-36, 43, 45, 48]. This approximation can be used to study transient (for example the time to fill a cache [START_REF] Gast | Transient and steady-state regime of a family of list-based cache replacement algorithms[END_REF]) or steady-state properties (for example the steady-state hit ratio [START_REF] Chen | A projection method and Kronecker product preconditioner for solving Sylvester tensor equations[END_REF][START_REF] Fagin | Asymptotic miss ratios over independent references[END_REF]).

One of the reasons of the success of mean field approximation is that it is often very accurate as soon as N , the number of objects in the system, exceeds a few hundreds. In fact, this approximation can be proven to be asymptotically exact as N goes to infinity, see for example [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF][START_REF] Thomas | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF][START_REF] Thomas | Strong approximation theorems for density dependent Markov chains[END_REF] and explicit bounds for the convergence rate exist [START_REF] Bortolussi | Bounds on the deviation of discretetime Markov chains from their mean-field model[END_REF][START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF][START_REF] Ying | Stein's Method for Mean Field Approximations in Light and Heavy Traffic Regimes[END_REF].

Recently, the authors of [START_REF] Gast | A Refined Mean Field Approximation[END_REF] proposed what they call a refined mean field approximation that can be used to characterize more precisely steady-state performance metrics. Their refinement uses that for many models, a steady-state expected performance metric of a system with N objects E[h(X )] is equal to its mean field approximation h(π ) plus a term in 1/N :

E [h(X )] = h(π ) + 1 N V (h) + o 1 N , (1)
where π is the fixed point of the ODE that describes the mean field approximation and V (h) is a constant that can be easily evaluated numerically.

By using a number of examples, they show that the refined approximation h(π ) + 1 N V (h) is much more accurate than the mean field approximation for moderate system sizes (i.e., a few tens of objects).

In this paper, we extend this method in two directions: First we generalize Equation (1) to the transient behavior; second we establish the existence of a second order term in 1/N 2 (both in transient and steady-state regimes). More precisely, we establish conditions such that for any smooth function h, there exist constants V (h) and A (h) such that for any time t ∈ [0; ∞) ∪ {∞} :

E [h(X (t))] = h(x(t)) + 1 N V (h) (t) + 1 N 2 A (h) (t) + o 1 N 2 . ( 2
)
We show that for the transient regime, V (h) (t) and A (h) (t) satisfy a linear time-inhomogeneous differential equation that can be easily integrated numerically (Theorem 1). The steady-state constants are directly computed from the fixed point of this linear differential equation (Theorem 3). We use Equation [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF] to propose two new approximations that depend on the system size N and that are expansions of the classical mean field approximation to the order 1/N and 1/N 2 , respectively. We then compare the following three approximations numerically on various examples :

• Mean field approximation: h(x(t)).

• 1/N -expansion: h(x(t)) + V (h) (t)/N ; • 1/N 2 -expansion: h(x(t)) + V (h) (t)/N + A (h) (t)/N 2 .

Our numerical results shows that the two expansions capture very accurately the transient behavior of such a system even when N ≈ 10. Moreover, they are generally much more accurate than the classical mean field approximation for small values of N (for transient and steady-state regimes). Our experiments also confirm that good accuracy of the 1/N -expansion approximation that was observed for steady-state values in [START_REF] Gast | A Refined Mean Field Approximation[END_REF] : In most cases, the largest gain in accuracy comes from the 1/N -term (both for the transient and steady-state values). The 1/N 2 -term does improve the accuracy but only marginally. We also study the limit of the method by studying an unstable mean field model that has an unstable fixed point. This last example has unique fixed point that is not an attractor which means that the classical mean field approximation cannot be use for steady-state approximation as shown in [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF]. We show that in this case, the 1/N and 1/N 2 expansions are not stable with time and are therefore inaccurate when the time becomes large.

To summarize, this paper makes theoretical contributions that are interesting from a practical perspective :

• Theoretical contributions -We show that the 1/N -expansion proposed in [START_REF] Gast | A Refined Mean Field Approximation[END_REF] for steady-state estimation can be extended to the transient regime and can be refined to the next order correction term in 1/N 2 . • Practical implications -We show that, despite the complexity of the formulas, it is relatively easy to compute the 1/N and 1/N 2 terms (in the transient and steady-state regimes) for realistic models such as the supermarket model. The developed method is generic and is implemented in a tool [START_REF] Gast | Refined Mean Field Tool[END_REF].

Roadmap. The rest of the paper is organized as follows. We discuss related work in Section 2. We describe the model in Section 3. We develop the main results in Section 4 where we also provide the proofs. We show a simple malware propagation model in Section 5 in order to illustrate the main concepts. We then study the supermarket model in more detail in Section 6. In Section 7 we show an example that illustrates the limitations of the approach. Finally, we conclude in Section 8.

Reproducibility. The code to reproduce the paper -including simulations, figures and text -is available at https://github.com/ ngast/sizeExpansionMeanField [START_REF] Gast | Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis[END_REF].

RELATED WORK

Our results apply to the classical density-dependent population process of Kurtz [START_REF] Thomas | Strong approximation theorems for density dependent Markov chains[END_REF] of which the supermarket model of [START_REF] David | The Power of Two Random Choices in Randomized Load Balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF] is an example.

Stein's Method

From a methodological point of view, our paper uses an approach similar to one of [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF][START_REF] Gast | A Refined Mean Field Approximation[END_REF][START_REF] Vassili N Kolokoltsov | Mean field games and nonlinear Markov processes[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF][START_REF] Ying | Stein's Method for Mean Field Approximations in Light and Heavy Traffic Regimes[END_REF] in which the key idea is to compare an asymptotic expansion of the generator of the stochastic process with the generator of the mean field approximation, by using ideas inspired by Stein's method. In the papers [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF][START_REF] Vassili N Kolokoltsov | Mean field games and nonlinear Markov processes[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF][START_REF] Ying | Stein's Method for Mean Field Approximations in Light and Heavy Traffic Regimes[END_REF], this is used to obtain the rate of convergence of mean field models to their limit. In [START_REF] Gast | A Refined Mean Field Approximation[END_REF], this idea is used to compute the 1/N -term for the steady-state behavior. The main theoretical contribution of the paper with respect to these is to show that this method can be pushed further to study transient regime and to obtain exact formula for the term in 1/N 2 . The work on Stein's method is not new [START_REF] Stein | Approximate computation of expectations[END_REF] but has seen a regain of interest in the stochastic networks' community in the recent years thanks to the work of [START_REF] Braverman | Stein's method for steady-state diffusion approximations of M /Ph/n + M systems[END_REF][START_REF] Braverman | Stein's method for steady-state diffusion approximations: an introduction through the Erlang-A and Erlang-C models[END_REF].

System Size Expansion

Our paper is also closely related to an approach developed in the theoretical biology literature, known as system size expansion (SSE). The core idea of SSE dates back to the work of Van Kampen [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF], and consists in working with the stochastic process expressing the fluctuations of the population model around the mean field limit, rescaled by N -1/2 , and approximating it by an absolute continuous process ξ (t) taking real values. Starting from the Kolmogorov equation of the population model, and relying on a perturbation expansion, Van Kampen obtains an Fokker-Plank (FP) equation for ξ (t) containing in the right hand side terms of order N -p/2 , for p = 0, 1, . . .. Keeping only lower order terms (i.e. of order 0 and -1/2) results in a linear FP equation, whose solution is known as the linear noise approximation, which is equivalent to the central limit theorem proved by Kurtz [START_REF] Stewart | Markov processes: characterization and convergence[END_REF].

Grima and coauthors, in [START_REF] Grima | An effective rate equation approach to reaction kinetics in small volumes: Theory and application to biochemical reactions in nonequilibrium steady-state conditions[END_REF] and following papers (see e.g. [START_REF] Grima | How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations[END_REF][START_REF] Thomas | Intrinsic Noise Analyzer: A Software Package for the Exploration of Stochastic Biochemical Kinetics Using the System Size Expansion[END_REF][START_REF] Thomas | Computation of biochemical pathway ?uctuations beyond the linear noise approximation using iNA[END_REF]), start from the FP and keep higher order terms of

1 √ N
, introducing non-linear corrections to the linear approximation. The resulting FP cannot be solved exactly, but it can be used to derive differential equations for the mean, covariance, and potentially higher order moments. As far as the mean of the populations is concerned, the equation derived in [START_REF] Grima | An effective rate equation approach to reaction kinetics in small volumes: Theory and application to biochemical reactions in nonequilibrium steady-state conditions[END_REF][START_REF] Grima | How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations[END_REF] shows an equivalent structure with the one obtained in this paper. The higher-order SSE equations, with corrections up to order N -2 , have been implemented in the tool iNA [START_REF] Thomas | Intrinsic Noise Analyzer: A Software Package for the Exploration of Stochastic Biochemical Kinetics Using the System Size Expansion[END_REF][START_REF] Thomas | Computation of biochemical pathway ?uctuations beyond the linear noise approximation using iNA[END_REF], and more recently in the Matlab toolbox CERENA [START_REF] Kazeroonian | CERENA: ChEmical REaction Network Analyzer?A Toolbox for the Simulation and Analysis of Stochastic Chemical Kinetics[END_REF], the only working implementation to the authors' knowledge.

Even if equations for the mean population and for covariance of SSE and our method coincide, our approach has some advantanges. First of all, its derivation is rigorous and does not rely on any approximation of the process ξ (t), being based on a perturbation expansion of the moment equations themselves. Secondly, it gives us an approximate equation for any function h of the population vector, which can be used to estimate higher order moments or hitting times. Finally, in this paper we validate our method with large-dimensional models : the 1/N -expansion can be computed for models with hundreds of dimensions and the 1/N 2 -expansion can be computed for models with a few tens of dimensions.

Moment-closure Approximation

Our way of deriving the equations is also related to moment closure techniques [START_REF] Gillespie | Moment-closure approximations for mass-action models[END_REF], which work by truncating, at a finite order of moments, the exact infinite dimensional system of ODEs which captures the evolution in time of all moments of the population process. The truncation strategy typically assumes some form of the distribution, and uses the relationship among moments implied by that assumption to express high-order moments as a function of lower order ones (e.g. a Gaussian distribution has odd centered moments of order 3 and more all equal to zero). These techniques are in theory applicable to higher order moment -see for example [START_REF] Ale | A general moment expansion method for stochastic kinetic models[END_REF] -but the approach presented in [START_REF] Ale | A general moment expansion method for stochastic kinetic models[END_REF] seems difficult to apply in high dimensional models, due to the exponential dependence on the order of moments of the number of moment equations. The accuracy of moment closure approximations was studied in [START_REF] Grima | A study of the accuracy of moment-closure approximations for stochastic chemical kinetics[END_REF], and more recently in [START_REF] Schnoerr | Validity conditions for moment closure approximations in stochastic chemical kinetics[END_REF][START_REF] Schnoerr | Comparison of different moment-closure approximations for stochastic chemical kinetics[END_REF]. These studies show that accuracy is subtle and hard to predict, and does not necessarily increase with the population size N . The method we present in this paper uses a more rigorous approach, rooted in convergence theorems, which guarantees exactness in the limit of large N , and can also be used to provide estimates of moments of any order without extra effort, by choosing proper functions h.

MODEL AND NOTATIONS 3.1 Density-Dependent Population Processes

We consider mean field models described by the classical model of density-dependent population process of [START_REF] Thomas | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF]. A density dependent population process is a sequence of continuous time Markov chains X (N ) , where the index N is called the size of the system. For each N , the Markov chain X (N ) evolves on a subset E ⊂ R d , where d is called the dimension of the model. We assume that there exists a set of vectors L ∈ E and a set of functions β ℓ : E → R + such that X (N ) jumps from x to x + ℓ/N at rate N β ℓ (x) for each ℓ ∈ L.

Note that we state all our results using the framework of densitydependent population processes. An alternative would have been to used a continuous-time version of the discrete-time model of [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF] for which our results can be adapted (see also the discussion in Section 2.3 of [START_REF] Gast | A Refined Mean Field Approximation[END_REF]).

Drift and Mean Field approximation

We define the drift f as

f (x) = ℓ ∈ L ℓβ ℓ (x)
The drift is the expected variation of X (N ) (t) when X (N ) (t) = x. By definition of the model, it is independent of N .

In all our results, we will assume that the ordinary differential equation (ODE) x = f (x) has a unique solution that starts in x(0) at time 0 that we denote t → Φ t x. It satisfies: Φ t x = x + ∫ t 0 Φ s xds. When it is not ambiguous, we will denote x(t) := Φ t x. The function t → x(t) is called the mean field approximation.

Tensors, Derivatives and Einstein Notations

Our results rely on tensor computation. To simplify the expression of the results and their derivations, we use Einstein notation (also known as Einstein summation convention) that we recall here.

All vectors (or tensors) are d-dimensional

(or of size d × d, d × • • • × d).
For a given vector or tensor, the upper indices denote the component. For example, X i denotes the ith component of a d-dimensional vector X , and C i jk denotes the (i, j, k) components of a d × d × d-dimensional tensor C. We use the symbol ⊗ for the Kronecker product between two tensors: for two

d-dimensional vectors X and Y , X ⊗ Y denotes a d × d-dimensional tensor whose component (i, j) is X i Y j . Also, Y ⊗3 = Y ⊗ Y ⊗ Y .
For a given function, the lower indices denote the variable on which we differentiate. Unless otherwise stated, the functions will always be evaluated at the mean field approximation x(t). We use uppercase letters to denote the function evaluated at x(t). To be more precise, this means that the quantity F i j 1 ...j k denotes the kth derivative of the ith component of f with respect to x j 1 . . . x j k evaluated at x(t):

F i j 1 ...j k = ∂ k f i ∂x j 1 . . . ∂x j k (x(t))
We use Einstein summation convention, which implies summation over a set of repeated indices: each index variable that appears twice implies the summation over all the values of the index. For example F i j V j := j F i j V j and F i j,k, ℓ B k, ℓ := k, ℓ F i j,k, ℓ B k, ℓ . This convention greatly compactifies and therefore simplifies the expression of our results.

For a given d ×k tensor T , we denote by Sym(T ) the symmetric part of a tensor, which is the summation of this tensor over all permutation of indices. Its (i 1 . . . i k )-component is:

Sym(T ) i 1 ...i k = 1 k! σ ∈S k T i σ 1 ...i σ k ,
where S k is the symmetric group on k elements.

Summary of the Assumptions

In order the prove our results for the transient regime, we will use the following assumptions.

(A1) The sequence of stochastic processes X (N ) is a density dependent process that evolves in a compact subset of E ⊂ R d . (A2) The drift function f (x) is well defined and continuously differentiable four times. The function q(x) = ℓ ∈L ℓ ⊗ ℓβ ℓ (x) is well defined and continuously differentiable twice.

The function r (x) = ℓ ∈L ℓ ⊗ ℓ ⊗ ℓβ ℓ (x) is well defined and continuous.

Note that assumption (A2) on the differentiability of the drift, combined with assumption (A1) on the compactness of E implies that the drift is Lipschitz-continuous and bounded and that therefore the differential equation x = f (x) has a unique solution. These assumptions are mainly technical and are verified by many of the mean field models of the literature.

For the steady-state analysis, we will assume in addition:

(A3) For each N , the stochastic process X (N ) has a unique stationary distribution. (A4) The differential equation x = f (x) has a unique fixed point π that is a globally exponentially stable attractor, meaning that there exists two constants a, b > 0 such that for all x ∈ E:

∥Φ t (x) -π ∥ ≤ ae -bt .
Assumption (A3) combined with the existence of a globally stable attractor is a natural condition when one wants to show that a stochastic model converges to the fixed point of its mean field approximation (this is often a necessary condition, as shown in [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802.11 MAC protocol[END_REF]). The exponential stability of this attractor is a natural condition to obtain rate of convergence for mean field models [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF]. Proving that a fixed point is an attractor is often difficult but showing that this attractor is exponentially stable is often much easier since it only depends on the eigenvalue properties of the Jacobian evaluated at the fixed point π .

MAIN RESULTS

In this section, we provide the main theoretical results. We start by stating the results for the transient case ( §4.1), and the steady-state case ( §4.2). We then comment on the numerical feasibility of the approach ( §4.3) and we finish with the proofs ( §4.4).

Transient Analysis

The main result of our analysis is Theorem 1, which characterizes how the moments of the difference between the stochastic system X (t) and its mean field approximation evolve with time. We show that each of these moments admits an expansion with a first term in 1/N and a second term in 1/N 2 . The constants of this asymptotic expansion are characterized by a system of linear ODEs. One of the direct consequence of this theorem is Corollary 2 that provides an asymptotic expansion of the mean and the variance of X (N ) .

Theorem 1. Under assumption (A1-A2), let x(t) denote the unique solution of the ODE x = f (x) starting in X N (0). There exists a series of time-dependent tensors V ,W , A, B, C and D such that, for any four-time differentiable function h : R d → R, we have:

E h(X (N ) (t)) = h(x(t)) + 1 N H i V i + 1 2 H i j W i j (3) + 1 N 2 H i A i + 1 2 H i j B i j + 1 6 H i jk C i jk + 1 24 H i jk ℓ D i jkl + o 1 N 2 ,
where the terms H i . . . H i jk ℓ denotes the first to fourth derivative of h evaluated at x(t).

The dimension of the tensors V and A is n; the dimension of W and B is n × n; the dimension of C is n × n × n; the dimension of D is n × n × n × n. For the 1/N -terms, these tensors satisfy the following ODE system (with the initial conditions V = 0 and W = 0):

V i = F i j V j + 1 2 F i j,k W j,k W i, j = F i k W k, j + F j k W k,i + Q i, j = Sym 2F i k W k j + Q i j
For the 1/N 2 -terms, the ODE system is as follows (with the initial conditions A = 0,B = 0, C = 0 and D = 0)

A i = F i j A j + 1 2 F i j,k B j,k + 1 6 F i j,k, ℓ C j,k, ℓ + 1 24 F i j,k, ℓ,m D j,k, ℓ,m B i j = Sym 2F i k B k j + F i k ℓ C k ℓj + 1 3 F i k ℓm D k ℓmj + Q i j k V k + 1 2 Q i j k ℓ W k ℓ C i jk = Sym 3F i ℓ C ℓjk + 3 2 F i ℓm D ℓmjk + 3Q i j V k + 3Q i j ℓ W ℓk + R i jk D i jk ℓ = Sym 4F i m D mjk ℓ + 6Q i j W k ℓ .
where the symmetric d × d tensor Q and d ×3 tensor R are:

Q = ℓ ∈ L (ℓ ⊗ ℓ)β ℓ (x(t)) (4)
R = ℓ ∈ L (ℓ ⊗ ℓ ⊗ ℓ)β ℓ (x(t)); (5)
The tensors Q k and Q k, ℓ correspond to the first and second derivatives of the function x → ℓ ∈L ℓ ⊗ ℓβ ℓ (x), evaluated in x(t):

Q k = ∂ ∂x k ℓ ∈L (ℓ ⊗ ℓ)β ℓ (x(t)) Q k, ℓ = ∂ 2 ∂x k ∂x ℓ ℓ ∈L (ℓ ⊗ ℓ)β ℓ (x(t)).
To prove this theorem, we will first prove the existence of the tensors and then will show that they satisfy the corresponding set of ODEs by computing how the moments evolve with time. In fact, an equivalent characterization of the tensors V , W ,... is to use these tensors to construct asymptotic expansions of the moments of X (N ) (t) -x(t). This is summarized in Corollary 2, which also has an interest in its own. This corollary also justifies why moment closure works: neglecting the first moment of X (N ) (t) -x(t) gives the mean field approximation, neglecting the moment three and above gives the expansion of order 1/N ; finally neglecting the moments five and above gives the expansion of order 1/N 2 . In theory, it should be possible to continue the asymptotic expansion but the at the price of a much higher complexity in the expressions. In the numerical examples, we will show that the asymptotic development of the expectation provides a very accurate estimation of the true expectation in many cases.

Corollary 2. Under the assumption of Theorem 1, we have

E X (N ) (t) -x(t) = 1 N V (t) + 1 N 2 A(t) + o(1/N 2 ) E (X (N ) (t) -x(t)) ⊗2 = 1 N W (t) + 1 N 2 B(t) + o(1/N 2 ) E (X (N ) (t) -x(t)) ⊗3 = 1 N 2 C(t) + o(1/N 2 ) E (X (N ) (t) -x(t)) ⊗4 = 1 N 2 D(t) + o(1/N 2 ) E (X (N ) (t) -x(t)) ⊗k = o(1/N 2 ) for k ≥ 5.
In particular:

cov(X (N ) (t), X (N ) (t)) = 1 N W (t) + 1 N 2 (B(t) -V (t) ⊗ V (t)) + o(1/N 2 ).
Proof. The first set of equation is a direct consequence of Theorem 1 applied to the functions h(X ) = (X -x) ⊗k for k = 1, 2 . . . . For the covariance, we have :

cov(X (N ) (t), X (N ) (t)) = E (X (N ) (t) -x(t) + x(t) -E X (N ) (t) ) ⊗2 = E (X (N ) (t) -x(t)) ⊗2 -(x(t) -E X (N ) (t) ) ⊗2 = 1 N W (t) + 1 N 2 (B(t) -V (t) ⊗ V (t)) + o(1/N 2 ). □

Steady-State Regime

We now turn our attention to the steady-state regime. The next theorem shows that when the system in the mean field approximation has a unique attractor, then the tensors of Theorem 1 have a limit as t goes to infinity, and this limit can be used to obtain an asymptotic expansion in 1/N and 1/N 2 in steady-state. For V and W , these equations are the same as ones developed in [START_REF] Gast | A Refined Mean Field Approximation[END_REF]. The novelty of this result is the 1/N 2 -expansion.

Theorem 3. In addition to the assumption of Theorem 1, assume (A3) and (A3). Then the ODE of Theorem 1 also has a unique attractor. Moreover, in steady state for any four times differentiable function h : R d → R, one has:

E h(X (N ) ) = h(π ) + 1 N H i V i + 1 2 H i j W i j + 1 N 2 H i A i + 1 2 H i j B i j + 1 6 H i jk C i jk + 1 24 H i jk ℓ D i jkl + o 1 N 2 ,
where the terms H i . . . H i jk ℓ denotes the first to fourth derivative of h evaluated at the fixed point π and where the tensors satisfy the following system of linear equations:

2Sym F i k W k j = -Q i j F i j V j = 1 2 F i jk W jk and 4Sym(F i m D mjk ℓ ) = -6Sym(Q i j W k ℓ ) 3Sym(F i ℓ C ℓjk ) = -Sym 3 2 F i ℓm D ℓmjk + 3Q i j V k + 3Q i j ℓ W ℓk + R i jk 2Sym(F i k B k j ) = -Sym F i k ℓ C k ℓj + 1 3 F i k ℓm D k ℓmj +Q i j k V k + 1 2 Q i j k ℓ W k ℓ F i j A j = - 1 2 F i jk B jk + 1 6 F i jk ℓ C jk ℓ + 1 24 F i jk ℓm D jk ℓm
where Q, R, Q k and Q k ℓ are evaluated at the fixed point π .

Also, as we will see in the proof, under the condition of Theorem 3, the convergence as N goes to infinity of Equation ( 3) is uniform in time. This is not necessarily the case when the mean field approximation does not have an attractor (see Section 7).

Computational Issues and Implementation

Transient Analysis.

For a given mean field model, the ODE x = f (x) is an ODE of dimension d. As the drift f is in general nonlinear, the solution x(t) can rarely be computed in closed form but can be easily computed numerically for high dimensional models. Once the solution x(t) is computed, the system of ODEs for V , W , A, B, C and D given by Theorem 1 is a system of linear ODEs with time-varying parameters.

The system of ODEs for V and W do not depend on A, B, C, D. It is therefore possible to compute the 1/N terms V (t) and W (t) by numerically integrating a system of O(d 2 ) variables. The computation of the 1/N 2 terms is more complicated because D has d 4 variables. This makes the computation of the 1/N 2 terms feasible for d of at most a few tens.

Fixed-Point Analysis.

The computation of the fixed point of Theorem 3 can also be solved by a numerical algorithm: The constants V to D are the solutions of a system of linear equations.

For the 1/N -term, these equations are the same as the ones developed in [START_REF] Gast | A Refined Mean Field Approximation[END_REF] and can therefore be solved in O(d 3 ) time in two steps:

• First, we obtain the matrix W from the solution of the Lyapunov equation MW + (MW ) T = Q for some matrix M. • Second, the vector V is the solution of a linear system of equations of dimension d.

The most costly step of the above is the computation of the solution of the Lyapunov equation, which can be done in O(d 3 ) time by using the Bartels-Stewart algorithm [START_REF] Bartels | Solution of the matrix equation AX+ XB= C [F4][END_REF].

Once the terms V and W have been computed, one can compute the tensors D, C, B, A (in this order) by exploiting the fact that the equation for D does not depend A, B, C (similarly, the equation for C does not depend on A and B; the equation for B does not depend on A). Each is a system of linear equations with respectively d 4 , d 3 , d 2 and d variables. For D and C, the system is a generalization of the classical Lyapunov equation MW + (MW ) T = Q to higher order tensors. Although the system of linear equations is large, in our numerical examples we were able to solve these equations for system as large as d = 50 dimensions in less than 20 seconds (which corresponds for D to a linear system with 50 4 = 6.25 × 10 6 unknowns).

Implementation.

To compute numerically the mean field expansions, we implemented a generic tool in Python that can construct and solve the above equation. The tool is available at https: //github.com/ngast/rmf_tool/ [START_REF] Gast | Refined Mean Field Tool[END_REF]. It takes as an input a description of the model and uses symbolic differentiation to construct the derivatives of the drift and of the functions Q and R.

The tool uses the function integrate.solve_ivp of the library scipy [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF] to numerically integrate the ODEs for computing V (t) and W (t) of Theorem 1. For the steady-state analysis, the tool uses the python library scipy.sparse to construct a sparse system of linear equations and the function scipy.sparse.linalg.lgmres to solve the sparse linear system.

Note that the use of symbolic differentiation makes the computation slow for large models. Hence, for the supermarket model, we directly implemented Python functions that compute the drift of the system and its derivative. All our specific implementation is available in the git repository of the paper [START_REF] Gast | Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis[END_REF].

Analysis of the computation time.

To give a flavor of the numerical complexity of the method, we report in Figure 1 We observe that, as expected, computing the time-varying constants of the transient regime is more costly than solving the fixed point equations because it requires solving an ODE: for a given time budget, one can compute the steady-state constants for a system of doubled size. Moreover, these results show that the computation of the 1/N -terms V (t) and W (t) can be done for models with hundreds of dimensions in 10 seconds. With the same constraints of 10 seconds, the 1/N 2 -terms can be computed for models with a few tens of dimensions.

Note that we only provide this figure for the supermarket model because, among our three examples, it is the only one for which we can vary the dimension by changing the maximal queue lengths. We believe that the computation time does not grow too much with the dimension because the tensors corresponding to the derivatives of the drift or of the matrix Q are relatively sparse. The computation time might be higher for a model with denser tensors.

Proofs

To simplify the notation, where it is not needed in the proofs, we drop the superscript N and denote X instead of X (N ) . 4.4.1 Proof of Theorem 1. The proof of Theorem 1 is divided in two parts. We first we show the existence of the constants A, B,. . . Second we show how to derive the ODE that they satisfy.

Existence of V , W -Here, we again use the notation Φ s x to denote the value at time s of the solution of the ODE x = f (x) that starts in x at time 0. According to [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF]Equation (19)] for any function h : E → R, we have

N E h(Φ t x) -h(X (N ) (t)) | X (N ) 0 = x = ∫ t 0 E ∆ (N ) h • Φ s (X (N ) (t -s)) | X (N ) 0 = x ds, (6)
where ∆ (N ) is the operator that, for a function д, gives the function ∆ (N ) д defined by:

(∆ (N ) д)(x) = N ℓ ∈ L β ℓ (x)(N (д(x + ℓ N ) -д(x)) -д j (x)ℓ j ), (7)
where we recall the use of Einstein summation convention :

д j (x)ℓ j = d j=1 (∂д(x))/(∂x j )ℓ j .
By using a Taylor expansion of д in the above equation, for a function д : R n → R that is twice differentiable, we have

∆ (N ) д(x) = 1 2 ℓ ∈L β ℓ (x)д i j (x)ℓ i ℓ j + o(1/N ), (8)
where the hidden constant in the o(1/N ) depends on the modulus of continuity of the second derivative of д. This shows that Equation ( 6) equals

E       ∫ t 0 1 2 ℓ ∈L β ℓ (X (N ) (t-s))(h • Φ s ) i j (X (N ) (t-s))ℓ i ℓ j ds       + o 1 N ,
where (h • Φ s ) i j (Φ t -s x) denotes the second derivative of h • Φ s with respect to x i and x j evaluated at Φ t -s x. Again, the hidden constant in the o(1/N ) depends on the modulus of continuity of the second derivative of (h • Φ s ) which is finite for any time t because of Assumption (A2). As X (N ) (t-s) converges weakly to Φ t -s x as N goes to infinity, the above quantity (to which Eq.( 6) is equal) yields Eq.( 6

) = ∫ t 0 1 2 ℓ ∈L β ℓ (Φ t -s x)(h • Φ s ) i j (Φ t -s x)ℓ i ℓ j ds + o(1/N ). (9)
In the quantity (h • Φ s ) i j (Φ t -s x), the only dependence in h is a linear combination of the first and second derivative of h evaluated at Φ t x. Indeed, by the chain rule, for two functions д and h, the first and second derivative of д • h evaluated in y is

(h • д) i = (h k • д)д k i (h • д) i j = (h k ℓ • д)д k i д ℓ j + (h k • д)д k i j
Replacing д by Φ s x and evaluating the function is Φ t -s x shows that the second derivative of h • Φ s evaluated in Φ t -s x is :

(h • Φ s ) i j (Φ t -s x) = h k ℓ (Φ t x)(Φ s ) k i (Φ t -s (x))(Φ s ) k j (Φ t -s (x)) + h k (Φ t x)(Φ s ) k i, j (Φ t -s (x)).
Plugging this into Equation [START_REF] Chaintreau | The Age of Gossip: Spatial Mean Field Regime[END_REF] shows that Equation ( 6) is equal to:

h k ℓ (Φ t x) 1 2 ∫ t 0 ℓ ∈L β ℓ (Φ t -s x)(Φ s ) k i (Φ t -s (x))(Φ s ) k j (Φ t -s (x))ℓ i ℓ j ds =:W k ℓ (t ) + h k (Φ t x) 1 2 ∫ t 0 ℓ ∈L β ℓ (Φ s ) k i, j (Φ t -s (x))ℓ i ℓ j ds =:V k (t ) +o(1/N ).
This implies the existence of V (t) and W (t) in Equation (3). Existence of A. . . D -The proof of the existence of the terms A to D is similar. Hence, for space constraints we only sketch the main differences. The first ideas is to refined the expansion [START_REF] Cecchi | Mean-field analysis of ultra-dense csma networks[END_REF] to This shows that Equation ( 6) equals

∆ (N ) д(x) = 1 2 ℓ ∈L β ℓ (x)д i j ℓ i ℓ j + 1 6N ℓ ∈L β ℓ (x)д i jk ℓ i ℓ j ℓ k + o( 1 N 2 ). (10
E       ∫ t 0 1 2 ℓ ∈ L β ℓ (X (N ) (t-s))(h • Φ s ) i j (X (N ) (t-s))ℓ i ℓ j ds       +E       1 6N ∫ t 0 ℓ ∈ L β ℓ (X (N ) (t-s))(h • Φ s ) i jk (X (N ) (t-s))ℓ i ℓ j ℓ k ds       +o 1 N 2 .
In the above equation, the second term is of order 1/N and involves the derivative up to order three of h. The first term is equal to (6) plus a correction term of order 1/N that involves the derivative up order four of h (evaluated at Φ t x).

Derivation of the ODEs -The evolution of the stochastic process X (t) -x(t) can be decomposed in two parts : a jump part due to the fact that X (t) jumps to X (t) + ℓ/N at rate N β ℓ (X (t)) and a drift part due to the fact x(t) satisfies the ODE x = f (x). This shows that for any function h, one has :

d dt E [h(X (t) -x(t))] = ℓ ∈ L E h(X (t) -x(t) + ℓ N ) -h(X (t) -x(t)) N β ℓ (X (t)) -E h j (X (t) -x(t))f j (x(t)) .
In the above equation, the first line corresponds to the stochastic jumps of X (t) while the second line corresponds the continuous variation of x(t).

Applying the above equation1 to the function h(X ) = (X -x) ⊗k shows that :

d dt E (X -x) ⊗k (11) = ℓ ∈ L E X -x + ℓ N ⊗k -(X -x) ⊗k N β ℓ (X ) -kSym f (x) ⊗ E (X -x) ⊗k -1 = k m=1 k m Sym E 1 N m-1 ℓ ⊗m β ℓ (X ) ⊗ (X -x) ⊗k -m -kSym f (x) ⊗ E (X -x) ⊗k -1 = kSym E (f (X ) -f (x)) ⊗ (X -x) ⊗k -1 + k m=2 k m Sym E 1 N m-1 ℓ ⊗m β ℓ (X ) ⊗ (X -x) ⊗k-m (12)
The the existence of the constants V , W , A. . . D combined with Equation [START_REF] Eldering | Normally Hyperbolic Invariant Manifolds-the Noncompact Case[END_REF] show that the derivative of E[(X (t) -x(t)) ⊗k ] admits an asymptotic expansion with a first term in 1/N and a second term in 1/N 2 . We are now ready to compute how the constants V , W , A... D evolve with time by computing the derivative with respect to time of E[(X -x) ⊗k ] for k ∈ {1 . . . 4} and identifying the 1/N and 1/N 2 terms.

1.

Case E[X -x] -By using Equation ( 12), we have :

d dt E [X -x] = E [f (X ) -f (x)]. Applying (3) to the function h(X ) = f i (X ) -f i (x(t)) implies that d dt E X i -x i = 1 N (F i j V j + 1 2 F i jk W jk ) + 1 N 2 (F i j A j + 1 2 F i jk B jk + 1 6 F i jkl C jkl + 1 24 F i jklm D jklm ) + o( 1 N 2 ) Using that d dt E[X i -x i ] = V i /N +A i /N 2 +o(1/N 2 )
and identifying the O(1/N ) and O(1/N 2 ) terms shows that:

d dt V i = F i j V j + 1 2 F i jk W jk d dt A i = F i j A j + 1 2 F i jk B jk + 1 6 F i jk ℓ C jk ℓ + 1 24 F i jk ℓm D jk ℓm 2. Case E[(X -x) ⊗2 ] -By using (12), we have d dt E (X -x) ⊗2 = 2Sym(E [(f (X ) -f (x)) ⊗ (X -x)]) + 1 N E [q(X )],
where q(X ) is a covariance matrix defined by

q(X ) = ℓ ∈L β ℓ (X )ℓ ⊗ ℓ
For the first term, we consider the function 2 h(X ) = ((f (X )-f (x))⊗ (X -x)) i j and we use Lemma 4(i). The first derivative of this function h evaluated at x is 0. The second derivative of h with respect to

x k and x ℓ is 2Sym(F k ⊗ J (ℓ) ), where J (ℓ) is the matrix whose (ℓ, ℓ)element is one, the others being zero. The third derivative with respect to x k , x ℓ and x m is equal to 3Sym(F k ℓ ⊗ J (m) ). The fourth derivative with respect to x k , x ℓ , x m and x n is 4Sym(F k ℓm ⊗ J (n) ).

Hence, applying Equation (3) to h hows that

d dt E Sym(((f (X ) -f (x)) ⊗ (X -x)) i j ) = Sym 2 2 F i k ( 1 N W k j + 1 N 2 B k j ) + 3 6N 2 F i k ℓ C k ℓj + 4 24N 2 F i k ℓm D k ℓmj
For the second term, applying (3) to the function q shows that E q i j (X )

= Q i j + 1 N Q i j k V k + 1 2N Q i j k ℓ W k ℓ + O(1/N 2 ).
This shows that:

W i j = 2Sym(F i k W k j ) + Q i j B i j = Sym 2F i k B k j + F i k ℓ C k ℓj + 1 3 F i k ℓm D k ℓmj + Q i j k V k + 1 2 Q i j k ℓ W k ℓ .
3. By using [START_REF] Eldering | Normally Hyperbolic Invariant Manifolds-the Noncompact Case[END_REF], with E[(X -x) ⊗3 ], we have :

d dt E (X -x) ⊗3 = 3Sym(E (f (X ) -f (x)) ⊗ (X -x) ⊗2 ) + 3 N Sym(E [q(X ) ⊗ (X -x)]) + 1 N 2 E [r (X )], (13)
where r (x) = ℓ ∈L ℓ ⊗3 β ℓ (x).

2 Recall that the exponent i j stands for the component (i j).

To study the first term of Equation ( 13), we consider the function h(X ) = ((f (X ) -f (x)) ⊗ (X -x) ⊗2 ) i jk . Applying Lemma 4(ii), the first two derivatives of this function evaluated at x are equal to 0. The third derivative of this function (with respect to x ℓ , x m , x n ) is equal to 6Sym(F ℓ ⊗ J (m) ⊗ J (n) ) and the fourth derivative is equal to 12Sym(F ℓm ⊗ J (m) ⊗ J (n) ⊗ J (o) ).

Hence, applying Equation ( 3) to h shows that

E Sym((f (X ) -f (x)) ⊗ (X -x) ⊗2 ) = 1 N 2 Sym 6 6 F i ℓ C ℓjk + 12 24 F i ℓm D ℓmjk + o(1/N 2 ).
The second term of Equation ( 13) can be treated by applying Equation (3) to h(X ) = q(X )(X kx k ), whose first derivative evaluated at x is Q i j and whose second derivative is 2Q

i j ℓ ⊗ J (ℓ) (see Lemma 4(i)). This shows that E Sym(q(X ) ⊗ (X -x) i jk ) = 1 N Sym Q i j V k + 2 2 Q i j ℓ W ℓk + o( 1 N ).
Finally, the last term of Equation ( 13) is equal to R/N 2 +o(1/N 2 ). This shows that Equation ( 13) has only terms in O(1/N 2 ) plus term of order o(1/N 2 ). By identifying the O(1/N 2 )-terms, we get

C i jk = 3Sym(F i ℓ C ℓjk ) + 3 2 Sym(F i ℓm D ℓmjk ) + 3Sym(Q i j V k ) + 3Sym(Q i j ℓ W ℓk ) + R i jk 4. The derivative is similar for E[Y ⊗4 t ]. Applying (12) shows that d dt E (X -x) ⊗4 = 4Sym(E (f (X ) -f (x)) ⊗ (X -x) ⊗3 ) + 6 N E q(X ) ⊗ (X -x) ⊗2 4 N 2 E r (X ) ⊗ (X -x) ⊗2 + 1 N 3 E       ℓ ∈ L β ℓ (X )ℓ ⊗4       By (3) with the function h(x) = (f (X ) -f (x))(X -x) ⊗3 the first term is equal to 4Sym(F i m D mjk ℓ )/N 2 + o(1/N 2
) (because the first three derivatives of this function h are equal to zero and the last one has a factor 4 × 3 × 2 = 24 by Lemma 4(iii)).

For the second term, we can again use Equation (3) with h(X ) = q(X )(X -x) 2 and Lemma 4(ii). The first derivative of h is zero and only the second term counts :

Sym(E q(X ) ⊗ (X -x) ⊗2 ) = 2 2N Sym(Q ⊗ W ) + o(1/N ).
Finally, the one before last is of order O(1/N 3 ) because of (3) and the last term is of order O(1/N 3 ). We therefore obtain :

D i jk ℓ = 4Sym(F i m D mjk ℓ ) + 6Sym(Q k ℓ W i j ).
In the above proof, we used the following lemma, whose proof is direct by using general Leibniz rule. Lemma 4. Let д : R → R be k-times differentiable. Then

(i) ∂ k (xд(x )) (∂x ) k = xд (k ) (x) + kд (k-1) (x) (ii) ∂ k (x 2 д(x )) (∂x ) k = x 2 д (k ) (x) + 2kxд (k -1) (x) + k(k -1)д (k-2) (iii) ∂ k (x 3 д(x )) (∂x ) k = x 3 д (k ) (x) + 3kx 2 д (k -1) (x) + 3k(k -1)xд (k -2) +k(k -1)(k -2)д (k -3)
4.4.2 Proof of Theorem 3. Most of the work needed to prove Theorem 3 was already done in the proof of Theorem 1. Indeed, it should be clear the linear equations of Theorem 3 correspond to the fixed point equation of the ODE of Theorem 1. Therefore, to prove Theorem 3, the only remaining steps are to prove that:

(1) These fixed point equations have a unique solution.

(2) The system of ODEs of Theorem 1 converges to this solution.

(3) One can exchange the limits lim t →∞ and lim N →∞ .

Uniqueness -the uniqueness of the solution, for V and W was already shown in [START_REF] Gast | A Refined Mean Field Approximation[END_REF]. For D, one can remark that its fixed point equation can be written as a matrix equation M (4) D = y where y is a vectorized version of -6Sym(Q ⊗ W ), and where the matrix M (4) is a d 4 × d 4 matrix that can be expressed as the Kronecker sum of four times the Jacobian of the drift evaluated at π :

M (4) i jk ℓ;abcd = F i a δ jb δ kc δ ℓd + δ ia F j b δ kc δ ℓd + δ ia δ jb F k c δ ℓd + δ ia δ jb δ kc F ℓ d , (14)
where δ i j is the Kronecker symbol that equals 1 if an only if i = j and 0 otherwise. Note that in the above equation, the lines and columns of the matrix M (4) are indexed by the tuples ijkℓ (for the lines) or abcd (for the columns). By property the Kronecker sum, an eigenvalue of M (4) is the sum of four eigenvalues of the Jacobian matrix (F i j ). As the system is exponentially stable, all the eigenvalues of the Jacobian matrix have negative real part. Therefore all eigenvalues of the matrix M (4) have negative real part and M (4) is invertible. This implies the existence and the uniqueness of the solution for D of the fixed point equation.

Once the D is fixed, the equation for C can be written is a similar way M (3) C = y where M (3) is the Kronecker sum of three times the Jacobian of the drift. A similar reasoning as the one for D shows that C is uniquely defined. This can be propagated to B and then A.

Convergence to the fixed point. The time-varying constant D(t) satisfies a time-inhomogeneous linear differential equation D = M (4) (t)D + y(t), where M (4) (t) is the Kronecker sum of four times the Jacobian of the drift evaluated in x(t) and y(t) (defined as in Equation ( 14)). As x(t) converges to an exponentially stable attractor π , all eigenvalues of the Jacobian of the drift f evaluated in π have negative real part. This implies that there exists a time after which all eigenvalues of the Jacobian of f have negative real part in which case all eigenvalues of the matrix M (4) (t) have negative real part. This implies that the ODE for D(t) is exponentially stable and that therefore D(t) converges to the unique fixed point of this system. The same reasoning applies for C, B and A.

Exchange of the limits. The above steps guarantee that the terms V (t) and A(t) of the development in 1/N and 1/N 2 converge as N goes to infinity. Informally, this shows that lim

t →∞ E X N (t) -π = lim t →∞ E X N (t) -x(t) = lim t →∞ 1 N V (t) + 1 N 2 A(t) + o(1/N 2 ) = 1 N V + 1 N 2 A + lim t →∞ o(1/N 2 ). (15)
In order to conclude the proof, we need to show that it is possible to exchange the limits, which is to show that the term lim

t →∞ o(1/N 2 ) is indeed a o(1/N 2 ) term.
To see that, we use Stein's method and the ideas developed in [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF] to show that, in steady-state,

E h(X (N ) ) -h(π ) = E ∆ (N ) ∫ ∞ 0 h(Φ s (X (N ) (s))) -h(π )ds ,
where ∆ (N ) is the operator defined in Equation [START_REF] Braverman | Stein's method for steady-state diffusion approximations: an introduction through the Erlang-A and Erlang-C models[END_REF]. Note that this equation is a consequence of Equation ( 10) of [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF] and is the analog of Equation ( 6) as t goes to infinity.

Concerning the exchangeability of the limits, for space constraints, we only sketch the main remaining ideas of the proofs. The first step is to show that the hidden constant of the o(1/N 2 ) of Theorem 1 depends on the modulus of continuity of the function

G (t ) (x) = ∫ t 0 h(Φ s (x)) -h(π )ds.
This comes from Equation [START_REF] Chen | A projection method and Kronecker product preconditioner for solving Sylvester tensor equations[END_REF]. The second idea is that the function G(x) = ∫ ∞ 0 h(Φ s (x)) -h(π )ds is four times differentiable and that the derivatives G (t ) converge uniformly to the derivatives of G as t goes to infinity. This comes from perturbation theory : by [START_REF] Eldering | Normally Hyperbolic Invariant Manifolds-the Noncompact Case[END_REF]Lemma C.1], if the flow Φ has an exponentially stable attractor and is four times differentiable, then the first four derivatives of Φ s (x) converge exponentially fast to 0. The same argument is used in the proof of Lemma 3.5 of [START_REF] Gast | A Refined Mean Field Approximation[END_REF]. These two arguments show that the modulus of continuity of the derivatives of G (t ) are uniformly bounded in time and that therefore the convergence is uniform in time.

EXAMPLE 1: MALWARE PROPAGATION

In this section we illustrate the above results with a simplified variant of the malware propagation model of [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Mhr Khouzani | Maximum damage malware attack in mobile wireless networks[END_REF]. It can be viewed as an instance of a basic infection model in epidemiology (e.g., [START_REF] Murray | Mathematical Biology I: An Introduction[END_REF]). We choose this model because of its simplicity: since it is a one-dimensional model, the constants of the 1/N and the 1/N 2 approximation can be computed in closed form and the stationary distribution can be evaluated numerically easily with high precision (it is a birth-death process). This allows us to assess the accuracy of the various approximations with high precision.

Model

We consider a model of malware propagation in a system composed of N agents. Each agent is either infected by the malware or not. Let X be the fraction of infected agents. We consider that each noninfected agent becomes infected at rate 1+X (the rate 1 corresponds to infection by an external source while the rate X corresponds an infection by a peer). An infected agent recovers at rate 1 due to some patching mechanism. This translates into the following transitions for X :

X → X + 1 N at rate N (1 -X )(1 + X ) X → X - 1 N at rate N X

Mean Field Approximations and Expansions

To apply Theorem 1 and 3, let us first compute the drift of the system, its derivative, the matrix Q and its derivative, and the tensor R. As the system is uni-dimensional, all tensors are in fact scalars. The drift is f

(x) = 1 -x 2 -x = r (x) and the function q(x) = 1 -x 2 + x.
The ODE of the mean field approximation x = f (x) is a Bernoulli type equation, hence, the mean field approximation has the closedform solution

x(t) = - 1 2 + √ 5 2 2 1 -αe - √ 5t -1 , (16)
where α = (4x(0

) + 1 - √ 5)/(4x(0) + 1 + √ 5) and x(0) is the initial condition.
As there is a close form solution for the mean field approximation, it might be doable to obtain a close form expression for the constants V (t), W (t),. . . but the expressions of such constant seem highly complex. Hence, in our illustrations, we use our tool [START_REF] Gast | Refined Mean Field Tool[END_REF] to compute numerically these constants.

The fixed point analysis is simpler. From Equation ( 16), it is clear that the ODE x = f (x) has a unique attractor π = ( √ 5 -1)/2 that is exponentially stable. Moreover, the derivatives of the drift (evaluated at π 4) (π ) = 0. Finally, the function q evaluated at π is q(π ) = √ 5 -1 and its derivatives are q ′ (π ) = 2 -√ 5, q ′′ (π ) = -2. Last, we have that r (π ) = 0.

) are f ′ (π ) = - √ 5, f ′′ (π ) = 2, f (3) (π ) = f (
After some algebra, it can be shown that the constants V and A that solve the fixed point equation of Theorem 3 are

V = √ 5 -1 10 and A = √ 5 -3 50
.
Plugging the above quantity into Theorem 3 shows that, in steadystate and as N goes to infinity, one has :

E [X ] = √ 5 -1 2 1 - 1 5N + √ 5 -3 50N 2 + o 1 N 2 .

Numerical Comparison

In this section, we propose a numerical comparison of the exact values, the mean field approximation and the two expansions (up to order 1/N and 1/N 2 ).

Transient regime.

To perform a numerical comparison of the various approximations with the exact values, we implemented two numerical procedures. For the mean field approximation and the expansions, we implemented a numerical integration of the system of ODEs of Theorem 1. For the exact values, we used the fact that for a given size N , the stochastic model is a continuous time Markov chain with N + 1 states ({0, 1/N , 2/N . . . , 1}). We again used a numerical integrator to integrate the Kolmogorov equations for this case. The results are reported in Figure 2 in which we compare the three approximations (mean field and the two expansions) with the exact values, for N = 5 and N = 10. At the beginning, we start in a system where X (0) = 0.6 (i.e. 3N /5 of the N agents are infected). We observe that the expansions provide a much better characterization of the transient regime that the classical mean field approximation. Note that for N = 5, the gain when going from the 1/N to the 1/N 2 is small. For N = 10, the gain is almost invisible.

N 2 (E[X] X V/N) N=5 N=10 N=20 N=40 Constant A
To observe more precisely what is the gain brought by the 1/N 2 approximation, we plot in Figure 3 the 1/N 2 -constant A(t) and compare it with the error of the 1/N -expansion rescaled by N 2 : N 2 (E[X (t)]-x(t)-V (t)/N ), for various values of N ∈ {5, 10, 20, 30}.

As shown by Theorem 1, the rescaled error of the 1/N -expansion converges to A(t) as N goes to infinity. This figures also shows that A(t) is of order 10 -2 . This explains why the gain in accuracy brought by the 1/N 2 -term is small: the error of the 1/N -approximation is only around 0.01/N 2 .

N E[X ]

1/N -expansion 1/N 2 -expansion Error Error 1 0.5000000 0.4944272 5.6e-03 0.4791486 2.1e-02 5 0.5929041 0.5933126 -4.1e-04 0.5927015 2.0e-04 10 0.6055449 0.6056733 -1.3e-04 0.6055205 2.4e-05 20 0.6118184 0.6118536 -3.5e-05 0.6118155 3.0e-06 30 0.6138977 0.6139138 -1.6e-05 0.6138968 8.7e-07 50 0.6155559 0.6155619 -5.9e-06 0.6155557 1.9e-07 ∞ 0.6180340 0.6180340 0 0.6180340 0 Table 1: Malware propagation model: comparison of the "true" expectation of X and the 1/N and 1/N 2 expansions. The "error" column is the difference between E[X ] and the expansion. Note that the classical mean field approximation is the value for N = ∞, which is π ≈ 0.6180340.

Steady-state.

We now verify the accuracy in steady-state. In Table 1, we verify the accuracy of the approximation for various values of N ∈ {1, 5, 10, 20, 30, 50}. We compare three values :

• E[X ] that we computed by using the fact that this model is a birth-death process whose stationary measure can therefore be easily computed numerically. • π + V /N , which is the refined approximation of [START_REF] Gast | Expected values estimated via mean-field approximation are 1/N-accurate[END_REF] and that we call the 1/N -expansion.

• π + V /N + A/N 2 that we call the 1/N 2 -expansion.
We observe that for this model, the 1/N and 1/N 2 expansions are already very accurate for N = 1 and they soon provide more than 4 digits of precision for N ≥ 10. For N ≥ 10, the error made by the 1/N 2 -expansion is an order of magnitude smaller than the error made by the 1/N -expansion (the ratio between the two errors is approximately 0.6N ). The high accuracy of the 1/N -expansion can be by the fact that the two constants are V ≈ 0.12 and A ≈ -0.015, hence, as for the transient regime, the difference between the two expansions is only 0.015/N 2 .

THE SUPERMARKET MODEL

We now focus on the classical supermarket model of [START_REF] David | The Power of Two Random Choices in Randomized Load Balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF]. We study the gain of the 1/N and 1/N 2 expansions for the transient and the steady-state regimes. As for the previous examples, the gain in accuracy of the 1/N -expansion over the mean field approximation is large but the gain of the 1/N 2 -expansion over the 1/N -expansion is smaller. Also, this model illustrates that it is possible to compute the 1/N and 1/N 2 terms for a realistic model.

The Model

We consider a queuing system composed of N identical servers. Jobs arrive at a central broker according to a Poisson process of rate ρN and are dispatched towards the servers by using the JSQ(k) policy: for each incoming job, the broker samples k servers at random and sends the jobs to the server that has the smallest number of jobs in its queue (ties are broken at random). The time to process a job is exponentially distributed with mean 1.

This system can be modeled as a density dependent population process defined in Section 3. To see that, we assume that the queue size is bounded by d and we denote by X i (t) the fraction of servers with queue size i or more at time t. X (t) is a Markov chain whose transitions are :

X → X -1 N e i at rate N (X i -X i+1 ) X → X + 1 N e i at rate N ρ(X k i-1 -X k i ), (17)
where e i is a vector whose ith component is 1 the other ones being 0. Also, note that we use the classical notation for indices : X i denotes the ith component of X and X k i denotes the kth power of X i . The explanation is as follows: A departure from a server with i ≥ 1 jobs modifies X into X -N -1 e i and occurs at rate N (X i -X i+1 ). An arrival at a server with i jobs modifies X into X + N -1 e i . Assuming that the k servers are picked with replacement, the least loaded among k servers has i -1 jobs with probability X k i-1 -X k i .

Mean Field Approximation and Expansions

To apply Theorems 1 and 3, we first compute the drift, the constants Q, R and the needed derivatives.

The ith component of the drift of this model evaluated at x is F i :

F i = ρ(x k i-1 -x k i ) + (x i+1 -x i ). (18)
The first derivative of the drift evaluated at a point x satisfies

F i i-1 = kρx k -1 i-1 ; F i i = -kρx k -1 i -1; F i i+1 = 1,
all other terms being equal to 0. Similarly, the second derivative satisfies

F i i-1,i-1 = k(k -1)ρx k -2 i-1 F i ii = -k(k -1)ρx k -2 i ,
all other terms being equal to 0. The expression is similar for the third and fourth derivatives. The tensors Q and R of Equation ( 4) and ( 5) satisfy:

Q ii = (β e i (x) + β -e i (x)) = ρ(x k i-1 -x k i ) + (x i -x i+1 ) R iii = F i = ρ(x k i-1 -x k i ) + (x i+1 -x i ).
Finally, the first and second derivatives of q evaluated in x satisfy

Q ii i-1 = kρx k -1 i-1 Q ii i = 1 -kρx k -1 i Q ii i+1 = -1 Q ii i-1,i-1 = k(k -1)ρx k -2 i-1 Q ii ii = -k(k -1)ρx k -2 i
To apply Theorem 3, the only technical condition to verify is that the fixed point is exponentially stable. This is done for example in [START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF][START_REF] Ying | On the Rate of Convergence of the Power-of-Two-Choices to its Mean-Field Limit[END_REF]. The constants for the steady-state approximation can be computed by evaluating the above equation in π .

Algorithmic Considerations

In order to perform numerical comparison of the refined approximations and an estimation of the true values, we implemented various numerical algorithms. For the expected values, we implemented a C++ simulator of the supermarket model that simulates a densitydependent population process whose transitions are exactly the ones of Equation [START_REF] Gast | Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis[END_REF]. For the transient analysis, to estimate the evolution of the expected queue length as a function of time, we performed an average of 10 5 (for N = 10) or 20000 (for N = 20) independent runs of simulations. This number of simulations is chosen as a compromise between computation time and accuracy. As we will observe in Figure 4, more simulations would give more accurate results but we choose to limit the computation time to 1h per panel. For the steady-state values, we compute the average of 1000 independent time-average of simulations after a warp-up period of 10000N events for each.

For the numerical analysis, we implemented a code to compute the parameters of the supermarket model and then use our tool [START_REF] Gast | Refined Mean Field Tool[END_REF] to solve numerically the ODEs of Theorem 1 or the fixed point equations of Theorem 3. As the size of the ODE for the 1/N 2approximation grows like d 4 , we choose to bound the queue length to d = 10 for the 1/N 2 -expansions. In practice, using a larger maximal queue length brings to the same numerical value. For the transient regime, the computation time of the 1/N 2 -term is around 10sec and the one of the 1/N -term less than one second. The computation of the fixed point is much faster than the one of the transient regime: it takes around 300ms for d = 20 and around 15s for d = 50 (on a 2013-laptop).

Numerical Comparisons

It is shown in [START_REF] Gast | A Refined Mean Field Approximation[END_REF] that the 1/N -expansion provides estimates of the steady-state average queue length that are much more accurate than the classical mean field approximation. In this section we show that the 1/N -expansion can also be used to improve the accuracy in the transient-regime and that the 1/N 2 -expansion improves on the 1/N -expansion (both for transient and steady-state analysis).

6.4.1 Transient regime. We first consider how the expected queue length evolves with time. We consider the supermarket model with k = 2 choices and ρ = 0.9. We start in a system where the expected queue length is 2.8 : out of the N queues, 0.2N queues start with 2 jobs and 0.8N queues start with 3 jobs. We choose this value as it is close to 2.75, the steady-state average queue length predicted by the 1/N -expansion for N = 10.

In Figure 4, we report how the expected queue length evolve with time compared to the three approximation (mean field, 1/Napproximation and 1/N 2 -approximation). We observe in this figure that both for N = 10 and N = 20, the expansions provide an estimation of the evolution of the expected queue length that is much more accurate than the one provided by the classical mean field approximation. Moreover, for N = 10, the 1/N 2 -expansion provides a better approximation than the 1/N -expansion. For N = 20, the two curves are almost indistinguishable.

For the simulation of the transient regime, the running time of simulation is approximately 0.1sec per run of our C++ simulator for N = 20 and 0.05sec for N = 10. This represents roughly 1h of computation for each of the two panels. As a comparison, the total time to compute the expansion of order 1/N 2 is about 10 seconds (and does not depend on N ), and the time to compute the expansion of order 1/N is around 1 second (using our python's implementation).

Note that we only present the results for k = 2 and ρ = 0.9. Similar results can be observed for other values of k and ρ with one difference: the smaller is ρ, the smaller is the difference between the approximations and the simulation (the difference between the 1/Nexpansion and the 1/N 2 -expansion can almost not be distinguished for ρ < 0.7). This is more visible in Table 2 complete results can be found in the git repository of the paper [START_REF] Gast | Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis[END_REF].

We observe that in all tested cases, the 1/N -expansion provides an estimation of the average queue length that is much more accurate than the one provided by the classical mean field approximation. The estimation provided by the 1/N 2 -expansion is generally more accurate but the gain brought by the 1/N 2 -term varies across the different parameters. The gain is the most visible for k = 2, in which case the 1/N 2 -expansion provides very accurate estimates, even for N = 10. This is less pronounced for k = 3 and k = 4, where the gain is more visible for higher values of N . Recall that in all cases, the mean field approximation provides estimates that do not depend on the system size N . They are systematically less accurate than the two expansions.

Theorem 3 can also be used to compute estimations of the queue length distribution. Indeed, for the supermarket model, E[X i ] is the probability that a given server has i jobs or more. In Table 3 Mean field 0.590 0.109 0.000 0.000 0.000 0.000 1/N -expansion 0.679 0.450 0.007 0.000 0.000 0.000 1/N 2 -expansion 0.652 0.341 0.131 0.000 0.000 0.000 Simulation 0.657 0.344 0.140 0.051 0.018 0.006 ρ=0.95, k=4, N =20 Mean field 0.774 0.341 0.013 0.000 0.000 0.000 1/N -expansion 0.802 0.600 0.178 0.000 0.000 0.000 1/N 2 -expansion 0.795 0.429 0.578 0.001 0.000 0.000 Simulation 0.798 0.509 0.236 0.092 0.034 0.012 Table 3: Supermarket : steady-state distribution.

which is the probability that a server is busy and is equal to ρ. We make two observations. First, for moderate values of ρ and k, the 1/N 2 -expansion provides a very accurate estimation of the "true" distribution that we estimate by using simulation. This is less clear for higher values such as k = 4 and ρ = 0.95 for which the 1/N 2 terms has a tendency to over-correct for small values of N . Also, in all tested cases, the values for moderate values of i are well approximated, but the tail of the distribution is less well approximated. Note that for a fixed set of parameters (ρ, d), the two expansions become more accurate as N grows. This is illustrated in the git repository of the paper [START_REF] Gast | Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis[END_REF].

LIMITATIONS OF THE APPROACH

In the previous examples, we concentrated on cases where the mean field approximation has a unique attractor, which implies that the mean field approximation and its expansions converge to the exact value of E[h(X )] uniformly in time (Theorem 3). In this section, we show that when the mean field approximation has a fixed point that is not an global attractor, this does not hold anymore. Moreover, in this setting, the two expansions do not work when t is too large compared to N .

An "Unstable" Malware Propagation Model

We consider a variation of the malware propagation example presented in Section 5 that is inspired by the model of [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF]. The system is composed of N nodes. Each node can be dormant (D), active (A) or susceptible (S). Let X D , X A , X S denote the proportion of dormant, active and susceptible nodes. A node that is dormant becomes active at rate 0.1 + 10X A . An active node becomes susceptible at rate 5 and a susceptible node becomes dormant at rate 1 + 10X A X D +δ , where δ is a parameter of the model. This translates into the following transitions:

(X D , X A , X S ) → (X D - 1 N , X A + 1 N , X S ) at rate N (0.1 + 10X A )X D (X D , X A , X S ) → (X D , X A - 1 N , X S + 1 N ) at rate N 5X A (X D , X A , X S ) → (X D + 1 N , X A , X S - 1 N ) at rate N (1 + 10X A X D + δ )X S
This model satisfies all the assumptions (A1-A2) needed to apply Theorem 1 that characterize the transient regime. It also satisfies (A3): There exists a unique stationary distribution because for each system size N , the stochastic model is a finite state irreducible Markov chain. This model, however, does not satisfy assumption (A4) for all possible values of the parameter δ . Indeed, there exists a parameter value δ * ≈ 0.18 such that the mean field limit has a unique attractor if and only if δ > δ * . For δ < δ * , the mean field approximation has a unique fixed point but unless the initial state is this fixed point, the limiting behavior of the solution of the ODE is an orbit. This is illustrated in Figure 5 where the two possible regimes are shown: for δ = 0.1 the system has a stable orbit and an unstable fixed point. For δ = 0.5 the system has a globally stable attractor.

It is known that when the mean field approximation has a globally stable attractor, then the sequence of stationary measures of the stochastic processes concentrates on this attractor as the system size N goes to infinity. On the other hand, when the mean field approximation has a (even unique) fixed point that is not an attractor (for example because there exist stable orbits), the sequence of stationary measures does not necessarily concentrate on this fixed point [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802.11 MAC protocol[END_REF].

When the stochastic model of size N is a finite-state irreducible continuous time Markov chain, it has a unique stationary distribution and X (N ) (t) converges in distribution to a variable X (N ) distributed according to this distribution. This shows that for any function h lim t →∞ E[h(X (N ) (t))] = E[h(X (N ) )]. Theorem 1 also shows that for any fixed time step t, lim N →∞ E[h(X (t))] = h(x(t)) where x(t) is the mean field approximation. These reasons explain why one cannot exchange the limits t → ∞ and N → ∞ : lim

N →∞ lim t →∞ E h(X (N ) (t)) lim t →∞ lim N →∞ E h(X (N ) (t)) = lim t →∞ h(x(t)),
because the limit on left hand side is independent of the initial condition of the Markov chain while the limit on the right-handside is not necessarily well defined if x(t) does not converge to a unique fixed point regardless of the initial condition.

Instability of the the Expansions

One may hope that the expansions could be able to correct the nonexchangeability of the limits or at least would be able to compensate for some of the deviation. We show in fact in Figure 6 that not only the expansions do not correct the error of the mean field approximation but they can even make it worse when the mean field approximation has a limiting cycle (case δ = 0.1).

To see that, we compare in Figure 6 the mean field approximation, the two expansions and an estimation of E[X (t)] obtained by simulation for the example described in Section 7.1 in the case where the fixed point is not an attractor (δ = 0.1). We observe that for N = 50, the mean field approximation provides an accurate approximation of E[X (t)] for t ≤ 1 and then starts oscillating for larger values of t whereas E[X (t)] stabilizes. The two expansions are slightly more accurate than the mean field approximation until t ≈ 1.2. After this time, they diverge quickly and are much less accurate than the mean field approximation. The main explanation for this fact is that when the mean field approximation does not have an attractor, the ODE of Theorem 1 are unstable and the oscillations of the constants V (t) and A(t) grow with time. Note that the larger is N , the later the mean field approximation and its expansions start diverging from the expectation estimated by simulation.

When δ = 0.5, the fixed point is an exponentially stable attractor. In this case, the error made by the mean field approximation (or by any of the two expansions) remains bounded with time, see Figure 7. Moreover in this case the expansions provide a more accurate estimate of the true value of : "Unstable" malware model : when the fixed point is not an attractor (δ = 0.1), the accuracy of the approximations is not uniform in time for a fixed system size N .

the sense that most of the mean field models studied in the queuing theory literature have a unique fixed point that is an attractor. This means that for these models it is more likely to observe a positive result like the one observed in Figure 7 rather than an oscillation like the one of Figure 6. This is no longer true when considering models from biochemistry [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF].

CONCLUSION

In this paper, we show how mean field approximation can be reby a term in 1/N and a second term 1/N 2 where N is the size of the system. We exhibit conditions that ensure that this asymptotic expansion can be applied for the transient as well as the steady-state regimes. In the transient regime, these constants satisfy ordinary differential equations that can be easily integrated numerically. We provide a few examples that show that the 1/N and 1/N 2 expansions are much more accurate than the classical mean field approximation. We also study the limitations of the approach and show that, when the mean field approximation does not have an attractor, these new approximations might be unstable for large time horizons. Obtaining a better approximation in this case remains a challenge that we leave for future work.

When we compare the accuracy of the classical mean field approximation to the one of the expansions of order 1/N and 1/N 2 , it seems that most of the gain in terms of accuracy are brought by the 1/N -term. As the 1/N 2 -term is much more expensive to compute than the 1/N term, we believe that when the 1/N 2 -expansion is too hard to compute, staying with the 1/N -expansion is already sufficient for many models. Finally, our derivation may also be exploited to obtain bounds on the error committed in the approximation of moments, which is something we aim at tackling as future work.

  the time taken by our algorithm to compute the expansions for the supermarket model described in Section 6. This figure shows the computation time as a function of the number of dimensions of the model d. It contains four panels that correspond to: (a) The time to compute V (t) and W (t) for t ∈ [0, 10]. (b) The time to compute V and W of Theorem 3. (c) The time to compute A(t), B(t), C(t) and D(t) for t ∈ [0, 10]. (d) The time to compute A, B, C and D of Theorem 3.

Figure 1 :

 1 Figure 1: Supermarket : time to compute the approximation as a function of the number of dimension d. We compare the 1/N -expansion (first line) and the 1/N 2 -expansions.

)

  Size Expansions of Mean Field Approximation: Transient and Steady-State Analysis Conference'17, July 2017, Washington, DC, USA

10 Figure 2 :

 102 Figure 2: Malware model, transient regime: comparison of the mean field approximation, the 1/N and 1/N 2 expansions and the exact value.

Figure 3 :

 3 Figure 3: Malware model, transient regime : we compare the error of the expansion of order 1/N with the constant A(t) of Theorem 1.

NN

  = 10, ρ = 0.9, k = 2. = 20, ρ = 0.9, k = 2.

Figure 4 :

 4 Figure 4: Supermarket model and transient regime: Comparison of the classical mean field approximation and the two expansions with data from simulations.

Figure 5 :

 5 Figure 5: The unstable malware model : illustration of the two possible regimes of the mean field approximation.

Figure 6

 6 Figure6: "Unstable" malware model : when the fixed point is not an attractor (δ = 0.1), the accuracy of the approximations is not uniform in time for a fixed system size N .

Figure 7 :

 7 Figure7: "Unstable" malware model in the stable case δ = 0.5 and N = 50 (complement of Figure6).

  .6.4.2 Steady-state.In Table2, we present results that illustrate the accuracy of the expansions compared to the one of the classical mean field approximation. We choose a few values of k and ρ. More

Table 2 :

 2 , we report the value of E[X i ] for various values of the parameters and i ∈ {2 . . . 7}. Note that we do not report the value E[X 1 ], Supermarket model, steady-state average queue length : comparison of the value computed by simulation with the three approximations.

	N	k	ρ Mean field1/N -expansion1/N 2 -expansionSimulation
	10	2	0.7 1.1301	1.2150		1.2191	1.2193
	20	2	0.7 1.1301	1.1726		1.1736	1.1737
	10	2	0.9 2.3527	2.7513		2.8045	2.8002
	20	2	0.9 2.3527	2.5520		2.5653	2.5662
	10	2 0.95 3.2139	4.1017		4.3265	4.2993
	20	2 0.95 3.2139	3.6578		3.7140	3.7124
	10	3	0.9 1.8251	2.2364		2.3322	2.3143
	20	3	0.9 1.8251	2.0307		2.0547	2.0517
	50	3	0.9 1.8251	1.9073		1.9112	1.9106
	100 3	0.9 1.8251	1.8662		1.8672	1.8672
	10	4 0.95 2.0771	2.9834		3.8704	3.3268
	20	4 0.95 2.0771	2.5303		2.7520	2.6376
	50	4 0.95 2.0771	2.2584		2.2939	2.2787
	100 4 0.95 2.0771	2.1678		2.1766	2.1732
				X 2	X 3	X 4	X 5	X 6	X 7
	ρ=0.9, k=2, N =10
	Mean field	0.729 0.478 0.206 0.038 0.001 0.000
	1/N -expansion 0.742 0.544 0.361 0.179 0.025 0.000
	1/N 2 -expansion 0.741 0.533 0.316 0.194 0.116 0.005
	Simulation	0.741 0.534 0.327 0.170 0.077 0.032
	ρ=0.95, k=2, N =20
	Mean field	0.857 0.698 0.463 0.204 0.039 0.001
	1/N -expansion 0.861 0.721 0.544 0.371 0.184 0.026
	1/N 2 -expansion 0.861 0.719 0.527 0.321 0.210 0.122
	Simulation	0.861 0.719 0.530 0.334 0.178 0.083
	ρ=0.9, k=4, N =10

In the remainder of the proof, we drop the dependence in t in most of the proof and write X instead of X (t ) and x instead of x (t ).


Carlos Martínez

Olivier Bernard

Francis Mairet

Maximizing microalgae productivity in a light-limited chemostat

Keywords: Dynamics and control, Industrial biotechnology, Photoinhibition, Biomass productivity, Microalgae, Photoacclimation, Chlorophyll

INTRODUCTION

Microalgae can be grown in photobioractors (PBR) illuminated by artificial light. [START_REF] Carvalho | Light requirements in microalgal photobioreactors: an overview of biophotonic aspects[END_REF] discusses the importance of providing light at the appropriate intensity to enhance growth. Low light intensities will become growth-limiting, notably in dense cultures due to selfshading. On the other hand, too high levels may lead to photoinhibition i.e. a decrease in the rate of photosynthesis due to high light intensities [START_REF] Long | Photoinhibition of photosynthesis in nature[END_REF]).

Despite the importance of the incident light, most of the mathematical works concerned with the optimization of PBRs focus only on the control of the dilution rate. See for example the works of [START_REF] Mairet | Driving species competition in a light-limited chemostat[END_REF], [START_REF] Masci | 11th ifac symposium on computer applications in biotechnology microalgal biomass surface productivity optimization based on a photobioreactor model[END_REF], or [START_REF] Grognard | Optimal strategies for biomass productivity maximization in a photobioreactor using natural light[END_REF].

Here, we determine optimal values for the dilution rate and the incident light intensity in order to maximize the biomass productivity in a continuous culture. In our approach, we follow the works of [START_REF] Gerla | Photoinhibition and the assembly of light-limited phytoplankton communities[END_REF] and [START_REF] Hsu | Dynamics of phytoplankton communities under photoinhibition[END_REF] for describing the dynamics of a microalgae culture. Then, we determine the optimal operational parameters at the steady state of the system.

This work was supported by the CONICYT doctoral grant (Carlos Martínez), and by the Phycover (ANR-14-CE04-0011) and IPL Algae in silico (INRIA) projects. 1 INRIA Sophia Antipolis, 2004, route des Lucioles BP 93, 06902 Sophia Antipolis Cedex, France (e-mail {carlos.martinez, olivier.bernard}@inria.fr) 2 francis.mairet@ifremer.fr Another route for enhancing microalgae productivity in PBRs can be to minimize light absorption by individual cells so that cells deeper inside the culture will receive enough light [START_REF] Carvalho | Light requirements in microalgal photobioreactors: an overview of biophotonic aspects[END_REF]). Light absorption is related to the chlorophyll content of microalgae. In this context, along the paper, we indicate how the chlorophyll content is related to some parameters. This allows us, in the last section, to show the existence of an optimal chlorophyll-carbon quota for maximizing productivity. This article is organized as follows. In section 2, we describe a model for light-limited cultures. In section 3, we analyze the dynamics of this model. Finally, in section 4, we present the criterion for maximizing productivity and discuss the influence of chlorophyll content on productivity.

MODELLING LIGHT-LIMITED GROWTH OF MICROALGAE.

Let us consider a perfectly mixed PBR of depth L where microalgae grow (see Figure 1). Let us assume that light is attenuated exponentially with depth according to the Lambert-Beer law i.e. at a distance z ∈ [0, L] from the illuminated surface, the corresponding light intensity

I(x, I in , z) satisfies I(x, I in , z) = I in e -axz , (1)
where x denotes the microalgae concentration, I in is the incident light intensity, and a > 0 is the specific light attenuation coefficient. According to [START_REF] Bernard | Modelling of Microalgae Culture Systems with Applications to Control and Optimization[END_REF], the extinction coefficient a is correlated to the chlorophyll quota θ by means of

= a 1 θ + a 2 , (2)
with a 1 and a 2 non-negative parameters. The chlorophyll quota represents the cellular chlorophyll-to-carbon ratio.

Let us denote by µ the specific growth rate function of the microalgae. We assume that nutrients and carbon dioxide are so that light is the single factor that limits growth i.e. µ : R + -→ R + is a function of the light intensity I perceived by the microalgae. Following [START_REF] Huisman | Light-limited growth and competition for light in well-mixed aquatic environments: An elementary model[END_REF] Eq. ( 6) shows that the AGR is determined by the incident light intensity and the light intensity at the bottom. The function g was first defined by [START_REF] Gerla | Photoinhibition and the assembly of light-limited phytoplankton communities[END_REF].

In the literature we find different models for µ, for example the model of [START_REF] Eilers | A model for the relationship between light intensity and the rate of photosynthesis in phytoplankton[END_REF] and the following reparametrization presented by [START_REF] Bernard | Validation of a simple model accounting for light and temperature effect on microalgal growth[END_REF];

µ(I) = µ max I I + µmax α I I * -1 2 , (8)
with α the initial slope of the light response curve, µ max the specific maximum growth rate, and I * the light intensity at which µ reaches µ max . Figure 2 shows the function µ given in ( 8) with kinetic parameters from Table 1. From now on, we assume that µ is given by Eq. ( 8).

According to [START_REF] Bernard | Modelling of Microalgae Culture Systems with Applications to Control and Optimization[END_REF], the parameter α can be written as

α = α 1 θ, (9)
with α 1 the initial slope of the chlorophyll-specific growth rate µ Chl := µ/θ which is independent of θ. Expression (9) shows that α is directly proportional to the pigment content, therefore lower values of α are related to more transparent microalgae.

Table 1. Kinetic parameters of µ in (8).

Parameter Value µmax 1.63 d -1 I * 87.2 µmol m -2 s -1 α 0.027µmol -1 m 2 s d -1
Let us remark that µ (and any Haldane-type model3 ) satisfies:

µ(I 1 ) = µ(I 2 ) if, and only if I 1 = I 2 or I 1 I 2 = I * 2 . ( 10
)
We will write σ(I) = I * 2 /I for any I > 0. Thus, according to (10), microalgae grow at same rate when they perceive the light intensity I or σ(I) (see Figure 2). Property ( 10) is reflected in the function g by the following property;

g(I in , σ(I in )) = g(I in , I * ), I in > I * . ( 11
)
To prove (11), by using Eq. ( 10), we can easily show that is enough to prove the following identity

I * I * 2 /Iin µ(I) I dI = Iin I * µ(I) I dI. (12)
This follows directly from doing the change of variables J = I * 2 /I on the right side.

To end this section, we give an explicit expression for g. Indeed, let

∆ := 4 µ max αI * -1 > 0 (13) then g(I in , I out ) = 2µ max / √ ∆ ln(I in /I out ) tan -1 (κ(I in )) -tan -1 (κ(I out )) , (14)
with

κ(I) = 2µmax αI * √ ∆ I I * -1 + 1 √ ∆ .

DYNAMICS OF A LIGHT-LIMITED CHEMOSTAT.

Considering a dilution rate D and a loss rate m (sum of respiration and mortality rates), the dynamic evolution of the microalgae concentration x is given by ẋ

= [g(I in , I out (x, I in )) -D -m]x. (15)
By standard arguments (i.e. ∂(xg(I in , I out (x, I in ))/∂x and g(I in , I out (x, I in )) are both bounded by µ max for all x ≥ 0) it can be shown that Eq.( 15) admits a unique global solution for any non-negative initial condition.

It is clear that µ (see Eq.

< m + D ≤ g(I in , γ(I in )). ( 16
)
Then, there exists a unique I s such that g(I in , I s ) = D + m, and

I s ≤ γ(I in ). ( 17
)
Thus, we have the positive equilibrium x s of Eq.( 15)

x s := 1 aL ln I in I s , (18)
If D + m ≤ µ(I in ), then x s is the unique positive equilibrium of Eq. ( 15) and is globally stable. If

µ(I in ) < D + m. ( 19
)
then there exist another equilibrium x u given by

x u = 1 aL ln I in I u , (20) b) a
In this situation x s attracts any solution of (15) with an initial concentration higher than x u , and x = 0 attracts any solution with an initial concentration lower than x u (see Fig. 3b)). In this case we say that the system faces bi-stability. This behavior is due to photoinhibition. If the initial biomass concentration is high enough, self-shading reduces this effect of inhibition. We note that if D + m = g(I in , γ(I in )) then x u = x s .

MAXIMIZING MICROALGAE PRODUCTIVITY.

Let us assume that condition (16) holds. We define the steady-state micro-algal biomass surface productivity as P := DLx s , with x s defined by Eq. ( 18). We will simply speak of productivity when referring to P . The productivity represents the quantity of microalgae that is produced per unit of area and time when the system reaches its steady state. Recalling Eq.( 17) we can rewrite the productivity as

P = D a ln I in I s . ( 22
)
From Eq.( 17) we note that I s is strictly increasing on m (see Figure 3 a)), therefore the productivity is decreasing on m. However, the dependence of P on D and I in is not trivial. By combining equations ( 17) and ( 22), we obtain that

P = 1 a (g(I in , I s ) -m) ln I in I s = 1 a Iin Is µ(I) -m I dI. (23)
Let Îin be such that

µ( Îin ) = m, Îin > I * (24)
then, the integral in ( 23) reaches its maximum value when

I s = σ( Îin ) and I in = Îin , (25)
because in that way the integral is calculated over the maximal interval where the function to be integrated is non-negative. In the following theorem we show that by taking

D = D with D = g( Îin , σ( Îin )) -m, (26)
I s satisfies Eq.( 25).

Theorem 1. The productivity is maximal when D = D and I in = Îin with D and Îin satisfying ( 24) and ( 26). The maximal productivity is given by

P max := 1 a Îin σ( Îin) µ(I) -m I dI. (27)
Under these conditions, Eq.( 15) faces bi-stability; the solutions of Eq.( 15) reach the stable equilibrium if and only if the initial microalgae concentration x 0 satisfies

x 0 > xu with xu = 1 aL ln I in I * . (28)
Proof. Let Îin > I * defined by Eq.( 24) and D defined by Eq.( 26). First, we note that D > 0, indeed

D = 1 ln Îin σ( Îin) Îin σ( Îin) µ(I) -m I dI > 0. ( 29
)
From the definition of D and Eq.( 11), we have that D + m = g( Îin , σ( Îin )) = g( Îin , I * ).

Since σ( Îin ) < I * , we conclude that I s = σ( Îin ) and that Eq.( 15) faces bi-stability with I u = I * . Finally, by using Eq.( 20) we conclude that the unstable equilibrium is given by Eq.( 28).

Theorem 1 not only gives the maximal productivity but an expression for the minimal initial microalgae concentration xu to avoid the washout. Eq.( 28) shows that xu increases by reducing the depth of the reactor. Thus, even if the productivity does not depend on the depth of the reactor, . Plot of the optimal incident light intensity, the optimal dilution rate, and the maximal productivity as functions of the mortality rate. The function µ is given by ( 8) with kinetic parameters of C. vulgaris given in Table 1, and a = 0.2m 2 g -1 .

the depth should be chosen such that xu is not too big. The condition x 0 > xu can be rewritten as I * > I out,0 := I in e -ax0L . Thus, the result related to bi-stability can be stated as follows: the system will washout only if at the beginning the light intensity at the bottom of the reactor is higher than I * (i.e. if all microalgae in the PBR are suffering from photoinhibition).

For evaluating the optimal parameters and the maximal productivity we take the kinetic parameters from Table 1. Figure 4 shows that for small values of m the optimal incident light intensity takes high values; if m = 0.1d -1 then Îin = 2151 µmol m -2 s -1 . In fact, as m approaches to zero, the optimal incident light intensity goes to infinity.

To see how Îin varies with other parameters, we note that it is the biggest solution of the following equation for I

(I -I * ) 2 I * 2 = α 1 m - 1 µ max I. (30)
The solutions of this equation correspond to the intercepts of a parabola (left side of Eq.( 30)) and a line (right side of Eq.( 30)). Figure 5 shows this situation. Thus, any increase and the line

y = α 1 m -1 µmax I.
in the slope of the line results in an increase of the value of Îin . In particular, this occurs when m decreases or when α or µ max increases. On the other hand, any increase on I * will open wider and translate to the right the parabola. Thus, any increase on I * increases the value of Îin .

Fig. 4 shows that the maximal productivity approaches a finite value, that we denote T , when m approaches to zero. The following proposition gives a simple expression for evaluating T . Proposition 2. The maximal productivity is bounded from above by

T = 2µ max a √ ∆ π 2 -arctan 1 √ ∆ 1 - 2µ max αI * , (31)
when ∆ > 0.

Proof. By taking derivative with respect to m in Eq.( 30) we can obtain

∂ Îin ∂m = -α Îin m 2 Îin-I * Î * 2 in 1 + I * Îin . (32)
By applying L'Hopital's rule and using Eq.( 32) we conclude that

lim m→0 + m ln Îin σ( Îin ) = 0. (33)
Now, by using Eq.( 27) and noting that Îin → ∞ and σ( Îin ) → 0 as m → 0 + , we obtain that

T := lim m→0 + P max = 1 a ∞ 0 µ(I) I dI. (34)
Finally, by using Eq.( 14) to evaluate the integral in Eq.( 34) we obtain the result of the Proposition.

In the example of Fig. 4, we have that T = 22.35g m -2 d -1 . Thus, according to the parameters in this paper, C.

= 1.44d -1 µmol -1 m 2 gC gChl -1 , a 1 = m 2 gChl -1 , m = 0.1 d -1
, and I * , µ max from Table 1.

vulgaris could never reach a higher productivity than 22.35g m -2 d -1 in a PBR operated at constant dilution rate and constant incident light.

With respect to the chlorophyll content in microalgae cells, according to Eq.( 9), α decreases with θ, and consequently Îin becomes closer to I * (see Fig. ( 5)). But how does this affect the productivity? Let us assume that µ max and I * do not vary with θ and that a and α vary according to equations ( 2) and ( 9) respectively. Fig. 6 shows how the maximal productivity P max varies with θ for different values of a 2 . We can see that if a 2 > 0 (i.e. if light is not only absorbed by chlorophyll), then there exists an optimal value for θ. However, if a 2 = 0 (i.e. light is absorbed only by chlorophyll), then the more transparent the microalgae, the higher the maximal productivity.

Experimental results of [START_REF] Anning | Photoacclimation in the marine diatom skeletonema costatum[END_REF] and [START_REF] Neidhardt | Photosystem-ii repair and chloroplast recovery from irradiance stress: relationship between chronic photoinhibition, light-harvesting chlorophyll antenna size and photosynthetic productivity in dunaliella salina (green algae)[END_REF] show that µ max and I * can also vary with θ. Let us assume that all the kinetic parameters (α, µ max , and I * ) and m depend on θ and let us define

M (θ) := α(θ)I * (θ) 1 m(θ) - 1 µ max (θ) . ( 35
)
The following theorem gives conditions over M such that for too high and too low values of θ the maximal productivity behaves as in Fig. 6. In particular, it states that for enhancing the maximal productivity when a 2 > 0, it is not convenient to reduce too much the chlorophyll quota in microalgae cells. Theorem 3. Assume that a varies according to Eq.( 2), that there exists N > 0 such that 0 < µ max (θ) -m(θ) ≤ N for all θ > 0, and that M (θ) (defined in Eq.( 35)) is a differentiable function of θ. We have that

(1) If lim θ→∞ M (θ) = ∞ and lim θ→∞ M (θ) < +∞, then lim θ→∞ P max = 0. (2) If lim θ→0 + M (θ) = 0 and lim θ→0 + M (θ) < +∞, then lim

θ→0 + P max = ∞ if a 2 = 0, 0 if a 2 > 0. ( 36
)
Proof. Since µ(I) ≤ µ max for all I ≥ 0, from Eq.( 27) and the definition of σ(I), we can easily obtain that

P max ≤ 2 a 1 θ + a 2 N ln Îin I * . (37)
Let J := Îin /I * . Then, from Eq.( 30), we have that

(J -1) 2 = M (θ)J. ( 38
)
From where we obtain that

∂J ∂θ = M (θ)J (J -1) 1 + 1 J . ( 39
)
From Eq.( 37) and the definition of J, we obtain that lim θ→∞ P max ≤ 2N lim θ→∞ ln(J)/(a 1 θ + a 2 ). After applying L'Hopital's rule and using Eq.( 39), we obtain that

lim θ→∞ P max ≤ lim θ→∞ M (θ) α 1 (J -1) 1 + 1 J . ( 40
)
Since M (θ) → ∞ as θ → ∞, from Eq.( 38) we conclude that J → ∞ as θ → ∞. Now, since lim θ→∞ M (θ) < +∞, we conclude that the right side of Eq.( 40) converges towards 0, and consequently lim θ→∞ P max = 0. Thus a) is proved. The part b) follows the same arguments.

CONCLUSIONS AND FUTURE WORK.

We determined the maximal microalgal productivity that can be reached in a light-limited chemostat operated at constant dilution rate and constant incident light intensity. A simple criterion is obtained: the biomass productivity is maximal when the incident light intensity and the light intensity at the bottom of the culture are such that the growth rates at the top and at the bottom equal the loss rate. Under these optimal conditions the system faces bistability i.e. the system washouts when the density is below a certain threshold. We provided a simple expression for determining this threshold. We also studied the effects of the chlorophyll content on the maximal productivity. We showed, that if light is not only absorbed by chlorophyll, then there is an intermediate value of the chlorophyll quota maximizing the productivity.

As a future work, we will deal with the optimization of productivity by considering the background turbidity of the medium (i.e. by adding terms to Eq.( 5)) and (or) by including a model of photo-acclimation. By adding turbidity the productivity becomes dependent of the depth of the reactor, which makes more tricky the choice of the depth. By adding photo-acclimation, in long term, the pigment content become dependent of the incident light which could change the compensation condition (25) for maximizing the productivity.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Scheme of the light-limited chemostat.

Following

  , we compute the average growth rate function (AGR) in the reactor, denoted μ, by integrating the local growth rate over all the reactor μGerla et al. (2011), we do the change of variables I = I(x, I in , z). From Eq. (1) dI/dz = -axI. So if x > 0, we can rewrite the AGR μwhere I out (x, I in ) is used for indicating the light intensity at the bottom of the reactor i.e. I out (x, I in ) = I in e -axL . If I in > 0, note from (1) that axL = ln(I in ) -ln(I out (x, I in )). (5) By substituting (5) in (4) we have μ(•) = g(I in , I out (x, I in )), (6) with g(I in , I out ) = 1 ln(I in ) -ln(I out )

  (8)) is differentiable, µ (I) > 0 for all I ∈ [0, I * ), µ (I) < 0 for all I ∈ (I * , ∞), µ(0) = 0 and lim I→∞ µ(I) = 0. Let us assume that I in > I * . Following the work of[START_REF] Hsu | Dynamics of phytoplankton communities under photoinhibition[END_REF], we know that g(I in , 0) = 0, g(I in , I in ) = µ(I in ) and that there exists γ(I in ) < I * such that the function g(I in , •) is strictly increasing on [0, γ(I in )] and strictly decreasing on [γ(I in ), I in ] (see Figure3 a)). Thus, if m + D > g(I in , γ(I in )), then the equation m + D = g(I in , I out (x, I in )) for x has no solution (see Figure3a)) and x s = 0 is the unique equilibrium of Eq. (15) which is globally stable. Let us assume that 0

)

  Fig. 3. a) Function g and its intersection with D + m. b) Graphical representation of bi-stability. The black lines represent the solution of Eq. (15) for different initial conditions. with I u satisfying g(I in , I u ) = D + m, and γ(I in ) ≤ I u . (21)

  Fig. 4. Plot of the optimal incident light intensity, the optimal dilution rate, and the maximal productivity as functions of the mortality rate. The function µ is given by (8) with kinetic parameters of C. vulgaris given in Table1, and a = 0.2m 2 g -1 .

Fig. 5 .

 5 Fig. 5. Plots of the parabola y = (I-I * ) 2 I * 2

Fig. 6 .

 6 Fig.6. Maximal productivity as a function of the chlorophyll quota for different values of a 2 . The parameters are taken to be α1 = 1.44d -1 µmol -1 m 2 gC gChl -1 , a 1 = m 2 gChl -1 , m = 0.1 d -1, and I * , µ max from Table1.

By a Haldane-type model, we mean µ(I) = I aI 2 +bI+c with a, c > 0.