Workload analysis of a two-queue fluid polling model

Abstract

In this paper, we analyze a two-queue random time-limited Markov-modulated polling model. In the first part of the paper, we investigate the fluid version: fluid arrives at the two queues as two independent flows with deterministic rate. There is a single server that serves both queues at constant speeds. The server spends an exponentially distributed amount of time in each queue. After the completion of such a visit time to one queue, the server instantly switches to the other queue, i.e., there is no switch-over time.

For this model, we first derive the Laplace–Stieltjes transform (LST) of the stationary marginal fluid content/workload at each queue. Subsequently, we derive a functional equation for the LST of the two-dimensional workload distribution that leads to a Riemann–Hilbert boundary value problem (BVP). After taking a heavy-traffic limit, and restricting ourselves to the symmetric case, the BVP simplifies and can be solved explicitly.

In the second part of the paper, allowing for more general (Lévy) input processes and server switching policies, we investigate the transient process limit of the joint workload in heavy traffic. Again solving a BVP, we determine the stationary distribution of the limiting process. We show that, in the symmetric case, this distribution coincides with our earlier solution of the BVP, implying that in this case the two limits (stationarity and heavy traffic) commute.

1. Introduction

The stationary analysis of multi-dimensional Markov processes associated with queueing models is often quite challenging. Even in the two-dimensional case, the characterization of the stationary distribution of fundamental queueing models (such as the shortest queue routing and the coupled processors [24, 21]) requires solving boundary value problems (BVPs). The intrinsic complexity of this analysis has led to the development of asymptotic techniques, studying the stationary distribution in some limiting regime of the model parameters; one prominent example is the heavy-traffic limit, first introduced by Kingman [35] for the single-server queue. In the heavy-traffic limit, a scaled version of the workload process is shown to have a nontrivial limit, which may serve as an approximation to the non-scaled process. The methodological contribution in this paper is to combine both approaches: for a specific fluid flow polling model with random time-limited service (which will be specified later), we first derive the BVP, which characterizes the stationary distribution, but for which no explicit solution is known. We then formulate the BVP obtained in the heavy-traffic limit, which in the symmetric case leads to an explicit solution for the two-dimensional stationary distribution (in heavy traffic). A second contribution of our paper is to investigate the heavy-traffic limit of a generalization of the polling model using process limits, allowing for Lévy input processes into the queues and a more general switching process for the server. Following [33], instead of directly focusing on the stationary distribution and deriving a functional equation for it, we characterize the entire scaled limit process as a two-dimensional reflected Brownian motion in the positive orthant. We show that in the earlier special (Markovian) symmetric case, the stationary distribution of the heavy-traffic process limit coincides with the heavy-traffic limit of the stationary distribution; thus the heavy-traffic limit and the time limit to stationarity commute.

Specifically, the model that we consider is a polling model with two queues and a single server that moves between the two queues to provide them with service. The policy that governs the switching is random time-limited (RTL): the duration of the service period at any queue is random, having an exponential distribution. All these service periods are independent, and the server always remains at a queue until the exponentially distributed time expires, even if that queue is empty and the other is not. The input and, when the served queue is not empty, the output processes for both queues are assumed to be deterministic fluid streams (with identical rates). Our motivation to study this RTL Markov-modulated fluid polling model comes from our earlier paper [40], in which the present fluid model emerged as an (asymptotic) approximation of a two-queue RTL polling model with Poisson arrivals and exponential service times. In the present paper, we show that even under the simplifying assumptions of fluid flows with constant inflow and outflow rates, and symmetric queues, determining the joint stationary workload distribution still requires solving a complicated BVP. In the heavy-traffic limit, we obtain and explicitly solve a BVP which is similar to that studied in [14] and belongs to a class of two-dimensional BVPs that is discussed in [21] (see also [24]). It is intuitive to recognize that in the asymmetric case with different loads on the two queues, the queue dynamics are easier to describe (compared to the symmetric case), as the workloads become independent in heavy traffic, reducing the analysis to that of the two marginals. The heaviest-loaded queue reaches the saturation point (and must be scaled), while the other queue remains stable (and needs not be scaled). For this reason, in the first part of the paper, we focus on the symmetric case: the two queues are entirely symmetric in terms of inflow and outflow rates, as well as the server visiting times. The symmetry assumption puts us in the most interesting case for the heavy-traffic setting that we consider in this paper: it ensures that the workloads in the two queues are of comparable magnitude in heavy traffic. In the second part of the paper, we extend the analysis both to a more general (Lévy input) model and to the study of the symmetric as well as the asymmetric case.

Related literature. Both fluid queueing models and polling models have received much attention in the literature on stochastic service systems; we refer to the surveys [37, 15] for overviews of the literature on fluid queues, and to the surveys [11, 12] for similar overviews on polling.

In contrast with the extensive literature on fluid queues and polling, there are very few studies focusing on polling systems with fluid input. Some exceptions are Czerniak and Yechiali [22], Boxma et al. [13], Remerova et al. [38], and Adan et al. [3]; see also [12, Section 6]. A recent heavy-traffic analysis of a fluid model with two queues in series is Koops et al. [36].

Polling models with time-limited service also have not been widely studied. Coffman, Fayolle and Mitrani [20] have analyzed a two-queue polling model with exponential visit periods; in their case (in contrast to the service protocol pertaining to the model studied in this paper), the server does not stay at an empty queue. They determine the probability generating function (PGF) of the joint stationary queue length distribution by solving a Riemann–Carleman BVP. In a series of papers, Al Hanbali et al. (see, e.g., [5]) consider a polling model with several queues and exponential visit periods. They relate the PGFs of the number of customers in a queue at the end of the server’s visit to that queue and at its beginning. This is used as input for a numerical scheme to approximate the joint queue length PGF at the instants of the server’s departure from the queues. Further references are provided in [39]; that paper and [40] also present a perturbation method for obtaining queue length PGFs in time-limited polling models.

Organization of the paper. In Section 2, we describe the RTL Markov-modulated fluid queue under consideration. In Section 3.1, we briefly present the Laplace–Stieltjes transforms (LSTs) of the model’s marginal stationary workload distributions and obtain their heavy-traffic limits. Section 3.2 is devoted to a discussion of the joint workload distribution analysis. In Section 4, restricting ourselves to the symmetric case, we derive an explicit expression for the LST of the joint stationary workload distribution in heavy traffic by solving a Riemann–Hilbert BVP. Several numerical experiments are performed in Section 5 in order to get more insight into the model. Section 6 is devoted to the computation of the scaled joint stationary distribution of an analogous model with a general Lévy input, generalizing the results obtained in Sections 3–4.

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

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 [32, Section 5], and also from the analysis performed in [17], 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 [39]), 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 [21]. 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 [21, 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. [21, Section I.3.3] and [14, 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. [14] and [21, 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}

Proof. The conformal mapping $\psi(z)$ is obtained by taking the composition of the following conformal mappings:

1. i. The conformal mapping $\eta(z) = z - \frac{c}{\mu}$ , where $z = x + i y$ , maps the parabola $y^2 = \frac{16c}{\mu} (x + \frac{3c}{\mu})$ onto the parabola $y^2 = \frac{16c}{\mu} (x + \frac{4c}{\mu})$ .

2. ii. From [9, p. 113], we have that the conformal mapping $\xi(z) = i \cosh\!\left(\frac{\pi}{4} \sqrt{\frac{\mu}{c}z}\right)$ maps the interior of the parabola $y^2 = \frac{16c}{\mu} (x + \frac{4c}{\mu})$ onto the interior of the upper half-plane $\textrm{Im}[\xi] > 0$ .

3. iii. As shown in [16, p. 326, Equation (6)], the conformal mapping $w(z) = \frac{1 + i \sqrt{2} z }{1 - i \sqrt{2} z}$ maps the upper half-plane (i.e., $\textrm{Im}[z] > 0$ ) onto the interior of the unit circle $|w| <1$ .

It follows from [9, Theorem III] that a composition of conformal mappings again is a conformal mapping. Hence the composition mapping $\psi(z) \;:\!=\; w(\xi(\eta(z)))$ conformally maps the interior of the parabola (20) onto the interior of the unit circle D.

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.

Lemma 3. The LST of the total scaled workload in heavy traffic, as given in (28), is analytic on the strip $-{3 c}/{\mu} < \textrm{Re}[s_2] < 0$ .

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 [16, 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. [8, 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.

Remark 2. To compute the limiting probability density function of the scaled total workload in heavy traffic, we need to invert the LST (28). The appearance of LSTs with a hyperbolic cosine and their probabilistic interpretation has a longstanding tradition in probability theory; see, e.g., [8] and the references therein. As we shall need these results for the proof of Lemma 4, we review them briefly 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. [8, 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. [19, 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. [8, 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 [8, Equation (3.11)].

As stated in [19, 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 [8, 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 [25, 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 [34, 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. [8, 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 [19]. 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 [30, 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 [19, 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$ .

Proof. By substituting $\tilde{\nu}^{\star}(s_j, s_j)$ , $j = 1, 2$ (obtained from the LST (28)), into Equation (14), we obtain $\tilde{\nu}^{\star}(s_1, s_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., [6]. 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}$ .

Remark 4. Notice from the simulation results in Table 1 that the correlation coefficient of the joint workload is not very sensitive to the traffic load.

Remark 5. The scaled two-dimensional workload LST $\tilde{\nu}^{\star}(s_1, s_2)$ can be inverted numerically; cf. [18, 23]. We have not been able to explicitly invert the LST. The scaled marginal distributions in heavy traffic are exponential (cf. Lemma 1), which suggests that the two-dimensional scaled workload distribution in heavy traffic might be a bivariate exponential distribution. It is observed in [10, Theorem 4.2] that the minimal correlation of any bivariate exponential distribution is $1 - {\pi^2}/{6} = -0.6449$ . This does not exclude the possibility that the joint workload in heavy traffic has a bivariate exponential distribution, as our correlation equals $-0.4203$ .

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 [1] 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 [26] and [42].

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 [41, 29]. This central limit version may also be concluded from [7, Theorem 3.2, p. 178]. The functional limit theorem may be concluded from, e.g., [27, 28]. 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 [33], 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., [31, 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 [21, 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.

Proof. The proof of the lemma is identical to that of Lemma 5, and as such it is omitted.

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

Proof. Setting $s_2=0$ yields on the one hand that the left-hand side of (61) is equal to $\hat{f}_{1}(0)=\hat{\theta}_1$ and on the other hand that the right-hand side of (61) is equal to $i \bar{{\beta}}_1 /\hat{\psi'}_{\!\!1}(0)$ . Substituting $\hat{\psi}_1(z)$ from Lemma 5, we obtain the value for ${\beta}_1$ . Moreover, since $\hat{f}_{1}(\infty)=0$ , we obtain the value for ${\alpha}_1$ . The same approach can also be used for the determination of $\hat{f}_2(s_1)$ . After tedious but straightforward computations, Equation (63) follows.

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

Lemma 6. For $j=1,2$ , the jth scaled stationary buffer content process LST in heavy traffic is analytic on the strip $ -{\hat{\theta}_j(2\hat{\theta}_{3-j}+\hat{\theta}_j)}/{2(\hat{\theta}_1+\hat{\theta}_2)}<\textrm{Re}[s_j] <0$ .

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.