Matthieu Jonckheere
email: matthieu.jonckheere@gmail.com

Balakrishna J Prabhu
email: balakrishna.prabhu@laas.fr

Asymptotics of Insensitive Load Balancing and Blocking Phases *

Keywords: Insensitive load balancing, blocking phases, mean-field scalings, QED-Jagerman-Halfin-Whitt regime

Load balancing with various types of load information has become a key component of modern communication and information systems. In many systems, characterizing precisely the blocking probability allows to establish a performance trade-off between delay and losses. We address here the problem of giving robust performance bounds based on the study of the asymptotic behavior of the insensitive load balancing schemes when the number of servers and the load scales jointly. These schemes have the desirable property that the stationary distribution of the resulting stochastic network depends on the distribution of job sizes only through its mean. It was shown that they give good estimates of performance indicators for systems with finite buffers, generalizing henceforth Erlang's formula whereas optimal policies are already theoretically and computationally out of reach for networks of moderate size.

We study a single class of traffic acting on a symmetric set of processor sharing queues with finite buffers and we consider the case where the load scales with the number of servers. We characterize the response of symmetric systems under those schemes at different scales and show that three amplitudes of deviations can be identified according to whether ρ < 1, ρ = 1, and ρ > 1. A central limit scaling takes place for a sub-critical load; for ρ = 1, the number of free servers scales like n θ θ+1 (θ being the buffer depth and n being the number of servers) and is of order 1 for super-critical loads. This further implies the existence of different phases for the blocking probability. Before a (refined) critical load ρ c (n) = 1 -an -θ θ+1 , the blocking is exponentially small and becomes of order n -θ θ+1 at ρ c (n). This generalizes the well-known Quality and Efficiency Driven (QED) regime or Halfin-Whitt regime for a one-dimensional queue, and leads to a generalized staffing rule for a given target blocking probability.

Introduction

Load balancing is a critical component in multi-servers systems such as call centers, server farms, as well as in distributed systems with many applications running on different servers (as a single example, see the load balancing needs of the CERN network [1]). Despite its intensive use, there are few efficient rules of thumb for dynamic load balancing schemes i.e. when the decision of the dispatcher depends on the instantaneous load (e.g. number of jobs) at each server, for which performance evaluation is within reach using currently known techniques. This fact becomes even more evident for large systems with asymmetric server speeds and blocking, where both the precise structure of optimal policies (for specific traffic descriptions) and their performance elude current knowledge and techniques.

Two types of ideas were largely employed to overcome this difficulty: considering large scale networks and obtain asymptotic results using propagation of chaos on the one hand; restricting the load balancing schemes to obtain more tractable reversible processes on the other hand. We aim here at combining both techniques and show that it may lead to very precise results which are typically out of reach with other techniques. In turn, these results allow to give universal (i.e. valid for all job-size distributions) lower bound of performance. In particular, different asymptotic scalings for the blocking probability can be identified: an exponentially small blocking probability for sub-critical loads; a polynomial order in the critical regime; and a constant level for super-critical loads. Before describing our contribution more precisely, let us recall the vast effort of research in the two mentioned directions.

A first way of overcoming the computing difficulty of load balancing problems for processor sharing systems (and more generally symmetric queues) with generic job-size distribution is by restricting the routing policies so that the stationary regime of the system becomes invariant to the job-size distribution (except for its mean), leading to the insensitive load balancing (see [START_REF] Bonald | Insensitive bandwidth sharing in data networks[END_REF] for more details). To understand the underlying principles, it is useful to come back to the very properties of the Erlang formula which was clearly a revolution for performance evaluation of telephone networks and arguably the true start of queuing theory. The Erlang formula, which gives the probability of loss for a set of telephone lines, bases its lasting success on simplicity and robustness: 1. the only assumptions which are required to apply the formula are Poisson arrivals and independent calls durations;

2. the formula is insensitive to the call duration distribution and depends on a unique parameter: the traffic intensity; and 3. it can efficiently be computed using a recursive formula.

At the mathematical level, the key property is the reversibility of the birth-and-death processes modeling the system under Markovian assumptions which implies the insensitivity property: the stationary measure of the system does not depend on the whole call distribution but only on its mean.

Those same principles were translated to performance models of best effort and voice traffic in [START_REF] Bonald | Insensitive bandwidth sharing in data networks[END_REF][START_REF] Bonald | A queueing analysis of max-min fairness, proportional fairness and balanced fairness[END_REF] for models with balanced allocations and in [START_REF] Bonald | Insensitive bandwidth sharing in data networks[END_REF][START_REF] Bonald | Insensitive load balancing[END_REF] for models with dynamic insensitive load balancing. For multi-class networks with insensitive load balancing, Markov Decision Programming techniques were employed in [START_REF] Leino | Insensitive load balancing in data networks[END_REF], structural results were provided in [START_REF] Jonckheere | Towards an Erlang formula for multiclass networks[END_REF] while extensive simulations were proposed in [START_REF] Pla | Optimal robust policies for bandwidth allocation and admission control in wireless networks[END_REF]. All this research dealt with networks of fixed size and allowed to draw the following conclusions. For networks with a unique class of traffic, the insensitive load balancing compares very accurately to optimal policies for a given job-size distribution, while delay estimation are a bit less accurate [START_REF] Bonald | Insensitive bandwidth sharing in data networks[END_REF][START_REF] Bonald | Insensitive load balancing[END_REF]. The penalization imposed by reversibility is greater for multi-class networks while the sensitivity (of optimal sensitive policies) also deteriorates [START_REF] Leino | Insensitive load balancing in data networks[END_REF][START_REF] Jonckheere | Towards an Erlang formula for multiclass networks[END_REF]. Hence a small to moderate price has to be paid for robustness and simplicity. It is perhaps counter-intuitive to notice that for models with infinite buffers, this price becomes very high. It was indeed proved that if the state space is infinite and in the absence of blocking, the optimal insensitive load balancing (for any reasonable criterion) is static (i.e., does not depend on the queue-lengths) and is hence much less efficient than a state-dependent sensitive load balancing [START_REF] Jonckheere | Insensitive versus efficient dynamic load balancing in networks without blocking[END_REF]. For sensitive schemes like join-the-shortest-queue, there is no characterization of the stationary measures in a general setting (see, for example, the approximations in [START_REF] Knessl | On the finite capacity shortest queue problem[END_REF]).

On the other side of the spectrum, a specific attention has been given in the last decades to mean-field type results for different type of networks with load balancing applications. In particular, a great deal of research, which started with the seminal works of [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An assymptotic approach[END_REF], has been dedicated to prove mean-fields limits for schemes like join-the-shortest-of-d (JSQ(d)) among n queues where n is large. For example, transient functional law of large numbers and propagation of chaos were obtained in [START_REF] Graham | Chaoticity on path space for a queueing network with selection of the shortest queue among several[END_REF] for FIFO scheduling, [START_REF] Stolyar | Pull-based load distribution among heterogeneous parallel servers: the case of multiple routers[END_REF] obtains mean-field limit for the join-the-idle queue (JIQ) while [START_REF] Mukherjee | Universality of Load Balancing Schemes on Diffusion Scale[END_REF] computes the diffusive limit in the Halfin-Whitt regime for a class of policies of which JIQ and JSQ(d) policies are a special case. For general service time distributions, propagation of chaos properties and asymptotic behaviour of the number of occupied servers were obtained for the JSQ policy in [START_REF] Bramson | Asymptotic independence of queues under randomized load balancing[END_REF]. All these results concern systems without blocking and are sensitive to the job-size distribution. For systems with blocking, the recent work of [START_REF] Mukhopadhyay | Mean field and propagation of chaos in multi-class heterogeneous loss models[END_REF] computes the mean-field limit for the JSQ(d) scheme and exponentially distributed job sizes.

It is hence natural and complementary to look at insensitive networks with a very large number of servers and given buffer depth in order to see if the results obtained for finite networks scale appropriately. As detailed below, this leads to very precise results which are qualitatively different from the case without blocking and are out of reach for sensitive policies with blocking. This in turn provides simple dimensioning rules.

Contributions

We study the asymptotics of a set of n processor sharing servers, each with buffer size θ, fed by a Poisson process of intensity ρn when n gets large, under the family of insensitive load balancing schemes shown to be optimal (in the class of insensitive load balancing) in [START_REF] Bonald | Insensitive load balancing[END_REF]. A more detailed description of the model is given in the next section.

Building on closed-form expressions for the stationary measure, we characterize precisely the asymptotics of the stationary measure and the blocking probability for various scalings of the load. Consequently, we provide universal benchmarks for achievable performance which have no known counterpart for sensitive policies.

We first obtain the stationary measure of the number of occupied servers and give its transient mean-field limit. Considering the symmetric version of the model, we show that the functional law of large numbers also holds for the stationary version of the system (limits in n and t commute). The existence and uniqueness of the limiting stationary probabilities are proved through a monotonicity argument involving the Erlang formula, while the stationary point is characterized through the Erlang formula. This implies simple conclusions on the asymptotic behavior of the blocking probability: the blocking is asymptotically vanishing for the sub-critical (ρ < 1) case and is equal to 1 -ρ -1 for the super-critical case ρ > 1. In both cases, this blocking probability corresponds to the optimal blocking probability achievable by any non anticipating policy. Of course this is far from being sufficiently informative and we are led to focus on a more detailed study of the stationary distribution for large n, establishing both large deviations principles for sub-and super-critical cases and moderate deviations results. We show that, when ρ < 1 is fixed, the blocking probability is exponentially small, and we characterize the most probable deviations from the mean-field limit. The large deviation cost is shown to be a sum of two terms: the "distance" to the stationary point from distributions with a given mean plus the cost of having a different mean from the true stationary mean. We also show that a central limit theorem is valid for the occupation numbers around the stationary point of the mean-field in the sub-critical regime. For the critical case ρ = 1, the right scaling is not anymore of order √ n. Using local limit theorems and exploiting the characterization of deviations from the mean-field limits, we show that the number of free servers scales like n θ θ+1 , the limiting distribution depending on θ and coinciding with the normal distribution only for θ = 1. In a third step, we study the critical case at a finer scale and show that a qualitative phase transitions occurs at the critical load ρ c (n) = 1 -an -θ θ+1 where θ is the buffer depth. The blocking probability is exponentially small until ρ c (n) and of order n -θ θ+1 at this critical load. This generalizes the Halfin-Whitt regime established for the M/M/n/n system, (in that case, the correct scaling for the moderate deviations stayed of order √ n), and show that the popular staffing rule established for the M/M/n/n system does actually change with the value of θ when load balancing is employed. The super-critical regime is simpler to characterize, the deviations being of order 1. We illustrate these findings on simple numerical experiments. Finally, we comment on how these results can be used for performance planning, for instance in trading delay for blocking, while controlling the level of blocking fixing the number of servers and how the level of information needed for a possible implementation can be reduced. We also give insights on possible future work.

2 Review of the optimal insensitive load balancing policy Notation. We use the following notations common to all sections. For any vector space (the exact one under consideration will be clear from the context), let e i be the point defined by (e i ) i = 1, (e i ) j = 0, j = i. We denote:

|x| = k i=1 x i and |x| x = |x|! x 1 ! . . . x k ! .
We denote by 1 S the indicator function of S, that is the map taking value 1 inside S and 0 outside and denote respectively by R + and R * + the set of non-negative and positive reals. This section is a review of the relevant definitions, merits and results known for the insenstive load-balancing policy investigated in this paper. The narrative here is for a more general model than the one we shall analyse. Nonetheless, it gives a flavour of the possible generalizations, some of which are elaborated upon in Section 6.

Consider a dispatcher and a set of n processor sharing servers with speed µ i for i = 1 . . . , n. Jobs with i.i.d. sizes sampled from a generic distribution of mean 1 arrive to the dispatcher according to a Poisson process of intensity λ. The dispatcher routes an incoming job to one of the servers according to the following insensitive load balancing rule. Let θ = (θ 1 , . . . , θ n ) a vector of natural numbers and X ⊂ N n a finite coordinate convex set describing the constraints on the number of jobs in each server (for instance X = {x : x i ≤ θ i , ∀i = 1 . . . n}). Then, the incoming job is routed to server i with probability

a θ i (x) = θ i -x i N j=1 (θ j -x j ) 1 x+ei∈X . (1)
Note that with this rule, the number of jobs in server i is smaller than θ i for all i = 1 . . . n.

One can view θ i as the buffer size of server i but it could be a smaller number chosen to guarantee a certain rate of service. Also note that this rule depends on the speeds µ i only through the vector θ. Nevertheless this load balancing rule was proved to be optimal 1 in the set of insensitive load balancing (for a unique class of traffic) in [START_REF] Bonald | Insensitive load balancing[END_REF], i.e., given the speeds µ i there exists an optimal vector θ such that this rule is optimal among all insensitive load balancing.

Let X be the stochastic process valued in X describing the number of on going jobs in each server. Under Poisson arrivals and exponentially distributed job sizes, X is a continuous-time jump Markov process, on the state space X , with infinitesimal generator Q = (q(x, y)) x,y∈X given by ∀x ∈ X ,

     q(x, x -e i ) = µ i if x -e i ∈ X q(x, x + e i ) = λa i (x) if x + e i ∈ X q(x, y) = 0 if y ∈ X , y = x -e i , x + e i . (2)
We recall that the family of insensitive load balancing corresponds to the routing rates λ i () = λa(•) such that there exists a balance function Λ :

X → R * + , ∀i, ∀x ∈ X , x + e i ∈ X , λ i (x) = Λ(x + e i )/Λ(x), (3)
which is equivalent to the detailed balance criterion. The relationship between this criterion and insensitivity was first formulated in [START_REF] Whittle | Partial balance and insensitivity[END_REF]. Under condition (3), the process X is reversible and the stationary distribution is given by

π(x) = Λ(x)Φ(x) y∈X Φ(y)Λ(y) , (4)
with

Φ(x) = n i=1 µ -xi i .
For the optimal insensitive load balancing corresponding to the mentionned rates λa i (•), the routing balance function Λ takes the form

Λ(x) = Λ θ (x) = |θ -x| θ -x λ |x| , (5)
where |θ-x| θ-x = |θ-x|! n i=1 (θi-xi)! are the multinomial coefficients. The blocking probability, B θ , of an arriving job can be determined using the PASTA property to be π(θ).

Model and preliminary results

In the rest of the paper, we shall assume that the servers are homogeneous, that is, they have the same speed and the same buffer size. (We shall comment later on the possibility of extending those results). Without loss of generality, let the common speed be 1. The common buffer size will be taken to be θ. (From now on, θ is a natural number and not a vector as in the previous section). For the asymptotic analysis we have in mind, it turns out to be more convenient to define the state as the number of servers processing a certain number of jobs instead of the number of jobs being processed in every server. Let S = {s ∈ {0, 1, . . . , n} θ+1 : θ i=0 s i = n} be the set of states where s i corresponds to the number of servers with i jobs. In state s ∈ S, the insensitive load balancing rule described in Section 2 will route an incoming job to a server with i jobs at rate

λ i (s) = λ (θ -i)s i nθ - s , (6)
where s = θ i=0 is i . Let {S (n) (t) ∈ S} t≥0 be a stochatic process denoting, at time t, the number of servers with i jobs, i = 0, . . . , θ. Under Poisson arrivals and exponentially distributed job sizes, S (n) (t) is a continuous-time jump Markov process on the state space S with the following transition rates

S (n) (t) → S (n) (t) + e i -e i-1 at rate λ i-1 (s), i ≥ 1; S (n) (t) + e i -e i+1 at rate s i+1 , (7)
assuming that the transitions take the process to a state within S.

The reversibility property of X is preserved by S. More precisely, Theorem 1. If the job-size distribution is exponential, the process S (n) (t) is a reversible Markov process and its stationary distribution is given by

π (n) (s) = π (n) 0 (nθ -s)! (nθ)! n s θ k=0 θ! (θ -k)! (nρ) k s k , ( 8
)
where is the total number of jobs in the system, and ρ = λ/n is the load per server, and π

(n) 0

corresponds to the probability of the state with all servers empty, that is, s = 0 and s = (n, 0, . . . , 0).

Proof.

A sufficient condition for a probability measure to be the stationary measure of a Markov chain is that it satisfy the local balance equations. Consider two states s and s + e i -e i-1 , both within S. From (8),

π (n) (s + e i -e i-1 ) π (n) (s) = λ(θ -(i -1))s i-1 nθ - s 1 (s i + 1)µ , (9)
= λ i-1 (s) (s i + 1) , (10)
which are in the same proportion as the local transition rates between these two states as computed from [START_REF] Bramson | Asymptotic independence of queues under randomized load balancing[END_REF].

Corollary 1. Using the PASTA property, the blocking probability is given by

B (n) θ = π (n) 0 (nρ) nθ (θ!) n (nθ)! . (11)
Instead of using π

(n) 0
as the normalizing constant, we can resort to B (n) for this purpose, and rewrite [START_REF] Cover | Elements of Information Theory[END_REF] as

π (n) (s) = B (n) θ (nθ -s)! n s θ k=0 1 (θ -k)! (nρ) k-θ s k , (12)
A special case that will reappear throughout this paper is the one with θ = 1, which corresponds to the classical M/M/n/n queue or the Erlang loss system. Upon setting θ = 1 in (8), we obtain:

π (n) (s 0 ) = (nρ) (n-s0) (n -s 0 )! π(0), (13)
where π

(n) 0 = k≤n (nρ) n-k (n -k)! . ( 14
)
where s 0 is the number of empty servers. We hence retrieve the formula corresponding to the M/M/n/n queue, as expected.

We remind the reader that for the JSQ policy, the stationary measure is quite intricate to compute even for the case of two servers [START_REF] Knessl | On the finite capacity shortest queue problem[END_REF], making it difficult to predict the performance of this policy.

Once the stationary measure is determined, the stationary performance measures such as the mean soujourn time and the blocking probability can be numerically computed for any given set of parameters such as n, θ, or ρ. However, for large n or ρ close to 1, the relationship between the performance measures and the parameters can be obtained in a more palatable (and exploitable) form using asymptotic analysis. The following sections will follow this path leading to a mean-field limit as well as the characterization of the large and moderate deviations.

A mean-field deterministic limit

In this section, we give the limiting transient and stationary behaviour in case of exponentially distributed job-sizes when n diverges. This limit called the mean-field limit has become a classical asymptotic regime for the analysis of large queuing systems and particle systems with a large number of servers or particles [START_REF] Boudec | The stationary behaviour of fluid limits of reversible processes is concnetrated on stationary points[END_REF][START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF][START_REF] Mukhopadhyay | Mean field and propagation of chaos in multi-class heterogeneous loss models[END_REF]. In load-balancing applications, this type of analysis has been used for several policies whose stationary measure is either unknown or known for relatively small values of the number of servers (for example, shorter of d choices [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF], joining the shortest queue [START_REF] Mukhopadhyay | Mean field and propagation of chaos in multi-class heterogeneous loss models[END_REF]). Indeed, for data-centers that can have hundreds to thousands of servers, the mean-field limit can give a first-order approximation to the system behaviour both in the transient and in the stationary phase.

In the mean-field limit, dynamics for the fraction of servers containing a certain number of jobs are as follows.

Theorem 2. Fix a ρ and a θ ≥ 1. For exponentially distributed job-sizes, for all fixed time,

S (n) (t)/n L 2
--→ y(t), which is the solution of the following set of differential equations:

dy j (t) dt = ρ θ -(j -1) θ -k ky k (t) y j-1 (t) + y j+1 (t) (15)
-ρ θ -j θ -k ky k (t) y j (t) -y j (t), 0 < j < θ, (16)
dy θ (t) dt = ρ 1 θ -k ky k (t) y θ-1 (t) -y θ (t), ( 17
)
dy 0 (t) dt = y 1 (t) -ρ θ θ -k ky k (t) y 0 (t). ( 18
)
with y(0) = lim n→∞

S (n) (0) n .
Proof. For Poisson arrivals and exponentially distributed job-sizes, the technical difficulty is much less compared to that for generic job-size distribution. We hence do not give the details of the proof but only sketch the argument. Fix a time interval [0, t]. The process is trivially tight and assuming exponential job-size makes it Markov. Dynkins formula allows to write this Markov process as a drift part plus a martingale and calculating the increasing process of the martingale, one proves that the martingale part goes to 0. As a consequence, the process converges along subsequences to a deterministic process in L 2 . It remains to prove that the limit is unique which is easy in this case since using the regularity of the rates of the process, the limit is a differential equation with a Lipschitz drift. 2

Remark 1. We expect the results to hold under generic job-sizes distribution but the proof becomes much more technical as one has to work with measure-valued processes and falls out of the scope of this paper. Results like asymptotic independence for randomized load balancing schemes such as join the shortest of d queues with generic job-sizes have been proved in [START_REF] Bramson | Asymptotic independence of queues under randomized load balancing[END_REF].

Theorem 3. For 0 < ρ ≤ 1, the unique steady-state solution of the system of equations ( 15)-( 18) is given by

pj = θ - ĉ ρ θ-j 1 (θ -j)! pθ , (19)
with pθ = 1

θ k=0 θ-ĉ ρ k 1 k! . ( 20
)
where

ĉ = θ -ρErl -1 θ (1 -ρ), (21)
with Erl -1 θ as the inverse function of the Erlang blocking viewed as a function of the traffic intensity for a fixed buffer depth θ.

If ρ > 1, the unique solution is ĉ = θ, pj = 0, for j ≤ θ -1 and pθ = 1 .

Proof. Suppose first ρ < 1. It can be easily verified that, with

ĉ = θ k=0 k pk , (22)
(19) and ( 20) is the steady-state solution of ( 15)- [START_REF] Leino | Insensitive load balancing in data networks[END_REF]. We now show that ĉ as defined in ( 22) verfies [START_REF] Mukherjee | Universality of Load Balancing Schemes on Diffusion Scale[END_REF].

After some simple algebraic manipulations, it can be verified that the fixed point equation ( 22) is equivalent to the equation

(1 -ρ) θ k=0 θ -x ρ k 1 k! = θ -x ρ θ 1 θ! (23)
in the set [0, θ]. Thus solving [START_REF] Mukhopadhyay | Mean field and propagation of chaos in multi-class heterogeneous loss models[END_REF] boils down to finding a traffic intensity a := θ-x ρ such that the Erlang blocking formula with intensity a gives 1 -ρ, i.e.:

Erl θ (a) = a θ 1 θ! θ k=0 a k 1 k! = (1 -ρ).
By a simple sample path argument, the Erlang formula is an increasing function of a that is 0 in 0 and 1 in +∞. Hence, it is invertible and there is a unique a > 0 such that Erl θ (a) = 1-ρ. Now observe that by a conservation argument (the traffic entering vs traffic outgoing) we have that for all a:

a(1 -Erl θ (a)) ≤ θ, which boils down (given the definition of a) to a ≤ θ ρ .

which in turns gives a unique solution to

θ -x ρ = a. Now consider ρ > 1.
It is straightforward to verify that ĉ = θ, pj = 0, for j ≤ θ -1 and pθ = 1 is a solution. If ĉ < θ, then the drift of the differential equation ( 19) cannot be 0 which implies that the solution given is unique.

Using the generic method developed in [START_REF] Boudec | The stationary behaviour of fluid limits of reversible processes is concnetrated on stationary points[END_REF] to inverse limits in n and t under the assumption of reversibility of the Markov process under study (which is indeed verified here), we can state that Proposition 1. For ρ < 1, π (n) converges point wise to p when n and t converge to infinity.

Proof. [START_REF] Boudec | The stationary behaviour of fluid limits of reversible processes is concnetrated on stationary points[END_REF] allows to state that if there is convergence for fixed time intervals, if the process is reversible for fixed n, and if there exists a unique stationary limiting point, (which we have verified) in the previous Proposition, then the conclusion of the Proposition hold. 2

Remark 2. By insensitivity, taking limit in time first define a sequence of limiting distribution π (n) which do not depend on the specific job-size distribution and which converges towards p.

Performance consequences

Let B θ denote the stationary blocking probability in the mean-field limit, that is, when n → ∞ and t → ∞. Using the PASTA property for fixed n, the blocking probability of a job is the probability that it finds all the servers in their blocking state i.e., that upon arrival all the servers have θ tasks. Before deriving B θ , we first give a lower bound on the blocking probability that could be achieved by any non-anticipating and size-unaware load balancing policy. Proposition 2. For θ > 0, the blocking probability of any non-anticipating and size-unaware load balancing policy is greater than max(0, 1 -ρ -1 ).

Proof. We give the argument for ρ > 1. The argument for ρ ≤ 1 is similar. Consider the system in which the resources are pooled, that is, there is one server of service rate of nµ and buffer size nθ. By a path-wise argument for Markovian versions of the systems (implying by insensitivity the result for all service times in stationary regime), the blocking probability of this system will be less than any system with a set of disjoints servers. In the pooled system, the scaled number of tasks in the system X(t) n will follow the differential equation:

ẋ(t) = ρ -1, 0 < x(t) < θ, (24)
ρ(1 -Bθ ) = 1, x(t) = θ. (25)
When x(t) is in the interior of the state space, all tasks will be accepted. On the boundary x(t) = θ, the tasks which cause overflow will be blocked. Hence, the blocking probability will be B θ = 1 -ρ -1 .

We now be shown that the insensitive load-balancing policy achieves this lower bound which is independent of θ. Proposition 3. The limiting blocking probability of the insensitive load balancing policy is given by

B θ = 0 if ρ < 1; 1 -ρ -1 otherwise. ( 26
)
Proof. For ρ < 1, it can be seen that pθ < 1, and therefore, B θ = 0. For ρ ≥ 1, we shall prove first this result for θ = 1. For θ = 1, from [START_REF] Mukhopadhyay | Mean field and propagation of chaos in multi-class heterogeneous loss models[END_REF],

c = p 1 ,
and from [START_REF] Lu | Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services[END_REF],

p 1 = ρ(1 -B 1 ) 1 -p 1 (1 -p 1 ) = ρ(1 -B 1 ). Since p 1 ≤ 1, B 1 is such that ρ(1 -B 1 ) = 1. That is, B 1 = 1 -ρ -1 .
Using Proposition 3 in [START_REF] Bonald | Insensitive load balancing[END_REF], B θ ≤ B 1 for all θ ≥ 1. Now, using the lower bound in proposition 2, this implies that B θ = B 1 .

The stationary blocking probability of the insensitive policy is thus minimal in the considered class of policies, and it is independent of θ. Hence, even a buffer of size 1 is sufficient to get the optimal stationary behavior. We would like to point out that the optimality is only valid in the limit n → ∞. In order to compute the blocking probability (or other performance measures) for values of n that are large but finite, one has to look at finer scales, which will be the objective of Section 5.

Finer scales and estimates

While the results of the previous section are interesting for some performance metrics like the mean number of customers (which gives the mean waiting time via Little's formula), they are too rough to be really informative in terms of blocking probabilities. Indeed, any reasonable dynamic load balancing may achieve the given bounds. To get useful and discriminative estimates, we hence need to investigate the process S (n) at finer scales. In particular, we aim at determining when blocking can be considered a large deviation event (with a probability exponentially small in n) and when it will be in the scale of the central limit theorem.

Large deviations

Let P = {p ∈ R θ + : θ i=0 p i = 1}. For c > 0, denote S (n) c
= {s ∈ S : s = nc}, and

P (n) c = {q ∈ P : nq ∈ S (n) c }. Since s = k ks k , we have k q k = c, ∀q ∈ S (n) c . Thus, P (n) c
is the set of discrete probability distributions taking values on a lattice of unit size 1/n and having a first moment of c.

Define p ∈ P (n) c by p k (c) := 1 (θ -k)! θ -c ρ θ-k 1 ψ(c) . ( 27
)
where

ψ(c) = θ k=0 1 k! θ -c ρ k , (28)
is a normalizing constant which ensures that p is a probability vector. There need not be a vector in

P (n) c
satisfying [START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An assymptotic approach[END_REF], in which case we define p to be the vector

2 in P (n) c
which is closest (say in norm l 1 ) to satisfying [START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An assymptotic approach[END_REF]. To simplify the notation, let π (n) (q; c) = π (n) (nq)1 {q∈Pc} be the stationary probability of observing q ∈ P (n)

c . Let c = arg max c∈[0,θ] e c ψ(c). (29)
We first characterize c. As a consequence of the definition of ψ, simple algebric computations show that c coincides (as it intuitively should) with the asymptotic mean value ĉ found in Theorem 3, i.e., Proposition 4. c = ĉ = θ -ρErl -1 θ (1 -ρ). The proof of the proposition is similar to the one for ĉ in Theorem 3. Let p = p(ĉ). The large deviation cost is shown to be a sum of two terms: the "distance" to the stationary point from distributions with equal means plus the cost of having a different mean from the stationary mean. More precisely, we have the following large deviation estimates for S (n) :

Theorem 4. For ρ < 1, lim n→∞ 1 n log π (n) (q; c) π (n) (p; ĉ) = (c -ĉ) + log ψ(c) ψ(ĉ) -D KL (q(c) p(c)). (30)
Proof. Applying Stirling's approximation in the term containing nθ in [START_REF] Janssen | Corrected asymptotics for a multi-server queue in the Halfin-Whitt regime[END_REF] and noting that

k kq k = c, we get π (n) (q; c) ∼ B (n) θ (2π(nθ -nc)) 1/2 e -nθ+nc • n nq θ k=0 1 (θ -k)! nθ -nc nρ θ-k nq k (31) = B (n) θ (2π(nθ -nc)) 1/2 e -nθ+nc ψ n • n nq θ k=0 p nq k k , (32)
Thus, π (n) (q; c) is proportional to the multinomial distribution with p k as the probability of success of the kth class. Using Stirling's approximation in (32), we get

π (n) (q; c) ∼ B (n) θ (2π(nθ -nc)) 1/2 e -nθ+nc ψ n (2πn) -θ/2 • θ k=0 p k q k nq k 1 q 1/2 k , (33)
from which the desired result can be deduced.

Corollary 2. For ρ < 1, lim n→∞ - 1 n log π (n) (q; c) π (n) (p; c) = D KL (q p). (34)
Corollary 2 says that, conditioned on observing nc jobs in the system, the probability of observing a certain distribution of jobs over the servers concentrates around p. The probability of observing any other q ∈ P (n) c decreases exponentially with rate nD KL (q p), that is the Kullback-Liebler distance serves as the large-deviations rate function. This result is akin to Sanov's theorem in information theory [START_REF] Cover | Elements of Information Theory[END_REF].

Corollary 3. For ρ < 1, lim n→∞ 1 n log π (n) (p; c) π (n) (p; ĉ) = (c -ĉ) + log ψ(c) ψ(ĉ) . ( 35
)
Corollary 3 states that the scaled number of tasks in the system concentrates around ĉ exponentially with rate(cĉ) + log(ψ(c)/ψ(ĉ)).

Asymptotics for the blocking probability

In this section, we shall look at the asymptotics for the blocking probability which is the main performance measure of interest for our system. Two different asymptotic regimes will be considered: (i) the number of servers, n scales, linearly with the arrival rate, nρ; and (ii) the Halfin-Whitt regime [START_REF] Halfin | Heavy-traffic limits for queues with many exponential servers[END_REF] which has a linear term as in (i) along with a sub-linear term that represents the safety margin.

The starting point for both these asymptotic regimes will be an integral characterization of B θ which is derived from the generating function of the stationary measure (see Theorem 13 in the supplementary material).

Theorem 5. For ρ ∈ (0, 1), the blocking B θ has the asymptotic form:

lim n→∞ B (n) θ exp(nR(γ θ,ρ )) 2πn α θ,ρ 1/2 = 1. ( 36
)
where

R(t) = log θ k=0 t k k! -ρt, (37)
γ θ,ρ = arg max t∈(0,∞) R(t) = θ - ĉ ρ , (38)
and

α θ,ρ = (1 -ρ) ρ θ ργ θ,ρ -1 . (39)
Corollary 4. For θ = 1, γ θ,ρ = 1-ρ ρ -1 and α θ,ρ = 1. Thus,

B (n) 1 ∼ e n(1-ρ) ρ n (2πn) -1/2 . ( 40
)
For the proof of Theorem 5, we shall need the following result whose proof follows along the same lines as that of Theorem 3.

Lemma 1. Let γ θ,ρ = arg max t∈(0,∞) R(t). Then, γ θ,ρ is the unique solution of the equation

(1 -ρ) θ k=0 x k k! = x θ θ! , (41)
is x = Erl -1 θ (1 -ρ).
Proof of Theorem 5. Using calculations on the generating functions for the stationary probability (see the supplementary material), we obtain that

lim n→∞ B (n) θ nρ ∞ 0 θ k=0 1 k! t k n e -tnρ dt = 1. ( 42
)
The asymptotic form of the integral can be determined using Laplace's method, which says that

∞ 0 e nf (t) dt ≈ e nf (t0) 2π n(-f (t 0 )) 1/2 , ( 43
)
where t 0 is the unique maximizer of f in (0, ∞).

Define:

R(t) = log θ k=0 t k k! -ρt. (44)
Then

lim n→∞ B (n) θ e nR(t0) 2πnρ 2 -R (t 0 ) 1/2 = 1. ( 45
)
where t 0 maximizes R and is characterized in Lemma 1. For Laplace's method to be applicable, one needs R (t 0 ) < 0. This is shown in Lemma 2 which appears in Appendix A.3. We shall now compute α θ,ρ = -ρ -2 R (γ θ,ρ ), and the main result then follows from (45).

Let

g θ (t) = θ k=0 t k k! . Since γ θ,ρ is the maximizer of R(t), we have R (γ θ,ρ ) = 0, which upon rearrangement gives g θ-1 (γ θ,ρ ) g θ (γ θ,ρ ) = ρ. ( 46
)
From the definition of g θ and Lemma 1, we have

g θ (γ θ,ρ ) = 1 1-ρ γ θ θ,ρ
θ! , and

g θ-2 (γ θ,ρ ) = g θ-1 (γ θ,ρ ) - γ θ-1 θ,ρ (θ -1)! (47) = γ θ θ,ρ θ! ρ 1 -ρ - θ γ θ,ρ (48)
Thus,

R (γ θ,ρ ) = g θ-2 (γ θ,ρ ) g θ (γ θ,ρ ) - g θ-1 (γ θ,ρ ) g θ (γ θ,ρ ) 2 (49) = ρ - (1 -ρ)θ γ θ,ρ -ρ 2 (50)
from which the expression for α θ,ρ follows.

For ρ > 1, we cannot use the directly use the technique that was used for ρ < 1 because the maximum of R(t) in the interval [0, ∞) occurs at t = 0, which is not an interior point of the support of R. So, we shall resort to a theorem due to Erdelyi that treats this case. Theorem 6. Let ρ > 1 and n(ρ -1) be bounded away from 0. As n → ∞,

B (n) θ ∼ 1 -ρ -1 + 1 (ρ -1) θ n θ + o(n -θ ). ( 51
)
Proof. We shall apply the Erdelyi theorem (see Theorem 1.1 in the arXiv preprint of [START_REF] Nemes | An explicit formula for the coefficients in Laplace's method[END_REF] for a precise statement with the notation relevant for our proof) with

f (t) = -R(t) and [a, b) = [0, ∞).
Let us verify that the four conditions of this theorem are satisfied by -R(t).

For ρ ≥ 1, the function -R(t) is increasing in the interval [0, ∞) with minimum at t = 0. To see this,

-R (t) = ρ - θ-1 k=0 t k k! θ k=0 t k k! ≥ 1 - θ-1 k=0 t k k! θ k=0 t k k! ≥ 0.
From Lemma 4, in a neighborhood of 0, -R(t) is analytic with expansion

-R(t) = (ρ -1)t + t θ+1 (θ + 1)! + o(t θ+1 ), (52)
and -R(t) is continuously differentiable with an analytic derivative. Finally, to show the absolute convergence of the integral, note that

∞ 0 e -n(-R(t)) dt ≤ ∞ 0 e -n(ρ-1)t dt = 1 n(ρ -1) < ∞, (53)
as long as n(ρ -1) is bounded away from 0. Thus, all the necessary conditions required by Erdelyi theorem are satisfied. The various parameters in the asymptotic expansion (1.5) in [START_REF] Nemes | An explicit formula for the coefficients in Laplace's method[END_REF] are: R(0) = 0, α = 1, a 0 = ρ -1, a 1 = . . . = a θ-1 = 0, a θ = 1 (θ+1)! , which gives β 0 = (ρ -1) -1 , β 1 = . . . β θ-1 = 0 and

β θ = -(ρ-1) -(θ+2) θ! , so that ∞ 0 e nR(t) dt ∼ 1 (ρ -1)n - 1 (ρ -1) θ+2 n θ+1 + o(n -(θ+1) ) (54)
Substituting the above asymptotic expansion in (42), we get the claimed result.

Since R(0) = 0, there is no exponential decay of the blocking probability when ρ > 1.

Corollary 5. Setting θ = 1 in the above theorem, we get the corresponding result for θ = 1 obtained in [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF] (see Theorem 13 in there).

The previous theorems give the asymptotics of the blocking probability for a fixed load per server for large number of servers. For ρ < 1, the blocking probability goes to 0 exponentially quickly in n while for ρ > 1 it goes to 1 -ρ -1 , a strictly positive quantity. The next theorem looks at the scaling law that results in a polynomial blocking probability. For the Erlang C model, that is, a system without blocking, this regime has the following interpretation: if the cost of servers is high, Halfin and Whitt [START_REF] Halfin | Heavy-traffic limits for queues with many exponential servers[END_REF] observed that it could be beneficial to reduce the number of servers in such a way such that the probability of waiting is no longer exponentially small but decays as n -1/2 . This increase in the waiting probability has the benefit of requiring λ + O(λ 1/2 ) instead of bλ, b > 1, servers. Thus, one gains in the cost and the utilization of servers at the expense of the waiting probability. In order to evoke this trade-off between these two quantities, this scaling regime is also called the Quality and Efficiency Driven (QED) regime. We note that this asymptotic regime was already studied for the Erlang B system in Jagerman [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF] (see Theorem 14) but the interpretation in terms of a trade-off is due to Halfin and Whitt [START_REF] Halfin | Heavy-traffic limits for queues with many exponential servers[END_REF] for systems without blocking and Whitt [START_REF] Whitt | Heavy traffic approximations for service systems with blocking[END_REF] for systems with blocking.

The following theorem gives the QED scaling for the balanced load-balancing policy and can be viewed as a generalization of the QED result for the Erlang loss model.

Theorem 7. For a ∈ (∞, ∞), let nρ = n + an 1/(θ+1) . ( 55
)
Then,

lim n→∞ B (n) θ n θ/(θ+1) ∞ 0 exp au - u (θ+1) (θ + 1)! du = 1. ( 56
)
Proof. The proof follows a similar reasoning as in the proof of Theorem 14 in [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF]. From (42) and using Lemma 4

1 = lim n→∞ B (n) θ nρ ∞ 0 θ k=0 1 k! t k n e -tnρ dt (57) ∼ lim n→∞ B (n) θ nρ ∞ 0 exp nt -n t (θ+1) (θ + 1)! -tnρ dt (58) ∼ lim n→∞ B (n) θ nρ ∞ 0 exp an 1/(θ+1) t -n t (θ+1) (θ + 1)! dt ( 59
)
Setting u = tn 1/(θ+1) in the integral gives:

∼ lim n→∞ B (n) θ (n θ/(θ+1) + a) ∞ 0 exp au - u (θ+1) (θ + 1)! du ( 60
)
Note that a can be positive or negative, which means that even with a total charge larger than the number of servers, the blocking probability can decay to 0 provided that (55) is satisfied asymptotically. When a = 0 and using simple computations, the Theorem leads to the following corrollary:

Corollary 6. If ρ = 1: B (n) θ ∼ (θ + 1)! 1 θ+1 θ + 1 Γ 1 θ + 1 n -θ/(θ+1) , ( 61
)
where Γ is the Gamma function.

Note that for θ = 1, we retrieve that:

B (n) 1 ∼ (0.5πn) -1/2 . ( 62
)

Moderate deviations

Using the previous estimates on the blocking probability (i.e. on the normalizing constant of the stationary distribution) we can now characterize the deviations around p of size smaller than O(n) for a fixed value of ρ. Three amplitudes of deviations will be identified according to whether ρ < 1, ρ = 1, and ρ > 1. The proof for the three results in this subsection appear in the appendix and supplementary material. The first result is for ρ < 1 and is a central-limit-theorem-type scaling when the deviations around the mean are of the order of √ n.

Theorem 8. For ρ < 1, 1 √ n S (n) (∞) 0≤i<θ -n(p) 0≤i<θ d ----→ n→∞ N (0, Σ), ( 63
)
where

Σ -1 = ψ(1, 1, . . . , 1) • (1, 1, . . . , 1) - 1 θ -ĉ (θ, θ -1, . . . , 1) • (θ, θ -1, . . . , 1) +      1/p 0 0 . . . 0 0 1/p 1 . . . 0 . . . . . . . . . . . . 0 0 . . . 1/p θ-1      (64)
Corollary 7. For ρ < 1, θ = 1, we have ĉ = ρ, p0 = 1 -ρ and ψ = ρ -1 leading to Σ -1 = ρ -1 and:

S (n) 0 (∞) -n(1 -ρ) √ nρ d ----→ n→∞ N (0, 1). ( 65
)
The next case corresponds to ρ = 1. For a, z ∈ R and θ ≥ 1, define

Φ θ (z; a) = ∞ z exp au - u (θ+1) (θ + 1)! du. ( 66
)
For θ = 1 and a = 0, (2π) -1/2 Φ θ reduces to the complementary cumulative distribution function of the standard normal distribution.

Theorem 9. For ρ = 1 and z ∈ R + ,

lim n→∞ P S (n) θ-1 (∞) n θ/(θ+1) > z = Φ θ (z; 0) Φ θ (0; 0) , ( 67
)
Corollary 8. For ρ = 1, θ = 1, and z > 0,

lim n→∞ P S (n) (∞) √ n > z = 2(1 -Φ(z)), ( 68
)
where Φ is the distribution function of the standard normal distribution.

Unlike in the ρ < 1 case, the deviations are now no longer of O( √ n) but are of higher order. On the other hand, the fluctuations take S (n) with high probability only to states with either θ -1 or θ jobs. All other configurations are on a scale lower than n θ/(θ+1) . This is in contrast with the behavior for ρ < 1 where the fluctuations can take the process to states with number of jobs ranging from 0 to θ. Thus, for ρ = 1, conditioned on being accepted, a customer has a high probability of being routed to a server θ -1 jobs. This property has a direct consequence on the state information the dispatcher needs to take routing decisions. We shall elaborate upon this in Section 6.

Finally, for ρ > 1, the following result shows that the deviations around nθ are of O(1) and are geometrically distributed. Moreover, the excursions take S (n) only to states with S (n) θ-1 > 0 and S (n) i = 0 for i < θ -1. That is, at a random time, there will be a geometrically distributed number of servers with θ -1 clients and there will be no servers with less than θ -1 clients. We give more precise on the blocking probability later on.

Theorem 10.

For ρ > 1, S (n) θ-1 (∞) d ----→ n→∞ Geo(ρ -1 ), ( 69
)
and the blocking probability is

B (n) θ ∼ 1 -ρ -1 . ( 70
)
The previous theorems give a more precise characterization of the system state as well as its performance for a fixed value of ρ (not depending on n) both in terms of blocking and waiting time. In particular, for ρ < 1, the blocking is exponentially small in n while the mean sojourn time is

1 ρ i θ -i θ - ĉ ip i ,
with a deviation of order 1 √ n .

Numerical experiments

We first provide a comparison in Figure 1, of the blocking probability obtained by the insensitive policy analysed in this paper with that of two other policies, namely JSQ and JIQ. The results were obtained through simulations. There are 20 servers each with a buffer size of 10. Two different job-size distributions were used: (i) exponential; and (ii) a discrete distribution, which we call custom, with point masses at 0.1 and 10. The probability of job-size being 0.1 (resp. 10) was 9/9.9 (resp. 0.9/9.9). The JSQ policy is known to be optimal for exponential job-size distributions and homogeneous server speeds, and thus gives a natural benchmark for comparison. The JIQ policy is an interesting policy from the practical point of view as it requires little state information compared to JSQ and the insensitive policy while at the same time it is optimal in the mean-field limit, that is its blocking probability goes to 0 when the number of servers goes to ∞. We have not included the JSQ(d) policy in our comparison as this policy is not optimal in the mean-field limit. We observed in the simulations that in the symmetric case and for high loads, the performance of JSQ and JIQ changes very little when changing the job-size distribution.

While JIQ requires less state information, it can be see that even for a load of 0.9, a few orders of magnitude of gains can be obtained by using the state information. The drawback of JIQ comes from the fact that at high loads, it behaves more and more like Bernoulli routing since there are fewer empty servers available. Thus, while JIQ is optimal in the mean-field limit, the number of servers required to get close this limit will be much higher than that of JSQ or the insensitive policy, motivating the asymptotics expressions for fixed n that we provided. For systems with a smaller number of servers in which state information can be obtained relatively cheaply, JSQ or balanced policies can give a considerable performance advantage over JIQ.

As mentioned in the Introduction, in the case of symmetric speeds, our motivation for studying the insensitive policy comes mainly from the fact that precise asymptotic estimates can be obtained which is obviously not the case with JSQ and JIQ. For asymmetric speed, insensitive policies might actually present performance gains over JSQ but this falls out of the scope of this paper.

We now illustrate the relationship between the blocking probability and the various parameter such as θ, n and λ for the insensitive policy. First, we evaluate the predictive abilities of some of the results obtained in this section by comparing them with the blocking probability obtained from simulating the Markov chain S (n) . In Figure 2, we plot the blocking probability for n = 200 servers and for different values of ρ and θ. The theoretical values were calculated using Theorem 5 for ρ < 1, Theorem 9 for ρ = 1, and Theorem 10 for ρ > 1.

We observe that even for n = 200 the prediction is already reasonably accurate for θ = 2 and θ = 3, except for loads very close to 1 where the accuracy is less (this comes from a singurality in the expression of the blocking probability at 1).

Theorem 7 says that the decay of the blocking probability goes from polynomial to exponential when the load per server ρ is below 1 + an -1/θ . As θ increases the frontier between the exponential and polynomial decay goes closer to ρ = 1. In other words, for a given n as θ increases, the blocking probability starts to decay exponentially from a value to ρ which is closer to 1. This phenomenon is shown in figure 3, in which the blocking probability was computed using Theorem 5.

For the final comparison, we illustrate the benefits of resource pooling. We shall compare three different systems, which we shall index by θ, each corresponding to θ = 1, 2, 3. In the system θ, there will be n/θ servers each of service rate θ and buffer size of θ. That is, the three systems have the same total service rate but differ in the buffer size. For the θ = 1 system we took 300 servers so that the θ = 2 (resp., θ = 3) system had 150 (resp., 100) each of service rate 2 (resp., 3). Figure 4 illustrates the blocking probabiilty as a function of ρ for the three systems. The data for this plot was obtained through simulations. For example, for a load of 0.9 the blocking probability for θ = 1 system is 4.4 × 10 -3 while for the system with θ = 2 the corresponding value is 6 × 10 -5 which represents a considerable reduction in 6 Engineering insights and future work

Performance planning

The asymptotic analysis of insensitive load balancing allows to give a conservative planning for managing the performance relationship engaged between delay guarantees depending on θ, blocking guarantees depending both on n, θ and ρ, and given levels of loads. Indeed, in many applications, a given level of quality of service in terms of delay has to be reached and this can be done by fixing θ. For a given buffer depth θ the mean delay of a job entering the system will be less than θ (the server speed and the mean job-size are fixed to 1). On the other hand, for a given θ and n, we have precisely characterized the asymptotics of the blocking probability, unveiling the critical load ρ c (n) as the frontier of the acceptable blocking probability for most applications. Hence, one can adapt the number of servers n to cope with a target blocking probability given the load or adapt the load given the number of servers. Note that this planning is completely out of reach for specific sensitive policies. Another way of looking at it is by considering the staffing rule which is the number of servers necessary to obtain a vanishing blocking probability in the limit when the total charge is large. In [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF] and [START_REF] Whitt | Heavy traffic approximations for service systems with blocking[END_REF], the staffing rule for θ = 1 was shown to be λ + O(λ 1/2 ), that is, at least these many servers are required to get a vanishing blocking probability when λ is large.

Theorem 7 generalizes the known results for θ = 1 to larger values of θ leading to the following staffing rule: Proposition 5. For a fixed target blocking probability, the number of servers has to scale as λ + aλ 1/(θ+1) , where a is determined by the target blocking probability and can computed using (56).

Practical schemes under the critical load

One of the major criticisms of state-dependent policies such as JSQ or the policy under study in this paper is that the dispatcher needs to know the state of every server in order to route an incoming job. The process of collecting state information can add significant delays and lead to lost revenue [START_REF] Lu | Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services[END_REF]. Practical policies such as the JSQ(d) [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF] or the JIQ [START_REF] Lu | Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services[END_REF] play on the trade-off between information and optimality, and aim to perform much better than state-independent policies while at the same time needing much less information than the whole set of servers. For example, JSQ(d), with the knowledge of the state of only d (which can be fixed number independent of n) servers, has a considerable gain at least in the case of exponentially distributed job-sizes and in the absence of blocking when d = 2 compared to d = 1. While at first glance, the insensitive load balancing policy seems to require full state information, Theorem 9 lends a helping hand in alleviating this need. Recall that this theorem has the following implication: for ρ = 1 and n large, most of the servers will have either θ or θ -1 jobs. One possible scheme to exploit this property is based on the idea first proposed for JIQ, which was motivated by the observation that collecting state information at arrival instants should be avoided in order to reduce delays for jobs. In JIQ, the servers inform the dispatcher (or leave information on a bulletin board) when they become idle. The dispatcher3 then knows which servers are idle, and it routes an incoming packet to one of these servers, if there is one, otherwise it routes based on no information. Thus, upon arrival a job can be routed immediately based on state information collected previously.

For the insensitive policy one can conceive a scheme in which servers inform the dispatcher whether they have θ -1 or fewer than θ -1 jobs (this scheme automatically implies that the dispatcher also knows which servers have θ jobs). When a job arrives, the dispatcher will need to determine the state of only those servers with less than θ -1 jobs. Since this number is expected to be on a smaller scale than n θ/(θ+1) (thanks to Theorem 9), one can expect to reduce the information flow between the servers and dispatchers at arrival instants. One of our future works will be to characterize precisely the variations in the number of servers with fewer than θ -1 jobs. A back of the envelope calculation based upon the proof of Theorem 9 leads one to believe that there will O(n (k+1)/(θ+1) ) servers with k jobs and hence O(n (θ-1)/(θ+1) ) servers with less that θ -1 jobs but this remains to be rigorously investigated.

Of course this reasoning is valid for a given blocking probability of order n -θ θ+1 and this could be significantly reduced for other blocking targets (and hence other loads).

Multi-speed servers

This planning is of course simplified by the fact that we considered a symmetric system depending only on three possibly inter-dependent parameters (n, ρ, θ). In a future work, we aim at generalizing the analysis to servers with different speeds or even to servers with state dependent speed. Though this generalization falls out of the scope of this paper, let us underline the possibility of this analysis by giving its first step, the expression of the stationary measure for the occupation of a multi-speed server farm.

Consider a server farm with n servers that are classified according to their speed into J different types. A server of type j has speed c j , buffer size of θ j , and there are n j servers of type j.

As for the symmetric system, it is convenient here to study the number of servers processing jobs instead of the number of jobs being processed. Let S j = {s ∈ {0, 1, . . . , n j } θj : θj -1 i=0 s i ≤ n j }, and let S (n) = J j=1 S j . Further, let {S (n) (t) ∈ S} t≥0 , where n = J i=1 n j . Let S (n) (t) be a random process defined on S (n) (t), where the component (i, j) of S (n) (t) denotes the number of servers of type j with i customers at time t.

We shall use a boldface font to denote an element of S j , and use calligraphic font to denote an element of S (n) . So, s j would be an vector in S j , and an element s ∈ S (n) can be written as s = (s 1 , s 2 , . . . , s J ).

In state s, the arrival rate to servers of type j and i tasks is given by

λ i,j (s) = (θ j -i)s i,j j (n j θ j -sj ) nρ, (71)
where sj = i is i,j is the total number of tasks in severs of type j.

Theorem 11. If the job-size distribution is exponential, the process S(t) is a reversible Markov process and its stationary distribution of S (n) (t) is given by

π(s) = π(0) (nθ -s)! (nθ)! J j=1 n j s j θ k=0 θ j ! (θ j -k)! (nρ j ) k s k,j , (72)
where s = j sj is the total number of tasks in the system, and ρ j = ρc -1 j . Given this first result established, all the steps presented in the presented analysis might (and should) be considered. This would in particular allow to characterize the optimal trunk reservation parameters (θ i ) for various trade-off of loads, delays and blocking.

Future research directions

Other than the directions described in the previous subsection, several open questions deserve attention. A natural related model would be the generalization of the Erlang C model, that is, the model in studied in this work but with a common waiting room where arrivals wait when all the servers are in the blocking phase. More fundamental questions that merit investigation are:

• Can similar results be established for sensitive policies (like join the shortest of d queues among n)? Are the meaningful scaling similar?

• Can we quantify the optimality gaps for specific families of jobs-size distributions?

• Can we obtain even finer estimates for the blocking probabilities in the QED regime, in the spirit of the body of work establishing precise asymptotics for Erlang's formula [START_REF] Janssen | Corrected asymptotics for a multi-server queue in the Halfin-Whitt regime[END_REF].

A Proof of Theorem 5

Proof. We first prove a local convergence. Let q = p+β/ √ n, and let c = k kq k , β = i iβ i .

Since k q k = k pk = 1, we have

k β k = 0, c = ĉ + β √ n . (73)
We remind the reader that in order to simplify notation, we shall use p instead of p(c). Starting from (33),

π(q) π(p) ∼ ψ(c) ψ(ĉ) n e √ n β k p k q k nq k (74) = ψ(c) ψ(ĉ) n e √ n β k pk q k p k pk nq k (75) = e √ n β k pk q k p k ψ(c) -1 pk ψ(ĉ) -1 nq k (76) = e √ n β k pk q k nq k k p k ψ(c) -1 pk ψ(ĉ) -1 nq k (77)
We shall compute the asymptotics of the two products separately. The first product gives

k pk q k nq k = k 1 + β k pk √ n -n pk - √ nβ k (78) ∼ k exp - √ nβ k - β 2 k 2p k (79) = k exp - β 2 k 2p k , (80)
where the last equality follows from (73). For the second product, from [START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An assymptotic approach[END_REF],

log p k ψ(c) pk ψ(ĉ) = log θ -ĉ -β/ √ n ρ θ-k θ - ĉ ρ -(θ-k) (81) ∼ log 1 - (θ -k) β (θ -ĉ) √ n + (θ -k)(θ -k -1) 2 β2 (θ -ĉ) 2 n (82) ∼ -(θ -k) β (θ -ĉ) √ n - (θ -k) 2 β2 (θ -ĉ) 2 n (83) Thus, k p k ψ(c) -1 pk ψ(ĉ) -1 nq k ∼ exp (np k + √ nβ k ) -(θ -k) β (θ -ĉ) √ n - (θ -k) β2 2(θ -ĉ) 2 n (84) ∼ - √ n β + β2 2(θ -ĉ) , (85)
where the equalities (73), k β k = β and k pk = ĉ helped in the simplification. Substituting the asymptotics of the two products in (77), we get

π(q) π(p) = exp β2 2(θ -ĉ) k exp - β 2 k 2p k (86)
Our starting point is again (33) which for the present case reduces to:

π(q) π(p) ∼ ψ(c) n e n β(n) k p k q k nq k (90) = ψ(c) q θ nq θ e n β(n) θ-1 k=0 p k ψ(c) β (n) k nq k (91) = ψ(c) q θ nq θ e n β(n) θ-1 k=0 1 (θ -k)! (θ -ĉ -β(n) ) θ-k β (n) k nq k (92) = ψ(c) q θ nq θ e n β(n) θ-1 k=0 1 (θ -k)! (-β(n) ) θ-k β (n) k nq k . ( 93
)
Since β (n) ∼ 0 and β(n) < 0, the value of k < θ that makes the largest contribution is θ -1. For all other values of k,

( β(n) ) θ-k β (n) k
→ 0 with respect to this fraction for k = θ -1. That is, fluctuations under this scaling will be visible only in S

θ-1 and S

(n) θ

and not in lower values of k, This further implies that, given the number of jobs in the system, there is only possible configuration of servers possible. In other words, given the number in the system, we can immediately deduce the configuration:

S (n) θ-1 = nθ -nc and S (n) θ = n -S (n)
θ-1 . Therefore, there in only one vector p in the set P (n) c . As a consequence, the only possible value of q in (90) is p, which then leads to:

π(q) π(p) ∼ ψ(c) n e n β(n) . (94)
Consider q θ-1 = β (n) ≥ 0, where β (n) is a scalar from now on. Since q θ = 1 -β (n) , we have β(n) = -β (n) . Let us compute the asymptotics of the term with ψ:

n log(ψ(c)) = n log   θ j=0 β (n) j j!   ∼ n β (n) - β (n) θ+1 (θ + 1)! ( 95
)
where the last asymptotic form is a consequence of Lemma 4. Substituting the above relation back in (94), we get π(q) π(p) ∼ exp -n

β (n) θ+1 (θ + 1)! , (96)
where we have used the identity β(n) = -β (n) which was noted previously. Consequently, the right scaling for β (n) is zn -1/(θ+1) , for z > 0, which means that S

(n)

θ-1 = nβ (n) lives on a scale of n θ/(θ+1) . As for the proof of the central-limit Theorem, we can pass from local to global convergence combining (96), the estimate on the blocking probabilities given in Theorem 7, and Theorem 1.29 in [START_REF] Sagitov | Weak Convergence of Probability Measures[END_REF] with δ 1 (n) = 1

√ n and k = 1.

A.2 Proof of Theorem 10

Proof. Following the same steps as in the proof of theorem 9 until (93), we can arrive at the conclusion S

(n)

θ-1 and S

(n) θ will be non-zero. Note that the only difference with the ρ = 1 case is that now

ψ(c) = θ j=0 (-β(n) ) j ρ j j! , (97)
and p k has a factor ρ -(θ-k) . Going further until (96) leads us to:

π(q) π(p) ∼ exp -nβ (n) + n β (n) ρ -n β (n) θ+1 ρ θ+1 (θ + 1)! . ( 98
)
The only possible scaling, is thus, β (n) = zn -1 , which means that the fluctuations of S We cannot carry on from this stage onwards in the same line as that in the proof of theorem 10 because to arrive at (94) we had assumed that the non-zero fluctuations we increasing with n (this was needed to apply Stirling's approximation). So, we shall work directly with the stationary distribution. From [START_REF] Janssen | Corrected asymptotics for a multi-server queue in the Halfin-Whitt regime[END_REF],

P(S (n) θ-1 = s) = B (n) θ s! n! s!(n -s)! (nρ) -s (99)
= P(S

(n) θ-1 = 0) n! (n -s)! (nρ) -s (100)
∼ P(S

(n) θ-1 = 0)ρ -s , (101)
which is a consequence of Stirling's approximation.

A.3 Concavity of R

Lemma 2. The function R : R + → R defined by

R(t) = log θ k=0 t k k! -ρt, (102)
is concave.

Proof. Recall that g θ (t) = θ k=0 t k k! . Rewrite g θ (t) in terms on the incomplete gamma function using the following steps:

g θ (t) = 1 Γ(θ + 1, 0) ∞ 0 (t + u) θ e -u du (103)
= e t Γ(θ + 1, t),

where Γ is the normalized incomplete gamma function, that is, Γ(m, x) = Γ(m,x) Γ(m,0) . To show the concavity of R, we shall show that its second derivative is negative. Note that g θ (t) = g θ-1 (t) so that

R (t) = g θ-1 (t) g θ (t) -ρ, (105)
and

R (t) = g θ (t)g θ-2 (t) -g θ-1 (t) 2 g θ (t) 2 (106) = Γ(θ + 1, t) Γ(θ -1, t) -Γ(θ, t) 2 Γ(θ + 1, t) 2 . ( 107
)
It is shown in [START_REF] Alzer | Functional inequalities for the incomplete gamma function[END_REF] that Γ viewed as a function of θ is log-concave for all t > 0. We can thus infer that R is concave in (0, ∞).

A.4 Generating functions

((nρ) -1 t) (θ-k) z k (θ -k)! s k e -t dt (112) = B θ ∞ 0 θ k=0 1 k! t k (nρ) -k z θ-k n e -t dt, (113)
where the last identity a consequence of the multinomial theorem followed by a relabelling of the index inside the sum. Finally, making the transformation t → tnρ inside the integral completes the proof. Proof. The result is a direct consequence of Theorem 12 and Lemma 3.

B Miscellaneous results

Lemma 4. For θ ≥ 1, log θ i=0 t i i! = t -1 (θ + 1)! t θ+1 + o(t θ+1 ), as t → 0. (122)

Proof. Let h θ (t) = log θ i=0 t i i! . The proof is based on computing the coefficients in Taylor series expansion of h around 0, that is, The first derivative of h is:

h (1) θ (t) = 1 - t θ θ! g θ (t) -1 , (123)
where g θ (t) = θ i=0 t i i! , which gives the coefficient of t as 1. For k ≤ θ, taking the kth derivative of h (1) and evaluating it using the product rule for higher order derivatives, we obtain the (k + 1)th derivative of h as:

h (k+1) θ (t) = - k j=0 k j t θ-j (θ -j)! (g θ (t) -1 ) (k-j) , (124)
where (g θ (t) -1 ) (k-j) is the (kj)th derivative of g θ (t) -1 . Assuming that the derivatives of g θ (t) -1 do not go to ∞ at t = 0 (which can be seen to be true), at t = 0, the only non-zero derivative is obtained for k = θ and j = k. That is,

h (k+1) θ (0) = 0 1 ≤ k < θ; -1 k = θ. (125)

Figure 1 :

 1 Figure 1: Comparison of the blocking probability for different load balancing policies. Number of servers is 20. Buffer size is 10.

Figure 2 :

 2 Figure 2: Comparison of the blocking probability computed from Theorems 5 and 10 with that obtained from simulations. Number of servers is 200.

Figure 3 :

 3 Figure 3: Decay of the blocking probability as predicted by theorem 5. Number of servers is 100.

Figure 4 :

 4 Figure 4: Benefits of resource pooling in terms of the blocking probability. Number of servers in the baseline system is 300.

  generating function of S(n) .Theorem 12.M (n) (z) = B θ (nρ) From (12) and using the fact that x! = ∞ 0 t x e -t dt and (nθs) = k (θk)s k , M(n) (z) = B θ s ∞ 0 t k (θ-k)s k e -t dt n s θ k=0 (nρ) -(θ-k) z k (θk)! ) -1 t) (θ-k) z k (θk)!

  generating function of the number of tasks in steady state.Lemma 3. M(n) (z) = M (n) (z 0 , z 1 , . . . , z θ ). )(z k ) s k (119) = M (n) (z 0 , z 1 , . . . , z θ ). (120)Theorem 13.M(n) (z) = B θ (nρ)

Optimal in the sense that it minimizes the blocking probability or any convex criterion.

When the value of c is clear from the context, we will use the notation p instead of p(c).

We are assuming a single dispatcher.

This work was partially supported by the Basque Center for Applied Mathematics BCAM and the Bizkaia Talent and European Commission through COFUND programme, under the project titled "High-dimensional stochastic networks and particles systems", awarded in the 2014 Aid Programme with request reference number AYD-000-273, and by the STIC-AmSud project N o 14STIC03.

Consider the exponent on the RHS. Since i β i = 0, we have β =

Since the multivariate Gaussian distribution has exponent -1 2 βΣ -1 β, we can deduce from the above equation the inverse of the covariance matrix to be one stated in the theorem and the local convergence of π(q) π(p) until the Gaussian density. Using the approximation in (32), combined with the blocking probability estimates obtained in Theorem 5, it can be easily seen that

which in turns implies that for any q = p + β/ √ n:

Generalizing slightly the previous computations, the same would hold for any q = p + (β + ε n )/ √ n, with ε n vanishing when n goes to infinity. Hence, to derive a global convergence result of the distribution function as stated in the Theorem, we can now appeal to a variant of Scheffé's lemma (see for instance Theorem 1.29 in [START_REF] Sagitov | Weak Convergence of Probability Measures[END_REF] with δ i (n) = 1 √ n , i = 1 . . . , k and k = θ).

A.1 Proof of Theorem 9

Proof. Instead of defining q according to a pre-defined scaling like in the previous proof, we shall this time define it with an arbitrary scaling which shall be made precise later. Let q = p + β (n) , where again we have k β (n) k = 0. For ρ = 1, p0 = . . . = pθ-1 = 0, pθ = 1, so we shall assume that β (n) k ≥ 0 for k < θ. Also, for ρ = 1, we have ĉ = θ and ψ(ĉ) = 1 so that

where β(n) = k kβ (n) k < 0.