1. Introduction
In the world of cloud-based computing, large data centers are often used for file storage. In these data centers, large networks of servers are used to store even larger sets of files. In order to improve reliability and retrieval speed, these files are often ‘coded’. By coded, we mean that the file is broken down into smaller pieces which are stored on multiple servers. Consider the situation in which there are four servers and one file. One can store the entire file on one server, but in such a configuration the file would be inaccessible if that server were to fail. In order to improve reliability, one can replicate the file across all four servers, but such a method would require much more memory. Suppose instead that we split the file into halves, A and B, and then store A, B, A + B, A − B in each of the four servers, respectively. Then the original file can be constructed from any two pieces. One can extend this idea to the case where equally sized pieces of a file are stored across L servers and any k pieces can reconstruct the original file. This can be accomplished using the maximum distance separable (MDS) code with parameters (L, k) [Reference Lin and Costello29]. The MDS code greatly improves reliability since L − k + 1 servers must fail before the file becomes irretrievable, while only requiring enough total memory to store L/k files. Given a coding scheme, one can consider load balancing mechanisms to improve file retrieval speed. In [Reference Li, Ramamoorthy and Srikant28], two routeing schemes, called batch sampling (BS) and redundant request with killing (RRK), were considered. In BS routeing, incoming jobs are routed to the k-shortest queues containing the file being requested, while in RRK routeing jobs are routed to all servers containing the requested file and then removed from the queue (killed) once k pieces of the file have been returned. Li [Reference Li, Ramamoorthy and Srikant28] formally calculated the steady state (T → ∞) queue length distribution in the large system limit (n → ∞) and gave simulation results for different values of L and k in both routeing schemes.
In this work we are interested in developing a rigorous limit theory for such load balancing schemes for systems with MDS coding as n becomes large. Specifically, we establish law of large numbers and diffusion approximations for such systems under an appropriate scaling, as n → ∞. Such limit theorems provide useful model simplifications that can then be employed for approximate simulation of the large and complex n-server systems (see Section 6 for some numerical results). These limit theorems are also the first steps towards making rigorous the program initiated in [Reference Li, Ramamoorthy and Srikant28] of developing steady state approximations for such systems, with provable convergence properties as n becomes large.
We will focus in this work only on BS routeing and leave analysis of the RRK scheme for future work. Specifically, we consider a system with n servers on which I(n) files are stored using MDS coding with parameters (L, k). A key assumption to our analysis is that the files are stored such that each combination of L servers has exactly c files. We further assume that jobs arrive in the system according to a Poisson process with rate nλ and request a file uniformly at random. This is another simplifying assumption on our model that roughly says that all files are in equal demand. These structural assumptions imply a convenient exchangeability property of the system which allows for the use of certain mean-field approximation techniques. A single file request spawns k jobs which are then routed into the k-shortest queues within the set of L servers containing the file being requested. Each server processes the jobs in their queue with exponential service rate k according to the first-in–first-out (FIFO) discipline and processing times are mutually independent. Regarding each server as a ‘particle’, the above formulation describes an interacting particle system with simultaneous jumps. Note that the symmetry structure introduced above implies that every time a file request arrives, it leads to a selection of L servers uniformly at random (from which the k servers with shortest queues are chosen). In particular, this says that the well-studied ‘power-of-d’ routeing scheme (also known as the ‘supermarket model’) is a special case of the scheme considered here on taking L = d and k = 1. Direct analysis of such large and complex n-server systems is challenging even by simulation methods as frequently the servers in networks of interest number in the hundreds of thousands with arrival rates of file requests of similar order. The goal of this work is to develop suitable approximate approaches to such systems.
The starting point of our analysis is to consider, as the state descriptor, the empirical measure of the n queue lengths rather than the individual values of the queue lengths. Thus, the state space for our system will be the space 
$\mathcal{P}_n(\mathbb{N}_0)$ of probability measures on ℕ0 that assign weights in ℕ0/n to sets in ℕ0 rather than the space 
$\mathbb R_+^n$. With this formulation, the state processes for all n-server systems can be regarded as taking values in a common space 
$\mathcal S\,{:{=}}\,\mathcal{P}(\mathbb{N}_0)$ (the space of probability measures on ℕ0). It follows from our symmetry assumptions that the state evolution of the n-server system describes a pure-jump Markov process with values in 
${\mathcal{P}}(\mathbb{N}_0)$ and thus one can bring to bear the theory of weak convergence of Markov processes to study the scaling limits as n becomes large. In particular, in Theorem 1 we prove a law of large numbers for the empirical measure process (πn(t))0≤t≤T as n → ∞. We show that πn converges to a deterministic function π in 𝔻([0, T]: 𝒮), where 𝔻([0, T]: 𝒮) is the space of functions from [0, T] to 𝒮 that are right continuous and have left limits, equipped with the usual Skorokhod topology. Next we consider the fluctuation process 
$X^n\,{:{=}}\,\sqrt{n}(\pi^n\,{-}\,\pi)$. This process can be regarded as taking values in the space of signed measures on ℕ0; however, for an asymptotics analysis, it is convenient to view it as taking values in the Hilbert space of square summable real sequences, ℓ 2. The study of the asymptotics of these fluctuations is the subject of Theorem 2, which shows that 
$X^n \,{:{=}}\, \sqrt{n}(\pi^n -\pi)$ converges in 𝔻([0, T]: ℓ 2) to an ℓ 2-valued diffusion process.
Limit theorems of the form studied in this work can be used for model simplification and for computing approximations for performance measures, e.g. through simulation methods. Direct simulation of the underlying n-server system would in general be prohibitively expensive for large n since the jumps in the system occur at a rate proportional to n. The asymptotic approximations given in this work allow a system manager to simulate performance metrics for the system at a coarser scale via numerical ordinary differential equation (ODE) solvers or Euler discretizations for stochastic differential equations (SDEs) (see Section 6 for an example). Although the systems considered here are required to satisfy certain symmetry conditions (all files are equally sized and all jobs are in equal demand), the simplified models given by the limiting ODE and SDE give useful qualitative insights into the behavior of large storage networks employing these types of coding schemes. The results obtained here are useful in analyzing the long-time behavior of such systems as well, e.g. providing information on the rate at which the queue lengths decay in steady state and how such a decay is impacted by different values of L and k [Reference Friedlander15]. Furthermore, techniques developed in this work can also be used for models satisfying weaker symmetry conditions (e.g. for multitype populations with a fixed finite number of types). The poisson representations used in the proofs of Theorems 1 and 2 crucially rely on the fact that interarrival and service times are exponentially distributed. Different methods will need to be used in order to treat the case of general distributions. We refer the reader to [Reference Bramson, Lu and Prabhakar6] and Section 6 of [Reference Stolyar35] for work in this direction.
Load balancing mechanisms similar to the type considered here have been studied in many works. Specifically, the join-the-shortest-queue (JSQ), join-the-idle-queue (JIQ), and power-of-d routeing schemes have garnered quite a bit of attention (see [Reference Bramson, Lu and Prabhakar6], [Reference Eschenfeldt and Gamarnik13], [Reference Graham16], [Reference Mitzenmacher31], [Reference Mukherjee, Borst, van Leeuwaarden and WHITING32], [Reference Stolyar35], [Reference Vvedenskaya, Dobrushin and Karpelevich37], and the references therein). Mitzenmacher [Reference Mitzenmacher31] and Vvedenskaya et al. [Reference Vvedenskaya, Dobrushin and Karpelevich37] first analyzed the power-of-d routeing scheme, showing that in steady state the fraction of queues with lengths exceeding m decay super exponentially in m, a large improvement over the exponential rate for the setting where jobs are routed to servers uniformly at random. Graham [Reference Graham16] established a functional law of large numbers for πn on 𝔻([0, T]: 𝒮) in the power-of-d routeing scheme using characterization results for nonlinear martingale problems. In [Reference Eschenfeldt and Gamarnik13], the authors derived a diffusion approximation for the JSQ routeing policy in the large-system limit under heavy-traffic scaling. It was shown that the limit can be characterized through a two-dimensional diffusion. In [Reference Mukherjee, Borst, van Leeuwaarden and WHITING32], it was shown that the JIQ routeing scheme yields the same diffusion approximation as the JSQ routeing scheme. In both these works, the diffusion approximations were derived under the same scaling regime as considered here. However, unlike for the JSQ and JIQ routeing regimes where the diffusion approximations can be described through two-dimensional processes, in the current work the limit is an infinite-dimensional diffusion described as an ℓ 2-valued process driven by a cylindrical Brownian motion. We refer the reader to [Reference Aghajani and Ramanan1], [Reference Decreusefond and Moyal12], [Reference Kaspi and Ramanan23], and [Reference Reed and Talreja33] for other queueing network systems where infinite-dimensional diffusions arise as approximate models. As noted earlier, the power-of-d scheme is a special case of our results, and, thus, Theorems 1 and 2 also provide law of large numbers and diffusion approximations for this classical load balancing scheme (see Corollaries 1 and 2). In particular, Corollary 1 recovers the law of large numbers established in [Reference Graham16]. Limit theorems giving fluctuation results for power-of-d schemes have not been studied previously.
Diffusion approximation methods have been used extensively in stochastic network theory. In particular, they have been very useful in the study of critically loaded stochastic processing networks (see [Reference Atar and Shifrin4], [Reference Bell and Williams5], [Reference Budhiraja and Ghosh8], [Reference Budhiraja, Ghosh and Lee9], [Reference Dai and Lin11], [Reference Harrison17], [Reference Kushner27], [Reference Whitt38], and the references therein). For such systems, the key mathematical tool is the functional central limit theorem for renewal processes which provides Brownian motion approximations for a finite collection of centered renewal processes with rates approaching ∞. The scaling regime and mathematical tools that are relevant for the analysis in the current context are quite different from those used in the above works. In particular, here the number of nodes approach ∞ and the tools for proving convergence come from martingale problems and Markov process theory. A similar scaling regime was considered in [Reference Budhiraja and Friedlander7] for certain systems motivated by ad-hoc wireless network models introduced in [Reference Antunes, Fricker, Robert and Tibi3]. A key simplifying feature there is that the state space of an individual queue is a finite set. Consequently, the limit diffusion in [Reference Budhiraja and Friedlander7] is finite dimensional and, thus, for diffusion approximations, classical convergence theorems from [Reference Joffe and Métivier20] and [Reference Kurtz24] can be invoked. In contrast, the queue length processes in this work are unbounded, taking values in ℕ0, and thus one needs to study diffusion approximations in an infinite-dimensional state space, namely the Hilbert space ℓ 2. The proofs employ appropriate criteria for tightness and characterization results for Hilbert space-valued stochastic processes.
A basic assumption in our analysis of the fluctuations around the law of large number limit (see statement of Theorem 2) is a uniform (in n) bound on the second moment of the empirical measure at time 0. This condition is not very stringent and allows for many types of initial configurations, e.g. one where no queue contains more than k jobs (where k is independent of n). We argue that these integrability properties at time 0 propagate through to any finite future time T. Tightness of the scaled fluctuation processes Xn which is shown by establishing, uniform in n, second moment bounds (on Xn) and by employing criteria for tightness of Hilbert space-valued semimartingales (cf. [Reference Joffe and Métivier20] and [Reference Métivier30]), relies on these integrability properties. Another ingredient in the proof of tightness is a suitable Lipschitz property of the map F introduced in (5) that enables the use of a Gronwall argument. For this argument, one needs a Lipschitz estimate in the ℓ 2 norm; however, it is not clear that F, as a map from ℓ 2 to ℓ 2, is Lipschitz. We instead restrict attention to a smaller space
\begin{align*}
\nu_M\ \!{:\!=} \Bigg\{r \in \ell_2\colon r_i\ge 0,\, \sum_{i=0}^\infty r_i=1,\,
\sum_{i=0}^\infty i r_i \le M\Bigg\},
\end{align*}
and argue that, for each M, the map F is Lipschitz from 𝒱M to ℓ 2. This ‘local’ Lipschitz property plays an important role in the proof of Proposition 4.
For characterization of limit points in the proof of the central limit theorem, one needs to argue that the associated SDE in ℓ 2 (see (11)) has a unique weak solution in an appropriate class of processes. It turns out that arguing this uniqueness among adapted processes with paths in ([0, T]: ℓ 2) (the space of continuous functions from ℂ[0, T] to ℓ 2) is not straightforward due to a lack of suitable regularity of the function G introduced in (16). In particular, once more, the Lipschitz property of the map x → G(x, π) (for a fixed 
$\pi \in\mathcal{P}(\mathbb{N}_0)$) from ℓ 2 to itself is not immediate. The key observation here is that this map is Lipschitz when restricted to the space
\begin{align*}
\tilde \ell_2 {:\!=}\Bigg\{x\in\ell_2\colon\sum_{j=0}^\infty j^2x_j^2<\infty,\,\sum_{j=0}^\infty x_j = 0\Bigg\}. \end{align*}
This observation, together with the property that the limit points X of 
$X^n = \sqrt{n}(\pi^n-\pi)$ satisfy X(t ∈ ℓ̃ 2) for all t ≥ 0 almost surely, is key to the characterization of the limit points as the unique solution of SDE (11) in a suitable class (see Proposition 2).
The paper is organized as follows. In Section 2 we give a precise mathematical formulation of our model and a statement of our main results. Specifically, Theorem 1 provides the convergence in probability of the empirical measure process in 𝔻([0, T]: 𝒮) to the unique solution of the ODE defined in (7). In Theorem 2 we present the main diffusion approximation result. This result says that the sequence of centered and scaled processes Xn, defined in (10), converges to the unique solution (in a suitable class) of the ℓ 2-valued SDE, driven by a cylindrical Brownian motion, given in (11). In Section 2.1 we record the corollaries of these results for the special setting of power-of-d schemes. We then proceed to the proofs of Theorems 1 and 2. In Section 3 we give a convenient representation of the state processes through a countable number of time-changed unit-rate Poisson processes. Such Poisson representations have been used extensively (cf. [Reference Anderson and Kurtz2], [Reference Kang, Kurtz and Popovic22], and [Reference Kurtz25]) in the study of diffusion approximations for pure-jump processes. Using this, we obtain a semimartingale decomposition (see (20)) that is central to our analysis. Section 4 is devoted to the proof of Theorem 1. In Section 4.1 we prove tightness of the sequence of state processes {πn}n∈ℕ (see Proposition 3) and the proof of Theorem 1 is completed in Section 4.2. In Section 5 we prove Theorem 2. In Section 5.1 we prove the propagation of integrability properties that was discussed earlier, and in Section 5.2 (see Proposition 4) we prove the key tightness property for the sequence of processes {Xn}n∈ℕ which relies on the Lipschitz property of F, in the ℓ 2 norm, on 𝒱M (Lemma 4). Theorem 2 is then proved in Section 5.3. Finally, in Section 6 we present some numerical results. In particular, we use our results to give numerical confidence intervals for several performance measures of interest and compare the results to those obtained from a direct simulation of the corresponding n-server systems.
1.1. Notation
The following notation will be used. Fix T < ∞. All stochastic processes will be considered over the time horizon [0, T]. We will use the notation (X t)0≤t≤T and (X(t))0≤t≤T interchangeably for stochastic processes. The space of probability measures on a Polish space 
$\mathbb S$, equipped with the topology of weak convergence, will be denoted by 
$\mathcal P(\mathbb S)$. When 𝕊 = ℕ0, we will metrize 
$\mathcal P(\mathbb S)$ with the metric d 0 defined as
\begin{align*}
d_0(\mu, \nu) \,{:{=}}\, \sum_{j=0}^{\infty} \frac{|\mu(j)-\nu(j)|}{2^j}, \qquad \mu, \nu \in \mathcal P(\mathbb N_0). \end{align*}
For 𝕊-valued random variables X, X n, n ≥ 1, convergence in distribution of X n to X as n → ∞ will be denoted as X n ⇒ X. The space of functions that are right continuous with left limits (RCLL) from [0, T] to 𝕊 will be denoted as 𝔻([0, T]: 𝕊) and equipped with the usual Skorokhod topology. Similarly ℂ([0, T]: 𝕊) will be the space of continuous functions from [0, T] to 𝕊, equipped with the uniform topology. We will usually denote by κ, κ 1, κ 2, … the constants that appear in various estimates within a proof. The values of these constants may change from one proof to another and, unless stated otherwise, will take values in the set (0, ∞). Let 
${\ell _2} = \{ ({a_j})_{j = 0}^\infty \mid \sum {_{j = 0}^\infty a_j^2} < \infty \}$ be the space of square summable real sequences. This space is a Hilbert space with inner product
\begin{align*}
\left\langle x,y\right\rangle_2 = \sum_{j=0}^\infty x_jy_j.
\end{align*}
We denote the corresponding norm as ‖ · ‖2. Similarly, 
${\ell _1} = \{ ({a_j})_{j = 0}^\infty |\sum\limits_{j = 0}^\infty | {a_j}|\, < \infty \}$ and ‖ · ‖1 is the norm on this Banach space. The Hilbert–Schmidt norm of a Hilbert–Schmidt operator A on ℓ 2 will be denoted by ‖A‖ HS (cf. Appendix A.3). We denote by I the identity operator. For a Hilbert Space 
$\mathbb {H}, \mathcal M^2_T(\mathbb H)$ will denote the space of all ℍ-valued continuous, square-integrable martingales M, such that M(0) = 0. For a real number a, (a)+ will denote the positive part of a.
2. Model description and main result
We consider a system with n servers each with its own infinite capacity queue. In all, there are I(n) equally sized files stored over the n servers. Each file is stored in equally sized pieces at L servers such that any k pieces can reconstruct the original file. The files are distributed such that each combination of L servers has exactly c files. This, in particular, implies that 
$I(n)=c\binom{n}{L}$. Jobs arrive from outside according to a Poisson process with rate nλ and request one of the I(n) files uniformly at random. Such a request corresponds to selection of one of the 
$\left({_L^n} \right)$ sets of L servers, uniformly at random, which is the set of servers containing the pieces of the requested file. The job is then routed to the k-shortest queues among this set of L servers. Each server processes queued jobs according to the FIFO discipline. Processing times at each server are mutually independent and exponentially distributed with mean k −1.
Let 
$Q^n(t)=\{Q^n_i(t)\}_{i=1}^n$, where 
$Q^n_i(t)$ represents the length of the ith queue at time t, and let 
$\pi^n(t)=\{\pi^n_i(t)\}_{i\in\mathbb N_0}$, where 
$\pi_i^n(t)$ represents the proportion of queues with length exactly i at time t. This can explicitly be written as
\begin{align}
\pi^n_i(t)=\frac{1}{n}\sum_{j=1}^n{\bf 1}_{\{Q^n_j(t)=i\}}.\label{eq:eq102}
\end{align}
We will assume for simplicity that Qn(0) = qn is nonrandom and, thus, 
${\pi ^n}(0) = (1/n)\sum {_{j = 1}^n} {{\bf{1}}_{\{ q_j^n = i\} }}$ is nonrandom as well. We identify 
$\mathcal P(\mathbb N_0)$ with the infinite-dimensional simplex 
${\cal S} = {\rm{\{ }}s \in\mathbb R_ + ^\infty \mid \sum {_{j = 0}^\infty } {s_i} = 1\}$, and let 
$\mathcal S_n=\mathbb N_0^\infty/n\cap\mathcal S$. It follows that π n(t) ∈ 𝒮n for all t ∈ [0, T]. Let 
$\Sigma=\{\ell=(\ell_i)_{i=1}^L\in\mathbb N^L_0\mid\ell_1\leq\ell_2\leq\cdots\leq\ell_L\}$ and, for ℓ ∈ Σ, define 
$\rho_i(\ell)\,{:{=}}\, \sum {_{j = 1}^L} {\bf 1}_{\{\ell_j=i\}}, i\in\mathbb N_0$. Roughly speaking, Σ will represent the set of possible states for L selected queues arranged by nondecreasing queue length. Note that each file will be stored at L servers and that at any given time t the queue lengths of these L servers (up to a reordering) will correspond to an element in Σ. We will refer to the elements of Σ as ‘queue length configurations’. Given a configuration ℓ ∈ Σ, ρ i(ℓ) gives the number of queues of length i (among the L selected). From the above description of the system, it follows that the empirical measure process, πn(t), is a continuous-time Markov chain with state space 𝒮n and generator
\begin{align}\label{eqn:generator}
\mathcal L^n f(r)
= \frac{n\lambda}{\binom{n}{L}}\sum_{\ell\in\Sigma}\bigg(\,\prod_{i=0}^\infty\binom{nr_i}{\rho_i(\ell)}\bigg)\Bigl[f\Big(r+\frac{1}{n}\Delta_\ell\Big)-f(r)\Bigr]
+ k\sum_{i=1}^\infty nr_i\Bigl[f\Bigl(r+\frac{1}{n}(e_{i-1}-e_i)\Bigr)-f(r)\Bigr]
\end{align}
for f: 𝒮n → ℝ, where
\begin{align}\label{eqn:deldef}
\Delta_\ell \,{:{=}}\, \sum_{i=1}^k e_{\ell_i+1}-\sum_{i=1}^k e_{\ell_i}
\end{align}
and, for y ∈ ℕ0, e y ∈ ℓ 2, is a vector with 1 at the yth coordinate and 0 elsewhere. Here we use the standard conventions that 
$0^0 = {\binom{0}{0} = 0\text{!}=1}$, and 
$\smash{\binom{a}{b}=0}$ when a < b. The above generator can be understood as follows. A typical term in the second expression corresponds to a jump as a result of a server, with exactly i jobs queued, completing a job. The term in the square brackets gives the change in value of f as a result of such a jump and the prefactor knr i corresponds to the fact that servers process jobs at rate k and there are in all nr i queues (prior to the jump) with exactly i jobs. The first expression in (2) corresponds to a jump resulting from an arrival of a job to the system. Typically, such an arrival makes a request for L servers with queue length configuration ℓ 1 ≤ ℓ 2 ≤ · · · ≤ ℓ L and results in the jump Δℓ/n. The sum in (3) only goes up to k (instead of L) since only the smallest k queues are affected by such a jump. Since, prior to the jump, there are nr i queues with exactly i jobs, the overall rate associated with the configuration ℓ = {ℓ 1 ≤ ℓ 2 ≤ · · · ≤ ℓ L} ∈ Σ equals
\begin{align*}\frac{n\lambda}{\binom{n}{L}}\bigg(\,\prod_{i=0}^\infty\binom{nr_i}{\rho_i(\ell)}\bigg).\end{align*}
In our setting the first entry in an element of ℓ 2 will typically correspond to the number of empty queues, and, thus, we refer to it as the ‘0th’ coordinate and any r ∈ ℓ 2 will correspond to a vector of the form (r 0, r 1, …). For notational convenience, for r ∈ ℓ 2, we set r −1 := 0.
The main results in this work provide scaling limits for πn. We first present the law of large numbers which describes the nominal state of the system for large n. Define, for r ∈ ℓ 1,
\begin{align*}
\bar{\zeta}^\delta(j,r) \,{:{=}}\, \bar{\zeta}(j-1,r)-\bar{\zeta}(j,r),
\end{align*}
where, adopting the convention that 
$\sum {_{i = b}^a} {x_i} = 0$ for a < b,
\begin{align}\label{eqn:zetabardef}
L\text{!}\bar{\zeta}(j,r)\,:\!=\! \sum_{i_1=0}^{k-1}\sum_{i_2=1}^{L-i_1}\frac{L\text{!}}{i_1\text{!}\,i_2\text{!}\,(L-i_1-i_2)\text{!}}\Biggl(\,\sum_{m=0}^{j-1}
r_m\!\Biggr)^{\!\!i_1}\!(r_j)^{i_2}\Biggl(\,\sum_{m=j+1}^{\infty}r_{m}\!\Biggr)^{\!\!L-i_1-i_2}[i_2\wedge(k-i_1)].
\end{align}
This allows us to define, for r ∈ ℓ 1,
\begin{align}\label{eqn:codingF}
F(r) \,{:{=}}\, \lambda L\text{!}\,\sum {_{j = 0}^\infty } \bar{\zeta}^\delta(j,r)e_j +
k\sum {_{j = 0}^\infty }[r_{j+1}-r_j]e_j+ kr_0e_0.
\end{align}
For j ≥ 0, the quantities k[r j+1 − r j] in (5) roughly represent the rate at which the jth coordinate of the state changes (in the limit) as a result of job completions. The final term in the display cancels the extra term in the second summation for the boundary case j = 0. The quantity λL! (ζ̄(j – 1, r) – ζ̄(j, r)) represents a similar quantity as a result of job arrivals. The various terms in (4) can be interpreted as follows. An arrival to a queue with j jobs implies that a queue length configuration vector ℓ = {ℓ 1 ≤ ℓ 2 ≤ · · · ≤ ℓ L} was selected which has the property that at least one of the k smallest i equals j, or, equivalently, exactly i 1 (i 1 = 0, 1, …, k − 1) of the smallest L selected are less than j, i 2 (i 2 = 1, …, L − i 1) of these are equal to j, and L − i 1 − i 2 are greater than j. The three ratios in (4) are contributions from these three types of queues. The term [i 2 ∧ (k − i 1)] arises from the fact that only the smallest k of the L queues are affected.
Also observe that, for some c ζ ∈ (0, ∞),
\begin{align}\label{eqn:czetabound}
\bar \zeta \left({j,r} \right) \le {c_\zeta }{r_j}\,\,\,\,{\rm{for}}\,\;{\rm{all}}\,j \in {{\mathbb{N}}_0}\,\,\,{\rm{and}}\,\,\,r = \left({r{_j}} \right)_{j = 0}^\infty \in {\mathcal{S}}
\end{align}
Thus, the infinite sum in (5) is well defined since 
$\sum {_{j = 0}^\infty {r_j}} = 1$, and, consequently, F is a well-defined map from 𝒮 to ℓ 1. A similar estimate shows that F is a well-defined map from ℓ 1 to ℓ 1 and 
$\sum {_{j = 0}^\infty {F_j}} (r) = 0$ for all r ∈ ℓ 1.
Consider the system of ODEs
\begin{align}\label{eqn:ODE}
\dot{\pi}(t) = F(\pi(t)),\qquad \pi(0) = \pi_0,
\end{align}
where F is defined in (5) and π 0 ∈ 𝒮. The solution of the equation is a continuous map π: [0, T] → 𝒮 such that
\begin{align}\label{eq:eq324}
\pi(t) = \pi_0 + \int_0^t F(\pi(s)) {\rm d} s, \qquad t \in [0,T],
\end{align}
where the integral on the right-hand side is the classical Bochner integral which is well defined since, from (5) and (6),
\begin{align} \label{eq:eq305}
\sup_{0\le s \le T}\|F(\pi(s))\|_1 \le \sup_{r \in \mathcal S} \|F(r)\|_{1} < \infty .
\end{align}
Equation (7) will characterize the law of large number limit of πn.
The following result on the well posedness of (7) will be shown in Section 4.2.
Proposition 1
Let π 0 ∈ 𝒮. Then there exists a π ∈ ℂ([0, T]: 𝒮) that solves (7). Furthermore, if π and π̃ are two elements of ℂ([0, T]: 𝒮) solving (7) with π(0) = π̃(0) = π 0, then π = π̃.
The next theorem gives a law of large numbers for the sequence {πn}n∈ℕ. Recall that we take πn(0) to be nonrandom.
Theorem 1
Suppose that πn(0) → π 0 in 𝒮 as n → ∞. Then πn → π in probability in 𝔻([0, T]: 𝒮), where π is the unique solution of (7) in ℂ([0, T]: 𝒮).
Remark 1
The occupancy process πn satisfies a natural monotonicity property. Let 
$\gamma _i^n(t){\mkern 1mu} : = \sum {_{j = i}^\infty \pi _j^n} (t)$, starting from some initial state γ n ,0 and let γ̃n be the corresponding process starting from an initial state γ̃n ,0. Suppose that γ̃n ,0i ≥ γ n ,0i for every i. Then γ n(t) is componentwise stochastically dominated by γ̃n(t) for every t ≥ 0. A deterministic analogue of this property will hold for the limiting trajectories π. Such a monotonicity can be used to deduce various qualitative properties of πn and the limit path π. Indeed, Friedlander [Reference Friedlander15] used this monotonicity behavior crucially in order to study the long-time behaviors of πn and π.
Our second main result studies the fluctuations of πn from its law of large number limit. Consider
\begin{align}\label{eqn:Xndef}
X^n(t) = \sqrt{n}[\pi^n(t)-\pi(t)], \qquad t \in [0,T],
\end{align}
where πn is the state process introduced in (1) and π is the unique solution of (7) in ℂ([0, T]: 𝒮).
We will show that, under conditions, Xn converges in distribution in 𝔻([0, T]: ℓ 2) to a stochastic process that can be characterized as the solution of an SDE of the form
\begin{align}\label{eqn:limitSDE}
{\rm d} X(t) = G(X(t),\pi(t)){\rm d} t+a(t){\rm d} W(t), \qquad X(0)=x_0.
\end{align}
The equation is again interpreted in the integrated form
\begin{align}\label{eqn:limitSDEb}
X(t) = x_0 + \int_0^t G(X(s),\pi(s)){\rm d} s + \int_0^t a(s){\rm d} W(s), \qquad t \in [0,T].
\end{align}
In the above equations, a is a measurable map from [0, T] to the space of Hilbert–Schmidt operators from ℓ 2 to ℓ 2 such that 
$\mathop \smallint \nolimits_0^T a(t)_{{\rm{HS}}}^2{\rm{ d}}t < \infty$, where ‖ · ‖HS denotes the Hilbert–Schmidt norm (see Appendix A.3), and W is an ℓ 2-cylindrical Brownian motion. Precise definitions are given in Appendix A.4, but roughly speaking, W can be identified with an independent and identically distributed sequence {β i}i∈ℕ0 of standard real Brownian motions over [0, T] and the stochastic integral 
$\int_0^t a (s){\rm{ d}}W(s)$ represents an ℓ 2-valued Gaussian martingale M(t) given as
\begin{align}\label{eq:eq410}
M_i(t) = \sum_{j=0}^{\infty} \int_0^t A_{ij}(s) {\rm d} \beta_j(s), \qquad t \in [0,T],\, i \in \mathbb N_0,
\end{align}
where A ij(s) = 〈e i, a(s)e j〉2, s ∈ [0, T], i, j ∈ ℕ0. We refer the reader to Chapter 4 of [Reference Da Prato and Zabczyk10] for construction and properties of the stochastic integral in (12). The Hilbert–Schmidt norm and integrability property of a ensure that the infinite sum in (13) converges. The operator a(t) is determined from the system parameters and the law of large number limit π in Theorem 1 as the symmetric square root of the following nonnegative trace class operator:
\begin{align}\label{eqn:squarematrix}
\Phi(t) \,{:{=}}\, \lambda
L\text{!}\,\sum_{\ell\in\Sigma}\Delta_\ell\Delta_\ell^T\prod_{i=0}^\infty\frac{\pi_i(t)^{\rho_i(\ell)}}{\rho_i(\ell)\text{!}}+k\sum_{i=1}^\infty(e_{i-1}-e_i)(e_{i-1}-e_i)^T\pi_i(t).
\end{align}
The trace class property of Φ(t) and the integrability of the squared Hilbert–Schmidt norm of a(t) are shown in Lemma 7. Define the space ℓ̃ 2 ⊂ ℓ 2 as
\begin{align}\label{eqn:tibfelldef}
\tilde \ell_2\,{:{=}}\,\Bigg\{x\in\ell_2\colon\sum {_{j = 0}^\infty } j^2x_j^2<\infty,\,\sum {_{j = 0}^\infty } x_j = 0\Bigg\}.
\end{align}
In (11) G is a map from ℓ̃ 2 × 𝒮 to ℓ 2 defined as
\begin{align}\label{eqn:opDFdef}
G_i(x,r)\,{:{=}}\, \frac{\partial}{\partial u}F_i(r+ux)\Big|_{u=0},\qquad i\in\mathbb N_0,\,\; u\in \mathbb R.
\end{align}
One of the difficulties in the analysis is that G as a map from ℓ 2 × 𝒮 to ℓ 2 is not well behaved and we need to restrict attention to the smaller space ℓ̃ 2 × 𝒮 in order to get unique olvability of (11). Note that, under the condition that 
$\sum {_{j = 0}^\infty } {j^2}x_j^2 < \infty $, the series 
$\sum {_{j = 0}^\infty } |{x_j}| < \infty $, and thus the series 
$\sum {_{j = 0}^\infty } {x_j}$, is convergent. Additionally, the right-hand side of (16) is well defined for every x ∈ ℓ̃ 2 and r ∈ 𝒮, since, for each j ∈ ℕ0 and r ∈ ℓ 1 with 
$\sum {_{i = 0}^\infty } {r_i} = 1,r \mapsto {F_j}(r)$ is a polynomial in (r 0, r 1, …, r j+1) given as
\begin{align*}
F_j(r)= \lambda L\text{!}\,[\bar{\zeta}(j-1,r)-\bar{\zeta}(j,r)]+k(r_{j+1}-r_j),
\end{align*}
where
\begin{align*}
\bar{\zeta}(j,r) \,=\,\sum_{i_1=0}^{k-1}\frac{(\sum\nolimits_{m=0}^{j-1}r_m)^{i_1}}{i_1\text{!}}\sum_{i_2=1}^{L-i_1}[i_2\wedge(k-i_1)]\frac{(r_j)^{i_2}}{i_2\text{!}}\frac{(1-\sum\nolimits_{m=0}^{j}r_{m})^{L-i_1-i_2}}{(L-i_1-i_2)\text{!}}.
\end{align*}
Also, from (4) and (5), it is easily checked that there is a c ∈ (0, ∞) such that, for all x ∈ ℓ̃ 2 and r ∈ 𝒮,
\begin{align*}
|G_i(x,r)|\leq c\bigg[|x_{i-1}|+|x_i|+|x_{i+1}|+(r_{i-1}+r_i)\sum_{m=0}^\infty|x_{m}|\bigg].
\end{align*}
This in particular implies that G(x, r) : = (G i(x, r))i∈ℕ0 ∈ ℓ 1 ⊂ ℓ 2 for all (x, r) ∈ ℓ̃ 2 × 𝒮.
The following result shows the well posedness of (12). The definition of an ℓ 2-cylindrical Brownian motion is given in Section A.4.
Proposition 2
There exists a filtered probability space (Ω, 
$\mathcal F$, ℙ, {
$\mathcal F_t$}) on which is given an ℓ 2-cylindrical Brownian motion W and a continuous 
$\{\mathcal F_t\}$-adapted process (X(t))0≤t≤T with sample paths in ℂ([0, T]: ℓ 2) that satisfies the integral equation (12) and is such that X(t) ∈ ℓ̃ 2 ⊂ ℓ 2for all t ∈ [0, T] almost surely. Furthermore, if {X̃ t}0≤t≤T is another such process then X̃ t = X t for all t ∈ [0, T] almost surely.
The above result establishes weak existence and pathwise uniqueness of (12). By a standard argument (cf. [Reference Ikeda and Watanabe18, Section IV.1]), it follows that (12) has a unique weak solution. We can now present our main result on the fluctuations of πn. Recall that 
$X^n(0)=\sqrt{n}(\pi^n(0)-\pi_0)$ is deterministic.
Theorem 2
Suppose that 
$\mathop {\sup }\limits_{n \in \mathbb N } \sum {_{j = 0}^\infty {j^2}} \pi _j^n(0)$ < ∞ and πn(0) → π 0 in 𝒮 as n → ∞. Let π be the unique solution of (7) and, with Xn defined as in (10), Xn(0) → x 0 in ℓ 2. In addition, suppose that
\begin{align}\label{eqn:XinitBound}
{\sup _{n\in\mathbb N}}\sum {_{j = 0}^\infty } j^2(X^n_j(0))^2<\infty.
\end{align}
Then Xn ⇒ X in 𝔻([0, T]: ℓ 2), where X is the unique weak solution to (11) given by Proposition 2.
Proposition 2 and Theorem 2 will be proved in Section 5. In Section 6 we will describe how Theorems 1 and 2 can be used for numerical computation of various performance measures using simulation of diffusion processes.
2.1. Supermarket model
Consider a system of n servers, each with its own queue. Jobs arrive in the system according to a Poisson process with rate nλ. When a job enters the system, d servers are chosen uniformly at random and the job is routed to the shortest of the d selected queues. All servers process jobs according to the FIFO discipline. Service times are mutually independent and exponentially distributed with mean 1. This model has been well studied and is known as power-of-d routeing or the ‘supermarket model’ (see [Reference Graham16], [Reference Mitzenmacher31], and [Reference Vvedenskaya, Dobrushin and Karpelevich37]). The model is a special case of the system considered in the current work, corresponding to L = d and k = 1. Theorems 1 and 2 then provide, as corollaries, the following law of large numbers and central limit theorem for the power-of-d routeing scheme.
Define by 
$\pi^n_d$ the empirical measure process of queue lengths in the power-of-d system. For r ∈ ℓ 1, define
\begin{align*}
F_d(r) \,{:{=}}\, \lambda\bigg[\sum_{i=1}^d\binom{d}{i}r_{j-1}^i\bigg(\sum_{m=j}^\infty r_m\bigg)^{d-i}-\sum_{i=1}^d\binom{d}{i}r_j^i\bigg(\sum_{m=j+1}^\infty r_m\bigg)^{d-i}\bigg]e_j + \sum {_{j = 0}^\infty }[r_{j+1}-r_j]e_j.
\end{align*}
The following is a direct corollary of Theorem 1.
Corollary 1
Suppose that 
$\pi_d^n(0) \to \pi_d(0)$ in 𝒮 as n → ∞. Then 
$\pi_d^n\to\pi_d$ in probability in 𝔻([0, T]: 𝒮), where πd is the unique solution in ℂ([0, T]: 𝒮) to the ODE
\begin{align*}
\dot{\pi}_d(t) = F_d(\pi_d(t)),\qquad \pi_d(0) = \pi_0.
\end{align*}
Remark 2
This result has been established in [Reference Graham16] (see Theorem 3.4 therein). In particular, it is easy to verify that 
${\mathcal V_m}(t){\mkern 1mu} : = {\mkern 1mu} \sum {_{j = m}^\infty {{({\pi _d}(t))}_j}}$ is the same function as in Equation (3.9) of [Reference Graham16] (see also [Reference Vvedenskaya, Dobrushin and Karpelevich37]).
Our second corollary studies the fluctuations of 
$\pi^n_d$ from its law of large number limit. Consider
\begin{align*}
X^n_d(t) = \sqrt{n}[\pi^n_d(t)-\pi_d(t)], \qquad t \in [0,T].
\end{align*}
Analogous to a(t) introduced in (11), let a d(t) be the symmetric square root of the nonnegative operator
\begin{align}\label{eqn:squarematrixd}
\Phi_d(t) \,{:{=}}\, \lambda
\sum {_{j = 0}^\infty }(e_{j+1}-e_{j})(e_{j+1}-e_{j})^T\bigg(\sum_{i=1}^d\binom{d}{i}[(\pi_d)_j(t)]^i\bigg(\sum_{m=j+1}^\infty(\pi_d)_m(t)\bigg)^{d-i}\bigg)
+\sum_{j=1}^\infty (e_{j-1}-e_j)(e_{j-1}-e_j)^T(\pi_d)_j(t).
\end{align}
Analogous to G in (16), let G d be a map from ℓ̃ 2 × 𝒮 to ℓ 2, where ℓ̃ 2 is as in (15), defined as
\begin{align}\label{eqn:opDFdefd}
(G_d)_i(x,r)\,{:{=}}\, \frac{\partial}{\partial u}(F_d)_i(r+ux)\Big|_{u=0},\qquad i\in\mathbb N_0, u\in \mathbb R.
\end{align}
In the special case that d = 2, this function simply reduces to
\begin{align*}
{\left({{G_2}} \right)_i}\left({x,r} \right) = 2\lambda \sum\limits_{m = i}^\infty {\left[ {{x_{i - 1}}{r_m} + {r_{i - 1}}{x_m} - {x_i}{r_{m + 1}} - {r_i}{x_{m + 1}}} \right]} + \left({{x_{i + 1}} - {x_i}} \right).
\end{align*}
The following result is immediate from Theorem 2.
Corollary 2
Suppose that 
$\mathop {\sup }\limits_{n \in \mathbb N} \sum {_{j = 0}^\infty {j^2}} {(\pi _d^n)_j}(0) < \infty$ and 
$\pi^n_d(0) \to \pi_0$ in 𝒮 as n → ∞. Also, suppose that 
$X_d^n(0)=\sqrt{n}[\pi_d^n(0)-\pi_0]\to x_0$ in probability in ℓ 2 and that
\begin{align*}
{\sup _{n\in\mathbb N}}\sum {_{j = 0}^\infty } j^2((X_d^n)_j(0))^2<\infty.
\end{align*}
Then 
$X_d^n\Rightarrow X_d$ in 𝔻([0, T]: ℓ 2), where Xd is the unique weak solution to (11) with values in 
${\tilde\ell _2}$, with G replaced by Gd defined by (19) and a(t) replaced by a d(t) which is given as the symmetric square root of the operator Φd(t) in (18).
3. Semimartingale representation
In this section we write the state processes using compensated time-changed Poisson processes to give a semimartingale representation for the system. Let {N ℓ, ℓ ∈ Σ} and {D i, i ∈ ℕ0} be collections of mutually independent unit-rate Poisson processes. The process N ℓ will be used to represent the stream of jobs requesting files which are stored at servers with queue length configuration (immediately before the time of arrival of the request) ℓ= (ℓ 1, …, ℓ L). Similarly, D i will represent the stream of jobs completed by servers whose queue length (immediately before the time of completion) is equal to i. From the form of the generator in (2) we see that the state process πn can be expressed as
\begin{align*}
\pi^n(t) = \pi^n(0) + \frac{1}{n}\sum_{\ell\in\Sigma}\Delta_\ell
N_\ell\biggl(\int_0^t\frac{n\lambda
}{\binom{n}{L}}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\biggr)
+\frac{1}{n}\sum_{i=1}^\infty
(e_{i-1}-e_i)D_i\Big(k\int_0^tn\pi^n_i(s){\rm d} s\Big).
\end{align*}
By adding and subtracting the compensators of the Poisson processes we can write the state process as a semimartingale. Namely,
\begin{align}\label{eqn:semimartrep}
\pi^n(t)=\pi^n(0)+A^n(t)+M^n(t),
\end{align}
where
\begin{align}\label{eqn:Adef}
A^n(t)
\,{:{=}}\,\sum_{\ell\in\Sigma}\Delta_\ell\int_0^t\frac{\lambda}{\binom{n}{L}}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}
{\rho_i(\ell)}{\rm d} s
+k\sum_{i=1}^\infty (e_{i-1}-e_i)\int_0^t\pi^n_i(s){\rm d} s
\end{align}
and
\begin{align}\label{eqn:MartRep}
M^n(t)\,{:{=}}\, \sum_{\ell\in\Sigma}\frac{1}{n}\Delta_\ell
N_\ell\biggl(\frac{n\lambda}{\binom{n}{L}}\int_0^t\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\biggr)-\sum_{\ell\in\Sigma}\Delta_\ell
\frac{\lambda}{\binom{n}{L}}\int_0^t\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s
+ \sum_{i=1}^\infty \frac{1}{n}(e_{i-1}-e_i)
D_i\bigg(k\int_0^tn\pi^n_i(s){\rm d} s\bigg)- k\sum_{i=1}^\infty (e_{i-1}-e_i)
\int_0^t\pi^n_i(s){\rm d} s.
\end{align}
It will follow from (45) and (54) that both Mn(t) and An(t) take values in ℓ 2. Specifically, (45) and (54) imply that, for some c ζ ∈ (0, ∞),
\begin{align*}
A_j^n(t) \leq \int_0^t\biggl(\frac{\lambda}{\binom{n}{L}}c_\zeta
n^L[\pi^n_{j-1}(s)+\pi_j^n(s)]+k[\pi_{j+1}^n(s)+\pi^n_j(s)]\biggr){\rm d} s
\end{align*}
for all t ∈ [0, T], n ∈ ℕ, and j ∈ ℕ0. Thus, there exists a κ ∈ (0, ∞) such that
\begin{align*}
\sum\limits_{j = 0}^\infty A_j^n{(t)^2} \le \kappa \sum\limits_{j = 0}^\infty{\int_0^t} [\pi _{j - 1}^n{(s)^2} + \pi _{j + 1}^n{(s)^2} + \pi _j^n{(s)^2}]{\kern 1pt} ds \le 3\kappa t
\end{align*}
for all t ∈ [0, T]. A similar argument shows that An(t) in fact takes values in ℓ 1. In the next section we establish tightness of {Mn} as a sequence of ℓ 2-valued processes.
Similarly, using (20) and (8) for π(t), we can express Xn as a semimartingale through the equation
\begin{align}\label{eqn:semiMart}
X^n(t) = X^n(0) + A^{n,c}(t) + M^{n,c}(t),
\end{align}
where
\begin{align}\label{eqn:Abardef}
A^{n,c}(t) = \sqrt{n}\Big[A^n(t)-\int_0^tF(\pi(s)){\rm d} s\Big]
\end{align}
and 
$M^{n,c}(t) = \sqrt{n}M^n(t)$. We note that there is a natural filtration 
$\{\mathcal F^n_t\}_{0\leq t\leq T}$ on the probability space where the processes N ℓ, D i, and πn are defined such that An, Mn, πn, Xn, Mn ,c , and An ,c are RCLL processes adapted to the filtration, and Mn and Mn ,c are 
$\{\mathcal F^n_t\}$-local martingales.
4. Law of large numbers
In this section we present the proof of Theorem 1. First, in Section 4.1 we use the semi-martingale representation from Section 3 to prove a key tightness property (see Proposition 3). Then, in Section 4.2 we prove the unique solvability of (7) and complete the proof of Theorem 1 by proving convergence of πn to the unique solution of (7) in ℂ([0, T]: 𝒮).
4.1. Tightness
In this section we prove tightness of {(πn, Mn)}n∈ℕ. We first recall the notion of ℂ-tightness. From [Reference Ethier and Kurtz14, Theorem 3.10.2], it can be seen that the definition presented below is equivalent to the more standard definition of ℂ-tightness (see [Reference Jacod and Shiryaev19, Definition VI.3.25]).
Definition 1
Let 
$(\mathcal Z, d_{\mathcal Z})$ be a Polish space. For 
$z\in\mathbb D([0,T]\colon\mathcal Z)$, let
\begin{align*}
j_T(z)\,{:{=}}\,\sup_{0\le t\leq T}d_{\mathcal Z}(z(t),z(t-)).
\end{align*}
A sequence {Z n}n∈ℕ of 
$\mathbb D([0,T]\colon\mathcal Z)$-valued random variables is said to be ℂ-tight if it is tight in 
$\mathbb D([0,T]\colon\mathcal Z)$ and j T(Z n) ⇒ 0.
If Z n and Z are 
$\mathbb D([0,T]\colon\mathcal Z)$-valued random variables and Z n ⇒ Z then 
$\mathbb P(Z\in\mathbb C([0,T]\colon\mathcal Z))=1$ if and only if {Z n}n∈ℕ is ℂ-tight [Reference Ethier and Kurtz14]. The following proposition proves the ℂ-tightness of {πn}n∈ℕ and convergence of Mn to the zero process.
Proposition 3
Suppose that πn(0) → π 0 in 𝒮 as n → ∞. Then {(πn, Mn)}n∈ℕ is a ℂ-tight sequence of 𝔻([0, T]: 𝒮 × ℓ 2)-valued random variables. Furthermore, Mn ⇒ 0 in 𝔻([0, T]: ℓ 2).
Proof. We first prove the second statement by arguing that 
$\mathbb E\sup_{0\leq s\leq T}\|M^n(s)\|_2^2\to 0$ as n → ∞. For this, from Doob’s inequality, it suffices to show that 𝔼|〈Mn〉(T)| → 0 as n → ∞, where
\begin{align*}
\left\langle M^n\right\rangle(s) \,{:{=}}\,\sum {_{j = 0}^\infty }\left\langle M^n_j\right\rangle(s),\qquad s\in[0,T].
\end{align*}
Note that 
${\left\langle {{e_i},{e_k}e_j^T{e_m}} \right\rangle _2} = {\delta _{ik}}{\delta _{jm}}$ for all i, j, k, l ∈ ℕ0, where δ xy is 1 if x = y and 0 otherwise. It then follows that
\begin{align}\label{eqn:processed1}
\sum_{i=1}^\infty \left\langle e_j,(e_{i-1}-e_i)(e_{i-1}-e_i)^Te_j\right\rangle_2 \pi^n_i(s)
=\sum_{i=1}^\infty \left\langle e_j,(e_{i-1}e_{i-1}^T+e_ie_i^T)e_j\right\rangle_2 \pi^n_i(s)\mathbb N
=\pi^n_{j+1}(s)+\pi^n_j(s).
\end{align}
Since {N ℓ, D i ℓ ∈ Σ, i ∈ ℕ0} are mutually independent Poisson processes, we now have, from (22),
\begin{align}\label{eqn:Mjs}
\left\langle {M_j^n} \right\rangle \left(t \right) = {\lambda \over {n\left({_L^n} \right)}}\int_0^t {Z\left({j,n{\pi ^n}\left(s \right)} \right)} {\rm{d}}s + {k \over n}\int_0^t {\left[ {\pi _{j + 1}^n\left(s \right) + \pi _j^n\left(s \right)} \right]} {\rm{d}}s
\end{align}
where
\begin{align}\label{eqn:incoming1}
Z(j,n\pi^n(s))=\sum_{\ell\in\Sigma}\left\langle e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2 \prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}.
\end{align}
The ℓth term in the sum on the right-hand side of (27) is the contribution from jobs that request servers with queue length configuration ℓ. A fixed ℓ ∈ Σ will make a nonzero contribution to 
$\left\langle e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2$ if j or j − 1 is one of the k-smallest coordinates in ℓ. Thus, for a fixed ℓ ∈ Σ, the ℓth term in (27) is nonzero only if j or j − 1 is a member of the set (ℓ 1, …, ℓ k). The contribution from all such in the sum (27) can be counted as follows. Suppose that 0 ≤ i 1 ≤ k − 1 servers are selected among those with queue length less than j − 1. This corresponds to 
$\left({n\sum {_{_{{i_1}}^{m = 0}}^{j - 2}} \pi _m^n\left(s \right)} \right)$ different choices of servers. In addition, suppose that i 2 ≤ L − i 1 and i 3 ≤ L − i 1 − i 2 servers are selected among those with queue length equal to j − 1 and j, respectively. This corresponds to 
$\Big(\begin{smallmatrix}{n\pi^n_{j-1}(s)}\\{i_2}\end{smallmatrix}\Big)$ and 
$\Big(\begin{smallmatrix}n\pi^n_{j}(s)\\{i_3}\end{smallmatrix}\Big)$ choices, respectively. It follows that L − i 1 − i 2 − i 3 servers must be selected which have queue length larger than j, which corresponds to 
$\left({n\sum {_{_{L - {i_1} - {i_2} - {i_3}}^{m = j + 1}}^\infty \pi _m^n\left(s \right)} } \right)$ possible choices. Since jobs are only routed to the k-shortest servers, we have, with
\begin{align}\label{eq:eq340}
\Sigma_j(i_1,i_2,i_3)\,{:{=}}\,\bigg\{\ell\in\Sigma|\sum_{i=1}^{j-2}\rho_i(\ell)=i_1,\,
\rho_{j-1}(\ell)=i_2,\,\rho_{j}(\ell)=i_3\bigg\},
\end{align}
for ℓ ∈ Σj(i 1, i 2, i 3),
\begin{align}\label{eqn:innerEval}
\left\langle e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2=[i_2\wedge (k-i_1)-i_3\wedge(k-i_1-i_2)_+]^2,
\end{align}
and, thus, for x ∈ n𝒮n,
\begin{align}\label{eqn:Zdef}
Z\left({j,x} \right) = \sum\limits_{{i_{1 = 0}}}^{k - 1} {\left({\sum {_{_{{i_1}}^{m = 0}}^{j - 2}{x_m}} } \right)} \sum\limits_{{i_2} = 0}^{L - {i_1}} {\left({_{{i_2}}^{{x_{j - 1}}}} \right)} \times \sum\limits_{{i_{3 = 0}}}^{L - {i_1} - {i_2}} {{{\left[ {{i_2} \wedge {{\left({k - {i_1} - {i_2}} \right)}_ + }} \right]}^2}\left({_{{i_3}}^{{x_j}}} \right)} \left({\sum\limits_{L - {i_1} - {i_2} - {i_3}} {_{m = j + 1{x_m}}^\infty } } \right),
\end{align}
recalling that we adopt the convention that, for a < b, 
$\sum {_{i = b}^a{x_i}} = 0$.
Note that, for nonnegative integers a, b, a ≥ b,
\begin{align}\label{eqn:binominequality}
\binom{a}{b}\leq \frac{a^b}{b\text{!}}.
\end{align}
his fact, combined with (30) and recalling the fact that πn(s) ∈ 𝒮 for s ∈ [0, T], gives the following bound on Z(j, nπn(s)):
\begin{align}\label{eqn:Zineq}
Z(j,n\pi^n(s))
\leq \sum_{i_1=0}^{k-1}\frac{(n\sum_{m=0}^{j-2}\pi^n_m(s))^{i_1}}{i_1\text{!}}\sum_{i_2=0}^{L-i_1}\frac{(n\pi^n_{j-1}(s))^{i_2}}{i_2\text{!}}
\mathbb N
\\
\times\sum_{i_3=0}^{L-i_1-i_2}k^2{\bf 1}_{\{i_2\vee i_3 > 0\}}\frac{(n\pi^n_j(s))^{i_3}}{i_3\text{!}}\frac{(n\sum_{m=j+1}^{\infty}\pi^n_{m}(s))^{L-i_1-i_2-i_3}}{(L-i_1-i_2-i_3)\text{!}}
\\
\leq n^L\sum_{i_1=0}^{k-1}\sum_{i_2=0}^{L-i_1}\sum_{i_3=0}^{L-i_1-i_2}k^2{\bf 1}_{\{i_2\vee i_3 > 0\}}(\pi^n_{j-1}(s))^{i_2}(\pi^n_j(s))^{i_3}
\leq c_Zn^L(\pi^n_{j-1}(s)+\pi^n_j(s))
\end{align}
for some c Z ∈ (0, ∞). Using (32) in (26) gives
\begin{align}\label{eqn:quadMRate}
\mathbb E|\left\langle M^n\right\rangle(t)|
\leq\mathbb E| \frac{2\lambda(n-L)\text{!}\,L\text{!}\,c_Zn^L}{n\times n\text{!}}\int_0^t\sum {_{j = 0}^\infty }\pi^n_j(s){\rm d} s|
+\mathbb E|\frac{2k}{n}\int_0^t\sum {_{j = 0}^\infty }\pi^n_j(s){\rm d} s|
\mathbb N
\leq \left| \frac{2\lambda(n-L)\text{!}\,L\text{!}\,c_Zn^L}{n\times n\text{!}}t\right|+\left|\frac{2k}{n}t\right|.
\end{align}
Thus, 𝔼|〈Mn〉T| → 0 and, consequently, 
$\mathbb E\mathop {\sup }\nolimits_{0 \le s \le T} \parallel {M^n}(s)\parallel _2^2 \to 0$ as n → ∞. It follows that Mn ⇒ 0 in 𝔻([0, T]: ℓ 2), which completes the proof of the second statement.
Tightness of {πn}n∈ℕ in 𝔻([0, T]: 𝒮) follows as in the proof of Theorem 3.4 of [Reference Graham16]. Namely, it suffices to show tightness of 
$\{Q_1^n\}_{n\in\mathbb N}$ in 𝔻([0, T]: 𝔻) (cf. [Reference Sznitman36]). However, this tightness is an immediate consequence of the fact that the jumps of 
$Q_1^n$ can be embedded in a Poisson process with rate λL + k.
Finally, in order to show that {πn}n∈ℕ is ℂ-tight, it suffices to show that
\begin{align*}%\label{eqn:ctight}
j_T(\pi^n) \,{:{=}}\,\sup{_{0 \le t \le T}}d_0(\pi^n(t),\pi^n(t-)) \to0\quad {\rm as}\,\;n\to\infty.
\end{align*}
Note that the arrivals to the nth system occur according to a Poisson process with rate nλ. When a job arrives in the system, the dispatcher assigns it to k different servers, causing the queue length of each of the k chosen servers to increase by one. The n servers in the system process jobs according to Poisson processes with rate 1/k. Any completion of job processing results in the corresponding queue length dropping by 1. These n + 1 Poisson processes are mutually independent from the construction in Section 3, which ensures that the compensated processes,
\begin{align*}
N_\ell\bigg(\int_0^t\frac{n\lambda}{\binom{n}{L}}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg) -\int_0^t\frac{n\lambda}{\binom{n}{L}}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s
\end{align*}
and
\begin{align*}
D_i\bigg(k\int_0^tn\pi^n_i(s){\rm d} s\bigg)- k\int_0^tn\pi^n_i(s){\rm d} s,
\end{align*}
are martingales with respect to a common filtration. Therefore, these Poisson processes do not have simultaneous jumps. Each arrival creates a jump for the vector πn(t) with
\begin{align*}
d_0(\pi^n(t),\pi^n(t-)) \le \frac{2k}{n},
\end{align*}
and each service completion event produces a jump for the vector πn(t) with
\begin{align*}
d_0(\pi^n(t),\pi^n(t-)) \le \frac{2}{n}.
\end{align*}
Thus, almost surely, at any instant the jump size d 0(πn(t), πn(t −)) is at most 2k/n. Therefore,
\begin{align*}
j_T(\pi^n) \leq \frac{2k}{n} \to 0,
\end{align*}
which completes the proof.
4.2. Convergence
In this section we provide the proof of Theorem 1. Since we have already proved tightness of {πn}n∈ℕ in Section 4.1, all that remains is to prove uniqueness of the solutions of (7) in an appropriate class and to characterize the limit of any weakly convergent subsequence as the unique solution to (7). We first present the following Lipschitz property for the map F: 𝒮 → ℓ 1, defined in (5), that will give uniqueness of the solutions to (7). We note that in the proof of Theorem 2 we will need a stronger Lipschitz property of F in the ℓ 2 norm. This Lipschitz property is not immediate on the space 𝒮, but, as shown in Lemma 4, is satisfied on a smaller class 𝒱M.
Lemma 1
The map F is a Lipschitz function from 𝒮 to ℓ 1. Namely, there exists a C 1 ∈ (0, ∞) such that, for any r, r̃ ∈ 𝒮,
\begin{align}\label{eqn:LipschitzDisplay1}
\|F(r)-F(\tilde r)\|_1\leq C_1\|r-\tilde r\|_1.
\end{align}
Proof. Let r, r̃ ∈ 𝒮 and, for i 1 ∈ ℕ0 and j, i 2 ∈ ℕ, define R j,i 1,i 2(r, r̃) as
\begin{align}\label{eqn:Rdef}
R_{j,i_1,i_2}(r,\tilde r)
\,{:{=}}\,\bigg(\sum_{m=0}^{j-1}r_m\bigg)^{i_1}\bigg(\sum_{m=j+1}^{\infty}r_{m}\bigg)^{L-i_1-i_2}r^{i_2}_j
- \bigg(\sum_{m=0}^{j-1}\tilde r_m\bigg)^{i_1}\bigg(\sum_{m=j+1}^{\infty}\tilde
r_{m}\bigg)^{L-i_1-i_2}\tilde r_j^{i_2}.
\end{align}
Note that, for any a, b, c, ã, b̃, c̃ ∈ ℝ+,
\begin{align}\label{eqn:tripleIneq}
abc-\tilde a \tilde b \tilde c= ab (c-\tilde c)+a(b-\tilde b)\tilde c+(a-\tilde a)\tilde b\tilde c.
\end{align}
Combining (35), (36), and the fact that r, r̃ ∈ 𝒮, we have
\begin{align}\label{eqn:Rbound1}
|R_{j,i_1,i_2}(r,\tilde r)| \leq |r_{j}^{i_2}-\tilde r_{j}^{i_2}| +\tilde
r_{j}^{i_2}\Bigg|\bigg(\sum_{m=j+1}^\infty
r_m\bigg)^{L-i_1-i_2}-\bigg(\sum_{m=j+1}^\infty \tilde
r_m\bigg)^{L-i_1-i_2}\Bigg|+\tilde r_{j}^{i_2}\Bigg|\bigg(\sum_{m=0}^{j-1}
r_m\bigg)^{i_1}-\bigg(\sum_{m=0}^{j-1} \tilde r_m\bigg)^{i_1}\Bigg|.
\end{align}
For any a, b ∈ ℝ and i ∈ ℕ, 
$\left({{a^i} - {b^i}} \right) = \left({a - b} \right)\sum {_{j = 1}^i} {a^{i - j}}{b^{j - 1}}$. Thus, if a, b ∈ [0, 1] and i ≤ L, |ai − bi| ≤ |a − b|L. This inequality along with (37) implies that there exist κ 1, 
${\kappa_1,\kappa '_1} > 0$ such that, for all i 1, i 2 ≤ L, i 2 > 0,
\begin{align}\label{eqn:Rineq}
|R_{j,i_1,i_2}(r,\tilde r)| \leq \kappa_1'\bigg(|r_{j}-\tilde r_{j}| +\tilde
r^{i_2}_{j}\sum_{m=j+1}^\infty|r_m-\tilde r_m|+\tilde
r^{i_2}_{j}\sum_{m=0}^{j-1}|r_m-\tilde r_m|\bigg)
\leq \kappa_1(|r_{j}-\tilde r_{j}|+\tilde r_{j}\|r-\tilde r\|_1).
\end{align}
The definition of F (see (5)) and the triangle inequality imply that
\begin{align}\label{eqn:Ftriangle}
\|F(r)-F(\tilde r)\|_1 \leq \lambda L\text{!}\sum {_{j = 0}^\infty }
|\bar{\zeta}^\delta(j,r)-\bar{\zeta}^\delta(j,\tilde r)| + k\sum {_{j = 0}^\infty }|(r-\tilde
r)_{j+1}-(r-\tilde r)_j|.
\end{align}
Noting that
\begin{align*}
\bar{\zeta}^\delta(j,r)-\bar{\zeta}^\delta(j,\tilde r) =
[\bar{\zeta}(j-1,r)-\bar{\zeta}(j-1,\tilde r)]-
[\bar{\zeta}(j,r)-\bar{\zeta}(j,\tilde r)],
\end{align*}
it follows that
\begin{align}\label{eqn:zetatoR}
\sum {_{j = 0}^\infty } |\bar{\zeta}^\delta(j,r)-\bar{\zeta}^\delta(j,\tilde r)| \leq
2\sum {_{j = 0}^\infty } |\bar{\zeta}(j,r)-\bar{\zeta}(j,\tilde r)|\leq
\kappa_2\sum {_{j = 0}^\infty }\sum_{i_1=0}^{k-1}\sum_{i_2=1}^{L-i_1}|R_{j,i_1,i_2}(r,\tilde
r)|,
\end{align}
where the second inequality follows from the the definitions of ζ̅ and R. Combining (40) with (39) and applying (38) yields, for some κ 3 < 0,
\begin{align*}
\|F(r)-F(\tilde r)\|_1 \leq \kappa_2\lambda
L\text{!}\sum {_{j = 0}^\infty }\sum_{i_1=0}^{k-1}\sum_{i_2=1}^{L-i_1}|R_{j,i_1,i_2}(r,\tilde
r)|+ 2k\sum {_{j = 0}^\infty }|r_j-\tilde r_j|
\leq \kappa_3\sum {_{j = 0}^\infty }[|r_{j}-\tilde r_{j}| +\tilde r_{j}\|r-\tilde
r\|_1]+ 2k\|r-\tilde r\|_1,
\end{align*}
and, thus, with C 1 := 2(κ 3 + k), (34) is satisfied for all r, r̃ ∈ 𝒮, which proves the result.
Using the above Lipschitz property of F, we can now complete the proof of Proposition 1.
Proof of Proposition 1. Existence of a π ∈ ℂ([0, T]: 𝒮) that solves (7) will be shown below in the proof of Theorem 1. We now argue uniqueness. Suppose that π and π̃ are two elements of ℂ([0, T]: 𝒮) satisfying (7) with π(0) = π̃(0) = π 0. The Lipschitz property of F proved in Lemma 1 implies that, for all t ∈ [0, T],
\begin{align*}
\|\pi(t)-\tilde\pi(t)\|_1 =
{\left\| {\int_0^t {[F\left({\pi \left(s \right)} \right) - F\left({\tilde \pi \left(s \right)} \right)]{\rm{d}}s} } \right\|_1}
\leq\int_0^t\|F(\pi(s))-F(\tilde\pi(s))\|_1{\rm d} s
\leq C_1\int_0^t\|\pi(s)-\tilde\pi(s)\|_1{\rm d} s.
\end{align*}
The result follows.
We now proceed to the proof of Theorem 1.
Proof of Theorem 1. From Proposition 3, {πn}n∈ℕ is a ℂ-tight sequence of 𝔻([0, T]: 𝒮)-valued random variables.
Note from (20) that, for all j ∈ ℕ0,
\begin{align}\label{eqn:piexpansion}
\pi^n(t) = \pi^n(0)+V^n(t)+M^n(t)+\int_0^tF(\pi^n(s)){\rm d} s,
\end{align}
where
\begin{align*}
V^n(t) \,{:{=}}\, A^n(t) - \int_0^tF(\pi^n(s)){\rm d} s.
\end{align*}
From the definition of An in (21) we see that
\begin{align}\label{eqn:Ajdef}
A^n_j(t) = \int_0^t\bigg(\,\sum_{\ell\in\Sigma}\left\langle\Delta_\ell,
e_j\right\rangle_2\frac{\lambda}{\binom{n}{L}}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}
+k [\pi_{j+1}^n(s)-\pi_{j}^n(s)]\bigg){\rm d} s.
\end{align}
By a similar argument (see the comments given below (44)) used to obtain the representation in (30),
\begin{align}\label{eqn:jexpansion5}
\sum_{\ell\in\Sigma}\left\langle\Delta_\ell,
e_j\right\rangle_2\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)} =
[\zeta(j-1,n\pi^n(s))-\zeta(j,n\pi^n(s))],
\end{align}
where, for x ∈ n𝒮n,
\begin{align}\label{eqn:zetadef}
\zeta \left({j,x} \right): = \sum\limits_{{i_{1 = 0}}}^{k - 1} {\left({\sum {_{_{{i_1}}^{m = 0}}^{j - 1}{x_m}} } \right)\sum\limits_{{i_{2 = 1}}}^{L - {i_1}} {\left[ {{i_2} \wedge \left({k - {i_1}} \right)} \right]} } \left({_{{i_2}}^{{x_j}}} \right)\left({\sum {_{_{L - {i_1} - {i_2}}^{m = j + 1}}^\infty {x_m}} } \right).
\end{align}
One can interpret ζ(j, x) as the rate at which jobs are being routed into queues of length j when the system is in state x. Recall that any incoming job corresponds to the selection of L queues. The term on the right-hand side of (44) then sums over all possible queue length configurations of this selection. In particular, i 1 represents the number of queues with lengths less than j, i 2 corresponds to the queues of length equal to j, and L − i 1 − i 2 are the queues of length greater than j. Since we are routeing jobs to the k-shortest queues, the rate must be multiplied by the factor [i 2 ∧ (k − i 1)] rather than i 2. From our convention that x −1 = 0, we see that ζ(− 1, x) = 0. In addition, recalling the conventions that, for a < b, 
$\sum {_{i = b}^a{x_i}} = 0$ and 
${\binom{0}{0}}=1$, we see that ζ(0, x) is well defined. Combining (42), (43), and (44) gives the following representation for 
$A^n_j$:
\begin{align}\label{eqn:AjzetaExp}
A_j^n\left(t \right) = {\lambda \over {\left({_L^n} \right)}}\int_0^t {[\zeta \left({j - 1,n{\pi ^n}\left(s \right)} \right) - \zeta \left({j,n{\pi ^n}\left(s \right)} \right)]ds} + k\int_0^t {[\pi _{j + 1}^n\left(s \right) - \pi _j^n\left(s \right)]ds}
\end{align}
For each fixed j, i 1 ∈ ℕ0 and i 2 ∈ ℕ with i 1, i 2 ≤ L, we have
\begin{align}\label{eqn:binomToExp}
\Bigg(\begin{array}{@{}c@{}}{n\sum_{m=0}^{j-1}\pi^n_m(s)}{i_1}\end{array}\Bigg)\,[i_2\wedge(k-i_1)]\binom{n\pi^n_j(s)}{i_2}\binom{n\sum_{m=j+1}^{\infty}\pi^n_{m}(s)}{L-i_1-i_2}
\\
=
n^L\frac{(\!\sum_{m=0}^{j-1}\pi^n_m(s))^{i_1}}{i_1\text{!}}[i_2\wedge(k-i_1)]
\frac{(\pi^n_j(s))^{i_2}}{i_2\text{!}}\frac{(\!\sum_{m=j+1}^{\infty}\pi^n_{m}(s))^{L-i_1-i_2}}{(L-i_1-i_2)\text{!}}
\\
+\hat{R}_n(j,i_1,i_2,s),
\end{align}
where
\begin{align*}
\sup_{i_1,i_2\leq L}|\hat{R}_n(j,i_1,i_2,s)|\leq \kappa_1n^{L-1}\pi^n_j(s),
\end{align*}
and, thus, from the definitions of ζ and 
$\bar{\zeta}$ in (44) and (4),
\begin{align}
\label{eqn:zetatobar}
\left|\zeta(j,n\pi^n(s))-
\frac{n\text{!}}{(n-L)\text{!}}\bar{\zeta}(j,\pi^n(s))\right|\leq
\kappa_2n^{L-1}\pi^n_j(s) \quad{\rm for\,\; all}\,\; s\in[0,T].
\end{align}
Furthermore, using the definitions of An in (45) and F in (5), (47) implies that
\begin{align}
\label{eqn:AtoF}
\sup{_{0 \le t \le T}}\|V^n(t)\|_2=\sup{_{0 \le t \le T}}\left\|A^n(t)-
\int_0^tF(\pi^n(s)){\rm d} s\right\|_2\leq \frac{\kappa_3}{n}.
\end{align}
Also, from Proposition 3, Mn ⇒ 0 in 𝔻([0, T]: ℓ 2). Combining these observations with tightness of πn, we have subsequential convergence of (πn, Mn, Vn) to (π, 0, 0), in distribution, in 𝔻([0, T]: 𝒮 × ℓ 2 × ℓ 2) for some ℂ([0, T]: 𝒮)-valued π. By appealing to the Skorokhod representation theorem we can assume that this convergence holds almost surely. Noting that r ↦ F j (r) is a continuous map from 𝒮 to ℝ for each j ∈ ℕ0, we have F j(πn(s)) → F j(π(s)) as ℕ0 and s ∈ [0, T]. Thus, upon sending n → ∞ in (41), (9) and the dominated convergence theorem imply that, almost surely,
\begin{align*}
\pi_j(t) = (\pi_0)_j+\int_0^tF_j(\pi(s)){\rm d} s \quad{\rm for\,\; all}\,\; t\in[0,T], j\in\mathbb N_0.
\end{align*}
This shows that π satisfies (7). The result now follows from the uniqueness property shown in Proposition 1.
5. Diffusion approximation
In this section we prove Theorem 2. Section 5.1 presents some moment estimates on πn which will be used in the proof of Theorem 2. Section 5.2 then proves tightness of the sequence of centered and scaled state processes {Xn}n∈ℕ. Section 5.3 completes the proof of Theorem 2 by proving unique solvability of the SDE (11) (Theorem 2) and characterizing limit points of Xn as this unique solution.
5.1. Moment bounds
The following elementary lemma will be useful in the proof of Lemma 3.
Lemma 2
For all t ≥ 0, k ∈ ℕ, and n ∈ ℕ, 
$\mathop {{{\lim }_{m \to \infty }}}\limits_{} {\mathbb Em^k}\mathop {{{\sup }_{0 \le s \le t}}}\limits_{} \pi _m^n(s) = 0$.
Proof. Fix n ℕ. Note that file requests arrive at rate nλ. Let N be a Poisson process representing the total flow of such file requests. Also, let 
$m^*=\sup\{m\colon\pi_m^n(0)>0\}$ be the length of the largest queue at time 0. Note that, since the system consists of n queues, m ∗ must be finite for any fixed n. Then, for m > m ∗,
\begin{align*}
\mathbb E m^k\sup_{0\leq s\leq t}\pi^n_m(s) = \mathbb E\sup_{0\leq s\leq t}{\bf
1}_{\{N(t)\geq m-m^*\}} m^k\pi^n_m(s)+\mathbb E\sup_{0\leq s\leq t}{\bf
1}_{\{N(t)< m-m^*\}}m^k\pi^n_m(s)
\leq m^k\mathbb P(N(t)\geq m-m^*).
\end{align*}
Thus, from Markov’s inequality, for m > m ∗,
\begin{align*}
\mathbb E m^k\sup_{0\leq s\leq t}\pi^n_m(s)
\leq m^k\rm e^{-(m-m^*)}\rm e^{n\lambda t(\rm e-1)}.
\end{align*}
The result follows.
In the next lemma we will we establish two key moment bounds that will be needed in the tightness proof (see proof of Proposition 4).
Lemma 3
Suppose that 
${\sup _{n \in\mathbb N}}\sum {_{j = 0}^\infty {j^2}} \pi _j^n(0) = :{\mkern 1mu} {c_{\pi (0)}} < \infty $. Then
\begin{align}\label{eqn:probmzrtight}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\Bigg(\sum {_{j = 0}^\infty } j\pi_j^n(t)\Bigg)^2<\infty
\end{align}
and
\begin{align}\label{eqn:secmomentbound}
{\sup _{n\in\mathbb N}}\mathbb E{\int_0^T}\sum {_{j = 0}^\infty } j^2\pi^n_j(t){\rm d} t<\infty.
\end{align}
Proof. Since πn(t) = πn(0) + An(t) + Mn(t), we can write, for fixed K ∈ ℕ,
\begin{align}\label{eqn:firstmomExpan}
\mathbb E\sup{_{0 \le t \le T}}\Bigg|\sum_{j=0}^K j\pi_j^n(t)\Bigg|^2\! \leq
3\Bigg|\sum_{j=0}^K j\pi_j^n(0)\Bigg|^2\!+3\mathbb E\sup_{0\leq t\leq
T}\Bigg|\sum_{j=0}^K jA^n_j(t)\Bigg|^2\!+3\mathbb E\sup_{0\leq t\leq
T}\Bigg|\sum_{j=0}^K jM^n_j(t)\Bigg|^2\!.
\end{align}
Using (43), for K ∈ ℕ, we can write
\begin{align}\label{eqn:jexpansion3}
\sum_{j=0}^K j\sum_{\ell\in\Sigma}\left\langle\Delta_\ell,
e_j\right\rangle_2\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)} =\sum_{j=1}^K
j[\zeta(j-1,n\pi^n(s))-\zeta(j,n\pi^n(s))] \mathbb N
=\sum_{j=0}^{K-1} \zeta(j,n\pi^n(s))-K\zeta(K,n\pi^n(s))
\end{align}
and
\begin{align}\label{eqn:jexpansion4}
k\sum_{j=0}^K j[\pi_{j+1}^n(s)-\pi_{j}^n(s)]
= -k\Bigg(\sum_{j=1}^K\pi_j^n(s)-K\pi_{K+1}^n(s)\Bigg).
\end{align}
Using similar bounds as in (32), for some c ζ ∈ (0, ∞),
\begin{align}\label{eqn:zetaineq}
\zeta(j,n\pi^n(s))\leq c_{\zeta}n^L \pi_j^n(s).
\end{align}
The above bound implies that, for some κ 1 ∈ (0, ∞) and all n, K ∈ ℕ,
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}\Bigg[\frac{\lambda
}{\binom{n}{L}}\int_0^t\sum_{j=1}^K\zeta(j-1,n\pi^n(s))+k\int_0^t\sum_{j=0}^K\pi_j^n(s){\rm d} s\Bigg]^2
\leq
\mathbb E\Bigg[\Bigg(c_\zeta n^L\frac{\lambda}{\binom{n}{L}}+k\Bigg)T\Bigg]^2
\leq \kappa_1.
\end{align*}
Combined with (45), (52), and (53), the above estimate gives, for all n, K ∈ ℕ,
\begin{align}\label{eq:eq1030}
\mathbb E\sup{_{0 \le t \le T}}\Bigg|\sum_{j=0}^K jA_j^n(t)\Bigg|^2 \le
\kappa_2\Big(1+ K \mathbb E\Big[\sup{_{0 \le t \le T}} (\pi^n_{K}(t) +
\pi^n_{K+1}(t))\Big]\Big).
\end{align}
We now consider 
$\mathbb E\rm sup{_{0 \le {t} \le {T}}}|\sum {_{j = 0}^Kj} M_j^n(t){|^2}$. Since 
$\sum {_{j = 0}^Kj} M_j^n(t)$ is a martingale, Doob’s inequality implies that
\begin{align}\label{eqn:BGD}
\mathbb E\sup{_{0 \le t \le T}}|\sum_{j=0}^K jM^n_j(t)|^2 \leq
4\mathbb E\left\langle\sum_{j=0}^K jM^n_j\right\rangle(T) =
4\mathbb E\sum_{j_1=0}^K\sum_{j_2=0}^K j_1j_2\left\langle M^n_{j_1},M^n_{j_2}\right\rangle(T).
\end{align}
The diagonal terms (j 1 = j 2) in the above sum are given by (26). We now consider the off-diagonal terms. Fix 0 ≤ j 1 < j 2 ≤ K, and note that in order to compute 
$\left\langle M^n_{j_1},M^n_{j_2}\right\rangle(T)$, we must expand
\begin{align}\label{eqn:incoming2}
Z(j_1,j_2,n\pi^n(s))\,{:{=}}\,\sum_{\ell\in\Sigma} \left\langle e_{j_1},\Delta_\ell
\Delta_\ell^T e_{j_2}\right\rangle_2\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}.
\end{align}
Similar to (27), the ℓth term in (57) is the contribution from jobs that request servers with queue length configuration ℓ. A fixed ℓ ∈ Σ will make a nonzero contribution to 〈e j1, 
$\Delta_\ell\Delta_\ell^Te_{j_2}$〉2 if (j 1 or j 1 − 1) and (j 2 or j 2 − 1) are among the k-smallest coordinates in ℓ. That is, for a fixed ℓ ∈ Σ, the ℓth term is nonzero only if (j 1 or j 1 − 1) is a member of the set (ℓ 1, …, ℓ k) and (j 2 or j 2 − 1) is also a member. The contribution from all such ℓ in the sum (57) can be counted in a method analogous to that used to obtain (30). Namely, we count the numbers of choices of servers with queue length less than j 1 − 1, equal to j 1 − 1, equal to j 1, between j 1 and j 2 − 1, equal to j 2 − 1, equal to j 2, and larger than j 2. One must be careful in the cases j 2 − 1 = j 1 and j 2 − 1 = j 1 + 1. In both cases there are no servers with length between j 1 and j 2 − 1. In the first case above (j 2 − 1 = j 1), we must also be careful not to double count. To ensure this, we include an indicator function 1{j2>j1+1} in the upper index of the binomial coefficient corresponding to the selection of servers with queue length equal to j 2 − 1. Combining these observations we see that, for x ∈ n𝒮n,
\begin{align}\label{eqn:Zdef4}
Z({j_1},{j_2},x) = \sum\limits_{\ell \in \Sigma } {\langle {e_{{j_1}}},{\Delta _\ell }\Delta _\ell ^T{e_{{j_2}}}\rangle _2}\prod\limits_{i = 0}^\infty \left({\matrix{
{{x_i}} \cr
{{\rho _i}(\ell)} \cr
} } \right)
= \sum\limits_{{i_1} = 0}^{k - 2} (\matrix{
{\sum\limits_{m = 0}^{{j_1} - 2} {x_m}} \cr
{{i_1}} \cr
})\sum\limits_{{i_2} = 0}^{k - {i_1} - 1} \left({\matrix{
{{x_{{j_1} - 1}}} \cr
{{i_2}} \cr
} } \right)\sum\limits_{{i_3} = 0}^{k - {i_1} - {i_2} - 1} [{i_2} - {i_3}]\left({\matrix{
{{x_{{j_1}}}} \cr
{{i_3}} \cr
} } \right)
\times \sum\limits_{{i_4} = 0}^{k - {i_1} - {i_2} - {i_3} - 1} (\matrix{
{\sum\limits_{m = {j_1} + 1}^{{j_2} - 2} {x_m}} \cr
{{i_4}} \cr
})\sum\limits_{{i_5} = 0}^{L - \sum\limits_{n = 1}^4 {i_n}} \left({\matrix{
{{x_{{j_2} - 1}}{1_{\{ {j_2} > {j_1} + 1\} }}} \cr
{{i_5}} \cr
} } \right)
\times \sum\limits_{{i_6} = 0}^{L - \sum\limits_{n = 1}^5 {i_n}} {[(1_{\{ {j_2} = {j_1} + 1\} }}({i_3} - {i_5}) + {i_5}) \wedge {(k - \sum\limits_{n = 1}^4 {i_n})_ + } - {i_6} \wedge {(k - \sum\limits_{n = 1}^5 {i_n})_ + }]
\times \left({\matrix{
{{x_{{j_2}}}} \cr
{{i_6}} \cr
} } \right)(\matrix{
{\sum\limits_{m = {j_2} + 1}^\infty {x_m}} \cr
{L - \sum\limits_{n = 1}^6 {i_n}} \cr
}).
\end{align}
For j 1 > j 2, we define Z(j 1, j 2, x) := Z(j 2, j 1, x). The contribution to 
$\left\langle M^n_{j_1},M^n_{j_2}\right\rangle(T)$ for j 1 ≠ j 2 from completed jobs is given by the term
\begin{align}\label{eqn:outgoing2}
\sum_{i=1}^\infty\left\langle
e_{j_1},(e_{i-1}-e_i)(e_{i-1}-e_i)^Te_{j_2}\right\rangle 2\pi_i^n(s) =-{\bf
1}_{\{j_1=j_2-1\}}\pi^n_{j_2}(s)-{\bf 1}_{\{j_1-1=j_2\}}\pi^n_{j_1}(s).
\end{align}
This follows on noting that if a job is completed from a queue of length j then its queue length becomes j − 1. This implies that the contribution is zero unless j 1 = j 2 − 1 or j 1 − 1 = j 2, which results in the above expression. Combining (58) and (59) gives, for j 1, j 2 ∈ ℕ0,
\begin{align}\label{eqn:Zdef2}
\left\langle M^n_{j_1},M^n_{j_2}\right\rangle(T)
=\frac{\lambda}{n\binom{n}{L}}\int_0^TZ(j_1,j_2,n\pi^n(s)){\rm d} s
+ {k \over n}{\int_0^T}[{\bf 1}_{\{j_1=j_2\}}[\pi^n_{j_1}(s)+\pi^n_{j_1+1}]-{\bf1}_{\{j_1=j_2-1\}}\pi^n_{j_2}(s)
{\bf 1}_{\{j_1-1=j_2\}}\pi^n_{j_1}(s)]{\rm d} s,
\end{align}
where, by convention, Z(j, j, x) : = Z(j, x). Referring to the definition of Z in (58), note that, for j 2 > j 1 + 1, Z(j 1, j 2, x) = 0 unless (i 2 or i 3) is greater than 0 and (i 5 or i 6) is greater than 0. In the case that j 2 = j 1 + 1, Z(j 1, j 2, x) = 0 unless (i 2 or i 3) is greater than 0 and (i 3 or i 6) is greater than 0. Therefore, (31) implies there exists a c̃ Z ∈ (0, ∞) such that, for r ∈ 𝒮n and j 1 < j 2,
\begin{align}\label{eqn:offDiagBound}
Z(j_1,j_2,nr)\leq\tilde
c_Zn^L[r_{j_1}r_{j_2}+r_{j_1-1}r_{j_2}+r_{j_1}r_{j_2-1}+r_{j_1-1}r_{j_2-1}+{\bf
1}_{\{j_2=j_1+1\}}r_{j_1}].
\end{align}
Combining this with (32) and (60), we have, for some 
${\kappa '_3}$, κ 3 ∈ (0, ∞) and all n, K ∈ ℕ,
\begin{align}\label{eqn:crossquad}
\sum_{j_1=0}^K\sum_{j_2=0}^K j_1j_2\left\langle M^n_{j_1},M^n_{j_2}\right\rangle(T)
\le \frac{\kappa_3'}{n}\bigg[\int_0^T\sum_{j_1=0}^K\sum_{j_2=0}^K
(j_1+1)(j_2+1)\pi^n_{j_1}(t)\pi^n_{j_2}(t){\rm d} t+\int_0^T\sum_{j=1}^{K+1}
j(j+1)\pi^n_j(t){\rm d} t\bigg]
\leq\frac{\kappa_3}{n}\bigg[\int_0^T\bigg(\sum_{j=0}^{K}
j^2\pi^n_j(t)+(K+1)^2\pi^n_{K+1}(t)\bigg){\rm d} t+1\bigg].
\end{align}
Recalling that πn(t) = πn(0) + An(t) + Mn(t), we have, for all K, n ∈ ℕ,
\begin{align*}
\mathbb E\int_0^T\sum_{j=0}^K j^2\pi^n_j(t){\rm d} t = \int_0^T\sum_{j=0}^K
j^2\pi^n_j(0){\rm d} t+\mathbb E\int_0^T\sum_{j=0}^K j^2A^n_j(t){\rm d} t+
\int_0^T\mathbb E\sum_{j=0}^K j^2M^n_j(t){\rm d} t
\le \mathbb E\int_0^T\sum_{j=0}^K j^2A^n_j(t){\rm d} t + \kappa_4,
\end{align*}
where κ 4 = c π(0)T and the last inequality follows on using the fact that 
$M^n_j(t)$ is a martingale. Thus, from (45), for some κ 5 ∈ (0, ∞) and all K, n ∈ ℕ,
\begin{align}\label{eqn:piSecMoment}
\mathbb E{\int_0^T}\sum_{j=0}^K j^2\pi^n_j(t){\rm d} t \le
\frac{\kappa_5}{n^L}\mathbb E{\int_0^T}\sum_{j=1}^K
j^2{\int_0^T}[\zeta(j-1,n\pi^n(s))-\zeta(j,n\pi^n(s))]{\rm d} s{\rm d} t \mathbb N
+\kappa_5\mathbb E{\int_0^T}\sum_{j=1}^K
j^2{\int_0^T}[\pi^n_{j+1}(s)-\pi^n_j(s)]{\rm d} s{\rm d} t+\kappa_5.
\end{align}
Using the facts that, for any a 0, …, a K ∈ ℝ,
\begin{align*}
\sum\limits_{j = 1}^K {{j^2}\left[ {{a_{j - 1}} - {a_j}} \right] = \sum\limits_{j = 1}^K {\left[ {{{\left({j - 1} \right)}^2}{a_{j - 1}} - {j^2}{a_j} + \left({{2_j} - 1} \right){a_{j - 1}}} \right]} }
= -K^2a_K+\sum_{j=1}^K(2j-1)a_{j-1}
\end{align*}
and
\begin{align*}
\sum_{j=0}^Kj^2[a_{j+1}-a_j] = \sum\limits_{j = 0}^K {\left[ {{{\left({j + 1} \right)}^2}{a_{j + 1}} - {j^2}{a_j} - \left({2j + 1} \right){a_{j + 1}}} \right]}
=(K+1)^2a_{K+1}-\sum_{j=0}^K(2j+1)a_{j+1}
\end{align*}
in (63) we have, for some κ 6 ∈ (0, ∞) and all K, n ∈ ℕ,
\begin{align}\label{eqn:piSecMoment2}
\mathbb E{\int_0^T}\sum_{j=0}^K j^2\pi^n_j(t){\rm d} t
\leq\frac{\kappa_5}{n^L}\mathbb E{\int_0^T}{\int_0^T}\sum_{j=0}^K (2j-1)\zeta(j-1,n\pi^n(s)){\rm d} s\,\;{\rm d} t
+\kappa_5\mathbb E{\int_0^T}\int_0^t (K+1)^2\pi^n_{K+1}(s){\rm d} s\,\;{\rm d} t+\kappa_5
\leq \kappa_6\mathbb E{\int_0^T} \bigg[K^2\sup_{0\le s \le t}\pi^n_{K+1}(s)+\sup_{0\le s \le t}\sum_{j=0}^K j\pi^n_{j}(s)\bigg]{\rm d} t+\kappa_6,
\end{align}
where the second inequality follows from (54). Thus, it follows from (56) and (62) that, for some κ 7 ∈ (0, ∞),
\begin{align}\label{eqn:MsecMoment}
\mathbb E\sup{_{0 \le t \le T}}|\sum_{j=0}^K jM^n(t)|^2 \leq
\frac{\kappa_3}{n}\bigg[{\int_0^T}\mathbb E\sum_{j=0}^K
j^2\pi^n_j(t){\rm d} t+\gamma^n_{K}T+1\bigg]
\leq \frac{\kappa_7}{n}\bigg[1 +\gamma^n_K + {\int_0^T}\mathbb E\sup_{0\leq u\leq
s}|\sum_{j=0}^K j\pi^n_j(u)|^2{\rm d} s\bigg],
\end{align}
where 
$\gamma _K^n = \mathbb E({K^2}{\sup _{0 \le s \le T}}\pi _{K + 1}^n(s))$. Combining (51), (55), and (65), and using the fact that 
$\big|\sum {_{j = 0}^\infty } j\pi^n_j(0)\big|\leq c_{\pi(0)}$,
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}|\sum_{j=0}^K j\pi^n_j(t)|^2\leq
\kappa_8 \bigg(1 + \mathbb E\sup{_{0 \le t \le T}} |\sum_{j=0}^K
jA^n_j(t)|^2 + \mathbb E\sup{_{0 \le t \le T}} |\sum_{j=0}^K
jM^n_j(t)|^2 \bigg)
\le \kappa_9\bigg(1 + \gamma^n_K + \frac{1}{n}{\int_0^T}\mathbb E\sup_{0\leq s\leq
t}|\sum_{j=0}^K j\pi^n_j(s)|^2 {\rm d} s\bigg).
\end{align*}
By Gronwall’s lemma (since the above inequality also holds for all T 1 ≤ T), there is a κ 10 ∈ (0, ∞) such that, for all n, K ∈ ℕ,
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}|\sum_{j=0}^K j\pi^n_j(t)|^2\leq \kappa_{10}(1 + \gamma^n_K).
\end{align*}
Sending K → ∞ and recalling from Lemma 2 that, for each fixed n, as K → ∞, 
$\gamma^n_K \to 0$ we have, for all n,
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}\bigg|\sum_{j=0}^{\infty} j\pi^n_j(t)\bigg|^2\leq\kappa_{10},
\end{align*}
where κ 10 is independent of n. This proves (49). Finally, (50) follows from (49) upon sending K → ∞ in (64).
5.2. Tightness
We now proceed with the proof of the tightness of {(Xn, Mn ,c )}n∈ℕ. Let, for M ∈ ℝ+,
\begin{align*}
\mathcal V_M \,{:{=}}\, \bigg\{r\in\mathcal S|\sum_{i=0}^\infty ir_i\leq M\bigg\},
\end{align*}
where 𝒱M is equipped with the topology inherited from ℓ 2. We begin by establishing the following Lipschitz property for F on 𝒱M.
Lemma 4
The map F is a Lipschitz function from 𝒱M to ℓ 2 for each M ∈ ℝ+. Namely, there exists an C(M) ∈ (0, ∞) such that, for any r, r̃ ∈ 𝒱M,
\begin{align}\label{eqn:LipschitzDisplay}
\|F(r)-F(\tilde r)\|_2\leq C(M)\|r-\tilde r\|_2.
\end{align}
Proof. Fix M ∈ ℝ+. Let r, r̃ ∈ 𝒱M and, for i 1 ∈ ℕ0 and j, i 2 ∈ ℕ, recall R j,i 1,i 2(r, r̃) from (35). Using (36) and the fact that r, r̃ ∈ 𝒮, we have
\begin{align*}
(R_{j,i_1,i_2}(r,\tilde r))^2 \leq 3[r_{j}^{i_2}-\tilde r_{j}^{i_2}]^2+3\tilde r_{j}^{2i_2}\bigg[\bigg(\sum_{m=j+1}^\infty
r_m\bigg)^{L-i_1-i_2}-\bigg(\sum_{m=j+1}^\infty \tilde
r_m\bigg)^{L-i_1-i_2}\bigg]^2
+3\tilde r_{j}^{2i_2}\bigg[\bigg(\sum_{m=0}^{j-1}
r_m\bigg)^{i_1}-\bigg(\sum_{m=0}^{j-1} \tilde r_m\bigg)^{i_1}\bigg]^2.
\end{align*}
By an argument similar to that used to derive (38) and an application of the Cauchy–Schwarz inequality we have the following inequality for all i 1, i 2 ≤ L, i 2 > 0:
\begin{align}\label{eqn:Rineq2}
(R_{j,i_1,i_2}(r,\tilde r))^2\leq \kappa_1([r_{j}-\tilde r_{j}]^2+(j+1)\tilde r_{j}\|r-\tilde r\|_2^2).
\end{align}
Using arguments analogous to those used in the proof of Lemma 1, we have
\begin{align}\label{eqn:Fineq2}
\|F(r)-F(\tilde r)\|_2 \leq \kappa_2\lambda
L\text{!}\bigg(\sum {_{j = 0}^\infty }\sum_{i_1=0}^{k-1}\sum_{i_2=1}^{L-i_1}
[R_{j,i_1,i_2}(r,\tilde r)]^2\bigg)^{1/2}+ 2k\bigg(\sum {_{j = 0}^\infty }(r-\tilde r)_j^2\bigg)^{1/2}
\leq \kappa_3\bigg(\sum {_{j = 0}^\infty }[[r_{j}-\tilde r_{j}]^2 +(j+1)\tilde r_{j}\|r-\tilde r\|_2^2]\bigg)^{1/2}+ 2k\|r-\tilde r\|_2
\leq \kappa_4\|r-\tilde r\|_2\bigg(1 +\sum {_{j = 0}^\infty } j\tilde r_{j}\bigg)^{1/2}+2k\|r-\tilde r\|_2.
\end{align}
Since r, r̃ ∈ 𝒱M, (68) gives
\begin{align*}
\|F(r)-F(\tilde r)\|_2 \leq \kappa_4(M+1)^{1/2}\|r-\tilde r\|_2+ 2k\|r-\tilde r\|_2,
\end{align*}
and, thus, with C(M) : = κ4(M + 1)1/2 + 2k, (66) is satisfied for all r, r̃ ∈ 𝒱M, which proves the result.
Recall the process Xn introduced in (10) and Mn,c defined below (24). The following proposition gives tightness of {(Xn, Mn,c)}n∈ℕ.
Proposition 4
Suppose that {πn}n∈ℕ is as in the statement of Theorem 1 with 
${\sup _{n \in\mathbb N}}\sum {_{j = 0}^\infty } {j^2}\pi _j^n(0) < \infty$. Let 
$X^n(0)=\sqrt{n}(\pi^n(0)-\pi_0)$ and suppose that (17) is satisfied. Then {(Xn, Mn,c)}n∈ℕ is a ℂ-tight sequence of 𝔻([0, T]: (ℓ̃ 2)2)-valued random variables.
Proof. We will make use of Theorem 4 in Appendix A.2. We first prove that {Mn,c}n∈ℕ is tight. In order to show that condition (A) in Theorem 4 is satisfied for {Mn,c}n∈ℕ it suffices (cf. Theorem 2.3.2 of [Reference Joffe and Métivier20]) to show that the condition is satisfied for the real-valued process 
$\left\langle M^{n,c}\right\rangle(t)\,{:{=}}\, \sum {_{j = 0}^\infty }\left\langle M^{n,c}_j\right\rangle(t)$. Fix ε ∈ (0, T] and T 0 ∈ [0, T − ε]. Let τ n ≤ T 0 be a sequence of 
$\{\mathcal F^n_t\}$-stopping times. Then, (43) and (54) imply that, for θ ∈ [0, ε],
\begin{align*}
|\left\langle M^{n,c}(\tau_n+\theta)\right\rangle-\left\langle M^{n,c}(\tau_n)\right\rangle|
=|\sum {_{j = 0}^\infty }\biggl[{\int_{\tau_n}^{\tau_n+\theta}}\sum_{\ell\in\Sigma}\left\langle\Delta_\ell,e_j\right\rangle_2
\frac{\lambda}{I(n)}\prod_{i=0}^\infty\binom{n\pi^n_i(s)}{\rho_i(\ell)}+k{\int_{\tau_n}^{\tau_n+\theta}}
[\pi_{j+1}^n(s)-\pi_j^n(s)]{\rm d} s\biggr]|
\leq\kappa_{1}\sum {_{j = 0}^\infty }{\int_{\tau_n}^{\tau_n+\theta}}[\pi_j^n(s)+\pi_{j-1}^n(s)
+ \pi_{j+1}^n(s)]{\rm d} s
\leq \kappa_{1}\varepsilon.
\end{align*}
The proof of (A) is now immediate.
We next show that {Mn,c}n∈ℕ satisfies condition (T1) of Theorem 4. For this, we will apply Theorem 3. We first verify that {Mn,c(t)}n∈ℕ satisfies Theorem 3(a) for all t ∈ [0, T]. It follows from (33) that
\begin{align}\label{eqn:24.9}
{\sup _{n\in\mathbb N}}\mathbb E\left\langle M^{n,c}\right\rangle(T)={\sup _{n\in\mathbb N}}n\mathbb E\left\langle M^n\right\rangle(T) \leq\kappa_2.
\end{align}
This, combined with Doob’s inequality, implies that, for each n 0 ∈ ℕ,
\begin{align*}
{\sup _{n\in\mathbb N}}\sum_{i=0}^{n_0}\mathbb E\sup{_{0 \le t \le T}}\left|M_i^{n,c}(t)\right|
\leq n_0+{\sup _{n\in\mathbb N}}\sum_{i=0}^{n_0}\mathbb E \Big(\sup{_{0 \le t \le T}}M_i^{n,c}(t)\Big)^2
\leq n_0+\kappa_3.
\end{align*}
Using Markov’s inequality, Theorem 3(a) follows.
We now verify Theorem 3(b) for {Mn,c(t)}n∈ℕ for each fixed t ∈ [0, T]. Note that 
$\left\langle M^{n,c}_j\right\rangle(t) = n\left\langle M^n_j\right\rangle(t)$ and, thus, from (26) and (32),
\begin{align}\label{eqn:tripleStar}
\left\langle M^{n,c}_j\right\rangle(t)\leq
\kappa_4\int_0^t(\pi^n_{j-1}(s)+\pi^n_{j}(s)+\pi^n_{j+1}(s)){\rm d} s.
\end{align}
It follows from (70) and the Cauchy–Schwarz inequality that
\begin{align}\label{eqn:MsqSum}
\sum_{j=n_0}^\infty \mathbb E(M_j^{n,c}(t))^2 &= \sum_{j=n_0}^\infty \mathbb E\left\langle
M_j^{n,c}(t)\right\rangle \mathbb N
\\
& \leq \kappa_5\mathbb E\int_0^t\sum_{j=n_0-1}^\infty\pi^n_j(s){\rm d} s\mathbb N
\\
\leq
\kappa_5\bigg(\sum_{j=n_0-1}^\infty\frac{1}{j^2}\bigg)^{1/2}\int_0^t\mathbb E\bigg(\sum_{j=n_0-1}^\infty
j^2(\pi^n_j(s))^2\bigg)^{1/2}{\rm d} s\kern-0.2pt.
\end{align}
From (50),
\begin{align}\label{eqn:piSecMoment3}
{\sup _{n\in\mathbb N}}\mathbb E{\int_0^T}\sum {_{j = 0}^\infty } j^2(\pi^n_j(s))^2{\rm d} s
\leq{\sup _{n\in\mathbb N}}\mathbb E{\int_0^T}\sum {_{j = 0}^\infty } j^2\pi_j^n(s){\rm d} s \,{{=}:}\,
\kappa_6 <\infty.
\end{align}
Using this observation in (71), we have
\begin{align*}
\sum_{j=n_0}^\infty \mathbb E(M_j^{n,c}(t))^2 \leq
\kappa_7\bigg(\sum_{j=n_0-1}^\infty\frac{1}{j^2}\bigg)^{1/2}\int_0^t\mathbb E\bigg(\sum_{j=n_0-1}^\infty
j^2\pi^n_j(s)\bigg)^{1/2}{\rm d} s
\leq \kappa_8\bigg(\sum_{j=n_0-1}^\infty\frac{1}{j^2}\bigg)^{1/2}.
\end{align*}
From Markov’s inequality we now see that, for any δ > 0,
\begin{align*}
\lim_{n_0\to\infty}{\sup _{n\in\mathbb N}}\mathbb P\bigg(\sum_{j=n_0}^\infty(M^{n,c}_j(t))^2>\delta\bigg)=0,
\end{align*}
which verifies Theorem 3(b). Thus, we have shown that {Mn,c(t)}n∈ℕ is a tight sequence of ℓ 2-valued random variables for all t ∈ [0, T]. From Theorem 4, it now follows that {Mn,c}n∈ℕ is a tight sequence of 𝔻([0, T]: ℓ 2)-valued random variables.
We will now argue that {Xn}n∈ℕ is a tight sequence of 𝔻([0, T]: ℓ 2)-valued random variables. Again, via Theorem 4, it suffices to show that {Xn(t)}n∈ℕ is tight for every t ∈ [0, T] (which will follow from verifying conditions (a) and (b) of Theorem 3) and that {Xn}n∈ℕ satisfies condition (A) of Theorem 4. We first show that, for all t ∈ [0, T], Theorem 3(a) holds for {Xn(t)}n∈ℕ. Namely, we show that, for each n 0 ∈ ℕ and t ∈ [0, T],
\begin{align}\label{eqn:25b}
\lim_{A\to\infty}{\sup _{n\in\mathbb N}}\mathbb P\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|>A\bigg)=0.
\end{align}
Fix ε > 0. From Lemma 3, there is a M ∈ (0, ∞) such that
\begin{align}\label{eqn:epsM}
{\sup _{n\in\mathbb N}}\mathbb E\bigg(\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty }j\pi_j^n(t)\bigg)\leq \frac{M\varepsilon}{2}.
\end{align}
Let 
$B_M^n\,{:{=}}\,\{\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j\pi^n_j(t)\leq M\}$. Then, for t ∈ [0, T] and n 0 ∈ ℕ,
\begin{align}\label{eqn:probaexpansion}
\mathbb P\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|>A\bigg)
\leq \mathbb P\bigg(\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j\pi^n_j(t)>M\bigg)
+\mathbb P\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|>A,B_M^n\bigg)\mathbb N
\leq \frac{\varepsilon}{2}+\mathbb P\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|>A,B_M^n\bigg).
\end{align}
The Cauchy–Schwarz inequality yields
\begin{align}\label{eqn:Xexpansion}
\sum_{i=0}^{n_0}|X^n_j(t)|
\leq \sqrt{n_0}\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|^2\bigg)^{1/2}
\leq \sqrt{n_0}\|X^n(t)\|_2.
\end{align}
Furthermore, from (20) and the triangle inequality,
\begin{align}\label{eqn:54b}
\|X^n(t)\|_2 \leq \|X^n(0)\|_2+\|A^{n,c}(t)\|_2+\|M^{n,c}(t)\|_2.
\end{align}
The definition of An,c in (24) gives
\begin{align*}
\|A^{n,c}(t)\|_2=\sqrt{n}\left\|A^n(t)-\int_0^tF(\pi(s)){\rm d} s\right\|_2.
\end{align*}
The moment bound (49) proved in Lemma 3 implies that
\begin{align}\label{eqn:firtMomSquaredBound}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}|\sum {_{j = 0}^\infty }
j\pi^n_j(t)|^2=\!:\,\kappa_7<\infty,
\end{align}
and, thus, for some κ 8 ∈ (0, ∞),
\begin{align}\label{eqn:firstMomSquaredBoundLim}
\sup{_{0 \le t \le T}}|\sum {_{j = 0}^\infty } j\pi_j(t)|^2\leq\kappa_8
\end{align}
as well. From (48) and the Lipschitz property proved in Lemma 4, with M ≥ κ 7 ˅ κ 8 on the set 
$B^n_M$,
\begin{align*}
\|A^{n,c}(t)\|_2 \leq \sqrt{n}\int_0^t\|F(\pi^n(s))-F(\pi(s))\|_2{\rm d} s +
\frac{\kappa_9}{\sqrt{n}\,} \leq C(M)\int_0^t\|X^n(s)\|_2{\rm d} s +
\frac{\kappa_9}{\sqrt{n}\,}.
\end{align*}
Thus, from (77) and Gronwall’s lemma, on the set 
$B_M^n$, for all n ≥ 1,
\begin{align}\label{eqn:Xngronwall}
\sup{_{0 \le t \le T}}\|X^n(t)\|_2\leq
\kappa_{10}\bigg(\frac{1}{\sqrt{n}\,}+\|X^n(0)\|_2+\sup_{0\leq t\leq
T}\|M^{n,c}(t)\|_2\bigg){\rm e}^{C(M)T}.
\end{align}
From (69) and Doob’s inequality,
\begin{align}\label{eqn:Mnunifsecmomentbound}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\|M^{n,c}(t)\|^2_2<\infty.
\end{align}
Also, by assumption, Xn(0) → x 0 in ℓ 2. Thus, for the given ε > 0, we can find α 0 such that, for all α ≥ α 0,
\begin{align*}
\mathbb P\bigg(\sup{_{0 \le t \le T}}\|X^n(t)\|_2\geq \frac{\alpha}{\sqrt{n_0}\,},
B^n_M\bigg)\leq\frac{\varepsilon}{2}.
\end{align*}
Therefore, from (75) and (76) we have, for all 
$A\geq {\alpha_0}/{\sqrt{n_0}}$,
\begin{align*}
{\sup _{n\in\mathbb N}}\mathbb P\bigg(\sum_{j=0}^{n_0}|X^n_j(t)|>A\bigg)\leq
\frac{\varepsilon}{2}+\frac{\varepsilon}{2}=\varepsilon.
\end{align*}
Since ε > 0 is arbitrary we get (73). Thus, we have verified Theorem 3(a) for {Xn(t)}n∈ℕ for each t ∈ [0, T].
We now consider Theorem 3(b). Namely, we show that, for every δ > 0 and t ∈ [0, T],
\begin{align*}
\lim_{n_0\to\infty}{\sup _{n\in\mathbb N}}\mathbb P\bigg(\sum_{j=n_0}^\infty(X^n_j(t))^2>\delta\bigg)=0.
\end{align*}
For this, it suffices to show that
\begin{align}\label{eqn:XnsecondMoment}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(X^n_j(t))^2<\infty.
\end{align}
Recalling that 
$X^n_j(t) = X^n_j(0)+A^{n,c}_j(t)+M^{n,c}_j(t)$ for each j ∈ ℕ, it follows that
\begin{align}\label{eqn:XsecMomentExp}
\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(X^n_j(t))^2 \leq
3{\sup _{n\in\mathbb N}}\sum {_{j = 0}^\infty } j^2(X^n_j(0))^2+3{\sup _{n\in\mathbb N}}\mathbb E\sup_{0\leq
t\leq T}\sum {_{j = 0}^\infty } j^2(A^{n,c}_j(t))^2 \mathbb N
+3{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty }
j^2(M^{n,c}_j(t))^2.
\end{align}
Using the definitions of An,c, An, and F in (24), (45), and (5), respectively, we can write
\begin{align}\label{eqn:AsecMoment1}
(A^{n,c}_j(t))^2 \leq
\kappa_{11}\biggl\{\int_0^tn\Big[\frac{(n-L)\text{!}}{n\text{!}}\zeta(j,n\pi^n(s))-\bar{\zeta}(j,\pi(s))\Big]^2{\rm d} s
+\int_0^tn\Big[\frac{(n-L)\text{!}}{n\text{!}}\zeta(j-1,n\pi^n(s))-\bar{\zeta}(j-1,\pi(s))\Big]^2{\rm d} s \mathbb N
+ \int_0^tn[\pi^n_j(s)-\pi_j(s)]^2{\rm d} s+\int_0^tn[\pi^n_{j+1}(s)-\pi_{j+1}(s)]^2{\rm d} s\biggr\}.
\end{align}
From (47) and in a similar manner as in (40) we have
\begin{align*}
n\Big[\frac{(n-L)\text{!}}{n\text{!}}\zeta(j,n\pi^n(s))-\bar{\zeta}(j,\pi(s))\Big]^2 \leq
\kappa_{12}\{(\pi^n_j(s))^2+n[\bar{\zeta}(j,\pi^n(s))-\bar{\zeta}(j,\pi(s))]^2\}
\leq
\kappa_{13}\Bigg\{(\pi^n_j(s))^2+n\sum_{i_1=0}^{k-1}\sum_{i_2=1}^{L-i_1}R_{j,i_1,i_2}(\pi^n(s),\pi(s))^2\Bigg\},
\end{align*}
where R j,i 1,i 2 is as in (35). By (67) and the Cauchy–Schwarz inequality we now have
\begin{align*}
nR_{j,i_1,i_2}(\pi^n(s),\pi(s))^2 \leq
\kappa_{14}\Bigg[(X^n_j(s))^2+\pi_j(s)\Bigg(\sum_{m=0}^\infty|X^n_m(s)|\Bigg)^2\Bigg]
\leq
\kappa_{14}\Bigg[(X^n_j(s))^2+\pi_j(s)\Bigg(\sum_{m=0}^\infty\frac{1}{m^2}\Bigg)\sum_{m=0}^\infty
m^2(X^n_m(s))^2\Bigg].
\end{align*}
Therefore,
\begin{align*}
n\Big[\frac{(n-L)\text{!}}{n\text{!}}\zeta(j,n\pi^n(s))-\bar{\zeta}(j,\pi(s))\Big]^2 \!\leq
\kappa_{15}\Bigg\{(\pi^n_j(s))^2+(X^n_j(s))^2+\pi_j(s)\sum_{m=0}^\infty
m^2(X^n_m(s))^2\Bigg\}.
\end{align*}
Combining this estimate with (72) and (84) yields
\begin{align}\label{eqn:AsecMoment2}
\mathbb E\sum {_{j = 0}^\infty } j^2(A^{n,c}_j(t))^2 \leq
\kappa_{16}\mathbb E\Bigg\{\int_0^t\sum_{j=1}^\infty
j^2\Bigg[(X^n_{j-1}(s))^2+(X^n_j(s))^2+(X^n_{j+1}(s))^2
+(\pi_j(s)+\pi_{j-1}(s))\sum_{m=0}^\infty
m^2(X^n_m(s))^2\Bigg]{\rm d} s\Bigg\}+\kappa_{16}
\leq \kappa_{17}\mathbb E\int_0^t\biggl(1+\sum_{j=1}^\infty
j^2\pi_j(s)\biggr)\biggl(\sum_{j=1}^\infty
j^2(X^n_j(s))^2\biggr){\rm d} s+\kappa_{17}.
\end{align}
Additionally, it follows from (70) and (50) that
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2M_j^{n,c}(t)^2 \leq
\kappa_{18}'\mathbb E\sum {_{j = 0}^\infty } j^2\left\langle M_j^{n,c}\right\rangle_T \leq
\kappa_{18}{\int_0^T}\biggl[1+\mathbb E\sum {_{j = 0}^\infty } j^2\pi_j^n(s)\biggr]{\rm d} s \leq
\kappa_{19}.
\end{align*}
Therefore, from (17), (83), and (85), for all t ∈ [0, T],
\begin{align*}
\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(X^n_j(t))^2 \leq \kappa_{20}+
\kappa_{20}{\int_0^T}\bigg(1+\sum_{j=1}^\infty j^2\pi_j(t)\bigg)\mathbb E\sup_{0\leq
s\leq t}\bigg(\sum_{j=1}^\infty j^2(X^n_j(s))^2\bigg){\rm d} t.
\end{align*}
From (50) and Fatou’s lemma, 
${\int_0^T} {\sum {_{j = 1}^\infty {j^2}{\pi _j}\left(s \right){\rm{d}}s} } < \infty$, and, thus, by Gronwall’s lemma,
\begin{align*}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(X^n_j(t))^2 \leq
\kappa_{19} \exp\bigg[\kappa_{20}{\int_0^T}\bigg(1+\sum_{j=1}^\infty
j^2\pi_j(s)\bigg){\rm d} s\bigg] <\infty.
\end{align*}
This proves (82) and verifies Theorem 3(b) for {Xn(t)}n∈ℕ for each t ∈ [0, T]. Thus, {Xn(t)}n∈ℕ is a tight sequence of ℓ 2-valued random variables for every t ∈ [0, T].
We now show that condition (A) of Theorem 4 holds for {Xn}n∈ℕ. Since Xn(t) = Xn(0) + An,c(t) + Mn,c(t) and we have shown the condition is satisfied by {Mn,c}n∈ℕ, it suffices to show that the condition holds for {An,c}n∈ℕ. Let T 0, η, ε, θ > 0, T 0 ≤ T − θ, and suppose that {τ n}n∈ℕ is a family of stopping times such that τ n ≤ T 0. From the definition of An,c (cf. (24)) and (48), we have
\begin{align}\label{eqn:AbarFluc}
\|A^{n,c}(\tau_n+\theta)-A^{n,c}(\tau_n)\|_2 \leq
\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}\|F(\pi^n(t))-F(\pi(t))\|_2{\rm d} t +\frac{\kappa_{21}}{\sqrt{n}\,},
\end{align}
where κ 21 is independent of the choices of τ n and T 0. Fix n 0 ∈ ℕ such that 
$\eta -{\kappa_{21}}/{\sqrt{n_0}}>0$ and let 
$\eta' = \eta -{\kappa_{21}}/{\sqrt{n_0}}$. Recall κ 7 and κ 8 introduced in (78) and (79), and 
$B_M^n$ introduced below (74). Select M ∈ (0, ∞) large enough that M > κ 7 ˅ κ 8 and (74) holds. Then, for all n ≥ n 0,
\begin{align}\label{eqn:Xnfluctuation1}
\mathbb P \biggl\{\biggl\|\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}[F(\pi^n(t))-F(\pi(t))]{\rm d} t\biggr\|_2>\eta'\biggr\} \mathbb N\leq \mathbb P
\biggl\{\biggl\|\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}[F(\pi^n(t))-F(\pi(t))]{\rm d} t\biggr\|_2>\eta',B_M^n\biggr\}+\mathbb P \biggl\{\sup_{0\leq t\leq
T}\sum {_{j = 0}^\infty } j\pi_j^n(t)>M\biggr\}\mathbb N\leq \mathbb P
\biggl\{\biggl\|\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}[F(\pi^n(t))-F(\pi(t))]{\rm d} t\biggr\|_2>\eta',B_M^n\biggr\}+\frac{\varepsilon}{2}.
\end{align}
It follows from the Lipschitz property of F proved in Lemma 4 that
\begin{align}\label{eqn:Xnfluctuation2}
\mathbb P \bigg\{\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}\|F(\pi^n(t))-F(\pi(t))]\|_2{\rm d} t
>\eta',B_M^n\bigg\}\leq \mathbb P \bigg\{C(M)\int_{\tau_n}^{\tau_n+\theta}\|X^n(t)\|_2{\rm d} t>\eta',B_M^n\bigg\}.
\end{align}
Recall from (80) that, for some C̃(M) ∈ (0, ∞) on the set 
$B_M^n$,
\begin{align*}
C(M)\sup{_{0 \le t \le T}}\|X^n(t)\|_2\leq \tilde C(M)\bigg(1+\sup{_{0 \le t \le T}}\|M^{n,c}(t)\|_2\bigg).
\end{align*}
Thus, from (88), Markov’s inequality, and (81), we have
\begin{align}\label{eqn:Xnfluc5}
\mathbb P
\Bigg\{\!\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}\|F(\pi^n(t))-F(\pi(t))]\|_2{\rm d} t>\eta',B_M^n\Bigg\}\leq \mathbb P\Big\{\theta\tilde C(M)\Big(1+\sup{_{0 \le t \le T}}\|M^{n,c}(t)\|_2\Big)>\eta'\Big\}
\mathbb N
\leq {\theta \tilde C(M)(1+\mathbb E\sup{_{0 \le t \le T}}\|M^{n,c}(t)\|_2)}{\eta'}\mathbb N
\leq \theta\tilde C(M)\kappa_{22}.
\end{align}
Combining (87) and (89) gives, whenever θ ≤ δ,
\begin{align*}
\sup_{0\leq\theta\leq\delta}\mathbb P
\Bigg\{\!\left\|\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}[F(\pi^n(t))-F(\pi(t))]{\rm d} t\right\|_2>\eta'\Bigg\} \leq C(M)\kappa_{22}\delta+ \frac{\varepsilon}{2}.
\end{align*}
Selecting δ small enough that the first term on the right-hand side is less than ε/2 we have
\begin{align}\label{eqn:Xnfluc4}
\sup_{0\leq\theta\leq\delta}\mathbb P
\Bigg\{\!\left\|\int_{\tau_n}^{\tau_n+\theta}\sqrt{n}[F(\pi^n(t))-F(\pi(t))]{\rm d} t\right\|_2>\eta'\Bigg\} \leq \frac{\varepsilon}{2}+\frac{\varepsilon}{2} =\varepsilon.
\end{align}
Therefore, combining (86) and (90) gives
\begin{align*}
\sup_{n\geq
n_0}\sup_{0\leq\theta\leq\delta}\mathbb P\{\|A^{n,c}(\tau_n+\theta)-A^{n,c}(\tau_n)\|_2>\eta\}
\leq\varepsilon
\end{align*}
which shows that condition (A) of Theorem 4 is satisfied for {An,c}n∈ℕ. Therefore, as discussed earlier, {Xn}n∈ℕ is a tight sequence of 𝔻([0, T]: ℓ 2)-valued random variables and, thus, {(Xn, Mn,c)}n∈ℕ is a tight sequence of 𝔻([0, T]: (ℓ 2)2)-valued random variables.
Finally, the ℂ-tightness of {(Xn, Mn,c)}n∈ℕ is immediate from the estimate
\begin{align*}
j_T(X^n) = j_T(M^{n,c}) \leq \frac{2+2k}{\sqrt{n}},\qquad n\in\mathbb N,
\end{align*}
and the fact that there are almost surely no simultaneous jumps, which follows as in the proof of Proposition 3.
5.3. Convergence
In this section we give the proofs of Proposition 2 and Theorem 2. Since we have shown tightness of {(Xn, Mn,c)}n∈ℕ in Section 5.2, all that remains in order to complete the proof of Theorem 2 is to characterize the weak limit points of this sequence of processes. This will be argued by showing that the limit point of any weakly convergent subsequence of {Xn}n∈ℕ will be a solution to SDE (11) and that uniqueness holds for (11) in an appropriate class, which will also prove Proposition 2. We begin by establishing a uniform integrability property for the sequence {Mn,c}n∈ℕ.
Lemma 5
Suppose that {πn}n∈ℕ satisfies the conditions of Proposition 4. Then the sequence 
$\{\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty }|M_j^{n,c}(t)|^2\}_{n\in\mathbb N}$ is uniformly integrable.
Proof. It follows from the Cauchy–Schwarz and Burkholder–Davis–Gundy inequalities that
\begin{align}\label{eqn:mbarnUI1}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\bigg(\sum {_{j = 0}^\infty }|M_j^{n,c}
(t)|^2\bigg)^2 \leq
{\sup _{n\in\mathbb N}}\bigg(\sum_{m=0}^\infty\frac{1}{m^2}\bigg)\sum {_{j = 0}^\infty }
\mathbb E\sup{_{0 \le t \le T}}j^2|M_j^{n,c}(t)|^4\mathbb N
\leq \kappa_1{\sup _{n\in\mathbb N}}\sum {_{j = 0}^\infty } j^2\mathbb E[M_j^{n,c}](T)^2.
\end{align}
Recalling the definition of Mn from (22), for each j, 
$\mathbb E[M_j^{n,c}](T)^2$ can be written as
\begin{align*}
\mathbb E[M_j^{n,c}](T)^2 =\mathbb E\bigg\{\!\sum_{\ell\in\Sigma}\frac{1}{n}\left\langle
e_j,\Delta_\ell \Delta_\ell^Te_j\right\rangle_2
N_\ell\Bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\Bigg)
+ \frac{1}{n}\biggl[D_j\biggl(k{\int_0^T}n\pi^n_j(s){\rm d} s\biggr)+D_{j+1} \biggl(k{\int_0^T}n\pi^n_{j+1}(s){\rm d} s\biggr)\biggr]\Bigg\}^2.
\end{align*}
We now consider the first term in the expectation on the right-hand side of the above equation corresponding to the stream of incoming jobs assigned to queues of length j,
\begin{align*}
\mathbb E\bigg(\,\sum_{\ell\in\Sigma}\frac{1}{n}\left\langle e_j,\Delta_\ell
\Delta_\ell^Te_j\right\rangle_2
N_\ell\bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)\bigg)^2.
\end{align*}
First, consider the diagonal terms in the above sum. It follows from Fubini’s theorem that
\begin{align}\label{eqn:diagTerm1}
\mathbb E\sum_{\ell\in\Sigma}\frac{1}{n^2}\left\langle e_j,\Delta_\ell\Delta_\ell^Te_j
\right\rangle_2^2 N_\ell\bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)^2
\leq \mathbb E\sum_{\ell\in\Sigma}\frac{1}{n^2}\left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2^2
\biggl(\biggl(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\biggr)^2
+\frac{n\lambda}{\binom{n}{L}} {\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\biggr).
\end{align}
The fact that 
$\pi^n_i(s)\in[0,1]$ and (31) imply that, for all i ∈ ℕ0,
\begin{align}\label{eqn:diagTerm2}
\frac{\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty \binom{n\pi^n_i(s)}
{\rho_i(\ell)}{\rm d} s \leq \frac{\lambda
n^L}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\frac{\pi^n_i(s)^{\rho_i(\ell)}}{\rho_i(\ell)\text{!}}{\rm d} s \leq \frac{\lambda
n^L}{\binom{n}{L}}T \leq \kappa_3.
\end{align}
Combining (93), (92), (27), and (32) yields
\begin{align}\label{eqn:diagTerm3}
\mathbb E\sum_{\ell\in\Sigma}\frac{1}{n^2}\left\langle e_j,\Delta_\ell\Delta
_\ell^Te_j\right\rangle_2^2 N_\ell\bigg(\frac{n\lambda}
{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)^2
\leq
\frac{\lambda(\kappa_3+1)}{\binom{n}{L}}\mathbb E{\int_0^T}Z(j,n\pi^n(s)){\rm d} s
\leq
\frac{\lambda(\kappa_3+1)c_Zn^L}{\binom{n}{L}}\mathbb E{\int_0^T}(\pi^n_{j-1}(s)
+\pi^n_j(s)){\rm d} s
\leq \kappa_4\mathbb E{\int_0^T}(\pi^n_{j-1}(s)+\pi^n_j(s)){\rm d} s.
\end{align}
We now analyze the cross terms from the sum. It follows from the independence of N ℓ and N ℓ′that
\begin{align}\label{eqn:crossTerm1}
\mathbb E\sum_{\ell,\ell'\in\Sigma\atop \ell\neq\ell'}\frac{1}{n^2} \left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\left\langle
e_j,\Delta_{\ell'}\Delta_{\ell'}^Te_j\right\rangle_2
N_\ell\bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)
\times N_{\ell'}\bigg(\frac{n\lambda}
{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell')}{\rm d} s\bigg)
= \mathbb E\sum_{\ell,\ell'\in\Sigma\atop
\ell\neq\ell'}\frac{1}{n^2}\left\langle e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\left\langle
e_j,\Delta_{\ell'}\Delta_{\ell'}^Te_j\right\rangle_2\bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)
\times \bigg(\frac{n\lambda}
{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell')}{\rm d} s\bigg).
\end{align}
Combining (95), (27), and (32) yields
\begin{align}\label{eqn:crossTerm2}
\mathbb E\sum_{\ell,\ell'\in\Sigma\atop \ell\neq\ell'}\frac{1}{n^2}\left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\left\langle
e_j,\Delta_{\ell'}\Delta_{\ell'}^Te_j\right\rangle_2N_\ell\bigg(\frac{n\lambda}{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell)}{\rm d} s\bigg)
\times N_{\ell'}\bigg(\frac{n\lambda}
{\binom{n}{L}}{\int_0^T}\prod_{i=0}^\infty
\binom{n\pi^n_i(s)}{\rho_i(\ell')}{\rm d} s\bigg)
\leq
\mathbb E\frac{\lambda^2}{\big(\binom{n}{L}\big)^2}{\int_0^T}Z(j,n\pi^n(s))^2{\rm d} s
\leq\mathbb E
\frac{\lambda^2n^{2L}}{\big(\binom{n}{L}\big)^2}{\int_0^T}(\pi^n_{j-1}(s)+
\pi^n_j(s))^2{\rm d} s
\leq \kappa_5\mathbb E{\int_0^T}(\pi^n_{j-1}(s)+\pi^n_j(s)){\rm d} s.
\end{align}
Similarly,
\begin{align*}
\mathbb E\bigg(\frac{1}{n}\bigg[D_j\bigg(k{\int_0^T}n\pi^n_j(s){\rm d} s\bigg)+D_{j+1} \bigg(k{\int_0^T}n\pi^n_{j+1}(s){\rm d} s\bigg)\bigg]\bigg)^2
\leq \kappa_6\mathbb E{\int_0^T}[\pi^n_j(s)+\pi^n_{j+1}(s)]{\rm d} s.
\end{align*}
Combining this estimate with (94), (96), and (50) gives
\begin{align*}
{\sup _{n\in\mathbb N}}\sum {_{j = 0}^\infty } j^2{\rm E}[M_j^{n,c}](T)^2
\leq\kappa_7{\sup _{n\in\mathbb N}}\mathbb E{\int_0^T}\sum_{j=1}^\infty
j^2(\pi_{j-1}^n(s)+\pi_{j}^n(s)+\pi_{j+1}^n(s)){\rm d} s<\infty,
\end{align*}
which, in view of (91), gives the desired uniform integrability.
The following lemma together with (82) shows that any weak limit point X of {Xn}n∈ℕ satisfies X(t) ∈ ℓ̃ 2 for all t ∈ [0, T] almost surely.
Lemma 6
Let zn and z be 𝔻([0, T]: ℓ 2)-valued random variables such that
\begin{align*}
\sup{_{0 \le t \le T}}\|z^n(t)-z(t)\|_2\to0 \quad{in\,\; probability\,\; as}\,\;n\to\infty.
\end{align*}
Suppose that 
${\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(z_j^n(t))^2<\infty$. Then 
$\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(z_j(t))^2< \infty$ almost surely and 
$\sup{_{0 \le t \le T}}|\sum {_{j = 0}^\infty } z_j^n(t)-\sum {_{j = 0}^\infty } z_j(t)|\to0$ in probability.
Proof. Let 
$\kappa = {\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2[z^n_j(t)]^2$. Note that
\begin{align*}
{\sup _{n\in\mathbb N}}\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } |z_j^n(t)|
\leq \bigg(\sum_{j=1}^\infty\frac{1}{j^2}\bigg)^{1/2}\sqrt{\kappa}
<\infty.
\end{align*}
Also, by Fatou’s lemma, 
$\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(z_j(t))^2\leq\kappa$ and so we have 
$\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } |z_j(t)| < \infty$ almost surely as well. Now
\begin{align*}
\mathbb E\bigg[\sup{_{0 \le t \le T}}|\sum {_{j = 0}^\infty } z^n_j(t)-
\sum {_{j = 0}^\infty } z_j(t)|\wedge 1\bigg]
\leq\mathbb E\bigg[\sup{_{0 \le t \le T}}|\sum_{j=0}^m
z^n_j(t)-\sum_{j=0}^m z_j(t)|\wedge 1\bigg]+ \mathbb E\bigg[\sup_{0\leq
t\leq T}|\sum_{j=m+1}^\infty z^n_j(t)|\wedge
1\bigg]
+\mathbb E\bigg[\sup{_{0 \le t \le T}}|\sum_{j=m+1}^\infty
z_j(t)|\wedge 1\bigg]
\equiv T_1^m(n)+T_2^m(n)+T_3^m(n).
\end{align*}
Then, for κ 1 ∈ (0, ∞),
\begin{align*}
(T_2^m(n))^2\leq \bigg(\sum_{j=m+1}^\infty\frac{1}{j^2}\bigg)\kappa_1
\quad\text{and}\quad (T_3^m(n))^2\leq \bigg(\sum_{j=m+1}^\infty
\frac{1}{j^2}\bigg)\kappa_1.
\end{align*}
The result now follows on first sending n → ∞ and then m → ∞.
The following result that shows that Φ(t) is a trace class operator will be useful in characterizing the martingale term in the limiting diffusion. Note that, from definition (14), Φ(t) is a nonnegative operator.
Lemma 7
For each t ∈ [0, T], Φ(t) is a nonnegative trace class operator. Denote by a(t) the nonnegative square root of Φ(t). Then 
${\int_0^T} {\left\| {a\left(s \right)} \right\|} _{{\rm{HS}}}^2{\rm{d}}s < \infty$.
Proof. We first show that Φ(t) is a trace class operator. Since Φ(t) is nonnegative (and hence self-adjoint), it suffices to show that
\begin{align*}
\sum {_{j = 0}^\infty }\left\langle e_j,\Phi(s)e_j\right\rangle_2<\infty.
\end{align*}
Using an argument similar to that used in the derivation of (30), we can write 〈e j, Φ(s)e j〉2 as
\begin{align}\label{eqn:limcrossquad}
\left\langle e_j,\Phi(s)e_j\right\rangle_2 = \lambda
L\text{!}\,\bar{Z}(j,\pi(s))+k(\pi_j(s)+\pi_{j+1}(s)),
\end{align}
where the definition of Z̄ is analogous to Z, given as
\begin{align}
\bar Z(j,\pi (s))\cr
\quad \quad {\kern 1pt} : = {\kern 1pt} \sum\limits_{{i_1} = 0}^{k - 2} {{{{(\sum\limits_{m = 0}^{j - 1} {\pi _m}(s))}^{{i_1}}}} \over {{i_1}{\rm{!}}}}\sum\limits_{{i_2} = 0}^{L - {i_1}} {{{\pi _{j - 1}}{{(s)}^{{i_2}}}} \over {{i_2}{\rm{!}}}}\sum\limits_{{i_3} = 0}^{L - {i_1} - {i_2}} {[{i_2} \wedge {(k - {i_1})_ + } - {i_3} \wedge {(k - {i_1} - {i_2})_ + }]^2}\cr
\quad \times {{{\pi _j}{{(s)}^{{i_3}}}} \over {{i_3}{\rm{!}}}}{{{{(\sum\limits_{m = j + 1}^\infty {\pi _m}(s))}^{L - {i_1} - {i_2} - {i_3}}}} \over {(L - {i_1} - {i_2} - {i_3}){\rm{!}}}}.
\end{align}
For completeness, the proof of (97) is provided in Appendix A.5. Using similar arguments as in (32) and (61), it is easy to see that there exists c Z̄ ∈ (0, ∞) such that, for all j ∈ ℕ0,
\begin{align}\label{eqn:zbarDiag}
\bar{Z}(j,\pi(s))
\leq c_{\bar{Z}} (\pi_{j-1}(s)+\pi_j(s)).
\end{align}
Once again, details on this step are given in Appendix A.5. From (97) and (99), it follows that there exists a κ 1 ∈ (0, M) such that
\begin{align*}
\sum {_{j = 0}^\infty }\left\langle e_j,\Phi(t) e_j\right\rangle_2
\leq \kappa_2\sum {_{j = 0}^\infty }[\pi_{j-1}(t)+\pi_{j}(t)+\pi_{j+1}(t)]\leq 3\kappa_1.
\end{align*}
Therefore, Φ(t) is a trace class operator. Finally, note that
\begin{align*}
{\int_0^T}\|a(s)\|_{\rm HS}^2{\rm d} s = {\int_0^T}\sum {_{j = 0}^\infty }\left\langle
a(s)e_j,a(s)e_j\right\rangle_2 {\rm d} s = {\int_0^T}\sum {_{j = 0}^\infty }\left\langle
e_j,\Phi(s)e_j\right\rangle_2 {\rm d} s\leq 3\kappa_1T,
\end{align*}
which completes the proof.
We now proceed with the proofs of Proposition 2 and Theorem 2.
Proof of Proposition 2. The existence of an (X(t))0≤t≤T as in the statement of Proposition 2 will be proved as part of Theorem 2. We now consider the second statement in Proposition 2, and let (X(t))0≤t≤T and (X̃(t))0≤t≤T be two 
$\{\mathcal F_t\}$-adapted processes solving (12) with sample paths in ℂ([0, T]: ℓ 2) such that X(t) ∈ ℓ̃ 2 and X̃(t) ∈ ℓ̃ 2 for all t almost surely. In order to show that X(t) = X̃(t) for all t ∈ [0, T] almost surely, it suffices to show the following Lipschitz property on G. There exists a C ∈ (0, ∞) such that, for all x, x̃ ∈ ℓ̃ 2,
\begin{align}\label{eqn:Glip}
\sup{_{0 \le t \le T}}\|G(x,\pi(t))-G(\tilde x,\pi(t))\|_2\leq C\|x-\tilde x\|_2.
\end{align}
Note that from (4), (5), and (16), for j ∈ ℕ0 and (x, r) ∈ ℓ̃ 2 × 𝒮,
\begin{align}\label{eqn:Gtoxi}
G_j(x,r) = \lambda L\text{!}\, [\xi^{\bf 1}_{j-1}(x,r)-\xi^{\bf
1}_{j}(x,r)+\xi^2_{j-1}(x,r)-\xi^2_{j}(x,r)+\xi^3_{j-1}(x,r)-
\xi^3_{j}(x,r)]+k\xi^4_{j}(x),
\end{align}
where
\begin{align*}
\xi^{\bf 1}_{j}(x,r)\,{:{=}}\, \sum_{i_1=0}^{k-1}i_1
{(\!\sum_{m=0}^{j-1}r_m)^{i_1-1}}{i_1{!}}
\sum\limits_{{i_2} = 1}^{L - {i_1}} {\left[ {{i_2} \wedge \left({k - {i_1}} \right)} \right]}
{(r_j)^{i_2}}{i_2\text{!}}
{(\!\sum_{m=j+1}^{\infty}r_{m})^{L-i_1-i_2}}{(L-i_1-i_2)\text{!}}
\sum_{m=0}^{j-1}x_m,
\xi^2_{j}(x,r)\,{:{=}}\, \sum_{i_1=0}^{k-1}{(\!\sum_{m=0}^{j-1}
r_m)^{i_1}}{i_1\text{!}}\sum_{i_2=1}^{L-i_1}i_2[i_2\wedge(k-i_1)]
{(r_j)^{i_2-1}}{i_2\text{!}}{(\!\sum_{m=j+1}^{\infty}
r_{m})^{L-i_1-i_2}}{(L-i_1-i_2)\text{!}}x_j,
\xi^3_{j}(x,r)\,{:{=}}\,\sum_{i_1=0}^{k-1}{(\!\sum_{m=0}^{j-1}
r_m)^{i_1}}{i_1\text{!}}\sum_{i_2=1}^{L-i_1}(L-i_1-i_2)[i_2
\wedge(k-i_1)]{(r_j)^{i_2}}{i_2\text{!}}{\big(\!\sum_{m=j+1}^{\infty}
r_{m}\big)^{L-i_1-i_2-1}}{(L-i_1-i_2)\text{!}}
\times\sum_{m=j+1}^{\infty}x_m,
\end{align*}
and
\begin{align*}
\xi^4_{j}(x)= [x_{j+1}-x_{j}].
\end{align*}
Also, let 
$\xi^i\,{:{=}}\,(\xi^i_{j})_{j=0}^\infty$ for i = 1, 2, 3, 4. Using the triangle inequality, it suffices to show that (100) holds with G replaced with ξi, i = 1, 2, 3, 4. Since π(t) ∈ 𝒮 for all t ∈ [0, T],
\begin{align}\label{eqn:xi1Bound}
\sup{_{0 \le t \le T}}\|\xi^1(x,\pi(t))-\xi^1(\tilde x,\pi(t))\|^2_2 \leq
\kappa_1'\sup_{0\leq t\leq
T}\sum {_{j = 0}^\infty }\pi_j(t)^2\bigg[\sum_{m=0}^{j-1}x_m-\sum_{m=0}^{j-1}\tilde
x_m\bigg]^2
\leq \kappa_1'\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j\pi_j(t)\|x-\tilde
x\|_2^2
\leq \kappa_1\|x-\tilde x\|_2^2,
\end{align}
where the last inequality is from (79). Also,
\begin{align*}
\sup{_{0 \le t \le T}}\|\xi^2(x,\pi(t))-\xi^2(\tilde x,\pi(t))\|^2_2 \leq
\kappa_2\sum {_{j = 0}^\infty }[x_j-\tilde x_j]^2 = \kappa_2\|x-\tilde x\|_2^2.
\end{align*}
Using the fact that 
$\sum {_{m = 0}^\infty } {x_m} = \sum {_{m = 0}^\infty {{\tilde x}_m} = 0}$ and the calculation in (102),
\begin{align*}
\sup{_{0 \le t \le T}}\|\xi^3(x,\pi(t))-\xi^3(\tilde x,\pi(t))\|^2_2 \leq
\kappa_3'\sup_{0\leq t\leq
T}\sum {_{j = 0}^\infty }\pi_j(t)^2\bigg[\sum_{m=j+1}^{\infty}x_m-\sum_{m=j+1}^{\infty}\tilde
x_m\bigg]^2
= \kappa_3'\sup_{0\leq t\leq
T}\sum {_{j = 0}^\infty }\pi_j(t)^2\bigg[\sum_{m=0}^{j}\tilde
x_m-\sum_{m=0}^{j}x_m\bigg]^2
\leq \kappa_3\|x-\tilde x\|_2^2.
\end{align*}
Finally,
\begin{align*}
\|\xi^4(x)-\xi^4(\tilde x)\|^2_2 \leq \sum {_{j = 0}^\infty }[x_j-\tilde
x_j]^2+\sum {_{j = 0}^\infty }[x_{j+1}-\tilde x_{j+1}]^2 \leq 2\|x-\tilde x\|_2^2.
\end{align*}
Combining the above Lipschitz estimates for ξi, i = 1, 2, 3, 4, we have (100) and the result follows.
We now proceed to the proof of Theorem 2.
Proof of Theorem 2. From Proposition 4, {(Xn, Mn,c)}n∈ℕ is ℂ-tight in 𝔻([0, T]: (ℓ 2)2). Suppose that (X, Mc) is a weak limit of a subsequence of {(Xn, Mn,c)}n∈ℕ (also indexed by {n}) given on some probability space (Ω 
$\mathcal F$, ℙ). Let m ∈ ℕ, and let 
$\mathcal H\colon(\ell_2\times\ell_2)^m\to\mathbb R$ be a bounded and continuous function. For s ≤ t ≤ T and 0 ≤ t 1 ≤ · · · ≤ t m ≤ s, we let 
$\xi_i^n=(X^n(t_i),M^{n,c}(t_i))$ and ξ i = (X(t i), Mc(t i)). Then, for all j ∈ ℕ0,
\begin{align*}
\mathbb E\mathcal H(\xi_1,\ldots,\xi_m)[M^c_j(t)-M^c_j(s)] =
\lim_{n\to\infty}\mathbb E\mathcal H(\xi_1^n,\ldots,\xi_m^n)[M^{n,c}_j(t)-M_j^{n,c}
(s)]=0,
\end{align*}
where the first equality follows from the uniform integrability property proved in Lemma 5 and the second equality follows from the fact that Mn,c is a martingale for each n ∈ ℕ. It follows that Mc is an 
$\{\mathcal F_t\}$-martingale, where 
$\mathcal F_t=\sigma\{X(s),M^c(s),\, s\leq t\}$.
As was shown in (60),
\begin{align}\label{eqn:crossquad2}
\left\langle M_i^{n,c},M_j^{n,c}\right\rangle(t) = n\left\langle M_i^n,M_j^n\right\rangle(t)
= \frac{\lambda}{\binom{n}{L}}\int_0^tZ(i,j,n\pi^n(s)){\rm d} s-
k\int_0^t{\bf 1}_{\{i = j+1\}}\pi_i^n(s){\rm d} s
- k\int_0^t{\bf 1}_{\{i+1 = j\}}\pi_j^n(s){\rm d} s+k\int_0^t{\bf
1}_{\{i=j\}}(\pi_j^n(s)+\pi_{j+1}^n(s)){\rm d} s
\end{align}
(see (30) and (58) for the definition of Z). Using similar arguments as in (58), we have the estimate
\begin{align*}
\left\langle e_i,\Phi(s)e_j\right\rangle_2 = \lambda L\text{!}\, \bar{Z}(i,j,\pi(s))-k{\bf
1}_{\{i+1=j\}}\pi_j(s)-k{\bf 1}_{\{i=j+1\}}\pi_i(s)+k{\bf
1}_{\{i=j\}}(\pi_j(s)+\pi_{j+1}(s)),
\end{align*}
where, for i < j,
\begin{align*}
\bar{Z}(i,j,\pi(s)) \,{:{=}}\, \sum_{i_1=0}^{k-2}
\frac{(\sum_{m=0}^{i-2}\pi_m(s))^{i_1}}{i_1\text{!}}
\sum_{i_2=0}^{k-i_1-1}{\pi_{i-1}(s)^{i_2}}
{i_2{!}}\sum_{i_3=0}^{k-i_1-i_2-1}[i_2-i_3]{\pi_{i}(s)^{i_3}}{i_3{!}}
\times\sum_{i_4=0}^{k-i_1-i_2-i_3-1}{(\sum_{m=i+1}^{j-2}
\pi_m(s))^{i_4}}{i_4\text{!}}\sum_{i_5=0}^{L-\sum_{n=1}^4i_n}
\frac{\pi_{j-1}(s)^{i_5}{\bf 1}_{\{j> i+1\}}}{i_5\text{!}}
\times\sum_{i_6=0}^{L-\sum_{n=1}^5i_n}\!\biggl[\!({\bf 1}_{\{j=
i+1\}}(i_3-i_5)+i_5)\wedge\!\biggl(k-\sum_{n=1}^4i_n\biggr)_+\! -i_6\wedge\!\biggl(k-\sum_{n=1}^5i_n\biggr)_+\biggr]
\times\frac{\pi_{j}(s)^{i_6}}{i_6\text{!}}\frac{(\sum_{m=j+1}^\infty
\pi_m(s))^{L-\sum_{n=1}^6i_n}}{ (L-\sum_{n=1}^6i_n)\text{!}}
\end{align*}
for i > j, Z̄(i, j, π(s)) := Z̄(j, i, π(s)), and, for i = j, Z̄(j, j, π(s)) := Z̄(j, π(s)), where Z̄(j, r) is defined in (98). Using arguments similar to those used in (46) and (47), we can write
\begin{align*}
\left|Z(i,j,n\pi^n(s)) -
\frac{n\text{!}}{(n-L)\text{!}}\bar{Z}(i,j,\pi^n(s))\right|\leq \kappa_1n^{L-1}.
\end{align*}
From this, (97), (103), and the fact that πn → π in probability, it follows that
\begin{align*}
\sup{_{0 \le t \le T}}\left|\left\langle M_i^{n,c}(t),M_j^{n,c}(t)\right\rangle-
\int_0^t\left\langle e_i,\Phi(s)e_j\right\rangle_2{\rm d} s\right|\to 0
\end{align*}
in probability. A similar argument as in Lemma 5 shows that 
$\left\{ {{{\left\langle {M_i^{n,c},M_j^{n,c}} \right\rangle }_t}} \right\}_{n \in\mathbb N}$ is uniformly integrable for each t ∈ [0, T] and i, j ∈ ℕ0. Applying the above convergence and uniform integrability properties,
\begin{align*}
\mathbb E\mathcal H(\xi_1,\ldots,\xi_m)[\left\langle M^c_i,M^c_j\right\rangle_t-\left\langle M^c_i,M^c_j
\right\rangle_s-\int_s^t\left\langle e_i,\Phi(u)e_j\right\rangle_2{\rm d} u]
= \lim_{n\to\infty}\mathbb E\mathcal H(\xi_1^n,\ldots,\xi_m^n)\Big[\left\langle
M_i^{n,c},M_j^{n,c}\right\rangle_t-\left\langle M^{n,c}_i,M^{n,c}_j\right\rangle_s-\int_s^t\left\langle
e_i,\Phi(u)e_j\right\rangle_2{\rm d} u\Big]
=0.
\end{align*}
Also, from Lemma 5 and Fatou’s lemma, 
$\mathbb E\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty }|M^c_j(t)|^2<\infty$. Thus, 
$M^c\,{:{=}}\, (M^c_j)_{j\in\mathbb N_0}$ is a collection of square-integrable 
$\{\mathcal F_t\}$-martingales with
\begin{align*}
\left\langle M^c_i,M^c_j\right\rangle (t) = \int_0^t\left\langle e_i,\Phi(s)e_j\right\rangle_2{\rm d} s,\qquad
t\in[0,T].
\end{align*}
From Theorem 8.2 of [Reference Da Prato and Zabczyk10], it now follows that there is an ℓ 2-cylindrical Brownian motion {(W t(h))0≤t≤T : h ∈ ℓ 2} on some extension 
$(\bar{\Omega},\bar{\mathcal F},\bar{\mathbb P},\{\bar{\mathcal F}_t\})$ of the filtered probability space (Ω, 
$\mathcal F$, ℙ, $\{\mathcal F_t\}$) such that
\begin{align}\label{eqn:limitMartingale}
M^c(t)=\int_0^ta(s){\rm d} W(s).
\end{align}
Recall the representation of Xn in terms of An,c and Mn,c from (23). We now argue that, together with Xn and Mn,c, An,c(·) converges to 
$\int_0^ \cdot G (X(s),\pi (s)){\rm{d}}s$ in 𝔻([0, T]: ℓ 2) in distribution as n → ∞ (along the chosen subsequence). The definition of An,c in (24) and the estimate in (48) imply that
\begin{align}\label{eqn:Aexp}
\sup{_{0 \le t \le T}}\bigg\|A^{n,c}(t) -
\int_0^t\sqrt{n}[F(\pi^n(s))-F(\pi(s))]{\rm d} s\bigg\|_2\leq\frac{\kappa_2}{\sqrt{n}\,}.
\end{align}
For r, r̃ ∈ 𝒮 such that (r – r̃) ∈ ℓ̃ 2, the ith component of F(r) − F(r̃) can be written as
\begin{align*}
{[}F(r)-F(\tilde r){]}_i = {\int_0^1}\frac{\partial}{\partial u}F_i(r u +
(1-u)\tilde r){\rm d} u\\ = {\int_0^1}G_i((r-\tilde r), ru+(1-u)\tilde r){\rm d} u\\ =
G_i(r - \tilde r, \tilde r)+{\int_0^1}[G_i((r-\tilde r, ru+(1-u)\tilde r)-G_i(r-\tilde r,\tilde r)]{\rm d} u.
\end{align*}
Therefore, observing that cG i(x, r) = G i(cx, r) for c ∈ ℝ and (x, r) ∈ ℓ̃ 2 × 𝒮, and noting from (82) that Xn(s) ∈ ℓ̃ 2 for every s ∈ [0, T] almost surely, we can write
\begin{align}\label{eqn:taylor}
\sqrt{n}[F(\pi^n(s))-F(\pi(s))]_i = G_i(X^n(s),\pi(s))+R_i^n(s),
\end{align}
where
\begin{align*}
R^n_i(s) = {\int_0^1}[G_i(X^n(s),\pi^n(s)u+(1-u)\pi(s))-G_i(X^n(s),\pi(s))]{\rm d} u.
\end{align*}
Thus,
\begin{align*}
\sqrt{n}[F(\pi^n(s))-F(\pi(s))] = G(X^n(s),\pi(s))+R^n(s),
\end{align*}
where 
$R^n(s)\,{:{=}}\, (R^n_i(s))_{i\in\mathbb N_0}$. We now show that 
${{\int_0^T} {\left\| {{R^n}\left(s \right)} \right\|} _2}{\rm{d}}s \to 0$ in probability as n → ∞. Since 
$\sum {_{m = 0}^j} X_m^n(s) = - \sum {_{m = j + 1}^\infty } X_m^n(s)$, it follows from (36) that, for r, r̃ ∈ 𝒮,
\begin{align*}
\|\xi^i(X^n(s),r)-\xi^i(X^n(s),\tilde r)\|_2^2
\leq \kappa_3'\sum {_{j = 0}^\infty }\bigg(\sum_{m=0}^{j}|X^n_m(s)
|\bigg)^2\bigg[[r-\tilde r]_j^2+\tilde r_j\bigg(\sum_{i=0}^{j-1}[r-\tilde r]_i\bigg)^2+\tilde r_j\bigg(\sum_{i=j+1}^{\infty}[r-\tilde r]_i\bigg)^2\bigg]
\leq \kappa_3\bigg(\sum {_{j = 0}^\infty }
j^2|X^n_j(s)|^2\bigg)\sum {_{j = 0}^\infty }[j\tilde r_j\|r-\tilde r\|_2^2+[r-\tilde r]_j^2]
\end{align*}
for i = 1, 2, 3. The triangle inequality, (101), and the observation that 
${\sup _{0 \le s \le T}}\sum {_{j = 0}^\infty j{\pi _j}\left(s \right) < \infty }$ (see (79)) then imply that
\begin{align*}
\|G(X^n(s),\pi^n(s)u+(1-u)\pi(s))-G(X^n(s),\pi(s))\|_2^2\leq
\kappa_3\bigg(\sum {_{j = 0}^\infty } j^2|X^n_j(s)|^2\bigg)\|
\pi^n(s)-\pi(s)\|_2^2.
\end{align*}
Since sup0≤s≤T ||πn(s) − π(s)||2 → 0 in probability and, from (82), 
${\sup _{n\in\mathbb N}}\mathbb E\sup_{0\leq s\leq T} \sum {_{j = 0}^\infty } j^2|X^n_j(s)|^2<\infty$, it follows that
\begin{align*}
\sup_{0\leq u\leq 1}\sup_{0\leq s\leq T}\|G(X^n(s),\pi^n(s)u+(1-u)\pi(s))-G(X^n(s),\pi(s))\|_2\to 0
\end{align*}
in probability as n → ∞, and, thus,
\begin{align}\label{eq:Eq504}
{\int_0^T}\|R_n(s)\|_2{\rm d} s \to 0 \quad\text{in probability.}
\end{align}
In view of (105), (106), and (107), it now suffices to show that, along the subsequence,
\begin{align*}
({X^n},{M^{n,c}},\int_0^ \cdot G({X^n}(s),\pi (s)){\kern 1pt} s) \Rightarrow (X,\bar M,\int_0^ \cdot G(X(s),\pi (s)){\kern 1pt} s)
\end{align*}
in 𝔻([0, T]: (ℓ 2)3). By appealing to the Skorokhod representation theorem we can assume without loss of generality that (Xn, Mn,c) converges almost surely in 𝔻([0, T]: (ℓ 2)2) to (X, M̄). From (82) and Fatou’s lemma we also have
\begin{align*}
\sup{_{0 \le t \le T}}\sum {_{j = 0}^\infty } j^2(X_j(t))^2<\infty \quad\text{almost surely.}
\end{align*}
Also, since 
$\sum\nolimits_{j = 0}^\infty X_j^n(t) = 0$ for all t ∈ [0, T] and n ∈ ℕ, by Lemma 6 and (82), we have 
${\sum {_{j = 0}^\infty }} X_j(t)=0$ for all t ∈ [0, T] almost surely as well. It then follows that Xn(t), Xn(t) ∈ ℓ̃ 2 for all t ∈ [0, T] almost surely for all n ∈ ℕ. From the Lipschitz property in (100), it now follows that, as n → ∞,
\begin{align*}%\label{eqn:Gconv}
{\int_0^T}\|G(X^n(s),\pi(s))-G(X(s),\pi(s))\|_2{\rm d} s \leq
C{\int_0^T}\|X^n(s)-X(s)\|_2{\rm d} s \to 0,
\end{align*}
which proves the desired convergence. Together with (23) and representation (104), it follows that the limit point (X, Mc) satisfies
\begin{align*}
X(t) = x_0 + \int_0^tG(X(s),\pi(s)){\rm d} s + \int_0^ta(s){\rm d} W(s)
\end{align*}
almost surely for all t ∈ [0, T]. Since Xn(t) ∈ ℓ̃ 2 for all t ∈ [0, T] almost surely, this in particular proves the existence part of Proposition 2. Finally, the uniqueness part of Proposition 2 (which was established earlier in this section) now says that Xn converges in distribution along the full sequence to the unique weak solution of (11) with values in ℓ̃ 2. The result follows.
6. Numerical results
In this section we present some simulation results comparing the prelimit n-server system with results of the corresponding law of large number and central limit approximations. We consider a system with n = 10 000 servers. For all combinations of L and k in the set {(L, k) ∈ ℕ × ℕ: 2 ≤ L ≤ 5, k < L}, we simulate 1000 realizations of both the n-server system and the diffusion approximation given in Theorem 2 using parameters T = 10, λ = 0.9, and c = 1. Note that, since the limiting processes are infinite-dimensional, we must truncate to a finite-dimensional approximation in order to perform simulations. In our numerical approximations, we truncate to the first 20 coordinates. All computations were performed in MATLAB®. A numerical ODE solver (ode45) was used to compute the ODE corresponding to the law of large number limit. The limit diffusion was simulated using Euler’s method with step sizes of 0.1. The realizations of the diffusion were used to create 95% confidence intervals for the following metrics at time T; the number of empty queues, the number of ‘large’ queues (queues with more than 5 jobs), and the mean queue length. The coverage rates (i.e. the proportion of the n-server system simulations which fall within the 95% confidence interval estimated by the diffusion approximation) can be found in Tables 1, 2, and 3. The diffusion approximation based confidence intervals for the first and third cases contain approximately 95% of the n-server simulated observations, as desired. However, the results given in Table 2 appear to be less satisfactory. This is not surprising since the events corresponding to large queues are rare, and, thus, their probabilities are harder to estimate.
TABLE 1: Empty queue coverage rate.
TABLE 2: Large queue coverage rate.
TABLE 3: Mean queue length coverage rate.
In Figure 1 we present two graphs showing that the diffusion approximation captures the fluctuations of the underlying processes around the limiting ODE given via the LLN. We consider the supermarket model, i.e. (L, k) = (2, 1) with λ = 0.9, T = 50, c = 1, and n = 10 000 servers. Figures 1(a) and 1(b) present the results for large queues (i.e. queues of length at least 5) and empty queues, respectively. In each, the solid line represents the numerical solution to the limiting ODE, the dashed line represents the a single simulation of the underlying system, and the dotted lines represent empirical 95% confidence intervals obtained from the diffusion approximation. Namely, 1000 realizations of the diffusion were computed and the 2.5 and 97.5 percentiles were taken at each time point. The figures show that both the LLN and diffusion approximation are doing a good job of approximating the dynamics of the finite system over time.
FIGURE 1: Empirical 95% confidence intervals obtained from the diffusion approximation, effectively capturing the fluctuations around the LLN ODE
The goal of this paper was to develop reliable approximations of the n-server system that are much quicker to simulate. In Table 4 we present the average time (in seconds) required to simulate one trial of the finite system and diffusion approximation. As is seen from these tables, the time required to simulate the diffusion approximations is substantially smaller than for the underlying n-server jump-Markov process. In addition, increasing n will further increase the amount of time required to simulate the n-server system. Indeed, n = 10 000 is a small number compared to the size of typical data centers and server farms that have machines which number in the hundreds of thousands. The point at which it becomes quicker to use the diffusion approximation will depend heavily on the system parameters. Indeed, simulation results indicate that this point occurs in the mid-hundreds for (L, k) = (2, 1) while it is in the mid-thousands for (L, k) = (5, 4). A further caveat is that this number will depend on the efficiency of implementations of both the numerical approximation and the simulation scheme.
TABLE 4: Average simulation times for (a) the finite system and (b) the limit diffusion.
Appendix A. Auxiliary results
A.1. Criterion for tightness of Hilbert-valued random variables
The following theorem gives sufficient conditions for tightness of a sequence of random variables taking values in a (possibly infinite-dimensional) Hilbert space. For a proof see Corollary 2.3.1 of [Reference Kallianpur and Xiong21].
Theorem 3
Let ℍ be a separable Hilbert space with inner product 〈·, ·〉 and complete orthonormal system 
$\{e_i\}_{i=1}^\infty$. Suppose that {Y n}n∈ℕ is a sequence of ℍ-valued random variables satisfying the following conditions:
(a) for each n 0 ∈ ℕ, limA→∞ supn∈ℕ ℙ(max1≤i≤n 0〈Y n, e i〉2 > A) = 0;
(b) for every δ > 0,
${\lim _{n0 \to \infty }}{\sup _{n \in {\rm{N}}}}{\rm{P}}\left({\sum {_{j = n0}^\infty } {{\left\langle {Y_{n,}}{e_j}\right\rangle }^2} > \delta } \right) = 0$.
Then {Y n}n∈ℕ is a tight sequence of ℍ-valued random variables.
A.2. Criterion for tightness of RCLL processes
The following theorem gives a criterion for tightness of a sequence of RCLL processes with values in a Polish space; see [Reference Kurtz26].
Theorem 4
Let 𝕊 be a Polish space and let {Y n}n∈ℕ be a sequence of 𝔻([0, T]: 𝕊)-valued 
$\{\mathcal F^n_t\}$-semimartingales satisfying the following conditions:
(T1) {Y n(t)}n∈ℕ is tight for every t in a dense subset of [0, T];
(A) for each ε > 0, η > 0, and T 0 ∈ [0, T − ε], there exists a δ > 0 and n 0 with the property that, for every collection of stopping times (τ n)n∈ℕ (τn being an
$\mathcal F^n_t\,{:\!=}\,\sigma\{Y_n(s)\colon s\leq t\}$-stopping time) with τ n ≤ T 0,
where d(·, ·) is the distance on 𝕊.\begin{align*}
\sup_{n\geq n_0}\sup_{0\leq\theta\leq\delta}\mathbb P\{d(Y_n(\tau_n+\theta),Y_n(\tau_n))\geq\eta\}\leq \varepsilon,
\end{align*}
Then {Y n}n∈ℕ is tight in 𝔻([0, T]: 𝕊).
A.3. Hilbert–Schmidt and trace class operators
Here we collect some elementary facts about trace class and Hilbert–Schmidt operators. We refer the reader to [Reference Reed and Simon34] for details. For a separable Hilbert space ℍ (with inner product 〈·, ·〉 and norm || · ||), let 
$\mathcal L(\mathbb H)$ be the collection of all bounded linear operators on ℍ. An operator 
$A\in\mathcal L(\mathbb H)$ is called nonnegative if 〈u, Au〉≥ 0 for all u ∈ ℍ. Such an operator is called trace class if, for a complete orthonormal system (CONS) {e i} in ℍ, Σi〈Ae i, e i〉 < ∞ in which case the quantity is finite (and is the same) for every CONS {e i}. An operator $A\in\mathcal L(\mathbb H)$ is called Hilbert–Schmidt if there exists a CONS {e i} in ℍ such that Σj〈Ae j, Ae j〉 = Σj ||Ae j||2 < ∞. In that case, this quantity is the same for all CONS {e i} and its square root is called the Hilbert–Schmidt norm of A, denoted as ||A||HS. For a nonnegative operator $A\in\mathcal L(\mathbb H)$, there is a unique nonnegative 
$B\in\mathcal L(\mathbb H)$ referred to as the nonnegative square root of A such that B 2 = A. If A is a trace class operator then B is a Hilbert–Schmidt operator.
A.4. Cylindrical Brownian motion
A collection of continuous real stochastic processes {(W t(h))0≤t≤T : h ∈ ℓ 2} given on a filtered probability space (Ω, 
$\mathcal F$, ℙ, 
$\{\mathcal F_t\}$) is called a ℓ 2-cylindrical Brownian motion if, for every h ∈ ℓ 2, (W t(h))0≤t≤T is a 
$\{\mathcal F_t\}$-Brownian motion with variance 
$\|h\|^2_2$ and, for h, k ∈ ℓ 2,
\begin{align*}
\left\langle W(h),W(k)\right\rangle_t = \left\langle h,k\right\rangle_2t,\qquad 0\leq t\leq T.
\end{align*}
For a measurable map a from [0, T] to the space of Hilbert–Schmidt operators from ℓ 2 to ℓ 2 such that 
${\int_0^T} {a(} s)_{{\rm{HS}}}^2{\rm{d}}s < \infty $, we denote by 
$\int_0^t a (s)dW(s)$ the ℓ 2-valued martingale defined as the limit of
\begin{align*}
\sum_{i=1}^n\sum_{j=1}^n\phi_i\int_0^t\left\langle\phi_i,a(s)\phi_j\right\rangle_2{\rm d}
W_s(\phi_j) \quad{{\rm as}\,\; n\to\infty,}
\end{align*}
where {φ i}i∈ℕ is a CONS in ℓ 2. We refer the reader to Chapter 4 of [Reference Da Prato and Zabczyk10] regarding the fact that the limit exists and is independent of the choice of the CONS.
A.5. Proofs of (97) and (99)
Recalling the definition of Φ(s) in (14) we can write
\begin{align}\label{eqn:incomingsquare1}
\left\langle e_j,\Phi(s)e_j\right\rangle_2 = \lambda L\text{!}\sum_{\ell\in\Sigma}\left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\prod_{i=0}^\infty
\frac{\pi_i(t)^{\rho_i(\ell)}}{\rho_i(\ell)\text{!}}
+k\sum_{i=1}^\infty
\left\langle e_j,(e_{i-1}-e_i)(e_{i-1}-e_i)^Te_j\right\rangle_2\pi_i(t).
\end{align}
Recalling from (28) that, for ℓ∈ Σj(i 1, i 2, i 3), (29) holds, we have from the decomposition,
\begin{align*}
\Sigma = \bigcup\limits_{{i_1} = 0}^{k - 2} \bigcup\limits_{{i_2} = 0}^{L - {i_1}} \bigcup\limits_{{i_3} = 0}^{L - {i_1} - {i_2}} {\Sigma _j}({i_1},{i_2},{i_3}),
\end{align*}
and
\begin{align*}
\lambda L\text{!}\,\sum_{\ell\in\Sigma}\left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\prod_{i=0}^\infty\frac{\pi_i(t)^{\rho_i
(\ell)}}{\rho_i(\ell)\text{!}}
=\lambda\sum_{i_1=0}^{k-2}\sum_{i_2=0}^{L-i_1}\sum_{i_3=0}^{L-i_1
-i_2}\sum_{\ell\in\Sigma_j(i_1,i_2,i_3)}\left\langle
e_j,\Delta_\ell\Delta_\ell^Te_j\right\rangle_2\binom{L}{i_1,i_2,i_3,L-i_1-i_2-i_3}
\times\binom{i_1}{\rho_1(\ell),\ldots,\rho_{j-2}(\ell)}
\binom{L-i_1-i_2-i_3}{\rho_{j+1}(\ell),\rho_{j+2}(\ell),\ldots}
\prod_{i=0}^\infty\pi_i(t)^{\rho_i(\ell)}.
\end{align*}
\begin{align*}
=\lambda\sum_{i_1=0}^{k-2}\sum_{i_2=0}^{L-i_1}\sum_{i_3=0}^{L-i_1
-i_2}\binom{L}{i_1,i_2,i_3,L-i_1-i_2-i_3}
\times\bigg(\sum_{m=0}^{j-2}\pi_m(s)\bigg)^{i_1}\pi_{j-1}
(s)^{i_2}\pi_j(s)^{i_3}\bigg(\sum_{m=j+1}^{\infty}\pi_{m}(s)\bigg)^{L-i_1-
i_2-i_3}
\times[i_2\wedge (k-i_1)_+-i_3\wedge (k-i_1-i_2)_+]^2
=\bar{Z}(j,\pi(s)),
\end{align*
where the second equality follows from the multinomial theorem. Futhermore, from (25), the second sum in (A.1) is
\begin{align*}
k\sum_{i=1}^\infty \left\langle e_j,(e_{i-1}-e_i)(e_{i-1}-e_i)^Te_j\right\rangle_2\pi_i(t)
= k(\pi_j(t)+\pi_{j+1}(t)).
\end{align*}
Combining the last two equations with (A.1) gives (97).
For (99), note that
\begin{align*}\label{eqn:Zbarineq}
\bar{Z}(j,\pi(s)) \leq
\sum_{i_1=0}^{k-1}\frac{(\sum_{m=0}^{j-2}\pi_m(s))^{i_1}}{i_1\text{!}}
\sum_{i_2=0}^{L-i_1}\frac{(\pi_{j-1}(s))^{i_2}}{i_2\text{!}}
\\
\times\sum_{i_3=0}^{L-i_1-i_2}k^2{\bf 1}_{\{i_2\vee i_3 >
0\}}\frac{(\pi_j(s))^{i_3}}{i_3\text{!}}\frac{(\sum_{m=j+1}^{\infty}
\pi_{m}(s))^{L-i_1-i_2-i_3}}{(L-i_1-i_2-i_3)\text{!}}
\\
\leq
\sum_{i_1=0}^{k-1}\sum_{i_2=0}^{L-i_1}\sum_{i_3=0}^{L-i_1-i_2}k^2{\bf
1}_{\{i_2\vee i_3 > 0\}}(\pi_{j-1}(s))^{i_2}(\pi_j(s))^{i_3}
\\
\leq
c_{\bar{Z}}(\pi_{j-1}(s)+\pi_j(s))
\end{align*}
for some c Z̄ ∈ (0, ∞), where the second inequality follows from the fact that all factorials are greater than 1 and π is a probability measure. This proves (99).
Acknowledgements
The research was supported in part by the National Science Foundation (DMS-1016441, DMS-1305120, DMS-1814894), the Army Research Office (W911NF-14-1-0331), and DARPA (W911NF-15-2-0122).
