2. Model and results

Let $(\Omega, {\mathcal{F}}, \{{\mathcal{F}}_t\}_{t\ge 0}, \mathbb{P})$ be a filtered probability space where all the random variables and processes are defined. We consider a queueing system with K service stations, each of which has its own dedicated arrival process and an infinite number of parallel servers. In addition, there is a flexible arrival stream, which can be served by any of the stations. Let $A_k = \{A_k(t)\colon t\ge 0\}$ be the dedicated arrival process at station $k=1,\dots,K$ , with arrival rate $\lambda_k$ and arrival times $\tau_{k,i}$ , $i\in {\mathbb N}$ , and $A_{0} = \{A_{0}(t)\colon t\ge0\}$ be the flexible arrival process, with arrival rate $\lambda_0$ and arrival times $\tau_{0,i}$ , $i\in {\mathbb N}$ . Assume that these arrival processes are mutually independent. Let $X_k = \{X_k(t)\colon t\ge 0\}$ be the process counting the number of jobs in station k for $k=1,\dots, K$ , and denote the counting process of jobs in the K service stations by $X= (X_1,\dots, X_K)$ .

For jobs initially in service in station k, let $\eta^0_{k,j}$ , $j=1,\dots, X_k(0)$ , be their remaining service times, and for the new arrivals from the stream $A_k$ , let $\eta_{k,i}$ , $i\in {\mathbb N}$ , be their service times. For the jobs from the arrival stream $A_{0}$ , let $\eta_{0,i}$ , $i\in {\mathbb N}$ , be their service times. We assume that the remaining service times $\{\eta^0_{k,j}\}_{k,j}$ are i.i.d. continuous random variables with cumulative distribution function (c.d.f.) $F_0$ , and the service times $\{\eta_{k,i}\}_{k,i}$ are all i.i.d. with c.d.f. F. Let us denote the upper tail probabilities by $F^{\mathrm{c}}_0=1-F_0$ and $F^{\mathrm{c}}=1-F$ for $F_{0}$ and F, respectively. Without loss of generality, assume that the mean of F is one.

We associate a ‘weight’ $\alpha_k>0$ with each station k in order to evaluate the status of station k by a score $\alpha_{k} X_{k}({\cdot})$ , and our routing policy depends on the scores. If the score $\alpha_{k} X_{k}$ of station k is lower than those of the other stations, then we route the newly arriving jobs to station k. More specifically, for the ith job from the arrival stream $A_{0}$ , at the arrival time $\tau_{0,i}$ , it is routed to station k if the score of station k is the lowest among the other scores, i.e.

(2.1) \begin{equation} \alpha_{k} X_{k}(\tau_{0,i}) < \min_{\ell \neq k } \alpha_\ell X_\ell(\tau_{0,i}).\end{equation}

If multiple stations have the lowest scores, then job i is routed to the station with the smallest number. For example, if stations 1 and 2 have the same score, i.e. $\alpha_{1} X_{1}(\tau_{0, i}) = \alpha_{2} X_{2}(\tau_{0, i})$ at the time $\tau_{0,i}$ of the arrival of job i, then job i is routed to station 1. Ties in the scores are resolved in this lexicographic way. For each $x= (x_1,\dots, x_K) \in {\mathbb R}^K_+$ , define an indicator function

(2.2) \begin{equation} \delta_k(x) = \begin{cases} 1 & \text{if} \ k = \min\{\,j \colon \alpha_{j} x_{j} = \min_{1\le \ell \le K} \alpha_{\ell} x_{\ell}\} , \\ 0 & \text{otherwise}, \end{cases}\end{equation}

and the set $\mathcal R_{k} \,:\!=\, \{ x \in \mathbb R^{K}_{+}\colon \delta_{k} (x) = 1\} $ for $k = 1, \ldots , K$ . Then we have $\delta_{k}(x) = \textbf{1}_{\mathcal R_{k}}(x)$ for $x \in \mathbb R_{+}^{K}$ , $k = 1, \ldots , K$ . Moreover, $\mathcal R_{j} \cap \mathcal R_{k} = \emptyset$ for $j \neq k$ , and $\cup_{k=1}^{K} \mathcal R_{k} = \mathbb R_{+}^{K}$ . Since $\delta_{k}(cx) = \delta_{k}(x) $ for every $x \in \mathbb R^{K}_{+}$ , $c > 0$ , we see that each set $\mathcal R_{k}$ forms a cone for $k = 1, \ldots , K$ .

Then, with these indicator functions, the process $X_k(t)$ can be described as

\begin{equation*} X_k(t) = \sum_{j=1}^{X_k(0)}{\mathbf 1}_{\{\eta^0_{k,j}>t\}} + \sum_{i=1}^{A_k(t)}{\mathbf 1}_{\{\tau_{k,i}+\eta_{k,i} >t\}} + \sum_{i=1}^{A_{0}(t)}\delta_{k}(X(\tau_{0,i}{-})){\mathbf 1}_{\{\tau_{0,i}+\eta_{0,i}>t\}}, \quad t\ge 0.\end{equation*}

We consider a sequence of such queueing systems, indexed by n, and the scaling limit. We first make the following assumption on the arrival rates.

Let $\bar{X}^n_k \,:\!=\, n^{-1}X^n_k$ for $k=1,\dots,K$ . The following functional law of large numbers (FLLN) holds. We use $D({\mathbb R}_+, {\mathbb R}^d)$ to denote ${\mathbb R}^d$ -valued càdlàg functions, endowed with the Skorokhod $J_1$ topology written as $(D({\mathbb R}_+, {\mathbb R}^d), J_1)$ (see [Reference Billingsley4] or [Reference Whitt26] for its definition). When $d=1$ , we write $D=D({\mathbb R}_+, {\mathbb R})$ and use $(D^d, J_1)$ to denote the d-fold product topology. When the space is restricted on the fixed time interval [0, T] for some $T > 0$ , we write $D_{[0, T]} \,:\!=\, D([0, T],{\mathbb R})$ and use $\big(D_{[0, T]}^{d}, J_{1}\big)$ to denote the d-fold product topology.

Theorem 2.1. Under Assumption 2.1, if there are deterministic constants $\bar{X}_k(0), k=1,\dots,K$ such that $\big(\bar{X}^n_1(0), \dots, \bar{X}^n_K(0)\big) \Rightarrow (\bar{X}_1(0), \dots, \bar{X}_K(0))$ in ${\mathbb R}_+^K$ as $n\to\infty$ , then

\begin{equation*} \big(\bar{X}^n_1,\dots,\bar{X}^n_K\big) \Rightarrow (\bar{X}_1, \dots, \bar{X}_K) \quad \mbox{in}\ (D^K, J_1) \ \mbox{as}\ n \to \infty, \end{equation*}

where

(2.3) \begin{equation} \bar{X}_k(t) = \bar{X}_k(0)F^{\mathrm{c}}_0(t) + \lambda_k \int_0^tF^{\mathrm{c}}(t-s)\,{\mathrm d} s, \quad t \ge 0. \end{equation}

In addition, $\bar{X}(t) \to \bar{X}^*= \big(\bar{X}^*_1, \dots, \bar{X}^*_K\big)$ as $t\to\infty$ , with $\bar{X}^*_k = \lambda_k$ for each $k=1,\dots,K$ .

Proof. Under the scaling in Assumption 2.1, the arrival rate $\lambda^n_0/\sqrt{n}\to \lambda_0$ implies that $\bar{A}^n_0(t) = A^n_0(t)/n \to 0$ in D in probability as $n\to\infty$ . Hence, the third term in $\bar{X}^n_k(t)$ ,

\begin{equation*} \frac{1}{n}\sum_{i=1}^{n\bar{A}^n_{0}(t)}\delta_{k}(\bar{X}^n(\tau_{0,i}{-})){\mathbf 1}_{\{\tau^n_{0,i}+\eta_{0,i}>t\}} \end{equation*}

tends to 0 in D in probability as $n\to\infty$ , since the summands are all bounded by one. Then, the convergence of $\bar{X}^n_k(t)$ reduces to the convergence of the first two terms, which is exactly the convergence of the LLN-scaled number of jobs in an infinite-server queueing model with an arrival process $A^n_k(t)$ and i.i.d. service times $\{\eta_{k,i}, i\ge 1\}$ . This follows from the existing result in the literature; see, e.g., [Reference Krichagina and Puhalskii16, Reference Pang, Talreja and Whitt21].

Observe that in the fluid limit, the ‘weighted’ shortest queue policy in (2.1) is irrelevant as indicated in (2.3), and the associated steady-state values $\bar{X}^*_k = \lambda_k$ . This is evident since the extra load $A^n_0(t)$ is of order $\sqrt{n}$ while the fluid scale is of order n. It is then necessary to consider the queueing dynamics in the diffusion scale.

Let $\hat{X}^n_k \,:\!=\, \sqrt{n} \big(\bar{X}^n_k - \bar{X}^*_k\big)$ for $k=1,\dots,K$ , and write $\hat{X}^n = \big(\hat{X}^n_1,\dots,\hat{X}^n_K\big)$ . Note that we center the process $\hat{X}^n_k$ by its equilibrium point. Then it is clear that if $\alpha_k = 1/\lambda_k$ , $\delta_k(X^n(t)) = \delta_k(\hat{X}^n(t))$ , $t \ge 0$ . We choose this specific value of $\alpha_k=1/\lambda_k$ to establish the following FCLT. Note that since all the service times are i.i.d. with mean one, the value $\lambda_k$ is the (offered) load at each station. Thus, the routing criterion chooses the station with minimum ratio of the current state and the steady state. When the $\lambda_k$ are equal for all k, the routing policy becomes the so-called ‘joining the shortest queue’ (JSQ) policy. Thus, the routing policy can be regarded as a ‘weighted’ JSQ policy with the weights being the reciprocal of the offered load.

Let $\hat{A}^n_k(t) \,:\!=\, (A^n_k(t) - \lambda^n_k t)/{\sqrt{n}}$ for $t\ge 0$ and $k=1,\dots, K$ , and let $\hat{A}^n_0(t) = A^n_0(t)/{\sqrt{n}}$ for $t\ge 0$ . We make the following assumption for these scaled arrival processes.

In the second condition, when the arrival processes $A_k$ are mutually independent renewal processes with the interarrival times having mean $\lambda^{-1}_k$ and variance $\sigma^2_k$ , if $A^n_k$ is defined by scaling the interarrival times by $n^{-1}$ , then the limit is given by $\hat{A}_k(t)= \sqrt{\lambda_k^3\sigma_k^2}\hat{B}_k(t)$ for some standard Brownian motion $\hat{B}_k(t)$ (see, e.g., [Reference Whitt26, Chapter 13.7]), $k = 1, \ldots , K$ , $t \ge 0$ .

We also make the following assumption on the initial condition.

Assumption 2.3. There exists a random vector $(\hat{X}_1(0), \dots, \hat{X}_K(0)) \in {\mathbb R}^K$ such that

\begin{equation*} \big(\hat{X}^n_1(0), \dots, \hat{X}^n_K(0)\big) \Rightarrow (\hat{X}_1(0), \dots, \hat{X}_K(0)) \quad \mbox{in}\ {\mathbb R}_+^K \ \mbox{as}\ n\to\infty. \end{equation*}

Under this condition, given the scaling of $\hat{X}_k$ , it is also clear that $\bar{X}_k(0) = \bar{X}^*_k=\lambda_k$ for each $k=1,\dots,K$ .

For the FCLT in Theorem 2.2, we also assume that $F_0(t) = F_{\mathrm{e}}(t)$ , where $F_{\mathrm{e}}(t)=\int_0^t F^{\mathrm{c}}(s)\,{\mathrm d}s$ for $t\ge 0$ is the equilibrium (stationary excess) distribution of F. Observe that for the fluid limit $\bar{X}(t)$ in (2.3), if $F_0=F_{\mathrm{e}}$ and $\bar{X}(0)=\bar{X}^*$ , then $\bar{X}(t) = \bar{X}^*$ for all $t\ge 0$ . Recall that we have assumed that the mean for the c.d.f. F is one. In the scaling limits, we have the following three mutually independent driving noises:

1. (i) some K-dimensional, independent Gaussian processes $ \hat{X}_{\cdot, 0} \,:\!=\, (\hat{X}_{1,0}, \ldots , \hat{X}_{K, 0})^{\prime}$ ;

2. (ii) another K-dimensional, independent Gaussian process, $\hat{X}_{\cdot, 1} \,:\!=\, (\hat{X}_{1,1}, \ldots , \hat{X}_{K, 1})^{\prime}$ ; and

3. (iii) the K-dimensional, independent Brownian motions $(\hat{A}_{1}, \ldots , \hat{A}_{K})$ with strictly positive variance rates $(c_{1}, \ldots , c_{K})$ from Assumption 2.2.

The processes $\hat{X}_{k,0}$ and $\hat{X}_{k,1}$ are independent continuous Gaussian processes, independent of $\hat{A}_{k}$ , with mean zero and covariance functions, for $t, t'\ge 0$ ,

(2.4) \begin{align} \text{Cov}(\hat{X}_{k,0}(t), \hat{X}_{k,0}(t')) & = \lambda_0\big(F_{\mathrm{e}}^{\mathrm{c}}(t\vee t') - F^{\mathrm{c}}_{\mathrm{e}}(t) F^{\mathrm{c}}_{\mathrm{e}}(t')\big),\qquad\qquad\qquad\qquad \end{align}

(2.5) \begin{align} \text{Cov}(\hat{X}_{k,1}(t), \hat{X}_{k,1}(t')) & = \lambda_k\int_0^{t\wedge t'}\big(F^{\mathrm{c}}(t\vee t'-s) - F^{\mathrm{c}}(t-s)F^{\mathrm{c}}(t'-s)\big)\,{\mathrm d}s.\end{align}

In addition, the processes $\hat{X}_{k,0}$ and $\hat{X}_{k,1}$ are independent of $\hat{X}_{k',0}$ and $\hat{X}_{k',1}$ , as well as $\hat{A}_{k'}$ , for every $k'\neq k$ . The process $\int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s)$ is also a continuous Gaussian process with covariance function, for $k=1,\dots,K$ ,

\begin{equation*}\text{Cov}\bigg(\int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s),\int_0^{t'} F^{\mathrm{c}}(t'-s)\,{\mathrm d}\hat{A}_k(s)\bigg) =c_k\int_0^{t\wedge t'}F^{\mathrm{c}}(t-s)F^{\mathrm{c}}(t'-s)\,{\mathrm d}s.\end{equation*}

As we shall see later in (3.2) and (3.5), respectively, $\hat{X}_{k,0}$ and $\hat{X}_{k,1}$ can be represented as a time-changed, Brownian bridge driven by another Brownian motion, and as a time-changed Kiefer process driven by a Brownian sheet, independent of Brownian motions.

Definition 2.1. (Weak solution to a system of stochastic Volterra integral equations.) We shall consider a weak solution to equation (2.6), that is, on some filtered probability space $(\Omega, \mathbb{P}, {\mathcal{F}},\{{\mathcal{F}}_t\}_{t\ge 0})$ , a continuous, adapted, K-dimensional process $\hat{\mathcal X}=(\hat{\mathcal X}_{1},\ldots,\hat{\mathcal X}_{K})^{\prime}$ , an independent, K-dimensional Gaussian process $\hat{\mathcal X}_{\cdot,0}\,:\!=\,(\hat{\mathcal X}_{1,0}, \dots, \hat{\mathcal X}_{K,0})'$ , another independent, K-dimensional Gaussian process $\hat{\mathcal X}_{\cdot,1}\,:\!=\,(\hat{\mathcal X}_{1,1}, \dots, \hat{\mathcal X}_{K,1})'$ determined by the covariance functions (2.4), and an independent, K-dimensional Brownian motion $\hat{\mathcal A}_{\cdot}$ with some variance rates $(c_{1}, \ldots, c_{K})$ satisfy almost surely, for $t \ge 0 $ ,

(2.6) \begin{align} \hat{\mathcal X}_k(t) & = \hat{\mathcal X}_k(0)F^{\mathrm{c}}_{\mathrm{e}}(t) + \hat{\lambda}_k F_{\mathrm{e}}(t) + \lambda_0\int_0^t\delta_k(\hat{X}(s))F^{\mathrm{c}}(t-s)\,{\mathrm d}s \nonumber \\ & \quad + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{\mathcal A}_k(s) + \hat{\mathcal X}_{k,0}(t) + \hat{\mathcal X}_{k,1}(t). \end{align}

Theorem 2.2 (FCLT.) Under Assumptions 2.1, 2.2, and 2.3, we have

\begin{equation*} \big(\hat{X}^n_1,\dots,\hat{X}^n_K\big) \Rightarrow (\hat{X}_1, \dots, \hat{X}_K) \quad \mbox{in}\ (D^K, J_1) \ \mbox{as}\ n \to \infty, \end{equation*}

where the limit $\hat{X} = (\hat{X}_1, \dots, \hat{X}_K)$ , together with the independent, continuous Gaussian processes $(\hat{X}_{1,0}, \ldots , \hat{X}_{K,0})$ , $(\hat{X}_{1,1}, \ldots , \hat{X}_{K,1})$ , and some Brownian motions $(\hat{\mathcal{A}}_{1}, \ldots , \hat{\mathcal{A}}_{K})$ (having the same law as that of $(\hat{A}_{1}, \ldots , \hat{A}_{K})$ in Assumption 2.2(ii), is a weak solution, unique in distribution, to the system of stochastic Volterra integral equations in (2.6).

For the convenience of notation, we write the limit $\hat{X} = (\hat{X}_1, \dots, \hat{X}_K)$ as

(2.7) \begin{align} \hat{X}_k(t) & = \hat{X}_k(0)F^{\mathrm{c}}_{\mathrm{e}}(t) + \hat{\lambda}_k F_{\mathrm{e}}(t) + \lambda_0\int_0^t\delta_k(\hat{X}(s))F^{\mathrm{c}}(t-s)\,{\mathrm d}s \nonumber \\ & \quad + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s) + \hat{X}_{k,0}(t) + \hat{X}_{k,1}(t). \end{align}

Remark 2.2. If we denote the last three Gaussian processes in (2.7) by

\begin{equation*} \hat{Y}_k(t) = \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s)+ \hat{X}_{k,0}(t) + \hat{X}_{k,1}(t)\end{equation*}

and assume F has a density f and that $F(0)=0$ , then we can write

(2.8) \begin{align} \hat{X}_k(t) & = \hat{X}_k(0)F^{\mathrm{c}}_{\mathrm{e}}(t) + \hat{\lambda}_k F_{\mathrm{e}}(t) + \lambda_0\int_0^t\delta_k(\hat{X}(s))F^{\mathrm{c}}(t-s)\,{\mathrm d}s + \hat{Y}_k(t) \nonumber \\ & = \int^{t}_{0}\bigg({-}\hat{X}_k(0)F^{\mathrm{c}}(s) + \hat{\lambda}_k F^{\mathrm{c}}(s) + \lambda_0\delta_k(\hat{X}(s)) - \lambda_0\int_0^s\delta_k(\hat{X}(s))f(s-u)\,{\mathrm{d}}u\bigg)\,{\mathrm{d}}s \nonumber \\ & \quad + \hat{Y}_k(t) \end{align}

for every $t \ge 0 $ . Observe that in the ‘drift’, there is not only a discontinuous term $\lambda_0 \delta_k(\hat{X}(t))$ , but also a memory of the process $\hat{X}(t)$ in the term $\lambda_0\int_0^t\delta_k(\hat{X}(s))f(t-s)\,{\mathrm{d}}s$ .

Remark 2.3. (Diffusion with discontinuous drift.) When the c.d.f. F is given by $F(t) =$ $1-{\mathrm{e}}^{-t}$ , $t \ge 0$ , it can be shown that the stochastic equation in (2.8) reduces to the diffusion with discontinuous drift studied in [Reference Chao5, Reference Krylov and Liptser18] (with some modification on the coefficients and also in the drift due to the absence of queue thresholds; see also the conjecture in [Reference Fleming and Simon10]), that is,

(2.9) \begin{equation} {\mathrm{d}}\hat{X}_k(t) = (\hat{\lambda}_k + \lambda_0\delta_k(\hat{X}(t)) - \hat{X}_k(t))\,{\mathrm{d}}t + \sqrt{\lambda_k + c_k^2}\,{\mathrm{d}}\tilde{B}_k(t) \end{equation}

for standard Brownian motion $(\tilde{B}_1, \dots, \tilde{B}_K)$ . The proof of this property is given in Appendix A.

3. Existence and uniqueness of a weak solution to the limiting stochastic Volterra integral equation with ranks

In this section we prove that the stochastic Volterra integral equation in (2.7) has a unique weak solution that has continuous paths.

Theorem 3.1. The stochastic Volterra integral equation (2.7) has a unique weak solution in $\mathbb{C}({\mathbb R}_+, {\mathbb R}^K)$ . Moreover, almost surely,

(3.1) \begin{equation} \int^{\infty}_{0}\sum_{k,\ell=1,\,k\neq\ell}^{K}\mathbf{1}_{\{\alpha_{k}\hat{X}_{k}(t)=\alpha_{\ell}\hat{X}_{\ell}(t)\}}\, {\mathrm d}t = 0, \end{equation}

where $\{\alpha_{k}\}_{1\le k \le K}$ are the associated weights in (2.1).

Before proceeding to the proof, we make the following observations on the Gaussian processes $\hat{X}_{\cdot,0}\,:\!=\,(\hat{X}_{1,0}, \dots, \hat{X}_{K,0})'$ and $\hat{X}_{\cdot,1}\,:\!=\,(\hat{X}_{1,1}, \dots, \hat{X}_{K,1})'$ . The process $\hat{X}_{k,0}$ is equivalent in distribution to a time-changed Brownian bridge $\hat{W}^0_k$ :

(3.2) \begin{equation} \hat{X}_{k,0}(t) = \lambda_k^{1/2}W^0(F_{\mathrm{e}}(t)) = \lambda_k^{1/2}\hat{W}^0_k(1-{\mathrm{e}}^{-t}).\end{equation}

Recall that the Brownian bridge $W^0_k$ is the unique strong solution to the one-dimensional stochastic differential equation (SDE) [Reference Karatzas and Shreve15]

\begin{equation*} {\mathrm{d}}\hat{W}^0_k(t) = -\frac{\hat{W}^0_k(t)}{1-t}\,{\mathrm{d}}t + {\mathrm{d}}\breve{B}_k(t),\end{equation*}

where $\breve{B}_k(t)$ is an independent standard Brownian motion. Thus we can write

(3.3) \begin{equation} \hat{X}_{k,0}(t) = \lambda_k^{1/2}\bigg({-}\int_0^{1-{\mathrm{e}}^{-t}}\frac{W^0_k(s)}{1-s}\,{\mathrm{d}}s + \breve{B}_k(1-{\mathrm{e}}^{-t})\bigg) = -\int_0^t\hat{X}_{k,0}(s)\,{\mathrm{d}}s + \lambda_k^{1/2}\breve{B}_k(1-{\mathrm{e}}^{-t})\end{equation}

for standard Brownian motions $\breve{B}_k(t)$ , independent of the $\hat{B}_k(t)$ in Assumption 2.2(ii). The second equality follows from a calculation using the change of variables for the integrable term.

Next, the Gaussian process $\hat{X}_{k,1}$ is equivalent in distribution to a time-changed Kiefer process, $\hat{\mathcal{K}}_k(t,x)$ , that is,

(3.4) \begin{equation} \hat{X}_{k,1}(t) = -\int_0^t\int_0^t\mathbf{1}_{s+x \le t}\,{\mathrm{d}}\hat{\mathcal{K}}_k(\lambda_k s,F(x)) = -\int_0^t\int_0^t\mathbf{1}_{s+x \le t}\,{\mathrm{d}}\hat{\mathcal{K}}_k(\lambda_k s,1-{\mathrm{e}}^{-x}).\end{equation}

Recall that, similar to a Brownian bridge, the Kiefer process $\hat{\mathcal{K}}_k(t,x)$ can be written in terms of Brownian sheets $\hat{W}_k(t,x)$ , i.e. $\hat{\mathcal{K}}_k(t,x)$ is the unique strong solution to the SDE [Reference Shorack and Wellner24]

\begin{equation*} \mathcal{K}(t,x) = -\int_0^x\frac{\mathcal{K}(t,y)}{1-y}\,{\mathrm d}y + \hat{W}_k(t,x).\end{equation*}

So we have

(3.5) \begin{align} \hat{X}_{k,1}(t) & = -\int_0^t\int_0^t\mathbf{1}_{s+x \le t}\,{\mathrm d}_{s,x} \bigg({-}\int_0^{1-{\mathrm{e}}^{-x}}\frac{\mathcal{K}(\lambda_k s,y)}{1-y}\,\mathrm{d}y + \hat{W}_k(\lambda_k s, 1-{\mathrm{e}}^{-x})\bigg) \nonumber \\ & = -\int_0^t\hat{X}_{k,1}(s)\,{\mathrm d}s + \int_0^t\int_0^t\mathbf{1}_{s+x \le t}\, {\mathrm d}\hat{W}_k(\lambda_k s, 1-{\mathrm{e}}^{-x}), \end{align}

where the $\hat{W}_k$ are mutually independent Brownian sheets, which are also independent of the standard Brownian motions $\hat{B}_k$ in (3.3) and $B_k$ in Assumption 2.2(ii).

Proof of Theorem 3.1.

Step 1: Construction. In order to construct the solution, let us analyze the solution to the stochastic Volterra integral equation (2.7) and derive equation (3.8) first. Consider a probability space $(\Omega, \mathbb{P}, {\mathcal{F}}, \{{\mathcal{F}}_t\}_{t\ge 0})$ under which the independent Brownian motions $(\hat{A}_{1},\ldots,\hat{A}_{K})$ in Assumption 2.2, independent Brownian bridges $(\hat{W}_{1}^{0},\ldots,\hat{W}_{K}^{0})$ in (3.2), and independent Kiefer processes $(\hat{\mathcal K}_{1},\ldots,\hat{\mathcal K}_{K})$ in (3.4) are defined. Assume for a moment that the solution $\hat{X}({\cdot})$ in (2.7) exists and, for $k=1,\dots, K$ , let us define

(3.6) \begin{align} \hat{R}_k(t) & \,:\!=\, \hat{\lambda}_{k}t + c_k\hat{B}_{k}(t) + \lambda_0\int^{t}_{0}\delta_{k}(\hat{X}(s))\,{\mathrm d}s, \end{align}

(3.7) \begin{align} \hat{S}_k(t) & \,:\!=\, \hat{X}_k(0)F^{\mathrm{c}}_{\mathrm{e}}(t) + \hat{X}_{k,0}(t) + \hat{X}_{k,1}(t). \ \ \quad \end{align}

Write $\hat{R}= (\hat{R}_1, \dots, \hat{R}_K)'$ and $\hat{S}= (\hat{S}_1, \dots, \hat{S}_K)'$ . Then $\hat{R}$ is a semimartingale with respect to the filtration $\{{\mathcal{F}}_t\}_{t\ge 0}$ , which is a Brownian motion with a random drift. The process $\hat{S}_{k}$ is a continuous Gaussian process that starts at $\hat{X}_{k}(0)$ for $k =1, \ldots , K$ . Recall the representations of $ \hat{X}_{k,0}(t)$ and $\hat{X}_{k,1}(t)$ in (3.2) and (3.4), respectively, using Brownian bridge $\hat{W}^0_k(t)$ and Kiefer process $\hat{\mathcal{K}}_k(t,x)$ (through the Brownian sheet $\hat{W}_k$ ). Also, let $\hat{X}_{\cdot, 0}=(\hat{X}_{1,0}, \dots, \hat{X}_{K,0})'$ and $\hat{X}_{\cdot, 1}=(\hat{X}_{1,1}, \dots, \hat{X}_{K,1})'$ . We similarly define, for $\hat{W}^0$ , $\hat{\mathcal{K}}$ and $\hat{W}$ . Note that since $\hat{B}(t)$ , $\hat{X}_{\cdot, 0}(t)$ , and $\hat{X}_{\cdot,1}(t)$ are mutually independent, we also have the mutual independence between $\hat{B}(t)$ , $\hat{W}^0$ , and $\hat{W}$ . Then, with (3.6) and (3.7), the stochastic integral equation (2.7) can be rewritten as

(3.8) \begin{equation} \hat{X}_{k}(t) = \int^{t}_{0}F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{R}_{k}(s) + \hat{S}_{k}(t) \end{equation}

for $k = 1, \ldots , K$ , $t \ge 0 $ , where the semimartingale $\hat{R}$ depends on $\hat{X}$ as in (3.6) and the continuous Gaussian process $\hat{S}$ is independent of $\hat{B}$ .

We shall first construct a weak solution. To construct a solution to this stochastic equation, we use the change of measure in such a way that the semimartingale $\hat{R}_{k}(s)$ becomes a Brownian motion under a new measure.

We consider a K-dimensional standard Brownian motion $\widetilde{\beta} \,:\!=\, (\widetilde{\beta}_{1},\ldots,\widetilde{\beta}_{K})$ and the independent Gaussian process $\hat{S}$ in (3.7) constructed from the independent Brownian bridges $\hat{W}_{\cdot}^{0}$ and the independent Kiefer processes $\hat{\mathcal K}_{\cdot}$ through the Brownian sheets, independent of $\widetilde{\beta}$ , on a filtered probability space $(\widetilde{\Omega},\widetilde{\mathcal F},\{\widetilde{\mathcal F}({t})\}_{t \ge 0},\widetilde{\mathbb P})$ , and we define $\xi \,:\!=\, (\xi_{1}, \ldots , \xi_{K})$ , $M \,:\!=\, (M_{1}, \ldots , M_{K})$ , with

(3.9) \begin{equation} \xi_{k}(t) \,:\!=\, \hat{S}_{k}(t) + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}\widetilde{\beta}_{k}(s), \quad M_{k}(t) \,:\!=\, \widetilde{\beta}_{k}(t) - \int^{t}_{0}\bigg(\frac{\hat{\lambda}_{k}}{c_{k}} + \frac{\lambda_0}{c_{k}}\delta_{k}(\xi(s))\bigg)\,{\mathrm d}s, \end{equation}

with the indicator function $\delta_{k}$ in (2.2) for $k=1,\ldots,K$ and $t \ge 0$ ,

\begin{equation*} Z(t) = \exp\Bigg(\sum_{k=1}^{K}\int^{t}_{0}\bigg(\frac{\hat{\lambda}_{k}}{c_{k}} + \frac{\lambda_0}{c_{k}}\delta_{k}(\xi(s))\bigg)\,{\mathrm d}\widetilde{\beta}_{k}(s) - \frac{1}{2}\sum_{k=1}^{K}\int^{t}_{0}\bigg(\frac{\hat{\lambda}_{k}}{c_{k}} + \frac{\lambda_0}{c_{k}}\delta_{k}(\xi(s))\bigg)^{2}\,{\mathrm d}s\Bigg). \end{equation*}

Here, M is a K-dimensional, drifted Brownian motion with at most linearly growing drifts.

Then the stochastic exponential Z is a continuous martingale under the probability measure $\widetilde{\mathbb P}$ , and hence, for a fixed $T > 0$ , we define a new probability measure $\widetilde{\mathbb Q}$ by

\begin{equation*} \frac{{\mathrm d}\,\widetilde{\mathbb Q}}{{\mathrm d}\,\widetilde{\mathbb P}} \bigg\vert_{\widetilde{\mathcal F}(T)} \,:\!=\, Z(T). \end{equation*}

Applying the Girsanov theorem ([Reference Karatzas and Shreve15, Theorem 3.5.1] for Brownian motions; see also, e.g., [Reference Nualart and Pardoux20, Proposition 1.6] for Brownian sheets), we see that M is a K-dimensional, standard Brownian motion, independent of $\hat{S}$ , under the probability measure $\widetilde{\mathbb Q}$ . Here, by the Girsanov theorem, we remove the drifts of the drifted Brownian motion but we do not shift the Brownian sheets that drive the independent Gaussian processes $\hat{S}$ . Thus, under $(\widetilde{\Omega},\widetilde{\mathcal F},\{\widetilde{\mathcal F}(t)\}_{t\ge0},\widetilde{\mathbb Q})$ , the adapted continuous process $\xi$ , the continuous Gaussian process $\hat{S}$ , and the Brownian motion M satisfy the equation

\begin{equation*} \xi_{k}(t) = \hat{S}_{k}(t) + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}M_{k}(s) + \int^{t}_{0}\hat{\lambda}_{k}F^{\mathrm{c}}(t-s)\,{\mathrm d}s + \lambda_0\int^{t}_{0}F^{\mathrm{c}}(t-s)\delta_{k}(\xi(s))\,{\mathrm d}s \end{equation*}

for $k = 1, \ldots, K$ , $0 \le t \le T$ . Thus, $\xi$ , M, and $\hat{S}$ in (3.9) satisfy the system (2.7) of stochastic Volterra integral equations for $0 \le t \le T$ under the new measure $\widetilde{\mathbb Q}$ . Since $T > 0$ is arbitrary, by the above construction there is a weak solution for $t \ge 0 $ to the system (2.7) of the stochastic Volterra integral equations.

Step 2: Uniqueness. The joint distribution of the weak solution $(\hat{X},\hat{A},\hat{X}_{\cdot,0},\hat{X}_{\cdot,1})$ to (2.7) is uniquely determined by the Girsanov change of measure as in the proof of [Reference Karatzas and Shreve15, Proposition 5.3.10], i.e. we firstly localize the problem by defining the sequence of the first passage times for $\hat{X}$ to the sphere of integer radii, centered at the origin; secondly we apply the Girsanov change of measure with those stopping times; and then we take the limits. Note that the Gaussian processes $(\hat{X}_{\cdot, 0}, \hat{X}_{\cdot, 1})$ are independent of the Brownian motion $\hat{A}$ . When the initial values $ \hat{X}(0) $ are randomized, we may determine the distribution as in [Reference Karatzas and Shreve15, Corollary 5.3.11].

Step 3: Proof of (3.1). To show (3.1), we first show that for the processes $\xi$ in (3.9) under $\widetilde{\mathbb P}$ ,

(3.10) \begin{equation} \int^{T}_{0}\sum_{k,\ell = 1,\, k\neq \ell}^{K}\mathbf{1}_{\{\alpha_{k}\xi_{k}(t)=\alpha_{\ell}\xi_{\ell}(t)\}}\, {\mathrm d}t = 0, \end{equation}

and then we once again apply the Girsanov theorem to show that (3.10) holds under $\widetilde{\mathbb Q}$ as in Step 1, and hence the weak solution $\hat{X}$ satisfies (3.1). Thus, it suffices to show (3.10) under $\widetilde{\mathbb P}$ where $\hat{S}$ and $\widetilde{\beta}$ are independent. Note that since the tail probability $F^{\mathrm{c}}$ and positive constants $c_{k}$ are deterministic, each integral $\int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}\widetilde{\beta}_{k}(s)$ is normally distributed with mean 0 and variance $\int^{t}_{0}(F^{\mathrm{c}}(t-s)c_{k})^{2}\,{\mathrm d}s$ for each $k=1,\ldots,K$ , independent of $\hat{S}$ under $\widetilde{\mathbb P}$ .

It follows that, for every $k \neq \ell$ and $t \ge 0$ , and for every fixed $(\theta_{1}, \ldots , \theta_{K}) \in \mathbb R^{K}$ ,

\begin{equation*} \widetilde{\mathbb P}\bigg(\alpha_{k}\bigg(\theta_{k} + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}\widetilde{\beta}_{k}(s)\bigg) = \alpha_{\ell}\bigg(\theta_{\ell} + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{\ell}\,{\mathrm d}\widetilde{\beta}_{\ell}(s)\bigg)\bigg) = 0, \end{equation*}

and hence, by the tower property of the conditional probability and by the independence, we have, for any $k \neq \ell$ and $t \in [0, T]$ ,

\begin{align*} & \widetilde{\mathbb P}(\alpha_{k}\xi_{k}(t) = \alpha_{\ell}\xi_{\ell}(t)) \\ & = \widetilde{\mathbb E}\bigg[\widetilde{\mathbb P}\bigg(\alpha_{k}\bigg(\hat{S}_{k}(t) + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}\widetilde{\beta}_{k}(s)\bigg) = \alpha_{\ell}\bigg(\hat{S}_{\ell}(t) + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{\ell}\,{\mathrm d}\widetilde{\beta}_{\ell}(s)\bigg) \mid \hat{S}(t)\bigg)\bigg] \\ & = \widetilde{\mathbb E}\bigg[\widetilde{\mathbb P}\bigg(\alpha_{k}\bigg(\theta_{k} + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{k}\,{\mathrm d}\widetilde{\beta}_{k}(s)\bigg) \\ & \qquad\qquad = \alpha_{\ell}\bigg(\theta_{\ell} + \int^{t}_{0}F^{\mathrm{c}}(t-s)c_{\ell}\,{\mathrm d}\widetilde{\beta}_{\ell}(s)\bigg)\bigg) \mid \theta_{k} = \hat{S}_{k}(t),\ \theta_{\ell} =\hat{S}_{\ell}(t)\bigg] = 0.\end{align*}

Here, $\widetilde{\mathbb E}$ represents the expectation under $\widetilde{\mathbb P}$ . Thus, we obtain (3.10) under $\widetilde{\mathbb P}$ , because

\begin{equation*} \widetilde{\mathbb E}\Bigg[\int^{T}_{0}\sum_{k,\ell=1,\,k\neq\ell}^{K} \mathbf{1}_{\{\alpha_{k}\xi_{k}(t)=\alpha_{\ell}\xi_{\ell}(t)\}}\,{\mathrm d}t\Bigg] = \int^{T}_{0}\sum_{k,\ell=1,\,k\neq\ell}^{K}\widetilde{\mathbb P}\big(\alpha_{k}\xi_{k}(t) = \alpha_{\ell}\xi_{\ell}(t)\big)\,{\mathrm d}t = 0 . \end{equation*}

By the reasoning in the previous paragraph, we claim the property (3.1).

Corollary 3.1. Recall the conic set $\mathcal R_{k}$ defined from the indicator function $\delta_{k}$ in (2.2) and $\delta_{k}({\cdot}) = \mathbf{1}_{\mathcal R_{k}}({\cdot})$ for $k = 1, \ldots , K$ . We denote the closure of $\mathcal R_{k}$ by $\overline{\mathcal R_{k}}$ for every k. The stochastic equation

\begin{align*} \hat{X}_k(t) & = \hat{X}_k(0)F^{\mathrm{c}}_{\mathrm{e}}(t) + \hat{\lambda}_k F_{\mathrm{e}}(t) + \lambda_0\int_0^t\mathbf{1}_{\overline{\mathcal R}_{k}}(\hat{X}(s))F^{\mathrm{c}}(t-s)\,{\mathrm d}s \\ & \quad + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s) + \hat{X}_{k,0}(t) + \hat{X}_{k,1}(t), \quad k = 1, \ldots , K \end{align*}

has a unique weak solution in $C(\mathbb R_{+}, \mathbb R^{K})$ .

4. Proof for the convergence to the limit

We have the representation

(4.1) \begin{align} \hat{X}^n_k(t) & = \hat{X}^n_k(0)F^{\mathrm{c}}_0(t) + \hat{\lambda}^n_k\int_0^tF^{\mathrm{c}}(t-s)\,{\mathrm d}s + \int_0^t\delta_k(\hat{X}^n(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^n_0(s) \nonumber \\ & \quad + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^n_k(s) + \hat{X}^n_{k,0}(t) + \hat{X}^n_{k,1}(t) + \hat{X}^n_{k,2}(t),\end{align}

where

\begin{align*} \hat{X}^n_{k,0}(t) & \,:\!=\, \frac{1}{\sqrt{n}\,}\sum_{j=1}^{X^n_k(0)}\Big({\mathbf 1}_{\eta^0_{k,j}>t} - F^{\mathrm{c}}_{\mathrm{e}}(t)\Big) = -\frac{1}{\sqrt{n}\,}\sum_{j=1}^{X^n_k(0)}\Big({\mathbf 1}_{\eta^0_{k,j}\le t} - F_{\mathrm{e}}(t)\Big), \\ \hat{X}^n_{k,1}(t) & \,:\!=\, \frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_k(t)}\Big({\mathbf 1}_{\tau^n_{k,i}+\eta_{k,i} >t} - F^{\mathrm{c}}\big(t-\tau^n_{k,i}\big)\Big), \\ \hat{X}^n_{k,2}(t) & \,:\!=\, \frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)}\delta_{k}\big(X^n\big(\tau^n_{0,i}{-}\big)\big) \Big({\mathbf 1}_{\tau^n_{0,i}+\eta_{0,i}>t} - F^{\mathrm{c}}\big(t-\tau^n_{0,i}\big)\Big).\end{align*}

Note that this assumption of $F_0=F_{\mathrm{e}}$ is essential since we are centering the process $\bar{X}^n_k$ by its equilibrium $\bar{X}_k^*$ , and in the derivation of the representation of $X^n_k(t)$ , the term $-\sqrt{n} \lambda_kF_0(t)$ cancels out the term $\sqrt{n}\lambda_k \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}s$ only if $F_0=F_{\mathrm{e}}$ .

The joint convergence in the following lemma follows directly from the existing results for each component in [Reference Krichagina and Puhalskii16, Reference Pang and Whitt22] for the heavy-traffic analysis of $\mathrm{G}/\mathrm{GI}/\infty$ queues and the mutual independence of the corresponding limits.

Lemma 4.1. Under Assumptions 2.1, 2.2, and 2.3,

\begin{align*} \bigg(\int_0^{\cdot}F^{\mathrm{c}}({\cdot}-s)\,{\mathrm d}\hat{A}^n_k(s), \hat{X}^n_{k,0}, \hat{X}^n_{k,1}\bigg)_{k=1,\dots,K} \Rightarrow \bigg(\int_0^{\cdot}F^{\mathrm{c}}({\cdot}-s)\,{\mathrm d}\hat{A}_k(s), \hat{X}_{k,0}, \hat{X}_{k,1}\bigg)_{k=1,\dots,K} \end{align*}

in $(D^{3K}, J_{1})$ as $n\to\infty$ , where the limits $\hat{X}_{k,0}$ and $\hat{X}_{k,1}$ are given in Theorem 2.2.

Proof. For each t, by conditioning on the filtration generated by the arrival process $A^n_0$ , we have

\begin{align*} \mathbb{E}\big[\big(\hat{X}^n_{k,2}(t)\big)^2\big] & = \frac{1}{n}\mathbb{E}\bigg[\sum_{i=1}^{A^n_0(t)}\delta_{k}\big(X^n\big(\tau^n_{0,i}-\big)\big) F\big(t-\tau^n_{0,i}\big)F^{\mathrm{c}}\big(t-\tau^n_{0,i}\big)\bigg] \\[4pt] & \le \frac{1}{\sqrt{n}\,}\mathbb{E}\bigg[\int_0^t F(t-s)F^{\mathrm{c}}(t-s)\,{\mathrm d} \frac{A^n_0(s)}{\sqrt{n}}\bigg] \to 0 \quad\mbox{as}\ n \to \infty, \end{align*}

where the convergence follows from Assumption 2.2 that ${A^n_0}/{\sqrt{n}} \Rightarrow \lambda_0e$ in D.

Next, we consider the increment, for $t, u\ge 0$ ,

(4.2) \begin{align} & \big|\hat{X}^n_{k,2}(t+u) - \hat{X}^n_{k,2}(t)\big| \nonumber \\[4pt] & \le \frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)}{\mathbf 1}_{t<\tau^n_{0,i}+\eta_{0,i}\le t+u} + \frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)}\big(F\big(t+u-\tau^n_{0,i}\big) - F\big(t-\tau^n_{0,i}\big)\big) \nonumber \\[4pt] & \quad + \frac{1}{\sqrt{n}\,}\sum_{i=A^n_{0}(t)+1}^{A^n_{0}(t+u)}\delta_{k}\big(X^n\big(\tau^n_{0,i}-\big)\big) \big|{\mathbf 1}_{\tau^n_{0,i}+\eta_{0,i}>t+u} - F^{\mathrm{c}}\big(t+u-\tau^n_{0,i}\big)\big|. \end{align}

For the first term, since it is non-decreasing in u, we have, for $\delta>0$ and $\varepsilon>0$ ,

\begin{align*} & \mathbb{P}\Bigg(\sup_{u\in[0,\delta]}\frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)} {\mathbf 1}_{t<\tau^n_{0,i}+\eta_{0,i}\le t+u} > \frac{\varepsilon}{3}\Bigg) \\[3pt] & \le \frac{9}{\varepsilon^2}\mathbb{E}\Bigg[\Bigg(\frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)} {\mathbf 1}_{t<\tau^n_{0,i}+\eta_{0,i}\le t+\delta}\Bigg)^2\Bigg] \\[3pt] & \le \frac{18}{\varepsilon^2}\mathbb{E}\Bigg[\Bigg(\frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)} {\mathbf 1}_{t<\tau^n_{0,i}+\eta_{0,i}\le t+\delta} - \big(F\big(t+\delta-\tau^n_{0,i}\big) - F\big(t-\tau^n_{0,i}\big)\big)\Bigg)^2\Bigg] \\[3pt] & \quad + \frac{18}{\varepsilon^2}\mathbb{E}\Bigg[\Bigg(\frac{1}{\sqrt{n}\,}\sum_{i=1}^{A^n_{0}(t)} \big(F\big(t+\delta-\tau^n_{0,i}\big) - F\big(t-\tau^n_{0,i}\big)\big)\Bigg)^2\Bigg] \\[3pt] & = \frac{18}{\varepsilon^2}\mathbb{E}\bigg[\frac{1}{\sqrt{n}\,}\int_0^t(F(t+\delta-s)-F(t-s)) (1 - (F(t+\delta-s)-F(t-s)))\,{\mathrm d}\frac{A^n_0(s)}{\sqrt{n}}\bigg] \\[3pt] & \quad + \frac{18}{\varepsilon^2}\mathbb{E}\bigg[\bigg(\int_0^t(F(t+\delta-s) - F(t-s))\, {\mathrm d}\frac{A^n_0(s)}{\sqrt{n}}\bigg)^2\bigg]. \end{align*}

Here, the first term converges to zero as $n \to \infty$ , and the second term satisfies

(4.3) \begin{align} & \frac{1}{\delta}\limsup_{N\to\infty}\sup_{t\in[0,T]}\mathbb{E}\bigg[\bigg(\int_0^t(F(t+\delta-s) - F(t-s))\, {\mathrm d}\frac{A^n_0(s)}{\sqrt{n}}\bigg)^2\bigg] \nonumber \\[3pt] & \qquad \le \frac{1}{\delta}\sup_{t\in[0,T]}\lambda_0^2\bigg(\int_0^t(F(t+\delta-s) - F(t-s))\, {\mathrm d}s\bigg)^2 \nonumber \\[3pt] & \qquad = \frac{1}{\delta}\sup_{t\in[0,T]}\lambda_0^2\bigg(\int_t^{t+\delta}F(s)\,{\mathrm d}s - \int_0^\delta F(s)\,{\mathrm d}s\bigg)^2 \le \lambda_0^2\delta, \end{align}

which converges to zero as $\delta \to 0$ .

The second term in (4.2) satisfies (4.3). The third term in (4.2) can be bounded by $\big(A^n_{0}(t+u) - A^n_{0}(t)\big)/{\sqrt{n}}$ , which is non-decreasing in u, so that the supremum over $u\in [0,\delta]$ is bounded by $\big(A^n_{0}(t+\delta) - A^n_{0}(t)\big)/{\sqrt{n}}$ . Then, by the convergence of ${A^n_0}/{\sqrt{n}} \Rightarrow \lambda_0 e$ in D, we obtain that, for small enough $\delta$ ,

\begin{equation*} \limsup_{N\to\infty}\mathbb{P}\bigg(\sup_{u\in[0,\delta]}\frac{1}{\sqrt{n}\,} \big(A^n_{0}(t+u) - A^n_{0}(t)\big) > \frac{\varepsilon}{3}\bigg) = 0. \end{equation*}

Thus, by [Reference Billingsley4, corollary on p. 83], we have shown that

\begin{align*} \frac{1}{\delta}\limsup_{N\to\infty}\sup_{t\in[0,T]}\mathbb{P}\bigg(\sup_{u\in[0,\delta]} \big|\hat{X}^n_{k,2}(t+u) - \hat{X}^n_{k,2}(t)\big| > \frac{\varepsilon}{3}\bigg) \to 0 \quad\mbox{as}\ \delta\to 0. \end{align*}

This completes the proof.

Proof of Theorem 2.2. We first observe that (4.1) can be rewritten as

(4.4) \begin{align} \hat{Z}^{n}_{k}(t) &\,:\!=\, \hat{X}^n_k(t) - \int_0^t\delta_k(\hat{X}^n(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^n_0(s) \nonumber \\ & = \hat{X}^n_k(0)F^{\mathrm{c}}_0(t) + \hat{\lambda}^n_k\int_0^tF^{\mathrm{c}}(t-s)\,{\mathrm d}s \nonumber \\ & \quad + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^n_k(s) + \hat{X}^n_{k,0}(t) + \hat{X}^n_{k,1}(t) + \hat{X}^n_{k,2}(t) \end{align}

for $k = 1,\ldots,K$ , $0 \le t \le T$ , and we write $\hat{Z}^{n} \,:\!=\, (\hat{Z}^{n}_{1},\ldots,\hat{Z}^{n}_{k})$ . Because of the tightness of the sequence $(\hat{X}^n_k(0),\hat{\lambda}^n_k)_{n\ge1}$ in $\mathbb R^{2}$ and the tightness of the sequence

\begin{equation*} \big(\big(\hat{A}^{n}_{0}({\cdot}), (\hat{A}^{n}_{k}({\cdot}), \hat{X}^{n}_{k,0}({\cdot}), \hat{X}^{n}_{k,1}({\cdot}), \hat{X}^{n}_{k,2}({\cdot})),\,k = 1,\ldots,K\big)\big)_{n\ge 1} \end{equation*}

in $\big(D_{[0, T]}^{4K+1}, J_{1}\big)$ , we claim the tightness of the sequence

(4.5) \begin{equation} \big(\hat{X}^{n}({\cdot}), \hat{Z}^{n}({\cdot}), \hat{A}_{0}^{n}({\cdot})\big)_{n\ge 1} = \big(\hat{X}^{n}_{1}({\cdot}),\ldots,\hat{X}^{n}_{K}({\cdot}),\hat{Z}^{n}_{1}({\cdot}),\ldots,\hat{Z}^{n}_{K}({\cdot}), \hat{A}^{n}_{0}({\cdot})\big)_{n\ge 1} \end{equation}

in $\big(D_{[0, T]}^{2K+1}, J_{1}\big)$ .

Now let us take a weak limit point $\big(\hat{X}^{\infty}({\cdot}),\hat{Z}^{\infty}({\cdot}),\hat{A}^{\infty}_{0}({\cdot})\big)$ of the sequence (4.5) in $\big(D_{[0, T]}^{2K+1}, J_{1}\big)$ . Without loss of generality, we may assume that the whole sequence may converge weakly to this limit point. By the Skorokhod representation theorem for the separable metric space $\big(D_{[0, T]}^{2K+1}, J_{1}\big)$ , we may take almost sure convergence

(4.6) \begin{equation} \lim_{n\to\infty}\big(\hat{X}^{n}(t),\hat{Z}^{n}(t),\hat{A}^{n}_{0}(t)\big) = \big(\hat{X}^{\infty}(t),\hat{Z}^{\infty}(t),\hat{A}^{\infty}_{0}(t)\big) \end{equation}

for all but countably many t on [0, T] by extending the probability space if necessary. The filtration for the corresponding probability space is taken to be the one generated by all these processes.

The limits of the right-hand side of (4.4), i.e. the limit of $\hat{Z}^{n}({\cdot}) = (\hat{Z}^{n}_{1}({\cdot}),\ldots,\hat{Z}^{n}_{K}({\cdot}))$ can be represented as

(4.7) \begin{equation} \hat{Z}^{\infty}_{k}(t) = \hat{X}_k(0)F^{\mathrm{c}}_0(t) + \hat{\lambda}_k\int_0^tF^{\mathrm{c}}(t-s)\,{\mathrm d}s + \int_0^t F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}_k(s) + \hat{X}_{k,0}(t) + \hat{X}_{k,1}(t) \end{equation}

for $1 \le k \le K$ , $0 \le t \le T$ , and it is continuous on [0, T]. Thus, the almost sure convergence of $\hat{Z}^{n}$ to $\hat{Z}^{\infty}$ in $\big(D_{[0, T]}^{K}, J_{1}\big)$ is uniform on [0, T]. With the same reasoning, the almost sure convergence $\lim_{n\to\infty}\hat{A}^{n}_{0}(t) = \lambda_{0}t$ is also uniform on [0, T], that is,

\begin{equation*} \lim_{n\to \infty}\sup_{0 \le t \le T}\big|\hat{A}^{n}_{0}(t) - \lambda_{0}t\big| = 0. \end{equation*}

Moreover, since the right-hand side of (4.7) is continuous almost surely in $t\in[0,T]$ , so is $\hat{Z}^{\infty}_{k}$ for each k. Note that the integrand $\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}({\cdot}-s)$ of $\int^{\cdot}_{0}\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}({\cdot}-s)$ ${\mathrm d}\hat{A}^{n}_{0}(s)$ is bounded and $F^{\mathrm{c}}$ is differentiable with a bounded derivative. Then, the sequence $\big\{\int^{t}_{0}\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^{n}_{0}(s),\, t \in [0, T], n \ge 1\big\}$ of absolutely continuous functions on [0, T] is uniformly equicontinuous. Thus, the almost sure limit

(4.8) \begin{equation} \Phi_{k}({\cdot}) \,:\!=\, \lim_{n\to\infty}\int^{\cdot}_{0}\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}({\cdot}-s)\, {\mathrm d}\hat{A}^{n}_{0}(s) = \lambda_{0}\lim_{n\to\infty}\int^{\cdot}_{0}\delta_{k}(\hat{X}^{n}(s{-})) F^{\mathrm{c}}({\cdot} - s)\,{\mathrm d}s \end{equation}

is uniformly converging by the Arzelà–Ascoli theorem, and hence it is also continuous in $t\in[0,T]$ . Thus, the limit $\hat{X}^{\infty}_{k}(t)$ is continuous in t, because of (4.4) and the continuity of $\hat{Z}^{\infty}$ . That is, $\hat{X}^{\infty}(t) = \hat{X}^{\infty}(t{-})$ for every $t\in[0,T]$ . Then we have the uniform convergence

\begin{equation*} \lim_{n\to\infty}\sup_{0 \le t \le T}|\hat{X}^{n}(t) - \hat{X}^{\infty}(t)| = 0, \end{equation*}

and $\hat{X}^{\infty}$ satisfies

(4.9) \begin{equation} \hat{X}^{\infty}_{k}(t) = \Phi_{k}(t) + \hat{Z}^{\infty}_{k}(t), \quad k=1,\ldots,K,\,t \ge 0, \end{equation}

where $\Phi_{k}$ and $\hat{Z}^{\infty}_{k}$ are given in (4.8) and (4.7), respectively.

Now we claim that each $\Phi_{k}$ is absolutely continuous with respect to Lebesgue measure by an application of the Riesz representation theorem for bounded linear functionals. Hence, by a similar argument of the change of measure as in the proof of Theorem 3.1, an analogue of (3.1) holds for $\hat{X}^{\infty}$ :

(4.10) \begin{equation} \int^{T}_{0}\sum_{k,\ell=1,k\neq\ell}^{K}\mathbf{1}_{\{\alpha_{k}\hat{X}^{\infty}_{k}(t)=\alpha_{\ell}\hat{X}^{\infty}_{\ell}(t)\}}\, {\mathrm d}t = \int^{T}_{0}\mathbf{1}_{\cup_{k=1}^{K}\partial\mathcal R_{k}}(\hat{X}^{\infty}(t))\,{\mathrm d}t = 0. \end{equation}

Here, $\partial\mathcal R_{k}$ is the boundary of $\mathcal R_{k}$ , i.e. $\partial\mathcal R_{k} = \{x\in\mathbb R_{+}^{K}\colon\alpha_{k}x_{k} = \alpha_{\ell}x_{\ell}$ for some $\ell\}$ for $k = 1, \ldots , K$ .

Since $\delta_{k}({\cdot})$ is an indicator function of a conic set in (2.2), the almost sure convergence (4.6) of $\hat{X}^{n}$ implies that the almost convergence

\begin{equation*} \lim_{n\to\infty}\delta_{k}(\hat{X}^{n}(s{-})) = \delta_{k}(\hat{X}^{\infty}(s{-})) = \delta_{k}(\hat{X}^{\infty}(s)) \end{equation*}

holds if $\hat{X}^{\infty}(s)$ is not on the boundary $\partial\mathcal R_{k}$ of the set $\mathcal R_{k}$ for $k=1,\ldots,K$ . Thus, thanks to (4.10), the almost sure limit $\Phi_{k} (t)$ of $\int^{t}_{0}\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^{n}_{0}(s)$ in (4.4) is

\begin{equation*} \lim_{n\to\infty}\int^{t}_{0}\delta_{k}(\hat{X}^{n}(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}\hat{A}^{n}_{0}(s) = \lambda_{0}\int^{t}_{0}\delta_{k}(\hat{X}^{\infty}(s{-}))F^{\mathrm{c}}(t-s)\,{\mathrm d}s, \quad 0 \le t \le T. \end{equation*}

Hence, we claim that (4.9) is reduced to

\begin{equation*} \hat{Z}^{\infty}_{k}(t) = \hat{X}_{k}^{\infty}(t) - \lambda_{0}\int^{t}_{0}\delta_{k}(\hat{X}^{\infty}(s))F^{\mathrm{c}}(t-s)\,{\mathrm d}s, \quad 0 \le t \le T,\,1\le k \le K. \end{equation*}

In other words, because of the representation in (4.7), the weak limit point $\hat{X}^{\infty}({\cdot})$ satisfies the stochastic equation (2.7). Thanks to the weak uniqueness in Theorem 3.1, the distribution of $\hat{X}^{\infty}({\cdot})$ is uniquely determined. Therefore, $\hat{X}^{n}({\cdot})$ converges weakly to $\hat{X}({\cdot})$ , which satisfies the stochastic equation (2.7) as $n\to\infty$ .

5. Concluding remarks

In this paper we have studied the non-Markovian parallel infinite-server model with a flexible stream under the scaling of its arrival rate being of order $\sqrt{n}$ . It will be interesting to investigate the case where the flexible stream has an arrival rate of order n, as also discussed in [Reference Fleming and Simon10]. It will also be interesting to consider such non-Markovian parallel many-server models with or without abandonment in the Halfin–Whitt regime [Reference Halfin and Whitt12, Reference Reed23]. The infinite-server model studied here can be used as an approximation for the many-server models in the underloaded regime, and can also be used as an approximation of the offered load process for the many-server models in the Halfin–Whitt regime. If the stationary distribution of the limiting process could be characterized or approximated, then it could be used for staffing decisions in each queue, as well as approximating the delay probabilities; see for example, the relevant work in [Reference Whitt27, Reference Whitt28].