#  Invariance principle and McKean-Vlasov limit for randomized load balancing
in heavy traffic

Rami Atar  Viterbi Faculty of Electrical and Computer Engineering
Technion – Israel Institute of Technology<EMAIL_ADDRESS>and  Gershon
Wolansky  Department of Mathematics
Technion – Israel Institute of Technology<EMAIL_ADDRESS>
(Date: February 28, 2024)

###### Abstract.

We consider a load balancing model where a Poisson stream of jobs arrive at a
system of many servers whose service time distribution possesses a finite
second moment. A small fraction of arrivals pass through the so called power-
of-choice algorithm, which assigns a job to the shortest among $\ell$,
$\ell\geq 2$, randomly chosen queues, and the remaining jobs are assigned to
queues chosen uniformly at random. The system is analyzed at critical load in
an asymptotic regime where both the number of servers and the usual heavy
traffic parameter associated with individual queue lengths grow to infinity.
The first main result is a hydrodynamic limit, where the empirical measure of
the diffusively normalized queue lengths is shown to converge to a path in
measure space whose density is given by the unique solution of a parabolic PDE
with nonlocal coefficients.

Further, two forms of an invariance principle are proved, corresponding to two
different assumptions on the initial distribution, where individual normalized
queue lengths converge weakly to solutions of SDE. In one of these results,
the limit is given by a McKean-Vlasov SDE, and propagation of chaos holds. The
McKean-Vlasov limit is closely related to limit results for Brownian particles
on ${\mathbb{R}}_{+}$ interacting through their rank (with a specific
interaction). However, an entirely different set of tools is required, as the
collection of $n$ prelimit particles does not obey a Markovian evolution on
${\mathbb{R}}_{+}^{n}$.

The PDE has a stationary solution expressed explicitly, which in particular
provides a quantification of the balance achieved by the algorithm. For fixed
$\ell$, as the intensity of the load balanced stream varies between its
limits, this solution varies from exponential distribution to a Dirac
distribution, demonstrating that the results cover the whole range from
independence to full coordination.

###### Key words and phrases:

Randomized load balancing, power of choice, low sampling rate, diffusion
limit, hydrodynamic limit, McKean-Vlasov limit, propagation of chaos,
parabolic initial boundary value problem, viscous scalar conservation law,
Daley-Miyazawa semimartingale representation.

###### 2010 Mathematics Subject Classification:

68M20, 90B22, 60F17, 35K20

##  1\. Introduction

This paper is concerned with a load balancing model where a Poisson stream of
arrivals faces a system of $n$ servers working in parallel, with a general
common service time distribution having a finite second moment. Routing is
based on the join the shortest of $\ell$ queues (abbreviated JSQ($\ell$))
algorithm, which assigns a job to the shortest among $\ell$, $\ell\geq 2$,
queues chosen uniformly at random. Following a setting proposed in [6],
motivated by maintaining low communication overhead, the stream undergoes
thinning, where a small fraction of the arrivals are routed via JSQ($\ell$)
while the other arrivals are assigned to a queue chosen uniformly at random
among all queues.

The focus is on the asymptotics where the number of servers, $n$, becomes
large, while individual queues are kept at their diffusive regime. This is
achieved by taking $n$ to also serve as the heavy traffic parameter of
individual queue lengths, letting the per-server processing rate and arrival
intensity scale like $n$, and critically loading the servers (to first order).
In order for the load balanced stream to have a non-negligible effect, its
average intensity per queue should scale (at least) like $n^{1/2}$. Hence it
is assumed that the intensity of the overall arrival stream, the overall load
balanced stream, and the per-server arrival stream scale like $n^{2}$,
$n^{3/2}$ and $n$, respectively (see Fig. 1). Specifically, we take the
overall load balanced intensity asymptotic to $bn^{3/2}$, and the average per-
server arrival asymptotic to $\lambda n$, where $b>0$ and $\lambda>0$ are
constants. Moreover, the average load on each queue, defined as the rate of
arrival of work per server divided by a server’s processing rate, is
asymptotic to $1+\rho n^{-1/2}$, where $\rho\in{\mathbb{R}}$, and stability is
expresses by the condition $\rho<0$. The $i$-th queue length process at the
$n$-th system is denoted by $X^{n}_{i}(t)$, and its normalized version by
$\hat{X}^{n}_{i}(t)=n^{-1/2}X^{n}_{i}(t)$, $i\in[n]:=\\{1,\ldots,n\\}$. The
empirical distribution at time $t$ is denoted by
$\bar{\xi}^{n}_{t}=n^{-1}\sum_{i\in[n]}\delta_{\hat{X}^{n}_{i}(t)}$.

The parameters $b$ and $\ell$ determine the volume of communication between
the servers and the dispatcher. They also have dramatic influence on the
degree to which the system is load balanced. One of the goals of this work is
to quantify the level of achieved load balancing, which could be defined e.g.
by the empirical standard deviation, or other functionals of
$\bar{\xi}^{n}_{t}$, as a function of these parameters. Assuming
$\bar{\xi}^{n}_{0}\to\xi_{0}$ in probability, $\xi_{0}$ a Borel probability
measure on ${\mathbb{R}}_{+}$, our first main result addresses this by
providing a hydrodynamic limit, which shows that $\bar{\xi}^{n}_{t}$ converges
in probability, as a measure-valued process, to a deterministic measure-valued
path given by $\xi_{t}(dx)=u(x,t)dx$, $t>0$. Here, the density $u$ is the
unique classical solution of the initial-boundary value problem

(1.1) $\begin{cases}u_{t}=-[(b\ell
v^{\ell-1}-c_{1})u]_{x}+au_{xx}&(x,t)\in(0,\infty)^{2}\\\ \displaystyle
v(x,t)=\int_{x}^{\infty}u(y,t)dy&(x,t)\in(0,\infty)^{2}\\\
(c_{1}-b\ell)u(0,t)+au_{x}(0,t)=0&t\in(0,\infty)\\\
u(\cdot,t)dx\to\xi_{0}(dx)&\text{weakly as $t\downarrow 0$.}\end{cases}$

Above, $c_{1}=-\lambda\rho+b$, and $a>0$ is a variance parameter. The equation
has a non-local drift coefficient that captures the algorithm’s instantaneous
effect on the macroscopic distribution, and a Robin boundary condition
expressing zero flux through the origin.

When $\rho<0$, (1.1) has a unique stationary solution given by

(1.2) $\begin{split}v_{\rm stat}(x)&=w(x)^{-1/(\ell-1)},\qquad u_{\rm
stat}(x)=\frac{c_{1}}{a}(1-\alpha)e^{\frac{c_{1}}{a}(\ell-1)x}w(x)^{-\ell/(\ell-1)},\\\
w(x)&=(1-\alpha)e^{\frac{c_{1}}{a}(\ell-1)x}+\alpha,\qquad\alpha=\frac{b}{c_{1}}=\Big{(}1-\lambda\frac{\rho}{b}\Big{)}^{-1}\in(0,1),\end{split}$

for $x\in{\mathbb{R}}_{+}$. It can be seen that for fixed
$(\ell,\rho,\lambda,a)$, as $b\to 0$ and, respectively, $b\to\infty$, the
measure $u_{\rm stat}(x)dx$ converges to the exponential distribution and to
$\delta_{0}(dx)$, showing that the setting covers the entire range from queues
operating independently to a state space collapse where the load balancing
algorithm achieves complete coordination.

Our second goal is to provide an invariance principle. Under the additional
assumption that, for a fixed $k$,
$(\hat{X}^{n}_{i}(0))_{i\in[k]}\Rightarrow(X_{i}(0))_{i\in[k]}$, it is shown
that $(\hat{X}^{n}_{i})_{i\in[k]}\Rightarrow(X_{i})_{i\in[k]}$ where the
latter is the solution to the SDE in ${\mathbb{R}}_{+}^{k}$

(1.3)
$X_{i}(t)=X_{i}(0)+b_{1}t+b_{0}\int_{0}^{t}v(X_{i}(s),s)^{\ell-1}ds+\sigma
W_{i}(t)+L_{i}(t),\qquad i\in[k],$

in which $L_{i}$ are boundary terms having continuous nondecreasing sample
paths, satisfying $\int_{[0,\infty)}X_{i}(t)dL_{i}(t)=0$. Alternatively, if
the initial queue lengths are assumed exchangeable, then the limit in
distribution of the $k$-tuple is given by $k$ independent copies of

(1.4) $X(t)=X(0)+b_{1}t+b_{0}\int_{0}^{t}v(X(s),s)^{\ell-1}ds+\sigma
W(t)+L(t),$

with $\int_{[0,\infty)}X(t)dL(t)=0$. In both cases, the diffusion limits
depend on the service time distribution only through its first two moments.

The last result, combined with the characterization of the limiting empirical
distribution in terms of (1.1), shows that the pair $(X,\bar{F})$, describing
the limiting stochastic dynamics of a rescaled queue length and its law,
satisfies

(1.5) $X(t)=X(0)+\int_{0}^{t}\mathbf{b}(X(s),\bar{F}(\cdot,s))ds+\sigma
W(t)+L(t),\qquad\bar{F}(x,t)={\mathbb{P}}(X(t)>x),$

where $\mathbf{b}(x,\bar{F})=b_{1}+b_{0}\bar{F}(x)^{\ell-1}$. Viewed this way,
the result is closely related to the literature on interacting diffusion
models, specifically diffusions interacting through their ranks, and is an
instance of a McKean-Vlasov limit. For background on McKean-Vlasov limits and
propagation of chaos we refer to [41, 21]. To explain the relation, consider a
parametric regime in which the heavy traffic parameter is taken to its limit
first. In this case, each normalized queue length behaves as a Brownian
particle on ${\mathbb{R}}_{+}$, with an interaction among the particles caused
by the load balancing algorithm. The algorithm effects the dynamics by
selecting a particle with probability depending on its rank. In the original
model, the selected queue increases by one job. In the particle system this
can be imitated by imposing a positive rank-dependent drift. It therefore
comes as no surprise that the McKean-Vlasov SDE (1.5) is a special case of the
one that arises in the study of rank-dependent diffusions. Indeed, existence
and uniqueness of solutions to (1.5), as well as convergence results, follow
from work on interacting diffusions, specifically [29, Sec. 2.4]. Clearly, the
above explanation is only heuristic. In the regime under consideration, the
prelimit objects are queue lengths, which do not follow a Markovian evolution
when considered as an $n$-tuple with state space ${\mathbb{R}}_{+}^{n}$. As
the basic assumptions of the interacting diffusion setting do not hold, the
treatment requires a completely different set of tools.

###  1.1. Motivation and related work

A large body of work has been devoted to subcritically loaded systems, where
on average queues are short. Other works analyzed critically loaded systems
for which either state space collapse or near optimal performance were
achieved; see references below. However, as far as many-server limits are
concerned, non perfect balancing has not been quantified before and very few
papers have addressed asymptotic behaviors in which individual queue lengths
exhibit diffusive fluctuations (in a setting with a fixed number of servers,
[6] considered both aspects). These two aspects are of great practical
relevance.

Invariance principles are desirable as they express robustness of performance
to underlying distributions. In the current context, an additional important
motivation for an invariance principle stems from the fact that, to the best
of our knowledge, all earlier diffusion limit results on load balancing in
many-server settings cover only the case of exponential service times, a
rather restrictive assumption from a modeling viewpoint.

This paper is not the first to study a regime in which the number of servers
scales like $n$ and individual queues are critically loaded, with each queue
length scaling like $n^{1/2}$. This precise parameterization was considered
before in [9], which studied a multi-agent game rather than a load balancing
model.

The first paper to study an asymptotic regime that combines heavy traffic
scaling with many-server scaling was [26], where a single-queue pooled
$n$-server system was considered. The regime introduced there, known as the
Halfin-Whitt regime, is one where the probability of a job to meet an
available server upon arrival is bounded away from $0$ and $1$, and is, in
particular, distinct from the regime studied here.

The literature on load balancing in asymptotic regimes is vast; see [18] for a
recent survey. The asymptotic regime where individual queues undergo heavy
traffic scaling and JSQ($\ell$) is applied on a small fraction of the arrival
stream was introduced in the aforementioned [6], in a setting where the number
of servers is fixed. This could model two different scenarios, that are also
relevant for our model. One is where a single stream is split by the
dispatcher. Another is where each server has a dedicated stream due, for
example, to geographical or compatibility constraints, and an additional
stream is shared by all servers. Under an exponential service assumption, it
was proved in [6] that the collection of rescaled queue lengths converges to a
rank-based diffusion.

Perhaps the most well known load balancing algorithm is the join-the-shortest-
queue (JSQ), which routes jobs the shortest among all queues, known to be
delay optimal under exponential service and asymptotically delay optimal in
various limiting regimes; see, for example, [15]. Its variation JSQ($\ell$),
that is typically applied with $\ell$ much smaller that $n$, significantly
reduces the communication overhead while maintaining good performance. It was
introduced in [43], where under the exponential service time assumption and
subcritical load, the $n\to\infty$ limit of the length of a typical queue in
stationarity was shown to have a doubly exponential tail decay, a dramatic
improvement over routing uniformly at random in which case decay is
exponential. Related results appeared in [34]. The result was extended in [23]
where empirical measures were considered and propagation of chaos was
obtained. Functional central limit theorems (CLT) were obtained in [24] and
strong approximation results, including law of large numbers (LLN) and CLT, in
[32]. (These CLT results were not invariance principles and were not concerned
with the asymptotic regime studied in this paper). Results on the mixing rate
and the size of the maximal queue length were obtained in [31], and aspects of
stability and performance under server heterogeneity were investigated in
[36]. The paper [35] studied how $\ell(n)$ should grow so that JSQ$(\ell(n))$
with exponential servers would perform like JSQ, and thus achieve asymptotic
delay optimality.

A series of papers [11, 12, 13] extended [43] to several families of general
service time distributions, based on an approach that first establishes
propagation of chaos and then uses it to compute queue length distribution in
equilibrium. Another line of research treating JSQ($\ell$) with general
service times is [2, 3] and [1]. In [2, 3], the dynamical behavior was studied
under general load and the hydrodynamic limit was shown to be given by an
infinite system of coupled PDE. A numerical method was developed to solve
these PDE. In [1], an infinite system of PDE was constructed and shown to
constitute the invariant state of the aforementioned system of PDE under a
subcriticality condition. Our setting, where individual queues undergo heavy
traffic scaling, as well as our results, are quite different from both these
lines of work.

The diffusion limit in heavy traffic of JSQ in a many server setting was
established in [19], where, assuming exponential servers, a rescaled empirical
measure of queue lengths was shown to converge to a process expressed in terms
of a 2d diffusion. Convergence of the steady state at the same scale was
proved in [14], and properties of the diffusion process were investigated in
[7, 8]. Under the regime considered in this line of work, the number of queues
that are of length $0$, $1$ and $2$ is of order $n^{1/2}$, $n$ and,
respectively, $n^{1/2}$, and only a negligible fraction exceeds length $2$.
Thus this regime captures quite a different behavior from what is described in
this paper. The papers [25] and [46] considered JSQ in other diffusive
regimes, namely the nondegenerate slowdown regime and, respectively, the
super-Halfin-Whitt regime.

Numerous load balancing algorithms besides JSQ and JSQ($\ell$) have been
proposed. Ones that emphasize sparse communication, where the messaging rate
often goes far below that of JSQ($\ell$), include [42, 45, 20, 4]. The
practical importance of sparse messaging is widely acknowledged; see e.g.
[33]. Further asymptotic results on load balancing with nonexponential service
times include join-the-idle-queue [20], pull-based load distribution [39, 40],
zero waiting algorithms [30], and join-the-shortest-estimated-queue [5].

McKean-Vlasov limits of diffusions interacting through their ranks were
addressed in [38]. Observing that in these models the dependence of each
particle’s coefficients on the empirical law is discontinuous, a situation not
covered by the classical treatment such as [21], this paper gave convergence
results assuming that the coefficients are merely measurable. The paper [29]
mentioned above addresses a setting of which rank-based interaction is a
special case, and attains convergence results in a stronger topology than weak
convergence.

###  1.2. Notation

Denote ${\mathbb{R}}_{+}=[0,\infty)$. Let
$\iota:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ denote the identity map. In
${\mathbb{R}}^{N}$, denote the Euclidean norm by $\|\cdot\|$. For
$({\mathbf{X}},d_{\mathbf{X}})$ a Polish space, let
$C({\mathbb{R}}_{+},{\mathbf{X}})$ and $D({\mathbb{R}}_{+},{\mathbf{X}})$
denote the space of continuous and, respectively, càdlàg paths, endowed with
the topology of uniform convergence on compacts and, respectively, the $J_{1}$
topology. Denote by ${\mathcal{M}}_{1}$ the space of probability measures on
${\mathbb{R}}_{+}$ equipped with the topology of weak convergence. Denote by
$C^{\uparrow}$ the set of members of $C({\mathbb{R}}_{+},{\mathbb{R}}_{+})$
that are nondecreasing and start at $0$. For $\xi\in
D({\mathbb{R}}_{+},{\mathbb{R}}^{N})$, an interval $I\subset{\mathbb{R}}_{+}$,
and $0\leq\delta\leq T$, denote

$\displaystyle\mathnormal{\Delta}\xi(t)$ $\displaystyle=\xi(t)-\xi(t-)\qquad
t>0,\qquad\mathnormal{\Delta}\xi(0)=\xi(0),$ $\displaystyle{\rm osc}(\xi,I)$
$\displaystyle=\sup\\{\|\xi(s)-\xi(t)\|:s,t\in I\\},$ $\displaystyle
w_{T}(\xi,\delta)$
$\displaystyle=\sup\\{\|\xi(t)-\xi(s)\|:s,t\in[0,T],|s-t|\leq\delta\\},$
$\displaystyle\|\xi\|^{*}_{T}$ $\displaystyle=\sup\\{\|\xi(t)\|:t\in[0,T]\\},$

and by $|\xi|(t)$ the total variation of $\xi$ in $[0,t]$.

For $\mathfrak{D}\subset{\mathbb{R}}^{N}$, denote by $\mathfrak{D}^{o}$ and
$\bar{\mathfrak{D}}$ its interior and, respectively, closure. Denote by
$C(\mathfrak{D})$ the set of continuous functions
$f:\mathfrak{D}\to{\mathbb{R}}$ and by $C_{0}(\mathfrak{D})$ the set of $f\in
C(\mathfrak{D})$ whose support is a compact subset of ${\mathbb{R}}^{N}$. For
$k,l\in{\mathbb{N}}$ and an open set $\mathfrak{D}\subset{\mathbb{R}}$ (resp.,
$\mathfrak{D}\subset{\mathbb{R}}\times{\mathbb{R}}_{+}$), denote by
$C^{k}(\mathfrak{D})$ (resp., $C^{k,l}(\mathfrak{D})$) the set of functions
$f:\mathfrak{D}\to{\mathbb{R}}$ possessing continuous derivatives up to and
including $k$ (resp., $(k,l)$). For $\mathfrak{D}\subset{\mathbb{R}}$,
$C^{k}(\mathfrak{D})$ denotes the set of functions
$f:\mathfrak{D}\to{\mathbb{R}}$ whose restriction to $\mathfrak{D}^{o}$ lies
in $C^{k}(\mathfrak{D}^{o})$, and whose derivatives up to and including $k$
have continuous extensions to $\mathfrak{D}$. Define $C^{k,l}(\mathfrak{D})$
analogously. For $\mathfrak{D}\in{\mathbb{R}}$ and
$\mathfrak{D}\in{\mathbb{R}}\times{\mathbb{R}}_{+}$ denote
$C^{\infty}(\mathfrak{D})=\bigcap_{k\in{\mathbb{N}}}C^{k}(\mathfrak{D})$ and
$C^{\infty}(\mathfrak{D})=\bigcap_{k,l\in{\mathbb{N}}}C^{k,l}(\mathfrak{D})$,
resp. A subscript $0$ denotes compact support, e.g.
$C^{k}_{0}(\mathfrak{D})=C^{k}(\mathfrak{D})\cap C_{0}(\mathfrak{D})$. A
subscript $b$ denotes boundedness, e.g. $C^{k}_{b}(\mathfrak{D})$ are
functions in $C^{k}(\mathfrak{D})$ whose derivatives of order $0\leq l\leq k$
are bounded. For $f$ defined on (a subset of) ${\mathbb{R}}$, $f^{\prime}$
denotes derivative, and for $f$ defined on (a subset of)
${\mathbb{R}}\times{\mathbb{R}}_{+}$, $f_{x}$ and $f_{t}$ denote spatial and
temporal derivatives, resp.

For an interval $I\subset{\mathbb{R}}_{+}$ and a normed space $S$, denote by
$\mathbb{L}^{p}(I;S)$, $1\leq p\leq\infty$, the usual $\mathbb{L}^{p}$ space
defined in terms of the Lebesgue measure on $I$. Denote by
$\mathbb{L}^{p}_{\rm loc}({\mathbb{R}}_{+};S)$ the set of
$f:{\mathbb{R}}_{+}\to S$ such that $f\in\mathbb{L}^{p}([0,T];S)$ for all
finite $T$. $\mathbb{L}^{p}({\mathbb{R}}_{+};{\mathbb{R}})$ is abbreviated to
$\mathbb{L}^{p}({\mathbb{R}}_{+})$.

For $f,g:{\mathbb{R}}_{+}\to{\mathbb{R}}$ and a measure $m$ on
${\mathbb{R}}_{+}$, denote $\langle
f,m\rangle=\int_{{\mathbb{R}}_{+}}f(x)m(dx)$ and $\langle
f,g\rangle=\int_{{\mathbb{R}}_{+}}f(x)g(x)dx$.

If $V\in{\mathbb{R}}^{n}$ then $V_{i}$, $i\in[n]$ denote its components in the
standard basis, and vice versa: Given $V_{i}$, $i\in[n]$, $V$ denotes the
vector $(V_{1},V_{2},\ldots,V_{n})$. Both these conventions hold also for
random variables $V_{i}$ and processes $V_{i}(\cdot)$. Occasionally, with a
slight abuse of standard terminology, a sequence of random elements (random
variables or processes) is referred to as tight when their laws form a tight
sequence of probability measures. The term with high probability (w.h.p.)
means “holds, for each $n$, on $\mathnormal{\Omega}_{n}$, where
$\lim_{n}{\mathbb{P}}(\mathnormal{\Omega}_{n})=1$”. Convergence in
distribution is denoted by $\Rightarrow$. $c$ denotes a positive constant
whose value may change from one expression to another.

Throughout the paper, superscript $n$ attached to scalars, random variables or
processes, denotes dependence on the index $n$ rather than power.

###  1.3. Paper organization

The load balancing model and the main results are presented in §2.1 and §2.2,
respectively. A discussion of the results appears in §2.3. An outline of the
proof is given in §2.4. Several tools required for the proof are developed in
§3. The proof of the main results appears in §4. For a detailed description of
the content of §§3–4 see §2.4.

Figure 1. Orders of magnitude of the various streams.

##  2\. Load balancing in heavy traffic

###  2.1. The load balancing model

####  2.1.1. Arrivals, queue lengths and basic relations

In the model there are $n$ servers and a queue in front of each. There is a
dedicated stream of arrivals into each queue and an additional stream of
arrivals, called the load balancing stream (LBS), that go through the
JSQ$(\ell,n)$ algorithm. These $n+1$ arrival streams are modeled as mutually
independent Poisson processes. Clearly, this could be recast as a single
Poisson stream, out of which $n+1$ thinned streams are created by means of
random selection.

In what follows, the queueing systems will be indexed by $n\in{\mathbb{N}}$,
the number of servers. The processes $X^{n}_{i}$, $E^{n}_{i}$, $D^{n}_{i}$ and
$T^{n}_{i}$ represent the $i$-th queue length process, dedicated arrival
process, departure process and busyness process, respectively. Denote by
$A^{n}_{0}$ the LBS arrival process, and by $A^{n}_{i}$ the process counting
LBS arrivals routed to server $i$. For each $i$, $E^{n}_{i}$ is a Poisson
process of parameter $\lambda^{n}$, and $A^{n}_{0}$ is Poisson of parameter
$\lambda^{n}_{0}$, all having right-continuous sample paths. We have

(2.1)
$X^{n}_{i}(t)=X^{n}_{i}(0-)+E^{n}_{i}(t)+A^{n}_{i}(t)-D^{n}_{i}(t),\qquad
i\in[n],\ t\in{\mathbb{R}}_{+}.$

Work conservation is assumed, hence

(2.2) $T^{n}_{i}(t)=\int_{0}^{t}1_{\\{X^{n}_{i}(s)>0\\}}ds.$

####  2.1.2. The load balancing algorithm

Upon each LBS arrival, $\ell$ out of the $n$ queues, chosen uniformly at
random, are sampled. The arrival is routed to the queue that is shortest among
the $\ell$; if there are ties, the queue with the smaller index is preferred.
Given $x\in{\mathbb{R}}^{n}$, let

(2.3) ${\rm rank}(i;x)=\\#\\{j:x_{j}<x_{i}\\}+\\#\\{j\leq i:x_{j}=x_{i}\\}.$

Then the probability that a LBS arrival is routed to the queue whose rank,
defined by (2.3), is $r$ is given by

(2.4) $p_{n,r}=\frac{{n-r\choose\ell-1}}{{n\choose\ell}},\qquad r\in[n],$

with ${k\choose j}=0$ when $j>k$. Note that

(2.5) $p_{n,r}=0\text{ for }r\geq
n-\ell+2,\qquad\max_{r\in[n]}p_{n,r}=p_{n,1}=\frac{\ell}{n}.$

Thus the randomization mechanism can equivalently be achieved by letting, for
each $n$, $\theta^{n}_{k}$, $k\in{\mathbb{N}}$ be IID random variables with
${\mathbb{P}}(\theta^{n}_{1}=r)=p_{n,r}$, $r\in[n]$, and routing the $k$-th
LBS arrival to the queue whose rank is $\theta^{n}_{k}$. That is,

(2.6)
$A^{n}_{i}(t)=\int_{[0,t]}1_{\\{{\mathcal{R}}^{n}_{i}(s-)=\theta^{n}_{A^{n}_{0}(s)}\\}}dA^{n}_{0}(s),\qquad{\mathcal{R}}^{n}_{i}(t)={\rm
rank}(i;X^{n}(t)),\qquad i\in[n],\ t\in{\mathbb{R}}_{+}.$

Of course, this randomization scheme is not used in practice as it requires
keeping track of the queue lengths of the entire system. However, it is
convenient for carrying out the analysis.

####  2.1.3. The initial condition.

We allow quite a general initial condition, where residual times of jobs
already being processed at time $0$ may be dependent and have unspecified
distributions. Denote the (random) set of queues that at time $0$ contain no
jobs and, respectively, at least one job, by
${\mathcal{N}}^{n}=\\{i\in[n]:X^{n}_{i}(0)=0\\}$ and
${\mathcal{P}}^{n}=\\{i\in[n]:X^{n}_{i}(0)>0\\}$. For $i\in{\mathcal{P}}^{n}$,
let $Z^{n}_{i}(0)$ denote the initial residual time of the head-of-the-line
job in queue $i$. For $i\in{\mathcal{N}}^{n}$ we add fictitious jobs having
zero processing time. To this end, rather than specifying $X^{n}(0)$ as the
initial queue length, $X^{n}(0-)$ is specified; and for each
$i\in{\mathcal{N}}^{n}$, the queue length is set to $X^{n}_{i}(0-)=1$ and the
residual processing time to $Z^{n}_{i}(0)=0$. Obviously, this results in
$X^{n}_{i}(0)=0$. This convention allows us to greatly simplify notation when
we later construct counting processes for service and departure. The initial
condition is thus a tuple

${\mathcal{I}}^{n}=(\\{X^{n}_{i}(0-),Z^{n}_{i}(0),i\in[n]\\},{\mathcal{N}}^{n},{\mathcal{P}}^{n}),$

where $({\mathcal{N}}^{n},{\mathcal{P}}^{n})$ partitions $[n]$, and

$\begin{split}X^{n}_{i}(0-)&=1,\ Z^{n}_{i}(0)=0,\ i\in{\mathcal{N}}^{n},\\\
X^{n}_{i}(0-)&\geq 1,\ Z^{n}_{i}(0)>0,\ i\in{\mathcal{P}}^{n}.\end{split}$

####  2.1.4. Service times

Let $\Phi_{\rm ser}$ be a Borel probability measure on $[0,\infty)$ with mean
1 and standard deviation $\sigma_{\rm ser}\in(0,\infty)$, such that $\Phi_{\rm
ser}(\\{0\\})=0$. Let $\Phi^{n}_{\rm ser}$ be defined as a scaled version of
$\Phi_{\rm ser}$ uniquely specified via $\Phi^{n}_{\rm ser}[0,x]=\Phi_{\rm
ser}[0,\mu^{n}x]$, $x\in{\mathbb{R}}_{+}$. Here, $\mu^{n}>0$ is the service
rate in the $n$-th system. For $k\geq 1$, let $Z^{n}_{i}(k)$ denote the
service time of the $k$-th job to be served by server $i$ after the head-of-
the-line job at time $0-$ there (for $i\in{\mathcal{N}}^{n}$ this means the
$k$-th job after the fictitious one). For every $i$,
${\mathcal{Z}}^{n}_{i}=(Z^{n}_{i}(k),k\geq 1)$ is an IID sequence with common
distribution $\Phi^{n}_{\rm ser}$.

Next, the potential service process $S^{n}_{i}$, evaluated at $t$, gives the
number of jobs completed by server $i$ by the time it has worked $t$ units of
time. It is given, with $\sum_{0}^{-1}=0$, by

$S^{n}_{i}(t)=\max\Big{\\{}k\in{\mathbb{Z}}_{+}:\sum_{j=0}^{k-1}Z^{n}_{i}(j)\leq
t\Big{\\}},\qquad t\geq 0.$

The departure processes are given by $D^{n}_{i}(t)=S^{n}_{i}(T^{n}_{i}(t))$.
This is the number of jobs completed by time $t$ by server $i$. Note that the
first departure counted by $D^{n}_{i}$ is the one initially processed if
$i\in{\mathcal{P}}^{n}$ and the fictitious job if $i\in{\mathcal{N}}^{n}$.

####  2.1.5. Dependence structure.

For each $n$, the $2n+3$ stochastic elements

(2.7) $E^{n}_{i}$, $i\in[n]$, ${\mathcal{Z}}^{n}_{i}$, $i\in[n]$,
${\mathcal{I}}^{n}$, $A^{n}_{0}$, and $\\{\theta^{n}_{k}\\}$ are mutually
independent.

####  2.1.6. Scaling and critical load condition.

The arrival and service rates are assumed to satisfy the following. There are
constants $\lambda>0$ and $\hat{\lambda}\in{\mathbb{R}}$ such that

(2.8) $\hat{\lambda}^{n}:=n^{-1/2}(\lambda^{n}-n\lambda)\to\hat{\lambda}\text{
as }n\to\infty,$

a constant $b>0$ such that

(2.9) $\hat{\lambda}^{n}_{0}:=n^{-3/2}\lambda^{n}_{0}\to b\text{ as
}n\to\infty,$

and constants $\mu>0$ and $\hat{\mu}\in{\mathbb{R}}$ such that

(2.10) $\hat{\mu}^{n}:=n^{-1/2}(\mu^{n}-n\mu)\to\hat{\mu}\text{ as
}n\to\infty.$

The critical load condition is assumed, namely

(2.11) $\lambda=\mu.$

Some further notation used throughout is

$\hat{b}^{n}_{1}=\hat{\lambda}^{n}-\hat{\mu}^{n},\qquad
b_{1}=\hat{\lambda}-\hat{\mu},\qquad c_{1}=-b_{1},\qquad
b_{0}=b\ell,\qquad\sigma^{2}=\lambda(1+\sigma_{\rm ser}^{2}),\qquad
a=\frac{\sigma^{2}}{2}.$

Note that the average load on a queue, defined as the rate of arrival of work
per server, $\lambda^{n}+n^{-1}\lambda^{n}_{0}$, divided by a server’s
processing rate, $\mu^{n}$, that is,

$\frac{n\lambda+n^{1/2}\hat{\lambda}+n^{1/2}b+o(n^{1/2})}{n\mu+n^{1/2}\hat{\mu}+o(n^{1/2})},$

is asymptotic to

$1+n^{-1/2}\rho,\qquad\rho:=\frac{\hat{\lambda}+b-\hat{\mu}}{\lambda}.$

Hence the load parameter $\rho$ is given by $\rho=(b_{1}+b)/\lambda$.

Let rescaled versions of the queue length processes and cumulative idle time
processes be defined by

(2.12)
$\hat{X}^{n}_{i}(t)=n^{-1/2}X^{n}_{i}(t),\qquad\hat{L}^{n}_{i}(t)=n^{-1/2}\mu^{n}(t-T^{n}_{i}(t)),$

and denote
$\bar{\xi}^{n}_{t}=n^{-1}\sum_{i\in[n]}\delta_{\hat{X}^{n}_{i}(t)}$. Further,
assume that $\bar{\xi}^{n}_{0-}\to\xi_{0}$ (equivalently
$\bar{\xi}^{n}_{0}\to\xi_{0}$) in probability, as $n\to\infty$, where
$\xi_{0}$ is a deterministic Borel probability measure on ${\mathbb{R}}_{+}$.
For the rescaled residual times $Z^{n}_{i}(0)$, assume that for every
$\varepsilon>0$,

(2.13)
$\lim_{n\to\infty}\max_{i\in[n]}{\mathbb{P}}(\tilde{Z}^{n}_{i}(0)>\varepsilon)=0\quad\text{
where}\quad\ \tilde{Z}^{n}_{i}(0):=\tilde{\mu}^{n}Z^{n}_{i}(0),\
\tilde{\mu}^{n}:=n^{-1/2}\mu^{n},$

and

(2.14) $\sup_{n}\max_{i\in[n]}{\mathbb{E}}[\tilde{Z}^{n}_{i}(0)^{2}]<\infty.$

Because $\tilde{\mu}^{n}$ scales like $n^{1/2}$, (2.13) and (2.14) impose the
condition that $Z^{n}_{i}(0)$ scale like $o(n^{-1/2})$. This is mild compared
to the condition on $Z^{n}_{i}(k)$, $k\geq 1$, which have been assumed to
scale like $n^{-1}$. Finally, assume that

(2.15) $\sup_{n}\max_{i\in[n]}{\mathbb{E}}[\hat{X}^{n}_{i}(0-)^{2}]<\infty.$

Note that, by Fatou’s lemma, this imposes a condition on $\xi_{0}$, namely
$\xi_{0}$ necessarily satisfies $\int x^{2}\xi_{0}(dx)<\infty$.

All assumptions made thus far are in force throughout the paper.

###  2.2. Main results

First we address well-posedness of the relevant PDE, starting with classical
solutions to (1.1).

######  Theorem 2.1.

Within the class of functions $u\in
C^{2,1}((0,\infty)\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$ satisfying, for
every $T<\infty$, $\sup_{t\in(0,T]}\int_{{\mathbb{R}}_{+}}xu(x,t)dx<\infty$,
there exists a unique solution to equation (1.1). Moreover, for each $t>0$,
$u(\cdot,t)$ is a probability density.

A function $u_{\rm stat}\in C^{2}({\mathbb{R}}_{+})$ is said to be a
stationary solution associated with (1.1), if it is a probability density
possessing a finite second moment, and, setting $\xi_{0}(dx)=u_{\rm
stat}(x)dx$, the solution to (1.1) is given by $u(x,t)=u_{\rm stat}(x)$ for
all $x,t$.

######  Proposition 2.2.

Assume $\rho<0$. Then there exists a unique stationary solution to (1.1). It
is given by (1.2).

Theorem 2.1 is a consequence of a result stated next, that gives uniqueness of
weak solutions to a class of viscous scalar conservation laws, and provides
one of the main tools used in this paper. This is a parabolic equation of the
type

(2.16)
$\begin{cases}v_{t}=\left(\mathfrak{f}(v)\right)_{x}+av_{xx},&(x,t)\in{\mathbb{R}}_{+}^{2},\\\
v(0,t)=1,&t>0,\\\ v(\cdot,0)=v_{0}.\end{cases}$

A key point is that this result does not require any regularity assumption in
the $x$ variable on the class of solutions. The definition of a weak solution
is as follows.

######  Definition 2.3.

A function $v\in\mathbb{L}_{\rm
loc}^{\infty}({\mathbb{R}}_{+};\mathbb{L}^{\infty}({\mathbb{R}}_{+}))\cap\mathbb{L}_{\rm
loc}^{\infty}({\mathbb{R}}_{+};\mathbb{L}^{1}({\mathbb{R}}_{+}))$ is a weak
solution of (2.16) if for any $t\in(0,\infty)$ and any $\phi\in
C^{\infty}_{0}({\mathbb{R}}_{+})$ satisfying $\phi(0)=0$,

(2.17) $\displaystyle\langle v(\cdot,t),\phi\rangle-\langle
v_{0},\phi\rangle=-\int_{0}^{t}\langle\mathfrak{f}(v(\cdot,s)),\phi^{\prime}\rangle
ds+a\int_{0}^{t}\langle v(\cdot,s),\phi^{\prime\prime}\rangle
ds+a\phi^{\prime}(0)t.$

The choice of $\mathfrak{f}$ that will be of interest is
$\mathfrak{f}(z)=c_{1}z-bz^{\ell}$. In particular, (1.1) is related (2.16)
when the latter takes the special form

(2.18)
$\begin{cases}v_{t}=(c_{1}v-bv^{\ell})_{x}+av_{xx},&(x,t)\in{\mathbb{R}}_{+}^{2},\\\
v(0,t)=1,&t>0,\\\ v(\cdot,0)=\xi_{0}(\cdot,\infty).\end{cases}$

The existence and uniqueness of a classical solution of (2.16) is well known,
and there is a vast literature on weak solutions in the $W^{1,p}_{\rm loc}$
sense. However, we found no uniqueness result of a weak solution as given in
Definition 2.3 under the mere assumption $v\in\mathbb{L}_{\rm
loc}^{\infty}({\mathbb{R}}_{+};\mathbb{L}^{\infty}({\mathbb{R}}_{+}))\cap\mathbb{L}_{\rm
loc}^{\infty}({\mathbb{R}}_{+};\mathbb{L}^{1}({\mathbb{R}}_{+}))$.

######  Theorem 2.4.

Assume $\mathfrak{f}\in C^{\infty}({\mathbb{R}})$ and
$v_{0}\in\mathbb{L}^{1}({\mathbb{R}}_{+})\cap\mathbb{L}^{\infty}({\mathbb{R}}_{+})$.
Then there exists a unique weak solution to (2.16) (in particular, to (2.18))
in the sense of Definition 2.3. This solution is in
$C^{\infty}({\mathbb{R}}_{+}\times(0,\infty))$. If $v$ denotes the weak
solution to (2.18) then $u=-v_{x}$ is the classical solution of (1.1).

Next is a hydrodynamic limit result.

######  Theorem 2.5.

Let $\xi_{0}$ be extended to a trajectory
$\xi=\\{\xi_{t},t\in{\mathbb{R}}_{+}\\}$ in
${\mathcal{M}}_{1}({\mathbb{R}}_{+})$ by setting

$\xi_{t}(dx)=u(x,t)dx,\qquad t>0,$

where $u$ is the unique solution to (1.1). Then $\xi\in
C({\mathbb{R}}_{+},{\mathcal{M}}_{1}({\mathbb{R}}_{+}))$, and one has
$\bar{\xi}^{n}\to\xi$ in probability in
$D({\mathbb{R}}_{+},{\mathcal{M}}_{1}({\mathbb{R}}_{+}))$, as $n\to\infty$.

We next provide two versions of an invariance principle under two different
assumptions on the initial conditions.

######  Theorem 2.6.

Fix $k\in{\mathbb{N}}$ and assume that
$(\hat{X}^{n}_{i}(0))_{i\in[k]}\Rightarrow(X_{i}(0))_{i\in[k]}$. Let $v$ be
the weak solution to (2.18). Then
$(\hat{X}^{n}_{i},\hat{L}^{n}_{i})_{i\in[k]}\Rightarrow(X_{i},L_{i})_{i\in[k]}$
in $(D({\mathbb{R}}_{+},{\mathbb{R}}_{+})\times
C({\mathbb{R}}_{+},{\mathbb{R}}_{+}))^{k}$, where the latter tuple is the
unique in law solution to the system (1.3), driven by a $k$-dimensional
standard BM $(W_{i})_{i\in[n]}$ independent of $(X_{i}(0))_{i\in[k]}$.

######  Theorem 2.7.

Assume that, for each $n$, the ${\mathbb{R}}_{+}$-valued random variables
$X^{n}_{i}(0-)$, $i\in[n]$, are exchangeable (necessitating that the limit law
of each $X^{n}_{i}(0-)$ is $\xi_{0}$). Let $v$ be the weak solution to (2.18).
Then for $k\in{\mathbb{N}}$,
$(\hat{X}^{n}_{i},\hat{L}^{n}_{i})_{i\in[k]}\Rightarrow(X_{i},L_{i})_{i\in[k]}$
in $(D({\mathbb{R}}_{+},{\mathbb{R}}_{+})\times
C({\mathbb{R}}_{+},{\mathbb{R}}_{+}))^{k}$, where the latter tuple is given by
$k$ independent copies of the unique in law solution $(X,L)$ to (1.4) driven
by a standard BM $W$ independent of $X(0)$, and the latter is distributed
according to $\xi_{0}$.

In view of Theorem 2.5, the SDE in Theorem 2.7 is a McKean-Vlasov SDE, as it
can be seen to take the form (1.5). Note that we do not assume that the entire
initial data is exchangeable, nor that the tie breaking rule is symmetric
w.r.t. exchanging $i$’s. The result implies that these issues have negligible
effect upon taking the limit.

###  2.3. Discussion

The hydrodynamic limit result establishes that the dynamics of the model are
given at the macroscopic scale by the PDE (1.1). One can use this to express
basic performance criteria in terms of the solution $u$. In particular,
denoting the macroscopic mean by

$m_{\rm mac}(t):=\int_{0}^{\infty}xu(x,t)dx,$

a first order approximation to the mean queue length is given by $m_{\rm
mac}(t)n^{1/2}$. Accordingly, the mean delay is, to first order, given by
$\lambda^{-1}m_{\rm mac}(t)n^{-1/2}$. Similarly, letting

$\sigma_{\rm mac}^{2}(t):=\int_{0}^{\infty}(x-m_{\rm mac}(t))^{2}u(x,t)dx,$

the macroscopic standard deviation, $\sigma_{\rm mac}(t)$, may be taken as an
index of balance. As a second order parabolic equation in one spatial
variable, there exist standard tools for numerically solving (1.1). Some
examples of solutions are plotted in Fig. 4.

Whereas the macroscopic dynamics can be solved numerically, the macroscopic
equilibrium state has an explicit formula, namely (1.2). One would like to
show rigorously that all solutions $u$ converge to $u_{\rm stat}$ as
$t\to\infty$, and further, that the invariant distribution of the stochastic
model’s dynamics converges to $u_{\rm stat}$ as $n\to\infty$; these will be
addressed in the future. However, even without establishing these limit
results, $u_{\rm stat}$ constitutes a legitimate solution to the PDE. In
particular, a combination of Theorem 2.5 and Proposition 2.2 implies that
$\bar{\xi}^{n}\to\xi$ in probability holds provided that $\xi_{0}(dx)=u_{\rm
stat}(x)dx$. Here, $\xi_{t}=\xi_{0}$ for all $t$. Fig. 4 shows graphs of
$u_{\rm stat}$ for different values of $\ell$ and $b$, based on (1.2).

We can also use $u_{\rm stat}$ in place of $u(\cdot,t)$ in the above
macroscopic performance criteria, and define analogously $m_{\rm mac}$ and
$\sigma_{\rm mac}$ which are now time independent quantities. Fig. 4 shows
their dependence on $\ell$ and $b$.

Figure 2. Solution of (1.1) for initial condition $\rm unif[0,10]$ (left) and
$\rm unif[15,30]$ (right) at different time instances, with $\rho=-.01$,
$c_{1}=.21$, $b=.2$, $\ell=4$, $a=1$.

Figure 3. Stationary solution for different values of $\ell$, with $b=.2$
(left), and for different values of $b$, with $\ell=4$ (right). In all cases,
$\rho=-.01,a=1,\lambda=1$.

Figure 4. Stationary mean and standard deviation for $2\leq\ell\leq 9$, with
$b=2$, $a=1$ (left) and for $0.1\leq b\leq 2$, with $\ell=4$, $a=1$ (right).

The nature of the drift in equation (1.1) is akin to the dynamics of a mass
distribution subjected to a gravitational field [44], where the force acting
on an infinitesimal mass depends on the amount of mass that piles above it.
Here it is the mass accumulated below it, specifically, a nonlinear function
of $v(x,t)$, that determines the force.

###  2.4. Outline of the proof

The approach used to proving the hydrodynamic limit result is based on PDE
uniqueness. Thus the first main tool required is the uniqueness of weak
solutions to (2.16). The crux of the argument is as follows. If $v_{1},v_{2}$
are two weak solutions then one can show that $w=v_{1}-v_{2}$ satisfies the
equality

$\langle
w(\cdot,t),\psi(\cdot,t)\rangle=\int_{0}^{t}\Big{\langle}w(\cdot,s),\Big{[}\psi_{s}(\cdot,s)-\Big{[}\frac{\mathfrak{f}(v_{1})-\mathfrak{f}(v_{2})}{v_{1}-v_{2}}\Big{]}\psi_{x}+\psi_{xx}\Big{]}\Big{\rangle}ds$

for any $t>0$ and test function $\psi$ satisfying $\psi\in
C_{0}^{\infty}(\mathbb{R}_{+}\times[0,t])$, $\psi(0,s)=0$ for $s\in[0,t]$. One
can then find an $\mathbb{L}^{\infty}$ approximation of the solution to the
linear, backward equation

$\psi_{s}(\cdot,s)-\Big{[}\frac{f(v_{1})-f(v_{2})}{v_{1}-v_{2}}\Big{]}\psi_{x}+\psi_{xx}=0$

on $\mathbb{R}_{+}\times[0,t]$ satisfying $\psi(\cdot,t)=\phi(\cdot)$, given
any smooth $\phi$. As a result, $\langle w(\cdot,t),\phi\rangle=0$ for any $t$
and any such $\phi$, hence $w\equiv 0$.

Next, $\bar{\xi}^{n}$ is shown to form a tight sequence, and then the next
main step is to prove that limit points satisfy (2.16) in weak sense.
Achieving this goal relies on two key elements. One is a semimartingale
decomposition for point processes, introduced in [17]. This decomposition is
particularly convenient in our setting, where the number of point processes
involved, given by the departure processes, grows to infinity. The other is an
estimate showing $C$-tightness of $\hat{X}^{n}_{i}$ uniformly in $i\in[n]$.
With this toolbox we can then show the following. Let $\tilde{\phi}$ be a test
function as in Definition 2.3 and $\phi$ its antiderivative. Then every limit
point $\xi$ of $\bar{\xi}^{n}$ satisfies

(2.19) $\displaystyle\langle\phi,\xi_{t}\rangle$
$\displaystyle=\langle\phi,\xi_{0}\rangle+\int_{0}^{t}\langle
b_{1}\phi^{\prime}+a\phi^{\prime\prime},\xi_{s}\rangle
ds+\frac{b_{0}}{\ell}\int_{0}^{t}\int_{{\mathbb{R}}_{+}}\phi^{\prime}(x)\mathfrak{S}(\xi_{s}[x,\infty),\xi_{s}(x,\infty))\xi_{s}(dx)ds,$

where

(2.20)
$\mathfrak{S}(a,b)=a^{\ell-1}+a^{\ell-2}b+\cdots+ab^{\ell-2}+b^{\ell-1},\qquad
a,b\in{\mathbb{R}}.$

To see the relation to equation (2.17) required by Definition 2.3, let
$v(x,t)=\xi_{t}(x,\infty)$. Then
$\langle\phi,\xi_{t}\rangle=\int_{{\mathbb{R}}_{+}}(v(x,t)-1)\tilde{\phi}(x)dx$.
Suppose that $\xi_{t}$ has no atoms for every $t>0$. Then, owing to the fact
that $\mathfrak{S}(a,a)=\ell a^{\ell-1}$, (2.19) reduces to (2.17) in the case
of interest $\mathfrak{f}(z)=c_{1}z-bz^{\ell}$. However, an a priori proof of
the atomless property is not required, as a calculation via integration by
parts shows that (2.19) and (2.17) are equivalent even in presence of atoms.
As a result, the existence of weak limit of $\bar{\xi}^{n}$ satisfying the PDE
follows, completing the proof of the hydrodynamic limit.

Equation (2.16) is well known to have $C^{\infty}$ solutions, and as a
consequence of the above result the relation of the hydrodynamic limit to the
PDE (1.1) for the density follows.

Finally, the two invariance principles follow by combining the $C$-tightness
of individual rescaled queue lengths with the hydrodynamic limit result, where
the atomless property, that by now has been proved, gives a control over the
interaction term.

The structure of §§3–4 is as follows. Preliminary steps are taken in §3,
starting from §3.1 where the uniqueness of weak solutions to the PDE (2.16) is
proved. The aforementioned semimartingale representations are provided in
§3.2. §3.3 gives uniform estimates on rescaled processes, including
$C$-tightness and second moment. Tightness of $\bar{\xi}^{n}$ is shown in
§3.4. With this set of tools, the proof is then carried out in §4, which
starts by writing down an equation for $\langle\phi,\bar{\xi}^{n}\rangle$. In
§4.1 it is shown that the linear term in this equation converges to the linear
term in (2.19), and in §4.2, that the interaction term converges to that in
(2.19). The reduction, alluded to above, of (2.19) to (2.17), is shown in
§4.3. Finally, all the main results are proved, based on the above, in §4.4.

##  3\. Preliminaries

###  3.1. PDE uniqueness

Here we prove the uniqueness part of Theorem 2.4, stated as Lemma 3.2 below.
In preparation for proving this result, we provide an extension of the class
of test functions allowed in Definition 2.3, which is relatively standard but
is given for completeness.

######  Lemma 3.1.

If $v$ is a weak solution of (2.16) (in the sense of Definition 2.3) then for
any $t\in(0,\infty)$ and any $\psi\in
C_{0}^{\infty}\left({\mathbb{R}}_{+}\times[0,t]\right)$ satisfying
$\psi(0,s)=0$, $0\leq s\leq t$, one has

$\displaystyle\langle v(\cdot,t),\psi(\cdot,t)\rangle-\langle
v_{0},\psi(\cdot,0)\rangle$ $\displaystyle=\int_{0}^{t}\langle
v(\cdot,s),\psi_{s}(\cdot,s)\rangle
ds-\int_{0}^{t}\langle\mathfrak{f}(v(\cdot,s)),\psi_{x}\rangle
ds+a\int_{0}^{t}\langle v(\cdot,s),\psi_{xx}\rangle ds$ (3.1)
$\displaystyle\qquad+a\int_{0}^{t}\psi_{x}(0,s)ds.$

Proof. Observe by (2.17) that for any $\phi$ satisfying Definition 2.3 and for
any $0\leq\tau_{1}<\tau_{2}\leq t$,

(3.2) $\langle v(\cdot,\tau_{2}),\phi\rangle=\langle
v(\cdot,\tau_{1}),\phi\rangle-\int_{\tau_{1}}^{\tau_{2}}\langle\mathfrak{f}(v(\cdot,s)),\phi^{\prime}\rangle
ds+a\int_{\tau_{1}}^{\tau_{2}}\left(\langle
v(\cdot,s),\phi^{\prime\prime}\rangle+\phi^{\prime}(0)\right)ds.$

Since $v$ and $\mathfrak{f}(v)$ are in $\mathbb{L}^{\infty}({\mathbb{R}}_{+})$
uniformly in $t$, it follows, in particular, that $t\mapsto\langle
v(\cdot,t),\phi\rangle$ is continuous and, moreover,

(3.3) $|\langle
v(\cdot,\tau_{2})-v(\cdot,\tau_{1}),\phi\rangle|\leq|\tau_{1}-\tau_{2}|K$

where $K=K(\phi,v)$.

Let now $\psi\in C_{0}^{\infty}({\mathbb{R}}_{+}\times[0,t])$ be a function
that satisfies the conditions of the lemma. Let $N\in\mathbb{N}$. Let
$t_{k}:=tk/N$ where $k=0,\ldots,N$. Define $\phi_{k}(x)=\psi(t_{k},x)$. Then,
from (3.2)

$\langle v(\cdot,t_{k+1}),\phi_{k}\rangle-\langle
v(\cdot,t_{k}),\phi_{k}\rangle=-\int_{t_{k}}^{t_{k+1}}\langle\mathfrak{f}(v(\cdot,s)),\phi^{\prime}_{k}\rangle
ds+a\int_{t_{k}}^{t_{k+1}}\left(\langle
v(\cdot,s),\phi^{\prime\prime}_{k}\rangle+\phi^{\prime}_{k}(0)\right)ds.$

Summing over $k=0,\ldots,N-1$ and using $\phi_{N}=\psi(\cdot,t)$,
$\phi_{0}=\psi(\cdot,0)$ we get

(3.4) $\langle v(\cdot,t),\psi(\cdot,t(1-1/N))\rangle-\langle
v_{0},\psi(\cdot,0)\rangle-\sum_{k=1}^{N-1}\langle
v(\cdot,t_{k+1}),\phi_{k+1}-\phi_{k}\rangle\\\
=\sum_{k=0}^{N-1}\left[-\int_{t_{k}}^{t_{k+1}}\langle\mathfrak{f}(v(\cdot,s)),\phi^{\prime}_{k}\rangle
ds+a\int_{t_{k}}^{t_{k+1}}\left(\langle
v(\cdot,s),\phi^{\prime\prime}_{k}\rangle+\phi^{\prime}_{k}(0)\right)ds\right].$

Evidently, $\langle v(\cdot,t),\psi(\cdot,t(1-1/N))\rangle\rightarrow\langle
v(\cdot,t),\psi(\cdot,t)\rangle$ as $N\rightarrow\infty$. Next, by (3.3)

$\langle
v(\cdot,t_{k+1}),\phi_{k+1}-\phi_{k}\rangle=\int_{t_{k}}^{t_{k+1}}\langle
v(\cdot,s),\psi_{s}(\cdot,s)\rangle ds+\int_{t_{k}}^{t_{k+1}}\langle
v(\cdot,t_{k+1})-v(\cdot,s),\psi_{s}(\cdot,s)\rangle ds$
$=\int_{t_{k}}^{t_{k+1}}\langle v(\cdot,s),\psi_{s}(\cdot,s)\rangle
ds+O(N^{-2}).$

Thus

$\sum_{k=1}^{N-1}\langle
v(\cdot,t_{k+1}),\phi_{k+1}-\phi_{k}\rangle\rightarrow\int_{0}^{t}\langle
v(\cdot,s),\psi_{s}(\cdot,s)\rangle ds$

as $N\rightarrow\infty$. Likewise, the right side of (3.4) converges to

$-\int_{0}^{t}\langle\mathfrak{f}(v(\cdot,s)),\psi_{x}\rangle
ds+a\int_{0}^{t}\left(\langle
v(\cdot,s),\psi_{xx}\rangle+\psi_{x}(0,s)\right)\rangle ds.$

This proves the lemma. ∎

######  Lemma 3.2.

Assume $\mathfrak{f}\in C^{\infty}({\mathbb{R}})$ and
$v_{0}\in\mathbb{L}^{1}({\mathbb{R}}_{+})\cap\mathbb{L}^{\infty}({\mathbb{R}}_{+})$.
Then there is at most one weak solution to (2.16) in the sense of Definition
2.3.

Proof. As can be seen by performing a change of variables $x\mapsto
a^{-1/2}x$, we may and will assume w.l.o.g. that $a=1$. Suppose $v_{1},v_{2}$
are two solutions. The goal is to show that $v_{1}=v_{2}$. Let $w=v_{1}-v_{2}$
and note that $w\in\mathbb{L}^{\infty}_{\rm
loc}({\mathbb{R}}_{+},\mathbb{L}^{\infty}({\mathbb{R}}_{+}))\cap\mathbb{L}^{\infty}_{\rm
loc}({\mathbb{R}}_{+},\mathbb{L}^{1}({\mathbb{R}}_{+}))$. Let

$J(x,t):=\begin{cases}\displaystyle\frac{\mathfrak{f}(v_{1}(x,t))-\mathfrak{f}(v_{2}(x,t))}{v_{1}(x,t)-v_{2}(x,t)}&\text{if}\
\ w(x,t)\not=0\\\ \mathfrak{f}^{\prime}(v_{1}(x,t))&\text{if}\ \
w(x,t)=0.\end{cases}$

Then for any test function $\psi$ satisfying the conditions of Lemma 3.1,

(3.5) $\langle w(\cdot,t),\psi(\cdot,t)\rangle=\int_{0}^{t}\langle
w(\cdot,s),\left[\psi_{s}(\cdot,s)-J(\cdot,s)\psi_{x}+\psi_{xx}\right]\rangle
ds.$

Fix $T$. By the assumptions on $\mathfrak{f}$ and the definition of a solution
we know
$J\in\mathbb{L}^{\infty}([0,T];\mathbb{L}^{\infty}({\mathbb{R}}_{+}))$.
Moreover, for any $x_{0}>0$,
$|\mathfrak{f}^{\prime}(x)-\mathfrak{f}^{\prime}(0)|\leq cx$ holds provided
that $|x|<x_{0}$, with $c=c(x_{0})$. Hence
$J-\mathfrak{f}^{\prime}(0)\in\mathbb{L}^{\infty}([0,T];\mathbb{L}^{1}({\mathbb{R}}_{+}))$.
Let now $J_{N}\in C^{\infty}({\mathbb{R}}_{+}\times[0,T])$ be a sequence that
is bounded uniformly in $\mathbb{L}^{\infty}({\mathbb{R}}_{+}\times[0,T])$ and
satisfies $\|J_{N}-J\|_{1}=O(N^{-1})$, where
$\|\cdot\|_{1}=\|\cdot\|_{\mathbb{L}^{1}({\mathbb{R}}_{+}\times[0,T])}$. Given
$t\in(0,T]$, let $\tilde{\psi}^{N}$ be the classical solution of the backward
linear problem on the time interval $[0,t]$:

(3.6)
$\tilde{\psi}^{N}_{s}(\cdot,s)-J_{N}(\cdot,s)\tilde{\psi}^{N}_{x}+\tilde{\psi}^{N}_{xx}=0,\
\ \tilde{\psi}^{N}(\cdot,t)=\phi,\ \ \tilde{\psi}^{N}(0,s)=0,\qquad 0\leq
s<t,$

where $\phi\in C^{\infty}_{0}({\mathbb{R}}_{+})$. Since $\tilde{\psi}^{N}$ and
all its derivatives decay to zero as $x\rightarrow\infty$, uniformly in
$s\in[0,t]$, we may replace them by $\psi^{N}$ which satisfy the conditions of
Lemma 3.1 and

$\|\psi^{N}-\tilde{\psi}^{N}\|_{\infty}+\|\psi^{N}_{x}-\tilde{\psi}^{N}_{x}\|_{\infty}+\|\psi^{N}_{xx}-\tilde{\psi}^{N}_{xx}\|_{\infty}+\|\psi^{N}_{t}-\tilde{\psi}^{N}_{t}\|_{\infty}=O(N^{-1})$

uniformly on $[0,t]$. It follows from (3.5) that

$\langle w(\cdot,t),\phi\rangle=\int_{0}^{t}\langle
w(\cdot,s),[\tilde{\psi}^{N}_{s}(\cdot,s)-J(\cdot,s)\tilde{\psi}^{N}_{x}+\tilde{\psi}^{N}_{xx}]\rangle
ds\ =\int_{0}^{t}\langle
w(\cdot,s),\left[J_{N}(\cdot,s)-J\right]\tilde{\psi}^{N}_{x}\rangle ds$
$=\int_{0}^{t}\langle
w(\cdot,s),\left[J_{N}(\cdot,s)-J\right]\psi^{N}_{x}\rangle ds+O(N^{-1})$

for any such $\phi$. Suppose

(3.7)
$\sup_{N}\sup_{s\in[0,t]}\|\tilde{\psi}^{N}_{x}(\cdot,s)\|_{\infty}<\infty$

if $t$ is small enough. Since
$w\in\mathbb{L}^{\infty}([0,T];\mathbb{L}^{\infty}({\mathbb{R}}_{+}))$ and
$\|J_{N}-J\|_{1}=O(N^{-1})$, this implies $\langle
w(\cdot,t),\phi\rangle=O(N^{-1})$ for any $\phi\in
C^{\infty}_{0}({\mathbb{R}}_{+})$ and any $N$, hence $w(\cdot,t)=0$. As a
consequence, $w(\cdot,t)=0$ for all $t\in[0,t_{0}]$, some $t_{0}>0$. Moreover,
$t_{0}$ does not depend on the initial condition $v_{0}$. Thus, iterating the
argument shows that $w=0$.

It thus suffices to show (3.7). To this end, denote
$M:=\sup_{N}\|J_{N}\|_{\infty}<\infty$. Let $\psi^{0}$ be the solution of

$\psi^{0}_{s}(\cdot,s)+\psi^{0}_{xx}=0,\ \ \psi^{0}(\cdot,t)=\phi,\ \
\psi^{0}(0,s)=0,\ \ 0\leq s\leq t.$

Let $m_{N}(y,\tau):=J_{N}(y,\tau)\tilde{\psi}^{N}_{y}(y,\tau)$. Since, as
mentioned earlier, $\|\tilde{\psi}^{N}_{y}\|_{\infty}<\infty$ for all $N$, one
also has $\|m_{N}\|_{\infty}<\infty$ for all $N$. Moreover,
$\tilde{\psi}^{N}=\psi_{0}+\hat{\psi}^{N}$ where $\hat{\psi}^{N}$ is the
solution of

$\hat{\psi}^{N}_{s}+\hat{\psi}^{N}_{xx}=m_{N},\ \ \
\hat{\psi}^{N}(\cdot,t)=0,\ \hat{\psi}^{N}(0,s)=0,\qquad 0\leq s\leq t.$

The solution $\hat{\psi}^{N}$ is given by Duhamel’s principle:

$\hat{\psi}^{N}(x,s)=\frac{1}{2\sqrt{\pi}}\int_{s}^{t}d\tau(t-\tau)^{-1/2}\left[\int_{0}^{\infty}m_{N}(y,t-\tau)e^{-\frac{(x-y)^{2}}{4(t-\tau)}}dy-\int_{0}^{\infty}m_{N}(y,t-\tau)e^{-\frac{(x+y)^{2}}{4(t-\tau)}}dy\right].$

Then $\hat{\psi}^{N}_{x}(x,s)=-A_{N}(x,s)+B_{N}(x,s)$, where

$\displaystyle A_{N}(x,s)$
$\displaystyle=\frac{1}{4\sqrt{\pi}}\int_{s}^{t}d\tau(t-\tau)^{-3/2}\int_{0}^{\infty}m_{N}(y,t-\tau)(x-y)e^{-\frac{(x-y)^{2}}{4(t-\tau)}}dy,$
$\displaystyle B_{N}(x,s)$
$\displaystyle=\frac{1}{4\sqrt{\pi}}\int_{s}^{t}d\tau(t-\tau)^{-3/2}\int_{0}^{\infty}m_{N}(y,t-\tau)(x+y)e^{-\frac{(x+y)^{2}}{4(t-\tau)}}dy.$

Both $|A_{N}(x,s)|$ and $|B_{N}(x,s)|$ are bounded by

$\frac{\|m_{N}\|_{\infty}}{4\sqrt{\pi}}\int_{s}^{t}d\tau(t-\tau)^{-3/2}\int_{-\infty}^{\infty}|x+y|e^{-\frac{(x+y)^{2}}{4(t-\tau)}}dy.$

Changing variables of integration,

$|A_{N}(x,s)|\vee|B_{N}(x,s)|\leq\frac{\|m_{N}\|_{\infty}}{2\sqrt{\pi}}\int_{s}^{t}d\tau(t-\tau)^{-1/2}\int_{-\infty}^{\infty}|z|e^{-\frac{z^{2}}{2}}dz\leq
C\|m_{N}\|_{\infty}t^{1/2},$

for a constant $C$. Recall that $\|m_{N}\|_{\infty}\leq
M\|\psi^{N}_{x}\|_{\infty}$ and $\tilde{\psi}^{N}=\psi^{0}+\hat{\psi}^{N}$.
Thus

$\|\tilde{\psi}^{N}_{x}(\cdot,s)\|_{\infty}\leq\|\psi^{0}_{x}\|_{\infty}+2C\|m_{N}\|_{\infty}t^{1/2}\leq\|\psi^{0}_{x}\|_{\infty}+2CM\|\psi^{N}_{x}(\cdot,s)\|_{\infty}t^{1/2}\
,\qquad 0<s\leq t.$

Hence $\|\tilde{\psi}^{N}_{x}(\cdot,s)\|_{\infty}\leq
2\|\psi^{0}_{x}\|_{\infty}$ for $0\leq s\leq(2CM)^{-2}/2$, completing the
proof that $v_{1}=v_{2}$. ∎

###  3.2. Martingale toolbox

In addition to the delayed renewal processes $S^{n}_{i}$ it will also be
useful to introduce the corresponding non-delayed renewal processes

$S^{0,n}_{i}(t)=\max\Big{\\{}k\in{\mathbb{Z}}_{+}:\sum_{j=1}^{k-1}Z^{n}_{i}(j)\leq
t\Big{\\}},\qquad t\geq 0.$

Note that $S^{0,n}_{i}(0)=1$ and that these processes are IID (unlike
$S^{n}_{i}$). Additional rescaled processes are denoted as follows

(3.8)
$\hat{E}^{n}_{i}(t)=n^{-1/2}(E^{n}_{i}(t)-\lambda^{n}t),\qquad\hat{S}^{n}_{i}(t)=n^{-1/2}(S^{n}_{i}(t)-\mu^{n}t),\qquad\hat{S}^{0,n}_{i}(t)=n^{-1/2}(S^{0,n}_{i}(t)-\mu^{n}t),$
(3.9)
$\hat{A}^{n}_{0}(t)=n^{-1/2}(A^{n}_{0}(t)-\lambda^{n}_{0}t),\qquad\hat{A}^{n}_{i}(t)=n^{-1/2}A^{n}_{i}(t).$

A semimartingale decomposition of counting processes was introduced in [17],
which deviates from the classical Doob-Meyer decomposition and is convenient
for our purposes. For the renewal process $S^{n}_{i}$, this decomposition is
constructed as follows. Denote by

$R^{n}_{i}(t)=\inf\\{s>0:S^{n}_{i}(t+s)>S^{n}_{i}(t)\\},\qquad t\geq 0,$

the residual time to the next counting instant (note that it is right-
continuous, hence at a time of counting it already shows the time until the
next counting). One has

(3.10) $t+R^{n}_{i}(t)=\sum_{k=0}^{S^{n}_{i}(t)}Z^{n}_{i}(k).$

Hence with

$M^{{\rm
ser},n}_{i}(t)=\sum_{k=1}^{S^{n}_{i}(t)}\zeta^{n}_{i}(k),\qquad\zeta^{n}_{i}(k)=1-\mu^{n}Z^{n}_{i}(k),\qquad
k\geq 1,$

and $\sum_{1}^{0}=0$, one has the Daley-Miyazawa semimartingale
representation,

(3.11) $S^{n}_{i}(t)=\mu^{n}(t-Z^{n}_{i}(0)+R^{n}_{i}(t))+M^{{\rm
ser},n}_{i}(t),\qquad t\geq 0.$

Denoting ${\mathcal{F}}^{{\rm
ser},n}_{i}(t)=\sigma\\{S^{n}_{i}(s),R^{n}_{i}(s),s\leq t\\}$, the first term
on the right is predictable on the filtration $\\{{\mathcal{F}}^{{\rm
ser},n}_{i}(t)\\}$ while $M^{{\rm ser},n}_{i}$ a martingale on it. A form of
this decomposition that will be useful here is based on a filtration
${\mathcal{F}}^{n}_{t}$ that, for each $t$, contains all relevant information
about the system by time $t$. It will be obtained once a time change
transformation of (3.11) is performed, to represent
$D^{n}_{i}(t)=S^{n}_{i}(T^{n}_{i}(t))$. Thus let

$\tilde{R}^{n}_{i}(t)=n^{-1/2}\mu^{n}R^{n}_{i}(t)=\tilde{\mu}^{n}R^{n}_{i}(t),\qquad\hat{M}^{{\rm
ser},n}_{i}(t)=n^{-1/2}M^{{\rm ser},n}_{i}(t),$

and

(3.12) $M^{{\rm
dep},n}_{i}(t)=\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k),\qquad\hat{M}^{{\rm
dep},n}_{i}(t)=n^{-1/2}M^{{\rm dep},n}_{i}(t).$

Then, by (3.11),

(3.13)
$D^{n}_{i}(t)=\mu^{n}(T^{n}_{i}(t)-Z^{n}_{i}(0)+R^{n}_{i}(T^{n}_{i}(t)))+M^{{\rm
dep},n}_{i}(t),$

and

(3.14)
$\hat{S}^{n}_{i}(T^{n}_{i}(t))=-\tilde{Z}^{n}_{i}(0)+\tilde{R}^{n}_{i}(T^{n}_{i}(t))+\hat{M}^{{\rm
dep},n}_{i}(t).$

The tuple

${\mathcal{S}}^{n}(t)=(E^{n}_{i}(t),A^{n}_{i}(t),D^{n}_{i}(t),X^{n}_{i}(t),T^{n}_{i}(t),R^{n}_{i}(T^{n}_{i}(t)),\,i\in[n],\,A^{n}_{0}(t),\theta^{n}_{A^{n}_{0}(t)})$

is referred to as the state of the system at time $t$. Denote the
corresponding filtration by

${\mathcal{F}}^{n}_{t}=\sigma\\{{\mathcal{I}}^{n},{\mathcal{S}}^{n}(s),s\in[0,t]\\}.$

Note that at times $t$ when server $i$ is active, $R^{n}_{i}(T^{n}_{i}(t))$ is
the residual time till the completion of service of the job being processed by
that server, whereas at times when the server is idle,
$R^{n}_{i}(T^{n}_{i}(t))$ gives the service duration of the job that will be
processed next by this server.

######  Lemma 3.3.

i. The processes $\hat{M}^{{\rm dep},n}_{i}$, $\hat{E}^{n}_{i}$, and
$\hat{A}^{n}_{0}$ are $\\{{\mathcal{F}}^{n}_{t}\\}$-martingales, with optional
quadratic variations given by

$[\hat{M}^{{\rm
dep},n}_{i}](t)=n^{-1}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)^{2},\qquad[\hat{E}^{n}_{i}](t)=n^{-1}E^{n}_{i}(t),\qquad[\hat{A}^{n}_{0}](t)=n^{-1}A^{n}_{0}(t),$

and one has

(3.15) ${\mathbb{E}}\\{[\hat{M}^{{\rm dep},n}_{i}](t)\\}=n^{-1}\sigma_{\rm
ser}^{2}{\mathbb{E}}\\{D^{n}_{i}(t)\\}<\infty.$

ii. For distinct $i,j\in[n]$,

(3.16) ${\mathbb{E}}\\{[\hat{M}^{{\rm dep},n}_{i},\hat{M}^{{\rm
dep},n}_{j}](t)\\}=0.$

iii. The process

$\hat{M}^{A,n}_{i}(t)=\hat{A}^{n}_{i}(t)-\hat{C}^{A,n}_{i}(t),\quad\text{
where
}\quad\hat{C}^{A,n}_{i}(t)=\lambda^{n}_{0}n^{-1/2}\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s)}ds,$

is an $\\{{\mathcal{F}}^{n}_{t}\\}$-martingale, whose optional quadratic
variation is $[\hat{M}^{A,n}_{i}](t)=n^{-1}A^{n}_{i}(t)$.

Proof. i. For adaptedness of $\hat{M}^{{\rm dep},n}_{i}$ it suffices to prove
that $Z^{n}_{i}(k)1_{k\leq D^{n}_{i}(t)}\in{\mathcal{F}}^{n}_{t}$ for all $k$.
This is shown as follows. By (3.10),
$T^{n}_{i}(s)+R^{n}_{i}(T^{n}_{i}(s))=\sum_{k=0}^{D^{n}_{i}(s)}Z^{n}_{i}(k)$.
Hence $\\{Z^{n}_{i}(k),k\leq D^{n}_{i}(t)\\}$ can all be recovered from the
tuple $\\{T^{n}_{i}(s),R^{n}_{i}(T^{n}_{i}(s)),D^{n}_{i}(s)\\}$, as $s$ varies
between $0$ and $t$. Since the latter is
$\\{{\mathcal{F}}^{n}_{t}\\}$-adapted, this proves the claim.

Next it is shown that $\hat{M}^{{\rm dep},n}_{i}(t)\in L_{1}(d{\mathbb{P}})$.
Note first that as a renewal process, $S^{n}_{i}(t)$ has finite expectation
for every $t$. Since $T^{n}_{i}(t)\leq t$ this gives
${\mathbb{E}}[D^{n}_{i}(t)]<\infty$. Let

$t^{n}_{i}(k)=\inf\\{t\geq 0:D^{n}_{i}(t)\geq k\\},\qquad k=1,2,\ldots.$

These are clearly stopping times on $\\{{\mathcal{F}}^{n}_{t}\\}$. Hence

(3.17) $t^{n}_{i}(k)\in{\mathcal{F}}^{n}_{t^{n}_{i}(k)-},\qquad k\geq 1,$

where we recall that for a stopping time $\tau$,

${\mathcal{F}}^{n}_{\tau-}={\mathcal{F}}^{n}_{0}\vee\sigma\\{A\cap\\{\tau<t\\}:A\in{\mathcal{F}}^{n}_{t},t\geq
0\\}$

(see [27, I.1.11 and I.1.14]). The state of the system up to $t^{n}_{i}(k)-$,
namely $\\{{\mathcal{S}}^{n}(t),t<t^{n}_{i}(k)\\}$, can be recovered from the
tuple ${\mathcal{I}}^{n}$, $(E^{n}_{i}(t),t\in{\mathbb{R}}_{+},i\in[n])$,
$(A^{n}_{0}(t),t\in{\mathbb{R}}_{+})$, $(\theta^{n}_{j},j\in{\mathbb{N}})$,
$(Z^{n}_{\ell}(j),j\in{\mathbb{N}})$, $\ell\in[n]\setminus\\{i\\}$ and
$(Z^{n}_{i}(j),j<k)$, as follows by the construction of the model. By our
assumptions, $Z^{n}_{i}(k)$ is independent of this tuple. As a result, it is
independent of ${\mathcal{F}}^{n}_{t^{n}_{i}(k)-}$. It follows from (3.17)
that

$\\{D^{n}_{i}(t)\leq k\\}=\\{t^{n}_{i}(k)\geq
t\\}\in{\mathcal{F}}^{n}_{t^{n}_{i}(k)-}.$

The structure that we have just proved, where

(3.18)
$Z^{n}_{i}(k^{\prime})\in{\mathcal{F}}^{n}_{t^{n}_{i}(k)-},k^{\prime}<k,\text{
whereas $Z^{n}_{i}(k)$ is independent of
${\mathcal{F}}^{n}_{t^{n}_{i}(k)-}$},$

$D^{n}_{i}(t)$ is a stopping time on the discrete parameter filtration
$\\{{\mathcal{F}}^{n}_{t^{n}_{i}(k)-},k\in{\mathbb{N}}\\}$,
${\mathbb{E}}[D^{n}_{i}(t)]<\infty$, along with the fact that $Z^{n}_{i}(k)$
are IID with ${\mathbb{E}}[Z^{n}_{i}(k)]<\infty$, allows us to use Wald’s
identity, showing

${\mathbb{E}}\Big{[}\sum_{k=1}^{D^{n}_{i}(t)}Z^{n}_{i}(k)\Big{]}={\mathbb{E}}\\{D^{n}_{i}(t)\\}{\mathbb{E}}\\{Z^{n}_{i}(1)\\}={\mathbb{E}}\\{D^{n}_{i}(t)\\}(\mu^{n})^{-1}<\infty.$

Since $|\zeta^{n}_{i}(k)|\leq 1+\mu^{n}Z^{n}_{i}(k)$, this shows that
${\mathbb{E}}\\{\hat{M}^{{\rm dep},n}_{i}(t)\\}<\infty$.

To show the martingale property, note that by the independence stated in
(3.18), we have

(3.19)
${\mathbb{E}}[\zeta^{n}_{i}(k)|{\mathcal{F}}^{n}_{t^{n}_{i}(k)-}]=0,\qquad
k\geq 1.$

Arguing now along the lines of the proof of [17, Lemma 2.1],

$\displaystyle\hat{M}^{{\rm dep},n}_{i}(t)$
$\displaystyle=n^{-1/2}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)=n^{-1/2}\sum_{k=1}^{\infty}\zeta^{n}_{i}(k)1_{\\{t^{n}_{i}(k)\leq
t\\}}.$

Hence for $s<t$,

$\displaystyle{\mathbb{E}}[\hat{M}^{{\rm
dep},n}_{i}(t)|{\mathcal{F}}^{n}_{s}]-\hat{M}^{{\rm dep},n}_{i}(s)$
$\displaystyle=n^{-1/2}\sum_{k=1}^{\infty}{\mathbb{E}}[\zeta^{n}_{i}(k)1_{\\{s<t^{n}_{i}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{s}]$
$\displaystyle=n^{-1/2}\sum_{k=1}^{\infty}{\mathbb{E}}[{\mathbb{E}}[\zeta^{n}_{i}(k)1_{\\{s<t^{n}_{i}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{t^{n}_{i}(k)-}]\,|{\mathcal{F}}^{n}_{s}]$
$\displaystyle=0,$

where we used (3.17) and (3.19).

For the martingale property of $\hat{E}^{n}_{i}$ one only needs to show that
$E^{n}_{i}(t)-E^{n}_{i}(s)$ is independent of ${\mathcal{F}}^{n}_{s}$ when
$s<t$. Again, this follows from the fact that all the processes comprising
${\mathcal{S}}^{n}_{u}$, $u\leq s$, can be recovered from the tuple
${\mathcal{I}}^{n}$, $(E^{n}_{i}(u),u\in[0,s],i\in[n])$,
$(A^{n}_{0}(u),u\in[0,s])$, $(\theta^{n}_{k},k\in{\mathbb{N}})$,
$(Z^{n}_{i}(k),k\in{\mathbb{N}},i\in[n])$; but the increment
$E^{n}_{i}(t)-E^{n}_{i}(s)$ is independent of this tuple.

A similar proof holds for $\hat{A}^{n}_{0}$.

The expressions for the quadratic variation are straightforward.

To show (3.15) we can use Wald’s identity as before, now with the IID sequence
$\zeta^{n}_{i}(k)^{2}$, which now gives

${\mathbb{E}}\Big{[}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)^{2}\Big{]}={\mathbb{E}}\\{D^{n}_{i}(t)\\}{\mathbb{E}}\\{\zeta^{n}_{i}(1)^{2}\\},$

and (3.15) follows.

ii. Similar to the argument following (3.17), one can recover the state of the
system up to $t^{n}_{ijkl}-$, where

$t^{n}_{ijkl}=\min(t^{n}_{i}(k),t^{n}_{j}(l)),$

namely $\\{{\mathcal{S}}^{n}(t),t<t^{n}_{ijkl}\\}$, from ${\mathcal{I}}^{n}$,
$(E^{n}_{i}(t),t\in{\mathbb{R}}_{+},i\in[n])$,
$(A^{n}_{0}(t),t\in{\mathbb{R}}_{+})$, $(\theta^{n}_{j},j\in{\mathbb{N}})$,
$(Z^{n}_{i^{\prime}}(k),k\in{\mathbb{N}})$,
$i^{\prime}\in[n]\setminus\\{i,j\\}$, $(Z^{n}_{i}(k^{\prime}),k^{\prime}<k)$
and $(Z^{n}_{j}(l^{\prime}),l^{\prime}<l)$. However, by our assumptions, the
pair $(Z^{n}_{i}(k),Z^{n}_{j}(l))$ is independent of this tuple, hence it is
independent of ${\mathcal{F}}^{n}_{t^{n}_{ijkl}-}$. Since $Z^{n}_{i}(k)$ and
$Z^{n}_{j}(l)$ are mutually independent, this gives

(3.20)
${\mathbb{E}}[\zeta^{n}_{i}(k)\zeta^{n}_{j}(l)|{\mathcal{F}}^{n}_{t^{n}_{ijkl}-}]=0,\qquad
k,l\geq 1.$

We have

$\displaystyle[\hat{M}^{{\rm dep},n}_{i},\hat{M}^{{\rm dep},n}_{j}](t)$
$\displaystyle=n^{-1}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)\sum_{l=1}^{D^{n}_{j}(t)}1_{\\{t^{n}_{j}(l)=t^{n}_{i}(k)\\}}\zeta^{n}_{j}(l)$
$\displaystyle=n^{-1}\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\zeta^{n}_{i}(k)\zeta^{n}_{j}(l)1_{\\{t^{n}_{j}(l)=t^{n}_{i}(k)\leq
t\\}}.$

In view of (3.17), $1_{\\{t^{n}_{j}(l)=t^{n}_{i}(k)\leq
t\\}}\in{\mathcal{F}}^{n}_{t^{n}_{ijkl}-}$. Hence by (3.20),

${\mathbb{E}}[\zeta^{n}_{i}(k)\zeta^{n}_{j}(l)1_{\\{t^{n}_{j}(l)=t^{n}_{i}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{t^{n}_{ijkl}-}]=0,\qquad k,l\geq 1,$

and (3.16) follows.

iii. Clearly $A^{n}_{i}$, defined in (2.6) is adapted and $A^{n}_{i}(t)$ is
integrable for all $t$. Hence the same is true for $\hat{M}^{n}_{i}$. Next,
let

$s^{n}(k)=\inf\\{t\geq 0:A^{n}_{0}(t)\geq k\\},\qquad k=1,2,\ldots.$

As in (i), these are stopping times and
$s^{n}(k)\in{\mathcal{F}}^{n}_{s^{n}(k)-}$. To show the martingale property,
we can write, using
$\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s-)}ds=\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s)}ds$,

$\displaystyle n^{1/2}\hat{M}^{A,n}_{i}(t)$
$\displaystyle=\int_{[0,t]}1_{\\{{\mathcal{R}}^{n}_{i}(s-)=\theta^{n}_{A^{n}_{0}(s)}\\}}dA^{n}_{0}(s)-\lambda^{n}_{0}\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s)}ds$
$\displaystyle=\int_{[0,t]}(1_{\\{{\mathcal{R}}^{n}_{i}(s-)=\theta^{n}_{A^{n}_{0}(s)}\\}}-p_{n,{\mathcal{R}}^{n}_{i}(s-)})dA^{n}_{0}(s)+\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s-)}(dA^{n}_{0}(s)-\lambda^{n}_{0}ds)$
$\displaystyle=:M^{n}_{i,1}(t)+M^{n}_{i,2}(t).$

For $M^{n}_{i,1}$, write

$A^{n}_{i}(t)=\sum_{k=1}^{A^{n}_{0}(t)}1_{\\{{\mathcal{R}}^{n}_{i}(s^{n}(k)-)=\theta^{n}_{k}\\}}.$

An argument similar to the one given before shows that $\theta^{n}_{k}$ is
independent of ${\mathcal{F}}^{n}_{s^{n}(k)-}$. Therefore, for $0\leq s<t$, we
have

$\displaystyle{\mathbb{E}}[A^{n}_{i}(t)|{\mathcal{F}}^{n}_{s}]-A^{n}_{i}(s)$
$\displaystyle=\sum_{k=1}^{\infty}{\mathbb{E}}[1_{\\{{\mathcal{R}}^{n}_{i}(s^{n}(k)-)=\theta^{n}_{k}\\}}1_{\\{s<s^{n}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{s}]$
$\displaystyle=\sum_{k=1}^{\infty}{\mathbb{E}}[{\mathbb{E}}[1_{\\{{\mathcal{R}}^{n}_{i}(s^{n}(k)-)=\theta^{n}_{k}\\}}1_{\\{s<s^{n}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{s^{n}(k)-}]|{\mathcal{F}}^{n}_{s}]$
$\displaystyle=\sum_{k=1}^{\infty}{\mathbb{E}}[p_{n,{\mathcal{R}}^{n}_{i}(s^{n}(k)-)}1_{\\{s<s^{n}(k)\leq
t\\}}|{\mathcal{F}}^{n}_{s}]$
$\displaystyle={\mathbb{E}}[C^{n}_{i}(t)|{\mathcal{F}}^{n}_{s}]-C^{n}_{i}(s),$

where

$C^{n}_{i}(t)=\sum_{k=1}^{A^{n}_{0}(t)}p_{n,{\mathcal{R}}^{n}_{i}(s^{n}(k)-)}=\int_{[0,t]}p_{n,{\mathcal{R}}^{n}_{i}(s-)}dA^{n}_{0}(s),$

showing that $A^{n}_{i}-C^{n}_{i}=M^{n}_{i,1}$ is a martingale.

In the expression for $M^{n}_{i,2}$, the integrand is
$\\{{\mathcal{F}}^{n}_{t}\\}$-adapted and has LCRL sample paths, while the
integrator is a martingale on the filtration. As a result, $M^{n}_{i,2}$ is a
local martingale [37, Theorem II.20]; using the estimate
$\|M^{n}_{i,2}\|^{*}_{t}\leq A^{n}_{0}(t)+c$ shows it is in fact a martingale.
As a result, so is $\hat{M}^{A,n}_{i}$. Finally, the expression for the
quadratic variation is straightforward. ∎

We will also need the following simple fact.

######  Lemma 3.4.

Let $M_{N}$, $N\in{\mathbb{Z}}_{+}$ be a martingale with $M_{0}=0$, for which
the increments $\mathnormal{\Delta}_{N}=M_{N}-M_{N-1}$ satisfy
${\mathbb{E}}(|\mathnormal{\Delta}_{N}|1_{\\{|\mathnormal{\Delta}_{N}|>a\\}})\leq\bar{r}(a)$,
$a\geq 0$, $N\in{\mathbb{N}}$, and $\bar{r}(a)\to 0$ as $a\to\infty$. Then
$N^{-1}{\mathbb{E}}\|M\|^{*}_{N}<A(N)\to 0$, where $\\{A(N)\\}$ depend only on
$\bar{r}$.

Proof. Let
$b_{N}={\mathbb{E}}[\mathnormal{\Delta}_{N}1_{\\{|\mathnormal{\Delta}_{N}|\leq
a\\}}]$ and note that
$b_{N}=-{\mathbb{E}}[\mathnormal{\Delta}_{N}1_{\\{|\mathnormal{\Delta}_{N}|>a\\}}]$
hence $|b_{N}|\leq\bar{r}(a)$. Write

$M_{N}=P_{N}+Q_{N},\qquad
P_{N}=\sum_{i=1}^{N}\\{\mathnormal{\Delta}_{i}1_{\\{|\mathnormal{\Delta}_{i}|>a\\}}+b_{i}\\},\qquad
Q_{N}=\sum_{i=1}^{N}\\{\mathnormal{\Delta}_{i}1_{\\{|\mathnormal{\Delta}_{i}|\leq
a\\}}-b_{i}\\}.$

The quadratic variation of the martingale $Q_{N}$ is bounded by
$(a+\bar{r}(a))^{2}n$, giving ${\mathbb{E}}[\|Q\|^{*}_{N}]\leq
c(a+\bar{r}(a))n^{1/2}$, where $c^{2}$ is the constant from the Burkholder-
Davis-Gundy (BDG) inequality with $p=2$. For $P_{N}$,

$|P_{N}|\leq\sum_{i=1}^{N}|\mathnormal{\Delta}_{i}|1_{\\{|\mathnormal{\Delta}_{i}|>a\\}}+N\bar{r}(a),$

thus ${\mathbb{E}}[\|P\|^{*}_{N}]\leq 2N\bar{r}(a)$. This gives
$n^{-1}{\mathbb{E}}[\|M\|^{*}_{N}]\leq c(a+\bar{r}(a))N^{-1/2}+2\bar{r}(a)$.
Taking $a=N^{1/4}$ completes the proof. ∎

###  3.3. Uniform estimates on individual queue length processes

The goal here is to calculate the dynamics of the individual rescaled
processes and develop estimates showing that they are $C$-tight uniformly in
$i\in[n]$.

In the following lemma, part (i) provides an equation, (3.21), for the
dynamics of individual queue lengths, and part (ii) shows that, at the cost of
introducing an error term, one can replace the term
$\hat{S}^{n}_{i}(T^{n}_{i})$ in that equation by a martingale. Both
representations (3.21) and (3.24) are used in this paper, where the former is
used for estimates on each $\hat{X}^{n}_{i}$ that are uniform in $i$, and the
latter is convenient for representations of $\bar{\xi}^{n}$, as the martingale
terms add up to a martingale. The last part of the lemma uses (3.21) and gives
uniform second moment and tightness estimates.

######  Lemma 3.5.

i. One has

(3.21) $\displaystyle\hat{X}^{n}_{i}(t)$
$\displaystyle=\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}(t)+\hat{A}^{n}_{i}(t)-\hat{S}^{n}_{i}(T^{n}_{i}(t))+\hat{b}^{n}_{1}t+\hat{L}^{n}_{i}(t).$

Moreover, the sample paths of $\hat{L}^{n}_{i}$ are in $C^{\uparrow}$ and

(3.22) $\int_{0}^{\infty}\hat{X}^{n}_{i}(t)d\hat{L}^{n}_{i}(t)=0.$

ii. One has

(3.23) $\hat{X}^{n}_{i}(t)=\hat{X}^{1,n}_{i}(t)+e^{1,n}_{i}(t),$

where

(3.24)
$\hat{X}^{1,n}_{i}(t)=\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}(t)+\hat{A}^{n}_{i}(t)-\hat{M}^{{\rm
dep},n}_{i}(t)+\hat{b}^{n}_{1}t+\hat{L}^{n}_{i}(t),$
$e^{1,n}_{i}(t)=\tilde{Z}^{n}_{i}(0)-\tilde{R}^{n}_{i}(T^{n}_{i}(t))=\tilde{\mu}^{n}(Z^{n}_{i}(0)-R^{n}_{i}(T^{n}_{i}(t))).$

iii. For
$H^{n}_{i}=\hat{S}^{n}_{i},\hat{E}^{n}_{i},\hat{A}^{n}_{i},\hat{L}^{n}_{i}$
and $\hat{X}^{n}_{i}$, one has

(3.25)
$\sup_{n}\sup_{i\in[n]}{\mathbb{E}}[(\|H^{n}_{i}\|^{*}_{t})^{2}]<\infty,\qquad
t\geq 0,$

and for every $t$, $\varepsilon>0$ and $\eta>0$ there is $\delta>0$ such that

(3.26)
$\limsup_{n}\max_{i\in[n]}{\mathbb{P}}(w_{t}(H^{n}_{i},\delta)>\varepsilon)<\eta.$

Proof. i. By (2.1) and (2.12),

$\displaystyle n^{-1/2}X^{n}_{i}(t)$
$\displaystyle=n^{-1/2}X^{n}_{i}(0-)+n^{-1/2}(E^{n}_{i}(t)-\lambda^{n}t)+n^{-1/2}(\lambda^{n}-n\lambda)t+n^{1/2}\lambda
t+n^{-1/2}A^{n}_{i}(t)$
$\displaystyle\quad-n^{-1/2}(S^{n}_{i}(T^{n}_{i}(t))-\mu^{n}T^{n}_{i}(t))-n^{-1/2}\mu^{n}T^{n}_{i}(t),$

and

$-n^{-1/2}\mu^{n}T^{n}_{i}(t)=-n^{-1/2}(\mu^{n}-n\mu)t-n^{1/2}\mu
t+\hat{L}^{n}_{i}(t).$

Using (2.11), (3.8) and (3.9) gives (3.21). The properties of
$\hat{L}^{n}_{i}$ and (3.22) follow from (2.2).

ii. Using (3.14) in (3.21) gives (3.23).

iii. By the central limit theorem for renewal processes [10, §17], and the
fact that, by (2.8) and (2.10), $n^{-1}(\lambda^{n},\mu^{n})\to(\lambda,\mu)$,
for each $i$, $(\hat{E}^{n}_{i},\hat{S}^{0,n}_{i})$ converge in law to
$(E,S)$, a pair of mutually independent BM starting at zero, with zero drift,
and diffusion coefficients $\lambda^{1/2}$ and $\mu^{1/2}\sigma_{\rm ser}$,
respectively (where we recall $\lambda=\mu$).

We prove that (3.25) and (3.26) hold for $\hat{S}^{n}_{i}$ by relating these
processes to $\hat{S}^{0,n}_{i}$, whose laws do not depend on $i$. To this
end, note that

$S^{n}_{i}(t)=S^{0,n}_{i}((t-Z^{n}_{i}(0))^{+})-1_{\\{t<Z^{n}_{i}(0)\\}}=S^{0,n}_{i}(t-Z^{\\#,n}_{i})-1_{\\{t<Z^{n}_{i}(0)\\}},$

where $Z^{\\#,n}_{i}=Z^{\\#,n}_{i}(t):=Z^{n}_{i}(0)\wedge t$. Hence by (3.8),

$\hat{S}^{n}_{i}(t)=\hat{S}^{0,n}_{i}(t-Z^{\\#,n}_{i})-\tilde{\mu}^{n}Z^{\\#,n}_{i}-n^{-1/2}1_{\\{t<Z^{n}_{i}(0)\\}}.$

As a result,

(3.27)
$\|\hat{S}^{n}_{i}\|^{*}_{t}\leq\|\hat{S}^{0,n}_{i}\|^{*}_{t}+\tilde{\mu}^{n}Z^{n}_{i}(0)+n^{-1/2}=\|\hat{S}^{0,n}_{i}\|^{*}_{t}+\tilde{Z}^{n}_{i}(0)+n^{-1/2},$

and

$w_{t}(\hat{S}^{n}_{i},\delta)\leq
w_{t}(\hat{S}^{0,n}_{i},\delta)+\tilde{\mu}^{n}Z^{n}_{i}(0)=w_{t}(\hat{S}^{0,n}_{i},\delta)+\tilde{Z}^{n}_{i}(0).$

The latter inequality and the fact that the law of $\hat{S}^{0,n}_{i}$ does
not depend on $i$ gives

$\max_{i\in[n]}{\mathbb{P}}(w_{t}(\hat{S}^{n}_{i},\delta)>\varepsilon)\leq{\mathbb{P}}\Big{(}w_{t}(\hat{S}^{0,n}_{1},\delta)>\frac{\varepsilon}{2}\Big{)}+\max_{i\in[n]}{\mathbb{P}}\Big{(}\tilde{Z}^{n}_{i}(0)>\frac{\varepsilon}{2}\Big{)}.$

Using (2.13) and the $C$-tightness of $\hat{S}^{0,n}_{1}$, $n\in{\mathbb{N}}$,
shows that (3.26) is satisfied by $\hat{S}^{n}_{i}$.

Next, under our second moment assumptions on the inter-renewal times, it is
well known that the rescaled non-delayed renewal processes satisfy

$\sup_{n}{\mathbb{E}}[(\|\hat{S}^{0,n}_{1}\|^{*}_{t})^{2}]<\infty,$

for every $t$ [28, Appendix 1]. Since the law of $\hat{S}^{0,n}_{i}$ does not
depend on $i$, it follows from (3.27) and our assumption (2.14) that
$\hat{S}^{n}_{i}$ satisfy (3.25).

It follows now that $\hat{E}^{n}_{i}$ satisfies both estimates, for the law of
$E^{n}_{i}$ is merely a special case of the law of the non-delayed renewal
processes $S^{0,n}_{i}$.

As for the processes $A^{n}_{i}$, recall from Lemma 3.3.iii that
$\hat{M}^{A,n}_{i}$ is a martingale and that
$[\hat{M}^{A,n}_{i}](t)=n^{-1}A^{n}_{i}(t)$. One has

${\mathbb{E}}[n^{-1}A^{n}_{i}(t)]=n^{-1}\lambda^{n}_{0}\int_{0}^{t}{\mathbb{E}}[p_{n,{\mathcal{R}}^{n}_{i}(s)}]ds\leq
cn^{-1}n^{3/2}n^{-1}t,$

where we used (2.5) and (2.9). Since this bound does not depend on $i$, it
follows by the BDG inequality that

(3.28)
$\lim_{n}\max_{i\in[n]}{\mathbb{E}}[(\|\hat{M}^{A,n}_{i}\|^{*}_{t})^{2}]=0.$

Again by the bound on $p_{n,r}$, the Lipschitz constant of $\hat{C}^{A,n}_{i}$
is bounded by $cn^{3/2}n^{-1/2}n^{-1}=c$. It follows that both estimates
(3.25) and (3.26) are satisfied by $\hat{A}^{n}_{i}$.

To treat $\hat{L}^{n}_{i}$, recall Skorohod’s lemma [16, §8], stating that for
a trajectory $y\in D({\mathbb{R}}_{+},{\mathbb{R}})$, if $(x,z)\in
D({\mathbb{R}}_{+},{\mathbb{R}}_{+})^{2}$ are such that $x=y+z$, $z$ is
nondecreasing and, with the convention $z(0-)=0$, $\int_{[0,\infty)}xdz=0$,
then $z$ is given by

$z(t)=\sup_{s\in[0,t]}y^{-}(s),\qquad t\geq 0.$

In particular,

(3.29) $z(t)\leq\|y\|_{t},\qquad w_{t}(z,\delta)\leq w_{t}(y,\delta),\qquad
t>0,\,\delta>0.$

In view of (3.21) and (3.22), it follows that

(3.30) $\hat{L}^{n}_{i}(t)\leq\|\hat{U}^{n}_{i}\|_{t},\qquad
w_{t}(\hat{L}^{n}_{i},\delta)\leq w_{t}(\hat{U}^{n}_{i},\delta),$

where

(3.31)
$\hat{U}^{n}_{i}(t)=\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}(t)+\hat{A}^{n}_{i}(t)-\hat{S}^{n}_{i}(T^{n}_{i}(t))+\hat{b}^{n}_{1}t,$

and

(3.32) $\hat{X}^{n}_{i}=\hat{U}^{n}_{i}+\hat{L}^{n}_{i}.$

Because $T^{n}_{i}$ are $1$-Lipschitz, we have that
$\hat{S}^{n}_{i}(T^{n}_{i}(\cdot))$ satisfy (3.25) and (3.26). Hence by the
second moment assumption on initial condition (2.15) and the boundedness of
$\hat{b}^{n}_{1}$, the same holds for $\hat{U}^{n}_{i}$. Finally, by (3.30)
and (3.32), this is also true for $\hat{L}^{n}_{i}$ and $\hat{X}^{n}_{i}$. ∎

###  3.4. Empirical process tightness

######  Lemma 3.6.

$\bar{\xi}^{n}$ are $C$-tight in
$D({\mathbb{R}}_{+},{\mathcal{M}}_{1}({\mathbb{R}}_{+}))$.

Proof. Denote by $C^{\varepsilon}$ the $\varepsilon$-neighborhood, in
${\mathbb{R}}_{+}$, of a set $C\in{\mathbb{R}}_{+}$, and let

$d_{\rm L}(p,q)=\inf\\{\varepsilon>0:p(C^{\varepsilon})+\varepsilon\geq
q(C)\text{ and }q(C^{\varepsilon})+\varepsilon\geq p(C)\text{ for all
}C\in{\mathcal{B}}({\mathbb{R}}_{+})\\},$

$p,q\in{\mathcal{M}}_{1}({\mathbb{R}}_{+})$, denote the Levy-Prohorov metric,
which induces the topology of weak convergence on
${\mathcal{M}}_{1}({\mathbb{R}}_{+})$. Since $\bar{\xi}^{n}_{0}$ converge in
probability, proving $C$-tightness of $\bar{\xi}^{n}$ amounts to showing that
for every $\varepsilon>0$ and $\eta>0$ there exists $\delta>0$ such that

(3.33)
$\limsup_{n}{\mathbb{P}}(w_{T}(\bar{\xi}^{n},\delta)>\varepsilon)<\eta.$

Here and below we use the same notation $w_{T}$ for the metric $d_{\rm L}$ and
for the usual metric on ${\mathbb{R}}$.

We show that (3.33) is a consequence of Lemma 3.5.iii. Fix $\varepsilon>0$ and
$\eta>0$. Given $n$ and $\delta>0$, consider the event
$\mathnormal{\Omega}^{n}_{\delta}=\\{w_{T}(\bar{\xi}^{n},\delta)>\varepsilon\\}$.
On this event there are $0\leq s\leq t\leq T$, $t-s\leq\delta$, and
$C\in{\mathcal{B}}({\mathbb{R}}_{+})$, such that

$\xi^{n}_{t}(C^{\varepsilon})+n\varepsilon<\xi^{n}_{s}(C)\quad\text{ or
}\quad\xi^{n}_{s}(C^{\varepsilon})+n\varepsilon<\xi^{n}_{t}(C).$

In both cases, the number of trajectories $\hat{X}^{n}_{i}$ whose displacement
between times $s$ and $t$ exceeds $\varepsilon$ is greater than
$n\varepsilon$. Therefore

$\\#\\{i\in[n]:w_{T}(\hat{X}^{n}_{i},\delta)\geq\varepsilon\\}>n\varepsilon.$

Hence by Chebychev’s inequality,

${\mathbb{P}}(\mathnormal{\Omega}^{n}_{\delta})\leq(n\varepsilon)^{-1}n\max_{i}{\mathbb{P}}(w_{T}(\hat{X}^{n}_{i},\delta)>\varepsilon).$

In view of (3.26), given any $\eta_{1}>0$ there is $\delta>0$ such that for
all large $n$, the maximum over $i$ in the above display is $<\eta_{1}$.
Hence, for such $\delta$, the above is $\leq\varepsilon^{-1}\eta_{1}$.
Choosing $\eta_{1}$ such that $\varepsilon^{-1}\eta_{1}<\eta$ and the
corresponding $\delta$ gives (3.33). ∎

##  4\. Proof

###  4.1. An equation for the empirical process

This section provides an equation for the empirical measure that is a
precursor to the parabolic PDE. Fix a test function $\tilde{\phi}$ as in
Definition 2.3, that is, $\tilde{\phi}\in C^{\infty}_{0}({\mathbb{R}}_{+})$,
$\tilde{\phi}(0)=0$. Let

(4.34) $\phi(x)=\int_{0}^{x}\tilde{\phi}(y)dy,\qquad x\in{\mathbb{R}}_{+}.$

Then $\phi\in C^{\infty}_{b}({\mathbb{R}}_{+})$ and
$\phi(0)=\phi^{\prime}(0)=0$. Moreover, $\phi,\phi^{\prime}$ and
$\phi^{\prime\prime}$ are uniformly continuous on ${\mathbb{R}}_{+}$. Apply
$\phi$ to the dynamics (3.24), noting that, unlike $X^{n}_{i}$, $X^{1,n}_{i}$
may assume negative values. On the r.h.s. of (3.24), the terms
$\hat{A}^{n}_{i}$ and $\hat{M}^{{\rm dep},n}_{i}$ are piecewise constant. Thus

$\displaystyle\phi(\hat{X}^{1,n}_{i}(t))$
$\displaystyle=\phi(\hat{X}^{n}_{i}(0-))+\int_{0}^{t}\phi^{\prime}(\hat{X}^{1,n}_{i}(s))(d\hat{E}^{n,c}_{i}(s)+\hat{b}^{n}_{1}ds)+e^{2,n}_{i}(t)$
(4.35) $\displaystyle\quad+\sum_{s\leq
t}(\phi(\hat{X}^{1,n}_{i}(s))-\phi(\hat{X}^{1,n}_{i}(s-))),$

where $\hat{E}^{n,c}_{i}$ is the continuous part of $\hat{E}^{n}_{i}$,

$e^{2,n}_{i}(t)=\int_{0}^{t}(\phi^{\prime}(\hat{X}^{1,n}_{i}(s))-\phi^{\prime}(\hat{X}^{n}_{i}(s)))d\hat{L}^{n}_{i}(s),$

and we have used (3.22) and the continuity of the sample paths of
$\hat{L}^{n}_{i}$ to write

$\int\phi^{\prime}(\hat{X}^{n}_{i}(s))d\hat{L}^{n}_{i}(s)=\phi^{\prime}(0)\hat{L}^{n}_{i}(t)=0.$

In the last term of (4.1), by Taylor’s expansion, jumps according to
$\hat{E}^{n}_{i}$ can be expressed as

$\phi^{\prime}(\hat{X}^{1,n}_{i}(s-))\mathnormal{\Delta}\hat{E}^{n}_{i}(s)+\frac{1}{2}\phi^{\prime\prime}(\hat{X}^{2,n}_{i}(s))\mathnormal{\Delta}\hat{E}^{n}_{i}(s)^{2},$

where $\hat{X}^{2,n}_{i}(s)$ is an intermediate value between
$\hat{X}^{1,n}_{i}(s-)$ and $\hat{X}^{1,n}_{i}(s)$ (we leave unspecified the
value of the processes $\hat{X}^{2,n}_{i}$ away from times of jumps of
$\hat{X}^{1,n}_{i}$). Similarly, jumps according to $\hat{M}^{{\rm
dep},n}_{i}$ are expressed as

$-\phi^{\prime}(\hat{X}^{1,n}_{i}(s-))\mathnormal{\Delta}\hat{M}^{{\rm
dep},n}_{i}(s)+\frac{1}{2}\phi^{\prime\prime}(\hat{X}^{2,n}_{i}(s))\mathnormal{\Delta}\hat{M}^{{\rm
dep},n}_{i}(s)^{2},$

where again $\hat{X}^{2,n}_{i}(s)$ are intermediate points, and because, a.s.,
jumps of $\hat{E}^{n}_{i}$ and $\hat{M}^{{\rm dep},n}_{i}$ do not occur
simultaneously (for the same $i$), we may express the intermediate values by
the same process. Finally, jumps according to $\hat{A}^{n}_{i}$ – here we only
need first order approximation – are expressed as

$\phi^{\prime}(\hat{X}^{2,n}_{i}(s))\mathnormal{\Delta}\hat{A}^{n}_{i}(s).$

Note that jumps of $\hat{E}^{n}_{i}$ and $\hat{A}^{n}_{0}$ are of size
$n^{-1/2}$. Hence in all three cases,

(4.36)
$|\hat{X}^{2,n}_{i}(t)-\hat{X}^{1,n}_{i}(t-)|\leq|\mathnormal{\Delta}\hat{X}^{1,n}_{i}(t)|\leq\max\\{|\mathnormal{\Delta}\hat{M}^{{\rm
dep},n}_{i}(t)|,n^{-1/2}\\},\qquad t\in{\mathcal{J}}^{n}_{i},$

where ${\mathcal{J}}^{n}_{i}$ is the set of jump times of $\hat{X}^{1,n}_{i}$.

By (4.1), using
$\hat{E}^{n}_{i}=\hat{E}^{n,c}_{i}+\mathnormal{\Delta}\hat{E}^{n}_{i}$ and
$\int_{0}^{t}\phi^{\prime}(\hat{X}^{1,n}_{i}(s))d\hat{E}^{n,c}_{i}(s)=\int_{0}^{t}\phi^{\prime}(\hat{X}^{1,n}_{i}(s-))d\hat{E}^{n,c}_{i}(s)$,
we have

(4.37)
$\displaystyle\langle\phi,\bar{\xi}^{n}_{t}\rangle=\langle\phi,\bar{\xi}^{n}_{0-}\rangle+\int_{0}^{t}\langle
b_{1}\phi^{\prime}+\frac{\sigma^{2}}{2}\phi^{\prime\prime},\bar{\xi}^{n}_{s}\rangle
ds+\mathnormal{\Gamma}^{n}(t)+\sum_{j=1}^{6}f^{j,n}(t),$

where

(4.38)
$\mathnormal{\Gamma}^{n}(t)=\frac{b_{0}}{n}\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))\Big{(}\frac{n-{\mathcal{R}}^{n}_{i}(s)}{n}\Big{)}^{\ell-1}ds$

is the interaction term,

$\displaystyle f^{1,n}(t)$ $\displaystyle=f^{1,n}_{1}(t)+f^{1,n}_{2}(t)$
$\displaystyle:=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{1,n}_{i}(s-))d\hat{E}^{n}_{i}(s)-\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{1,n}_{i}(s-))d\hat{M}^{{\rm
dep},n}_{i}(s),$ $\displaystyle f^{2,n}(t)$
$\displaystyle=\langle\phi,\bar{\xi}^{n}_{t}\rangle-\langle\phi,\bar{\xi}^{1,n}_{t}\rangle,$
$\displaystyle f^{3,n}(t)$ $\displaystyle=\frac{1}{n}\sum_{i}e^{2,n}_{i}(t),$
$\displaystyle f^{4,n}(t)$
$\displaystyle=\hat{b}^{n}_{1}\int_{0}^{t}\langle\phi^{\prime},\bar{\xi}^{1,n}_{s}\rangle
ds-b_{1}\int_{0}^{t}\langle\phi^{\prime},\bar{\xi}^{n}_{s}\rangle ds,$
$\displaystyle f^{5,n}(t)$
$\displaystyle=\frac{1}{2n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{2,n}_{i}(s))\\{d[\hat{E}^{n}_{i}](s)+d[\hat{M}^{{\rm
dep},n}_{i}](s)\\}-\frac{\sigma^{2}}{2}\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{n}_{s}\rangle
ds,$ $\displaystyle f^{6,n}(t)$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{2,n}_{i}(s))d\hat{A}^{n}_{i}(s)-\mathnormal{\Gamma}^{n}(t)$

are “error” terms, and

$\bar{\xi}^{1,n}_{t}=n^{-1}\sum_{i\in[n]}\delta_{\hat{X}^{n,1}_{i}(t)}.$

Note that the terms $f^{5,n}$ and $f^{6,n}$, which involve
$\hat{X}^{2,n}_{i}$, are well defined, for their evaluation requires the
values of $\hat{X}^{2,n}_{i}(t)$ only at ${\mathcal{J}}^{n}_{i}$.

######  Lemma 4.1.

With $\phi$ as in (4.34), for $1\leq j\leq 6$, $f^{j,n}\to 0$ in probability
in $D({\mathbb{R}}_{+},{\mathbb{R}})$.

We note that one of the reasons the proof is somewhat involved is that we work
under minimal moment assumptions. For example, estimates on the second order
term $f^{5,n}$ use martingales that are driven by squares of the primitive
data $\zeta^{n}_{i}(k)$. Since the latter are assumed only to possess a second
moment, estimates based on the corresponding quadratic variation are not
applicable, for it does not, in general, have finite expectation.

Proof. Step 1. Estimating $f^{1,n}$. The cross variation between independent
Poisson processes is zero, hence $[\hat{E}^{n}_{i},\hat{E}^{n}_{j}](t)=0$ for
$i\neq j$. Hence by Lemma 3.3,

$[f^{1,n}_{1}](t)\leq cn^{-2}\sum_{i}[\hat{E}^{n}_{i}](t)\leq
cn^{-3}\sum_{i}E^{n}_{i}(t).$

Since each $E^{n}_{i}$ is a Poisson process of intensity $\leq cn$,
${\mathbb{E}}\\{[f^{1,n}_{1}](t)\\}\leq ctn^{-1}$, and it follows by the BDG
inequality that $f^{1,n}_{1}\to 0$ in probability. As for $f^{1,n}_{2}$, we
have by Lemma 3.3.ii that $[\hat{M}^{{\rm dep},n}_{i},\hat{M}^{{\rm
dep},n}_{j}](t)$, $i\neq j$, has zero mean. Hence by Lemma 3.3.i,

${\mathbb{E}}\\{[f^{1,n}_{2}](t)\\}\leq
cn^{-2}\sum_{i}{\mathbb{E}}\\{[\hat{M}^{{\rm dep},n}_{i}](t)\\}\leq
cn^{-3}\sum_{i}{\mathbb{E}}\\{D^{n}_{i}(t)\\}\leq cn^{-1}t.$

Thus $f^{1,n}_{2}\to 0$ in probability, and we conclude that $f^{1,n}\to 0$ in
probability.

Step 2. We show that for all $t>0$ and $\varepsilon>0$,

(4.39)
$\lim_{n\to\infty}\max_{i\in[n]}{\mathbb{P}}(\|e^{1,n}_{i}\|^{*}_{t}>\varepsilon)=0.$

In accordance with (2.13) we denote
$\tilde{Z}^{n}_{i}(k)=\tilde{\mu}^{n}Z^{n}_{i}(k)$, $k\geq 1$. By the
definition of $e^{1,n}_{i}(t)$ and $T^{n}_{i}(t)\leq t$,

$\displaystyle|e^{1,n}_{i}(t)|\leq\tilde{Z}^{n}_{i}(0)\vee\|\tilde{R}^{n}_{i}\|^{*}_{t}$
$\displaystyle\leq\max\\{\tilde{Z}^{n}_{i}(k):0\leq k\leq S^{n}_{i}(t)\\}$
(4.40) $\displaystyle\leq\max\\{\tilde{Z}^{n}_{i}(k):0\leq k\leq
S^{0,n}_{i}(t)\\}.$

Thus, by (2.13) and the fact that the law of $((\tilde{Z}^{n}_{i}(k),k\geq
1),S^{0,n}_{i}(t))$ does not depend on $i$, it suffices to prove that

$\tilde{Y}^{n}:=\max\\{\tilde{Z}^{n}_{1}(k):1\leq k\leq S^{0,n}_{1}(t)\\}\to
0\text{ in probability.}$

This can be argued via the $C$-tightness of $\\{\hat{S}^{0,n}_{1}\\}$ as
follows. Given $n$ and $t_{1}<t_{2}$, we have

(4.41) $\text{if $S^{0,n}_{1}(t_{1})=S^{0,n}_{1}(t_{2})$ then
}\hat{S}^{0,n}_{1}(t_{1})-\hat{S}^{0,n}_{1}(t_{2})=\tilde{\mu}^{n}(t_{2}-t_{1}).$

This shows that

$\max\\{Z^{n}_{1}(k):1\leq k\leq S^{0,n}_{1}(t)-1\\}\leq
2(\tilde{\mu}^{n})^{-1}\|\hat{S}^{0,n}_{1}\|^{*}_{t}.$

To include also $k=S^{0,n}_{1}(t)$, let the residual time process
$R^{0,n}_{1}$ be defined analogously to $R^{n}_{1}$, for $S^{0,n}_{1}$ in
place of $S^{n}_{1}$. Note that if $R^{0,n}_{1}(t)>1$ holds then (4.41) holds
with $t_{1}=t$ and $t_{2}=t+1$ hence
$\hat{S}^{0,n}_{1}(t)-\hat{S}^{0,n}_{1}(t+1)=\tilde{\mu}^{n}$. Because
$\tilde{\mu}^{n}\to\infty$ and $\|\hat{S}^{0,n}_{1}\|^{*}_{t+1}$ are tight,
this shows that w.h.p., $R^{0,n}_{1}(t)\leq 1$. As a result, again by (4.41),
w.h.p.,

$Y^{n}:=\max\\{Z^{n}_{1}(k):1\leq k\leq S^{0,n}_{1}(t)\\}\leq
2(\tilde{\mu}^{n})^{-1}\|\hat{S}^{0,n}_{1}\|^{*}_{t+1}.$

This shows that $Y^{n}\to 0$ in probability. Next, given $\delta>0$, on the
event $Y^{n}<\delta$, we have by (4.41),

$\displaystyle w_{t+1}(\hat{S}^{0,n}_{1},\delta)$
$\displaystyle\geq\sup\\{|\hat{S}^{0,n}_{1}(t_{1})-\hat{S}^{0,n}_{1}(t_{2})|:0\leq
t_{1}<t_{2}\leq t+1,\,S^{0,n}_{1}(t_{1})=S^{0,n}_{1}(t_{2})\\}$
$\displaystyle\geq\tilde{\mu}^{n}\sup\\{t_{2}-t_{1}:0\leq t_{1}<t_{2}\leq
t+1,\,S^{0,n}_{1}(t_{1})=S^{0,n}_{1}(t_{2})\\}$
$\displaystyle\geq\tilde{\mu}^{n}Y^{n}.$

Hence, w.h.p., $\tilde{Y}^{n}=\tilde{\mu}^{n}Y^{n}\leq
w_{t+1}(\hat{S}^{0,n}_{1},\delta)$. Using the $C$-tightness of
$\hat{S}^{0,n}_{1}$, sending $n\to\infty$ and then $\delta\to 0$ shows that
$\tilde{Y}^{n}\to 0$ in probability. This proves (4.39).

Step 3. We can now control $f^{2,n}$ and $f^{4,n}$. Recall (3.23). Then,
denoting by $m_{\phi}(\cdot)$ the modulus of continuity of $\phi$, for any
$\delta>0$,

$\|f^{2,n}\|^{*}_{t}\leq\frac{1}{n}\sum_{i}\|\phi(\hat{X}^{n}_{i})-\phi(\hat{X}^{1,n}_{i})\|^{*}_{t}\leq
m_{\phi}(\delta)+\frac{2\|\phi\|_{\infty}}{n}\\#\\{i:\|e^{1,n}_{i}\|^{*}_{t}>\delta\\}.$

Given $\varepsilon>0$ let $\delta>0$ be such that
$m_{\phi}(\delta)<\varepsilon/2$. Then

${\mathbb{P}}(\|f^{2,n}\|^{*}_{t}>\varepsilon)\leq
4\|\phi\|_{\infty}\varepsilon^{-1}\max_{i}{\mathbb{P}}(\|e^{1,n}_{i}\|^{*}_{t}>\delta).$

Sending $n\to\infty$, the above converges to $0$ by (4.39), proving that
$f^{2,n}\to 0$ in probability.

Next, similarly,
$\langle\phi^{\prime},\bar{\xi}^{1,n}\rangle-\langle\phi^{\prime},\bar{\xi}^{n}\rangle\to
0$ in probability, and

$\|f^{4,n}\|^{*}_{t}\leq
t\,|\hat{b}^{n}_{1}|\,\|\langle\phi^{\prime},\bar{\xi}^{1,n}\rangle-\langle\phi^{\prime},\bar{\xi}^{n}\rangle\|^{*}_{t}+|\hat{b}^{n}_{1}-b_{1}|\int_{0}^{\cdot}\langle\phi^{\prime},\bar{\xi}^{n}_{s}\rangle
ds.$

The first term on the right converges to $0$ in probability by the argument
shown for $f^{2,n}$ and the boundedness of $\hat{b}^{n}_{1}$. In the second
term, the integral is bounded by $t\|\phi^{\prime}\|_{\infty}$, hence the
convergence of this term to zero follows from $\hat{b}^{n}_{1}\to b_{1}$. This
proves $f^{4,n}\to 0$ in probability.

Step 4. To estimate $f^{3,n}$, note that, for any $\delta>0$,

$\|f^{3,n}\|^{*}_{t}\leq\frac{m_{\phi^{\prime}}(\delta)}{n}\sum_{i}\hat{L}^{n}_{i}(t)+\frac{2\|\phi^{\prime}\|_{\infty}}{n}\sum_{i}1_{\\{\|e^{1,n}_{i}\|^{*}_{t}>\delta\\}}\hat{L}^{n}_{i}(t).$

Hence

${\mathbb{E}}[\|f^{3,n}\|^{*}_{t}]\leq\frac{m_{\phi^{\prime}}(\delta)}{n}\sum_{i}{\mathbb{E}}[\hat{L}^{n}_{i}(t)]+\frac{2\|\phi^{\prime}\|_{\infty}}{n}\sum_{i}{\mathbb{P}}(\|e^{1,n}_{i}\|^{*}_{t}>\delta)^{1/2}({\mathbb{E}}[\hat{L}^{n}_{i}(t)^{2}])^{1/2}.$

Since by Lemma 3.5.iii,
$\max_{i\in[n]}{\mathbb{E}}[\hat{L}^{n}_{i}(t)^{2}]<c$, it follows that

${\mathbb{E}}[\|f^{3,n}\|^{*}_{t}]\leq
cm_{\phi^{\prime}}(\delta)+c\max_{i\in[n]}{\mathbb{P}}(\|e^{1,n}_{i}\|^{*}_{t}>\delta).$

In view of (4.39), if we take $n\to\infty$ and then $\delta\to 0$, the
expression on the right converges to $0$, hence $f^{3,n}\to 0$ in probability.

Step 5. Recall

$f^{5,n}(t)=\frac{1}{2n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{2,n}_{i}(s))\\{d[\hat{E}^{n}_{i}](s)+d[\hat{M}^{{\rm
dep},n}_{i}](s)\\}-\frac{\sigma^{2}}{2}\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{n}_{s}\rangle
ds,$ $[\hat{E}^{n}_{i}](t)=n^{-1}E^{n}_{i}(t)\qquad[\hat{M}^{{\rm
dep},n}_{i}](t)=n^{-1}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)^{2}.$

Let

$\displaystyle f^{5,n}_{1}$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s-))d[\hat{E}^{n}_{i}](s)-\lambda\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{1,n}_{s}\rangle
ds,$ $\displaystyle f^{5,n}_{2}$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s-))d[\hat{M}^{{\rm
dep},n}_{i}](s)-\lambda\sigma_{\rm
ser}^{2}\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{1,n}_{s}\rangle ds,$
$\displaystyle f^{5,n}_{3}$
$\displaystyle=\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{n}_{s}\rangle
ds-\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{1,n}_{s}\rangle ds$
$\displaystyle f^{5,n}_{4}$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}(\phi^{\prime\prime}(\hat{X}^{2,n}_{i}(s))-\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s-)))\\{d[\hat{E}^{n}_{i}](s)+d[\hat{M}^{{\rm
dep},n}_{i}](s)\\}.$

Then $f^{5,n}=\sum_{j}f^{5,n}_{j}$. For $f^{5,n}_{1}$, write

$\displaystyle f^{5,n}_{1}$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s-))n^{-1}dE^{n}_{i}(s)-\lambda\int_{0}^{t}\langle\phi^{\prime\prime},\bar{\xi}^{1,n}_{s}\rangle
ds$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s-))n^{-1/2}d\hat{E}^{n}_{i}(s)+\frac{1}{n}\sum_{i}\int_{0}^{t}\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(s))n^{-1/2}\hat{\lambda}^{n}ds.$

In the first sum, each term is a martingale, with expected quadratic variation
$\leq cn^{-2}\lambda^{n}t\leq cn^{-1}t$ where we use Lemma 3.3.i, and $c$ does
not depend on $i$ or $n$. Hence the first normalized sum is a martingale whose
expected quadratic variation is $\leq cn^{-2}t$. In particular, it converges
to zero in probability. In the second sum, each term is bounded in absolute
value by $ctn^{-1/2}$, hence the second normalized sum converges to zero. This
shows $f^{5,n}_{1}\to 0$ in probability.

Using Lemma 3.3 for the expression for $[\hat{M}^{{\rm dep},n}_{i}]$ and
denoting

$P^{n}_{i}(t)=\phi^{\prime\prime}(\hat{X}^{1,n}_{i}(t)),\qquad
Q^{n}_{i}(t)=\sum_{k=1}^{D^{n}_{i}(t)}q^{n}_{i}(k),\qquad
q^{n}_{i}(k)=\zeta^{n}_{i}(k)^{2}-\sigma_{\rm ser}^{2},$

write $f^{5,n}_{2}$ as

(4.42) $\displaystyle f^{5,n}_{2}$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}P^{n}_{i}(s-)n^{-1}dQ^{n}_{i}(s)+\frac{\sigma_{\rm
ser}^{2}}{n}\sum_{i}\int_{[0,t]}P^{n}_{i}(s-)\\{n^{-1}dD^{n}_{i}(s)-\lambda
ds\\}.$

To estimate the first term above, note that, by (3.8),

$D^{n}_{i}(t)=S^{n}_{i}(T^{n}_{i}(t))\leq
S^{n}_{i}(t)=n^{1/2}\hat{S}^{n}_{i}(t)+\mu^{n}t\leq
n^{1/2}\hat{S}^{n}_{i}(t)+c_{1}nt,$

for a constant $c_{1}$. For any $\beta>2c_{1}t$,

$\max_{i\in[n]}{\mathbb{P}}(D^{n}_{i}(t)>\beta
n)\leq\max_{i\in[n]}{\mathbb{P}}(\hat{S}^{n}_{i}(t)>n^{1/2}\beta/2)\leq
cn^{-1}\beta^{-2},$

by Lemma 3.5.iii, where $c=c(t)$. Denote the jump times of $D^{n}_{i}$ by
$0\leq t^{n}_{i}(1)<t^{n}_{i}(2)<\cdots$, and
$m^{n}_{i}(N)=\sum_{k=1}^{N}P^{n}_{i}(t^{n}_{i}(k)-)q^{n}_{i}(k)$. Then the
first term in (4.42) is given by $n^{-2}\sum_{i}m^{n}_{i}(D^{n}_{i}(t))$,
which we write as

(4.43) $\frac{1}{n^{2}}\sum_{i}m^{n}_{i}(D^{n}_{i}(t)\wedge(\beta
n))+\frac{1}{n^{2}}\sum_{i}\\{m^{n}_{i}(D^{n}_{i}(t)-m^{n}_{i}(D^{n}_{i}(t)\wedge(\beta
n))\\}.$

Now, $m^{n}_{i}(\cdot)$ is a discrete parameter martingale with increments
bounded by $\|\phi^{\prime\prime}\|_{\infty}|q^{n}_{i}(k)|$. Moreover,
$q^{n}_{i}(k)$ are IID in $k$ and $i$, and possesses a first moment in view of
our second moment assumption on the service times. Hence Lemma 3.4 is
applicable, showing that

$N^{-1}{\mathbb{E}}\|m^{n}_{i}\|^{*}_{N}<A(N)\to 0\qquad\text{as }N\to\infty,$

where $A(\cdot)$ does not depend on $n$ ot $i$. Hence

${\mathbb{E}}\Big{|}\frac{1}{n^{2}}\sum_{i}m^{n}_{i}(D^{n}_{i}(t)\wedge(\beta
n))\Big{|}\leq\frac{1}{n}\sum_{i}n^{-1}{\mathbb{E}}[\|m^{n}_{i}\|^{*}_{\beta
n}]\leq n^{-1}\beta nA(\beta n)=\beta A(\beta n),$

where throughout this paragraph $\beta n$ should be read as $\lfloor\beta
n\rfloor$. Hence the limit of the above expression is $0$ for any $\beta$.
Next, using $D^{n}_{i}\leq S^{n}_{i}$,

$\displaystyle{\mathbb{E}}\Big{|}\frac{1}{n^{2}}\sum_{i}\\{m^{n}_{i}(D^{n}_{i}(t)-m^{n}_{i}(D^{n}_{i}(t)\wedge(\beta
n))\\}\Big{|}$
$\displaystyle\leq\|\phi^{\prime\prime}\|_{\infty}n^{-2}\sum_{i}1_{\\{S^{n}_{i}(t)>\beta
n\\}}\sum_{k=\beta n+1}^{S^{n}_{i}(t)}|q^{n}_{i}(k)|$
$\displaystyle=cn^{-2}\sum_{i}\sum_{k=1}^{(S^{n}_{i}(t)-\beta
n)^{+}}|q^{n}_{i}(k-\beta n)|$
$\displaystyle=cn^{-2}\sum_{i}{\mathbb{E}}[(S^{n}_{i}(t)-\beta
n)^{+}]{\mathbb{E}}[|q^{n}_{i}(1)|]$ $\displaystyle\leq
cn^{-1}\max_{i\in[n]}{\mathbb{E}}[(S^{n}_{i}(t)-\beta n)^{+}],$

where Wald’s identity is used on the third line. Since it follows from Lemma
3.5.iii that $n^{-1}S^{n}_{i}(t)$ are uniformly integrable in $i$ and $n$,
$\lim_{\beta\to\infty}\limsup_{n}n^{-1}\max_{i\in[n]}{\mathbb{E}}[(S^{n}_{i}(t)-\beta
n)^{+}]=0$. As a result, the expression in (4.43) converges to zero in
probability.

As for the second term in (4.42). In view of Lemma 3.5.iii and (2.8), one has

$\max_{i\in[n]}\|n^{-1}S^{n}_{i}-\lambda\iota\|^{*}_{t}\to 0\qquad\text{in
probability}.$

Also from Lemma 3.5.iii and the relation (2.12) between $\hat{L}^{n}_{i}$ and
$T^{n}_{i}$, $\max_{i\in[n]}\|T^{n}_{i}-\iota\|^{*}_{t}\to 0$ in probability.
It follows that $D^{n}_{i}=S^{n}_{i}(T^{n}_{i})$ also satisfies
$\kappa_{n}:=\max_{i\in[n]}\|n^{-1}D^{n}_{i}-\lambda\iota\|^{*}_{t}\to 0$ in
probability. Denote $\varepsilon=t/j$, $s_{k}=k\varepsilon$ and
$\tilde{P}^{n}_{i}(s)=P^{n}_{i}(s-)$. Then

$\displaystyle\int_{(0,t]}\tilde{P}^{n}_{i}(s)\\{n^{-1}dD^{n}_{i}(s)-\lambda
ds\\}$
$\displaystyle=\sum_{k=0}^{j-1}\Big{\\{}\int_{(s_{k},s_{k+1}]}(\tilde{P}^{n}_{i}(s)-\tilde{P}^{n}_{i}(s_{k}))\\{n^{-1}dD^{n}_{i}(s)-\lambda
ds\\}$
$\displaystyle\qquad\qquad+\tilde{P}^{n}_{i}(s_{k})\\{n^{-1}D^{n}_{i}(s_{k+1})-n^{-1}D^{n}_{i}(s_{k})-\lambda\varepsilon\\}\Big{\\}}.$

This gives

(4.44)
$\displaystyle\Big{|}\int_{(0,t]}\tilde{P}^{n}_{i}(s)\\{n^{-1}dD^{n}_{i}(s)-\lambda
ds\\}\Big{|}$ $\displaystyle\leq
w_{T}(\tilde{P}^{n}_{i},\varepsilon)(n^{-1}D^{n}_{i}(t)+\lambda
t)+\|\tilde{P}^{n}_{i}\|^{*}_{t}2j\kappa_{n}$ $\displaystyle\leq
w_{T}(\tilde{P}^{n}_{i},\varepsilon)(2\lambda
t+\kappa_{n})+\|\phi^{\prime\prime}\|_{\infty}2j\kappa_{n}.$

We recall from Lemma 3.5.iii that $\hat{X}^{n}_{i}$ are uniformly $C$-tight.
Along with the uniform estimate (4.39) on $e^{1,n}_{i}$, this shows that so
are $\tilde{P}^{n}_{i}$. In particular, these processes satisfy (3.26). Hence,
on sending $n\to\infty$ and then $\varepsilon\to 0$, it follows that the right
side of (4.44) converges to zero in probability, uniformly in $i\in[n]$. Hence
the second term in $f^{5,n}_{2}$ converges to zero in probability, and we
conclude that $f^{5,n}_{2}\to 0$ in probability.

Next, recalling the result regarding $f^{2,n}$ and replacing $\phi$ by
$\phi^{\prime\prime}$ shows that $f^{5,n}_{3}\to 0$ in probability.

To bound $f^{5,n}_{4}$, fix $\varepsilon>0$ and let $\delta>0$ be such that
$|x-y|<\delta$ implies
$|\phi^{\prime\prime}(x)-\phi^{\prime\prime}(y)|<\varepsilon$. Let
${\mathcal{J}}^{n,\delta}_{i}=\\{s\in{\mathcal{J}}^{n}_{i}:|\hat{X}^{2,n}_{i}(s)-\hat{X}^{1,n}_{i}(s-)|>\delta\\}$.
For $n$ sufficiently large, the size of jumps of $\hat{E}^{n}$ (which is
$n^{-1/2}$) is smaller than $\delta$, hence, by (4.36), the corresponding jump
times are not members of ${\mathcal{J}}^{n,\delta}_{i}$. Thus the $i$-th term
in $f^{5,n}_{4}$ is bounded, in absolute value, by

$\displaystyle\varepsilon([\hat{E}^{n}_{i}](t)+[\hat{M}^{{\rm
dep},n}_{i}](t))+2\|\phi^{\prime\prime}\|_{\infty}\sum_{s\in{\mathcal{J}}^{n,\delta}_{i},s\leq
t}\mathnormal{\Delta}[\hat{M}^{{\rm dep},n}_{i}](s)$
$\displaystyle\leq\varepsilon
n^{-1}\Big{(}E^{n}_{i}(t)+\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)^{2}\Big{)}+cn^{-1}\sum_{k=1}^{D^{n}_{i}(t)}\zeta^{n}_{i}(k)^{2}1_{\\{n^{-1/2}|\zeta^{n}_{i}(k)|>\delta\\}}$

where we used the fact that for $s\in{\mathcal{J}}^{n,\delta}_{i}$,
$|\mathnormal{\Delta}\hat{X}^{1,n}_{i}(s)|=|\mathnormal{\Delta}\hat{M}^{{\rm
dep},n}_{i}(s)|$ by (3.24). Using (3.15), the expected value of the above
expression is

$\leq\varepsilon n^{-1}(\lambda^{n}_{i}t+\sigma_{\rm
ser}^{2}{\mathbb{E}}[D^{n}_{i}(t)])+cn^{-1}{\mathbb{E}}[\zeta^{n}_{1}(1)^{2}1_{\\{n^{-1/2}|\zeta^{n}_{1}(1)|>\delta\\}}]{\mathbb{E}}[D^{n}_{i}(t)],$

where for the last term we used Wald’s identity in exactly the same way as in
the proof of Lemma 3.3. In view of Lemma 3.5.iii,
${\mathbb{E}}[S^{n}_{i}(t)]<cn$, where $c$ does not depend on $i$ or $n$.
Since $D^{n}_{i}(t)\leq S^{n}_{i}(t)$, it follows that

${\mathbb{E}}[|f^{5,n}_{4}|]\leq
c\varepsilon+c{\mathbb{E}}[\zeta^{n}_{1}(1)^{2}1_{\\{n^{-1/2}|\zeta^{n}_{1}(1)|>\delta\\}}].$

Taking $n\to\infty$ then $\varepsilon\to 0$ shows that $f^{5,n}_{4}\to 0$ in
probability.

Step 6. Finally,

$f^{6,n}(t)=\frac{1}{n}\sum_{i}\int_{[0,t]}[\phi^{\prime}(\hat{X}^{2,n}_{i}(s))-\phi^{\prime}(\hat{X}^{n}_{i}(s))]d\hat{A}^{n}_{i}(s)+\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{n}_{i}(s))d\hat{A}^{n}_{i}(s)-\mathnormal{\Gamma}^{n}(t).$

By Lemma 3.3.iii,

$\displaystyle\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{n}_{i}(s-))d\hat{A}^{n}_{i}(s)$
$\displaystyle=\lambda^{n}_{0}n^{-3/2}\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))p_{n,{\mathcal{R}}^{n}_{i}(s)}ds+\bar{M}^{A,n}(t),$
$\displaystyle\bar{M}^{A,n}(t)$
$\displaystyle=\frac{1}{n}\sum_{i}\int_{[0,t]}\phi^{\prime}(\hat{X}^{n}_{i}(s-))d\hat{M}^{A,n}_{i}(s).$

By Lemma 3.3.iii, $[\bar{M}^{A,n}](t)\leq
cn^{-3}\sum_{i}A^{n}_{i}(t)=cn^{-3}A^{n}_{0}(t)$. Since $A^{n}_{0}$ is a
Poisson process of intensity $\lambda^{n}_{0}\leq cn^{3/2}$ (by (2.9)), it
follows that $\bar{M}^{A,n}\to 0$ in probability.

Now, writing

$p_{n,r}=\frac{\ell}{n}\,\frac{(n-r)(n-r-1)\cdots(n-r-\ell+2)}{n(n-1)\cdots(n-\ell+1)}$

shows that

(4.45)
$p_{n,r}=\frac{\ell}{n}\Big{[}\Big{(}1-\frac{r}{n}\Big{)}^{\ell-1}+\alpha_{n,r}\Big{]},\qquad\alpha_{n}:=\max_{r\in[n]}|\alpha_{n,r}|\to
0.$

Hence

$\displaystyle\lambda^{n}_{0}n^{-3/2}\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))p_{n,{\mathcal{R}}^{n}_{i}(s)}ds-\mathnormal{\Gamma}^{n}(t)$
$\displaystyle\qquad=\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))\Big{[}\lambda^{n}_{0}n^{-3/2}p_{n,{\mathcal{R}}^{n}_{i}(s)}-\frac{b_{0}}{n}\Big{(}\frac{n-{\mathcal{R}}^{n}_{i}(s)}{n}\Big{)}^{\ell-1}\Big{]}ds$
$\displaystyle\qquad=\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))\Big{(}\frac{n-{\mathcal{R}}^{n}_{i}(s)}{n}\Big{)}^{\ell-1}\Big{[}\lambda^{n}_{0}n^{-3/2}\frac{\ell}{n}-\frac{b_{0}}{n}\Big{]}ds+\sum_{i}\int_{0}^{t}\phi^{\prime}(\hat{X}^{n}_{i}(s))\lambda^{n}_{0}n^{-3/2}\frac{\ell}{n}\alpha_{n,{\mathcal{R}}^{n}_{i}(s)}ds.$

We have $\lambda^{n}_{0}n^{-3/2}\ell\to b\ell=b_{0}$, which shows that the
first sum converges to $0$. From the fact $\alpha_{n}\to 0$ we also have that
the last sum converges to $0$.

To bound the first term in $f^{6,n}$, note that $\hat{A}^{n}_{i}$ are
nondecreasing. Also, since the jumps of $\hat{A}^{n}_{i}$ are of size
$n^{-1/2}$, one has
$|\phi^{\prime}(\hat{X}^{2,n}_{i}(s)-\phi^{\prime}(\hat{X}^{1,n}_{i}(s))|\leq
m_{\phi^{\prime}}(n^{-1/2})$ at any jump time $s$ of $\hat{A}^{n}_{i}$.
Moreover,
$\max_{i\in[n]}|\phi^{\prime}(\hat{X}^{1,n}_{i}(s)-\phi^{\prime}(\hat{X}^{n}_{i}(s))|\leq\max_{i\in[n]}m_{\phi^{\prime}}(\|e^{1,n}_{i}\|^{*}_{t})\to
0$ in probability, by (4.39). Hence

$\max_{i\in[n]}\Big{|}\int_{[0,t]}[\phi^{\prime}(\hat{X}^{2,n}_{i}(s))-\phi^{\prime}(\hat{X}^{n}_{i}(s))]d\hat{A}^{n}_{i}(s)\Big{|}\leq\max_{i\in[n]}\\{(m_{\phi^{\prime}}(n^{-1/2})+m_{\phi^{\prime}}(\|e^{1,n}_{i}\|^{*}_{t}))\hat{A}^{n}_{i}(t)\\}\to
0$

in probability, where Lemma 3.5.iii is used for a uniform estimate on
$\hat{A}^{n}_{i}(t)$. It follows that $f^{6,n}\to 0$ in probability. ∎

###  4.2. Interaction term under limit

In this subsection we prove the following.

######  Lemma 4.2.

Given a subsequence along which $\bar{\xi}^{n}\Rightarrow\xi$, one has
$\mathnormal{\Gamma}^{n}\Rightarrow\mathnormal{\Gamma}$ in
$C({\mathbb{R}}_{+},{\mathbb{R}})$ along the subsequence, where

$\mathnormal{\Gamma}(t)=\frac{b_{0}}{\ell}\int_{0}^{t}\int_{{\mathbb{R}}_{+}}\phi^{\prime}(x)\mathfrak{S}(\xi_{s}[x,\infty),\xi_{s}(x,\infty))\xi_{s}(dx)ds.$

Recall that if $\xi_{t}$ is atomless for all $t>0$ then the integrand in the
above expression simplifies to
$\ell\phi^{\prime}(x)\xi_{s}[x,\infty)^{\ell-1}$, in which case
$\mathnormal{\Gamma}$ is directly related to the form of the interaction term
in equations (1.1), (1.3) and (2.16). However, the atomless property has not
been established at this stage of the proof.

Proof. Invoking the Skorohod representation theorem, assume without loss of
generality that $\bar{\xi}^{n}\to\xi$ a.s. Because
$0\leq{\mathcal{R}}^{n}_{i}\leq n$, the integrands in (4.38) are bounded by
$\|\phi^{\prime}\|_{\infty}$. Hence by bounded convergence, it suffices to
prove that, a.s., for every $t$,

(4.46)
$\gamma^{n}(t):=\frac{1}{n}\sum_{i}\phi^{\prime}(\hat{X}^{n}_{i}(t))\bar{\mathcal{R}}^{c,n}_{i}(t)^{\ell-1}\to\frac{1}{\ell}\int_{{\mathbb{R}}_{+}}\phi^{\prime}(x)\mathfrak{S}(\xi_{t}[x,\infty),\xi_{t}(x,\infty))\xi_{t}(dx),$

where

${\mathcal{R}}^{c,n}_{i}(t)=n-{\mathcal{R}}^{n}_{i}(t),\qquad\bar{\mathcal{R}}^{c,n}_{i}(t)=n^{-1}{\mathcal{R}}^{c,n}_{i}(t).$

Fix $t$ and $\varepsilon>0$. The function $x\mapsto\xi_{t}(x,\infty)$ has at
most countably many discontinuities. Hence one can find a finite sequence
$x_{*}=y_{0}<y_{1}<\ldots<y_{K}=x^{*}$ such that $[x_{*},x^{*}]$ contains the
compact support of $\phi^{\prime}$ (recall (4.34)),
$y_{k}-y_{k-1}<\varepsilon$, $k=1,\ldots,K$, and $\xi_{t}(\\{y_{k}\\})=0$,
$k=0,\ldots,K$. Because $y_{k}$ are not charged, both
$\bar{\xi}^{n}_{t}(y_{k-1},y_{k})$ and $\bar{\xi}^{n}_{t}[y_{k-1},y_{k}]$
converge a.s. to $\xi_{t}(y_{k-1},y_{k})$. Define the complementary rank by
${\rm rank}^{c}(i;x)=n-{\rm rank}(i;x)$, $x\in{\mathbb{R}}^{n}$, $i\in[n]$.
Then $0\leq{\rm rank}^{c}(i;x)\leq n-1$ and by (2.3),

${\rm rank}^{c}(i;x)=\\#\\{j:x_{j}>x_{i}\\}+\\#\\{j>i:x_{j}=x_{i}\\}.$

Denote $x_{i}=\hat{X}^{n}_{i}(t)$, $x=(x_{1},\ldots,x_{n})$, and for
$k=1,\ldots,K$,

$I_{k}=[y_{k-1},y_{k}),\qquad V_{k}=\\{i:x_{i}\in I_{k}\\},\qquad
P_{k}=\xi^{n}_{t}(I_{k}),\qquad Q_{k}=\xi^{n}_{t}[y_{k},\infty).$

Fix $k$, and for $i\in V_{k}$ let $j(i)={\rm rank}(i;(x_{l})_{l\in V_{k}})$.
This is a relabeling of $V_{k}$ according to the rank of its members within
$V_{k}$. With this notation we have

${\rm rank}^{c}(i;x)=Q_{k}+P_{k}-j(i),\qquad i\in V_{k}.$

Because $j(i)$ take all values between $1$ and $P_{k}$ as $i$ varies in
$V_{k}$, we have

$\sum_{i\in V_{k}}{\rm
rank}^{c}(i;x)^{\ell-1}=\sum_{j=1}^{P_{k}}(Q_{k}+P_{k}-j)^{\ell-1}=\sum_{j=0}^{P_{k}-1}(Q_{k}+j)^{\ell-1}.$

If $p_{n}$ is a sequence satisfying $n^{-1}p_{n}\to p\geq 0$ then
$n^{-1}\sum_{j=0}^{p_{n}}(j/n)^{\ell-1}\to\int_{0}^{p}z^{\ell-1}dz=\ell^{-1}p^{\ell}$.
Noting, as mentioned, that $n^{-1}P_{k}\to\xi_{t}(I_{k})$, we get

(4.47) $\frac{1}{n}\sum_{i\in
V_{k}}\bar{\mathcal{R}}^{c,n}_{i}(t)^{\ell-1}\to\frac{1}{\ell}(\xi_{t}(y_{k-1},\infty)^{\ell}-\xi_{t}(y_{k},\infty)^{\ell})=\frac{1}{\ell}\xi_{t}(I_{k})\mathfrak{S}(\xi_{t}(y_{k-1},\infty),\xi_{t}(y_{k},\infty)).$

Fix a sequence $\varepsilon_{m}>0$, $\varepsilon_{m}\to 0$. For each $m$, let
points as the above $y_{k}$, corresponding to $\varepsilon_{m}$, be denoted by
$y_{k}^{m}$. For $x\in[x_{*},x^{*}]$ let $x_{(m)}$ and $x^{(m)}$ be the unique
two points $y_{k-1}^{m}$ and $y_{k}^{m}$, respectively, for which
$x\in[y_{k-1}^{m},y_{k}^{m})$. This gives $x^{(m)}\to x$, $x^{(m)}>x$. By
right-continuity of $x\mapsto\xi_{t}(x,\infty)$, it follows that
$\xi_{t}(x^{(m)},\infty)\to\xi_{t}(x,\infty)$ for every $x\in[x_{*},x^{*}]$.
Similarly, by left continuity of $x\mapsto\xi_{t}[x,\infty)$, and $x_{(m)}\to
x$, $x_{(m)}\leq x$, we obtain $\xi_{t}(x_{(m)},\infty)\to\xi_{t}[x,\infty)$.
Let

$h_{(m)}(x)=\inf\\{\phi^{\prime}(z):z\in[x_{(m)},x^{(m)}]\\},\qquad
h^{(m)}(x)=\sup\\{\phi^{\prime}(z):z\in[x_{(m)},x^{(m)}]\\}.$

Summing over $k$ in (4.47), we have for every $m$,

(4.48) $\displaystyle\liminf_{n}\gamma^{n}(t)$
$\displaystyle\geq\frac{1}{\ell}\int_{[x_{*},x^{*}]}h_{(m)}(x)\mathfrak{S}(\xi_{t}(x_{(m)},\infty),\xi_{t}(x^{(m)},\infty))\xi_{t}(dx)$
(4.49) $\displaystyle\limsup_{n}\gamma^{n}(t)$
$\displaystyle\leq\frac{1}{\ell}\int_{[x_{*},x^{*}]}h^{(m)}(x)\mathfrak{S}(\xi_{t}(x_{(m)},\infty),\xi_{t}(x^{(m)},\infty))\xi_{t}(dx).$

It remains to show the expressions of the r.h.s. of (4.48) and (4.49) both
converge to the r.h.s. of (4.46) as $m\to\infty$. By bounded convergence, it
suffices that the integrands in both these integrals converge, for every $x$,
to $\phi^{\prime}(x)\mathfrak{S}(\xi_{t}[x,\infty),\xi_{t}(x,\infty))$. The
latter convergence holds because, by continuity of $\phi^{\prime}$,
$h^{(m)}(x)\to\phi^{\prime}(x)$ and $h^{(m)}(x)\to\phi^{\prime}(x)$, and
moreover $\xi_{t}(x_{(m)},\infty)\to\xi_{t}[x,\infty)$ and
$\xi_{t}(x^{(m)},\infty)\to\xi_{t}(x,\infty)$. This completes the proof. ∎

###  4.3. PDE satisfied by the limit

We take limits along an arbitrary convergent subsequence. By (4.37), (4.38),
Lemma 4.1 and Lemma 4.2, we have now established that every limit point $\xi$
satisfies (2.19).

Again, fix a convergent subsequence and denote its limit by $\xi$. Some
notation used below is as follows. For a right-continuous nondecreasing
function $V:{\mathbb{R}}\to[0,1]$ with $V(0-)=0$, let $V(dx)$ denote the
Stieltjes measure induced on ${\mathbb{R}}_{+}$. Let

$v(x,t)=\xi_{t}(x,\infty),\qquad\tilde{v}(x,t)=1-v(x,t)=\xi_{t}[0,x].$

These are right continuous functions for every $t$. Thus
$\tilde{v}(dx,t)=\xi_{t}(dx)$. Next, for $V$ as above, let the pure jump part
and continuous part be denoted, respectively, by

$V^{\rm jmp}(x)=\sum_{y\in[0,x]}\mathnormal{\Delta}V(y),\qquad V^{\rm
cts}(x)=V(x)-V^{\rm jmp}(x).$

######  Lemma 4.3.

$v$ is a weak solution to (2.18).

Proof. By equations (4.37) and (4.38) and Lemmas 4.1 and 4.2,

$\langle\phi,\xi_{t}\rangle=\langle\phi,\xi_{0}\rangle+\int_{0}^{t}\langle
b_{1}\phi^{\prime}+a\phi^{\prime\prime},\xi_{s}\rangle
ds+\mathnormal{\Gamma}(t).$

Write $\mathnormal{\Gamma}(t)$ as $\frac{b_{0}}{\ell}\int_{0}^{t}A(s)ds$ where

$A(t)=\int\phi^{\prime}(x)\mathfrak{S}(v(x-,t),v(x,t))\tilde{v}(dx,t).$

Let $\tilde{v}(x,t)=\tilde{v}^{\rm cts}(x,t)+\tilde{v}^{\rm jmp}(x,t)$ denote
the decomposition alluded to above. Let $U(x,t)=1-v(x,t)^{\ell}$. Then

$U(dx,t)=\ell v(x,t)^{\ell-1}\tilde{v}(dx,t).$

Hence

$U^{\rm cts}(dx,t)=\ell v(x,t)^{\ell-1}\tilde{v}^{\rm cts}(dx,t),\qquad U^{\rm
jmp}(dx,t)=\ell v(x,t)^{\ell-1}\tilde{v}^{\rm jmp}(dx,t).$

Now,

$A(t)=\int\phi^{\prime}(x)\ell v(x,t)^{\ell-1}\tilde{v}^{\rm
cts}(dx,t)+\int\phi^{\prime}(x)\frac{v(x-,t)^{\ell}-v(x,t)^{\ell}}{v(x-,t)-v(x,t)}\tilde{v}^{\rm
jmp}(dx,t).$

The second integral on the right is

$\sum_{x\in[0,\infty)}\phi^{\prime}(x)(v(x-,t)^{\ell}-v(x,t)^{\ell})=\sum_{x\in[0,\infty)}\phi^{\prime}(x)(U(x,t)-U(x-,t))=\int\phi^{\prime}(x)U^{\rm
jmp}(dx,t).$

Therefore, using $\phi^{\prime}(0)=0$,

$\displaystyle A(t)$ $\displaystyle=\int\phi^{\prime}(x)(U^{\rm
cts}(dx,t)+U^{\rm jmp}(dx,t))=\int\phi^{\prime}(x)U(dx,t)$
$\displaystyle=-\int\phi^{\prime\prime}(x)U(x,t)dx=\int\phi^{\prime\prime}(x)v(x,t)^{\ell}dx.$

Using integration by parts we have
$\langle\phi,\xi_{t}\rangle=\int_{{\mathbb{R}}_{+}}(v(x,t)-1)\phi^{\prime}(x)dx$.
Recalling $\phi^{\prime}=\tilde{\phi}$ and combining the above calculations,
we obtain

$\displaystyle\int\tilde{\phi}(x)(v(x,t)-1)dx$
$\displaystyle=\int\tilde{\phi}(x)(v(x,0)-1)dx+\int_{0}^{t}\int(b_{1}\tilde{\phi}^{\prime}(x)+a\tilde{\phi}^{\prime\prime}(x))(v(x,s)-1)dxds$
$\displaystyle\qquad+\frac{b_{0}}{\ell}\int_{0}^{t}\int\tilde{\phi}^{\prime}(x)v(x,t)^{\ell}dxds.$

Thus

$\displaystyle\int\tilde{\phi}(x)v(x,t)dx$
$\displaystyle=\int\tilde{\phi}(x)v(x,0)dx+\int_{0}^{t}\int(b_{1}\tilde{\phi}^{\prime}(x)+a\tilde{\phi}^{\prime\prime}(x))v(x,s)dxds+a\int_{0}^{t}\tilde{\phi}^{\prime}(0)ds$
$\displaystyle\qquad+\frac{b_{0}}{\ell}\int_{0}^{t}\int\tilde{\phi}^{\prime}(x)v(x,t)^{\ell}dxds.$

According to Definition 2.3, this shows that $v$ is a weak solution of (2.18)
once it is verified that $v\in\mathbb{L}^{\infty}_{\rm
loc}({\mathbb{R}}_{+},\mathbb{L}^{1}({\mathbb{R}}_{+}))$.

This property can otherwise be stated as
$\sup_{t\in(0,T]}\int_{0}^{\infty}\xi_{t}[x,\infty)dx<\infty$ a.s., for every
$T<\infty$. Invoking Skorohod’s representation theorem we may assume w.l.o.g.
that, along the chosen convergent subsequence, one has $\bar{\xi}^{n}\to\xi$
a.s. Arguing by contradiction, let there exist a sequence $t_{N}\in(0,T]$ and
an event of positive ${\mathbb{P}}$ measure on which
$\int_{0}^{\infty}\xi_{t_{N}}[x,\infty)dx\to\infty$. Now, for each fixed
$t\in(0,T]$, by Fatou’s lemma,

$\displaystyle{\mathbb{E}}\int_{0}^{\infty}\xi_{t}[x,\infty)dx$
$\displaystyle\leq\liminf_{n}{\mathbb{E}}\int_{0}^{\infty}\bar{\xi}^{n}_{t}[x,\infty)dx$
$\displaystyle\leq\liminf_{n}\frac{1}{n}\sum_{i\in[n]}{\mathbb{E}}\hat{X}^{n}_{i}(t)$
$\displaystyle\leq\sup_{n}\max_{i\in[n]}{\mathbb{E}}\|\hat{X}^{n}_{i}\|^{*}_{T}\leq
c<\infty,$

where $c=c(T)$ and the last assertion follows from Lemma 3.5.iii. This
contradicts the assumption, and hence the uniform $\mathbb{L}^{1}$ property
follows. ∎

###  4.4. Proof of main results

We can now prove Theorems 2.1, 2.4, 2.5, 2.6, 2.7 and Proposition 2.2.

Proof of Theorem 2.4. Having proved uniqueness in Lemma 3.2, we proceed to the
remaining assertions. Existence of a
$C^{\infty}({\mathbb{R}}_{+}\times(0,\infty))$ solution is well known; see
e.g. [22, Lemma 3.2, p. 68].

In the case $\mathfrak{f}(z)=c_{1}z-bz^{\ell}$ and
$v_{0}(x)=\xi_{0}(x,\infty)$, we have proved in Lemma 4.3 that the function
$v$ induced by any subsequential limit $\xi$ is a weak solution. By uniqueness
it must be equal to the smooth function mentioned above. Since $x\mapsto
v(x,t)$ is automatically nonincreasing, the function $u:=-v_{x}$ is
nonnegative and smooth. Using this in (2.17), an integration by parts gives
the first 3 parts of (1.1) as well as that $u$ integrates to $1$ for each
$t>0$.

As for the initial condition $u(\cdot,t)dx\to\xi_{0}(dx)$ in (1.1), again
using (2.17) with $-v_{x}=u$ gives
$\int\phi^{\prime}(x)u(x,t)dx=\int\phi^{\prime}(x)\xi_{0}(dx)+O(t)$ as $t\to
0$, which proves that the initial condition is satisfied. This shows that $u$
thus defined forms a classical solution to (1.1). ∎

Proof of Proposition 2.2. The fact that $u_{\rm stat}$ given in (1.2) is a
stationary solution is verified directly. As for uniqueness, for any solution
$u$, the function $v=\int_{\cdot}^{\infty}u(x)dx$ takes values in $[0,1]$, is
monotone, and satisfies

$\mathfrak{f}(v)^{\prime}+av^{\prime\prime}=0,\qquad
v(0)=1,\qquad\lim_{x\to\infty}v(x)=0.$

Integrating we obtain $\mathfrak{f}(v)+av^{\prime}=c_{1}$, where $c_{1}$ is a
constant. Since $\mathfrak{f}(0)=0$, we obtain by the condition at infinity
that $\lim_{x\to\infty}v^{\prime}(x)=c_{1}/a$, which can only hold if
$c_{1}=0$. Thus $v$ must satisfy the ODE

$v^{\prime}=-a^{-1}\mathfrak{f}(v),\qquad\text{ on }{\mathbb{R}}_{+},\qquad
v(0)=1.$

Since $\mathfrak{f}$ is Lipschitz on $[0,1]$, the result follows by ODE
uniqueness. ∎

Proof of Theorem 2.1. The existence of a classical solution has been shown in
the proof of Theorem 2.4 above. The property $\sup_{t\in(0,T]}\int
xu(x,t)dx<\infty$ follows from the fact that $v\in\mathbb{L}^{\infty}_{\rm
loc}({\mathbb{R}}_{+};\mathbb{L}^{1}({\mathbb{R}}_{+}))$. The uniqueness in
the class of $C^{2,1}$ functions satisfying the above uniform integrability
follows directly from uniqueness of weak solutions to (2.16), via integration
by parts. ∎

Proof of Theorem 2.5. The $C$-tightness of $\bar{\xi}^{n}$ stated in Lemma
3.6, the fact stated in Lemma 4.3 that for any limit point $\xi$,
$v(x,t):=\xi_{t}(x,\infty)$ is a weak solution of (2.18), and the uniqueness
stated in Lemma 3.1 imply that $\bar{\xi}^{n}$ possesses a deterministic weak
limit $\xi$, and that $\xi\in
C({\mathbb{R}}_{+},{\mathcal{M}}_{1}({\mathbb{R}}_{+}))$. The assertion that
$\xi$ is given by $\xi_{t}(dx)=u(x,t)dx$, where the latter is the unique
classical solution to (1.1), follows now from Lemma 4.3 and Theorem 2.4. ∎

Proof of Theorem 2.6. By Lemma 3.3.iii and Lemma 3.5.i, for $i\in[k]$,

(4.50)
$\hat{X}^{n}_{i}(t)=\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}(t)-\hat{S}^{n}_{i}(T^{n}_{i}(t))+\hat{b}^{n}_{1}t+\hat{\lambda}^{n}_{0}n^{-1/2}\int_{0}^{t}p_{n,{\mathcal{R}}^{n}_{i}(s)}ds+\hat{M}^{A,n}_{i}(t)+\hat{L}^{n}_{i}(t),$

and $\int\hat{X}^{n}_{i}(t)d\hat{L}^{n}_{i}(t)=0$. From Step 5 in the proof of
Lemma 4.1 we have that $T^{n}_{i}\to\iota$ in probability. Recall that
$(\hat{E}^{n}_{i},\hat{S}^{n}_{i})\Rightarrow(E_{i},S_{i})$ and that
$E_{i}-S_{i}$ are equal in law to $\sigma_{i}W_{i}$, where $W_{i}$ are
mutually independent standard BMs. Since, by assumption,
$(\hat{X}^{n}_{i}(0-))_{i\in[k]}\Rightarrow(X_{i}(0))_{i\in[k]}$, using the
dependence structure (2.7), we obtain
$(\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}-\hat{S}^{n}_{i}(T^{n}_{i}))\Rightarrow(X_{i}(0)+\sigma_{i}W_{i})_{i\in[k]}$,
where $(X_{i}(0))_{i\in[k]})$ is independent of $(W_{i})_{i\in[k]}$.

From (3.28) in the proof of Lemma 3.5 we have that $\hat{M}^{A,n}_{i}\to 0$ in
probability. From Lemma 3.5.iii, we have $C$-tightness of
$(\hat{X}^{n}_{i},\hat{L}^{n}_{i})_{i\in[k]}$. If we denote the integral term
in (4.50) by $\hat{K}^{n}_{i}$, we have tightness of the tuple
$(\hat{X}^{n}_{i},\hat{X}^{n}_{i}(0-)+\hat{E}^{n}_{i}-\hat{S}^{n}_{i}(T^{n}_{i}),\hat{K}^{n}_{i},\hat{L}^{n}_{i})_{i\in[k]}$,
and denoting a subsequential weak limit point by $(X_{i},X_{i}(0)+\sigma
W_{i},K_{i},L_{i})_{i\in[k]}$, one has, for $i\in[k]$,

$X_{i}(t)=X_{i}(0)+\sigma W_{i}(t)+b_{1}t+K_{i}(t)+L_{i}(t),\qquad\int
X_{i}dL_{i}=0.$

By uniqueness in law of the system of SDE (1.3), the result will be proved
once it is shown that $K_{i}(t)=b_{0}\int_{0}^{t}v(X_{i}(s))^{\ell-1}ds$.
Using Skorohod’s representation theorem we assume without loss that the
convergence along the subsequence is a.s. Thus, in view of (4.45), it suffices
to show that, along the subsequence, one has

(4.51)
$(\bar{\mathcal{R}}^{c,n}_{i})_{i\in[k]}\to(v(X_{i}(\cdot),\cdot))_{i\in[k]},$

a.s., in the uniform topology on $[t_{0},T]$, for any $0<t_{0}<T$.

To this end, recall that Theorem 2.5 establishes, in particular, that
$\bar{\xi}^{n}\to\xi$ in $D({\mathbb{R}}_{+},{\mathcal{M}}_{1})$, where
$\xi\in C({\mathbb{R}}_{+},{\mathcal{M}}_{1})$, and $\xi_{t}$ is atomless for
every $t>0$. Again it may be assumed that the convergence is a.s. Hence, with
$d_{\rm L}$ the Levy-Prohorov metric, $\|d_{\rm
L}(\bar{\xi}^{n},\xi)\|^{*}_{T}\to 0$ a.s. Because for any $t_{0}>0$ and
$x_{0}>0$, $v$ is uniformly continuous on $[0,x_{0}]\times[t_{0},T]$, this
gives

(4.52) $\sup_{t_{0}\leq t\leq T}\sup_{x\leq
x_{0}}|\bar{\xi}^{n}_{t}[x,\infty)-\xi_{t}[x,\infty)|\to 0,\quad\text{a.s.}$

By the atomless property of $\xi_{t}$, $t>0$, we know that

$\bar{k}_{n}:=\sup\\{\bar{\xi}^{n}_{t}(\\{x\\}):(x,t)\in[0,x_{0}]\times[t_{0},T]\\}\to
0,\quad\text{a.s.}$

Now, $\bar{\mathcal{R}}^{c,n}_{i}=n^{-1}(n-{\mathcal{R}}^{n}_{i})$ satisfies

$\bar{\xi}^{n}_{t}[\hat{X}^{n}_{i}(t),\infty)-\bar{k}_{n}\leq\bar{\mathcal{R}}^{c,n}_{i}(t)\leq\bar{\xi}^{n}_{t}[\hat{X}^{n}_{i}(t),\infty)+\bar{k}_{n},\qquad(x,t)\in[0,x_{0}]\times[t_{0},T],\
i\in[k].$

Using this in (4.52) shows

$\sup_{t_{0}\leq t\leq
T}|\bar{\mathcal{R}}^{c,n}_{i}(t)-v(\hat{X}^{n}_{i}(t),t)|1_{\\{\hat{X}^{n}_{i}(t)\leq
x_{0}\\}}\to 0,\quad\text{a.s.}$

Hence (4.51) follows by the a.s. convergence $\hat{X}^{n}_{i}\to X^{n}$ and
a.s. finiteness of $\|X_{i}\|^{*}_{T}$, $i\in[k]$, completing the proof. ∎

Proof of Theorem 2.7. The assumed exchangeability of $\hat{X}^{n}_{i}(0-)$ for
every $n$, along with the convergence in probability of their empirical law
$\bar{\xi}^{n}_{0-}$ to the deterministic limit $\xi_{0}$, imply the
convergence $(\hat{X}^{n}_{i}(0-))_{i\in[k]}\Rightarrow(X_{i}(0))_{i\in[k]}$
where the latter are mutually independent and each is $\xi_{0}$-distributed
[41, Proposition 2.2]. Hence the hypotheses of Theorem 2.6 hold. Because
Theorem 2.6 asserts that $(\hat{X}^{n}_{i},\hat{L}^{n}_{i})_{i\in[k]}$
converge in law to the solution of the system (1.3), in which $(W_{i})$ are
mutually independent and are independent of $(X_{i}(0))$, the additional
mutual independence of $(X_{i}(0))$ that has just been shown completes the
proof. ∎

Acknowledgment. Research of RA supported by ISF (grant 1035/20).

## References

  * [1] P. Agarwal and K. Ramanan.  Invariant states of hydrodynamic limits of randomized load balancing networks.  arXiv preprint arXiv:2008.08510, 2020.
  * [2] R. Aghajani, X. Li, and K. Ramanan.  The PDE method for the analysis of randomized load balancing networks.  Proceedings of the ACM on Measurement and Analysis of Computing Systems, 1(2):1–28, 2017.
  * [3] R. Aghajani and K. Ramanan.  The hydrodynamic limit of a randomized load balancing network.  The Annals of Applied Probability, 29(4):2114–2174, 2019.
  * [4] R. Atar, I. Keslassy, G. Mendelson, A. Orda, and S. Vargaftik.  Persistent-idle load-distribution.  Stochastic Systems, 10(2):152–169, 2020.
  * [5] R. Atar and D. Lipshutz.  Heavy traffic limits for join-the-shortest-estimated-queue policy using delayed information.  Mathematics of Operations Research, 46(1):268–300, 2021.
  * [6] S. Banerjee, A. Budhiraja, and B. Estevez.  Load balancing in parallel queues and rank-based diffusions.  arXiv preprint arXiv:2302.10317, 2023.
  * [7] S. Banerjee and D. Mukherjee.  Join-the-shortest queue diffusion limit in Halfin–Whitt regime: Tail asymptotics and scaling of extrema.  The Annals of Applied Probability, 29(2):1262–1309, 2019.
  * [8] S. Banerjee and D. Mukherjee.  Join-the-shortest queue diffusion limit in Halfin–Whitt regime: Sensitivity on the heavy-traffic parameter.  The Annals of Applied Probability, 30(1):80–144, 2020.
  * [9] E. Bayraktar, A. Budhiraja, and A. Cohen.  Rate control under heavy traffic with strategic servers.  The Annals of Applied Probability, 29(1):1–35, 2019.
  * [10] P. Billingsley.  Convergence of Probability Measures.  John Wiley & Sons, 2013.
  * [11] M. Bramson, Y. Lu, and B. Prabhakar.  Randomized load balancing with general service time distributions.  ACM SIGMETRICS performance evaluation review, 38(1):275–286, 2010.
  * [12] M. Bramson, Y. Lu, and B. Prabhakar.  Asymptotic independence of queues under randomized load balancing.  Queueing Systems, 71:247–292, 2012.
  * [13] M. Bramson, Y. Lu, and B. Prabhakar.  Decay of tails at equilibrium for FIFO join the shortest queue networks.  The Annals of Applied Probability, 23(5):1841–1878, 2013.
  * [14] A. Braverman.  Steady-state analysis of the join-the-shortest-queue model in the Halfin–Whitt regime.  Mathematics of Operations Research, 45(3):1069–1103, 2020.
  * [15] H. Chen and H.-Q. Ye.  Asymptotic optimality of balanced routing.  Operations Research, 60(1):163–179, 2012.
  * [16] K. L. Chung and R. J. Williams.  Introduction to Stochastic Integration, volume 2.  Springer, 1990.
  * [17] D. J. Daley and M. Miyazawa.  A martingale view of Blackwell’s renewal theorem and its extensions to a general counting process.  Journal of Applied Probability, 56(2):602–623, 2019.
  * [18] M. V. der Boor, S. C. Borst, J. S. Van Leeuwaarden, and D. Mukherjee.  Scalable load balancing in networked systems: A survey of recent advances.  SIAM Review, 64(3):554–622, 2022.
  * [19] P. Eschenfeldt and D. Gamarnik.  Join the shortest queue with many servers. The heavy-traffic asymptotics.  Mathematics of Operations Research, 43(3):867–886, 2018.
  * [20] S. Foss and A. L. Stolyar.  Large-scale join-idle-queue system with general service times.  Journal of Applied Probability, 54(4):995–1007, 2017.
  * [21] J. Gärtner.  On the McKean-Vlasov limit for interacting diffusions.  Mathematische Nachrichten, 137(1):197–248, 1988\.
  * [22] E. Godlewski and P.-A. Raviart.  Hyperbolic Systems of Conservation Laws.  Number 3-4. Ellipses, 1991.
  * [23] C. Graham.  Chaoticity on path space for a queueing network with selection of the shortest queue among several.  Journal of Applied Probability, 37(1):198–211, 2000.
  * [24] C. Graham.  Functional central limit theorems for a large network in which customers join the shortest of several queues.  Probability Theory and Related Fields, 131:97–120, 2005\.
  * [25] V. Gupta and N. Walton.  Load balancing in the nondegenerate slowdown regime.  Operations Research, 67(1):281–294, 2019.
  * [26] S. Halfin and W. Whitt.  Heavy-traffic limits for queues with many exponential servers.  Operations Research, 29(3):567–588, 1981.
  * [27] J. Jacod and A. Shiryaev.  Limit Theorems for Stochastic Processes.  Springer-Verlag, Berlin, 1987.
  * [28] E. V. Krichagina and M. I. Taksar.  Diffusion approximation for GI/G/1 controlled queues.  Queueing Systems, 12:333–367, 1992.
  * [29] D. Lacker.  On a strong form of propagation of chaos for McKean-Vlasov equations.  Electron. Commun. Probab., 23(45):1–11, 2018\.
  * [30] X. Liu, K. Gong, and L. Ying.  Large-system insensitivity of zero-waiting load balancing algorithms.  ACM SIGMETRICS Performance Evaluation Review, 50(1):101–102, 2022.
  * [31] M. J. Luczak and C. McDiarmid.  On the maximum queue length in the supermarket model.  The Annals of Probability, 34(2):493–527, 2006\.
  * [32] M. J. Luczak and J. Norris.  Strong approximation for the supermarket model.  The Annals of Applied Probability, 15(3):2038–2061, 2005.
  * [33] G. Mendelson and K. Xu.  Care: Resource allocation using sparse communication.  arXiv preprint arXiv:2206.02410, 2022.
  * [34] M. Mitzenmacher.  The power of two choices in randomized load balancing.  IEEE Transactions on Parallel and Distributed Systems, 12(10):1094–1104, 2001.
  * [35] D. Mukherjee, S. C. Borst, J. S. Van Leeuwaarden, and P. A. Whiting.  Universality of power-of-d load balancing in many-server systems.  Stochastic Systems, 8(4):265–292, 2018.
  * [36] A. Mukhopadhyay and R. R. Mazumdar.  Analysis of randomized join-the-shortest-queue (JSQ) schemes in large heterogeneous processor-sharing systems.  IEEE Transactions on Control of Network Systems, 3(2):116–126, 2015.
  * [37] P. Protter.  Stochastic Integration and Differential Equations. A New Approach. Springer, Berlin, 1990.
  * [38] M. Shkolnikov.  Large systems of diffusions interacting through their ranks.  Stochastic Processes and their Applications, 122(4):1730–1747, 2012.
  * [39] A. L. Stolyar.  Pull-based load distribution in large-scale heterogeneous service systems.  Queueing Systems, 80:341–361, 2015.
  * [40] A. L. Stolyar.  Pull-based load distribution among heterogeneous parallel servers: the case of multiple routers.  Queueing Systems, 85:31–65, 2017.
  * [41] A.-S. Sznitman.  Topics in propagation of chaos.  Lecture Notes in Mathematics, (1464):165–251, 1991.
  * [42] M. van der Boor, M. Zubeldia, and S. Borst.  Zero-wait load balancing with sparse messaging.  Operations Research Letters, 48(3):368–375, 2020.
  * [43] N. D. Vvedenskaya, R. L. Dobrushin, and F. I. Karpelevich.  Queueing system with selection of the shortest of two queues: An asymptotic approach.  Problemy Peredachi Informatsii, 32(1):20–34, 1996.
  * [44] G. Wolansky.  Comparison between two models of self-gravitating clusters: Conditions for gravitational collapse.  Nonlinear Analysis: Theory, Methods & Applications, 24(7):1119–1129, 1995.
  * [45] L. Ying, R. Srikant, and X. Kang.  The power of slightly more than one sample in randomized load balancing.  Mathematics of Operations Research, 42(3):692–722, 2017.
  * [46] Z. Zhao, S. Banerjee, and D. Mukherjee.  Many-server asymptotics for join-the-shortest queue in the super-Halfin-Whitt scaling window.  arXiv preprint arXiv:2106.00121, 2021.