2. Random splitting of Hawkes processes

A one-dimensional Hawkes process, $N=\{N(t),t\ge0\}$ , is a simple counting process with conditional intensity

(2.1) \begin{equation}\lambda(t)=\lambda_{0}+\sum_{j=1}^{N(t)}H(t-\tau_{j})=\lambda_{0}+\int_{0}^{t}H(t-s)\,{\textrm{d}} N(s),\end{equation}

where $\lambda_{0}>0$ is a constant, called the baseline intensity, $H\;:\; {\mathbb R}_{+}\to{\mathbb R}_{+}$ is the self-exciting function, and $\tau_{j}$ denotes the jth event time of N.

We consider the randomly split Hawkes processes $N_k=\{N_k(t),t\ge0\}$ , $k=1,\dots, d$ , defined as follows:

\begin{align*}N_{k}(t)=\sum_{j\ge1}\textbf{1}(\xi_{j}=k,\tau_{j}\le t),\end{align*}

where $\{\xi_{j},j\ge1\}$ is a sequence of i.i.d. random variables, independent of N, with the distribution

\begin{align*}\mathbb{P}(\xi_{j}=k)=p_{k} \in (0,1) \quad\text{and}\quad \sum_{k=1}^{d}p_{k}=1.\end{align*}

It is shown in [Reference Li and Pang23] that the splitting Hawkes process $(N_{k})_{k}$ is a multivariate Hawkes process with conditional intensity

(2.2) \begin{equation}\lambda_{k}(t)=p_{k}\lambda(t)=p_{k}\lambda_{0}+\sum_{k'=1}^{d}\int_{0}^{t}(p_{k}H(t-s))\,{\textrm{d}} N_{k'}(s).\end{equation}

Notice that the split processes $N_{k}$ are not independent.

We consider a sequence of Hawkes processes indexed by n, that is, $N^{(n)}$ is a Hawkes process for the nth system with intensity process $\lambda^{(n)}(\!\cdot\!)$ whose baseline intensity in (2.1) is $\lambda_{0}^{(n)}$ while the self-exciting function H stays the same. The splitting variables are denoted by $\{\xi_{j}^{(n)},j\ge1\}$ with distribution $(p^{(n)}_{k})_{k}$ .

Define the LLN and CLT-scaled processes for the splitting $\bigl(N^{(n)}_{k}\bigr)_{k}$ by

(2.3) \begin{equation}\bar{N}^{(n)}_{k}(t)=\dfrac{1}{n}N^{(n)}_{k}(nt) \quad\text{and}\quad\hat{N}^{(n)}_{k}(t)=\sqrt{n}\bigl(\bar{N}^{(n)}_{k}(t)-\mathbb{E}\bigl[\bar{N}^{(n)}_{k}(t)\bigr]\bigr),\end{equation}

respectively. It is easy to see (see e.g. [Reference Bacry, Delattre, Hoffmann and Muzy3]) that

(2.4) \begin{equation}\mathbb{E}\bigl[\bar{N}^{(n)}_{k}(t)\bigr]=\lambda_{0}^{(n)}p^{(n)}_{k}\int_{0}^{t}(1+\varphi*\mathfrak{e}(ns))\,{\textrm{d}} s ,\end{equation}

where $\varphi=\sum_{j\ge1} H^{*j}$ is the renewal function of H,

\begin{align*}f*g(x)=\int_{0}^{x}f(y)g(x-y)\,{\textrm{d}} y\end{align*}

denotes the convolution of f and g on ${\mathbb R}_{+}$ , and $\mathfrak{e}$ is the identity function.

The following FCLT for the splitting processes $\bigl(N^{(n)}_{k}\bigr)_{k}$ is proved in [Reference Li and Pang23].

Proposition 2.1. Suppose that

(2.5) \begin{equation}\lambda_{0}^{(n)}\to\lambda_{0} \quad\textit{and}\quad \bigl(p^{(n)}_{k}\bigr)\to (p_{k})_{k}\quad\textit{as}\;\text{$n\to\infty$}\quad \textit{and}\quad\|H\|_{1}=\int_{0}^{\infty}H(t)\,{\textrm{d}} t\in(0,1).\end{equation}

We have

\begin{align*}\bigl(\hat{N}^{(n)}_{k}\bigr)_{k}\Rightarrow (\hat{N}_{k})_{k}\quad\textit{in}\;\text{$({\mathbb D}^{d},J_{1})$} \quad \textit{as}\;\text{$ n\to\infty $,}\end{align*}

where $(\hat{N}_{k})_{k}$ is a d-dimensional Brownian motion with mean zero and covariance matrix

\begin{align*}\textrm{cov}(\hat{N}_{k}(t),\hat{N}_{k'}(s))=(t\wedge s)\cdot\biggl(\dfrac{\lambda_{0}p_{k}(\delta_{kk'}-p_{k'})}{1-\|H\|_{1}} + \dfrac{\lambda_{0}p_{k}p_{k'}}{(1-\|H\|_{1})^{3}}\biggr),\end{align*}

where $\delta_{kk'}=1$ if $k=k'$ and $\delta_{kk'}=0$ if $k\neq k'$ .

The process $\hat{N}_k$ admits the representations

(2.6) \begin{align}\hat{N}_{k}&= \dfrac{\lambda_{0}^{1/2}}{(1-\|H\|_{1})^{1/2}} \sqrt{p_{k}}W_{k}+p_{k}\dfrac{\lambda_{0}^{1/2}\|H\|_{1}}{(1-\|H\|_{1})^{3/2}} W \notag \\[5pt]&= \dfrac{\lambda_{0}^{1/2}}{(1-\|H\|_{1})^{1/2}} \hat{S}_{k}+ p_{k}\dfrac{\lambda_{0}^{1/2}}{(1-\|H\|_{1})^{3/2}}W,\end{align}

where $(W_{k})_k$ is a standard d-dimensional Brownian motion, $W=\sum_{k=1}^{d}\sqrt{p_{k}}W_{k}$ , and $\hat{S}_{k}=\sqrt{p_{k}}W_{k}- p_{k}W$ .

The condition $\|H\|_{1}\in(0,1)$ is referred to as the stability condition in the literature of Hawkes processes, under which a stationary version of Hawkes process exists and $\varphi*\mathfrak{e}(t)\to {{\|H\|_{1}}/{(1-\|H\|_{1})}}$ as $t\to\infty$ . The first representation in (2.6) is a direct consequence of [Reference Bacry, Delattre, Hoffmann and Muzy3, Theorem 2] for the vector-valued Hawkes process $\bigl(N^{(n)}_{k}\bigr)_{k}$ characterized by (2.2), while the second representation follows from [Reference Whitt32, Chapter 9.5], where $\hat{S}$ is also a d-dimensional Brownian motion, independent of W, with covariance function

\begin{align*} \textrm{cov}(\hat{S}_{k}(t),\hat{S}_{k'}(s))=p_{k}(\delta_{kk'}-p_{k'})(t\wedge s).\end{align*}

It is necessary that $\sum_{k}\hat{S}_{k}=0$ .

In our study of parallel server queues with a split Hawkes process as arrivals, we will also need the following alternative diffusion-scaled process by replacing the centering term $\mathbb{E}\bigl[\bar{N}^{(n)}_{k}(t)\bigr]$ in (2.3) with a linear function:

(2.7) \begin{equation} \check{N}^{(n)}_{k}(t) \;:\;eq \sqrt{n}\biggl(\bar{N}^{(n)}_{k}(t)-\dfrac{\lambda_{0}^{(n)}p^{(n)}_{k}t}{1-\|H\|_{1}}\biggr), \quad t \ge 0. \end{equation}

The proof of the proposition is given in the Appendix.

Proposition 2.2. Under the conditions in Proposition 2.1, if, in addition,

(2.8) \begin{equation}t^{1/2}\int_{t}^{\infty}H(s)\,{\textrm{d}} s\to0 \quad \textit{as}\;\text{$ t \to \infty$,}\end{equation}

then

(2.9) \begin{equation}\bigl(\check{N}^{(n)}_{k}\bigr)_{k} \Rightarrow (\hat{N}_{k})_{k}\quad\textit{in}\;\text{$ ({\mathbb D}^{d},J_{1})$} \quad \textit{as}\;\text{$ n\to\infty $,}\end{equation}

where $(\hat{N}_{k})_{k}$ is the same limit as given in Proposition 2.1.

3. Parallel single-server queues with split Hawkes arrival processes

In this model, there are d parallel servers and each server has its own queue. Arrivals in each queue come from the randomly split Hawkes process and are served in the FCFS discipline. Let $\{N(t)\;:\; t \ge 0 \}$ be the Hawkes process, arriving in the system, and let $(N_{k})_{k}$ be the splitting Hawkes arrival processes with the splitting mechanism as described in Section 2. Recall the baseline rate $\lambda_0$ , splitting probability $(p_k)_{k}$ , and self-exciting function H. For every $k=1,2,\ldots,d$ , let $\{\eta_{j,k}\}_j$ represent the service times of customers in the kth queue. We assume that $\{\eta_{j,k}\}_j$ are i.i.d. with finite mean $ m_{k}$ and variance $\sigma_{k}^2$ , and are independent of the Hawkes arrival process and the random splitting process. Denote the service rate of the kth queue by $\mu_{k}\;:\!=\; 1/m_{k}$ , and write $\mu = (\mu_k)_k$ .

We now give a definition of the model. For every $k=1,\ldots,d$ , let

\begin{equation*}V_{k}(m)=\sum_{j=1}^{m} \eta_{j,k}\end{equation*}

be the partial sum with the i.i.d. service times, and let $U_{k}(t)\;:\!=\; \sup\{m\ge 0\;:\; V_{k}(m)\le t\}$ be its right-continuous inverse. Then $U_{k}(t)$ is a renewal process representing the number of jobs that the server k can potentially complete by time t. Let $Z=(Z_{k})_{k}$ be the workload processes, and let $Q=(Q_{k})_{k}$ be the queue length processes, with $Q(0)=(Q_{k}(0))_{k}$ being the number of initial jobs in each queue. Under our assumptions, the processes N and (V, U) are independent, and we further assume that Q(0) is independent of the new arrivals $N=(N_{k})_{k}$ and the service processes. Then we have the following flow-balance equations for the dynamics at each queue: for $k=1,\ldots,d$ ,

(3.1) \begin{equation}\begin{aligned}Z_{k}(t)&= V_{k}(Q_{k}(0)+N_{k}(t))-B_{k}(t),\\[5pt]Q_{k}(t)&= Q_{k}(0)+N_{k}(t)-U_{k}(B_{k}(t)),\end{aligned}\end{equation}

where

\begin{equation*}B_{k}(t)=\int_{0}^{t}\textbf{1}(Q_{k}(s)>0)\,{\textrm{d}} s=\int_{0}^{t}\textbf{1}(Z_{k}(s)>0)\,{\textrm{d}} s\end{equation*}

is the busy-time process, and its paired idle-time process is given by

\begin{equation*}I_{k}(t)=t- B_{k}(t)=\int_{0}^{t}\textbf{1}(Q_{k}(s)=0)\,{\textrm{d}} s=\int_{0}^{t}\textbf{1}(Z_{k}(s)=0)\,{\textrm{d}} s.\end{equation*}

The processes $(Q_{k},Z_{k})_{k}$ take values in the non-negative orthant ${\mathbb R}_{+}^{2d}$ , and can be characterized by the following version of reflection maps (see [Reference Whitt32, Chapter 14.2]).

Definition 3.1. For any $x\in{\mathbb D}^{d}$ and any reflection matrix, i.e. a matrix with $\texttt{Q}_{ij}\ge0$ and $\sum_{i}\texttt{Q}_{ij}\le 1$ such that $(\texttt{Q})^{n}\to0$ as $n\to\infty$ , let the feasible regulator set be

\begin{align*}\Psi_{\texttt{Q}}(x)\equiv \bigl\{w\in {\mathbb D}_{\uparrow}^{d}\;:\; x+(\texttt{I}-\texttt{Q})w\ge0\bigr\},\end{align*}

where $\texttt{I}$ is the identity matrix and ${\mathbb D}_{\uparrow}^{d}$ is the subset of functions in ${\mathbb D}^{d}$ that are non-decreasing and non-negative in each coordinate. Then the reflection mapping is defined as

\begin{align*}(z,y)\;:\!=\; (\phi_{\texttt{Q}},\psi_{\texttt{Q}})(x)\;:\; {\mathbb D}^{d}\to {\mathbb D}^{2d},\end{align*}

where y is called the regulator component, given by

\begin{align*}y\;:\!=\; \psi_{\texttt{Q}}(x)\;:\!=\; \inf \Psi_{\texttt{Q}}(x)=\inf\{w\;:\; w\in\Psi_{\texttt{Q}}(x)\},\end{align*}

that is, for all i and t,

\begin{align*}y_{i}(t)\;:\!=\; \inf\{w_{i}(t)\in{\mathbb R}\;:\; w\in\Psi_{\texttt{Q}}(x)\},\end{align*}

and where z is called the content component, given by

\begin{align*}z=\phi_{\texttt{Q}}(x)=x+(\texttt{I}-\texttt{Q})y.\end{align*}

In the definition above, the matrix $\texttt{I}-\texttt{Q}$ is the direction of reflection (see [Reference Dai and Harrison9, Reference Harrison and Reiman18]), that is, whenever the boundary face $\{z\in{\mathbb R}^{d}_{+}\;:\; z_{j}=0\}$ is hit for some j, the process $w_{j}$ increases and causes an instantaneous displacement of z in the direction given by $\text{col}_{j}(\texttt{I}-\texttt{Q})$ , the jth column of $(\texttt{I}-\texttt{Q})$ ; $y=\psi_{\texttt{Q}}$ is the minimal element in $\Psi_{\texttt{Q}}$ such that $z\ge0$ . Moreover, the complementarity property holds ([Reference Whitt32, Theorem 14.2.3]), that is,

\begin{align*}\int_{0}^{\infty} z_{i} \,{\textrm{d}} y_{i}=0\quad\text{for all}\, \textit{i},\end{align*}

which also characterizes the regulator uniquely. It is proved in [Reference Whitt32, Theorems 14.2.5 and 14.2.7] that $\psi_{\texttt{Q}}$ is well-defined and Lipschitz under both the uniform norm and the Skorokhod $J_1$ topology.

In particular, if the reflection matrix $\texttt{Q}=0$ , then the angle of reflection is 0 (see [Reference Harrison, Landau and Shepp19]), and the reflection direction is orthogonal to the boundary, which is the case in our paper. Thus the reflection is also referred to as orthonormal. Letting $(\phi,\psi)\equiv(\phi_{0},\psi_{0})$ denote the operator in this case, and recalling that $\mathfrak{e}(t)\equiv t$ is the identity function on ${\mathbb R}_+$ , we can rewrite the processes in (3.1) in terms of Definition 3.1 as follows:

(3.2) \begin{equation}(Z,I)=(\phi,\psi)(V\circ(Q(0)+N)-\mathfrak{e}),\end{equation}

and the queue length process can be rewritten similarly.

We consider a sequence of such queueing systems in heavy traffic and indexed by the parameter n with $n\to \infty$ . The traffic intensity for the kth queue is given by

(3.3) \begin{equation}\rho^{(n)}_{k}=\dfrac{\lambda_{0}^{(n)}p^{(n)}_{k}}{(1-\|H\|_{1})\mu^{(n)}_{k}}, \quad k=1,\dots,d.\end{equation}

Define the following scaled processes:

(3.4) \begin{equation}\begin{gathered}\bar{V}^{(n)}(t)=\dfrac{1}{n}V^{(n)}([nt]),\quad\bar{U}^{(n)}(t)=\dfrac{1}{n}U^{(n)}(nt), \quad\bar{B}^{(n)}(t)=\dfrac{1}{n}B^{(n)}(nt), \\[5pt]\hat{V}^{(n)}(t)=\sqrt{n}\bigl(\bar{V}^{(n)}(t)-t m^{(n)}\bigr), \quad\hat{U}^{(n)}(t)=\sqrt{n}\bigl(\bar{U}^{(n)}(t)-t\mu^{(n)}\bigr),\end{gathered}\end{equation}

and

\begin{align*}\hat{Q}^{(n)}(t)=\dfrac{1}{\sqrt{n}} Q^{(n)}(nt),\quad \hat{Z}^{(n)}(t)=\dfrac{1}{\sqrt{n}} Z^{(n)}(nt), \quad t\ge 0.\end{align*}

We assume the following conditions on the service processes.

Assumption 3.1. Assume that the service times $\{\eta_{j,k}\}_j$ for $k=1,2,\ldots,d$ , satisfy

(3.5) \begin{equation} m^{(n)}_{k}=\mathbb{E}\bigl[ \eta_{1,k}^{(n)}\bigr]\to m_{k}\;=\!:\; \mu_{k}^{-1}\quad\textit{and}\quad\sigma^{(n)}_{k}=\sqrt{\textrm{var}\bigl( \eta_{1,k}^{(n)}\bigr)}\,\to\, \sigma_{k} \quad\textit{as}\;\text{$ n \to \infty$,}\end{equation}

and letting $F^{(n)}_{k}$ be the cumulative distribution function (CDF) of $ \eta_{1,k}^{(n)}$ , we have for every $\varepsilon>0$

(3.6)

The traffic intensity at each queue satisfies the following: for $k=1,\dots,d$ ,

(3.7)

Observe that by the definition in (3.3), under (3.7), we have

(3.8) \begin{equation} \lambda_{0}p_{k}m_{k}=\dfrac{\lambda_{0}p_{k}}{\mu_{k}}=1-\|H\|_{1}\quad\text{for each $k=1,2,\ldots,d$}.\end{equation}

Remark 3.1. Condition (3.6) is known as Lindeberg’s condition for the triangular array $\{\eta^{(n)}_{j,k}\}_{j,k}$ (see [Reference Billingsley5, Theorem 27.2]), under which one can show that as $n\to\infty$ ,

\begin{align*} \sup_{j\le n}\dfrac{1}{\sqrt{n}} \eta^{(n)}_{j,k}\Rightarrow 0, \end{align*}

and the processes in (3.4) satisfies

(3.9) \begin{equation}\begin{gathered}\bigl(\bar{V}^{(n)},\bar{U}^{(n)}\bigr)\rightarrow (\mu^{-1}, \mu) \mathfrak{e}\quad\text{u.o.c.} \, {in probability},\\[5pt]\bigl(\hat{V}^{(n)},\hat{U}^{(n)}\bigr)\Rightarrow(\hat{V},\hat{U}) \quad\text{in $ ({\mathbb D}^{2d}, J_{1})$,}\end{gathered}\end{equation}

where u.o.c. is short for ‘uniformly on every compact set on ${\mathbb R}_{+}$ ’, $\hat{V}$ is a d-dimensional Brownian motion with mean zero and covariance matrix

\begin{align*}\textrm{cov}(\hat{V}_{k}(t),\hat{V}_{k'}(s))=(t\wedge s)\sigma_{k}^{2}\delta_{kk'},\end{align*}

and $\hat{U}_{k}(t)=- \mu_{k}\hat{V}_{k}(\mu_{k}t)$ for $t\ge 0$ .

Remark 3.2. In addition to (3.5), if we assume

\begin{align*}\sqrt{n}\bigl(\lambda_{0}-\lambda_{0}^{(n)}\bigr)\to\hat{\lambda}_{0}, \quad \sqrt{n}\bigl(\mu_{k}- \mu^{(n)}_{k}\bigr)\to\hat{\mu}_{k},\quad \sqrt{n}\bigl(p_{k}-p^{(n)}_{k}\bigr)\to\hat{p}_{k}\end{align*}

for some $\hat{\lambda}_{0}, \hat{\mu}_{k},\hat{p}_{k}\in{\mathbb R}$ , as $n\to\infty$ , where $p_k$ is given in (2.5), then (3.7) holds with

\begin{align*}\hat{\rho}_{k}=\dfrac{\hat{\lambda}_{0}}{\lambda_{0}}+\dfrac{\hat{p}_{k}}{p_{k}}-\dfrac{\hat{\mu}_{k}}{\mu_{k}}.\end{align*}

The process limits remain the same.

We have the following FCLT for the diffusion-scaled processes $\bigl(\hat{Q}^{(n)},\hat{Z}^{(n)}\bigr)$ .

Theorem 3.1. Suppose that Assumption 3.1 and the conditions in (2.5) and (2.8) hold, and that $\hat{Q}^{(n)}(0)\Rightarrow\hat{Q}(0)\geq0$ . Then

(3.10) \begin{equation} \bigl(\hat{Q}^{(n)},\hat{Z}^{(n)}\bigr)\Rightarrow(\hat{Q},\hat{Z}) \quad\textit{in $ ({\mathbb D}^{2d}, J_1)$} \quad\textit{as}\;\text{$ n \to \infty$,}\end{equation}

with the limit

\begin{equation*}\hat{Q}= \phi(\hat{Q}(0) + \hat{W}-\hat{\theta} \mathfrak{e} ) \quad\textit{and}\quad\hat{Z}=m \hat{Q}=(m_{k}\hat{Q}_{k})_{k},\end{equation*}

where

\begin{align*}\hat{\theta}=\mu\hat{\rho}=(\mu_{k}\hat{\rho}_{k})_{k}\end{align*}

with $\hat{\rho}_{k}$ in (3.7), and $\hat{W}$ is a d-dimensional Brownian motion with covariance function

(3.11) \begin{equation}\begin{gathered}\textrm{cov}(\hat{W}_{k}(t),\hat{W}_{k'}(s))=(t\wedge s)\cdot\hat{R}_{kk'}, \\[5pt]\begin{aligned}\hat{R}_{kk'}&\;:\!=\; \dfrac{\lambda_{0}p_{k}p_{k'}}{(1-\|H\|_{1})^{3}}-\dfrac{\lambda_{0}p_{k}p_{k'}}{1-\|H\|_{1}}+\mu_k \bigl(1+ \mu_k^2\sigma_{k}^{2}\bigr) \delta_{kk'} \\[5pt]&= \mu_{k}\biggl(\dfrac{p_{k'}}{(1-\|H\|_{1})^{2}}-p_{k'}+\bigl(1+ \mu_k^2\sigma_{k}^{2}\bigr) \delta_{kk'} \biggr).\end{aligned}\end{gathered}\end{equation}

3.1. The stationary distribution of $\;\hat{\!Q}$

Observe that $\hat{Q}$ in Theorem 3.1 is a reflected Brownian motion on the orthant ${\mathbb R}_{+}^{d}$ with drift vector $\hat{\theta}$ , covariance matrix $\hat{R}$ , and reflection matrix $\texttt{I}$ . Since $\hat{\theta}_{k}>0$ for every $ k=1,2,\ldots,d$ , there exists a unique stationary distribution of $\hat{Q}$ .

Recall that a probability measure $\pi$ on ${\mathbb R}_{+}^{d}$ is called a stationary distribution for the reflected Brownian motion $\hat{Q}$ if, for every bounded Borel function f on ${\mathbb R}_{+}^{d}$ and every $t>0$ ,

\begin{align*}\int_{{\mathbb R}_{+}^{d}}\mathbb{E}_{x}[f(\hat{Q}(t))]\pi({\textrm{d}} x)=\int_{{\mathbb R}_{+}^{d}}f(x)\pi({\textrm{d}} x).\end{align*}

It is shown in [Reference Harrison and Reiman18] and [Reference Williams33] that the stationary distribution for a reflected Brownian motion in the orthant is equivalent to Lebesgue measure on ${\mathbb R}_{+}^{d}$ , and the occupation measures on faces are absolutely continuous with respect to the $(d-1)$ -dimensional Lebesgue measure $\nu_{k}(\!\cdot\!)$ on the face $F_{k}\;:\!=\; \{z \in {\mathbb R}_{+}^{d}\;:\; z_{k}=0\}$ , that is,

\begin{equation*}\pi({\textrm{d}} x)=q_{0}(x)\,{\textrm{d}} x\quad\text{and}\quad\mathbb{E}_{\pi}\biggl[\int_{0}^{t}\textbf{1}_{A}(\hat{Q}(s))\,{\textrm{d}}\hat{Y}_{k}(s)\biggr]=\dfrac{1}{2} t \int_{A} q_{k}(y)\nu_{k}({\textrm{d}} y)\quad\text{for all $A\in\mathscr{B}(F_{k})$.}\end{equation*}

Furthermore, $(q_{0};\, q_{1},\ldots, q_{d})$ in this paper jointly satisfy the basic adjoint relationship (BAR):

\begin{equation*} \int_{{\mathbb R}_{+}^{d}}\Biggl(\dfrac{1}{2}\sum_{i,j}\hat{R}_{ij}\dfrac{\partial^{2} f(x)}{\partial x_{i}\partial x_{j}}-\sum_{k}\hat{\theta}_{k}\dfrac{\partial f(x)}{\partial x_{k}}\Biggr) q_{0}(x)\,{\textrm{d}} x+\dfrac{1}{2}\sum_{k}\int_{F_{k}}\dfrac{\partial f}{\partial x_{k}}(x) q_{k}(x)\nu_{k}({\textrm{d}} x)=0\end{equation*}

for all $f\in C_{b}^{2}({\mathbb R}_{+}^{d})$ . It can be checked that the coefficients for $\hat{Q}$ do not satisfy the conditions of Propositions 8 and 9 in [Reference Dai and Harrison9] (see also [Reference Williams33]) which means the stationary density cannot be a product function. The numerical approach developed in [Reference Dai and Harrison9] to compute the stationary distribution is applicable to our model.

Noticing that the marginal distribution $\hat{Q}_{k}(\infty)$ for each k is exponentially distributed from the one-dimensional reflected Brownian motion results, we readily obtain

\begin{align*}\mathbb{E}[\hat{Q}_{k}(\infty)]=\dfrac{\hat{R}_{kk}}{2\hat{\theta}_{k}}\quad\text{and}\quad\textrm{var}(\hat{Q}_{k}(\infty))=\biggl(\dfrac{\hat{R}_{kk}}{2\hat{\theta}_{k}}\biggr)^{2}.\end{align*}

We apply the numerical algorithm in [Reference Dai and Harrison9] to a model with two queues and compute some characteristics of the stationary distribution. Similar to [Reference Dai and Harrison9], we compute the mean of the queueing limit and the correlation coefficients of $(\hat{W}_{1},\hat{W}_{2})$ and $(\hat{Q}_{1}(\infty),\hat{Q}_{2}(\infty))$ , and examine the effect of the splitting probability parameter $p_{k}$ . We fix $\lambda_{0}=1,\sigma_{1}=\sigma_{2}=1,\|H\|_{1}=0.3$ , and change the values of $p_{1}$ . Observe that by equation (3.8), the value of $\mu_{k}$ changes accordingly as $p_{k}$ changes. We also consider two scenarios of $(\hat{\rho}_{1},\hat{\rho}_{2})$ , taking values $(0.3,0.6)$ and $(1,1.5)$ . Notice that the correlation of $(\hat{W}_{1},\hat{W}_{2})$ from (3.11) is independent of the choice of $(\hat{\rho}_{1},\hat{\rho}_{2})$ by definition. Table 1 shows how the splitting probability $p_{1}$ and the parameters $(\hat{\rho}_{1},\hat{\rho}_{2})$ affect the correlation of $(\hat{Q}_{1}(\infty),\hat{Q}_{2}(\infty))$ .

Table 1. A numerical example to illustrate the stationary distribution of the limiting queueing process in a model with two queues.

4. Parallel single-server queues with split Hawkes arrival processes and abandonment

In this section we consider the parallel single-server queues as described in the previous section but with an additional feature of abandonment, that is, each customer joining the queue has a patience time and will leave the queue if the patience time runs out before entering service. The arrival and service processes are modeled in the same way as the previous section. Patience times are assumed to be independent of the arrival and service processes as well as the splitting process. Let $\{\vartheta_{j,k,p}\}_{j\in{\mathbb Z}}$ represent the patience times for every customers with CDF $G_{k}$ , $k=1,\dots,d$ . The dependence on k may be interpreted as the impact of each queue upon patience, since their services may have different distributions. Of course, one may also consider the homogeneous scenario with patience times having the same distribution in all the queues. In this model, variables and processes will be indexed with an additional subscript p to distinguish from the model in the previous section whenever necessary.

To give a rigorous definition of the model, we adapt the definition from [Reference Ward and Glynn30] and introduce the offered workload process

(4.1) \begin{equation}\tilde{Z}_{k,p}(t)=\Theta_{k,p}(Q_{k}(0))+\sum_{j\ge1}\eta_{j,k}\textbf{1}(\tau_{j,k}\le t, \tilde{Z}_{k,p}(\tau_{j,k}-)<\vartheta_{j,k,p}) - B_{k,p}(t)\end{equation}

for $t>0$ , where $\{\tau_{j,k}\;:\; j \ge 1\}$ are the arrival times of the split Hawkes process $N_k(t)$ ,

\begin{align*}\Theta_{k,p}(j)= \Theta_{k,p}(j-1)+ \eta_{-j,k}\textbf{1}(\Theta_{k,p}(j-1)<\vartheta_{-j,k,p})\quad\text{for $j\ge1$},\end{align*}

with $\Theta_{k,p}(0)=0$ representing the waiting time for the $(j+1)$ th customers initially in the kth queue. Note that $\tilde{Z}_{k,p}(0)=\Theta_{k,p}(Q_{k}(0))$ . The cumulative busy-time process is defined by

(4.2) \begin{equation}B_{k,p}(t)=\int_{0}^{t}\textbf{1}(\tilde{Z}_{k,p}(s)>0)\,{\textrm{d}} s.\end{equation}

Different from the offered workload process, the observed workload process is defined by

(4.3) \begin{align}Z_{k,p}(t)& = \tilde{Z}_{k,p}(t)+\sum_{j\ge1}\eta_{-j,k}\textbf{1}(j\le Q_{k}(0), \Theta_{k,p}(j-1)\ge \vartheta_{-j,k,p}>t )\notag \\[5pt]&\quad + \sum_{j\ge1} \eta_{j,k}\textbf{1}(\tilde{Z}_{k,p}(\tau_{j,k}-)\ge\vartheta_{j,k,p}>t-\tau_{j,k}\ge0),\end{align}

which retrospectively counts both the workload from customers that eventually receive service, $\tilde{Z}_{k,p}$ , and from those that are currently in a queue but eventually renege, $\textbf{1}(\cdot\ge\vartheta_{j,k,p}>\cdot)$ , in contrast to the prospective definition for $Z_{k}$ in (3.1).

Similarly, we define the observed queue length process as follows:

\begin{align*} Q_{k,p}(t) &=\sum_{j=1}^{Q_{k}(0)}\bigl(\textbf{1}(\Theta_{k,p}(j-1)<\vartheta_{-j,k,p}, t<\Theta_{k,p}(j))+\textbf{1}(\Theta_{k,p}(j-1)\ge\vartheta_{-j,k,p}> t)\bigr)\notag \\[5pt]&\quad +\sum_{j=1}^{N_{k}(t)}\bigl(\textbf{1}(\tilde{Z}_{k,p}(\tau_{j,k}-)<\vartheta_{j,k,p}, t<\tilde{Z}_{k,p}(\tau_{j,k}))+\textbf{1}(\tilde{Z}_{k,p}(\tau_{j,k}-)\ge\vartheta_{j,k,p}>t-\tau_{j,k})\bigr),\end{align*}

which counts the number of customers currently in the queue. Observe that the busy-time process in (4.2) also satisfies

\begin{align*}B_{k,p}(t)=\int_{0}^{t}\textbf{1}(Z_{k,p}(s)>0)\,{\textrm{d}} s=\int_{0}^{t}\textbf{1}(\tilde{Z}_{k,p}(s)>0)\,{\textrm{d}} s=\int_{0}^{t}\textbf{1}(Q_{k,p}(s)>0)\,{\textrm{d}} s\end{align*}

and the idle-time process is given by

\begin{equation*}I_{k,p}(t)\;:\!=\; t-B_{k,p}(t)=\int_{0}^{t}\textbf{1}(Q_{k,p}(s)=0)\,{\textrm{d}} s=\int_{0}^{t}\textbf{1}(\tilde{Z}_{k,p}(s)=0)\,{\textrm{d}} s.\end{equation*}

We also need the following definition generalizing the multidimensional reflection mapping in Definition 3.1, which is adapted from Appendix A.1 in [Reference Ward and Glynn30].

Furthermore, let u be the unique solution to the integral equation

(4.4) \begin{equation} u(t)+ \int_{0}^{t}\Gamma u(s)\,{\textrm{d}} s=x(t).\end{equation}

Define the mapping $\mathcal{M}_{\Gamma}\;:\; {\mathbb D}^{d}\to{\mathbb D}^{d}$ by $\mathcal{M}_{\Gamma}(x)=u$ . Then

\begin{align*}(z,y)=(\phi_{\texttt{Q}},\psi_{\texttt{Q}})(\mathcal{M}_{\Gamma}(x)),\end{align*}

where $\phi_{\texttt{Q}}$ and $\psi_{\texttt{Q}}$ are the content and reflection operators, respectively, in Definition 3.1.

By [Reference Ward and Glynn30, Lemma 3], $\mathcal{M}_{\Gamma}$ is Lipschitz-continuous with respect to uniform topology on every compact set. One can find from the integral equation (4.4) that

\begin{align*}u(t)&= \mathcal{M}_{\Gamma}(x)(t)\\[-1pt]&=x(0)+\int_{0}^{t}\,{\textrm{e}}^{\Gamma(s-t)}x({\textrm{d}} s)\\[-1pt]&=x(t)-\Gamma\int_{0}^{t}\,{\textrm{e}}^{\Gamma(s-t)}x(s)\,{\textrm{d}} s\\[-1pt]&= x(t)-\sum_{j\ge0}\Gamma^{j+1}\dfrac{(-1)^{j}}{j!}\int_{0}^{t}(t-s)^{j}x(s)\,{\textrm{d}} s.\end{align*}

In our model, $\Gamma=\text{diag} (\gamma)$ is a diagonal matrix with $\gamma=(\gamma_{k})_{k}\in{\mathbb R}^{d}$ on the diagonal, and we have

\begin{align*}u_{k}(t)=x_{k}(t)-\gamma_{k}\int_{0}^{t}\,{\textrm{e}}^{\gamma_{k}(s-t)}x_{k}(s)\,{\textrm{d}} s\end{align*}

for each coordinate process. In this case, we simply write $u=\mathcal{M}_{\gamma}(x)$ . If $x(t)=x_{0}+ a t+\sigma B(t)$ in Definition 4.1, where B is a d-dimensional standard Brownian motion, $a\in{\mathbb R}^{d}$ and $\sigma\in{\mathbb R}^{d\times d}$ , then $u=\mathcal{M}_{\gamma}(x)$ is the unique strong solution to the stochastic differential equation

\begin{equation*} {\textrm{d}} u(t)=(a-\gamma u(t))\,{\textrm{d}} t+\sigma \,{\textrm{d}} B_{t}=-\gamma u(t)\,{\textrm{d}} t+(a\,{\textrm{d}} t+\sigma \,{\textrm{d}} B_{t}),\quad u(0)=x_{0}.\end{equation*}

This is well-defined and called an Ornstein–Uhlenbeck (OU) process. Moreover, $z=\phi(\mathcal{M}_{\gamma}(x))$ is the regulated (reflected) OU process. Recall that $\phi=\phi_0$ with the matrix $\texttt{Q}\equiv 0$ .

As in the previous section, let $\bigl(Q^{(n)}_{p},Z^{(n)}_{p}\bigr)=\bigl(Q^{(n)}_{k,p}, Z^{(n)}_{k,p}\bigr)_{k}$ be the associated observed queue length and workload processes for the nth queueing system. We are interested in the diffusion-scaled observed processes in the heavy-traffic regime, and define

\begin{align*}\hat{Q}_{k,p}^{(n)}(t)\;:\!=\; \dfrac{1}{\sqrt{n}}Q_{k,p}^{(n)}(nt)\quad\text{and}\quad\hat{\tilde{Z}}_{k,p}^{(n)}(t)\;:\!=\; \dfrac{1}{\sqrt{n}}\tilde{Z}_{k,p}^{(n)}(nt).\end{align*}

We make the following assumption on $G^{(n)}_{k}$ for the variables $\vartheta^{(n)}_{j,k,p}$ .

Assumption 4.1. Assume that $G^{(n)}_{k}$ is continuous, and for some $\gamma_{k}>0$ ,

\begin{align*}\hat{G}^{(n)}_{k}(t)\;:\!=\; \sqrt{n}G^{(n)}_{k}(\sqrt{n}t)\rightarrow \gamma_{k}t\quad\textit{as}\;\text{$ n \to \infty$,}\end{align*}

where the convergence holds uniformly on compacts (u.o.c.) over ${\mathbb R}_+$ .

Note that in the case $\vartheta^{(n)}_{j,k,p}\;:\!=\; n\times\vartheta_{j,k,p}$ , that is, $G^{(n)}_{k}(t)=G_{k}(t/n)$ for some common CDF $G_{k}$ , the assumption above is equivalent to $G_{k}$ being differentiable at 0 with $\gamma_{k}=G'_{k}(0)$ . This is the so-called critical case studied in [Reference Ward and Glynn30].

Theorem 4.1. Suppose that Assumptions 3.1–4.1 and the conditions in (2.5) and (2.8) hold, and that $\hat{Q}^{(n)}(0)\Rightarrow\hat{Q}(0)\geq0$ . Then

\begin{equation*} \bigl(\hat{Q}^{(n)}_p,\hat{Z}^{(n)}_p\bigr)\Rightarrow(\hat{Q}_p,\hat{Z}_p) \quad\textit{in}\;\text{$ ({\mathbb D}^{2d}, J_1)$} \quad\textit{as}\;\text{$ n \to \infty$,}\end{equation*}

with the limit

\begin{align*}\hat{Q}_p=\phi (\mathcal{M}_{\gamma}(\hat{Q}(0) + \hat{W}-\hat{\theta} \mathfrak{e} ))\quad\textit{and}\quad \hat{Z}_p= \mu^{-1} \hat{Q}_p,\end{align*}

where $\hat{\theta}$ and $\hat{W}$ are as given in Theorem 3.1.

Here $\hat{Q}_p$ is a d-dimensional reflected OU process with initial value $\hat{Q}(0)$ . Although the limit in Theorem 4.1 formally reduces to that in Theorem 3.1 if $\gamma=0$ in Assumption 4.1, it is worth mentioning that the proof of Theorem 4.1 relies on Theorem 3.1, and in the case $\gamma>0$ ,

\begin{align*}(Z_{p},I_{p})\neq(\phi,\psi)(\mathcal{M}_{\gamma}(V\circ(Q(0)+N)-\mathfrak{e})),\end{align*}

in contrast to (Z, I) in (3.2).

5. A parallel server model with multivariate Hawkes arrivals

As noted in Section 2, a splitting Hawkes process is a multivariate Hawkes process with a particular conditional intensity process given by (2.2). In this section we present the limits for a parallel server model with a multivariate Hawkes arrival process, defined as a simple point process with conditional intensity process given by

\begin{equation*}\lambda_{i}(t)=\lambda_{i,0}+\sum_{j=1}^{d}\int_{0}^{t}H_{ij}(t-s)\,{\textrm{d}} N_{j}(s), \quad i =1,\dots, d, \quad t\ge 0,\end{equation*}

where $\lambda_{0}=(\lambda_{i,0})_{i}$ is the baseline intensity vector, and $H=(H_{ij})_{ij}$ is the kernel matrix function with $H_{ij}\;:\; {\mathbb R}_+\to {\mathbb R}_+$ . The non-explosion criterion in [Reference Bacry, Delattre, Hoffmann and Muzy3] is given by $ \int_{0}^{t}H_{ij}(s)\,{\textrm{d}} s<\infty$ , for all $ t>0$ , under which the point process is well-defined. Given the multivariate Hawkes process as arrivals and the service process as defined in Section 3 in the parallel single-server queues, the queue length process and workload process can be defined in the same way.

We consider a sequence of such queueing systems in heavy traffic with index n, where the baseline intensity is $\lambda^{(n)}$ and the kernel function is H.

We note in particular that the notation differs from Section 3: $\lambda_0$ is a vector and H is a matrix.

Under Assumption 5.1, it is shown (see [Reference Bacry, Delattre, Hoffmann and Muzy3, Theorem 2] and also Proposition 2.2) that

(5.1) \begin{equation}\check{N}^{(n)}(t)=\sqrt{n}\biggl(\dfrac{1}{n}N^{(n)}(nt)-(\texttt{I}-\|H\|_{1})^{-1}\lambda^{(n)}_{0}t\biggr)\Rightarrow \hat{N}(t)\end{equation}

in $({\mathbb D}^{d},J_{1})$ , where $\hat{N}$ is a d-dimensional Brownian motion with covariance function

\begin{equation*}\textrm{cov}(\hat{N}_{i}(t),\hat{N}_{j}(s))=(t\wedge s)\cdot\text{ent}_{ij}\bigl((\texttt{I}-\|H\|_{1})^{-1}\cdot \text{diag} ((\texttt{I}-\|H\|_{1})^{-1}\lambda_{0})\cdot(\texttt{I}-\|H\|_{1}^{\top})^{-1}\bigr)\end{equation*}

for $t, s \ge 0$ , where the superscript $\top$ denotes the transpose of a matrix, comparing it with Proposition 2.1 Note that in this model setup, the kernel matrix H is assumed to be independent of n. The traffic intensity for the kth queue in the nth system is given by (compare with (3.3))

\begin{equation*}\rho^{(n)}_{k}=\text{ent}_{k}\bigl((\texttt{I}-\|H\|_{1})^{-1}\lambda^{(n)}_{0}\bigr)\big/ \mu^{(n)}_{k}.\end{equation*}

The heavy-traffic condition will then imply that

\begin{align*}(\texttt{I}-\|H\|_{1})^{-1}\lambda_{0}=\mu.\end{align*}

Under Assumptions 3.1 and 5.1 together with the new traffic intensity above, and $\hat{Q}^{(n)}(0)\Rightarrow\hat{Q}(0)\ge0$ , we obtain

\begin{equation*}\bigl(\hat{Q}^{(n)},\hat{Z}^{(n)}\bigr)\Rightarrow(\hat{Q},\hat{Z}) \quad\text{in $ ({\mathbb D}^{2d}, J_1)$} \quad\text{as $ n \to \infty$,}\end{equation*}

with the limit

\begin{equation*}\hat{Q}= \phi(\hat{Q}(0) + \hat{W}-\hat{\theta} \mathfrak{e} ) \quad\text{and}\quad\hat{Z}=\mu^{-1} \hat{Q}=\bigl(\mu^{-1}_{k}\hat{Q}_{k}\bigr)_{k},\end{equation*}

where

\begin{align*}\hat{\theta}=\mu\hat{\rho}=(\mu_{k}\hat{\rho}_{k})_{k},\end{align*}

and $\hat{W}$ is a d-dimensional Brownian motion with the following covariance function: for $t,s \ge 0$ ,

(5.2) \begin{equation}\begin{gathered}\textrm{cov}(\hat{W}(t),\hat{W}(s))=(t\wedge s)\cdot \hat{R}_{kk'}, \\[5pt](\hat{R}_{kk'})_{kk'}=(\texttt{I}-\|H\|_{1})^{-1}\cdot\text{diag} (\mu)\cdot\bigl(\texttt{I}-\|H\|_{1}^{\top}\bigr)^{-1}+\text{diag}\bigl(\mu^{3}\sigma^{2}\bigr).\end{gathered}\end{equation}

For the parallel server queueing model with abandonment as described in Section 4, suppose the arrival process is also a multivariate Hawkes process as described above. Let $Q^{(n)}_{p}$ and $Z^{(n)}_{p}$ be the queue length process and the workload process for the model, respectively, and suppose that Assumptions 3.1, 4.1, and 5.1 hold, and that $\hat{Q}^{(n)}(0)\Rightarrow\hat{Q}(0)\geq0$ ; then

(5.3) \begin{equation}\bigl(\hat{Q}^{(n)}_p,\hat{Z}^{(n)}_p\bigr)\Rightarrow\bigl(\hat{Q}_p,\hat{Z}_p\bigr) \quad\text{in $ ({\mathbb D}^{2d}, J_1)$} \quad\text{as $ n \to \infty$,}\end{equation}

with the limit

\begin{align*}\hat{Q}_p=\phi (\mathcal{M}_{\gamma}(\hat{Q}(0) + \hat{W}-\hat{\theta} \mathfrak{e} ))\quad\text{and}\quad \hat{Z}_p= \mu^{-1} \hat{Q}_p,\end{align*}

where $\hat{W}$ is as given in (5.2).

We remark that the proofs for both models follow similar arguments by adapting those with a split Hawkes process, which we omit for brevity. However, we note that the split Hawkes processes as arrivals to the parallel server queues have a particular structure due to the splitting mechanism, while the multivariate Hawkes process has a matrix kernel $H=(H_{ij})$ (independent of n). It is important to highlight that the limit processes in both parallel server queueing models with a split Hawkes arrival process and a general multivariate Hawkes process are essentially of the same nature (with orthonormal reflections), except that the driving Brownian motions have different covariance functions (comparing (3.11) and (5.2)).

6. Proofs

6.1. Proof of Theorem 3.1

The proof follows from some modification of the arguments for the heavy-traffic convergence in single-server queues; see e.g. [Reference Chen and Yao7, Chapter 6]. We highlight the differences here for the parallel single-server queueing model.

Recall the diffusion-scaled processes $\hat{Q}^{(n)}_k$ and $\hat{Z}^{(n)}_k$ for the nth system from (3.1):

(6.1) \begin{equation}\begin{aligned}\hat{Z}^{(n)}_{k}(t)&= \dfrac{1}{\sqrt{n}} Z^{(n)}_{k}(nt)=\sqrt{n}\bigl(\bar{V}^{(n)}_{k}\bigl(\bar{Q}^{(n)}_{k}(0)+\bar{N}^{(n)}_{k}(t)\bigr)-\bar{B}^{(n)}_{k}(t)\bigr), \\[5pt]\hat{Q}^{(n)}_{k}(t)&= \dfrac{1}{\sqrt{n}} Q^{(n)}_{k}(nt)=\sqrt{n}\bigl(\bar{Q}^{(n)}_{k}(0)+\bar{N}^{(n)}_{k}(t)-\bar{U}^{(n)}_{k}\bigl(\bar{B}^{(n)}_{k}(t)\bigr)\bigr),\end{aligned}\end{equation}

where

\begin{align*}\bar{Q}^{(n)}_{k}(0)=\frac{1}{n}Q^{(n)}_{k}(0)\end{align*}

and $\bar{N}^{(n)}_{k}$ is the process in (2.3). We will apply the FCLT established for the split Hawkes processes and renewal process, Proposition 2.2 and Remark 3.1 respectively, and the continuous mapping approach to the multidimensional reflection mapping; see e.g. [Reference Chen and Yao7, Theorems 6.1 and 7.2] and [Reference Whitt32, Chapter 14.2]. The following lemma is a modification of the so-called random change of time result from [Reference Billingsley4, page 151], which is also called the continuity of composition in [Reference Whitt32, Theorem 13.2.2]. In this version, each sub-counting process has a different time scaling, and all the limit processes are assumed to have continuous sample paths. The conditions are slightly different from those in [Reference Whitt32, Chapter 14.2]. Recall that ${\mathbb C}_{\uparrow}^{d}={\mathbb C}^{d}\cap{\mathbb D}_{\uparrow}^{d}$ , where ${\mathbb C}^{d}$ denotes the space of ${\mathbb R}^{d}$ -valued continuous functions.

Lemma 6.1. Let $\bigl(x^{(n)},\chi^{(n)}\bigr)=\bigl(x^{(n)}_{k},\chi^{(n)}_{k}\bigr)_{k}\in{\mathbb D}^{d}\times{\mathbb D}_{\uparrow}^{d}$ and $(x,\chi)=(x_{k},\chi_{k})_{k}\in{\mathbb C}^{d}\times{\mathbb C}_{\uparrow}^{d}$ . If

\begin{align*}\bigl(x^{(n)},\chi^{(n)}\bigr)\rightarrow (x,\chi)\quad\textit{in}\;\text{$ ({\mathbb D}^{2d}, J_{1})$} \quad\textit{as}\;\text{$ n \to \infty$,}\end{align*}

then

\begin{align*}x^{(n)}\circ\chi^{(n)}=\bigl(x^{(n)}_{k}\circ\chi^{(n)}_{k}\bigr)_{k}\rightarrow (x_{k}\circ\chi_{k})_{k}=x\circ\chi\quad\textit{u.o.c.} \quad\textit{as}\;\text{$ n \to \infty$.}\end{align*}

Proof of Theorem 3.1. We have from (6.1), for every $ k=1,2,\ldots,d$ ,

\begin{align*}\hat{Q}^{(n)}_{k}(t)&= \sqrt{n}\bar{Q}^{(n)}_{k}(0) +\mu^{(n)}_{k}\sqrt{n}\biggl(\dfrac{\lambda_{0}^{(n)}p^{(n)}_{k} }{(1-\|H\|_{1}) \mu^{(n)}_{k}}-1\biggr)t \\[5pt]&\quad + \sqrt{n}\biggl(\bar{N}^{(n)}_{k}(t)-\dfrac{\lambda_{0}^{(n)}p^{(n)}_{k}t}{1-\|H\|_{1}}\biggr)-\sqrt{n}\bigl(\bar{U}^{(n)}_{k}(s)-\mu^{(n)}_{k}s\bigr)\bigr|_{s=\bar{B}^{(n)}_{k}(t)}\\[5pt]&\quad + \mu^{(n)}_{k}\sqrt{n}\bigl(t-\bar{B}^{(n)}_{k}(t)\bigr)\\[5pt]&= \hat{Q}^{(n)}_{k}(0)-\hat{\theta}^{(n)}_{k} t+ \check{N}^{(n)}_k(t)-\hat{U}^{(n)}_k\circ\bar{B}^{(n)}_k(t)+\mu^{(n)}_{k}\hat{I}^{(n)}_{k}(t),\end{align*}

where

\begin{align*}\hat{\theta}^{(n)}_k=\mu^{(n)}_{k}\cdot \sqrt{n}\bigl(1-\rho^{(n)}_k\bigr)\quad\text{and}\quad\hat{I}^{(n)}_{k}(t) = \sqrt{n}\bigl(t-\bar{B}^{(n)}_{k}(t)\bigr),\end{align*}

$\rho^{(n)}_k$ is defined in (3.3), and $\check{N}^{(n)}_k(t)$ is defined in (2.7). Observe that what differs from a single-server queue is that the $\check{N}^{(n)}_k(t)$ are correlated Hawkes processes, which introduces dependence among the different queueing processes.

By the definition of $\hat{Z}^{(n)}$ in (6.1) and $\bigl(\hat{V}^{(n)},\hat{U}^{(n)}\bigr)$ in (3.4), we also obtain

(6.2) \begin{equation}\hat{Z}^{(n)}- m^{(n)}\hat{Q}^{(n)}=\hat{V}^{(n)}\circ\bigl(\bar{Q}^{(n)}(0)+\bar{N}^{(n)}\bigr)+ m^{(n)} \hat{U}^{(n)}\circ\bar{B}^{(n)}.\end{equation}

We now consider the convergence of various components in these representations.

Under Assumption 3.1 and the conditions in Proposition 2.2 for Hawkes processes, from (2.9) for $\check{N}^{(n)}$ , (3.9) for $(\hat{V}^{(n)},\hat{U}^{(n)})$ and their independences, we have the joint weak convergence

\begin{align*}\bigl(\check{N}^{(n)},\hat{V}^{(n)},\hat{U}^{(n)}\bigr)\Rightarrow(\hat{N},\hat{V},\hat{U})\quad\text{in $ ({\mathbb D}^{3d},J_{1})$} \quad\text{as $ n \to \infty$,}\end{align*}

where $\hat{N}$ is the Brownian motion in Proposition 2.1, $(\hat{V},\hat{U})$ is the Brownian motion in (3.9) and independent of $\hat{N}$ .

Moreover, it is easy to see that

(6.3) \begin{equation}\bigl(\bar{Q}^{(n)}(0), \bar{N}^{(n)},\bar{B}^{(n)}\bigr)\rightarrow\biggl(0, \dfrac{\lambda_{0}(p_{k})_k} {1-\|H\|_{1}}\mathfrak{e}, \mathfrak{e}\biggr)=(0, \mu\mathfrak{e}, \mathfrak{e})\quad\text{u.o.c.},\end{equation}

in probability as $n\to\infty$ , where the identity (3.8) is used in the equality, and abusing notation, $\mathfrak{e}$ is a vector of identity functions.

Thus, applying the continuous mapping theorem and the composition mapping (Lemma 6.1), we obtain the joint convergence

(6.4) \begin{equation} \bigl(\check{N}^{(n)}_k,\hat{V}^{(n)}_{k}\bigl(\bar{Q}^{(n)}_{k}(0)+\bar{N}^{(n)}_{k}\bigr),\hat{U}^{(n)}_{k}\bigl(\bar{B}^{(n)}_{k}\bigr)\bigr)_{k}\Rightarrow\bigl(\hat{N}_k,\hat{V}_{k}(\mu_{k}\cdot), \hat{U}_{k}\bigr)_{k}\end{equation}

in $({\mathbb D}^{3d},J_{1})$ as $n\to\infty$ .

Now, applying the continuous mapping theorem to the multi-dimensional reflection mapping in Definition 3.1 and using the convergence results in (6.3) and (6.4), together with the fact $ \hat{U}_{k}(t)=-\mu_{k}\hat{V}_{k}(\mu_{k}t)$ from (3.9), we obtain the joint convergence of $\bigl(\hat{Q}^{(n)},\hat{Z}^{(n)}\bigr)$ in (3.10).

Finally, since $\hat{N}$ and $\hat{U}$ are independent, it is easy to check that $\hat{W}=\hat{N}-\hat{U}$ is a Brownian motion with the covariance function given in (3.11). This completes the proof.

6.2. Proof of Theorem 4.1

We adapt the proof idea in [Reference Ward and Glynn30] and highlight the main differences in the proof. For the system indexed by n, recall that the diffusion-scaled processes are defined by

\begin{equation*}\hat{Q}_{k,p}^{(n)}(t)=\dfrac{1}{\sqrt{n}}Q_{k,p}^{(n)}(nt),\quad\hat{Z}_{k,p}^{(n)}(t)=\dfrac{1}{\sqrt{n}}Z_{k,p}^{(n)}(nt), \quad\hat{\tilde{Z}}_{k,p}^{(n)}(t)=\dfrac{1}{\sqrt{n}}\tilde{Z}_{k,p}^{(n)}(nt),\end{equation*}

and

\begin{align*}\hat{I}^{(n)}_{k,p}(t) = \sqrt{n}\bigl(t-\bar{B}^{(n)}_{k,p}(t)\bigr).\end{align*}

Theorem 4.1 is proved following the procedure from [Reference Ward and Glynn30], where we start from the analysis of the offered workload process in (4.1). Specifically, the proof proceeds in the following steps.

1. Step 1. The weak convergence of the diffusion-scaled process $\bigl(\hat{\tilde{Z}}^{(n)}_{p},\hat{I}^{(n)}_{p}\bigr)$ in Theorem 6.1, by relating to $\bigl(\hat{Z}^{(n)},\hat{I}^{(n)}\bigr)$ in the corresponding model without abandonment.

2. Step 2. The following asymptotic equivalence properties: for every $T>0$ and k, as $n\to \infty$ , we have

   \begin{align*}\sup_{t\le T}\bigl|\hat{Z}_{k,p}^{(n)}(t)-\hat{\tilde{Z}}_{k,p}^{(n)}(t)\bigr|\Rightarrow 0\quad\text{and}\quad\sup_{t\le T}\bigl|\hat{\tilde{Z}}_{k,p}^{(n)}(t)- m^{(n)}_{k}\hat{Q}_{k,p}^{(n)}(t)\bigr|\Rightarrow 0.\end{align*}
   They follow the same procedure as [Reference Ward and Glynn30], so their proofs are given in the Appendix for completeness.

3. Step 3. Completing the proof: given the convergence results in the two steps, the joint convergence in Theorem 4.1 follows immediately.

Therefore we focus only on the proof of the following theorem.

Theorem 6.1. Under the assumptions of Theorem 4.1,

(6.5) \begin{equation}\bigl(\hat{\tilde{Z}}_{p}^{(n)},\hat{I}_{p}^{(n)}\bigr)\Rightarrow \bigl(\hat{Z}_{p},\hat{I}_{p}\bigr) \quad\textit{in}\;\text{$ ({\mathbb D}^{2d}, J_1)$} \quad\textit{as}\;\text{$ n \to \infty$,}\end{equation}

with the limit

\begin{align*} \bigl(\hat{Z}_{p},\hat{I}_{p}\bigr) \;:\!=\;(\phi_{0},\psi_{0}) \mathcal{M}_{\gamma}(\hat{Z}+\hat{I})= m \, (\phi_{0},\psi_{0}) \mathcal{M}_{\gamma}(\hat{Q}(0)+\hat{W} -\hat{\theta} \mathfrak{e} )\end{align*}

where $\hat{W}$ and $\hat{\theta}$ are as given in Theorem 3.1.

The proof of Theorem 6.1 uses the following two lemmas. In the proof we need Theorem 3.1, and it helps to rewrite $Z^{(n)}_{k}(t)$ in (3.1) in the corresponding queueing model without abandonment as

\begin{align*}Z^{(n)}_{k}(t)=Z^{(n)}_{k}(0)+\sum_{j\ge1} \eta^{(n)}_{j,k}\textbf{1}\bigl(\tau^{(n)}_{j,k}\le t\bigr)-B^{(n)}_{k}(t),\end{align*}

with

\begin{align*}Z^{(n)}_{k}(0)=\sum_{j\ge1} \eta_{-j,k}^{(n)}\textbf{1}\bigl(\,j\le Q^{(n)}_{k}(0)\bigr).\end{align*}

The diffusion-scaled process $\hat{Z}^{(n)}$ is defined in the same way. The corresponding convergence results for $\bigl(\hat{Q}^{(n)},\hat{Z}^{(n)}\bigr)$ in Theorem 3.1 will be used.

Define

\begin{align*}\hat{M}_{k,p,1}^{(n)}(t)=\dfrac{1}{\sqrt{n}}M_{k,p,1}^{(n)}([nt])\quad\text{and}\quad \hat{M}_{k,p,2}^{(n)}(t)=\dfrac{1}{\sqrt{n}}M_{k,p,2}^{(n)}([nt])\end{align*}

with

\begin{align*}&M_{k,p,1}^{(n)}(m)=\sum_{j=1}^{m}\bigl( \eta^{(n)}_{j,k}- m^{(n)}_{k}\bigr)\textbf{1}\bigl(\vartheta^{(n)}_{j,k,p}\le u\bigr)\big|_{u=\tilde{Z}_{k,p}^{(n)}(\tau^{(n)}_{j,k}-)},\\[5pt]&M_{k,p,2}^{(n)}(m)=\sum_{j=1}^{m}\bigl(\textbf{1}\bigl(\vartheta^{(n)}_{j,k,p}\le u\bigr)-G^{(n)}_{k}(u)\bigr)\big|_{u=\tilde{Z}_{k,p}^{(n)}(\tau^{(n)}_{j,k}-)}.\end{align*}

Lemma 6.2. For every $T>0$ and every k,

(6.6) \begin{equation}\sup_{t\le T}\bigl|\hat{M}_{k,p,1}^{(n)}(t)\bigr|\Rightarrow 0\quad\textit{and}\quad\sup_{t\le T}\bigl|\hat{M}_{k,p,2}^{(n)}(t)\bigr|\Rightarrow 0 \quad\textit{as}\;\text{$ n \to \infty$.}\end{equation}

Proof. It is easy to see that $\bigl\{M_{k,p,1}^{(n)}(m)\bigr\}_{m\ge1}$ is an $\bigl\{\mathscr{F}_{k,p,m}^{(n)}\bigr\}_{m\ge1}$ -adapted square-integrable martingale, where for $m\ge1$

\begin{align*}\mathscr{F}_{k,p,m}^{(n)}=\tilde{Z}_{k,p}^{(n)}(0)\vee\sigma\bigl\{\eta^{(n)}_{j,k},\vartheta^{(n)}_{j,k,p},\tau^{(n)}_{j,k}\bigr\}_{1\le j\le m}\vee\sigma\bigl\{\tau_{m+1,k}^{(n)}\bigr\}.\end{align*}

Applying Doob’s maximal inequality [Reference Shiryaev27, Theorem VII.3.3], we have

\begin{align*} \mathbb{E}\biggl[\sup_{t\le T}\bigl(\hat{M}_{k,p,1}^{(n)}(t)\bigr)^{2}\biggr]& \le 4\mathbb{E}\bigl[\bigl(\hat{M}_{k,p,1}^{(n)}(T)\bigr)^{2}\bigr]\\[5pt]& = \dfrac{4}{n}\textrm{var}\bigl( \eta^{(n)}_{k}\bigr)\sum_{j=1}^{[nT]}\mathbb{E}\bigl[G^{(n)}_{k}\bigl(\tilde{Z}_{k,p}^{(n)}\bigl(\tau_{j,,k}^{(n)}-\bigr)\bigr)\bigr]\\[5pt]& \le 4\textrm{var}\bigl( \eta^{(n)}_{k}\bigr) \dfrac{[nT]}{n}\biggl(\dfrac{1}{\sqrt{n}}\hat{G}^{(n)}_{k}(K_{0})+ \mathbb{P}\biggl(\sup_{\bar{N}^{(n)}_{k}(t)\le T}\hat{\tilde{Z}}_{k,p}^{(n)}(t)\ge K_{0}\biggr)\biggr)\end{align*}

for every T and $K_{0}$ . Together with Proposition 2.2 for $\bar{N}^{(n)}_{k}$ , the fact that

(6.7) \begin{equation}0\le \hat{\tilde{Z}}_{k,p}^{(n)}(t)\le \hat{Z}_{k,p}^{(n)}(t)\le \hat{Z}^{(n)}_{k}(t)\quad\text{for all $t\ge0$}\end{equation}

and Theorem 3.1 is established for $\hat{Z}^{(n)}_{k}$ , we can establish the weak convergence result: $\sup_{t\le T}\bigl|\hat{M}_{k,p,1}^{(n)}(t)\bigr|\Rightarrow 0$ as $n \to \infty$ . A similar procedure can be used to prove the convergence for $\hat{M}_{k,p,2}^{(n)}$ .

Recall the following martingale representations associated with Hawkes processes [Reference Bacry, Delattre, Hoffmann and Muzy3, Reference Li and Pang23]:

(6.8) \begin{equation}X^{(n)}(t)\;:\!=\; N^{(n)}(t)-\int_{0}^{t}\lambda^{(n)}(s)\,{\textrm{d}} s \quad\text{and}\quad X^{(n)}_{k}(t)\;:\!=\; N^{(n)}_{k}(t)-\int_{0}^{t}\lambda^{(n)}_{k}(s)\,{\textrm{d}} s\end{equation}

are martingales adapted to the filtrations generated by $N^{(n)}$ and $N^{(n)}_{k}$ , respectively. Denote

\begin{align*}\bar{X}^{(n)}_{k}(t)\;:\!=\; \dfrac{1}{n}X^{n}_{k}(nt)=\bar{N}^{(n)}_{k}(t)-\int_{0}^{t}\bar{\lambda}^{(n)}_{k}(s)\,{\textrm{d}} s\quad\text{and}\quad\bar{X}^{(n)}(t)=\dfrac{1}{n}X^{(n)}(nt).\end{align*}

Lemma 6.3. For every $T>0$ and every k,

\begin{align*}\sup_{t\le T}\biggl|\int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigl({\textrm{d}} \bar{N}^{(n)}_{k}(s)-\mathbb{E}\bigl[\bar{\lambda}^{(n)}_{k}(s)\bigr]\,{\textrm{d}} s\bigr)\biggr|\Rightarrow 0 \quad\textit{as}\;\text{$ n \to \infty$.}\end{align*}

Proof. Using the martingales in (6.8), the integral in the lemma can be split into two martingale integrals and a third component with bounded variation, and we show that each term converges to 0 uniformly on [0, T]. Specifically, we have

(6.9) \begin{align}& \int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigl({\textrm{d}} \bar{N}^{(n)}_{k}(s)-\mathbb{E}\bigl[\bar{\lambda}^{(n)}_{k}(s)\bigr]\,{\textrm{d}} s\bigr)\notag \\[-1pt]&\quad = \int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigl({\textrm{d}} \bar{X}^{(n)}_{k}(s)+p^{(n)}_{k}\bigl(\bar{\lambda}^{(n)}(s)-\mathbb{E}\bigl[\bar{\lambda}^{(n)}(s)\bigr]\bigr)\,{\textrm{d}} s\bigr)\notag \\[-1pt]&\quad = \int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr) \,{\textrm{d}} \bar{X}^{(n)}_{k}(s)+ p^{(n)}_{k}\int_{0}^{t}\,{\textrm{d}} \bar{X}^{(n)}(r)\biggl(n\int_{r}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\varphi(n(s-r))\,{\textrm{d}} s\biggr)\notag \\[-1pt]&\quad = \int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigl({\textrm{d}} \bar{X}^{(n)}_{k}(s)+p^{(n)}_{k}\|\varphi\|_{1}\,{\textrm{d}} \bar{X}^{(n)}(s)\bigr)\notag \\[-1pt]&\quad\quad +p^{(n)}_{k}\int_{0}^{t}\,{\textrm{d}} \bar{X}^{(n)}(s)\biggl(n\int_{s}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(r-)\bigr)\varphi(n(r-s))\,{\textrm{d}} r-\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\|\varphi\|_{1}\biggr),\end{align}

where we make use of (2.2) and the following identity from [Reference Bacry, Delattre, Hoffmann and Muzy3]:

\begin{align*}\lambda^{(n)}(t)-\mathbb{E}\bigl[\lambda^{(n)}(t)\bigr]=\int_{0}^{t}\varphi(t-s) \,{\textrm{d}} X^{(n)}(s),\end{align*}

and recall that $\varphi=\sum_{j\ge1} H^{*j}$ is the renewal function of H.

Since $X^{(n)}_{k}, X^{(n)}$ are martingales with quadratic variation $N^{(n)}_{k}$ and $N^{(n)}$ , respectively, we have from Doob’s maximal inequality that

\begin{align*}& \mathbb{E}\biggl[\sup_{t\in[0,T]}\biggl(\int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\,{\textrm{d}} \bar{X}^{(n)}_{k}(s)\biggr)^{2}\biggr]\\[-1pt]&\quad \le 4 \mathbb{E}\biggl[\biggl(\int_{0}^{T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\,{\textrm{d}} \bar{X}^{(n)}_{k}(s)\biggr)^{2}\biggr]= \dfrac{4}{n} \mathbb{E}\biggl[\int_{0}^{T}\bigl(\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigr)^{2}\bar{\lambda}^{(n)}_{k}(s)\,{\textrm{d}} s\biggr]\\[-1pt]&\quad \le \dfrac{4}{n}\bigl(\hat{G}^{(n)}_{k}(K_{0})\bigr)^{2}\mathbb{E}\biggl[\int_{0}^{T}\bar{\lambda}^{(n)}_{k}(s)\,{\textrm{d}} s\biggr]+ 4\mathbb{E}\biggl[\int_{0}^{T}\bar{\lambda}^{(n)}_{k}(s)\,{\textrm{d}} s;\ \sup_{s\le T}\hat{\tilde{Z}}_{k,p}^{(n)}(s)>K_{0} \biggr]\end{align*}

for every $K_{0}>0$ , where $\hat{G}^{(n)}_{k}(z)\le\sqrt{n}$ by definition. Given (6.7), we have

\begin{align*}\sup_{t\le T}\biggl|\int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\,{\textrm{d}} \bar{X}^{(n)}_{k}(s)\biggr|\Rightarrow 0 \quad\text{as $ n \to \infty$.}\end{align*}

One can prove the uniform convergence of the second term in (6.9) similarly.

For the last integral in (6.9), notice that $n\varphi(n(s-r))\,{\textrm{d}} s$ degenerates to a Dirac measure at r as $n\to\infty$ and $\bar{X}^{(n)}$ has bounded variation on [0, T]. For arbitrary $\delta_{n}>0$ , if $t-s>\delta_{n}$ , we have

\begin{align*}& n\int_{s}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(r-)\bigr)\varphi(n(r-s))\,{\textrm{d}} r-\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\|\varphi\|_{1} \\[-1pt]&\quad \le n\int_{s}^{s+\delta_{n}}\bigl|\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(r-)\bigr)-\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigr|\varphi(n(r-s))\,{\textrm{d}} r \\[-1pt]& \qquad+ 2\sup_{u\le T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr) \int_{s+\delta_{n}}^{\infty}n\varphi(n(r-s))\,{\textrm{d}} r \\[-1pt]&\quad \le \sup_{\substack{0<v<u\le T\\[5pt] u-v\le\delta_{n}}}\bigl|\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)-\hat{G}_{k,p}^{(n)}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(v)\bigr)\bigr|\cdot \|\varphi\|_{1}+ 2\sup_{u\le T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)\int_{n\delta_{n}}^{\infty}\varphi(u)\,{\textrm{d}} u.\end{align*}

On the other hand, if $t-s<\delta_{n}$ , we have

\begin{align*}n\int_{s}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(r-)\bigr)\varphi(n(r-s))\,{\textrm{d}} r-\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\|\varphi\|_{1} \le \sup_{u\le T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)\cdot \|\varphi\|_{1}.\end{align*}

Plugging these into the last integral in (6.9), we obtain the following bound:

\begin{align*}& \sup_{\substack{0<v<u\le T\\[5pt] u-v\le\delta_{n}}}\bigl|\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)-\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(v)\bigr)\bigr|\cdot \|\varphi\|_{1}\cdot\int_{0}^{T}\bigl(\bar{N}^{(n)}({\textrm{d}} s)+\bar{\lambda}^{(n)}({\textrm{d}} s)\bigr)\\[5pt]&\quad + 2 \int_{n\delta_{n}}^{\infty}\varphi(u)\cdot \sup_{u\le T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)\cdot \int_{0}^{T}\bigl(\bar{N}^{(n)}({\textrm{d}} s)+\bar{\lambda}^{(n)}({\textrm{d}} s)\bigr)\\[5pt]&\quad + \sup_{u\le T}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(u)\bigr)\cdot\|\varphi\|_{1}\cdot \sup_{\substack{0<v<u\le T\\[5pt] u-v\le\delta_{n}}}\biggl(\bigl|\bar{N}^{(n)}(u)-\bar{N}^{(n)}(v)\bigr|+\int_{v}^{u}\bar{\lambda}^{(n)}(s)\,{\textrm{d}} s\biggr).\end{align*}

We next need the following result. For every $T>0$ and $\delta_{n}\to0$ with $\sqrt{n}\delta_{n}\to0$ , we have for every k

(6.10) \begin{equation} \sup_{\substack{0<s<t<T\\[5pt] t-s\le \delta_{n}}}\bigl|\hat{\tilde{Z}}_{k,p}^{(n)}(t)-\hat{\tilde{Z}}_{k,p}^{(n)}(s)\bigr|\Rightarrow 0 \quad\text{as $ n \to \infty$.}\end{equation}

It follows from their definitions that

\begin{align*}\bigl|\hat{\tilde{Z}}_{k,p}^{(n)}(t)-\hat{\tilde{Z}}_{k,p}^{(n)}(s)\bigr|\le \bigl|\hat{Z}^{(n)}_{k}(t)-\hat{Z}^{(n)}_{k}(s)\bigr|+2\sqrt{n}\delta_{n}\end{align*}

for every $t>s>0$ with $t-s\le \delta_{n}$ . Therefore the convergence of $\hat{Z}^{(n)}_{k}$ in Theorem 3.1 and the fact that $\hat{Z}^{(n)}_{k}\in {\mathbb C}$ can be applied to show the convergence property in (6.10).

Now we continue with the last integral in (6.9), taking $\delta_{n}=n^{-2/3}$ and then $\sqrt{n}\delta_{n}\to0$ , $n\delta_{n}\to\infty$ , by (6.10); we prove the uniform convergence of the last piece and finish the proof.

Now we are ready to prove Theorem 6.1 by making use of Lemmas 6.2 and 6.3.

Proof of Theorem 6.1. We claim that as $n\to\infty$ ,

(6.11) \begin{equation}\bigl(\hat{Z}^{(n)}_{k}(t)-\hat{I}^{(n)}_{k}(t)\bigr)-\biggl(\hat{\tilde{Z}}_{k,p}^{(n)}(t)+\gamma_{k}\int_{0}^{t}\hat{\tilde{Z}}_{k,p}^{(n)}(s)\,{\textrm{d}} s-\hat{I}_{k,p}^{(n)}(t)\biggr)\Rightarrow 0\quad\text{u.o.c.}\end{equation}

for every k. The claim can also be rephrased as

(6.12) \begin{equation}\hat{\tilde{Z}}_{p}^{(n)}(t)+ \text{diag} (\gamma) \int_{0}^{t}\hat{\tilde{Z}}_{p}^{(n)}(s)\,{\textrm{d}} s=\bigl(\hat{Z}^{(n)}(t)-\hat{I}^{(n)}(t)+\hat{\varepsilon}_{p}^{(n)}(t)\bigr)+\hat{I}_{p}^{(n)}(t),\end{equation}

where $\hat{I}_{p}^{(n)}$ is the regulator of $\hat{\tilde{Z}}_{p}^{(n)}$ satisfying the conditions of Definition 4.1 with the reflection matrix $\texttt{Q}\equiv 0$ , and $\hat{\varepsilon}_{p}^{(n)}$ is the error term which converges weakly to 0 u.o.c. In other words,

\begin{align*}\bigl(\hat{\tilde{Z}}_{p}^{(n)},\hat{I}_{p}^{(n)}\bigr)=(\phi_{0},\psi_{0})\bigl(\mathcal{M}_{\gamma}\bigl(\hat{Z}^{(n)}-\hat{I}^{(n)}+\hat{\varepsilon}_{p}^{(n)}\bigr)\bigr).\end{align*}

We have shown in Section 6.1 that $\hat{Z}^{(n)}-\hat{I}^{(n)}\Rightarrow m(\hat{Q}(0)-\hat{\theta}\mathfrak{e}+W)$ . By Theorem 3.1 we obtain joint convergence in (6.5).

To obtain (6.11), denote

\begin{align*}\hat{M}_{k,p,0}^{(n)} \;:\!=\; \dfrac{1}{\sqrt{n}}\sum_{j\ge1} \eta_{-j,k}^{(n)}\textbf{1}\bigl(\,j\le Q^{(n)}_{k}(0), \vartheta_{-j,k,p}^{(n)}\le \Theta_{k,p}^{(n)}(\,j-1)\bigr).\end{align*}

A simple calculation gives

(6.13) \begin{align}\text{LHS of (6.11)}&=\hat{M}_{k,p,0}^{(n)}+ \dfrac{1}{\sqrt{n}}\sum_{j=1}^{N^{(n)}_{k}(nt)} \eta^{(n)}_{j,k}\textbf{1}\bigl(\vartheta^{(n)}_{j,k,p}\le\tilde{Z}_{k,p}^{(n)}\bigl(\tau^{(n)}_{j,k}-\bigr)\bigr)-\gamma_{k}\int_{0}^{t}\hat{\tilde{Z}}_{k,p}^{(n)}(s)\,{\textrm{d}} s\notag \\[5pt]&= \hat{M}_{k,p,0}^{(n)}+\hat{M}_{k,p,1}^{(n)}\bigl(\bar{N}^{(n)}_{k}(t)\bigr)+\hat{M}_{k,p,2}^{(n)}\bigl(\bar{N}^{(n)}_{k}(t)\bigr)\notag \\[5pt]&\quad + m^{(n)}_{k}\int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s-)\bigr)\bigl({\textrm{d}} \bar{N}^{(n)}_{k}(s)-\mathbb{E}\bigl[\bar{\lambda}^{(n)}_{k}(s)\bigr]\,{\textrm{d}} s\bigr)\notag \\[5pt]&\quad - \lambda_{0}^{(n)} m^{(n)}_{k} p^{(n)}_{k}\int_{0}^{t}\hat{G}^{(n)}_{k}\bigl(\hat{\tilde{Z}}_{k,p}^{(n)}(s)\bigr)\biggl(\int_{ns}^{\infty}\varphi(u)\,{\textrm{d}} u\biggr)\,{\textrm{d}} s\notag \\[5pt]&\quad +\int_{0}^{t}\bigl(\rho^{(n)}_{k}\hat{G}^{(n)}_{k}(u)-\gamma_{k}u\bigr)\bigr|_{u=\hat{\tilde{Z}}_{k,p}^{(n)}(s)}\,{\textrm{d}} s,\end{align}

where $\bar{\lambda}^{(n)}_{k}(t)=\lambda^{(n)}_{k}(nt)$ and we use the expression of $\mathbb{E}\bigl[\bar{\lambda}^{(n)}_{k}(t)\bigr]$ in (2.4). It is thus sufficient to show that all terms on the right-hand side of (6.13) converge weakly to 0 u.o.c.

For the term $\hat{M}_{k,p,0}^{(n)}$ , since $\Theta_{k,p}^{(n)}(j)\le Z^{(n)}_{k}(0)$ for all $j\le Q^{(n)}_{k}(0)$ , we have

(6.14) \begin{equation} \mathbb{E}\bigl[\hat{M}_{k,p,0}^{(n)};\;\hat{Z}^{(n)}_{k}(0)\le K_{0}\bigm|Q^{(n)}_{k}(0)\bigr]\leq\bar{Q}^{(n)}_{k}(0) m^{(n)}_{k} \hat{G}^{(n)}_{k}(K_{0})\end{equation}

for every $K_{0}>0$ . The Markov inequality can be applied to show that $ \hat{M}_{k,p,0}^{(n)}\Rightarrow 0$ as $n\to \infty$ .

Now, given Theorem 3.1 for $\hat{Z}^{(n)}_{k}$ and (6.7), the desired convergences of the last two terms in (6.13) can be checked directly. Further, applying Lemma 6.2, Lemma 6.3, and Proposition 2.1, we can prove the remaining terms and finish the proof.