Mark Van Der Boor

Céline Comte

Load Balancing in Heterogeneous Server Clusters: Insights From a Product-Form Queueing Model *

Keywords: Load balancing, performance analysis, product-form queueing model, Jackson network, insensitivity

Efficiently exploiting servers in data centers requires performance analysis methods that account not only for the stochastic nature of demand but also for server heterogeneity. Although several recent works proved optimality results for heterogeneity-aware variants of classical load-balancing algorithms in the many-server regime, we still lack a fundamental understanding of the impact of heterogeneity on performance in finite-size systems. In this paper, we consider a load-balancing algorithm that leads to a product-form queueing model and can therefore be analyzed exactly even when the number of servers is finite. We develop new analytical methods that exploit its product-form stationary distribution to understand the joint impact of the speeds and buffer lengths of servers on performance. These analytical results are supported and complemented by numerical evaluations that cover a large variety of scenarios.

Introduction

Distributing a stochastic demand across a set of heterogeneous servers is not only a fundamental problem in queueing theory but also an essential building block of applications like parallel computing and production systems. Besides static approaches that do not require communication between the dispatcher and the servers, classical solutions include join-the-shortest-queue, power-of-d-choices [START_REF] David | The power of two choices in randomized load balancing[END_REF], and join-idle-queue [START_REF] Lu | Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services[END_REF]. These solutions were originally designed for service systems that are homogeneous in the sense that all servers have the same speed [START_REF] Van Der Boor | Scalable load balancing in networked systems: A survey of recent advances[END_REF]. Although these solutions successfully cope with the stochastic nature of demand, they are not always suitable when servers have unequal speeds [START_REF] Gardner | Scalable load balancing in the presence of heterogeneous servers[END_REF].

Several heterogeneity-aware approaches were introduced to improve the performance of these classical solutions, for instance by using information on the job sizes and server speeds or by delaying the assignment decision. For instance, join-the-shortest-workload yields optimal performance if the service times of all jobs over all servers are known, an assumption that is rarely satisfied in practice. Redundancy scheduling achieves the same performance gain [START_REF] Ayesta | On a unifying product form framework for redundancy models[END_REF], this time by delaying the assignment decision, which may again be practically infeasible if the communication time between the dispatcher and servers is not negligible. To achieve good performance without these strong assumptions, more recent works introduced speed-aware variants of the above-mentioned well-known algorithms. More specifically, [START_REF] Weng | Optimal load balancing in bipartite graphs[END_REF] introduced variants of join-the-shortest-queue and join-idle-queue where the server speeds are used as a tie-breaking rule, and proved that these variants minimize the mean response time in the many-server regime; [START_REF] Gardner | Scalable load balancing in the presence of heterogeneous servers[END_REF] proposed variants of power-of-d-choices and join-idle-queue for service systems with two server types (fast and slow) by adapting the degree of diversity and assignment probabilities to the server speeds, and proved stability, again in the many-server regime. Despite these advances, we still lack a fundamental understanding of the impact of heterogeneity on performance in service systems with a finite number of servers.

In this paper, we make one step further into this direction by considering a load-balancing algorithm [START_REF] Bonald | Insensitive load balancing[END_REF][START_REF] Comte | Dynamic load balancing with tokens[END_REF][START_REF] Comte | Resource management in computer clusters: algorithm design and performance analysis[END_REF][START_REF] Jonckheere | Asymptotics of insensitive load balancing and blocking phases[END_REF]] that leads to a product-form queueing model (assuming that jobs arrive according to a Poisson process), and can therefore be analyzed exactly even when the number of servers is finite. Assuming that each server has a finitelength buffer, this algorithm assigns each incoming job to a server chosen at random, with a probability proportional to the number of available slots in the server's buffer; an incoming job that finds the buffers of all servers full is rejected and considered permanently lost. Although this algorithm does not account for the server speeds in the online assignment decision, these speeds can be used offline to adjust the buffer lengths. Besides its analytical tractability, this algorithm has the advantage of making performance insensitive to the job size distribution beyond its mean, provided that servers apply the processor-sharing policy. This insensitivity property, which contributed to the success of the Erlang-B formula for dimensioning circuit-switched networks, guarantees that the long-run performance metrics are not impacted by fine-grained traffic characteristics. In this paper, we use the product-form queueing model to better understand the impact of parameters on performance.

Related work on insensitive load balancing

The works [START_REF] Adan | A loss system with skill-based servers under assign to longest idle server policy[END_REF][START_REF] Bonald | Insensitive load balancing[END_REF][START_REF] Jonckheere | Towards an erlang formula for multiclass networks[END_REF] analyze the performance of variants and generalizations of this load-balancing algorithm in heterogeneous service systems where jobs have constraints that restrict their assignment to resources. However, the objective of these works is to develop methods to calculate performance metrics, which is a different goal than understanding the impact of parameters on performance. Some of the formulas derived in these works are used in the numerical-evaluation section to assess our analytical results. The objectives of the related work [START_REF] Jonckheere | Asymptotics of insensitive load balancing and blocking phases[END_REF] are closer to our work. This work analyzes the performance of the same algorithm, but it focuses on the many-server and heavy-traffic regimes and assumes, for the most part, that servers are homogeneous. The product-form stationary distribution of load-balancing models has also been studied in [START_REF] Van Der Boor | Load balancing in large-scale systems with multiple dispatchers[END_REF][START_REF] Comte | Dynamic load balancing with tokens[END_REF][START_REF] Comte | Resource management in computer clusters: algorithm design and performance analysis[END_REF][START_REF] Gardner | Product forms for FCFS queueing models with arbitrary server-job compatibilities: an overview[END_REF] for systems with multiple dispatcher or arbitrary server-job compatibilities (also see references therein).

Additional related work on load balancing If all servers have unit-length buffers, our algorithm can be seen as a loss variant of join-idle-queue [START_REF] Lu | Join-idle-queue: A novel load balancing algorithm for dynamically scalable web services[END_REF], whereby a job is rejected if all servers are busy upon its arrival. With arbitrary buffer lengths, our algorithm is related to idle-one-queue [START_REF] Gupta | Load balancing in the nondegenerate slowdown regime[END_REF] and join-below-threshold [START_REF] Goldsztajn | Self-learning threshold-based load balancing[END_REF][START_REF] Illés | Mean field analysis of join-below-threshold load balancing for resource sharing servers[END_REF][START_REF] Zhou | Flexible load balancing with multi-dimensional state-space collapse: Throughput and heavy-traffic delay optimality[END_REF][START_REF] Zhou | Heavy-traffic delay optimality in pull-based load balancing systems: Necessary and sufficient conditions[END_REF][START_REF] Zhou | Designing low-complexity heavy-traffic delay-optimal load balancing schemes: Theory to algorithms[END_REF], two generalizations of join-idle-queue introduced to improve performance in the heavy-traffic regime or when servers have unequal speeds. The idea is that servers notify the dispatcher when the number of jobs in their buffer falls below a threshold, so that the dispatcher assigns incoming jobs to lightly-loaded servers if possible. The threshold of a server, equal to one under join-idle-queue, offers a trade-off between performance improvement and communication overhead; in case of unequal server speeds, it can also be used to favor faster servers. In our algorithm, these thresholds correspond to the buffer lengths, and the trade-off between performance and communication overhead is materialized by the overall buffer length (that is, the sum of the lengths of the buffers at all servers). Assuming that this overall buffer length is fixed, we would like to understand how to optimally choose the buffer length of each server depending on the job arrival rate and the server speeds.

Contributions Our contributions can be summarized as follows. We first show the following results for a cluster of two servers in which the overall buffer size across the two servers is fixed. When the arrival rate is low, the optimal buffer lengths (to minimize the loss probability) are proportional to the server speeds, meaning that the buffer of a server is longer if this server is faster. On the contrary, when the arrival rate is large, the optimal buffer lengths are uniform, and the server speeds only intervene to break ties if the overall buffer length is odd. We also show that, between these two limiting regimes, the optimal buffer lengths evolve monotonically with the arrival rate. Besides the practical implications of these results in terms of system design, the analytical methods that we develop, based on an analogy with weighted paths in the two-dimensional lattice, are of independent interest. Afterwards, we explain how these results extend to clusters of more than two servers. We finally turn to numerical evaluations to assess the validity of these results in practice and understand the impact of parameters on the mean response time.

Organization of the paper

The remainder of the paper is organized as follows. Section 2 introduces the cluster model and the equivalent closed queueing model. In Section 3, we use this model to derive closed-form expressions for several performance metrics, such as the loss probability. Sections 4 to 6 contain our main contributions for clusters of two servers. In Section 4, we prove that the optimal buffer lengths in terms of the loss probability are proportional to the server speeds when the arrival rate tends to zero, while, in Section 5, we prove that the optimal buffer lengths are uniform when the arrival rate tends to infinity. Section 6 fills the gap between these two limiting regimes by showing a monotonicity result. These results are generalized to clusters of more than two servers in Section 7. Section 8 gives numerical evaluations and Section 9 concludes the paper.

Heterogeneous server cluster

We consider a cluster that consists of a single dispatcher and two servers. Incoming jobs arrive at the dispatcher according to a Poisson process with rate λ. Each server has a finite buffer, of length ℓ1 for server 1 and ℓ2 for server 2, that contains all jobs assigned to this server (either waiting or in service), with the total buffer size being denoted by L = ℓ1 + ℓ2. An incoming job is permanently lost if the buffers of both servers are full upon its arrival, otherwise the dispatcher immediately assigns the job to one of the servers, as will be elaborated further. Two equivalent state descriptors are the vector n = (n1, n2) that counts the jobs in the buffer of each server, and the vector x = (x1, x2) that counts the available slots in the buffers of each server. The generalization of this model to clusters of more than two servers will be considered in Section 7.

µ 1 Server 1 n 1 = 4 µ 2 Server 2 n 2 = 2 Dispatcher λ ?
Figure 1: A heterogeneous cluster with two servers.

Load balancing and scheduling

The dispatcher applies the following randomized load-balancing algorithm, considered in [START_REF] Bonald | Insensitive load balancing[END_REF][START_REF] Comte | Dynamic load balancing with tokens[END_REF][START_REF] Comte | Resource management in computer clusters: algorithm design and performance analysis[END_REF][START_REF] Jonckheere | Asymptotics of insensitive load balancing and blocking phases[END_REF]. When a new job arrives, the dispatcher chooses a server at random, with a probability proportional to the number of available slots in the buffer of the server, and assigns the job to this server. In Figure 1 for instance, there are one and two available slots in the buffers of servers 1 and 2, respectively, so that an incoming job would be assigned to server 1 with probability 1 3 and to server 2 with probability 2 3 . In general, if there are x1 and x2 available slots in the buffers of servers 1 and 2, respectively, then server 1 is chosen with probability x 1

x 1 +x 2 and server 2 with probability x 2

x 1 +x 2 if x1 + x2 ≥ 1, otherwise the job is lost. This algorithm assumes that the dispatcher always knows the number of available slots in the buffer of each server; this happens, for instance, if jobs go through the dispatcher again when they leave the system, or if servers notify the dispatcher upon service completions. This load-balancing algorithm was considered in [START_REF] Bonald | Insensitive load balancing[END_REF][START_REF] Comte | Dynamic load balancing with tokens[END_REF][START_REF] Comte | Resource management in computer clusters: algorithm design and performance analysis[END_REF][START_REF] Jonckheere | Asymptotics of insensitive load balancing and blocking phases[END_REF] for clusters with an arbitrary number of servers. We assume that each server processes the jobs in its buffer according to a non-anticipating work-conserving scheduling algorithm, such as processor-sharing or first-come-first-served. The service rates of the servers, assumed to be constant for simplicity, are denoted by µ1 and µ2. We assume without loss of generality that 1 > µ1 > µ2 > 0 and µ1 + µ2 = 1, so that server 1 is the fastest. The job service requirements are independent and exponentially distributed with unit mean, so that the remaining service time of a job is exponentially distributed with mean 1 µ if this job is currently served at rate µ. This memoryless assumption is actually not required to perform the subsequent analysis if each server applies processor-sharing or preemptive-resume last-come-first-served. Each job leaves the system immediately upon service completion.

Even if the dispatcher does not take the service rates µ1 and µ2 into account when making the assignment decision, we will see later that varying the buffer lengths ℓ1 and ℓ2 allows us to optimize the load distribution with respect to the service rates.

Queueing model

The dynamics can be described by a closed Jackson network [START_REF] Serfozo | Introduction to Stochastic Networks[END_REF] of three stations with two customer (or token) classes [START_REF] Van Der Boor | Load balancing in large-scale systems with multiple dispatchers[END_REF][START_REF] Comte | Dynamic load balancing with tokens[END_REF]. In Section 3, we will use this observation to derive closed-form expressions for several performance metrics.

Instead of counting how many jobs are present in the cluster, we keep track to how tokens evolve in a network, where every token corresponds to a specific slot in a server's buffer. A token corresponding to a slot in the buffer of server i is said to be of class i, for i ∈ {1, 2}. When a slot in the buffer of server i is not occupied, the corresponding class-i token is located at the dispatcher station; when this slot becomes occupied by a job, the class-i token moves to the server-i station. The routing mechanism of tokens is deterministic, and each class-i token moves only between the dispatcher station and the server-i station. More specifically, whenever a job completes service at server i (this happens with rate µi), the corresponding slot becomes available and the corresponding token moves to the dispatcher station. Whenever a job is assigned to server i, a slot in the buffer of this server becomes occupied by this job, so that the corresponding class-i token moves to server-i station.

The corresponding closed Jackson network consists of the dispatcher station and the server-1 and 2 stations. There are ℓi tokens of class i, for i ∈ {1, 2}, and these tokens are either in the dispatcher station or in the server-i station. The service policy in the dispatcher station is processor-sharing, while the service policies in the server stations match those applied by the servers. The state of this network is described by the vector x = (x1, x2) that counts the number of tokens in the dispatcher station, corresponding to available slots in the servers' buffers. In particular, x = 0 means that both buffers are full, while x = ℓ, with ℓ = (ℓ1, ℓ2), means that the cluster is empty. We let e1 = (1, 0) and e2 = (0, 1) denote the states corresponding to a single (slot occupied by a) job at server i. Known results on closed Jackson networks [20, Chapter 1 and Section 3.1] show that the stationary distribution of the Markov process defined by the evolution of the network state is given by

π(x) = 1 G(ℓ) x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 , x ≤ ℓ, (1)
where G(ℓ) is a normalization constant.

Performance analysis

As in Jackson networks with a single class of customers [START_REF] Buzen | Computational algorithms for closed queueing networks with exponential servers[END_REF], the closed-form expressions of several metrics of interest stem directly from the normalization constant.

Normalization constant

The normalization constant follows from the normalization equation x≤ℓ π(x) = 1. Using (1), we obtain

G(ℓ) = x≤ℓ x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 , ℓ ∈ N 2 . ( 2
)
The following geometric interpretation guides the proofs of Theorems 1 and 2 in Sections 4 and 5. Consider the square lattice consisting of the 2-dimensional vectors x = (x1, x2) with integer components. We are interested in the direct paths going from the origin 0 to some vector x ≤ ℓ, where by direct we mean a path that consists only of increasing unit steps in the horizontal or vertical direction. If each horizontal step has weight µ 1 λ and each vertical step weight

µ 2 λ , then ( µ 1 λ ) x 1 ( µ 2 λ )
x 2 is the multiplicative weight of any direct path going from the origin 0 to the vector x. Since there are x 1 +x 2 x 1 such paths, it follows that the normalization constant G(ℓ) is the sum of the weights of all direct paths going from the origin to a vector x ≤ ℓ. Intuitively, this suggests that:

• If λ ≪ µ1 + µ2(= 1)
, the paths that have the most weight are the longest, meaning that the states with large values of x1 + x2 (that is, with many available slots) are the most likely. This intuition guides the proof of Theorem 1.

• If λ ≫ 1, the paths that have the most weight are the shortest, meaning that the states with small values of x1 + x2 (that is, with few available slots) are the most likely. This intuition guides the proof of Theorem 2.

Although [START_REF] Ayesta | On a unifying product form framework for redundancy models[END_REF] unveils interesting properties of the constant G(ℓ), this expression is inconvenient when it comes to compute this constant because it leads to numerical instability. For the numerical results, we use instead the recursive expression

G(ℓ) = 1 + 2 i=1 ℓ i ≥1 µi λ G(ℓ -ei), ℓ ∈ N 2 \ {0},
with the base case G(0) = 1. This expression, which is a special case of the formula derived in [3, Proposition 2], can be seen as a generalization of the Erlang-B formula [START_REF] Krarup | Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges[END_REF]. It will often be convenient to consider the quantity

δG(ℓ) = G(ℓ + e1 -e2) -G(ℓ), ℓ ∈ N 2 : ℓ2 ≥ 1, (3)
that gives the variation of the normalization constant obtained by replacing a server-2 slot with a server-1 slot. By injecting (2) into this definition and making simplifications, we obtain

δG(ℓ) = ℓ 1 +ℓ 2 n=ℓ 1 +1 n ℓ1 + 1 µ1 λ ℓ 1 +1 µ2 λ n-ℓ 1 -1 - ℓ 1 +ℓ 2 n=ℓ 2 n ℓ2 µ1 λ n-ℓ 2 µ2 λ ℓ 2 . ( 4
)

Long-term performance metrics

We now consider three performance metrics called the loss probability, occupation rate, and mean response time.

The formulas below are simple extensions of formulas derived for closed Jackson networks with a single class of customers [START_REF] Buzen | Computational algorithms for closed queueing networks with exponential servers[END_REF]. However, to the best of our knowledge, the results regarding the occupation rate and mean response time have never been derived in the literature on insensitive load balancing.

Loss probability

The loss probability β(ℓ) is defined as the probability that an incoming job is rejected, which happens when the buffers of all servers are full upon its arrival. According to the PASTA property [START_REF] Wolff | Poisson arrivals see time averages[END_REF], the loss probability is equal to the stationary probability π(0) that the buffers of all servers are full, so that then by [START_REF] Adan | A loss system with skill-based servers under assign to longest idle server policy[END_REF] we obtain

β(ℓ) = 1 G(ℓ) . (5)
By (2), the loss probability decreases when the number of slots in a given buffer increases (as intuition suggests).

Occupation rate For each i ∈ {1, 2}, the occupation rate of server i is defined as the fraction of time that this server is busy. According to [START_REF] Adan | A loss system with skill-based servers under assign to longest idle server policy[END_REF], this quantity ρi(ℓ) is given by

1 G(ℓ) x≤ℓ-e i x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 = G(ℓ -ei) G(ℓ) .
We have ρ1(ℓ) > ρ2(ℓ) if and only if G(ℓ -e1) > G(ℓ -e2), which by ( 5) also means that the loss probability increases less when we remove a slot from server 1 than when we remove a slot from server 2.

Mean response time

The response time of a job is defined as the duration between its arrival in the cluster and its departure. We can show that the mean numbers of jobs in the buffers of servers 1 and 2 are given by

α1(ℓ) = ℓ 1 -1 x 1 =0 G(x1, ℓ2) G(ℓ) , α2(ℓ) = ℓ 2 -1 x 2 =0 G(ℓ1, x2) G(ℓ) .
According to Little's law, the mean response time of a job is given by ∆(ℓ) = α 1 (ℓ)+α 2 (ℓ) λ(1-β(ℓ)) .

Problem statement

We would like to understand the joint impact of the arrival rate λ, service rates µ1 and µ2, and buffer lengths ℓ1 and ℓ2 on these performance metrics. Despite the apparent simplicity of these expressions, using them to gain intuition on the impact of parameters on performance is a well-known difficult problem [START_REF] Harel | Convexity results for the erlang delay and loss formulae when the server utilization is held constant[END_REF]. The results of Section 3.2 suggest however that the loss probability is the easiest of these metrics to analyze, as it is simply the inverse of the normalization constant.

We observed earlier that, although the load-balancing algorithm does not account for the service rates µ1 and µ2 to make the assignment decision, the buffer lengths ℓ1 and ℓ2 can be adjusted to optimize performance. Therefore, Sections 4 to 6, which contain our main contributions, will focus more specifically on the following question:

Given the arrival rate λ, service rates µ1 and µ2, and overall buffer length L = ℓ1 + ℓ2, which value(s) of ℓ1 and ℓ2 minimize(s) the loss probability? Since µ1 > µ2, a natural strategy consists of allocating (almost) all slots to server 1. However, the more jobs are in service on server 1, the smaller the fraction of the service rate of this server they receive, and the longer they stay in the system, thus preventing other jobs from entering. This simple observation suggests that it would be better to maximize the minimum service rate received by any job by allocating slots to servers proportionally to their service rates, that is, allocate a fraction µ1 of the slots to server 1 and a fraction µ2 of the slots to server 2. Theorem 1 shows that this allocation is indeed optimal when the arrival rate is small. Theorems 2 and 3 show that, as the arrival rate increases, the optimal allocation evolves monotonically towards a uniform allocation, in which both servers have approximately the same number of slots.

Extending these results to understand the impact of parameters on the mean response time is not straightforward. This impact will be assessed numerically in Section 8.

Low-traffic regime

Our first main contribution is a theorem showing that, when the arrival rate is small, the loss probability is minimized by allocating slots to servers in proportion to their service rate.

Theorem 1. There is a λ * > 0 such that, for each λ ∈ (0, λ * ], the loss probability is minimized when ℓ = (ℓ1, L -ℓ1) with ℓ1 = ⌈µ1L -µ2⌉ or ⌊µ1L + µ1⌋ .

Proof. We first show that, as λ → 0, the monotonicity of G(ℓ) as a function of ℓ is entirely dictated by that of the term corresponding to x = ℓ in (2). We will then study this term. Let ℓ = (ℓ1, ℓ2) ∈ N 2 with ℓ1 + ℓ2 = L. First observe that

x≤ℓ:x =ℓ x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 ≤ 1 λ L -1 1 λ -1
.

By injecting this inequality into (2), we obtain

G(ℓ) = L ℓ1 µ1 λ ℓ 1 µ2 λ ℓ 2 + O λ→0 1 λ L-1 . (7)
If ℓ2 ≥ 1, we can apply this result to both ℓ and ℓ + e1 -e2, so that, by (3), we obtain

δG(ℓ) = ℓ2 ℓ1 + 1 µ1 µ2 -1 • L ℓ1 µ1 λ ℓ 1 µ2 λ ℓ 2 + O λ→0 1 λ L-1
.

As λ tends to zero, the first term tends to +∞ like ( 1 λ ) L , while the second term tends to +∞ at most like ( 1 λ ) L-1 . Therefore, there is λ * > 0 such that, for each λ ∈ (0, λ * ] and each ℓ ∈ N 2 such that ℓ1 + ℓ2 = L and ℓ2 ≥ 1, δG(ℓ) is of the same sign as ( ℓ 2 ℓ 1 +1 µ 1 µ 2 -1), provided that this quantity is nonzero. Now let λ ∈ (0, λ * ] and assume that servers 1 and 2 have ℓ1 and ℓ2 slots, respectively, with ℓ1 + ℓ2 = L. According to the above equality, replacing a slot of server 2 (if any) with a slot of server 1 reduces the loss probability whenever

ℓ 2 ℓ 1 +1 µ 1 µ 2 < 1, that is, ℓ1 < µ1L -µ2.
By taking the symmetrical statement and rearranging the terms, we obtain that replacing a slot of server 1 (if any) with a slot of server 2 reduces the loss probability whenever ℓ1 > µ1L + µ1. Therefore, the slot allocation can only be optimal when µ1L -µ2 ≤ ℓ1 ≤ µ1L + µ1, which is equivalent to [START_REF] Buzen | Computational algorithms for closed queueing networks with exponential servers[END_REF].

Remark 1. The ceil and floor values in [START_REF] Buzen | Computational algorithms for closed queueing networks with exponential servers[END_REF] are different only in the pathological case where µ1L + µ1 = µ1L -µ2 + 1 is an integer, which means that µ1 and µ2 can be written as fractions with denominator L + 1. In this case, the term corresponding to x = ℓ is zero in the expression of δG(ℓ), and which of ℓ1 = ⌈µ1L -µ2⌉ or ℓ1 = ⌊µ1L + µ1⌋ is optimal depends on the terms corresponding to x ≤ ℓ and x = ℓ.

Heavy-traffic regime

The following theorem shows that, when the arrival rate λ is large, the loss probability is minimized by taking ℓ1 ≃ ℓ2. The arrival rates µ1 and µ2 only serve to allocate the remaining slot when the overall buffer length L is odd.

Theorem 2. There is λ * > 0 such that, for each λ ∈ [λ * , ∞), the loss probability is minimized by choosing:

• ℓ1 = ℓ2 = L 2 if L is even; • ℓ1 = L+1 2 and ℓ2 = L-1 2 if L is odd.
Proof. Let ℓ ∈ N 2 such that ℓ1 + ℓ2 = L and ℓ1 ≤ L -1. Using (4), we proceed by exhaustion, distinguishing several cases depending on the values of ℓ1 and ℓ2.

Case 1 (ℓ 1 ≥ ℓ 2 ) We can split the second sum of (4) into two parts, corresponding to n ∈ {ℓ1 + 1, . . . , ℓ1 + ℓ2} and n ∈ {ℓ2, . . . , ℓ1}, respectively. We obtain

δG(ℓ) = ℓ 1 +ℓ 2 n=ℓ 1 +1 n ℓ1 + 1 µ1 λ ℓ 1 +1 µ2 λ n-ℓ 1 -1 - n ℓ2 µ1 λ n-ℓ 2 µ2 λ ℓ 2 - ℓ 1 n=ℓ 2 n ℓ2 µ1 λ n-ℓ 2 µ2 λ ℓ 2 .
As λ tends to +∞, each term in the first sum tends to zero at least as fast as ( 1 λ ) ℓ 1 +1 , while each term in the second sum tends to zero at most as fast as ( 1λ ) ℓ 1 . Therefore, when λ is sufficiently large, we have δG(ℓ) < 0 whenever ℓ1 ≥ ℓ2.

Case 2 (ℓ 1 ≤ ℓ 2 -2) We can split the first sum of (4) into two parts, corresponding to n ∈ {ℓ1 + 1, . . . , ℓ2 -1} and to n ∈ {ℓ2, . . . , ℓ1 + ℓ2}, respectively. We obtain

δG(ℓ) = ℓ 2 -1 n=ℓ 1 +1 n ℓ1 + 1 µ1 λ ℓ 1 +1 µ2 λ n-ℓ 1 -1 + ℓ 1 +ℓ 2 n=ℓ 2 n ℓ1 + 1 µ1 λ ℓ 1 +1 µ2 λ n-ℓ 1 -1 - n ℓ2 µ1 λ n-ℓ 2 µ2 λ ℓ 2 .
As λ tends to +∞, each term in the first sum tends to zero at most as fast as ( 1 λ ) ℓ 2 -1 , while each term in the second sum tends to zero at least as fast as ( 1 λ ) ℓ 2 . Therefore, when λ is sufficiently large, we have δG(ℓ) > 0 whenever ℓ1 ≤ ℓ2 -2.

Case 3 (ℓ 1 = ℓ 2 -1, assuming that L is odd) The two sums in (4) contain the same number of terms, and we obtain:

δG(ℓ) = ℓ 1 +ℓ 2 n=ℓ 1 +1 n ℓ1 + 1 µ1 λ n-ℓ 1 -1 µ2 λ n-ℓ 1 -1 × µ1 λ ℓ 1 +ℓ 2 +1-n - µ2 λ ℓ 1 +ℓ 2 +1-n .
Since µ1 > µ2, it follows that δG(ℓ) > 0.

Conclusion We now gather the three cases. If L is even, the sequence ℓ1 → β(ℓ1, L -ℓ1) is decreasing on {0, 1, . . . , L 2 } and increasing on { L 2 , . . . , L -1, L}. Therefore, the loss probability is minimal when ℓ1 = ℓ2 = L 2 . If L is odd, the sequence ℓ1 → β(ℓ1, L -ℓ1) is decreasing on {0, 1, . . . , L+1 2 } and increasing on { L+1 2 , . . . , L -1, L}. Therefore, the loss probability is minimal when ℓ1 = L+1 2 and ℓ2 = L-1 2 .

Monotonicity

To make the connection between the results of the last two sections, we now show that the optimal slot allocation is monotonic with respect to the arrival rate λ. In the following theorem, with "optimal number of slots allocated to the fastest server", we mean the smallest number of slots allocated to the fastest server that minimizes the loss probability. This will be discussed again in Remark 2.

Theorem 3. The optimal number of slots allocated to the fastest server, in terms of the loss probability, is decreasing with the arrival rate λ.

Proof. It will be convenient to write the loss probability as a function β(λ, ℓ) of both the arrival rate λ and buffer lengths ℓ = (ℓ1, ℓ2). The proof relies on two monotonicity results that are presented in the following two propositions.

Proposition 4. Let λ > 0 and ℓ ∈ N 2 such that ℓ1 ≥ 1 and ℓ2 ≥ 1. If β(λ, ℓ + e1 -e2) ≤ β(λ, ℓ), then β(λ, ℓ -xe1 + xe2) < β(λ, ℓ -(x + 1)e1 + (x + 1)e2) for each x ∈ {0, 1, 2, . . . , ℓ1 -1}. Proposition 5. Let λ * > 0 and ℓ ∈ N 2 such that ℓ2 ≥ 1. If β(λ * , ℓ+e1 -e2) < β(λ * , ℓ), then β(λ, ℓ+e1 -e2) < β(λ, ℓ)
for each λ ∈ (0, λ * ).

Propositions 4 and 5 are proven in Appendices A and B, respectively. These propositions can be rephrased as follows. Assume that, for a given arrival rate λ * and overall buffer length L, we know that allocating ℓ1 + 1 slots to the fastest server yields better performance than allocating ℓ1 slots to this server. Then Proposition 4 shows that, with the same arrival rate λ * , allocating even fewer slots than ℓ1 to the fastest server is even worse in terms of performance. Proposition 5 shows that allocating ℓ1 + 1 slots to the fastest server remains better than ℓ1 slots for any arrival rate λ ∈ (0, λ * ).

Now consider two arrival rates λ * ∈ (0, ∞) and λ ∈ (0, λ * ). Let ℓ λ * and ℓ λ denote the optimal slot allocations under these arrival rates, with

L = ℓ λ * 1 + ℓ λ * 2 = ℓ λ 1 + ℓ λ 2 . Our objective is to prove that ℓ λ * 1 ≤ ℓ λ 1 . The optimality of ℓ λ * implies that β(λ * , ℓ λ * ) < β(λ * , ℓ λ * -e1 + e2). Therefore, Proposition 4 gives β(λ * , ℓ λ * -xe1 + xe2) < β(λ * , ℓ λ * -(x + 1)e1 + (x + 1)e2) for each x ∈ {0, 1, 2, . . . , ℓ λ * 1 -1}. By Proposition 5, each of these inequalities yields β(λ, ℓ λ * -xe1 + xe2) < β(λ, ℓ λ * -(x + 1)e1 + (x + 1)e2) for each x ∈ {0, 1, 2, . . . , ℓ λ * 1 -1}. This in turn implies that ℓ λ * 1 ≤ ℓ λ 1 .
Remark 2. Proposition 4 shows that, among all vectors ℓ ∈ N 2 such that L = ℓ1 + ℓ2, at most two can minimize the loss probability, and these are separated by only one slot. Therefore, in Theorem 3, the "optimal number of slots allocated to the fastest server" is defined up to plus or minus one, and we systematically choose the smallest value for convenience.

Extension to more than two servers

We now consider a cluster that consists of a dispatcher and a set I = {1, 2, . . . , N } of servers. For each i ∈ I, we let µi denote the service rate of server i, ℓi the length of its buffer, and xi the number of available slots in this buffer. We assume without loss of generality that 1 > µ1 ≥ µ2 ≥ . . . ≥ µN > 0 and i∈I µi = 1, and we let L = i∈I ℓi denote the overall buffer length. All definitions of Sections 2 and 3 are generalized in a natural way to this N -server cluster. In particular, the load-balancing algorithm is generalized as follows: an incoming job is assigned to server i with probability

x i j∈I x j for each i ∈ I if j∈I xj ≥ 1, otherwise the job is rejected. The corresponding closed Jackson network consists of N + 1 stations that correspond to the dispatcher and the N servers, respectively. The set of token classes is I and there are ℓi class-i tokens, for each i ∈ I.

Stationary distribution

The stationary distribution of the Markov process defined by the evolution of the network state x = (x1, x2, . . . , xN ) over time is given by

π(x) = 1 G(ℓ) x1 + x2 + . . . + xN x1, x2, . . . , xN N i=1 µi λ x i , x ≤ ℓ,
where

x 1 +x 2 +...+x N x 1 ,x 2 ,...,x N = (x 1 +x 2 +...+x N )! x 1 !x 2 !•••x N !
is a multinomial coefficient, and the constant G(ℓ) is obtained by normalization. For each i, j ∈ I and ℓ ∈ N N with ℓj ≥ 1, the variation of the normalization constant obtained by replacing a server-j slot with a server-i slot is denoted by δj→iG(ℓ) = G(ℓ + ei -ej) -G(ℓ). We now explain how to generalize our results.

Theorem 1 Following the same approach as in the proof of Theorem 1, we obtain the following generalization of [START_REF] Comte | Dynamic load balancing with tokens[END_REF]:

G(ℓ) = L ℓ1, ℓ2, . . . , ℓN N i=1 µi λ ℓ i + O λ→0 1 λ L-1
.

As λ tends to zero, the variations of the second term become negligible compared to those of the first. Therefore, an optimal slot allocation is a mode of the multinomial distribution with parameters L and µ1, µ2, . . . , µN .

Theorem 2 To generalize the proof of this theorem and the next, it is helpful to rewrite the stationary distribution as

π(x) = 1 G(ℓ) ϕ1(x1, x2)ϕ2(x1 + x2, x-1,2), x ≤ ℓ, (8)
where x-1,2 = (x3, . . . , xN ), and

ϕ1(x1, x2) = x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 , ϕ2(x1 + x2, x-1,2) = x1 + x2 + . . . + xN x1 + x2, x3, . . . , xN N i=3 µi λ x i .
The first factor is equal (up to a multiplicative constant) to the stationary distribution obtained in a two-server cluster, while the second factor depends on x1 and x2 only via their sum x1 + x2. Using this expression, we can rewrite G(ℓ), and then δ2→1G(ℓ), as a nested sum over x1 ∈ {0, 1, . . . , ℓ1} and x2 ∈ {0, 1, . . . , ℓ2}, where the term corresponding to x1 and x2 is the product of ϕ1(x1, x2) and a factor that depends on x1 and x2 only via their sum. Upon observing that this factor is Θ λ→0 (( 1 λ ) ℓ 3 +...+ℓ N ), we conclude in a similar way as in the proof of Theorem 1 that, when λ is small enough, we have δ2→1G(ℓ) < 0 if ℓ1 ≥ ℓ2 and δ2→1G(ℓ) > 0 if ℓ1 ≤ ℓ2 -1. A similar result holds for Gi→j (ℓ) for each i, j ∈ I such that ℓj ≥ 1, so that the loss probability is again minimized by choosing the buffer lengths approximately equal to each other.

Theorem 3

The monotonicity result of this theorem can only be generalized to the fastest and slowest servers.

More specifically, we can show that the optimal buffer length of the fastest server(s) decreases with the arrival rate, while the optimal buffer length of the slowest server(s) increases with the arrival rate. Indeed, all intermediary results in Appendices A and B can be generalized to two arbitrary servers i, j ∈ I using the same approach as in the previous paragraph. The only restriction is that µi ≥ µj (the assumption µi > µj is used in Case 1 in the proof of Lemma 7 in Appendix B, and one can verify that it can be alleviated to µi ≥ µj ). Therefore, we can only conclude for a server i ∈ I such that either µi ≥ µj for each j ∈ I \ {i} or µi ≤ µj for each j ∈ I \ {i}.

Numerical results

We now proceed to some case studies, in which we numerically evaluate the performance measures from Section 3.2 and show how the loss probability and mean response time are dependent on the arrival rate, service times and slot allocation.

Cluster of two servers

We first consider a cluster of two servers. The overall buffer length is kept constant equal to L = 20 in the interest of space, but we observed a similar behavior for other values of L. As before, we assume that µ1 + µ2 = 1 and µ1 > µ2.

Loss probability In Figure 3, the loss probability is shown as a function of the number of slots allocated to the fast server, for several values of the arrival rate λ, first in a linear plot and then in a log-linear plot. As intuition suggests, a larger arrival rate yields a larger loss probability. Since µ1 is large, the loss probability tends to be lower when the first server has more than half of the slots. The log-linear plot reveals that, even in this range, the loss probability can still be reduced by several orders of magnitude by correctly choosing the buffer lengths.

To make more decisive statements, we consider, in Figure 4, the optimal buffer length of the fast server (called the optimal buffer length for brevity) as a function of the arrival rate, for several values of the service rates. The optimal buffer length is always at least ⌈L/2⌉, as it does not make sense to allocate more slots to a slower server, and it increases with µ1.

When the arrival rate λ is small, the optimal buffer length is approximately µ1L, which is consistent with Theorem 1. Although it may at first seem that the larger the difference in server speeds, the later the optimal number of slots changes when increasing λ, this is a coincidence, as can be seen by comparing the cases µ1 = 0.9 and µ1 = 0.95. As λ increases, the optimal buffer length converges to L/2 = 10, as predicted by Theorem 2, but the convergence is slower when µ1 is larger. Figure 4 lastly shows that the optimal buffer length decreases as the arrival rate increases, as proven in Theorem 3.

Mean response time

We now turn to the mean response time, for which we only show numerical results. The mean response time is shown in Figure 5 as a function of the arrival rate λ, for several values of the buffer length of the first server. We first turn to the paradoxical behavior of the system when the majority of slots is allocated to the slowest server. In this case, the mean response time is not necessarily monotonous in the arrival rate λ, as can be seen with ℓ1 = 4. Indeed, when the arrival rate is low, most jobs are served by the slowest server; as the arrival rate increases, a larger fraction of jobs is sent to the fastest server (even if the buffer of this server is shorter), so that the mean response time decreases. Note that this scenario is suboptimal anyway, as switching around the slots (i.e. giving the majority of slots to the fast server) leads to better performance. As a side remark, the mean response times always converges to the same value as λ increases, as long as at least one slot is allocated to each server.

Figure 6 shows the optimal buffer length of the fast server (to minimize the mean response time) as a function of the arrival rate λ. The same monotonicity property as for the loss probability seems to hold. Furthermore, when the arrival rate is large, the optimal buffer length also seems converge to L/2. The main difference with the loss probability appears in the low-traffic regime: instead of converging to ⌈µ1L⌉, the optimal buffer length of the fastest server seems to converge to L.

Cluster of four servers

We finally consider a cluster with four servers to illustrate the absence of monotonicity of the optimal buffer length in larger clusters. The server speeds are µ1 = 0.45, µ2 = 0.3, µ3 = 0.2, and µ4 = 0.05 and overall buffer size L = 40. The optimal buffer lengths for each server, in terms of either the loss probability or the mean response time, are shown in Figure 7 and 8 as functions of the arrival rate.

As mentioned in Section 7, all theorems for the loss probability can be generalized to more than two servers: the optimal buffer lengths are proportional to the server speeds when the arrival rate is low and uniform when the arrival rate is large, and the optimal buffer length of the fastest and slowest servers evolve monotonically with the arrival rate. The optimal buffer lengths of servers 2 and 3 are however not monotonic. From extensive numerical experiments conducted for clusters with four to ten servers, we conjecture that, for clusters with N servers, the total optimal buffer size of the n = 1, 2, . . . , N fastest servers is decreasing, e.g., with N = 4, ℓ1, ℓ1 + ℓ2, and ℓ1 + ℓ2 + ℓ3 are decreasing in λ.

We observe in Figure 8 that when the arrival rate is low, the optimal buffer size of the fastest server (in terms of mean response time) is L, as in the two-server case. For a large arrival rate, the optimal buffer lengths are uniform. A similar sense of monotonicity is observed for the mean response time: the optimal buffer lengths of the fastest and slowest servers are monotonous. The total optimal buffer length of the n = 1, 2, . . . , L fastest servers also seems to be monotonous.

Conclusion

In this paper, we considered a load-balancing algorithm that leads to a product-form stationary distribution in clusters of servers with unequal service speeds. We developed analytical methods to understand the joint impact of the speeds and buffer lengths of the servers on performance. For a two-server cluster with arbitrary service speeds, we proved analytically that the optimal slot allocation in terms of the loss probability evolves monotonically from proportional to the server speeds to uniform as the arrival rate increases. We then generalized these results to clusters of more than two servers and assessed their validity on numerical examples. For the future works, we would like to generalize these analytical results to the other performance metrics mentioned in the introduction. The main difficulty is that their expressions involve fractions of two sums that both span all states. We would also like to generalize these results in other directions, for instance by considering open variants of the algorithm or by accounting for assignment constraints [START_REF] Comte | Dynamic load balancing with tokens[END_REF][START_REF] Weng | Optimal load balancing in bipartite graphs[END_REF].

Appendix

In the following two proofs, the arrival λ will be mentioned explicitly as a parameter of the quantities involved, so that in particular we will let G(λ, ℓ) denote the normalization constant under arrival rate λ and buffer length vector ℓ, and δG(λ, ℓ) = G(λ, ℓ + e1 -e2) -G(λ, ℓ). Also, it will be convenient to rewrite the normalization constant G(λ, ℓ) as

G(λ, ℓ) = x≤ℓ π(x), ℓ ∈ N 2 , ( 9
)
where π is a stationary measure such that π(0) = 1, that is,

π(x) = x1 + x2 x1 µ1 λ x 1 µ2 λ x 2 , x ∈ N 2 . ( 10
)
Injecting this definition into the definition of δG(ℓ) yields δG(λ, ℓ) = ℓ 2 -1

x 2 =0
π(ℓ1 + 1, x2) -

ℓ 1 x 1 =0 π(x1, ℓ2). (11)
Appendices A and B give the proofs of Propositions 4 and 5, respectively.

A Proof of Proposition 4 Consider an arrival rate λ > 0. We will prove that, for each ℓ ∈ N 2 such that ℓ1 ≥ 1 and ℓ2 ≥ 1, β(λ, ℓ + e1 -e2) ≤ β(λ, ℓ) implies β(λ, ℓ) < β(λ, ℓ -e1 + e2), that is, δG(λ, ℓ) > 0 implies δG(λ, ℓ -e1 + e2) > 0. Proposition 4 then follows from an induction argument that we omit. The following lemma, which follows directly from [START_REF] Gardner | Scalable load balancing in the presence of heterogeneous servers[END_REF], will be useful. Let ℓ ∈ N 2 such that ℓ1 ≥ 1, ℓ2 ≥ 1, and δG(λ, ℓ) > 0. By applying [START_REF] Gardner | Product forms for FCFS queueing models with arbitrary server-job compatibilities: an overview[END_REF] to k = ℓ -e1 + e2 and removing a (positive) term from the first sum, we obtain δG(λ, k) > ℓ 2

x 2 =1 π(ℓ1, x2) -

ℓ 1 -1 x 1 =0
π(x1, ℓ2 + 1).

Applying Lemma 6 to each term, using the fact that x2 < ℓ2 + 1 in the first sum and x1 < ℓ1 in the second sum, and making a change of variable yields

ℓ2 + 1 ℓ1 + 1 µ1 µ2 δG(λ, k) > ℓ 2 -1 y 2 =0
π(ℓ1 + 1, y2) -

ℓ 1 y 1 =1
π(y1, ℓ2).

Lastly, we add a (positive) term in the last sum and apply [START_REF] Gardner | Product forms for FCFS queueing models with arbitrary server-job compatibilities: an overview[END_REF], so as to obtain ℓ2 + 1 ℓ1 + 1 µ1 µ2 δG(λ, k) > δG(ℓ).

Therefore, δG(λ, ℓ) > 0 implies that δG(λ, k) > 0.

Conclusion

The sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } satisfies Conditions (1), (2), or (3) of the proposition.

We now prove Proposition 5. Let ℓ ∈ N 2 such that ℓ2 ≥ 1. Using notation from Lemma 7, we distinguish two cases. If the sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } satisfies Conditions (1) or (2), the result is immediate because the sign of δG(λ, ℓ) does not depend on λ. Now assume that the sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } satisfies Condition (3). If n * = ℓ1 + ℓ2 and cn * = 0, the conclusion is the same as under Condition [START_REF] Ayesta | On a unifying product form framework for redundancy models[END_REF]. Otherwise, we study the variations of the function f : λ → δG(λ, ℓ) on (0, +∞). This function is infinitely differentiable on this interval and, for each λ > 0, we have Hence, if there is λ > 0 such that f (λ) ≤ 0, then f ′ (λ) < 0. Lemma 7 also implies lim λ→0 + f (λ) = +∞. Therefore, either f is positive on (0, +∞), or there is λ0 > 0 such that f is positive on (0, λ0) and negative on (λ0, +∞), with f (λ0) = 0. Both cases lead to the conclusion.

f ′ (λ) = -

Figure 2 :

 2 Figure 2: Jackson network associated with the cluster of Figure 1.

Figure 3 :

 3 Figure 3: Loss probability in a two-server cluster, as a function of the buffer length of the first server and for several values of the arrival rate, with L = 20, µ 1 = 0.9, and µ 2 = 0.1.

Figure 4 :

 4 Figure 4: Optimal buffer length of the fast server (to minimize the loss probability) in a two-server cluster, as a function of the arrival rate and for several values of the service rates, with L = 20.

Figure 5 :

 5 Figure 5: Mean response time in a two-server cluster, as a function of the arrival rate and for several values of the buffer length of the first server, with L = 20 and µ 1 = 0.75.

Figure 6 :

 6 Figure 6: Optimal buffer length of the fast server in a two-server cluster (to minimize the mean response time), as a function of the arrival rate and for several values of the service rates, with L = 20.

4 Figure 7 : 4 Figure 8 :

 4748 Figure 7: Optimal buffer length of each server (to minimize the loss probability) in a four-server cluster, as function of λ, with L = 40, µ 1 = 0.45, µ 2 = 0.3, µ 3 = 0.2, µ 4 = 0.05.

Proposition 4 .

 4 Let λ > 0 and ℓ ∈ N 2 such that ℓ1 ≥ 1 and ℓ2 ≥ 1. If β(λ, ℓ + e1 -e2) ≤ β(λ, ℓ), then β(λ, ℓ -xe1 + xe2) < β(λ, ℓ -(x + 1)e1 + (x + 1)e2) for each x ∈ {0, 1, 2, . . . , ℓ1 -1}.

Lemma 6 .

 6 For each x ∈ N 2 such that x2 ≥ 1, we have π(x) = x1 + 1 x2 µ2 µ1 π(x + e1 -e2).

ℓ 1 +ℓ 2 n=min(ℓ 1

 21 +1,ℓ 2 ) (n * -1))cn 1 λ n+1 , = -(n * -1)f (λ) -ℓ 1 +ℓ 2 n=n * (n -(n * -1))cn 1 λ n+1 , < -(n * -1)f (λ).

Acknowledgments

The research of Mark van der Boor is supported by the NWO Gravitation Networks grant 024.002.003. The authors thank Sem Borst for his valuable comments on an earlier draft of the paper, and in particular for suggesting the generalization of Theorem 1 to clusters with more than two servers. The authors are also grateful to Mor Harchol-Balter for an insightful discussion on Theorems 1 and 2.

B Proof of Proposition 5

Proposition 5. Let λ * > 0 and ℓ ∈ N 2 such that ℓ2 ≥ 1. If β(λ * , ℓ+e1 -e2) < β(λ * , ℓ), then β(λ, ℓ+e1 -e2) < β(λ, ℓ) for each λ ∈ (0, λ * ).

We will prove equivalently that, for each λ * > 0 and ℓ ∈ N 2 such that ℓ2 ≥ 1, δG(λ * , ℓ) > 0 implies δG(λ, ℓ) > 0 for each λ ∈ (0, λ * ). The following lemma will be useful.

The sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } depends on the buffer lengths ℓ1 and ℓ2 and service rates µ1 and µ2 but not on the arrival rate λ, and satisfies one of the following conditions:

(1) cn > 0 for each n ∈ {min(ℓ1 + 1, ℓ2), . . . , ℓ1 + ℓ2},

(2) cn < 0 for each n ∈ {min(ℓ1 + 1, ℓ2), . . . , ℓ1 + ℓ2}, or

(3) there is n * ∈ {min(ℓ1 + 1, ℓ2), . . . , ℓ1 + ℓ2} such that cn < 0 for each n ∈ {min(ℓ1 + 1, ℓ2), . . . , n * -1}, cn * ≥ 0, and cn > 0 for each n ∈ {n * + 1, . . . , ℓ1 + ℓ2}.

Proof of the lemma. As in the proof of Theorem 2, we use (4) and distinguish two cases depending on the values of ℓ1 and ℓ2.

Case 1 (ℓ 1 + 1 ≤ ℓ 2 ) We can rewrite (4) as [START_REF] Goldsztajn | Self-learning threshold-based load balancing[END_REF], where the sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } is given by (13a) in Figure 9. It follows immediately that cn > 0 for each n ∈ {ℓ1 + 1, . . . , ℓ2 -1}. Additionally, we verify that cn > 0 for each n ∈ {ℓ2, . . . , ℓ1 + ℓ2} by calculating the ratio between the first and second term of the subtraction and showing, thanks to the assumption µ1 > µ2, that this ratio is (strictly) larger than one. Therefore, the sequence satisfies Condition (1).

Case 2 (ℓ 1 ≥ ℓ 2 ) We can rewrite (4) as [START_REF] Goldsztajn | Self-learning threshold-based load balancing[END_REF], where the sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } is given by (13b) in Figure 9. It follows directly that cn < 0 for n ∈ {ℓ2, . . . , ℓ1}. For n ∈ {ℓ1 + 1, . . . , ℓ1 + ℓ2}, cn is of the same sign as gn, where

We conclude that the sequence {cn} n∈{min(ℓ 1 +1,ℓ 2 ),...,ℓ 1 +ℓ 2 } satisfies Conditions (2) or (3) by observing that g ℓ 1 +ℓ 2 +1 = 0, and then showing that the sequence {gn} n∈{ℓ 1 +1,...,ℓ 1 +ℓ 2 +1} is strictly concave in the sense that gn -gn+1 > gn+1 -gn+2 for each n ∈ {ℓ1 + 1, . . . , ℓ1 + ℓ2 -1}. The details of the calculation are omitted due to space constraints.

, n ∈ {ℓ 1 + 1, . . . , ℓ 1 + ℓ 2 }.

(13b)

Figure 9: Definitions of the sequence {c n } n∈{min(ℓ1+1,ℓ2),...,ℓ1+ℓ2} depending on the value of ℓ = (ℓ 1 , ℓ 2 ).


Sylvain Mayolle

Roger Soliva

Stéphane Dominguez

Christopher Wibberley

Yannick Caniven

Nonlinear fault damage zone scaling revealed through analog modeling

come

INTRODUCTION

Fracture damage zones around faults control many aspects of fluid flow properties and the seismogenic behavior of faults. Their width is a key parameter controlling the hydraulic behavior of the fault, energy dissipation during earthquake ruptures, and stress distribution (e.g., [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF][START_REF] Perrin | Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth[END_REF]. Several studies over the past decade have evidenced a linear scaling between fault displacement (D) and damage zone thickness (T) up to a scale limit of a few hundred meters [START_REF] Micarelli | Structural evolution and permeability of normal fault zones in highly porous carbonate rocks[END_REF][START_REF] Savage | Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones[END_REF][START_REF] Mitchell | Towards quantifying the matrix permeability of fault damage zones in low porosity rocks[END_REF][START_REF] Johri | A scaling law to characterize fault-damage zones at reservoir depths[END_REF][START_REF] Solum | Toward the creation of models to predict static and dynamic fault-seal potential in carbonates[END_REF][START_REF] Mayolle | Scaling of fault damage zones in carbonate rocks[END_REF][START_REF] Torabi | Fault zone architecture and its scaling laws: Where does the damage zone start and stop[END_REF], defining a threshold in the scaling. Although discussed in terms of seismogenic depth for earthquake-related damage [START_REF] Ampuero | Upper limit on damage zone thickness controlled by seismogenic depth[END_REF], the existence of such a threshold for faults is not yet supported by a physical explanation. In this study, we developed a rheologically consistent and scaled analog model of Earth's brittle crust to experimentally investigate the existence of such a threshold in D-T scaling and to identify the main controlling parameters. With these aims, we measured the strain field evolution, D-T scaling, and volumetric strain using high-resolution image correlation for models of different layer thickness.

MODEL SETUP

We present the two most illustrative tests from a series of 16 experiments designed to study reproducibility, the role of material composition, and the relevant boundary conditions (see the Supplemental Material1 ). These two experiments were performed in a sandbox device 500 mm long, 300 mm high, and 100 mm wide, equipped with two lubricated low-friction glass walls and containing two layers of granular materials (Fig. 1). The upper layer was composed of 97.5% dry sand of 150 μm mean grain size and 2.5% pumice powder of similar packing to ensure reproducibility. The static friction angle and the cohesion were 38° ± 1° and 185 ± 20 Pa, respectively, and dilation angle decreased with increasing plastic strain. These physical properties are analogous to those of upper-crustal rocks, and cohesion was scaled to satisfy a model-to-nature scaling of 1 cm equivalent to 0.5-1 km in nature (e.g., [START_REF] Schellart | Shear test results for cohesion and friction coefficients for different granular materials: Scaling implications for their usage in analogue modelling[END_REF][START_REF] Arzúa | Dilation in granite during servo-controlled triaxial strength tests[END_REF][START_REF] Reber | What model material to use? A review on rock analogs for structural geology and tectonics[END_REF]. Such materials consistently reproduce the main geometries of faults and their segmentation, interaction, displacement profiles, and gradients, but also displacement-length and size distribution scaling laws (e.g., [START_REF] Bellahsen | Fault reactivation control on normal fault growth: an experimental study[END_REF][START_REF] Schlagenhauf | Incremental growth of normal faults: Insights from a laser-equipped analog experiment[END_REF]. The lower layer was a mixture of 98% dry sand of 150 μm mean grain size and 2% polydimethylsiloxane silicone (kinetic sand), inducing a visco-plastic mechanical behavior simulating that of the lower crust or ductile shales. These two layers lay on a pre-compressed elastic polyurethane foam which, once relaxed by movement of a lateral backstop, induces distributed extensional strain in the layers above (Fig. 1). The extension velocity (30 mm/h) was regulated by a computer-controlled motor and allowed a total extension of 100 mm (20%). We performed two sets of experiments with different frictional upper-layer thicknesses: 50 mm and 200 mm (i.e., ∼4 km and ∼15 km thickness in nature). The kinetic sand thickness of the lower layer was 20 mm for the thinner model and 40 mm for the thicker model. All the experiments were monitored through the lateral glass wall using a Sony α7R II camera taking one photo every 30 s with a resolution of 1 pixel = 54 μm.

The acquired images were processed to analyze the displacement field using a subpixel spectral correlation algorithm from [START_REF] Van Puymbroeck | Measuring earthquakes from optical satellite images[END_REF] (Fig. 1) and were also used to provide movies of the strain field *E-mail: sylvain.mayolle@umontpellier.fr (see the Supplemental Material). The cumulative displacement field (U; i.e., the sum of displacement fields derived from each image correlation) and cumulative net shear displacement along the faults (D) were calculated using the horizontal and vertical displacement fields. Cumulative damage zone thickness (T) was determined using the distribution of displacement gradient across the fault (dU/dx; i.e., grdgradient function in Generic Mapping Tools software (GMT, https://www.genericmapping-tools.org/), the amplitude of the gradient) (Fig. 2A). This displacement gradient field is the most relevant field attribute for measuring fault thickness because it accounts for both the absolute displacement (any direction) and displacement modes (shear, dilation, and compaction). For each model, strain evolution is analyzed along six profiles normal to two selected fault planes. The Gaussian distribution of the displacement gradient along these profiles is compared to background deformation outside the Gaussian curve to define damage zone thickness (Fig. 2B). In an area devoid of faults, an upper limit of the background displacement gradient is used at each stage of image correlation, defined as 3σ of the modal value of the background displacement gradient (horizontal dashed line in Fig. 2B). The distance between the two points where the Gaussian distribution reaches this upper background value defines T, the cumulative damage zone thickness (Fig. 2B, stage 2). Thus, T is measured in relation to background deformation as is commonly done along scanlines from field data (e.g., [START_REF] Beach | Reservoir damage around faults; Outcrop examples from the Suez Rift[END_REF][START_REF] Choi | Definition and classification of fault damage zones: A review and a new methodological approach[END_REF][START_REF] Ostermeijer | Damage zone heterogeneity on seismogenic faults in crystalline rock: A field study of the Borrego Fault[END_REF]. Areas perturbed by boundary conditions are excluded from measurements (see the Supplemental Material).

We also measure the incremental volumetric strain (dε v ) along the faults to investigate its evolution in relation to the scaling behavior. To quantify its evolution relative to the incremental shear displacement during fault growth, we also calculate its derivative in the form, dε v /dD.

RESULTS

Fault Geometry

Experimental faults initiate as shear zones dipping ∼60° toward the moving backstop, and cumulate, at the end of the experiment, a net shear displacement of as much as 4 mm for the model with 50 mm upper-layer thickness (8%) and as much as 7 mm for the model with 200 mm upper-layer thickness (3.5%) (Fig. 2A). They die out at depth by diffuse deformation into the visco-plastic kinetic sand layer. In the 200 mm model, individual dilational-shear segments are observed in the very early stages of shear zone development, linking up as the faults grow, leading to irregular fault surface geometries (Fig. 2C, stages 1 and 2). During the experiment, all the faults grow wider as displacement increases (Figs. 2A and2B, stages 1 and 2). With increasing displacement, the fault deformation concentrates on narrow zones, which become straight and continuous through the sand layer (Fig. 2C, transition from stage 2 to stage 3). This strain localization is also shown by the change in shape of the Gaussian curve of displacement gradient profiles across the faults, which increases its amplitude faster than its wavelength (Fig. 2B).

The incremental volumetric strain field reveals that in the frictional layer, the first stages of fault thickening occur under a dominant dilational-shear mode (Fig. 2C, stage 1). The deriva-tive of this incremental volumetric strain with respect to the net shear displacement (dε v /dD) increases but never exceeds 0.15. Fault widening continues to evolve by a combination of dilational-shear and compactional-shear modes (Fig. 2C, stage 2). Dilational shear occurs in a wide and diffuse volume of the fault, especially at the steeper parts of the fault plane. Compactional shear initiates in limited regions of the fault, especially at contractional bends along the irregular fault surface (Fig. 2C, stage 2). Following this stage, a very thin zone of compactionalshear deformation, visible in the incremental volumetric strain field, develops over a larger extent along the fault and can reach a value of -0.14 (Fig. 2C, stage 3). Although its geometry is variable from one fault to another, this very thin compactional-shear zone is surrounded by synchronous dilational-shear deformation, spreading away from the fault at specific locations such as zones of segment linkage (Fig. 2B, stage 3). The cumulative effect of these two synchronous deformation modes results in net dε v /dD values ranging from -0.05 to 0.05 over the fault width, commonly close to 0 because incremental compaction and dilation are similar in magnitude.

D-T Scaling

Figure 3A shows a plot of D versus T measured in the two sets of experiments along profiles both across fault segments with simple geometries (i.e., no intersections, no branching) and across more complex zones such as fault branches and relay zones. The dε v /dD ratio is indicated for each stage measured. This diagram reveals that damage zone thickness grows nonlinearly with increasing fault displacement (Figs. 3A and3B). The first stage of fault growth shows a nonlinear rapid increase of fault thickness with fault displacement, in which the fault develops in a hybrid dilationalshear failure mode (blue symbols in Figs. 3A and3B). In both models, this nonlinear increase occurs up to a value of D = ∼0.8 mm, where the dT/dD ratio stabilizes around a very low mean value of 0.18 (with 2σ = 0.6), reflecting very little increase of damage zone thickness with displacement (thickness threshold). This stabilization of the dT/dD ratio occurs nearly synchronously with a change of failure mode from dominant dilational-shear behavior (dark blue symbols) to isochoric shear (white, light blue, or pink symbols). The dε v /dD value then remains nearly stable until the end. Note that this behavior is observed in each of the different model layer thicknesses, but the damage zone thickness at which the fault stabilizes is ∼20% higher for the thicker model in cases of simple fault segments. In the thicker model, T values measured at relay zones are ∼30% larger than along segments. At relay zones, the main curvature of the D-T data trend is also smoother (Figs. 2A and2B, green profile; Fig. 3A), showing a later and more gradual transition of failure mode from dilational shear to isochoric shear. Finally, due to this thickness threshold, D becomes larger than T in the later stages of the models (see the Supplemental Material), consistent with natural faults when D > 100 m.

DISCUSSION

Our experiments show that ( 1) nonlinear D-T scaling occurs toward a plateau at a threshold value of damage zone thickness dependent on layer thickness, (2) a change of dominant failure mode occurs when the dT/dD ratio stabilizes toward a value close to zero, and (3) the threshold damage zone thickness increases significantly at fault segment relay zones, for which the width is large in the thicker model. This latter effect promotes a delay in the stabilization of the dT/dD ratio and allows a more progressive transition in failure mode.

A hybrid dilational-shear mode is prominent in the first stage of fault growth during fault segment initiation and propagation. This

B C

A dilation is preserved later in extensional relay zones and around the fault as shear localization occurs (Fig. 2C). Hybrid dilational-shear failure modes are recorded in low-cohesion granular materials as shear-enhanced dilation due to grain disaggregation and rotation during shear (e.g., [START_REF] Rowe | The stress-dilatancy relation for static equilibrium of an assembly of particles in contact[END_REF][START_REF] Antonellini | Microstructure of deformation bands in porous sandstones at Arches National Park, Utah[END_REF]. Although the granular material does not produce transgranular fractures, it reproduces the main features of fault network geometry such as tip propagation, segmentation, linkage at relays, and fault bends, all of which contribute to generating the deformation volume around the fault as observed in nature irrespective of the deformation micro-mechanisms (e.g., [START_REF] Childs | A geometric model of fault zone and fault rock thickness variations[END_REF][START_REF] Wibberley | Recent advances in the understanding of fault zone internal structure: A review[END_REF]Figure 2C; also see the Supplemental Material). This is because the mechanical behavior is consistent with stiffer rocks, and the grain size in our models is much lower than the fault thickness at saturation (by one to two orders of magnitude). For the same reason, initial grain-scale heterogeneity prob-ably has little impact (<10%) on measured threshold thickness. The final dilational-shear thickness depends on initial fault-segment spacing and relay zone width, which in turn depend on frictional layer thickness, probably because the distance at which segments can interact increases with layer thickness [START_REF] Soliva | Spacing and linkage of confined normal faults: Importance of mechanical thickness[END_REF]. Consequently, larger relay zones provide larger damage volumes, themselves needing higher displacements to evolve into smooth isochoric shear surfaces.

When segment linkage is mature, the fault clearly cuts the entire thickness of the frictional layer and compactional-shear initiates on favorably oriented parts of the irregularities related to the earlier segmentation [START_REF] Mayolle | Scaling of fault damage zones in carbonate rocks[END_REF]. In a latter stage, compaction develops through the entire fault, probably due to the low strength of the dilated zones [START_REF] Wong | Mechanical compaction of porous sandstone: Oil & Gas Science and Technology[END_REF][START_REF] Rutter | The deformation of porous sandstones: Are Byerlee friction and the critical state line equivalent?[END_REF], which helps suppress further thickening. This process reflects strain localization and failure of the frictional layer and progressively inhibits dilational shear to promote isochoric shear within a fault of constant thickness.

COMPARISON TO NATURE AND CONCLUSION

The models show first-order similarities with natural data such as fault dip angles, conjugate geometries, segmentation, and D-T trend. Although precise scaling cannot be achieved, the ratio of layer thickness to mean threshold thickness found in the thick model (200 mm / 4.5 mm) is close to the ratio expected for natural faults in the seismogenic crust (15 km / 0.3 km; [START_REF] Savage | Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones[END_REF][START_REF] Solum | Toward the creation of models to predict static and dynamic fault-seal potential in carbonates[END_REF]. Also, the observed volumetric strain can be considered, to a certain extent, analogous to fracture dilation and compaction occurring within faults (mode I fractures, dilatant breccias versus cataclasites, gouge smearing, S-C structures, cleavage, and pressure-solution seams; e.g., [START_REF] Billi | The damage zone-fault core transition in carbonate rocks: Implications for fault growth, structure and permeability[END_REF][START_REF] Gratier | Aseismic sliding of active faults by pressure solution creep: Evidence from the San Andreas Fault Observatory at Depth[END_REF][START_REF] Ishii | Far-field stress dependency of the failure mode of damage-zone fractures in fault zones: Results from laboratory tests and field observations of siliceous mudstone[END_REF][START_REF] Vrolijk | Clay smear: Review of mechanisms and applications[END_REF]. In addition, it is well known that mode I fractures generally develop early in the growth history of a fault in natural rocks (e.g., [START_REF] Petit | Crack-seal', slip: A new fault valve mechanism?[END_REF][START_REF] Crider | Initiation of brittle faults in the upper crust: A review of field observations[END_REF]. Larger faults generally show a combination of preferentially dilational-shear damage zones and compactional-shear core zones (e.g., [START_REF] Wibberley | Internal structure and permeability of major strikeslip fault zones: The Median Tectonic Line in Mie Prefecture[END_REF][START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF][START_REF] Holdsworth | Fault rocks from the SAFOD core samples: Implications for weakening at shallow depths along the San Andreas Fault[END_REF], with an average spatial arrangement consistent with the volumetric strain distribution presented in this study.

Finally, as the overall mechanical controls on the failure-mode transition govern the threshold of damage zone growth, this threshold is likely to be specific to the geological context of any one fault, probably accounting for the large scatter observed in the natural data [START_REF] Savage | Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones[END_REF]. However, this study suggests that mechanical layering is probably one of the most important parameters governing this threshold and scattering in nature, and such layering exists at several scales, with the ultimate scale limit of layer thickness controlling fault damage zone width being the thickness of the brittle upper crust. Other parameters like lithology, tectonic setting (e.g., [START_REF] Balsamo | Tectonic control on vein attributes and deformation intensity in fault damage zones affecting Natih platform carbonates, Jabal Qusaybah[END_REF], and the physico-chemical behavior of fluids in the fault are probably also important for scattering. For this reason, the specific context of any one fault must be taken into account when comparing different data sets for establishing overall scaling laws.

ACKNOWLEDGMENTS

TOTAL SA and Geoscience Montpellier provided financial support of this research developed during the DAMAGE research program hosted at the University of Montpellier. We thank C. Romano for technical assistance. We thank A. Braathen, F. Balsamo, R Schultz, and anonymous reviewers for their helpful reviews.

Figure 1

 1 Figure 1. Experimental setup. At the base of the model, elastic foam is uncompacted by a computerized motor, imposing extension in the upper visco-plastic kinetic sand and frictional layer. The vertical displacement and vertical strain field derived from sub-pixel image correlation are shown on the right-hand side.

Figure 2

 2 Figure 2. Example of fault growth in the experiment with a 200-mm-thick upper layer at three different stages. (A) Displacement gradient field. Damage zone thickness is measured along five profiles crossing simple fault geometry (in black) and one profile in the relay zone (green). (B) Displacement gradient distribution along profiles shown in A. Gray dashed line is upper limit of background deformation. (C) Incremental volumetric strain in the black dashed rectangle shown in A. Pixels represent image correlation resolution. Strain field movies are presented in the Supplemental Material (see footnote 1).

Figure

  Figure 3. (A) Displacement (D)-damage zone thickness (T) data from our experiments for two different frictional layer thicknesses (200 mm and 50 mm) measured along 24 profiles crossing four selected faults. Three examples of segment interactions are shown with corresponding displacement gradient fields. Green dashed lines show positions of measurement profiles. (B) dT/dD versus displacement. Symbol color scale is dε v /dD in both graphs, where dε v is incremental volumetric strain.

Supplemental Material.[START_REF] Holdsworth | Fault rocks from the SAFOD core samples: Implications for weakening at shallow depths along the San Andreas Fault[END_REF] Movie and description of displacement gradient field evolution; (2) movie of strain fields evolution; (3) cummulative displacement field; (4) displacement gradient profiles through faults; (5) volumetric strain profiles through faults; (6) mechanical properties of the frictional layer; (7) images and parametric data of all the models; (8) a graph of all the measured D-T data; and additional information on the method and scaling to nature. Please visit https://doi .org/10.1130/XXXXX to access the supplemental material, and contact editing@geosociety.org with any questions.