2. Model description and notation

We consider an RTL Markov-modulated fluid polling model with two queues. In our initial description we will not make any assumptions of symmetry between the two queues, to facilitate later presentation and discussions regarding these assumptions in Section 3.2. As alluded to in the introduction, our main contributions in the first part of the paper (Sections 3–4) concern the heavy-traffic limit for identical parameter settings for the two queues. Arguably, this is the most interesting case under heavy-traffic conditions, because—as we will make precise later—it ensures that the workload processes of the two queues obey a similar scaling when approaching the heavy-traffic saturation point, and consequently exhibit a nontrivial correlation. By contrast, in an asymmetric setting, one of the two queues will approach heavy traffic while the other remains bounded. In that case, the two workloads are asymptotically independent, and their joint heavy-traffic distribution can be obtained from the marginal scaled limit for the queue with heaviest load and the ordinary (non-heavy-traffic) marginal distribution of the lighter-loaded queue.

In the asymmetric setting, fluid enters queue j (say $Q_j$ ) at a constant rate of $\lambda_j>0$ , $j = 1, 2$ . There is a single server that serves both queues with constant rate $\mu_j>0$ , $j = 1, 2$ . A special feature of the model is that the server spends random amounts of time at each queue; these times are independent of the fluid content levels (workloads at the queues). In particular, when a queue becomes empty, the server will remain at that queue (although not providing any service), even if the other queue is not empty, until the expiration of the random visit time. We denote the length of a generic time interval for which the server resides at $Q_j$ by $T_j$ , $j = 1, 2$ . The periods $T_j$ are exponentially distributed with rate $c_j>0$ , $j = 1, 2$ . Upon completion of the residence time at $Q_j$ , the server instantaneously switches to the other queue $Q_{3-j}$ , $j=1,2$ ; i.e., there is no switch-over time. All residence times are independent.

To analyze the model under consideration, we let $V_j(t)$ denote the workload at $Q_j$ at time t, $t\geq0$ . Assuming $V_1(0) = u_0$ , $V_2(0) = l_0$ , and the server being at $Q_1$ at time 0, we can describe the workload at time $T_1$ and $T_1 + T_2$ as follows:

• In the interval $(0, T_1]$ the server is serving $Q_1$ ; therefore, the workload (the fluid content) at $Q_1$ decreases linearly as long as it is positive: $V_1(T_1) = \max\{0, u_0 + (\lambda_1 - \mu_1) T_1\}$ . During this time period, the workload at $Q_2$ increases linearly: $V_2(T_1) = l_0 + \lambda_2 T_1$ .

• Analogously, as explained above, in the interval $(T_1, T_1 + T_2]$ the server is serving $Q_2$ , so the workload at $Q_2$ decreases linearly as long as it is positive; hence $V_2(T_1 + T_2) = \max\{0, V_2(T_1) + (\lambda_2 - \mu_2) T_2\}$ . The workload at $Q_1$ increases linearly, so $V_1(T_1 + T_2) = V_1(T_1) + \lambda_1 T_2$ .

Stability condition. For the model under consideration, the stability condition states that both queues must have larger capacities than the respective loads imposed on them:

(1) \begin{equation}\rho_1< \frac{ c_2}{ c_1 + c_2}\; \text{ and }\; \rho_2< \frac{c_1}{c_1 + c_2},\end{equation}

with $\rho_j=\frac{\lambda_j}{\mu_j}$ , $j=1,2$ ; cf. [Reference Saxena, Boxma, Kapodistria and Núñez-Queija39].

3. Workload analysis

In Section 3.1, we first briefly focus on the stationary workload of the marginal queues (initially of $Q_1$ , and hence by identical arguments also of $Q_2$ ). In Section 3.2, we proceed with the joint workload analysis.

3.1. Marginal workload analysis

Let $V_1$ denote the stationary workload at an arbitrary epoch. From a special case of [Reference Kella and Whitt32, Section 5], and also from the analysis performed in [Reference Chen and Yao17], the marginal queue length distributions of the model under consideration are known. We include them here for completeness.

Theorem 1. The LST of the (marginal) workload of the first queue in stationarity under the stability condition (1) is given by

(2) \begin{align}\mathbb{E}\!\left(e^{- s V_1}\right) = \frac{1 + \frac{\rho_1 \mu_1}{c_1 + c_2} s}{1 + \frac{\rho_1 \mu_1}{c_2\!\left(1 - \frac{c_1}{c_2} \frac{\rho_1}{1 - \rho_1}\right)} s}.\end{align}

An equivalent formula holds for the LST of $V_2$ under the stability condition (1).

Remark 1. From the result of Theorem 1, it is evident that, for $\theta_1=\frac{c_1+c_2}{\rho_1\mu_1}$ , $\theta_2=\frac{c_2}{\rho_1\mu_1}\!\left(1-\frac{c_1}{c_2}\frac{\rho_1}{1-\rho_1}\right)$ , and with $\mathcal{E}_\theta\sim\text{exp}(\theta)$ , the following hold:

1. 1. $V_1+\mathcal{E}_{\theta_1}$ is distributed like $\mathcal{E}_{\theta_2}$ , with $ V_1,\mathcal{E}_{\theta_1}$ independent.

2. 2. The distribution of $V_1$ is a mixture of zero and $\mathcal{E}_{\theta_2}$ .

With the result of Theorem 1, we can study the behavior of the workload $V_1$ in heavy traffic, i.e., when $\rho_1 \uparrow \frac{c_2}{c_1 + c_2}$ .

Lemma 1. For $\rho_1 \uparrow \frac{c_2}{c_1 + c_2}$ ,

(3) \begin{eqnarray}\!\left(\frac{c_2}{c_1 + c_2}-\rho_1\right) V_1 \overset{d}{\longrightarrow} Z,\end{eqnarray}

where Z is an exponentially distributed random variable with mean $\frac{c_1c_2 \mu_1}{(c_1 + c_2)^3}$ .

Proof. Replacing s by $\!\left(\frac{c_2}{c_1 + c_2}-\rho_1 \right)s$ in (2) and taking the limit as $\rho_1 \uparrow \frac{c_2}{c_1 + c_2}$ yields

(4) \begin{equation} \lim\limits_{\rho_1 \uparrow \frac{c_2}{c_1 + c_2}}\mathbb{E}\!\left(e^{-\!\left(\frac{c_2}{c_1 + c_2}-\rho_1\right)sV_1}\right) = \frac{1}{ 1 + \frac{c_1c_2 \mu_1}{(c_1 + c_2)^3} s}.\end{equation}

Note that the right-hand side in (4) corresponds to the LST of an exponentially distributed random variable with mean $\frac{c_1c_2 \mu_1}{(c_1 + c_2)^3}$ .

3.2. Joint workload analysis

We now focus on the joint workload distribution, restricting ourselves to the symmetric case, i.e., $c_1 = c_2 = c, \lambda_1 = \lambda_2 = \lambda$ , and $\mu_1 = \mu_2 = \mu$ . A main stepping-stone in our analysis is the functional equation in (13) below. A corresponding functional equation can be derived for the asymmetric case (see also [Reference Saxena, Boxma, Kapodistria and Núñez-Queija39]), but for the purposes of this paper it suffices to show that the symmetric case leads to a complicated BVP that, although it can be solved, provides little probabilistic insight into the problem at hand. Our next step is to analyze it under heavy traffic, and, as explained earlier, the symmetric case is then the interesting one.

As a side remark, note that for both queues to reach heavy traffic simultaneously, it suffices to have $\lambda_1/\mu_1=\lambda_2/\mu_2$ if $c_1=c_2$ ; additionally demanding that $\lambda_1=\lambda_2$ (and hence $\mu_1=\mu_2$ ) amounts to choosing a different scaling unit for the workloads.

Calculation of $\mathbb{E}\!\left(e ^{-s_1 V_1(T_1) - s_2 V_2(T_1)}|V_1(0)=u_0,V_2(0)= l_0\right)\!.$

In this part, we calculate the LST of the joint workload distribution at time $T_1$ . From the observations listed above the stability condition (1) in Section 2, we obtain

(5) \begin{align} &\mathbb{E}\!\left(e^{-s_1 V_1(T_1) - s_2 V_2(T_1)}|V_1(0)=u_0,V_2(0)=l_0\right)\nonumber\\[5pt] &= c e^{-s_2 l_0} \!\left(\int_{t = 0}^{\frac{u_0}{\mu - \lambda}} e^{- ( s_2 \lambda + c )t} e^{-s_1 (u_0 + (\lambda - \mu)t)} \textrm{d}t + \int_{t=\frac{u_0}{\mu - \lambda}}^{\infty} e^{- ( s_2 \lambda + c)t} \textrm{d}t\right) \nonumber \\[5pt] &= \frac{c}{ c + s_2 \lambda} \left[\frac{s_2 \lambda + c}{s_2 \lambda + c + s_1(\lambda - \mu)} e^{-s_1 u_0 - s_2 l_0} + \frac{ s_1(\lambda - \mu)}{s_2 \lambda + c + s_1(\lambda - \mu)} e^{- \frac{c + s_2 \lambda}{\mu - \lambda}u_0 -s_2 l_0}\right] .\end{align}

Calculation of $\mathbb{E}\!\left(e^{-s_1 V_1(T_1) - s_2 V_2(T_1)}\right).$

Unconditioning, we obtain from (5) the following:

(6) \begin{align} &\mathbb{E}\!\left(e^{-s_1 V_1(T_1) -s_2 V_2(T_1)}\right) \nonumber \\[5pt] &=\frac{c}{ c + s_2 \lambda} \Big[\frac{s_2 \lambda + c}{s_2 \lambda + c + s_1(\lambda - \mu)} \mathbb{E}\!\left(e^{-s_1 V_1(0) - s_2 V_2(0)}\right) \nonumber \\[5pt] & \quad + \frac{ s_1(\lambda - \mu)}{s_2 \lambda + c+ s_1(\lambda - \mu)} \mathbb{E}\!\left(e^{- \frac{c + s_2 \lambda}{\mu - \lambda} V_1(0) - s_2 V_2(0)}\right) \Big].\end{align}

Formulation of a functional equation. Since we are interested in the symmetric case, we can formulate a functional equation corresponding to (6) by defining

\begin{equation*}\tilde{\nu}(s_1, s_2) \;:\!=\; \mathbb{E}\!\left(e^{-s_1 V_1(0) -s_2 V_2(0)}\right),\end{equation*}

and then by symmetry,

\begin{equation*}\tilde{\nu}(s_2, s_1) = \mathbb{E}\!\left(e^{-s_1 V_1(T_1) -s_2 V_2(T_1)}\right).\end{equation*}

Further, defining

(7) \begin{equation}f(s_1, s_2) \;:\!=\; s_1 \lambda + c+ s_2(\lambda - \mu),\end{equation}

we obtain

(8) \begin{align} \tilde{\nu}(s_2, s_1) = \frac{c}{f(s_2, s_1) } \tilde{\nu}(s_1, s_2) + \frac{c}{ f(s_2, s_1)} \frac{ s_1(\lambda - \mu)}{s_2 \lambda + c} \tilde{\nu}\!\left(\frac{s_2 \lambda + c}{\mu - \lambda}, s_2\right).\end{align}

Now substituting $s_1 = s_2$ gives

(9) \begin{align}\tilde{\nu}\!\left(\frac{s_2 \lambda + c}{\mu - \lambda}, s_2\right) = \frac{(2 \lambda - \mu)(s_2 \lambda + c)}{c (\lambda - \mu)} \tilde{\nu}(s_2, s_2).\end{align}

Combining (8) and (9) yields

(10) \begin{align}\tilde{\nu}(s_2, s_1) = \frac{c}{f(s_2, s_1) } \tilde{\nu}(s_1, s_2) + \frac{(2\lambda - \mu)s_1}{f(s_2, s_1)} \tilde{\nu}(s_2, s_2).\end{align}

By symmetry (after interchanging the indexes),

(11) \begin{align}\tilde{\nu}(s_1, s_2) = \frac{c}{f(s_1, s_2) } \tilde{\nu}(s_2, s_1) + \frac{(2\lambda - \mu)s_2}{f(s_1, s_2)} \tilde{\nu}(s_1, s_1).\end{align}

Combining (10) and (11), it follows that

(12) \begin{align}\tilde{\nu}(s_1, s_2) & = \frac{c^2}{f(s_1, s_2)f(s_2, s_1) } \tilde{\nu}(s_1, s_2) + \frac{c(2\lambda - \mu)s_1}{ f(s_1, s_2) f(s_2, s_1)} \tilde{\nu}(s_2, s_2) \nonumber \\[5pt] & \quad + \frac{(2\lambda - \mu)s_2}{f(s_1, s_2)} \tilde{\nu}(s_1, s_1),\end{align}

so that finally

(13) \begin{align} \frac{\tilde{k}(s_1, s_2)}{s_1s_2} \tilde{\nu}(s_1, s_2) = -2 \mu \!\left(\frac{1}{2} - \frac{\lambda}{\mu}\right)\left[\frac{c}{s_2}\ \tilde{\nu}(s_2, s_2) + \frac{f(s_1, s_2)}{s_1} \ \tilde{\nu}(s_1, s_1)\right],\end{align}

with $\tilde{k}(s_1, s_2) = f(s_1, s_2) f(s_2, s_1) - c^2$ and with $f(s_1, s_2)$ defined in Equation (7).

Equations of this type have been studied in the monograph [Reference Cohen and Boxma21]. There a solution procedure for the present problem is outlined, which amounts to the following global steps:

1. Step A. Consider the zeros of the kernel equation $\tilde{k}(s_1, s_2)$ that have $\textrm{Re}[s_1]$ , $\textrm{Re}[s_2] \geq 0$ . For such pairs $(s_1, s_2)$ , $\tilde{\nu}(s_1, s_2)$ is analytic, and hence, for those pairs, the left-hand side of (13) is equal to zero.

2. Step B. A suitable set of these zeros of the kernel may form a contour. The fact that the right-hand side of (13) is zero on that contour (the ‘boundary’), in combination with analyticity properties of $\tilde{\nu}(s_1,s_1)$ and $\tilde{\nu}(s_2,s_2)$ inside and/or outside that contour, can be used to formulate a Riemann or Riemann–Hilbert BVP. The solution of such a problem yields $\tilde{\nu}(s_1,s_1)$ and $\tilde{\nu}(s_2,s_2)$ . Then $\tilde{\nu}(s_1,s_2)$ follows via (13).

Unfortunately, the above steps do not constitute a simple, straightforward recipe. For example, several choices of zero pairs are possible in the present problem, and it is not a priori clear what is the best choice. Therefore, to successfully employ this boundary value method (BVM) requires a detailed investigation of the zeros of the kernel $\tilde{k}(s_1, s_2)$ of the functional equation.

In the appendix, we demonstrate the steps for a set of zero pairs that lie on an ellipse. For such pairs, the functional equation can be transformed (through a conformal mapping, which is explicitly expressed in the Jacobi elliptic function) to a Riemann–Hilbert BVP on the unit circle. Solving the reduced BVP and using the inverse of the conformal mapping produces the solution to the problem at hand. From the analysis presented in the appendix, it is made clear that this choice of the zero pairs leads to an intricate LST expression involving a complex integral on a smooth contour (ellipse). To further invert this LST expression requires cumbersome computational procedures. In addition to the numerical complications, because of the nature of the solution of the BVP, it is difficult to gain probabilistic insight into the problem at hand. For all of the aforementioned reasons, we instead focus on the heavy-traffic setting of the model and solve the resulting simpler BVP.

4. Heavy-traffic analysis of the joint workload distribution

In this section, we shall determine the heavy-traffic limit of the LST of the scaled joint workload distribution of the symmetric model in stationarity. In what follows, we use the functional equation (13). Let $\rho$ be the load on each of the two queues ( $\rho = {\lambda}/{\mu}$ ). We scale the functional equation by replacing $s_1$ by $\!\left({1}/{2} - \rho\right)s_1$ , and $s_2$ by $\!\left({1}/{2} - \rho\right)s_2$ . After dividing by $-2 \mu c$ in (13) and taking the limit $\rho \uparrow {1}/{2}$ , we obtain the following functional equation:

(14) \begin{align}\frac{\tilde{k}^{\star}(s_1, s_2)}{s_1 s_2} \tilde{\nu}^{\star}(s_1, s_2) = \frac{\tilde{\nu}^{\star}(s_1, s_1)}{s_1} + \frac{\tilde{\nu}^{\star}(s_2, s_2)}{s_2},\end{align}

where $\tilde{\nu}^{\star}(s_1, s_2) = \lim\limits_{\rho \uparrow {1}/{2}} \mathbb{E}(e^{-s_1\!\left({1}/{2} - \rho\right)V_1 -s_2\!\left({1}/{2} - \rho\right)V_2 })$ and

(15) \begin{align}\tilde{k}^{\star}(s_1, s_2) & = -\lim\limits_{\rho \uparrow {1}/{2}} \frac{1}{2 \mu c\!\left({1}/{2} - \rho\right)^2}\;\; \tilde{k}\!\left(\left({1}/{2} - \rho\right)s_1,\!\left({1}/{2} - \rho\right)s_2\right) \nonumber \\[5pt] & = s_1 + s_2 + \frac{\mu}{8c}\big(s_1 - s_2\big)^2 .\end{align}

There is one unknown function in the right-hand side of (14): $\tilde{\nu}^{\star}(s, s)$ . We calculate this unknown function using the BVM by applying Step A and Step B as discussed in Section 3.2.

Kernel analysis. To apply the BVM, one needs to investigate the zeros of the kernel $\tilde{k}^{\star}(s_1, s_2)$ . By setting $\tilde{k}^{\star}(s_1, s_2) = 0$ , we obtain

(16) \begin{align}s_2^{\pm}(s_1) = \frac{-1 + \frac{\mu}{4c} s_1 \pm \sqrt{1 - \frac{\mu}{c} s_1}}{\frac{\mu}{4c}}.\end{align}

Note that $s_2^{\pm}(s_1)$ has a single branching point at $s_1 = \frac{c}{\mu}$ . For real-valued $s_1$ with $s_1 > {c}/{\mu}$ , the function $s_2^{\pm}(s_1)$ is complex-valued. Letting $s_2^{\pm}(s_1) = u \pm i v$ , we obtain

(17) \begin{align}u^2 + v^2 = s_2^{+}(s_1) s_2^{-}(s_1) = \frac{\!\left(\!-\!1 + \frac{\mu}{4c} s_1\right)^2 - 1 + \frac{\mu}{c} s_1}{\!\left( \frac{\mu}{4c} \right)^2}\end{align}

and

(18) \begin{align}2 u = s_2^{+}(s_1) + s_2^{-}(s_1)= \frac{-1 + \frac{\mu}{4c} s_1}{\frac{\mu}{8c}}.\end{align}

Computing $s_1 = u + {3c}/{\mu}$ from the above equation and substituting it into (17), we have

(19) \begin{align}u^2 + v^2 = \frac{\!\left(\frac{\mu}{4c}\right)^2 u^2 - 1 + \frac{\mu}{c}u + 4}{\!\left(\frac{\mu}{4c}\right)^2}.\end{align}

Simplifying the above equation yields

(20) \begin{align}v^2 = \frac{16 c}{\mu} \!\left(u + \frac{3 c}{\mu}\right),\end{align}

which describes a parabola in the complex plane. We will restrict ourselves to the following set:

\begin{align*}E = \left\{(u, v) \in \!\left(-\frac{3c}{\mu}, \infty\right) \times \mathbb{R} \mid v^2 = \frac{16 c}{\mu} \!\left(u + \frac{3 c}{\mu}\right) \right\}.\end{align*}

BVM: solution of the functional equation (14). Now we take $s_1$ with $s_1 > {c}/{\mu}$ and $s_2^{\pm}(s_1) = u \pm iv$ , with $(u, v) \in E$ . For all such $(s_1, s_2^{\pm}(s_1))$ pairs, the right-hand side of (14) becomes zero, and hence, for all $s_2 = s_2^{\pm}(s_1)$ , we have

(21) \begin{equation}\frac{\tilde{\nu}^{\star}(s_2, s_2)}{s_2} = - \frac{\tilde{\nu}^{\star}(s_1, s_1)}{s_1}.\end{equation}

For $s_1 > {c}/{\mu}$ , the right-hand side of the above equation is real, thus yielding

(22) \begin{equation}\textrm{Re}\left[g(s_2) \right] = 0, \, \, \mbox{for} \, \, s_2 = s_2^{\pm}(s_1) = u \pm iv, \;\mbox{with}\; (u, v) \in E \backslash \{(0, 0)\},\end{equation}

where $g(s_2) = i \frac{\tilde{\nu}^{\star}(s_2, s_2)}{s_2}$ . Notice that $\tilde{\nu}^{\star}(s_2, s_2)$ is analytic for $\textrm{Re}[s_2] \geq 0$ . Below, in Lemma 3, we prove that $\tilde{\nu}^{\star}(s_2, s_2)$ is analytic on the strip $- {3 c}/{\mu} < \textrm{Re}[s_2] < 0$ . For clarity of exposition, we postpone the proof of this lemma until after Theorem 2, at which point we will have introduced all necessary notation.

We thus see that $g(s_2)$ is analytic inside the contour E, say on $E^+$ , except for $s_2=0$ , which is a pole in $E^+$ . The above problem now reduces to a Riemann–Hilbert problem with a pole, and with boundary E; see [Reference Cohen and Boxma21, Section I.3.3]. To transform it into a (standard) Riemann–Hilbert problem on the unit circle D, we define $\phi$ (with inverse $\psi$ ) to be a conformal mapping of the interior of the unit circle D onto the region bounded by E with normalization conditions $\phi\!\left(\!-\!1\right)=\infty$ , $\phi(0)=0$ , and $\phi\!\left(1\right)=-{3 c}/{\mu}$ . That allows us to translate the Riemann–Hilbert BVP on and inside E to the following simple Riemann–Hilbert BVP with a pole (cf. [Reference Cohen and Boxma21, Section I.3.3] and [Reference Boxma and van Houtum14, Section 6]): defining $h(w) \;:\!=\; g(\phi(w))$ , we obtain for $h(\!\cdot\!)$ on the unit circle D

(23) \begin{equation}\textrm{Re} [h(w)] = 0, \, \, w \in D \backslash \{0\},\end{equation}

with $h(\!\cdot\!)$ analytic in $D^+\backslash \{0\}$ and continuous in $D^+\cup D\backslash \{0\}$ . The solution of the BVP (23) is as follows; cf. [Reference Boxma and van Houtum14] and [Reference Cohen and Boxma21, Section I.3.3]:

(24) \begin{equation}h(w) = i \alpha + \beta w - \frac{\bar{\beta } }{w}, \, \, w \in D^+\cup D\backslash \{0\},\end{equation}

where $\alpha , \beta $ are constants that we will calculate explicitly in Theorem 2. This determines $g(x) = h(\psi(x))$ ; substituting it in the above equation, we obtain

(25) \begin{equation}g(s_2) = i \alpha + \beta \psi(s_2) - \frac{\bar{\beta } }{\psi(s_2)}, \, \, s_2 \in E^+\cup E\backslash\{0\},\end{equation}

where $\psi(\!\cdot\!)$ is the conformal mapping from the parabola E to the unit circle D. Since $g(s_2) = i {\tilde{\nu}^{\star}(s_2, s_2)}/{s_2}$ , we obtain for $\textrm{Re}[s_2] > - {3 c}/{\mu}$ the following:

(26) \begin{equation}\tilde{\nu}^{\star}(s_2, s_2) = \alpha s_2 - i\;\beta \; \psi(s_2) s_2 + i \frac{\bar{\beta } s_2}{\psi(s_2)}.\end{equation}

With that we have calculated the LST of the total workload in heavy traffic, based on the conformal mapping $ \psi(s_2)$ . Before materializing this in Theorem 2, we give an explicit expression for $\psi(s_2)$ in the next lemma. That will enable us to explicitly determine $\tilde{\nu}^{\star}(s_2, s_2)$ in Theorem 2.

Lemma 2. For $z \in \mathbb{C}$ , we have a conformal map $\psi(z)$ which maps the interior of the parabola (20) onto the interior of the unit circle D, and it is given explicitly as follows:

(27) \begin{equation}\psi(z) = \frac{1 - \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} z - 1}\right)}{1 + \sqrt{2} \cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} z - 1}\right)}.\end{equation}

Now we state the first main theorem of this section, in which we obtain an explicit expression for the total stationary workload LST in heavy traffic.

Theorem 2. The scaled total workload LST in heavy traffic is given by, for $\textrm{Re}[s] > -{3 c}/{\mu}$ ,

(28) \begin{equation}\lim\limits_{\rho \uparrow {1}/{2}} \mathbb{E}\!\left( e^{-s\left({1}/{2} - \rho\right)(V_1 +V_2) }\right ) = \frac{\pi}{4}\frac{\mu}{c} \frac{s}{\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s - 1}\right)}.\end{equation}

Proof. Substituting $\psi(z)$ from Lemma 2 in (26) yields

(29) \begin{align}\tilde{\nu}^{\star}(s, s) = \alpha s - i \;\beta \;\!\left(\frac{1 - \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}{1 + \sqrt{2} \cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}\right)s \nonumber \\[5pt] \quad + i \;\bar{\beta }\;\!\left(\frac{1 + \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}{1 - \sqrt{2} \cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}\right)s.\end{align}

Since $\tilde{\nu}^{\star}(0, 0) = 1$ , we obtain from the above equation $\bar{\beta } = \frac{\pi}{16}\frac{\mu}{c} i$ , and since $\tilde{\nu}^{\star}(\infty, \infty) = 0$ , we obtain $\alpha = -\frac{\pi}{8}\frac{\mu}{c}$ . Substituting the values of $\alpha $ , $\beta $ , and $\bar{\beta }$ into the above equation and thereafter simplifying it, we obtain

(30) \begin{align} & \tilde{\nu}^{\star}(s, s) \nonumber \\[5pt] &= -\frac{\pi}{8}\frac{\mu}{c} \left[1 + \frac{1}{2} \frac{1 - \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}{1 + \sqrt{2} \cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)} + \frac{1}{2} \frac{1 + \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}{1 - \sqrt{2} \cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s - 1}\right)}\right] s.\end{align}

The theorem follows after some further simplifications and using the trigonometric square formula $\cosh^2 x=\!\left(\cosh (2x)+1\right)/2$ .

It is now convenient to formulate and prove the postponed Lemma 3. As we discussed in Step A of Section 3.2, we are interested in finding the LST in the domain $\textrm{Re}[s_2] \geq 0$ . In the previous theorem, we calculated the LST expression $\tilde{\nu}^{\star}(s_2, s_2) $ in $\textrm{Re}[s_2] > -{3 c}/{\mu}$ . We want to show that $\tilde{\nu}^{\star}(s_2, s_2) $ is analytic in the strip $-{3 c}/{\mu} < \textrm{Re}[s_2] < 0$ . From (28), we have an explicit expression, and it is sufficient to show that the denominator $\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right)$ has no zeros on that strip.

Proof. In the proof of Lemma 2, we observed that $\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right)$ is a conformal mapping for $\textrm{Re}[s_2] > -{3 c}/{\mu}$ , and hence it is an analytic function for $\textrm{Re}[s_2] > -{3 c}/{\mu}$ . Moreover, the reciprocal of this analytic function is also analytic (see [Reference Brown and Churchill16, p. 74]) if the denominator has no zeros in that domain. To show that the denominator $\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right)$ has no zeros in $-{3 c}/{\mu} < \textrm{Re}[s_2] < 0$ , we solve

(31) \begin{equation}0 = \cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right) = \frac{e^{\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}} + e^{-\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}}}{2},\end{equation}

so

(32) \begin{equation}e^{\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}} = e^{\pi i -\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}},\end{equation}

and hence

\begin{align*}\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1} = \pi i -\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1} + 2 \pi n i, \qquad n \in \mathbb{Z}.\end{align*}

This implies that the zeros of the function $\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right)$ are $s_2 = \frac{c}{\mu}\!\left(1 - (2n + 1)^2\right)$ , $n \in \mathbb{Z}$ . There are two different cases for the zeros that we need to discuss: (i) when $n = 0$ or $n = -1$ , we have $s_2 = 0$ , and in this case we know $\tilde{\nu}^{\star}(s_2, s_2)$ is 1; (ii) when $n \in \mathbb{Z} \backslash \{0, -1\}$ , we have $s_2 = \frac{c}{\mu} \!\left(1 - (2n + 1)^2\right) \leq - {8c}/{\mu} < -{3 c}/{\mu}$ , which concludes the proof of the claim of the lemma.

Note that instead of working directly with the roots appearing in the simplified Equation (28), one could consider the roots of the two denominators appearing in Equation (30), i.e., the zeros of $1\pm \sqrt{2}\cosh\!\left(\frac{\pi}{4}\sqrt{\frac{\mu}{c} s_2 - 1}\right) $ .

Equivalently, one can prove the analytic continuation using Euler’s formula (cf. [Reference Biane, Pitman and Yor8, Equation (3.3)]), which converts the hyperbolic cosine into an infinite product (we use this approach to rewrite Equation (28) as an infinite product expansion; cf. Equation (41)). Using the infinite product expansion, it becomes evident that, in the domain $ \textrm{Re}[s_2]>-{8 c}/{\mu}$ , the denominator has no roots.

The LST of the total workload lends itself to explicitly determining the heavy-traffic stationary workload distribution, as shown in the following lemma.

Lemma 4. With $f_{({1}/{2}- \rho)(V_1 + V_2)}(\!\cdot\!)$ the probability density function of the scaled total workload $({1}/{2}- \rho)(V_1 + V_2)$ , we have

(33) \begin{align} \lim\limits_{\rho \uparrow {1}/{2}} f_{({1}/{2}- \rho)(V_1 + V_2)}(x) =\sum_{n=1}^\infty(\!-\!1)^{n+1}(2n+1)\frac{c}{\mu}\!\left((2n+1)^2-1\right) e^{-\frac{c}{\mu}\!\left((2n+1)^2-1\right)x},\qquad x>0.\end{align}

Moreover, the limiting distribution as $\rho\uparrow 1/2$ of the scaled total workload $({1}/{2}- \rho)(V_1 + V_2)$ is infinitely divisible and is distributed like

(34) \begin{align}\sum_{n=1}^{\infty}\frac{\mathcal{E}_n}{\frac{c}{\mu}\!\left((2n+1)^2-1\right)},\end{align}

where $\{\mathcal{E}_n\}_{n\in\mathbb{N}}$ is a sequence of independent and identically exponentially distributed random variables with rate 1.

The infinite divisibility of the scaled total workload distribution is a consequence of the infinite divisibility of the exponential distribution.

Before proceeding with the proof of Lemma 4, we review the relevant results needed in the remark below.

Consider a random variable defined as

(35) \begin{align}C=\frac{2}{\pi^2}\sum_{n=1}^{\infty}\frac{\mathcal{E}_n}{(n-{1}/{2})^2},\end{align}

with $\{\mathcal{E}_n\}_{n\in\mathbb{N}}$ a sequence of independent and identically exponentially distributed random variables with rate 1. Then

(36) \begin{align}\mathbb{E}\!\left(e^{-sC}\right)=\mathbb{E}\!\left(\prod_{n=1}^{\infty}e^{-\frac{2s}{\pi^2 (n-{1}/{2})^2}\mathcal{E}_n}\right)=\frac{1}{\prod_{n=1}^{\infty}\!\left(1+\frac{2s}{\pi^2 (n-{1}/{2})^2}\right)}=\frac{1}{\cosh \sqrt{2s}},\end{align}

where the last equality is known as Euler’s formula; cf. [Reference Biane, Pitman and Yor8, Equation (3.3)]. Moreover, using the Mittag-Leffler expansion, based on the poles of the right-hand side of Equation (36), yields

\begin{align*}\mathbb{E}\!\left(e^{-sC}\right)=\pi\sum_{n=1}^\infty \frac{(\!-\!1)^n(n-{1}/{2})}{s+(n-{1}/{2})^2\pi^2/2};\end{align*}

cf. [Reference Ciesielski and Taylor19, Equation (2.21)]. Noting that

\begin{align*}\frac{1}{s+(n-{1}/{2})^2\pi^2/2}=\int_{x=0}^{\infty} e^{-sx}e^{- (n-{1}/{2})^2\pi^2x/2}\textrm{d}x,\end{align*}

this last expression yields the density function of the random variable C; more concretely,

(37) \begin{align}f_{C}(x)=\pi \sum_{n=1}^\infty (\!-\!1)^n (n-{1}/{2})e^{-(n-{1}/{2})^2\pi^2 x/2},\qquad x>0.\end{align}

Moreover, expressions equivalent to (37) can be produced using the reciprocal relation $f_C(x)=\!\left(\frac{2}{\pi x}\right)^{3/2}f_C\!\left(\frac{4}{\pi^2 x}\right)$ ; cf. [Reference Biane, Pitman and Yor8, Table 1 (continued), Row 5]. This immediately implies that

\begin{align*}f_{C}(x)=\sqrt{\frac{2}{\pi x^3}}\sum_{n=1}^\infty (\!-\!1)^n (2n-1)e^{-(2n-1)^2/2x},\qquad x>0;\end{align*}

see [Reference Biane, Pitman and Yor8, Equation (3.11)].

As stated in [Reference Ciesielski and Taylor19, p. 441], this turns out to be the density for the maximum displacement of a one-dimensional standard Brownian motion in a fixed time interval, or, as stated in [Reference Biane, Pitman and Yor8, Table 2, Row 3], the density of the hitting time of 1 of the one-dimensional standard Brownian motion with reflection at 0.

To further understand the infinite divisibility of the hitting time, the interested reader may refer to [Reference Feller25, p. 550], where the idea relies on the fact that the hitting time from 0 to 1 can be divided into the hitting time from 0 to any point in the interval (0, 1) plus the independent (by the strong Markov property) hitting time from that point to 1. By putting more and more points between 0 and 1, one can express the hitting time as the limit of a null triangular array, which gives rise to the infinite divisibility property expressed in (35).

A similar approach can be applied for the random hitting time of a one-dimensional standard Brownian motion with drift $\mu\geq 0$ to $\{\pm 1\}$ , say C ^′. As shown in [Reference Kent34, Theorem 7.1], the random hitting time has the representation

(38) \begin{align}C'=2\sum_{n=1}^{\infty}\frac{\mathcal{E}_n}{\mu^2+\pi^2(n-{1}/{2})^2},\end{align}

and one shows, by performing the same computations as in (36), that it has the following LST:

\begin{align*}\mathbb{E}\!\left(e^{-sC'}\right)=\frac{\cosh \mu}{\cosh \sqrt{2s+\mu^2}}.\end{align*}

Proof of Lemma 4.We express the LST (28) as an infinite product of LSTs of independent exponentially distributed random variables. For this purpose, we need the following two identities. First,

(39) \begin{align}\frac{1}{\cosh \sqrt{s}}=\prod_{n=1}^\infty \!\left(1+\frac{s}{\pi^2(n-{1}/{2})^2}\right),\end{align}

which is Euler’s formula; cf. [Reference Biane, Pitman and Yor8, Equation (3.3)]. Moreover,

\begin{align*}\cos \pi s =\prod_{n=0}^\infty \!\left(1-\!\left(\frac{s}{n+{1}/{2}}\right)^2\right).\end{align*}

From this last equation, by taking out the $n=0$ term, we can show that

(40) \begin{align}\prod_{n=1}^\infty \!\left(1-\!\left(\frac{{1}/{2}}{n+{1}/{2}}\right)^2\right)=\lim_{s\to1/2}\frac{\cos \pi s}{1-4s^2}=\frac{\pi}{4}.\end{align}

Using (39) and (40) yields, after straightforward computations, that

(41) \begin{align}\frac{\pi}{4}\frac{\mu}{c} \frac{s}{\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s - 1}\right)} &= \frac{\mu s}{c}\frac{\prod_{n=1}^\infty \!\left(1-\!\left(\frac{{1}/{2}}{n+{1}/{2}}\right)^2\right)}{\prod_{n=1}^\infty \!\left(1+\frac{\frac{\pi^2 \mu s}{4c}-\frac{\pi^2}{4}}{\pi^2(n-{1}/{2})^2}\right)}\nonumber\\[5pt] &= \frac{\prod_{n=2}^\infty \!\left(1-\!\left(\frac{{1}/{2}}{n-1/2}\right)^2\right)}{\prod_{n=2}^\infty \!\left(1+\frac{\frac{ \mu s}{4c}-\frac{1}{4}}{(n-{1}/{2})^2}\right)}\nonumber\\[5pt] &=\prod_{n=1}^\infty\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)}.\end{align}

Note that the last equality reveals that the LST at hand is associated with the random variable of Equation (34). For convenience, we shall denote the random variable of Equation (34) by $\tilde{C}$ .

We now turn our attention to the computation of the density function. Note that the conventional approach to producing the density function (based on a meromorphic expansion) does not work, as the corresponding (meromorphic) series diverges. We shall overcome this by following the approach of [Reference Ciesielski and Taylor19]. More concretely, we consider, for $N>0$ ,

(42) \begin{align}\prod_{n=1}^N\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)} &=\sum_{n=1}^N\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)}\prod\limits_{k=1,k\neq n}^N\frac{k(k+1)}{k(k+1)-n(n+1)}\nonumber\\[5pt] &=\sum_{n=1}^N\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)}\frac{(\!-\!1)^{n+1} (2 n+1) N! (N+1)!}{(N-n)! (N+n+1)!}\nonumber\\[5pt] &=\sum_{n=1}^\infty(\!-\!1)^{n+1} (2 n+1)\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)}\nonumber\\[5pt] & \quad\quad\quad\times\frac{ (N-n+1)\cdots N }{(N+2)\cdots (N+n+1)}\unicode{x1d7d9}_{\{n\leq N\}},\end{align}

by taking partial fractions and noting that $(2n+1)^2-1=4n(n+1)$ . Note that, as $N\to\infty$ , the left-hand side of Equation (42) converges to (41). The Laplace transform on the right-hand side of Equation (42) can easily be inverted, from which we obtain that the density function of (41) is given by

(43) \begin{align}\lim_{N\to\infty}\sum_{n=1}^\infty (\!-\!1)^{n+1} (2 n+1) \frac{c}{\mu}\!\left((2n+1)^2-1\right) e^{-\frac{c}{\mu}\!\left((2n+1)^2-1\right)x}\nonumber\\[5pt] \times\frac{ (N-n+1)\cdots N }{(N+2)\cdots (N+n+1)}\unicode{x1d7d9}_{\{n\leq N\}},\qquad x>0.\end{align}

Applying the dominated convergence theorem immediately yields Equation (33), as the terms (with respect to N) inside the series are bounded,

\begin{align*}\frac{ (N-n+1)\cdots N }{(N+2)\cdots (N+n+1)}\leq 1,\qquad \forall \ 1\leq n\leq N,\end{align*}

as

\begin{align*}&\lim_{N\to\infty} (\!-\!1)^{n+1} (2 n+1) \frac{c}{\mu}\!\left((2n+1)^2-1\right) e^{-\frac{c}{\mu}\!\left((2n+1)^2-1\right)x}\frac{ (N-n+1)\cdots N }{(N+2)\cdots (N+n+1)}\unicode{x1d7d9}_{\{n\leq N\}}\\[5pt] &\quad =(\!-\!1)^{n+1} (2 n+1) \frac{c}{\mu}\!\left((2n+1)^2-1\right) e^{-\frac{c}{\mu}\!\left((2n+1)^2-1\right)x},\end{align*}

and the series

\begin{align*}\sum_{n=1}^\infty (\!-\!1)^{n+1} (2 n+1) \frac{c}{\mu}\!\left((2n+1)^2-1\right) e^{-\frac{c}{\mu}\!\left((2n+1)^2-1\right)x}\end{align*}

converges for $x>0$ . From [Reference Hirschman and Widder30, Theorem 8.2, p. 60], it follows that (33) is indeed the density function in question. This is intuitively validated by noting that, if we were to directly use the Mittag-Leffler expansion, based on the poles of the right-hand side of Equation (41), this would yield

\begin{align*}\frac{\pi}{4}\frac{\mu}{c} \frac{s}{\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s - 1}\right)} =\sum_{n=1}^\infty (\!-\!1)^{n+1}(2n+1)\frac{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}{s+\frac{c}{\mu}\!\left((2n+1)^2-1\right)}.\end{align*}

However, the right-hand side of the above equation does not converge, but still yields the same result as in (33). This was noticed and commented upon in [Reference Ciesielski and Taylor19, p. 441].

We now state the most important result of this section. In (28), we have computed an explicit expression for the scaled total workload LST in heavy traffic. In the following theorem, we give an explicit expression for the scaled joint workload LST in heavy traffic.

Theorem 3. For $\textrm{Re}[s_j] > - {3c}/{\mu}$ , $j=1,2$ , the scaled joint workload LST in heavy traffic is given by

(44) \begin{align} &\lim\limits_{\rho \uparrow {1}/{2}} \mathbb{E}(e^{-s_1\!\left({1}/{2} - \rho\right)V_1 -s_2\!\left({1}/{2} - \rho\right)V_2 }) \nonumber\\[5pt] & \qquad = \frac{\pi}{4}\frac{\mu}{c} \frac{s_1 s_2}{\tilde{k}^{\star}(s_1, s_2)}\left[ \frac{1}{\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_1 - 1}\right)} + \frac{1}{\cosh\!\left(\frac{\pi}{2}\sqrt{\frac{\mu}{c} s_2 - 1}\right)} \right],\end{align}

where $\tilde{k}^{\star}(s_1, s_2) = s_1 + s_2 + \frac{\mu}{8c}\big(s_1 - s_2\big)^2$ .

Remark 3. Notice that if we let $s_2 \to 0$ in (44), the right-hand side tends to

\begin{align*}\frac{s_1}{\tilde{k}^{\star}(s_1, 0)} = \frac{1}{1 + \frac{\mu}{8c} s_1},\end{align*}

which is the heavy-traffic limit LST of the marginal workload as given in Lemma 1.

As a corollary we compute the first and second stationary moments of the joint workload in heavy traffic.

Corollary 1. For $j = 1, 2$ , it holds that

\begin{align*}\mathbb{E}\!\left(\textrm{lim}_{\rho \uparrow 1/2} \!\left({1}/{2} - \rho\right) V_j\right) & =\frac{\mu}{8 c}, \\[5pt] \mathbb{E}\!\left(\textrm{lim}_{\rho \uparrow 1/2} \!\left({1}/{2} - \rho\right)^2 V_j^2\right) &=\frac{\mu^2}{32 c^2}, \\[5pt] \mathbb{E}\!\left(\textrm{lim}_{\rho \uparrow 1/2}\!\left({1}/{2} - \rho\right)^2 V_1 V_2\right) &= \frac{\mu^2}{32 c^2} \frac{\pi^2 - 9}{3}, \\[5pt] \mathbb{R}\!\left(\textrm{lim}_{\rho \uparrow 1/2}\!\left(\left({1}/{2} - \rho\right)V_1, \!\left({1}/{2} - \rho\right)V_2\right) \right) &= \frac{2}{3} \pi^2 - 7 \approx -0.4203,\end{align*}

where $\mathbb{R}(\cdot, \cdot)$ is the correlation coefficient.

Proof. The marginal moments of the workload $V_j$ , $j = 1, 2$ , in heavy traffic are computed directly from Lemma 1. Equivalently, the expression (34) can be used to compute the moments, namely

\begin{align*}\mathbb{E}(\tilde{C}) &=\sum_{n=1}^{\infty}\frac{1}{\frac{c}{\mu}\!\left((2n+1)^2-1\right)}=\frac{\mu}{4c},\\[5pt] \mathbb{V}\textrm{ar}(\tilde{C}) &=\sum_{n=1}^{\infty}\frac{1}{\frac{c^2}{\mu^2}\!\left((2n+1)^2-1\right)^2}=\frac{\mu^2}{16c^2}\frac{\pi^2-9}{3},\end{align*}

with $\tilde{C}$ denoting the scaled total workload $\lim\limits_{\rho \uparrow {1}/{2}} ({1}/{2}- \rho)(V_1 + V_2)$ . The joint moment of $\lim\limits_{\rho \uparrow {1}/{2}} ({1}/{2} - \rho)^2 V_1 V_2$ is computed by differentiating the LST expression (44) with respect to $s_1$ and $s_2$ .

5. Numerical results

In this section, we verify the heavy-traffic results obtained above via simulations. Note that there are situations where simulation is not very efficient, and one such scenario appears in the heavy-traffic analysis of queueing models; see, e.g., [Reference Asmussen6]. Here it has been noted repeatedly that the standard simulation methods do not perform satisfactorily, one main problem being that the run lengths need to be exceedingly large to obtain even moderate precision. We have conducted simulations to validate our findings. One expects that as $\rho \uparrow {1}/{2}$ , the correlation coefficient tends to the exact correlation coefficient $\mathbb{R}\!\left(\textrm{lim}_{\rho \uparrow 1/2}\!\left(({1}/{2} - \rho)V_1,({1}/{2} - \rho)V_2\right)\right) = -0.4203$ . For the parameters $c = 0.1$ and $\rho = 0.49$ , we perform 1000 batches of MaxTime ( $2\times10^7$ ) simulations and calculate the correlation coefficient, the lower limit (LL), and the upper limit (UL) of the 95% confidence interval using the 1000 samples of the correlation coefficients. The runtime of each simulation is approximately 2 hours.

Table 1. Simulated correlation coefficient. The theoretical value for $\rho \to {1}/{2}$ is $-0.4203$ . By properties of the correlation coefficient, $\mathbb{R}(V_1, V_2)$ equals $\mathbb{R}(({1}/{2} - \rho)V_1, ({1}/{2} - \rho)V_2)$ .

From Table 1, we observe that as $\rho$ approaches $0.5$ from below, the simulation result approaches $ \mathbb{R}\!\left(\textrm{lim}_{\rho \uparrow 1/2}\!\left(({1}/{2} - \rho)V_1,({1}/{2} - \rho)V_2\right)\right) = -0.4203$ . We also observe that the upper and lower limits of the confidence interval increase as $\rho$ approaches ${1}/{2}$ .

For further validation of our heavy-traffic results, we plot the empirical cumulative distribution function of the scaled total workload in heavy traffic. For the parameters mentioned above, we first compute the inverse Laplace transform of the expression given in (28) by using Talbot’s method [Reference Abate and Whitt1] in MATLAB, then compare it with the simulation results. The simulations are performed for the loads $\rho = 0.2, 0.3, 0.4,$ and $0.49$ . Each simulation is performed for MaxTime ( $1\times10^9$ ), which takes approximately 1 hour. In Figure 1, one can see that as $\rho$ approaches $0.49$ the simulation results also approach the results obtained from the empirical cumulative distribution computed numerically from the inverse LST of the expression given in (28), i.e., ECDF_invLST.

Figure 1. Empirical cumulative distribution of the scaled total workload in heavy traffic.

6. Process limit in heavy traffic

Our main result so far, in Theorem 3, established the heavy-traffic limit of the stationary joint distribution of the scaled workloads in the symmetric case. In this section, we investigate the heavy-traffic limit of the entire process of scaled workloads, under less restrictive assumptions on the input processes and the server switching process. We show that the stationary distribution of this limit process corroborates the limit distribution of Theorem 3, establishing that the time-stationary limit and the heavy-traffic limit can be interchanged. Similar interchanges of limits have been previously demonstrated for different models in [Reference Gamarnik and Zeevi26] and [Reference Zhang and Zwart42].

As mentioned above, for the analysis in this section, we relax our assumptions regarding the input processes and the server switching process between the queues. For the two queues, we assume Lévy subordinator inputs instead of constant fluid flows, and the server visit periods form an alternating renewal process with possibly dependent consecutive visiting periods to server 1 and server 2. In the following, we make our assumptions precise.

We start with the server switching process. Specifically, we consider a sequence of independent and identically distributed nonnegative random pairs $\{(T_1(k),T_2(k)),\,k\ge 1\}$ distributed like $(T_1,T_2)$ , where $\mathbb{E}(T_j)^2$ , $j=1,2,$ are assumed finite (the marginal distributions of the $T_j$ are no longer assumed to be exponential). As before, ${1}/{c_j}=\mathbb{E} (T_j)$ , and we set $\sigma_j^2=\mathbb{V}\textrm{ar}[T_j]$ . The covariance between consecutive visit periods to queue 1 and queue 2 is denoted by $\zeta=\mathbb{C}\textrm{ov}(T_1,T_2)$ .

Let $S_0=0$ , $S_n=\sum_{k=1}^n(T_1(k)+T_2(k))$ for $n\ge 1$ , and set $I(t)=1$ if $t\in\bigcup_{n=0}^\infty [S_n,S_n+T_1(n+1))$ and $I(t)=0$ otherwise. Assuming that $T_1+T_2$ is not almost surely zero, with $p_1\;:\!=\;{c_2}/{(c_1+c_2)}$ it is well known that $\frac{1}{t}\int_0^tI(u)\textrm{d}u\to p_1$ almost surely, and it is also known that

(45) \begin{equation}\frac{1}{\sqrt{n}}\int_0^{nt}(I(u)-p_1)\textrm{d}u\end{equation}

converges weakly (in $D[0,\infty)$ endowed with the Skorokhod $J_1$ -topology) to a zero-drift Brownian motion with variance given by

(46) \begin{align}\sigma^2 &=\frac{\mathbb{V}\textrm{ar}\!\left(\int_0^{T_1+T_2}I(u)\textrm{d}u-p_1(T_1+T_2)\right)}{\mathbb{E}(T_1+T_2)}=\frac{\mathbb{V}\textrm{ar}(T_1-p_1(T_1+T_2))}{\mathbb{E}(T_1+T_2)}\nonumber\\[5pt] &=\frac{c_1^2\sigma_1^2-2c_1c_2\zeta+c_2^2\sigma_2^2}{{(c_1+c_2)^3}}c_1c_2.\end{align}

For a central limit theorem version of this, see [Reference Takács41, Reference Hew29]. This central limit version may also be concluded from [Reference Asmussen7, Theorem 3.2, p. 178]. The functional limit theorem may be concluded from, e.g., [Reference Glynn and Whitt27, Reference Glynn and Whitt28]. Let us denote this Brownian motion by $\sigma W(t)$ , where $\{W(t),t\ge t\}$ denotes a Wiener process (standard Brownian motion).

Next we describe the input processes into the two queues, which we assume to be independent of the server switching process just described. We no longer assume that the input processes are constant fluid flows, but instead let the input into $Q_j$ be a Lévy process $\{J_j(t), t\geq0\}$ , $j=1,2$ . To be precise, we assume that $\{J(t)\equiv (J_1(t),J_2(t)),\,t\ge 0\}$ is a bivariate subordinator with Laplace exponent $-\eta(s)$ , where, for $(s_1,s_2) \in\mathbb{R}_+^2$ ,

(47) \begin{align}\eta(s_1,s_2)=b_1s_1+b_2s_2+\int_{\mathbb{R}_+^2}(1-e^{-s_1x_1-s_2x_2})\Pi(\textrm{d}x_1, \textrm{d}x_2).\end{align}

Here $(b_1,b_2)\in\mathbb{R}^2_+$ , and $\Pi$ is the Lévy measure satisfying $\int_{\mathbb{R}_+^2}x_j\wedge 1\,\Pi(\textrm{d}x_1, \textrm{d}x_2)<\infty$ for $j=1,2$ . However, here we actually assume that $\int_{\mathbb{R}^2_+}x_j^2\Pi(\textrm{d}x_1, \textrm{d}x_2)<\infty$ , which is equivalent to the assumption that $\mathbb{E} (J_j^2(1))<\infty$ for $j=1,2$ . Consistent with our earlier notation, for $j,i\in\{1,2\}$ we write

(48) \begin{align}\lambda_j =\mathbb{E}(J_j(1))&=b_j+\int_{\mathbb{R}_+^2}x_j\Pi(\textrm{d}x_1,\textrm{d}x_2)=\frac{\partial\eta}{\partial s_j}(0\!+,0\!+\!), \nonumber\\[5pt] \sigma_{ji} &=\mathbb{C}\textrm{ov}(J_j(1),J_i(1))=\int_{\mathbb{R}_+^2}x_jx_i\Pi(\textrm{d}x_1, \textrm{d}x_2)=\frac{-\partial^2\eta}{\partial s_j\partial s_i}(0\!+,0\!+\!)\ ,\end{align}

and $\Sigma=(\sigma_{ji})_{j,i\in\{1,2\}}$ is the covariance matrix. Then $n^{-1/2}\left ( J_1(nt)-\lambda_1 nt,J_2(nt)-\lambda_2 nt\right)$ converges weakly to a zero-mean (two-dimensional) Brownian motion with covariance matrix $\Sigma$ . Let us denote this Brownian motion by $\{B(t)\equiv (B_1(t),B_2(t)),\,t\geq0\}$ . Since we have assumed that the processes $\{(T_1(k),T_2(k)),\,k\ge 1\}$ and $\{J(t), \, t \geq 0\}$ are independent, the Brownian motions W (one-dimensional) and B (two-dimensional) are independent as well.

We are now ready to describe the buffer content process of the two queues. The cumulative input to $Q_j$ up to time t is $J_j(t)$ , $j=1,2$ . At instances where $I(t)=1$ (resp., $I(t)=0$ ), the server is working at $Q_1$ (resp., $Q_2$ ) at a rate of $\mu_1$ (resp., $\mu_2$ ). If we let $p_2=1-p_1$ and define the free processes

(49) \begin{align}X_1(t) &=J_1(t)-\mu_1\int_0^tI(u)\textrm{d}u\nonumber\\[5pt] &=J_1(t)-\lambda_1t+(\lambda_1-p_1\mu_1)t-\mu_1\int_0^t(I(u)-p_1)\textrm{d}u, \nonumber\\[5pt] X_2(t) &=J_2(t)-\mu_2\int_0^t(1-I(u))\textrm{d}u\\[5pt] &=J_2(t)-\lambda_2t+(\lambda_2-p_2\mu_2)t+\mu_2\int_0^t(I(u)-p_1)\textrm{d}u,\nonumber\end{align}

then the buffer content process associated with the jth station ( $j=1,2$ ) is given by the (continuous) functional

(50) \begin{equation}V_j(t)=X_j(t)-\inf_{0\le u\le t}X_j(u).\end{equation}

As is natural in our model, we assume that $\lambda_j>0$ for $j=1,2$ and that $0<p_1<1$ . Let us replace $(\mu_1,\mu_2)$ by a sequence $(\mu_1^n,\mu_2^n)$ such that, as $n\to\infty$ ,

(51) \begin{equation}\sqrt{n}(p_1\mu_1^n-\lambda_1,p_2\mu_2^n-\lambda_2)\to (\theta_1,\theta_2).\end{equation}

Although not necessary at this point, for later considerations we will assume that $\theta_j>0$ , $j=1,2$ . For each value of n, $X_j^n(t)$ is the resulting free process with service rates $\mu_j^n$ , and $V_j^n(t)=X_j^n(t)-\inf_{0\le u\le t}X_j^n(u)$ is the corresponding buffer content process of $Q_j$ , $j=1,2$ . Observing that $\mu_j^n\to {\lambda_j}/{p_j}$ , it follows that $n^{-1/2}X_j^n(nt)$ converges weakly to

\begin{align*}X_1^\star(t) =-\theta_1t+B_1(t)-\frac{\lambda_1\sigma}{p_1}W(t), \\[5pt] X_2^\star(t) =-\theta_2t+B_2(t)+\frac{\lambda_2\sigma}{p_2}W(t).\end{align*}

In particular, the covariance matrix of the limiting Brownian motion is given by

\begin{align*}\Sigma^\star =\begin{pmatrix}\sigma_{11}+\frac{\lambda_1^2\sigma^2}{p_1^2} & \quad\sigma_{12}-\frac{\lambda_1\lambda_2\sigma^2}{p_1p_2}\\[10pt] \sigma_{12}-\frac{\lambda_1\lambda_2\sigma^2}{p_1p_2} & \quad\sigma_{22}+\frac{\lambda_2^2\sigma^2}{p_2^2}\end{pmatrix} .\nonumber\end{align*}

By the continuous mapping theorem it also follows that $n^{-1/2}V_j^n(nt)$ converges weakly to $V_j^\star(t)$ with $V_j^\star(t)=X_j^\star(t)-\inf_{0\le u\le t}X_j^\star(u)$ , $j=1,2$ .

In the previous sections, we considered the special case $J_j(t)=\lambda_j t$ , so that $\sigma_{ji}=0$ for $j,i=1,2$ . If in addition $\sigma>0$ (note that this assumption only excludes the case in which $T_1(k)/T_2(k)$ is a fixed constant), we can define $\hat X_j^\star=\frac{p_j}{\lambda_j\sigma}X_j^\star$ and $\hat \theta_j=\frac{p_j}{\lambda_j\sigma}\theta_j$ , $j=1,2$ . This results in

(52) \begin{align}\hat X_1^\star(t)=-\hat \theta_1 t-W(t), \qquad\hat X_2^\star(t)=-\hat \theta_2 t+W(t).\end{align}

Finally, defining $\hat V_j^\star(t)=\hat X_j^\star(t)-\inf_{0\le u\le t}\hat X_j^\star(u)$ , we observe that $\hat V_j^\star(t)=\frac{p_j}{\lambda_j\sigma}V_j^\star(t)$ , so that to study the stationary behavior of $V_j^\star(t)$ it suffices to study that of $\hat V_j^\star(t)$ . From now on it will be necessary that $\hat \theta_j>0$ for $j=1,2$ , which is ensured by our earlier assumption that $\theta_j>0$ .

Let us first observe that for $s \in\mathbb{R}_+^2$ (actually for all $s \in\mathbb{R}^2$ ), from (52), it follows after straightforward computations that

\begin{align*}\hat{k}(s_1,s_2)\equiv\log \mathbb{E}\left[e^{-s_1\hat X_1^\star(1)-s_2\hat X_2^\star(1)}\right]=\hat \theta_1 s_1+\hat \theta_2 s_2+\frac{1}{2}(s_1-s_2)^2.\end{align*}

With $\hat L_j^\star(t)=-\inf_{0\le u\le t}\hat X_j^\star(u)$ , $j=1,2$ , we know that $V_j^\star(t)=0$ , for every point of (right) increase of $\hat L_j^\star(t)$ . From this and the martingale of [Reference Kella and Whitt33], it may be concluded that the following is a zero-mean martingale:

(53) \begin{align}\hat{k}(s_1,s_2)\int_0^t e^{-s_1\hat V_1^\star(u)-s_2\hat V_2^\star(u)}\textrm{d}u &-e^{-s_1\hat V_1^\star(t)-s_2\hat V_2^\star(t)}+e^{-s_1\hat V_1^\star(0)-s_2\hat V_2^\star(0)}\nonumber\\[5pt] &-s_1\int_0^t e^{-s_2 \tilde{V}_2^\star(u)}\textrm{d}\hat L^\star_1(u)\nonumber\\[5pt] &-s_2\int_0^te^{-s_1 \tilde{V}_1^\star(u)}\textrm{d}\hat L_2^\star(u)\, .\end{align}

It has become standard by now (see, e.g., [Reference Kella31, Corollary 2.3]; it also follows from the theory of multivariate reflected Brownian motions on the nonnegative orthant) that if $\hat{v}^\star(s_1,s_2)$ is the LST of the stationary version of $\hat V^\star$ , then taking expectations in Equation (53) yields

(54) \begin{align}0 &=t\hat{k}(s_1,s_2)\hat{\nu}^\star(s_1,s_2)-\hat{\nu}^\star(s_1,s_2)+\hat{\nu}^\star(s_1,s_2)\nonumber \\[5pt] &\quad -s_1\mathbb{E}\int_0^t e^{-s_2 \tilde{V}_2^\star(u)}\textrm{d}\hat L^\star_1(u)-s_2\mathbb{E}\int_0^t e^{-s_1 \tilde{V}_1^\star(u)}\textrm{d}\hat L^\star_2(u),\end{align}

and in particular for $t=1$ we have

(55) \begin{equation}\hat{k}(s_1,s_2)\hat{\nu}^\star(s_1,s_2)=s_1\hat{f}_1(s_2)+s_2 \hat{f}_2(s_1)\,,\end{equation}

where

\begin{align*}\hat{f}_1(s_2) =\mathbb{E}\int_0^1 e^{-s_2 \tilde{V}_2^\star(u)}\textrm{d}\hat L^\star_1(u),\\[5pt] \hat{f}_2(s_1) =\mathbb{E}\int_0^1 e^{-s_1 \tilde{V}_1^\star(u)}\textrm{d}\hat L^\star_2(u)\,.\end{align*}

Our objective is to determine the unknown function in the left-hand side of (55): $\hat{\nu}^{\star}(s_1, s_2)$ .

The present setting is in several respects much more general: a two-dimensional Lévy input process, non-exponential visit periods, and asymmetry.

In the symmetric case, viz. for $\hat{\theta}_1=\hat{\theta}_2=\hat{\theta}$ , the key functional equation (55) reduces to

(56) \begin{equation}\frac{s_1+ s_2+\frac{1}{2 \hat \theta}(s_1-s_2)^2}{s_1s_2}\hat{\nu}^\star(s_1,s_2)=\frac{\hat{f}_1(s_2)}{s_2}+\frac{ \hat{f}_2(s_1)}{s_1},\end{equation}

which is in essence identical to (14) for $\hat{\theta}=4c/\mu$ . In this case, the starting point of the analysis matches, revealing that the results also match. It is important to note that, in the symmetric case, although the analysis is identical to the one performed in Section 4, the setting of this section is much broader than that of Section 4.

In the analysis that follows, we do not restrict ourselves to a symmetric system as we did in Section 4; instead we consider general $\theta_1$ , $\theta_2$ . For this general setting, we calculate the unknown function in the left-hand side of (55) using the BVM by applying Step A and Step B in an analogous manner as in Section 4. Unfortunately, several of the convenient simplifications that transpire in the symmetric case and that eventually lead to the elegant result of Theorem 2 are not allowed in the asymmetric case $\theta_1\neq\theta_2$ , as can be seen in the analysis that follows and in the result of Theorem 4.

Kernel analysis. To apply the BVM, one needs to investigate the zeros of the kernel $\hat{k}(s_1, s_2)$ . By setting $\hat{k}(s_1, s_2) = 0$ , we obtain

(57) \begin{align}\hat{s}_2^{\pm}(s_1) = s_1-\hat{\theta}_2\pm\sqrt{\hat{\theta}_2^2-2s_1( \hat{\theta}_1 + \hat{\theta}_2)}.\end{align}

Note that $\hat{s}_2^{\pm}(s_1)$ has a single branching point at $s_1 = { \hat{\theta}_2^2}/{2( \hat{\theta}_1+ \hat{\theta}_2)}$ . For real-valued $s_1$ with $s_1 > { \hat{\theta}_2^2}/{2( \hat{\theta}_1+ \hat{\theta}_2)}$ , the function $\hat{s}_2^{\pm}(s_1)$ is complex-valued. Letting $\hat{s}_2^{\pm}(s_1) = u \pm i v$ , we obtain, after straightforward computations, that

(58) \begin{align}v^2 = 2(\hat{\theta}_1+\hat{\theta}_2)\!\left(u+\frac{\hat{\theta}_2(2\hat{\theta}_1+\hat{\theta}_2)}{2(\hat{\theta}_1+\hat{\theta}_2)}\right)\!,\end{align}

which describes a parabola in the complex plane. We shall restrict ourselves to the following set:

\begin{align*}\hat{E}_1 = \left\{(u, v) \in (\!-\!{\hat{\theta}_2(2\hat{\theta}_1+\hat{\theta}_2)}/{2(\hat{\theta}_1+\hat{\theta}_2)}, \infty) \times \mathbb{R} \mid v^2 = 2(\hat{\theta}_1+\hat{\theta}_2)\!\left(u+\frac{\hat{\theta}_2(2\hat{\theta}_1+\hat{\theta}_2)}{2(\hat{\theta}_1+\hat{\theta}_2)}\right) \right\}\!.\end{align*}

This domain will allow us to determine $\hat{f}_1(\!\cdot\!)$ , while the symmetric domain obtained by considering the roots $\hat{s}_1^{\pm}(s_2)$ (which will result in a symmetric parabola with $\hat{\theta}_1$ and $\hat{\theta}_2$ interchanged) will allow us to determine $\hat{f}_2(\!\cdot\!)$ .

BVM: solution of the functional equation (55). Notice that, by definition, $\hat{f}_1(s_2)$ is analytic for $\textrm{Re}[s_2] \geq 0$ . It still remains to show that $\hat{f}_1(s_2)$ is analytic on the strip $-{\hat{\theta}_2(2\hat{\theta}_1+\hat{\theta}_2)}/{2(\hat{\theta}_1+\hat{\theta}_2)}< \textrm{Re}[s_2] < 0$ . We shall return to this point at a later stage; cf. Lemma 6.

Now we take $s_1$ with $s_1 > { \hat{\theta}_2^2}/{2( \hat{\theta}_1+ \hat{\theta}_2)}$ and $\hat{s}_2^{\pm}(s_1) = u \pm iv$ , with $(u, v) \in \hat{E}_1$ . For all such pairs $(s_1, \hat{s}_2^{\pm}(s_1))$ , the left-hand side of (55) becomes zero, and hence, for all $s_2 = \hat{s}_2^{\pm}(s_1)$ , we have

(59) \begin{equation}\frac{\hat{f}_1(s_2)}{s_2}=-\frac{ \hat{f}_2(s_1)}{s_1}.\end{equation}

For $s_1 > { \hat{\theta}_2^2}/{2( \hat{\theta}_1+ \hat{\theta}_2)}$ , the right-hand side of the above equation is real, thus yielding

(60) \begin{equation}\textrm{Re}\left[i\,\frac{\hat{f}_1(s_2)}{s_2} \right] = 0, \, \, \mbox{for} \, \, s_2 = \hat{s}_2^{\pm}(s_1) = u \pm iv, \;\mbox{with}\; (u, v) \in \hat{E}_1 \backslash \{(0, 0)\}.\end{equation}

We thus see that ${\hat{f}_1(s_2)}/{s_2} $ is analytic inside the contour $\hat{E}_1$ , say on $\hat{E}_1^+$ , except for $s_2=0$ , which is a pole in $\hat{E}_1^+$ . The above problem now reduces to a Riemann–Hilbert problem with a pole and with boundary $\hat{E}_1$ ; see [Reference Cohen and Boxma21, Section I.3.3]. To transform it into a (standard) Riemann–Hilbert problem on the unit circle D, we define $\hat\phi_1$ (with inverse $\hat\psi_1$ ) to be a conformal mapping of the interior of the unit circle D onto $\hat{E}_1^+$ with normalization conditions

\begin{align*} \hat{\phi}_1\!\left(\!-\!1\right) &=\infty, \\[5pt] \hat{\phi}_1(0) &=\frac{1 - \sqrt{2}\cos\!\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}{1 + \sqrt{2}\cos\!\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}, \quad \text{and} \\[5pt] \hat{\phi}_1\!\left(1\right) &=-\frac{\hat{\theta}_2(2\hat{\theta}_1+\hat{\theta}_2)}{2(\hat{\theta}_1+\hat{\theta}_2)}.\end{align*}

Following the same steps as in Section 4, leading to Theorem 2, we again translate the Riemann–Hilbert BVP on and inside $\hat{E}_1$ to the simple Riemann–Hilbert BVP with a pole. The solution of the BVP (60) is

(61) \begin{align}\hat{f}_1(s_2) &= \alpha_1 s_2- i \beta_1 \!\left(\hat{\psi}_1(s_2)-\frac{1 - \sqrt{2}\cos\!\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}{1 + \sqrt{2}\cos\!\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}\right)s_2 \nonumber\\[5pt] &\quad +i \frac{\bar{\beta} _1}{\hat{\psi}_1(s_2)-\frac{1 - \sqrt{2}\cos\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}{1 + \sqrt{2}\cos\left(\frac{\pi \hat{\theta}_1}{2(\hat{\theta}_1+\hat{\theta}_2)} \right)}} s_2, \qquad s_2 \in \hat{E}_1^+\cup\hat{E}_1\backslash\{0\},\end{align}

where $\hat{\psi}_1(\!\cdot\!)$ is the conformal mapping from the parabola $\hat{E}_1$ to the unit circle D given in the following lemma, and the constants $ \alpha_1$ and $ \beta_1$ , together with the full solution for the scaled buffer content processes, are given in Theorem 4.

Lemma 5. For $z \in \mathbb{C}$ and for $j=1,2$ , the conformal map

(62) \begin{equation}\hat{\psi}_j(z) =\frac{1 - \sqrt{2}\cosh\!\left(\frac{\pi}{2(\hat{\theta}_1+\hat{\theta}_2)} \sqrt{2(\hat{\theta}_1+\hat{\theta}_2)z-{\hat{\theta}_j^2}}\right)}{1 + \sqrt{2}\cosh\!\left(\frac{\pi}{2(\hat{\theta}_1+\hat{\theta}_2)} \sqrt{2(\hat{\theta}_1+\hat{\theta}_2)z-\hat{\theta}_j^2}\right)}\end{equation}

maps the interior of parabola

\begin{align*}v^2= 2(\hat{\theta}_1+\hat{\theta}_{2})\!\left(u+\frac{\hat{\theta}_{3-j}(2\hat{\theta}_j+\hat{\theta}_{3-j})}{2(\hat{\theta}_1+\hat{\theta}_2)}\right)\end{align*}

onto the interior of the unit circle D.

Now we are in position to state the main theorem of this section, in which we obtain an explicit expression for the scaled stationary buffer content process LST in heavy traffic.

Theorem 4. For $j=1,2$ , the scaled stationary buffer content process LST in heavy traffic is given by, for $\textrm{Re}[s_j] > -{\hat{\theta}_j(2\hat{\theta}_{3-j}+\hat{\theta}_j)}/{2(\hat{\theta}_1+\hat{\theta}_2)}$ ,

(63) \begin{align}\hat{f}_j(s_{3-j}) &=\mathbb{E}\left[\int_0^1 e^{-s_{3-j} \hat{V}_{3-j}^\star(u)}\textrm{d}\hat L^\star_j(u)\right]\nonumber\\[5pt] &= s_{3-j} \frac{\pi \sin \!\left(\frac{\pi \hat{\theta} _j}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right) }{\!\left(\sqrt{2} \sin \!\left(\frac{\pi \hat{\theta} _{3-j}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1\right)^2}\Bigg[-\frac{\cos \!\left(\frac{\pi \hat{\theta} _j}{ \hat{\theta} _1+ \hat{\theta} _2}\right)+4}{2 \sin \!\left(\frac{\pi \hat{\theta} _{3-j}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+\sqrt{2}}\nonumber\\[5pt] & \quad+2 \sqrt{2} \Bigg(\frac{1}{\sqrt{2} \sin \!\left(\frac{\pi \hat{\theta} _{3-j}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1}-\frac{1}{\sqrt{2} \cosh \!\left(\frac{\pi \sqrt{2\left( \hat{\theta} _1+ \hat{\theta} _2\right) s_{3-j} - \hat{\theta} _j^2}}{2\left( \hat{\theta} _1+\theta _2\right)}\right)+1}\Bigg)\nonumber\\[5pt] & \quad+\frac{1}{\sqrt{2} \Bigg(\frac{1}{\sqrt{2} \sin \!\left(\frac{\pi \hat{\theta} _{3-j}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1}-\frac{1}{\sqrt{2} \cosh \!\left(\frac{\pi \sqrt{2\left( \hat{\theta} _1+ \hat{\theta} _2\right) s_{3-j} - \hat{\theta} _j^2}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1}\Bigg)}-\sqrt{2}\Bigg].\end{align}

For $j=1,2$ , the scaled joint stationary buffer content process LST in heavy traffic is given by, for $\textrm{Re}[s_j] > -{\hat{\theta}_j(2\hat{\theta}_{3-j}+\hat{\theta}_j)}/{2(\hat{\theta}_1+\hat{\theta}_2)}$ ,

(64) \begin{align}\hat{\nu}^\star(s_1,s_2) &=\mathbb{E}\left[\int_0^1 e^{-s_1\hat V_1^\star(u)-s_2\hat V_2^\star(u)}\textrm{d}u\right]=\nonumber\\[5pt] &= \frac{s_1s_2}{\hat{k}(s_1, s_2)}\!\left( \frac{\hat{f}_1(s_2)}{s_2}+\frac{\hat{f}_2(s_1)}{s_1}\right),\end{align}

where $\hat{k}(s_1, s_2) = \hat \theta_1 s_1+\hat \theta_2 s_2+\frac{1}{2}(s_1-s_2)^2$ .

It is now convenient to formulate and prove the postponed Lemma 6.

Proof. Similarly to the proof of Lemma 3, we need to show that $\hat{f}_1(s_2)$ has no poles in $ -{\hat{\theta}_j(2\hat{\theta}_{3-j}+\hat{\theta}_j)}/{2(\hat{\theta}_1+\hat{\theta}_2)}<\textrm{Re}[s_j] <0$ . This is equivalent to considering the roots of the two denominators appearing in Equation (63), i.e., the zeros of

\begin{align*}1+ \sqrt{2}\cosh\!\left(\frac{\pi}{2(\hat{\theta}_1+\hat{\theta}_2)} \sqrt{2(\hat{\theta}_1+\hat{\theta}_2)s_{3-j}-\hat{\theta}_{j}^2}\right)\end{align*}

and the zeros of

\begin{align*} \frac{1}{\sqrt{2} \sin \!\left(\frac{\pi \hat{\theta} _{3-j}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1}-\frac{1}{\sqrt{2} \cosh \!\left(\frac{\pi \sqrt{2\left( \hat{\theta} _1+ \hat{\theta} _2\right) s_{3-j} - \hat{\theta} _j^2}}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)+1}.\end{align*}

For the former zeros, note that these are

(65) \begin{align}s_{3-j}=\frac{1}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right) }\!\left(\hat{\theta} _j^2-4(\hat{\theta} _1+ \hat{\theta} _2)^2({3}/{4}+2n)^2\right),\qquad n\in\mathbb{Z}.\end{align}

For the latter zeros, straightforward computations reveal that these are

(66) \begin{align}s_{3-j}=\frac{1}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right) }\!\left(\hat{\theta} _j^2-(\!-\hat{\theta}_j+4(\hat{\theta} _1+ \hat{\theta} _2)n)^2\right),\qquad n\in\mathbb{Z}\setminus\{0\},\end{align}

where we needed to exclude the case $s_{3-j}=0$ (which is equivalent to $n=0$ in the last expression), as this is not a pole for Equation (63). In both cases, (65) and (66), it is straightforward to show that $s_{3-j}<-{\hat{\theta}_j(2\hat{\theta}_{3-j}+\hat{\theta}_j)}/{2(\hat{\theta}_1+\hat{\theta}_2)}$ for all n.

Concluding this section, we would like to remark that in the case $\hat{\theta}_1=\hat{\theta}_2=\hat{\theta}$ , the result of Theorem 4 reduces exactly to that of Theorem 3 for $\hat{\theta}={4c}/{\mu}$ . This proves that the two limits (stationarity and heavy traffic) commute. Moreover, one can easily verify that, in the asymmetric case, taking the limit $\hat\theta_j\downarrow0$ while $\hat\theta_{3-j}>0$ yields

\begin{align*}\lim_{\hat\theta_j\downarrow0}\hat{f}_j(s_{3-j}) =\lim_{\hat\theta_j\downarrow0}\mathbb{E}\left[\int_0^1 e^{-s_{3-j} \hat{V}_{3-j}^\star(u)}\textrm{d}\hat L^\star_j(u)\right]=0,\end{align*}

as the $\sin \!\left(\frac{\pi \hat{\theta} _j}{2\left( \hat{\theta} _1+ \hat{\theta} _2\right)}\right)$ becomes zero and all other quantities are bounded.