Heavy-traffic limits for parallel single-server queues with randomly split Hawkes arrival processes

Abstract

We consider parallel single-server queues in heavy traffic with randomly split Hawkes arrival processes. The service times are assumed to be independent and identically distributed (i.i.d.) in each queue and are independent in different queues. In the critically loaded regime at each queue, it is shown that the diffusion-scaled queueing and workload processes converge to a multidimensional reflected Brownian motion in the non-negative orthant with orthonormal reflections. For the model with abandonment, we also show that the corresponding limit is a multidimensional reflected Ornstein–Uhlenbeck diffusion in the non-negative orthant.

1. Introduction

Single-server queues are fundamental models in applied probability and queueing theory, and heavy-traffic limits have been well studied; see the recent surveys in [7], [31], and [32], and the recent work on the model with abandonment in [28] and [30]. In these studies, the arrival processes are usually assumed to be Poisson or renewal processes. Recently, to account for the excessive burstiness, clustering effects, and path-dependence, other point processes such as Hawkes and Pólya processes have been used to model arrivals in queueing models. They can be used to model, for example, internet/social media traffic flows, patient flows during a pandemic, neuron interaction processes, and high-frequency transaction processes in limit order books. Queues with Hawkes input are studied in [8], [11], [16], [22], and [26], and single-server queues with Pólya arrival processes are considered in [14] and [15]. On the other hand, many systems have parallel servers, such as data centers and internet processors, and various randomized routeing schemes have been developed in the literature (see e.g. the recent survey [12]). A simple scheme is to send incoming jobs/customers randomly to any of the servers upon arrival.

In this paper we consider parallel single-server queues in heavy traffic, where the arrival process for each queue is randomly split from one Hawkes process. That is, upon each arrival, a customer or job is randomly assigned to one of the queues. The service times in the different queues are mutually independent and may have different distributions, and within a queue the service times are assumed to be independent and identically distributed (i.i.d.). The service discipline in each queue is first-come first-served (FCFS). We assume that each queue is critically loaded in heavy traffic, that is, the traffic intensity gets close to one. We also consider the queueing model with abandonment, where customers joining each queue may abandon the system while waiting in the queue and before receiving service. We focus on the joint queue length and workload processes, and prove the functional central limit theorem (FCLT) for the associated diffusion-scaled processes.

Randomly split Hawkes processes have recently been studied by the present authors [23]. Unlike Poisson processes, the split Hawkes processes become a multivariate Hawkes process with a particular dependence structure (see the covariance function of scaling limits in Proposition 2.1). As a consequence, the vector of queueing (workload) processes of the parallel server queues will be correlated. In contrast, in the case of randomly split Poisson processes, because of the independence property of Poisson thinning/sampling, the queues at the parallel servers are mutually independent and each queue can be analyzed as an M/G/1 queue. However, our model becomes more challenging to analyze directly. This is in a similar flavor to the open problem recently posed by Mandjes [24] about multivariate M/G/1 queues with coupled input and parallel service. Our model presents another example of multivariate single-server queues with correlated input, and our results provide heavy-traffic diffusion approximations for such a system.

For the queueing models without abandonment, the queueing limit process is a multidimensional reflected Brownian motion (RBM) in the non-negative orthant with orthonormal reflections. In many queueing networks, RBMs arise as the heavy-traffic limits of the queueing and workload processes; see e.g. the recent surveys [10] and [34]. In the very original article by Harrison [17], as a diffusion approximation for a pair of single-server queues in series, a two-dimensional RBM in the non-negative quadrant was introduced, which has a normal reflection at one axis and a tangential reflection at the other axis. On the other hand, our model with two queues provides an example of RBMs in the non-negative quadrant with normal reflections at both axes. A multidimensional RBM limit with orthonormal reflections was derived for a network of parallel single-server queues with Markov-modulated service speeds in [13]. This work contributes to the literature of queueing network scaling limits by providing a concrete example of RBMs in a non-negative orthant with orthonormal reflections.

Although the correlation structure of the split Hawkes processes can be explicitly characterized (as can the Brownian motion without reflection), it is challenging to compute the covariance functions for the limiting RBM in the parallel server queueing model. Moreover, one can check that the conditions of Propositions 8 and 9 in [9] (see also [33]) are not satisfied so that the limit process does not possess a product-form stationary distribution. Fortunately, the numerical approach developed in [9] to compute the stationary distribution can be used for the limiting RBM in our model. We implement that algorithm in a two-queue example to illustrate how the parameters in the Hawkes process and the splitting probabilities as well as the heavy-traffic scalings affect the correlation between the steady state of the two queues in the limit (see Section 3.1).

For the queueing models with abandonment, the queueing limit process is a multidimensional reflected Ornstein–Uhlenbeck (OU) diffusion in the non-negative orthant with orthonormal reflections. This extends the one-dimensional results in [30]. This result is of interest in two respects. First, there is a very limited result in the queueing network literature that gives a multidimensional reflected OU limit. Huang and Zhang [20] study open Jackson networks with reneging, and obtain a multidimensional reflected OU diffusion limit with a reflection being characterized via the routeing probability matrix (extending the RBM approximations for open Jackson networks without abandonment in [25]). The reflection is very complex in that model, but the reflections in the limits of our model are orthonormal, which is of interest in its own right. It is worth noting that in the studies of dynamic scheduling in multiclass GI/GI/1+GI queues [1, 21], although the queueing processes are multidimensional, the Brownian control problem is solved via the so-called workload process as a one-dimensional controlled OU process with reflection. Second, there have been studies of the stationary distributions and ergodic properties of reflected OU processes in one dimension; see e.g. [29] and [35]. However, such results beyond one dimension are wide open. Our limit process provides a concrete example for which an explicit characterization of the stationary distribution could potentially be obtained. We leave this as an open problem for future work. (We also refer the readers to a relevant article [2] for the recurrence properties of multidimensional reflected diffusions in convex polyhedral cones.) It would also be interesting to explore whether the numerical scheme in [9] could be further developed for the multidimensional reflected OU processes such as those arising from our model.

The proof for the model without abandonment uses the continuous mapping approach for the multidimensional reflection maps (see [31]). This is standard, but the particular scaling involved in the Hawkes process and its splitting scheme must be taken into account. For the queueing models with abandonment, we adapt the approach in [30] by comparing with the model without abandonment. It is therefore also important to first study the model without abandonment. In the meantime, the comparison approach is non-trivial, since this requires us to use the martingale representations for the Hawkes process and the randomly split Hawkes processes [3, 23], as well as some martingale inequalities and properties.

We also discuss two related models in which the arrival processes for the parallel servers come from a multivariate Hawkes process in Section 5. The limits for the models with and without abandonment are of the same type with orthonormal reflections but the driving Brownian motions have a different covariance function. Although the split Hawkes process is also a multivariate Hawkes process, the splitting scheme introduces a particular structure together with only one self-exciting function. On the other hand, a general multivariate Hawkes process has a matrix of self-exciting functions. This is yet another example of the open problems posed in [24], and also complements the parallel server queueing model with correlated services due to Markov modulation in [13]. The limits for these models will also be of interest in their own right.

1.1. Organization of the paper

In Section 2 we discuss random splitting of Hawkes processes and review the FCLT results. In Section 3 we present the parallel single-server queueing model with split Hawkes arrival processes and the FCLT on the joint queueing and workload processes. A numerical example is provided in Section 3.1 to illustrate the stationary distribution of the limiting queueing process in a model with two queues. In Section 4 we discuss the model with abandonment and state the corresponding FCLT. In Section 5 we consider the parallel server queueing model with a general multivariate Hawkes process and state the corresponding scaling limits. The proofs for these models are given in Section 6 and the Appendix.

1.2. Notation

All random variables and processes are defined in a common complete probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ . Throughout the paper, ${\mathbb N}$ denotes the set of natural numbers. ${\mathbb R}^d ({\mathbb R}^d_{+})$ denotes the space of d-dimensional real (non-negative) vectors; we write ${\mathbb R} ({\mathbb R}_{+})$ when $d=1$ . Let ${\mathbb D}^{d}={\mathbb D}^{d}({\mathbb R}_{+},{\mathbb R}^{d})$ denote the space of ${\mathbb R}^{d}$ -valued càdlàg functions on ${\mathbb R}_{+}$ . $({\mathbb D},J_{1})$ denotes the space ${\mathbb D}$ equipped with Skorokhod $J_{1}$ topology (see [4]), which is complete and separable. ${\mathbb D}^{d}$ denotes the space of ${\mathbb R}^{d}$ -valued càdlàg functions endowed with the weak Skorokhod $J_{1}$ topology [32], for which we write $({\mathbb D}^{d},J_{1})$ . For an integrable function $f\;:\; {\mathbb R}\to{\mathbb R}$ , its $L^{1}$ norm is denoted by $\|f\|_{1}$ . Notations $\to$ and $\Rightarrow$ mean convergence of real numbers and convergence in distribution, respectively. For a vector a, $\text{diag} (a)$ denotes the diagonal matrix with the elements of vector a on the main diagonal. For two vectors $a=(a_{k})_{k}$ and $b=(b_{k})_{k}$ , $a b=(a_{k}b_{k})_{k}=\text{diag} (a) b$ denotes their elementwise product. We use $\mathfrak{e}$ to denote the identity function $\mathfrak{e}(t)=t$ for $t \in {\mathbb R}_+$ . Additional notation is introduced in the paper whenever necessary.

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 [23] 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. [3]) 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 [23].

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 [3, Theorem 2] for the vector-valued Hawkes process $\bigl(N^{(n)}_{k}\bigr)_{k}$ characterized by (2.2), while the second representation follows from [32, 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 [32, 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 [9, 18]), 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 ([32, 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 [32, 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 [19]), 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 [5, 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 [18] and [33] 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 [9] (see also [33]) which means the stationary density cannot be a product function. The numerical approach developed in [9] 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 [9] to a model with two queues and compute some characteristics of the stationary distribution. Similar to [9], 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 [30] 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 [30].

Definition 4.1. Let $\Gamma$ be a $d\times d$ matrix of real entries and let $\texttt{Q}$ be a reflection matrix in Definition 3.1. Then, for each $x\in{\mathbb D}^{d}$ with $x(0)\ge0$ , let $(z,y)\in{\mathbb D}^{2d}$ be a unique pair satisfying

1. (i) $z(t)=x(t)- \int_{0}^{t}\Gamma z(s)\,{\textrm{d}} s+ (\texttt{I}-\texttt{Q})y(t) \ge 0$ for all $t\ge0$ ,

2. (ii) $y(0)=0$ , $y\in{\mathbb D}^{d}_{\uparrow}$ and $\int_{0}^{\infty}z_{k}(t) \,{\textrm{d}} y_{k}(t)=0$ for every k.

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 [30, 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 [30].

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 [3] 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.

Assumption 5.1. Suppose the following conditions hold:

1. (i) $\lambda^{(n)}_{0}\to\lambda_{0}$ as $n\to\infty$ ,

2. (ii) $ \int_{0}^{\infty}H_{ij}(t)\,{\textrm{d}} t<\infty$ and the spectral radius of the matrix $\|H\|_{1}=\bigl(\int_{0}^{\infty}H_{ij}(t)\,{\textrm{d}} t\bigr)_{ij}$ is less than 1,

3. (iii) $\lim_{t\to\infty}\sqrt{t}\int_{t}^{\infty}H_{ij}(s)\,{\textrm{d}} s=0$ for every i and j.

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 [3, 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. [7, 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. [7, Theorems 6.1 and 7.2] and [32, Chapter 14.2]. The following lemma is a modification of the so-called random change of time result from [4, page 151], which is also called the continuity of composition in [32, 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 [32, 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 [30] 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 [30], 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 [30], 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 [27, 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 [3, 23]:

(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 [3]:

\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.