The generalized join the shortest orbit queue system: stability, exact tail asymptotics and stationary approximations

Abstract

We introduce the generalized join the shortest queue model with retrials and two infinite capacity orbit queues. Three independent Poisson streams of jobs, namely a smart, and two dedicated streams, flow into a single-server system, which can hold at most one job. Arriving jobs that find the server occupied are routed to the orbits as follows: Blocked jobs from the smart stream are routed to the shortest orbit queue, and in case of a tie, they choose an orbit randomly. Blocked jobs from the dedicated streams are routed directly to their orbits. Orbiting jobs retry to connect with the server at different retrial rates, i.e., heterogeneous orbit queues. Applications of such a system are found in the modeling of wireless cooperative networks. We are interested in the asymptotic behavior of the stationary distribution of this model, provided that the system is stable. More precisely, we investigate the conditions under which the tail asymptotic of the minimum orbit queue length is exactly geometric. Moreover, we apply a heuristic asymptotic approach to obtain approximations of the steady-state joint orbit queue-length distribution. Useful numerical examples are presented and shown that the results obtained through the asymptotic analysis and the heuristic approach are agreed.

1. Introduction

In this work, we focus on the asymptotic stationary behavior of the generalized join the shortest orbit queue (GJSOQ) policy with retrials. This model is a natural generalization of the join the shortest orbit queue system, recently introduced in [11], by considering both non-identical retrial rates, and additional dedicated arrival streams that route jobs directly to the orbit queues if an arriving job finds the server busy; for more information about developments on the analysis of retrial queues, see the seminal books in [7,12].

We consider a single-server retrial system with two infinite capacity orbit queues accepting three independent arrival streams. The service station can handle at most one job. Arriving jobs are directed initially to the service station. An arriving job (of either stream) that finds the server idle starts service immediately. In case the server is busy upon a job's arrival, the blocked arriving job is routed to an orbit queue as follows: two of the arrival streams are dedicated to each orbit queue, i.e., an arriving job of stream $m$ that finds the server busy, it joins orbit queue $m$, $m=1,2$. An arriving job of the third stream, i.e., the smart stream, that finds the server busy, it joins the least loaded (i.e., the shorter) orbit queue, and in case of a tie, the smart job joins either orbit queue with probability $1/2$. Orbiting jobs retry independently to connect with the service station after a random time period (and according to the constant retrial policy) that depends on the type of the orbit queue. Note that the model at hand is described by a non-homogeneous Markov modulated two-dimensional random walk; see Figure 1. Our main concern is to investigate the asymptotic behavior of the GJSOQ system and to derive stationary approximations of its joint orbit queue-length distribution.

Figure 1. Matrix transition diagram.

1.1. Related work

The join-the-shortest-queue (JSQ) policy is widely used for load-balancing in stochastic networks. Yet the analytic derivation of the stationary distribution is known to be far from trivial.

The stationary behavior of the standard (i.e., non-modulated, without retrials) two-dimensional JSQ problem had initially studied in [18] and further investigated in [21], which was shown that the minimum queue length has exactly geometric asymptotics using the generating function approach. In [9,14], the authors presented a robust mathematical approach through generating functions and complex variables arguments to study its stationary behavior. However, this approach does not provide explicit expressions for the equilibrium distribution, and it is not so useful for numerical computations. The compensation method (CM), introduced in [1–3], provides an elegant and direct method to obtain explicitly the equilibrium joint queue-length distribution as an infinite series of product form terms, by solving directly the equilibrium equations.

Numerical/approximation methods were also applied: see the power series algorithm (PSA), e.g., [8], and the matrix geometric method; see e.g., [17,37], for which connections with CM were recently reported in [19]. PSA is numerically satisfactory for relatively lower-dimensional models, although, the theoretical foundation of this method is still incomplete. By expressing the equilibrium distribution as a power series in some variable based on the model parameters (usually the load), PSA transforms the balance equations into a recursively solvable set of equations by adding one dimension to the state space. For the (non-modulated) multidimensional JSQ model, the authors in [4] constructed upper and lower bounds for any performance measure based on two related systems that were easier to analyze. For a comparative analysis of the methods used for the analysis of multidimensional queueing models (including JSQ), see [5].

The stationary behavior of a two-queue system under the JSQ policy with Erlang arrivals was investigated in [6] by using the CM. The queueing model in [6] is described by a multilayer random walk in the quarter plane. Quite recently, in [11], the CM was applied to investigate the joint stationary distribution for a Markov-modulated random walk in the quarter plane, which describes a standard (i.e., without dedicated traffic streams) symmetric join the shortest orbit queue system with retrials. For such a model, it was also shown that the tail decay rate of the minimum orbit queue has an exactly geometric asymptotic.

Since the exact solutions discussed above are extremely complicated, it is useful to evaluate these expressions in certain limiting cases, in order to gain more insight into the qualitative structure of the particular model. Asymptotic formulas often clearly show the dependence of the solutions on the various variables/parameters in the problem, whereas the full exact expressions may be difficult to interpret in terms of the underlying model. Clearly, an asymptotic formula can never contain as much quantitative (numerical) information as an exact answer, but it provides reasonably accurate numerical results at a greatly reduced computational cost.

The tail asymptotic behavior for the two-queue case has been extensively studied. For the standard JSQ, i.e., without retrials, the problem for the case of homogeneous servers was answered in [21], while for the case of heterogeneous servers in [40]. The behavior of the standard generalized JSQ (GJSQ) problem was investigated by using a random walk structure in [15], and by using a quasi birth-death (QBD) formulation in [26]; see also [25] which extends Malyshev's approach [27] to the GJSQ paradigm. However, those two papers have not completely solved the tail asymptotic problem, since they focus on the so called strongly pooled condition. By using the Markov additive approach and the optimization technique developed in [29], the author in [30] completely characterized the weak decay rates in terms of the transition probabilities and provided a complete solution. The decay rate of stationary probabilities was also analyzed in heavy traffic and via large deviations in [41] and [42], respectively; see also [36]. We further mention [22], where the authors obtain heuristics from the balance equations.

In [39], the authors studied the tail decay problem of JSQ system with two queues and PH-type arrivals; see also in [38] for the case of more than two queues and threshold jockeying. These works investigate the tail behavior of non-homogeneous multilayered random walks. We further mention the recent work in [33], which focused on a thorough investigation of the tail behavior of space homogeneous two-dimensional skip-free Markov modulated reflecting random walk; see also [31]. Using a QBD and a reflecting random walk formulation, the authors in [23] showed that the tail asymptotic of the minimum queue length in a Markovian system of $k\geq 3$ parallel queues is exactly geometric.

1.2. Contribution

This work considers for the first time the generalized join the shortest queue policy with retrials. We use this model as a vehicle to investigate the asymptotic behavior of non-homogeneous Markov-modulated two-dimensional random walks, and with a particular interest on a modulation that allows a completely tractable analysis. We first investigate the stability conditions. Then, we focus on the tail asymptotic behavior of the stationary distribution and show that under specific conditions, the tail asymptotic of the minimum orbit queue length is exactly geometric. Finally, by using directly the equilibrium equations, we provide an asymptotic approach to obtain approximations to the steady-state joint orbit queue-length distribution for both states of the server. The approximation results agreed with the results obtained by the asymptotic analysis.

1.2.1. Fundamental contribution

In [11], we provided exact expressions for the stationary distribution of the standard (i.e., with smart arrival stream and no dedicated arrival streams), symmetric (i.e., identical retrial rates) JSOQ system by using the CM. To our best knowledge, the work in [11] is the only work in the related literature that provided analytic results regarding the stationary behavior of the JSQ policy with retrials. It is well known from the related literature on the standard (i.e., without retrials) JSQ model with two queues that under the presence of the dedicated streams, the CM collapses. This is because, under the presence of dedicated traffic, the corresponding two-dimensional random walk allows transitions to the North, and to the East. A similar situation arises also in our case, where in the corresponding Markov-modulated two-dimensional random walk, we have a similar behavior; see Figures 1 and 7. Thus, this work provides the only analytic results for the generalized JSQ problem with retrials available in the literature.

Note that our system is described by a non-homogeneous two-dimensional random walk modulated by a two-state Markov process. In such a case, the phase process represents the state of the server and affects the evolution of the level process, i.e., the orbit queue lengths in two ways: (i) The rates at which certain transitions in the level process occur depend on the state of the phase process; see Figure 2. Thus, a change in the phase might not immediately trigger a transition of the level process but changes its dynamics (indirect interaction). (ii) A phase change does trigger an immediate transition of the level process (direct interaction).

Figure 2. An instance of state transitions from states belonging to angle $r_{1}$ (left), and ray $d$ (right), given the state of the server.

For such a queueing system, we show that the tail asymptotic of the minimum orbit queue length for fixed values of their difference and server's state is exactly geometric. To accomplish this task, we transform the original process to a Markov-modulated random walk in the half plane. Then, we focused on the censored Markov chain referring to the busy states, which is now a standard two-dimensional random walk in the half plane and studied the tail asymptotic behavior. Using a relation among idle and busy states, we also study the tail behavior for the idle states. To our best knowledge, there is no other analytic study on the tail asymptotics properties of Markov-modulated random walks in the half plane.

Moreover, we provide a simple heuristic approach to approximate the stationary distribution, valid when one of the orbit queue lengths is large, by distinguishing the analysis between the symmetric and asymmetric cases (this is because the asymmetric case reveals additional technical difficulties compared with the symmetric case; see subsections 5.2 and 5.3). Our derived theoretical results through the heuristic approach agreed with those derived by the tail asymptotic analysis. Moreover, the advantage of our approach is that we cope directly with the equilibrium equations, without using generating functions or diffusion approximations.

Stability conditions are also investigated. In particular, having in mind that the orbit queues grow only when the server is busy, we focused on the censored chain at busy states and provide necessary and sufficient conditions for the ergodicity of the censored chain by using Foster–Lyapunov arguments. We conjecture that these conditions are also necessary and sufficient for the original process. Simulation experiments indicate that our conjecture is true but a more formal justification is needed. We postponed the formal proof in a future work.

1.2.2. Application-oriented contribution

Besides its theoretical interest, the queueing model at hand has interesting applications in the modeling of wireless relay-assisted cooperative networks. Such systems operate as follows: There is a finite number of source users that transmit packets (i.e., the arrival streams) to a common destination node (i.e., the single service station), and a finite number of relay nodes (i.e., the orbit queues) that assist the source users by retransmitting their blocked packets, i.e., the packets that cannot access upon arrival of the destination node; see e.g., [10,35]. We may consider here a smart source user (possibly the use of the highest importance that transmits priority packets) that employs the JSQ protocol (i.e., the cooperation strategy among the smart source and the relays, under which, the smart user chooses to forward its blocked packet to the least loaded relay node), and two dedicated source users (that transmit packets of lower priority or packets that can only be handled by the specific relay). This work serves as a major step towards the analysis of even general retrial models operating under the JSQ policy.

The rest of the paper is organized as follows. In Section 2, we present the mathematical model in detail. Stability conditions are investigated in Section 3. The main results of this work are presented in Section 4. More precisely, in subsection 4.1, we present the exact geometric behavior in the minimum orbit queue-length direction under specific conditions (see Theorem 4.1 and Corollary 1), while in subsection 4.2, we present explicit asymptotic expressions for approximating the stationary distribution of the model; see Lemmas 2 and 3. The proofs of the main results are presented in Section 5. Our theoretical findings are validated through detailed numerical results that are presented in Section 6.

2. System model

We consider a single-server retrial system with two infinite capacity orbit queues. Three independent Poisson streams of jobs, say $S_{m}$, $m=0,1,2$, flow into the single-server service system, namely the smart stream $S_{0}$, and the dedicated streams $S_{m}$, $m=1,2$. The arrival rate of stream $S_{m}$ is $\lambda _{m}$, $m=0,1,2,$, and let $\lambda :=\lambda _{0}+\lambda _{1}+\lambda _{2}$. The service system can hold at most one job. The required service time of each job is independent of its type and is exponentially distributed with rate $\mu$. If an arriving type-$m$, $m=1,2,$ job finds the (main) server busy, it is routed to a dedicated retrial (orbit) queue that operates as an $\cdot /M/1/\infty$ queue. An arriving type-$0$ job (i.e., the smart job) that finds the server busy, it is routed to the shortest orbit queue (i.e., the least loaded orbit), while in case of a tie, it is routed to either orbit with probability 1/2. Orbiting jobs try to access the server according to a constant retrial policy (i.e., orbits behave as $\cdot /M/1/\infty$ queues). In particular, the orbit queue $m$ attempts to retransmit a job (if any) to the main service system at a Poisson rate of $\alpha _{m}$, $m=1,2$.

Let $N_{m}(t)$ be the number of jobs in orbit queue $m$, $m=1,2$, at time $t$, and $C(t)$ be the state of the server, i.e., $C(t)=1$, when the server is busy, and $C(t)=0$ when it is idle at time $t$, respectively. Then, $X(t)=\{(N_{1}(t),N_{2}(t),C(t));t\geq 0\}$ is an irreducible Markov process with state space $S=\mathbb {N}_{0}\times \mathbb {N}_{0}\times \{0,1\}$, where $\mathbb {N}_{0}=\{0,1,2,\ldots \}$. Let $X=\{(N_{1},N_{2},C)\}$ the stationary version of $\{X(t);t\geq 0\}$, and define the stationary probabilities

$$p_{i,j}(k)=\lim_{t\to\infty}\mathbb{P}((N_{1}(t),N_{2}(t),C(t)) =(i,j,k))=\mathbb{P}((N_{1},N_{2},C)=(i,j,k)),$$

for $(i,j,k)\in S$. The equilibrium equations are

(2.1)\begin{equation} (\lambda+\alpha_{1}1_{\{i>0\}}+\alpha_{2}1_{\{j>0\}})p_{i,j}(0)=\mu p_{i,j}(1),\quad i,j\geq 0,\end{equation}

(2.2)\begin{align} & (\lambda+\mu)p_{i,j}(1)=[\lambda_{0}H(i-j+1)+\lambda_{2}]p_{i,j-1}(1) +[\lambda_{0}H(j-i+1)+\lambda_{1}]p_{i-1,j}(1)\nonumber\\ & \quad + \lambda p_{i,j}(0)+\alpha_{1}p_{i+1,j}(0)+\alpha_{2}p_{i,j+1}(0),\quad i,j\geq 1,\end{align}

(2.3)\begin{equation}(\lambda+\mu)p_{i,0}(1)=\lambda p_{i,0}(0)+\lambda_{1}p_{i-1,0}(1)+\alpha_{1}p_{i+1,0}(0)+\alpha_{2}p_{i,1}(0),\quad i\geq 2, \end{equation}

(2.4)\begin{align} & (\lambda+\mu)p_{0,j}(1)=\lambda p_{0,j}(0)+\lambda_{2}p_{0,j-1}(1)+\alpha_{1}p_{1,j}(0)+\alpha_{2}p_{0,j+1}(0),\quad j\geq 2, \end{align}

(2.5)\begin{align} & (\lambda+\mu)p_{1,0}(1)=\lambda p_{1,0}(0)+\alpha_{1}p_{2,0}(0)+\alpha_{2}p_{1,1}(0)+\left(\lambda_{1}+\frac{\lambda_{0}}{2}\right)p_{0,0}(1), \end{align}

(2.6)\begin{align} & (\lambda+\mu)p_{0,1}(1)=\lambda p_{0,1}(0)+\alpha_{1}p_{1,1}(0)+\alpha_{2}p_{0,2}(0)+\left(\lambda_{2}+\frac{\lambda_{0}}{2}\right)p_{0,0}(1), \end{align}

(2.7)\begin{align} & (\lambda+\mu)p_{0,0}(1)=\lambda p_{0,0}(0)+\alpha_{1}p_{1,0}(0)+\alpha_{2}p_{0,1}(0), \end{align}

with the normalization condition $\sum _{i=0}^{\infty }\sum _{j=0}^{\infty }(p_{i,j}(0)+p_{i,j}(1))=1$, and where

$$H(n)=\left\{\begin{array}{ll} 1, & n\geq 1,\\ \dfrac{1}{2}, & n=0,\\ 0, & n\leq{-}1. \end{array}\right.$$

For reasons that will become clear in the following sections, we consider the corresponding uniformized discrete time Markov chain through the uniformization technique. Since the total transition rate from each state is bounded by $\lambda +\mu +\alpha _{1}+\alpha _{2}$, we can construct by uniformization a discrete time Markov chain with the same stationary distribution as that of $\{X(t);t\geq 0\}$. Without loss of generality, let the uniformization parameter $\theta =\lambda +\mu +\alpha _{1}+\alpha _{2}= 1$. The uniformized Markov chain $X(n)=\{(N_{1,n},N_{2,n},C_{n})\}$ of $\{X(t);t\geq 0\}$ has six regions of spatial homogeneity: two angles, say $r_{1}=\{i>j> 0,k=0,1\}$ and $r_{2}=\{j>i> 0,k=0,1\}$, three rays, say $h=\{j=0,i>0,k=0,1\}$, $v=\{i=0,j>0,k=0,1\}$, $d=\{i=j>0,k=0,1\}$ and the points $\{(0,0,k), k=0,1\}$. Then, the matrix transition diagram partitioned according to the state of the server, is depicted in Figure 1, where

\begin{align*} A_{1,0}^{(r_1)}& =\begin{pmatrix} 0 & 0\\ 0 & \lambda_{1} \end{pmatrix}=A_{1,0}^{(h)}=A_{1,0}^{(0)},\quad A_{1,0}^{(r_2)}=\begin{pmatrix} 0 & 0\\ 0 & \lambda_{1}+\lambda_{0} \end{pmatrix}=A_{1,0}^{(v)},\quad A_{1,0}^{(d)}=\begin{pmatrix} 0 & 0\\ 0 & \lambda_{1}+\lambda_{0}/2 \end{pmatrix},\\ A_{0,1}^{(r_2)}& =\begin{pmatrix} 0 & 0\\ 0 & \lambda_{2} \end{pmatrix}=A_{0,1}^{(v)}=A_{0,1}^{(0)},\quad A_{0,1}^{(r_1)}=\begin{pmatrix} 0 & 0\\ 0 & \lambda_{0}+\lambda_{2} \end{pmatrix}=A_{0,1}^{(h)},\quad A_{0,1}^{(d)}=\begin{pmatrix} 0 & 0\\ 0 & \lambda_{0}/2+\lambda_{2} \end{pmatrix},\\ A_{0,-1}^{(r_2)}& =\begin{pmatrix} 0 & \alpha_{2}\\ 0 & 0 \end{pmatrix}=A_{0,-1}^{(v)}=A_{0,-1}^{(d)}=A_{0,-1}^{(r_1)},\quad A_{0,0}^{(v)}=\begin{pmatrix} \mu+\alpha_{1} & \lambda\\ \mu & \alpha_{1}+\alpha_{2} \end{pmatrix},\\ A_{{-}1,0}^{(r_1)}& =\begin{pmatrix} 0 & \alpha_{1}\\ 0 & 0 \end{pmatrix}=A_{{-}1,0}^{(h)}=A_{{-}1,0}^{(d)}=A_{{-}1,0}^{(r_2)},\quad A_{0,0}^{(h)}=\begin{pmatrix} \mu+\alpha_{2} & \lambda\\ \mu & \alpha_{1}+\alpha_{2} \end{pmatrix},\\ A_{0,0}^{(r_1)}& =\begin{pmatrix} \mu & \lambda\\ \mu & \alpha_{1}+\alpha_{2} \end{pmatrix}=A_{0,0}^{(r_2)}=A_{0,0}^{(d)},\quad A_{0,0}^{(0)}=\begin{pmatrix} \mu+\alpha_{1}+\alpha_{2} & \lambda\\ \mu & \alpha_{1}+\alpha_{2} \end{pmatrix}. \end{align*}

An illustration of the transitions from a state belonging to the angle $r_{1}$ (e.g., from the state $(4,2,k)$, $k=0,1$), and from a state belonging to the ray $d$ (e.g., from the state $(3,3,k)$, $k=0,1$), is given in Figure 2, left and right, respectively.

Lets say some more words regarding the derivation of $A_{i,j}^{(q)}$, $-1\leq i$, $j\leq 1$, $q=r_{1},r_{2},h,v,0$. For example, $A_{1,0}^{(r_1)}$ contains the transition probabilities from a state that belongs to the angle $r_{1}$ and results in an increase by one at the orbit queue 1, i.e., from $(i,j,k)$ to $(i+1,j,l)$, where $i>j$, $k,l=0,1$. Clearly, such a transition occurs only when we have a dedicated arrival of type 1, and the server is busy. Similarly, given a state in the ray $d$, $A_{1,0}^{(d)}$ contains transition probabilities that result in an increase by one at orbit queue 1, given that both orbits have the same occupancy. Such a transition is done either with the arrival of a dedicated job or with the arrival of a smart job who sees both orbits with the same number of jobs and is routed with probability $1/2$ to orbit queue 1. The rest of $A_{i,j}^{(q)}$, $-1\leq i,j\leq 1$, $q=r_{1},r_{2},h,v,0$, are constructed similarly.

In the following sections, we provide our main results that refer to the asymptotic behavior of $\{X(t);t\geq 0\}$. We first investigate the stability conditions. Then, our first main result refers to the investigation of the decay rate of the tail probabilities for the shortest queue length in a steady state; see Theorem 4.1. We cope with this task in a number of steps, after considering the corresponding discrete time Markov chain $\{X(n);n\geq 0\}$ through the uniformization technique:

1. 1. We transform the uniformized Markov chain to create the minimum and the difference of the queue states process (i.e., the Markov modulated random walk in the half plane), and then consider the censored process at busy states.

2. 2. We investigate the stationary tail decay rate of the censored process at busy states.

3. 3. Using a relation among busy and idle states, we show that the asymptotic properties of the shortest queue for the idle states are the same with its asymptotic properties for the busy states.

Our second main result relies on a heuristic approach to construct approximations for the stationary distribution of the original process. In particular, our aim (see Lemmas 2 and 3) is to solve (2.1)–(2.6) as either $i$ or $j\to \infty$. This task is accomplished by constructing a solution of (2.1)–(2.6) separately in each of the regions defined as follows: Region I: $i\gg 1$, $j\gg 1$, Region II: $i=0$ or $1$, $j\gg 1$, Region III: $i\gg 1$, $j=0$ or $1$.

3. On the stability condition

In this section, we investigate the stability condition of the model at hand. Stability condition for the standard generalized join the shortest queue (GJSQ) model without retrials was recently investigated; see e.g., [15,24]. Clearly, the presence of retrials along with the presence of the dedicated arrival flows to each orbit when the server is busy, complicates the problem considerably. Note that our model is described by a non-homogeneous Markov-modulated two-dimensional nearest-neighbor random walk, and its stability condition, to our best knowledge, is still an open problem. We mention here that the stability condition of a homogeneous Markov-modulated two-dimensional nearest-neighbor random walk was recently investigated in [34] by using the concept of induced Markov chains.

Here on, by noting that the orbit queues grow only when the server is busy, we construct a new discrete time Markov chain embedded at epochs in which the server is busy, i.e., the censored Markov chain at busy states (censored Markov chains have been widely used for proving the uniqueness of the stationary vector for a recurrent countable-state Markov chain [20,43]). Then, using standard Foster–Lyapunov arguments, we provide its stability conditions and having in mind that orbits grow only when the server is busy (it is natural to assume that the behavior of the original process at the busy states heavily affects its convergence), we conjecture that the obtained stability conditions for the censored Markov chain on the busy states coincide with those of the original one. Simulation experiments indicate that our conjecture is true. The formal justification of our conjecture is an interesting open problem and we let it be a future study.

We first consider the uniformized discrete time Markov chain, say $X(n)=\{(N_{1,n},N_{2,n},C_{n});n\geq 0\}$, of $\{X(t);t\geq 0\}$ with transition diagram given in (1) and state space $S=E\cup E^{c}$, where $E=\mathbb {N}_{0}\times \mathbb {N}_{0}\times \{1\}$, $E^{c}=\mathbb {N}_{0}\times \mathbb {N}_{0}\times \{0\}$. Then, we partition the transition matrix $P$ of $X(n)$ according to $E$, $E^{c}$ into:

where

\begin{align*} P_{0,0}& ={\rm diag}(A_{0},A_{1},A_{1},\ldots),\quad P_{1,0}={\rm diag}(L_{0},L_{0},L_{0},\ldots), \\ P_{0,1}& =\begin{pmatrix} B_{0} & & & & \\ B_{1} & B_{0} & & & \\ & B_{1} & B_{0} & & \\ & & \ddots & \ddots & \end{pmatrix},\quad P_{1,1}=\begin{pmatrix} D_{0,0} & D_{0,1} & & & \\ & D_{1,1} & D_{1,2} & & \\ & & D_{2,2} & D_{2,3} & \\ & & & \ddots & \ddots \end{pmatrix},\\ A_{0}& ={\rm diag}(\mu+\alpha_{1}+\alpha_{2},\mu+\alpha_{1},\mu+\alpha_{1},\ldots),\quad A_{1}={\rm diag}(\mu+\alpha_{2},\mu,\mu,\ldots),\\ L_{0}& =\mu I_{\infty},\quad B_{1}=\alpha_{1}I_{\infty},\\ B_{0}& =\begin{pmatrix} \lambda & & & & \\ \alpha_{2} & \lambda & & & \\ & \alpha_{2} & \lambda & & \\ & & \ddots & \ddots & \end{pmatrix},\quad D_{i,i}=(\alpha_{1}+\alpha_{2})I_{\infty}+Q_{i+1,i+2},\quad i=0,1,\ldots, \end{align*}

where $I_{\infty }$ is the identity matrix of infinite dimension and $Q_{i+1,i+2}$ is a infinite dimension matrix with only non-zero entries at the superdiagonal and such that the $(k,k+1)-$entry equal to $\lambda _{0}+\lambda _{2}$, $k=1,2,\ldots,i$, the $(i+1,i+2)-$entry equals ${\lambda _{0}}/{2}+\lambda _{2}$, and the $(i+k,i+k+1)-$entry equals $\lambda _{2}$, $k=2,3,\ldots$. Finally, the matrix $D_{i,i+1}$, $i=0,1,\ldots$ is a diagonal matrix of infinite dimension with $(k,k)-$entry equal to $\lambda _{1}$, $k=1,2,\ldots,i$, the $(i+1,i+1)-$entry equals ${\lambda _{0}}/{2}+\lambda _{1}$, and the $(i+k,i+k)-$entry equals $\lambda _{0}+\lambda _{1}$, $k=2,3,\ldots$. In what follows, denote $\hat {\lambda }:=\lambda (\lambda +\alpha _{1}+\alpha _{2})$, $\hat {\lambda }_{k}:=\lambda _{k}(\lambda +\alpha _{1}+\alpha _{2})$, $k=0,1,2,$ and $\hat {\mu }_{k}=\mu \alpha _{k}$, $k=1,2$.

Since $P_{0,0}$ is a diagonal matrix, its fundamental matrix, say $\tilde {P}_{0,0}$ has the form

\begin{align*} \tilde{P}_{0,0}& =\sum_{n=0}^{\infty}P_{0,0}^{n}={\rm diag}(\tilde{A}_{0},\tilde{A}_{1},\tilde{A}_{1},\ldots),\\ \tilde{A}_{0}& ={\rm diag}\left(\frac{1}{\lambda},\frac{1}{\lambda+\alpha_{2}},\frac{1}{\lambda+\alpha_{2}},\ldots\right)\\ \tilde{A}_{1}&={\rm diag}\left(\frac{1}{\lambda+\alpha_{1}},\frac{1}{\lambda+\alpha_{1}+\alpha_{2}},\frac{1}{\lambda+\alpha_{1}+\alpha_{2}},\ldots\right). \end{align*}

The censored chain $\{\tilde {X}(n);n\geq 0\}$ at busy states has six regions of spatial homogeneity: two angles, say $r_{1}=\{i>j> 0\}$ and $r_{2}=\{j>i> 0\}$, three rays, say $h=\{j=0,i>0\}$, $v=\{i=0,j>0\}$, $d=\{i=j>0\}$ and the point $(0,0)$. Then, the one-step transition probability matrix of the censored chain $\{\tilde {X}(n);n\geq 0\}$, given by $P^{(E)}=P_{1,1}+P_{1,0}\tilde {P}_{0,0}P_{0,1}$ is as follows:

• • In region $i>j>0$ (labeled as $r_{1}$),

  \begin{align*} p_{1,0}^{(r_{1})}& =\lambda_{1},\quad p_{0,1}^{(r_{1})}=\lambda_{0}+\lambda_{2},\quad p_{{-}1,0}^{(r_{1})}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\\ p_{0,-1}^{(r_{1})}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\quad p_{0,0}^{(r_{1})}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}}, \end{align*}

• • In region $j>i>0$ (labeled as $r_{2}$),

  \begin{align*} p_{1,0}^{(r_{2})}& =\lambda_{0}+\lambda_{1},\quad p_{0,1}^{(r_{2})}=\lambda_{2}, \quad p_{{-}1,0}^{(r_{2})}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\\ p_{0,-1}^{(r_{2})}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\quad p_{0,0}^{(r_{2})}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}}, \end{align*}

• • In region $i=j>0$ (labeled as $d$),

  \begin{align*} p_{1,0}^{(d)}& =\frac{\lambda_{0}}{2}+\lambda_{1},\quad p_{0,1}^{(d)}=\frac{\lambda_{0}}{2}+\lambda_{2},\quad p_{{-}1,0}^{(d)}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\\ p_{0,-1}^{(d)}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\quad p_{0,0}^{(d)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}}, \end{align*}

• • In region $i=0$, $j>0$ (labeled as $v$),

  $$p_{1,0}^{(v)}=\lambda_{0}+\lambda_{1},\quad p_{0,1}^{(v)}=\lambda_{2},\quad p_{0,-1}^{(v)}=\frac{\hat{\mu}_{2}}{\lambda+\alpha_{2}},\quad p_{0,0}^{(v)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{2}},$$

• • In region $j=0$, $i>0$ (labeled as $h$),

  $$p_{1,0}^{(h)}=\lambda_{1},\quad p_{0,1}^{(h)}=\lambda_{0}+\lambda_{2},\quad p_{{-}1,0}^{(h)}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}},\quad p_{0,0}^{(h)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}},$$

• • For $i=j=0$ (labeled as $O$),

  $$p_{1,0}^{(O)}=\lambda_{1}+\frac{\lambda_{0}}{2},\quad p_{0,1}^{(O)}={+}\frac{\lambda_{0}}{2}+\lambda_{2},\quad p_{0,0}^{(O)}=\alpha_{1}+\alpha_{2}+\mu.$$

Remark 1 Note that the censored chain $\{\tilde {X}(n);n\geq 0\}$ describes a new discrete time queueing system consisting of two parallel coupled queues with three job arrival streams, where one of them joins the shortest queue, and the other two are dedicated to each queue. The coupling feature of the two queues is easily realized by noting that, e.g., $p_{-1,0}^{(h)}={\hat {\mu }_{2}}/{(\lambda +\alpha _{1})}>{\hat {\mu }_{2}}/{(\lambda +\alpha _{1}+\alpha _{2})}=p_{-1,0}^{(r_{2})}$. The combination of the JSQ feature, along with the coupled processors feature considerably complicates the analysis.

Lemma 1 The censored chain $\{\tilde {X}(n);n\geq 0\}$ is positive recurrent if and only if one of the following conditions hold:

1. 1. $\rho _{1}:={\hat {\lambda }_{1}}/{\hat {\mu }_{1}}<1$, $\rho _{2}:={\hat {\lambda }_{2}}/{\hat {\mu }_{2}}<1$, $\rho ={\hat {\lambda }}/{(\hat {\mu }_{1}+\hat {\mu }_{2})}<1$,

2. 2. $\rho _{1}\geq 1$, $f_{1}:={\lambda (\lambda _{1}+\alpha _{1})}/{\hat {\mu }_{1}}-1<0$,

3. 3. $\rho _{2}\geq 1$, $f_{2}:={\lambda (\lambda _{2}+\alpha _{2})}/{\hat {\mu }_{2}}-1<0$.

Proof. The proof is given in Appendix A.

We now study the dynamics of the orbits in the original process $\{X(t);t\geq 0\}$ based on the stability criteria of the censored process $\{\tilde {X}(n);n\geq 0\}$ given in Lemma 1.

Figure 3 refers to the stability criteria 1. In Figure 3 (left), we have $\rho =0.78>\max \{\rho _{1}=0.433,\rho _{2}=0.325\}$ and $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|=19.084>0$. It is seen that in such a case, the original process is stable, and that both orbits are well balanced. This case corresponds to the strongly pooled case mentioned in [15] Theorem 2, i.e., the proportion of jobs that is routed to the least loaded orbit queue in case the server is busy is large. In Figure 3 (middle), we have the case $1>\rho >\max \{\rho _{1},\rho _{2}\}$, with $\rho _{1}>\rho _{2}$ but now $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|<0$. The original process is still stable but we can observe that the orbit lengths are not so close any more. This means that the stream that joins the shortest orbit failed to keep the orbit queue lengths close enough to each other. In case $\rho >1$, both orbits are unstable, even though $\rho _{1}<1$, $\rho _{2}<1$ (Figure 3 (right)).

Figure 3. Orbit dynamics for case $1>\rho >\max \{\rho _{1},\rho _{2}\}$, $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|>0$ (left), $1>\rho >\max \{\rho _{1},\rho _{2}\}$, $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|<0$ (middle), and for the case $\rho >1$, $\rho _{1}<1$, $\rho _{2}<1$ (right).

In Figure 4 (left), we have $\rho _{1}=0.428>\max \{\rho =0.4042,\rho _{2}=0.2675\}$, and $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|=-2.393<0$. In such a case, the system is still stable, but the smart stream failed to keep the orbit queue length very close. Similar observations can be also deduced from Figure 4 (right) for which $\rho _{1}=0.0944>\max \{\rho =0.061,\rho _{2}=0.021\}$, and $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|=0.7663>0$.

Figure 4. Orbit dynamics for case $1>\rho _{1}>\max \{\rho,\rho _{2}\}$, $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|<0$ (left), $1>\rho _{1}>\max \{\rho,\rho _{2}\}$, $\hat {\lambda }_{0}-|\hat {\lambda }_{1}-\hat {\lambda }_{2}+\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})|<0$ (right).

In Figure 5 (left), although $\rho _{1}>1$ ($\rho <1$, $\rho _{2}<1$), the fact that $f_{1}<0$ keeps the network stable. On the other hand when $f_{1}>0$ (Figure 5, middle), orbit queue 1 becomes unstable, and so it is the network. Note that in such a case orbit queue 2 remains stable. In Figure 5 (right), when $\rho >1$ both orbits becomes unstable, even though $\rho _{2}<1$. Similar observations can be deduced from Figure 6 that refers to the stability criteria 3.

Figure 5. Orbit dynamics for case $\rho <1$, $\rho _{1}>1$, $\rho _{2}<1$, $f_{1}<0$ (left), for the case $\rho <1$, $\rho _{1}>1$, $\rho _{2}<1$, $f_{1}>0$ (middle), and for the case $\rho >1$, $\rho _{1}>1$, $\rho _{2}<1$, $f_{1}>0$ (right).

Figure 6. Orbit dynamics for case $\rho <1$, $\rho _{1}<1$, $\rho _{2}>1$, $f_{2}<0$ (left), for the case $\rho <1$, $\rho _{1}<1$, $\rho _{2}>1$, $f_{2}>0$ (right).

Therefore, simulations experiments indicate that the stability criteria given in Lemma 1 constitute stability criteria for the GJSOQ system. A formal justification of this result is postponed in a future work.

4. Main results

4.1. Geometric decay in the minimum direction

In the following, we investigated the tail asymptotic of the stationary distribution of our model provided that $\rho _{1}<1$, $\rho _{2}<1$, $\rho <1$. We show that the tail asymptotic for the minimum orbit queue lengths for a fixed value of their difference, and server's state is exactly geometric. Moreover, we show that its decay rate is the square of the decay rate, i.e., $\rho ^{2}$, of the corresponding modified single-orbit queue, say the reference system, which is formally described at the end of this section.

Clearly, when $\lambda _{1}=\lambda _{2}=0$, both orbit queues are well balanced by the smart stream, and we again expect that the decay rate of the minimum of orbit queues equals $\rho ^{2}$. Note that this result is motivated by the standard Markovian JSQ model (without retrials), in which the tail decay rate of the minimum of queue lengths is the square of the tail decay rate of the M/M/2. Theorem 4.1 states that our model has a similar behavior.

Tail asymptotic properties for random walks in the half plane are still very limited including a recent paper by [30], in which the method developed in [29] was extended for studying the standard generalized join the-shortest-queue model. Our goal here is to cope with the tail asymptotic properties of a Markov modulated random walk in the half plane, with particular interest on a modulation that results in a tractable analysis.

Our approach is summarized in the following steps:

1. Step 1: We focus on the uniformized discrete time Markov chain $\{X(n);n\geq 0\}$, and perform the transformation $Z_{1,n}=\min \{N_{1,n},N_{2,n}\}$, $Z_{2,n}=N_{2,n}-N_{1,n}$.

2. Step 2: The discrete time Markov chain $Z(n)=\{(Z_{1,n},Z_{2,n},C_{n});n\geq 0\}$ is a Markov modulated random walk in the half plane with transition diagram as shown in Figure 7, where the matrices $A_{i,j}^{(m)}$, $i,j=0,\pm 1$, $m=r_{1},r_{2},d,v,h,0$ are given in Section 2. Then, we consider the censored Markov chain of $\{Z(n);n\geq 0\}$ on the busy states, which can be expressed explicitly. This censored Markov chain is a random walk in the half plane.

3. Step 3: We investigate the tail asymptotic behavior of this censored Markov chain at the busy states following [26], and show that the minimum orbit queue length decays geometrically; see Theorem 4.1.

4. Step 4: By exploiting a relation between idle and busy states (see (4.5)), we further show that the minimum orbit queue length at the idle states decays also geometrically; see Corollary 1.

Figure 7. The matrix transition diagram of $\{Z(n);n\geq 0\}$.

Note that the presence of dedicated streams that are routed to the shortest orbit in case of a busy server complicates considerably the problem of decay rates. We are also interested in when the two orbit queues are balanced. Intuitively, the two orbit queues being balanced when the smart arrival stream, which routes the job to the least loaded orbit queue, keeps the two orbit queue lengths close. Since the difference of the two orbit queues is the background state, the balancing of the two orbit queues is characterized by how the background process behaves as the level process (i.e., the shortest orbit queue) goes up.

By applying a similar procedure as in Section 3, the censored Markov chain on the busy states of $\{Z(n);n\geq 0\}$ is a two-dimensional random walk on the half plane $\{(m,l)\in \mathbb {Z}^{2};l\geq 0\}$, say $Y(n)=\{(Y_{1,n},Y_{2,n},1);n\geq 0\}$, where $Y_{1,n}=N_{2,n}-N_{1,n}$, $Y_{2,n}=\min \{N_{1,n},N_{2,n}\}$. Its transition rate diagram is given in Figure 8, where its one-step transition probabilities are:

(4.1)\begin{equation} \begin{aligned} q_{0,-1}^{(-)} & =\lambda_{1},\quad q_{{-}1,-1}^{(-)}=\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{0,1}^{(-)}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{1,1}^{(-)}=\lambda_{0}+\lambda_{2},\\ q_{0,0}^{(-)} & =\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{0,1}^{(+)}=\lambda_{2},\quad q_{0,-1}^{(+)}=\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\\ q_{{-}1,1}^{(+)} & =\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{1,-1}^{(+)}=\lambda_{0}+\lambda_{1},\quad q_{0,0}^{(+)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}},\\ q_{0,1}^{(2)} & =\lambda_{2}+\frac{\lambda_{0}}{2},\quad q_{{-}1,-1}^{(2)}=\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{{-}1,1}^{(2)}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}},\\ q_{0,-1}^{(2)} & =\frac{\lambda_{0}}{2}+\lambda_{1},\quad q_{0,0}^{(2)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}},\quad q_{0,-1}^{(1-)}=\lambda_{1},\\ q_{0,1}^{(1-)} & =\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}},\quad q_{1,1}^{(1-)}=\lambda_{0}+\lambda_{2},\quad q_{0,0}^{(1-)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}},\quad q_{0,1}^{(1+)}=\lambda_{2},\\ q_{0,-1}^{(1+)} & =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{2}},\quad q_{1,-1}^{(1+)}=\lambda_{0}+\lambda_{1},\quad q_{0,0}^{(1+)}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{2}},\\ q_{0,1}^{(0)} & =\lambda_{2}+\frac{\lambda_{0}}{2},\quad q_{0,0}^{(0)}=\alpha_{1}+\alpha_{2}+\mu,\quad q_{0,-1}^{(0)}=\lambda_{1}+\frac{\lambda_{0}}{2}. \end{aligned} \end{equation}

Remark 2 Note that $\{Y(n);n\geq 0\}$ can be also obtained by the censored Markov chain $\{\tilde {X}(n);n\geq 0\}$ (on the busy states of $\{X(n);n\geq 0\}$), by applying the transformation $Y_{1,n}=\min \{N_{1,n},N_{2,n}\}$, $Y_{2,n}=N_{2,n}-N_{1,n}$.

Figure 8. Transition rate diagram of the censored chain $\{Y(n);n\geq 0\}$.

The main result is summarized in the following theorem. Let $\boldsymbol {\pi }\equiv \{\boldsymbol {\pi }_{m},m=0,1,\ldots \}$ be the stationary joint distribution of $\min \{N_{1,n},N_{2,n}\}$ and $N_{2,n}-N_{1,n}$ when the server is busy, where $\boldsymbol {\pi }_{m}=\{\pi _{m,l}(1),l=0,\pm 1,\pm 2,\ldots \}$ is the subvector for level $m$, i.e.,

$$\pi_{m,l}(1)=\lim_{n\to\infty}\mathbb{P}(\min\{N_{1,n},N_{2,n}\}=m,N_{2,n}-N_{1,n}=l,C_{n}=1),\quad m\geq 0,l\in\mathbb{Z}.$$

Theorem 4.1 For the generalized join the shortest orbit queue system satisfying $\max \{\rho _{1},\rho _{2}\}<\rho <1$, and $\hat {\lambda }_{0}>|\hat {\lambda }_{2}-\hat {\lambda }_{1}+\rho ^{2}(\hat {\mu }_{1}-\hat {\mu }_{2})|$, the stationary probability vector $\boldsymbol {\pi }_{m}$ decays geometrically as $m\to \infty$ with decay rate $\rho ^{2}$, i.e.,

(4.2)\begin{equation} \lim_{m\to\infty}\rho^{{-}2m}\pi_{m,l}(1)=x_{l}(\rho^{{-}2}),\end{equation}

where

(4.3)\begin{equation} x_{l}(\rho^{{-}2})=\left\{\begin{array}{ll} \dfrac{\gamma_{2}+\hat{\lambda}_{0}/2}{\gamma_{2}}x_{0}\left(\rho\dfrac{\gamma_{2}}{\gamma_{1}+\hat{\lambda}_{0}}\right)^{|l|}, & l\leq{-}1,\\ x_{0}, & l=0,\\ \dfrac{\gamma_{1}+\hat{\lambda}_{0}/2}{\gamma_{1}}x_{0}\left(\rho\dfrac{\gamma_{1}}{\gamma_{2}+\hat{\lambda}_{0}}\right)^{l}, & l\geq 1, \end{array}\right. \end{equation}

where $x_{0}$ is a constant and $\gamma _{1}=\hat {\mu }_{1}\rho ^{2}+\hat {\lambda }_{2}$, $\gamma _{2}=\hat {\mu }_{2}\rho ^{2}+\hat {\lambda }_{1}$. Moreover, the two orbit queues are strongly balanced if and only if $\gamma _{2}<\rho (\gamma _{1}+\hat {\lambda }_{0})$, $\gamma _{1}<\rho (\gamma _{2}+\hat {\lambda }_{0})$.

Proof. The proof is based on a series of results and is given in Section 5.1.

Having known the exact tail asymptotic properties for $C(t) = 1$ (i.e., when the server is busy), we can investigate the tail asymptotic properties for $C(t) = 0$ (i.e., when the server is idle) based on the relationship given in (4.5) (which is similar to (2.1)).

Corollary 1 Based on the relationship given in (4.5),

(4.4)\begin{equation} \rho^{{-}2m}\pi_{m,l}(0)\sim \frac{\mu}{\lambda+\alpha_{1}+\alpha_{2}}x_{l}(\rho^{{-}2}),\quad l\in\mathbb{Z},\end{equation}

where $x_{l}(\rho ^{-2})$, $l\in \mathbb {Z}$ as given in Theorem 4.1.

Proof. We consider the Markov chain $Z(n)=\{(Z_{1,n},Z_{2,n},C_{n});n\geq 0\}$ where $Z_{1,n}=\min \{N_{1,n},N_{2,n}\}$, $Z_{2,n}=N_{2,n}-N_{1,n}$; see Figure 7. Knowing the asymptotic properties of $\pi _{m,l}(1)$, we will investigate those of $\pi _{m,l}(0)=\lim _{n\to \infty }\mathbb {P}(\min \{N_{1,n},N_{2,n}=m\},N_{2,n}-N_{1,n}=l,C_{n}=0)$. Note that the equilibrium equations that deal with the idle states are for $m>0$,

(4.5)\begin{equation} \pi_{m,l}(0)(\lambda+\alpha_{1}+\alpha_{2})=\mu\pi_{m,l}(1),\quad l\in\mathbb{Z},\end{equation}

which obviously leads to (4.4).

Therefore, we shown that under the conditions $\max \{\rho _{1},\rho _{2}\}<\rho <1$, $|\gamma _{1}-\gamma _{2}|<\hat {\lambda }_{0}$, the tail decay rate for the shortest orbit queue in our GJSOQ system equals $\rho ^{2}$. It is easy to see that this is the square of the tail decay rate for the orbit queue length of the reference system mentioned at the beginning of this section. The reference system operates as follows: Jobs arrive according to a Poisson process with rate $\lambda =\lambda _{0}+\lambda _{1}+\lambda _{2}$, and if they find the server available, they occupy it and get served after an exponentially distributed time with rate $\mu$. Otherwise they enter an infinite capacity orbit queue. Orbiting jobs retry to access the server according to the following rule: If there is only one job in orbit it retries after exponentially distributed time with rate $\alpha _{1}$ (note that this is not a restriction and whatever the retrial rate is, in this case, does not affect the final result). If there are more than one orbiting jobs, the first job in orbit queue retries with rate $\alpha _{1}+\alpha _{2}$ (This situation is similar to the case where there is a basic server in orbit queue, which transmit jobs to the service station at a rate $\alpha _{1}$, and in case there are more than one orbiting jobs, an additional server, which transmits jobs at a rate $\alpha _{2}$ helps the former one). Details on the stationary orbit queue-length distribution of the reference system as well as its decay rate are given in Appendix B.

Remark 3 Another point of interest refers to the determination of the decay rate of the marginal distribution, since Theorem 4.1 refers to the decay rate of the stationary minimum orbit queue length distribution for a fixed value of the orbit queue difference and server state. Clearly, the states of the server can be aggregated, since they are finite. However, the aggregation on the difference of queue sizes is not a trivial task and extra conditions maybe needed. One may expect that

(4.6)\begin{equation} \lim_{m\to\infty}\rho^{{-}2m}\mathbb{P}(\min\{N_{1},N_{2}\}=m) =\frac{1}{1-\mu}\sum_{l\in\mathbb{Z}}x_{l}(\rho^{{-}2}), \end{equation}

where $\{(N_{1},N_{2})\}$, the stationary version of $\{(N_{1,n},N_{2,n});n\geq 0\}$. This can be obtained using Theorem 4.1 and summing (4.2), (4.4) over the difference of the two orbit queues. Note here that the sum in the right-hand side of (4.6) definitely converges under the strongly balanced condition (see Theorem 4.1), since in such a case ${\rho \gamma _{2}}/{(\hat {\lambda }_{0}+\gamma _{1})}<\rho ^{2}<1$, ${\rho \gamma _{1}}/{(\hat {\lambda }_{0}+\gamma _{2})}<\rho ^{2}<1$. More importantly, under the conditions of Theorem 4.1, we always have ${\rho \gamma _{1}}/{(\gamma _{2}+\hat {\lambda }_{0})}<1$, ${\rho \gamma _{1}}/{(\gamma _{2}+\hat {\lambda }_{0})}<1$. Indeed, under conditions of Theorem 1, $\rho <1$, and $-\hat {\lambda }_{0}<\gamma _{1}-\gamma _{2}<\hat {\lambda }_{0}$. Thus, ${\rho \gamma _{1}}/{(\gamma _{2}+\hat {\lambda }_{0})}<1\Leftrightarrow \rho \gamma _{1}-\gamma _{2}<\hat {\lambda }_{0}$, which is true since $\rho \gamma _{1}-\gamma _{2}<\gamma _{1}-\gamma _{2}<\hat {\lambda }_{0}$. Similarly, ${\rho \gamma _{2}}/{(\gamma _{1}+\hat {\lambda }_{0})}<1\Leftrightarrow \rho \gamma _{2}-\gamma _{1}<\hat {\lambda }_{0}$, which is true since $\rho \gamma _{2}-\gamma _{1}<\gamma _{2}-\gamma _{1}<\hat {\lambda }_{0}$. The major problem to be resolved comes from the fact that since the difference of the two orbit queues is unbounded, in order to obtain the left-hand side of (4.6), it requires to verify the exchange of the summation and the limit. It seems that following the lines in [30] Theorem 1.5 by using a different framework based on [29], and focusing on the weak decay rates, we can show that the marginal distribution of the $\min \{N_{1},N_{2}\}$, under the same assumptions as given in Theorem 4.1, also decays with rate $\rho ^{2}$. Moreover, for the standard GJSQ, the authors in [26] also focused on the decay rate for the marginal probabilities of the $\min \{Q_{1},Q_{2}\}$; see Corollary 3.2.1, Theorem 2.2.1 and Corollary 2.2.1 in [26]. Following their lines, under the conditions of Theorem 4.1, let

$$\lim_{l\to\infty}\frac{\pi_{0,l,1}}{x_{l}(\rho^{{-}2})}=s,$$

where $\pi _{0,l,1}$, $l\geq 0$ is the stationary probability of an empty orbit queue for the censored Markov chain $\{Y(n);n\geq 0\}$, and $x_{l}(\rho ^{-2})$ as given in (4.3). Following [26] Corollary 2.2.1, if $0< s<\infty$, the marginal distribution of $\min \{Q_{1},Q_{2}\}$ for the busy states decays geometrically with rate $\rho ^{2}$. Now, using the relation in (4.5), the marginal distribution of $\min \{Q_{1},Q_{2}\}$ for the idle states decays also geometrically with the same rate. Thus, (4.6) has been proved. Therefore, the boundedness of $\pi _{0,l,1}/x_{l}(\rho ^{-2})$ is very crucial in proving that the marginal probabilities of $\min \{Q_{1},Q_{2}\}$ decay geometrically with rate $\rho ^{2}$. We postpone this verification in a future work.

Remark 4 In [39], the authors proved that the tail decay rate of the stationary minimum queue length in the standard symmetric JSQ system fed by PH arrivals (no dedicated traffic) is the square of the decay rate for the queue length in the standard PH/M/2. Such a system was described by a Markov-modulated two-dimensional random walk in the quarter plane. They proved this result based on matrix analytic approach and standard results on the decay rates. We believe that it is possible to extend their approach to cope with Markov-modulated two-dimensional random walk in the half plane, and get the same result as in Theorem 4.1.

4.2. A heuristic approach for stationary approximations

Our aim here is to develop a scheme to obtain approximations for the joint orbit queue-length distribution for each server state, valid when one of the queue lengths is large. The main results are summarized in Lemmas 2 and 3, where we distinguish the analysis to the symmetric and the asymmetric case, respectively. With the term symmetric, we mean that $\lambda _{1}=\lambda _{2}\equiv \lambda _{+}$, $\alpha _{1}=\alpha _{2}\equiv \alpha$.

We have to note that although the proofs of Lemmas 2 and 3 (see subsections 5.2 and 5.3, respectively) are similar, the asymmetric case reveals additional difficulties in the sense that we need further conditions to ensure that the approximated probabilities are well defined.

The main result for the symmetric case is the following.

Lemma 2 For $\hat {\lambda }<2\hat {\mu }$,

(4.7)\begin{equation} \begin{aligned} & p_{i,j}(1)\sim c\left\{\begin{array}{ll} \left(\dfrac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})}\right)^{i} \left(\dfrac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu})}{2\hat{\mu}(\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu})}\right)^{j} & \\ \quad\times \left[1+A_{-}\left(\dfrac{4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+})(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})^{2}}{\hat{\lambda}^{2}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))^{2}}\right)^{i}\right], & i< j,\ j\gg 1,\\ \left(\dfrac{\hat{\lambda}}{2\hat{\mu}}\right)^{i+j} \left(\dfrac{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}{\hat{\lambda}(\hat{\lambda}+2\hat{\mu})}\right), & i=j\gg 1,\\ \left(\dfrac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})}\right)^{j} \left(\dfrac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu})}{2\hat{\mu}(\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu}}\right)^{i} & \\ \quad\times \left[1+A_{-}\left(\dfrac{4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+})(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})^{2}}{\hat{\lambda}^{2}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))^{2}}\right)^{j}\right], & j< i,\ i\gg 1, \end{array}\right.\\ & p_{i,j}(0)=\frac{\mu}{\lambda+\alpha (1_{\{i\geq 1\}}+1_{\{j\geq 1\}})}p_{i,j}(1), \end{aligned} \end{equation}

where $c$ be the normalization constant, and

\begin{align*} A_{-}& =\frac{\lambda(\lambda+\alpha)+\hat{\mu}(1-x_{-}(1+\frac{\lambda+\alpha}{\lambda+2\alpha}))+\frac{\lambda_{+}(\lambda+\alpha)}{x_{+}}}{\hat{\mu}(\frac{\lambda+\alpha}{\lambda+2\alpha}x_{-}+x_{+}-1)-\lambda(\lambda+\alpha)-\frac{\lambda_{+}(\lambda+\alpha)}{x_{+}}},\\ x_{+}& =\frac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})},\quad x_{-}=\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\hat{\mu}x_{+}}. \end{align*}

Proof. The proof is given in Section 5.2.

The main result for the asymmetric case is the following.

Lemma 3 When, $\rho <1$, and $\hat {\lambda }_{0}>|\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$,

(4.8)\begin{equation} \begin{aligned} & p_{i,j}(1)\sim c\left\{\begin{array}{ll} \sqrt{\dfrac{\epsilon_{-}}{\epsilon_{+}}} \left(\dfrac{\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{i} \left(\dfrac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}\right)^{j} & \\ \quad \times \left[1+A_{-}\left(\dfrac{x_{-}}{x_{+}}\right)^{i}\right], & j>i,\ j\gg 1,\\ \rho^{i+j} \left[\dfrac{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\hat{\lambda}(1+\rho)} \sqrt{\dfrac{\epsilon_{+}}{\epsilon_{-}}}+ \dfrac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}(1+\rho)} \sqrt{\dfrac{\epsilon_{-}}{\epsilon_{+}}}\right], & i=j\gg 1,\\ \sqrt{\dfrac{\epsilon_{+}}{\epsilon_{-}}}\left(\dfrac{\rho(\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}\right)^{i} \left(\dfrac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{j} & \\ \quad \times \left[1+B_{-}\left(\dfrac{y_{-}}{y_{+}}\right)^{j}\right], & i>j,\ i\gg 1, \end{array}\right.\\ & p_{i,j}(0)=\frac{\mu}{\lambda+\alpha_{1}1_{\{i>0\}}+\alpha_{2}1_{j>0}}p_{i,j}(1), \end{aligned} \end{equation}

where $c$ be the normalization constant, $\epsilon _{\pm }$ are given in (5.38), and

\begin{align*} x_{+}& =\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}},\quad x_{-}=\frac{(\hat{\lambda}_{0}+\hat{\lambda}_{1})(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\mu}_{1}\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})},\\ y_{-}& =\frac{(\hat{\lambda}_{0}+\hat{\lambda}_{2})(\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\mu}_{2}\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})},\quad y_{+}=\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}} \end{align*}

\begin{align*} A_{-}& =\frac{x_{+}^{2}(\hat{\mu}_{1}\frac{\lambda+\alpha_{2}}{\lambda+\alpha_{1}+\alpha_{2}} +\frac{\lambda_{2}(\lambda+\alpha_{2})}{\rho^{2}})-x_{+}(\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2}) +\rho^{2}\hat{\mu}_{2}}{x_{+}(\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2})-[(\lambda_{0}+\lambda_{1})(\lambda+\alpha_{2})+\rho^{2}\hat{\mu}_{2}] -\frac{\lambda_{2}(\lambda+\alpha_{2})}{\rho^{2}}x_{+}^{2}},\\ B_{-}& =\frac{y_{+}^{2}(\hat{\mu}_{2}\frac{\lambda+\alpha_{1}}{\lambda+\alpha_{1}+\alpha_{2}} +\frac{\lambda_{1}(\lambda+\alpha_{1})}{\rho^{2}})-y_{+}(\lambda(\lambda+\alpha_{1})+\hat{\mu}_{1}) +\rho^{2}\hat{\mu}_{1}}{y_{+}(\lambda(\lambda+\alpha_{1})+\hat{\mu}_{1}) -[(\lambda_{0}+\lambda_{2})(\lambda+\alpha_{1})+\rho^{2}\hat{\mu}_{1}] -\frac{\lambda_{1}(\lambda+\alpha_{1})}{\rho^{2}}y_{+}^{2}}. \end{align*}

Proof. The proof is given in Section 5.3.

Remark 5 Note that in proving Lemma 3, we need to have the condition $\hat {\lambda }_{0}>|\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$. This condition was also needed in Theorem 4.1 and corresponds to the so called strongly pooled case in [15]. This condition refers to the case where the proportion of smart jobs is large compared with the dedicated traffic and results in balancing the orbit queue lengths. In such a case, the orbit queues overload in tandem; see Figure 9 (left). More importantly, by technical point of view, this condition is crucial in proving that the approximated join orbit queue-length probabilities are well defined as given in Lemma 3; see subsection 5.3.

Figure 9. Orbit dynamics when $\rho <1$, $\alpha _{1}>\alpha _{2}$, and $\hat {\lambda }_{0}> |\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$ (left), and when $\hat {\lambda }_{0}< |\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$, and $\lambda _{1}>\lambda _{0}>\lambda _{2}$ (right).

In case $\hat {\lambda }_{0}\leq |\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$, i.e., the weakly pooled case, there is a spectrum of possible drift directions, depending on how large is the proportion of the dedicated traffic. Moreover, the approximated joint orbit queue-length distribution cannot any more be written in such an elegant form and is not well defined. An illustration of both scenarios is given in Figure 9, where we performed some simulations. Note that in the strongly pooled case, orbit queue lengths grow in tandem as expected; Figure 9(left). In Figure 9(right), $\hat {\lambda }_{0}< |\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$, and $\lambda _{1}>\lambda _{0}>\lambda _{2}$, $\alpha _{1}>\alpha _{2}$, we observe that $N_{2}(t)$ lags behind $N_{1}(t)$.

5. Proofs of main results

To improve the readability of the paper we group in this section the proofs of the main results presented in section 4.

5.1. Proof of Theorem 4.1

The proof follows the lines in [26], and a series of results are needed to be proven. Note that $\{Y(n);n\geq 0\}$ is a discrete time QBD process with the following block structured transition matrix

$$P_{Y}=\begin{pmatrix} L_{0} & A_{1} & & & \\ A_{{-}1} & A_{0} & A_{1} & & \\ & A_{{-}1} & A_{0} & A_{1} & \\ & & A_{{-}1} & A_{0} & A_{1}\\ & & & \ddots & \ddots & \ddots \end{pmatrix},$$

where the infinite dimensional matrices $L$, $A_{i}$, $i=0,\pm 1$ are

\begin{align*} (L)_{l,l^{\prime}}& =\mathbb{P}(Y(n+1)=(0,l^{\prime},1)|Y(n)=(0,l,1)),\quad l,l^{\prime}\in\mathbb{Z},\\ (A_{i})_{l,l^{\prime}}& =\mathbb{P}(Y(n+1)=(m+i,l^{\prime},1)\,|\,Y(n)=(m,l,1)),\quad l,l^{\prime}\in\mathbb{Z},\ m\geq 1, \end{align*}

with their elements (i.e., the one-step transition probabilities) are given in (4.1). Under the stability condition, the stationary distribution exists. Denote it by the row vector $\boldsymbol {\pi }=(\boldsymbol {\pi _{0}},\boldsymbol {\pi _{1}},\ldots )$. It is well known that the stationary distribution has the following matrix geometric form:

(5.1)\begin{equation} \boldsymbol{\pi_{m}}=\boldsymbol{\pi_{0}}R^{m},\quad m\geq 1,\end{equation}

where $R$ is the minimal nonnegative solution of the equation

$$R=A_{1}+RA_{0}+R^{2}A_{{-}1}.$$

Since the size of $R$ is infinite, conditions for geometric tail decay rate of (5.1) are given in [26]:

Proposition 1 (Theorem 2.1 in [26])

Let $A\equiv A_{-1}+A_{0}+A_{1}$ is irreducible and aperiodic, and assume that the Markov additive process generated by $\{A_{i},i=0,\pm 1\}$ is 1-arithmetic. If there exists $z>1$ and positive vectors $\mathbf {x}$, $\mathbf {y}$ such that

(5.2)\begin{align} & \mathbf{x}=\mathbf{x}A_{*}(z),\quad \mathbf{y}=A_{*}(z)\mathbf{y}, \end{align}

(5.3)\begin{align} & \mathbf{x}\mathbf{y}<\infty, \end{align}

where $A_{*}(w)=w^{-1}A_{-1}+A_{0}+wA_{1}$, $w\neq 0$. Then $R$ has left and right eigenvectors $\mathbf {x}$ and $\mathbf {r}\equiv (I-A_{0}-RA_{-1}-z^{-1}A_{-1})\mathbf {y}$, respectively, with eigenvalue $z^{-1}$. Moreover, if

(5.4)\begin{equation} \boldsymbol{\pi}_{0}\mathbf{y}<\infty,\end{equation}

then

(5.5)\begin{equation} \lim_{m\to\infty}z^{n}\boldsymbol{\pi}_{m}=\frac{\boldsymbol{\pi}_{0}\mathbf{r}}{\mathbf{x}\mathbf{r}}\mathbf{x}.\end{equation}

Remark 6 It is easy to see that for our queueing model, both the aperiodicity and the irreducibility of $A$ are easily verified. Moreover, the 1-arithmetic property of the Markov additive process generated by $\{A_{i};i=0,\pm 1\}$ is readily satisfied by its transition structure.

For our case, the matrix $A_{*}(z)$ is given by

where

\begin{align*} a_{{-}1}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}}z^{{-}1}+\lambda_{1},\quad a_{0}=\alpha_{1}+\alpha_{2}+\frac{\lambda\mu}{\lambda+\alpha_{1}+\alpha_{2}},\quad a_{1}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}}+(\lambda_{0}+\lambda_{2})z,\\ c_{{-}1}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}}z^{{-}1}+\lambda_{1}+\frac{\lambda_{0}}{2},\quad c_{0}=a_{0},\quad c_{1}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}}z^{{-}1} +\frac{\lambda_{0}}{2}+\lambda_{2},\\ b_{{-}1}& =\frac{\hat{\mu}_{2}}{\lambda+\alpha_{1}+\alpha_{2}}+(\lambda_{0}+\lambda_{1})z,\quad b_{0}=a_{0},\quad b_{1}=\frac{\hat{\mu}_{1}}{\lambda+\alpha_{1}+\alpha_{2}}z^{{-}1}+\lambda_{2}. \end{align*}

Denote the left and the right invariant vectors of $A_{*}(z)$, $z>1$ by

$$\mathbf{x}(z)=(\ldots,x_{{-}1}(z),x_{0},x_{1}(z),\ldots),\quad \mathbf{y}(z)=(\ldots,y_{{-}1}(z),y_{0},y_{1}(z),\ldots)^{{\rm T}}.$$

We are left to examine conditions (5.2), (5.3), (5.4). To this end, let

\begin{align*} \beta_{i,j}(z)& =(\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})^{2}-4(\hat{\mu}_{i} +(\hat{\lambda}_{0}+\hat{\lambda}_{j})z)(\hat{\mu}_{j}z^{{-}1}+\hat{\lambda}_{i}),\quad i,j=1,2,\\ f(z)& =\frac{\hat{\mu}_{2}z^{{-}1}+\hat{\lambda}_{1}+ \frac{\hat{\lambda}_{0}}{2}}{2(\hat{\mu}_{2}z^{{-}1}+\hat{\lambda}_{1})} \left(1-\frac{\sqrt{\beta_{1,2}(z)}}{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}}\right) +\frac{\hat{\mu}_{1}z^{{-}1}+\hat{\lambda}_{2} +\frac{\hat{\lambda}_{0}}{2}}{2(\hat{\mu}_{1}z^{{-}1}+\hat{\lambda}_{2})} \left(1-\frac{\sqrt{\beta_{2,1}(z)}}{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}}\right). \end{align*}

Lemma 4 For the censored chain $\{Y(n);n\geq 0\}$ on the busy states of our GJSOQ, the conditions (5.2) and (5.3) hold if and only if there exists a $z>1$ such that

(5.6)\begin{equation} \beta_{1,2}(z)>0,\quad \beta_{2,1}(z)>0,\end{equation}

and $f(z)=1$. In such a case, the left and right invariant vectors $\mathbf {x}(z)$ and $\mathbf {y}(z)$ of $A_{*}(z)$ satisfying $\mathbf {x}(z)\mathbf {y}(z)<\infty$ have the following forms:

(5.7)\begin{equation} \begin{aligned} x_{l}(z) & =\left\{\begin{array}{ll} \left(1+\dfrac{\hat{\lambda}_{0}}{2(\hat{\mu}_{2}z^{{-}1}+\hat{\lambda}_{1})}\right)x_{0}\eta_{{\rm min}}^{|l|}, & l\leq{-}1,\\ \left(1+\dfrac{\hat{\lambda}_{0}}{2(\hat{\mu}_{1}z^{{-}1}+\hat{\lambda}_{2})}\right)x_{0}\theta_{{\rm min}}^{l},\quad l\geq 1, \end{array}\right.,\\ y_{l}(z) & =\left\{\begin{array}{ll} y_{0}\eta_{{\rm max}}^{-|l|}, & l\leq{-}1,\\ y_{0}\theta_{{\rm max}}^{{-}l}, & l\geq 1, \end{array}\right. \end{aligned} \end{equation}

where $\eta _{{\rm min}}<\eta _{{\rm max}}$ are the roots of the equation:

(5.8)\begin{equation} \phi_{-}(\eta):=\eta^{2}(\hat{\mu}_{1}+(\hat{\lambda}_{0}+\hat{\lambda}_{2})z) -\eta(\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})+\hat{\mu}_{2}z^{{-}1}+\hat{\lambda}_{1}=0,\end{equation}

and $\theta _{{\rm min}}<\theta _{{\rm max}}$ are the roots of the equation:

(5.9)\begin{equation} \phi_{+}(\theta):=\theta^{2}(\hat{\mu}_{2}+(\hat{\lambda}_{0}+\hat{\lambda}_{1})z) -\theta(\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})+\hat{\mu}_{1}z^{{-}1}+\hat{\lambda}_{2}=0.\end{equation}

Proof. See Appendix C.

Remark 7 Note that from (5.8), (5.9),

\begin{align*} \eta_{{\rm min}}& =\frac{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}-\sqrt{\beta_{1,2}(z)}}{2(\hat{\mu}_{1}+(\hat{\lambda}_{2}+\hat{\lambda}_{0})z)},\quad \eta_{{\rm max}}=\frac{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}+\sqrt{\beta_{1,2}(z)}}{2(\hat{\mu}_{1}+(\hat{\lambda}_{2}+\hat{\lambda}_{0})z)},\\ \theta_{{\rm min}}& =\frac{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}-\sqrt{\beta_{2,1}(z)}}{2(\hat{\mu}_{2}+(\hat{\lambda}_{1}+\hat{\lambda}_{0})z)},\quad \theta_{{\rm max}}=\frac{\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}+\sqrt{\beta_{2,1}(z)}}{2(\hat{\mu}_{2}+(\hat{\lambda}_{1}+\hat{\lambda}_{0})z)}. \end{align*}

We now turn our attention in the solution of the equation $f(z)=1$. To improve the readability, set $\gamma _{1}=\hat {\mu }_{1}\rho ^{2}+\hat {\lambda }_{2}$, $\gamma _{2}=\hat {\mu }_{2}\rho ^{2}+\hat {\lambda }_{1}$. Note that $\gamma _{1}+\gamma _{2}+\hat {\lambda }_{0}=\rho (\hat {\lambda }+\hat {\mu }_{1}+\hat {\mu }_{2})$.

Lemma 5 $z=\rho ^{-2}$ is the only solution of $f(z)=1$, $z>1$, and

(5.10)\begin{align} (\eta_{{\rm min}},\eta_{{\rm max}})& =\left(\rho\frac{\gamma_{2}}{\gamma_{1}+\hat{\lambda}_{0}},\rho\right)\text{ or }\left(\rho,\rho\frac{\gamma_{2}}{\gamma_{1}+\hat{\lambda}_{0}}\right), \end{align}

(5.11)\begin{align} (\theta_{{\rm min}},\theta_{{\rm max}})& =\left(\rho\frac{\gamma_{1}}{\gamma_{2}+\hat{\lambda}_{0}},\rho\right)\text{ or }\left(\rho,\rho\frac{\gamma_{1}}{\gamma_{2}+\hat{\lambda}_{0}}\right). \end{align}

Proof. See Appendix D.

The results given in Lemmas 4 and 5 provide the basic ingredients for the proof of Theorem 4.1. Note that based on Lemma 5, there are four cases for the expression of $\mathbf {x}(\rho ^{-2})$ based on the possible values of $\eta _{{\rm min}}$, $\theta _{{\rm min}}$. In particular,

\begin{equation*}\begin{pmatrix} \eta_{{\rm min}}\\ \theta_{{\rm min}} \end{pmatrix}=\begin{pmatrix} \rho\frac{\gamma_{2}}{\gamma_{1}+\hat{\lambda}_{0}}\\ \rho\frac{\gamma_{1}}{\gamma_{2}+\hat{\lambda}_{0}} \end{pmatrix}\text{ or }\begin{pmatrix} \rho\\ \rho\frac{\gamma_{1}}{\gamma_{2}+\hat{\lambda}_{0}} \end{pmatrix}\text{ or }\begin{pmatrix} \rho\frac{\gamma_{2}}{\gamma_{1}+\hat{\lambda}_{0}}\\ \rho\end{pmatrix}\text{ or }\begin{pmatrix} \rho\\ \rho\end{pmatrix}.\end{equation*}

The last one is rejected since corresponds to the case $\gamma _{1}+\hat {\lambda }_{0}\leq \gamma _{2}$, $\gamma _{2}+\hat {\lambda }_{0}\leq \gamma _{1}$, which is impossible since $\lambda _{0}>0$. It is easily seen that only the first pair satisfies $\mathbf {x}(\rho ^{-2})=\mathbf {x}(\rho ^{-2})A_{*}(\rho ^{-2})$, which means that $A_{*}(\rho ^{-2})$ is 1-positive. Note that in such a case $\rho ({\gamma _{2}}/{(\gamma _{1}+\hat {\lambda }_{0})})<\rho \Rightarrow \gamma _{2}-\gamma _{1}<\hat {\lambda }_{0}$, and $\rho ({\gamma _{1}}/{(\gamma _{2}+\hat {\lambda }_{0})})<\rho \Rightarrow \gamma _{1}-\gamma _{2}<\hat {\lambda }_{0}$, which is equivalent to $|\gamma _{1}-\gamma _{2}|<\hat {\lambda }_{0}$ or to $|\hat {\lambda }_{2}-\hat {\lambda }_{1}+\rho ^{2}(\mu _{1}-\mu _{2})|<\hat {\lambda }_{0}$ (note that this case is equivalent to the so-called strongly pooled case in [15]).

Now it remains to verify the boundary condition (5.4). From Lemma 5, we know that $\eta _{{\rm max}}=\theta _{{\rm max}}=\rho$. Thus, $\boldsymbol {\pi }_{0}\mathbf {y}<\infty$ implies

$$\sum_{l=0}^{\infty}\rho^{{-}l}(\mathbb{P}(N_{1}=l,N_{2}=0,C=1)+\mathbb{P}(N_{1}=0,N_{2}=l,C=1))<\infty,$$

where $\{(N_{1},N_{2},C)\}$ is the stationary version of $\{X(n);n\geq 0\}$. The last inequality is the same as condition C.10 in [15], which was proved under the corresponding conditions $\max \{\rho _{1},\rho _{2}\}<\rho <1$, and $|\hat {\lambda }_{2}-\hat {\lambda }_{1}+\rho ^{2}(\mu _{1}-\mu _{2})|<\hat {\lambda }_{0}$ in [15] Proposition 1, following [28] Theorem 14.3.7, so further details are omitted.

Following [26] Definition 2.2, the two orbit queue are strongly balanced, if the difference between orbit queue lengths (i.e., the background process) is dominated by the level process, i.e., $\min \{N_{1,n},N_{2,n}\}$. Intuitively, it means that the level process has more influence when the same scaling is applied to the level and background process. This is equivalent with the condition

$$\lim_{l\to+\infty}\rho^{{-}2l}x_{l}(\rho^{{-}2})=\lim_{l\to-\infty}\rho^{2l}x_{l}(\rho^{{-}2})=0.$$

Noting (4.3), it is easily realized that this is the case if and only if $\gamma _{2}<\rho (\gamma _{1}+\hat {\lambda }_{0})$, $\gamma _{1}<\rho (\gamma _{2}+\hat {\lambda }_{0})$.

5.2. Proof of Lemma 2

Starting with region I ($i\gg 1$, $j\gg 1$), Eqs. (2.1) and (2.2) are the only relevant. By introducing the new variables $m=i+j$, $n=j-i$, and setting $p_{i,j}(k)\equiv R_{m,n}(k)$, (2.1) and (2.2) become

(5.12)\begin{equation} R_{m,n}(0)=\frac{\mu}{\lambda+2\alpha}R_{m,n}(1),\end{equation}

(5.13)\begin{align}(\lambda+2\alpha)R_{m,n}(1)& =[\lambda_{0}H(1-n)+\lambda_{+}]R_{m-1,n-1}(1)\nonumber\\ & \quad +[\lambda_{0}H(n+1)+\lambda_{+}]R_{m-1,n+1}(1)+\lambda R_{m,n}(0)\nonumber\\ & \quad +\alpha[R_{m+1,n-1}(0)+R_{m+1,n+1}(0)], \end{align}

Substituting (5.12) to (5.13) yields

(5.14)\begin{align} & (\hat{\lambda}+2\hat{\mu})R_{m,n}(1)=[\hat{\lambda}_{0}H(1-n)+\hat{\lambda}_{+}]R_{m-1,n-1}(1) \nonumber\\ & \quad +[\hat{\lambda}_{0}H(n+1)+\hat{\lambda}_{+}]R_{m-1,n+1}(1) +\hat{\mu}[R_{m+1,n-1}(1)+R_{m,n}(1)], \end{align}

where $\hat {\lambda }=\lambda (\lambda +2\alpha )$, $\lambda =\lambda _{0}+2\lambda _{+}$, $\hat {\lambda }_{l}=\lambda _{l}(\lambda +2\alpha )$, $l=0,+$, $\hat {\mu }=\alpha \mu$. Note that regarding variable $m$, (5.14) is second-order difference equation with constant coefficients, which admits solutions of the form $R(m,n,1)=\gamma ^{m}Q(n)$, where $\gamma$ a constant to be determined. Substituting in (5.14), we realize that $Q(n)$ satisfies

(5.15)\begin{equation} (\hat{\lambda}+2\hat{\mu})Q(n)=\left[\frac{\hat{\lambda}_{0}H(1-n)+\hat{\lambda}_{+}}{\gamma} +\hat{\mu}\gamma\right]Q(n-1)+\left[\frac{\hat{\lambda}_{0}H(1+n)+\hat{\lambda}_{+}}{\gamma} +\hat{\mu}\gamma\right]Q(n+1). \end{equation}

Note that symmetry implies that $Q(n)=Q(-n)$, and thus, it is sufficient to only consider $n\geq 0$. Then, for $n=0$, (5.14) gives

(5.16)\begin{equation} (\hat{\lambda}+2\hat{\mu})Q(0)=2\left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma} +\hat{\mu}\gamma\right]Q(1).\end{equation}

Setting $n=1$ in (5.15), and using (5.16) we obtain,

(5.17)\begin{equation} \left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]Q(2) =\left[\hat{\lambda}+2\hat{\mu}-\frac{2}{\hat{\lambda}+2\hat{\mu}} \left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right] \left[\frac{\hat{\lambda}}{2\gamma}+\hat{\mu}\gamma\right]\right]Q(1), \end{equation}

while, for $n\geq 2$,

(5.18)\begin{equation} (\hat{\lambda}+2\hat{\mu})Q(n)=\left[\frac{\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]Q(n-1) +\left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]Q(n+1).\end{equation}

The form of (5.18) implies that $Q(n)=c\delta ^{n}$, $n\geq 1$. Thus, by substituting in (5.17), (5.18) we obtain

(5.19)\begin{align} \hat{\lambda}+2\hat{\mu}& =\left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]\delta +\frac{2}{\hat{\lambda}+2\hat{\mu}}\left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right] \left[\frac{\hat{\lambda}}{2\gamma}+\hat{\mu}\gamma\right], \end{align}

(5.20)\begin{align} \hat{\lambda}+2\hat{\mu}& =\left[\frac{\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]\frac{1}{\delta} +\left[\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\gamma}+\hat{\mu}\gamma\right]\delta. \end{align}

Eqs. (5.19) and (5.20) constitute a pair of algebraic curves in $(\gamma,\delta )-plane$ with four intersection points. From these points, the only possible candidates are: $(1,{(\hat {\lambda }_{+}+\hat {\mu })}/{(\hat {\lambda }_{0}+\hat {\lambda }_{+}+\hat {\mu })})$, $({\hat {\lambda }}/{2\hat {\mu }},{(\hat {\lambda }^{2}+4\hat {\lambda }_{+}\hat {\mu })}/ {(\hat {\lambda }^{2}+4(\hat {\lambda }_{0}+\hat {\lambda }_{+})\hat {\mu })})$. Note that $\gamma =1$ yields $\sum _{i,j\geq 0}p_{i,j}(1)=\infty$, i.e., corresponds to an un-normalizable solution for $p_{i,j}(1)$, thus the pair $({\hat {\lambda }}/{2\hat {\mu }},{(\hat {\lambda }^{2}+4\hat {\lambda }_{+}\hat {\mu })}/ {(\hat {\lambda }^{2}+4(\hat {\lambda }_{0}+\hat {\lambda }_{+})\hat {\mu })})$, with $\hat {\lambda }<2\hat {\mu }$ corresponds to the feasible pair. Note that $\hat {\lambda }<2\hat {\mu }$ corresponds to the stability condition for our system. Finally, from (5.16) by substituting $(\gamma,\delta )=({\hat {\lambda }}/{2\hat {\mu }}, {(\hat {\lambda }^{2}+4\hat {\lambda }_{+}\hat {\mu })}/{(\hat {\lambda }^{2}+4(\hat {\lambda }_{0}+\hat {\lambda }_{+})\hat {\mu })})$ yields

$$Q(0)=c\frac{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}{\hat{\lambda}(\hat{\lambda}+2\hat{\mu})}.$$

To summarize, the asymptotic expansions of $p_{i,j}(k)$, in region I ($i\gg 1$, $j\gg 1$) are

(5.21)\begin{equation} \begin{aligned} & p_{i,j}(1)\sim c\left\{\begin{array}{ll} \left(\dfrac{\hat{\lambda}}{2\hat{\mu}}\right)^{i+j}\left(\dfrac{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}{\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu}})\right)^{|i-j|}, & i\neq j,\\ \left(\dfrac{\hat{\lambda}}{2\hat{\mu}}\right)^{i+j}\left(\dfrac{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}{\hat{\lambda}(\hat{\lambda}+2\hat{\mu})}\right), & i=j.\end{array}\right.\\ & p_{i,j}(0)=\frac{\mu}{\lambda+2\alpha}p_{i,j}(1), \end{aligned} \end{equation}

where the multiplicative constant $c$ is the normalization constant.

We proceed with region II for $i=0$ or 1, and $j\gg 1$. Eqs. (2.1) and (2.2) for $j>i+1$, and (2.3) are the only relevant. Using these equations, we have to solve

(5.22)\begin{align} & (\lambda+\alpha 1_{\{i>0\}}+\alpha)p_{i,j}(0)=\mu p_{i,j}(1), \end{align}

(5.23)\begin{align} & (\hat{\lambda}+2\hat{\mu})p_{i,j}(1)=\hat{\lambda}_{+}p_{i,j-1}(1) +[\hat{\lambda}_{0}+\hat{\lambda}_{+}]p_{i-1,j}(1)+\hat{\mu}[p_{i+1,j}(1)+p_{i,j+1}(1)], \end{align}

(5.24)\begin{align} & (\lambda(\lambda+\alpha)+\hat{\mu})p_{0,j}(1)=\frac{\hat{\mu}(\lambda+\alpha)}{\lambda+2\alpha}p_{1,j}(1) +\hat{\mu}p_{0,j+1}(1)+\lambda_{+}(\lambda+\alpha)p_{0,j-1}(1). \end{align}

Clearly, the solution of (5.22)–(5.24) should agree with the expansion (5.21) for $i\gg 1$. Thus, we should add the condition

$$p_{i,j}(1)\sim c\left(\frac{\hat{\lambda}}{2\hat{\mu}}\right)^{i+j} \left(\frac{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}{\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu}}\right)^{j-i},\quad i\gg 1.$$

This condition implies that we seek for solutions of (5.22)–(5.24) of the form

(5.25)\begin{equation} p_{i,j}(1)=c\left(\frac{\hat{\lambda}}{2\hat{\mu}}\times \frac{\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu}}{\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu}}\right)^{j}F(i), \end{equation}

where $F(i)\sim [{\hat {\lambda }(\hat {\lambda }^{2}+4\hat {\lambda }_{+}\hat {\mu })}/{(2\hat {\mu }(\hat {\lambda }^{2}+4(\hat {\lambda }_{0}+\hat {\lambda }_{+})\hat {\mu }))}]^{i}$, $i\gg 1$. Substitute (5.25) in (5.23)and (5.24), we obtain the following set of equations that $F(i)$ should satisfy:

(5.26)\begin{align} (\lambda(\lambda+\alpha)+\hat{\mu})F(0)& =\hat{\mu}\left[\frac{\lambda+\alpha}{\lambda+2\alpha}F(1) +\frac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\lambda}_{+}\hat{\mu})} {2\hat{\mu}(\hat{\lambda}^{2}+4(\hat{\lambda}_{0}+\hat{\lambda}_{+})\hat{\mu})}F(0)\right]\nonumber\\ & \quad +\lambda_{+}(\lambda+\alpha)\frac{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))} {\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})}F(0), \end{align}

(5.27)\begin{align}(\hat{\lambda}+2\hat{\mu})F(i)& =(\hat{\lambda}_{0}+\hat{\lambda}_{+})F(i-1)+\hat{\lambda}_{+} \frac{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})}F(i) \nonumber\\ & \quad +\hat{\mu}F(i+1)+\frac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+})}{2(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}F(i),\quad i\geq 1. \end{align}

The form of (5.27) implies that $F(i)=A_{+}x_{+}^{i}+A_{-}x_{-}^{i}$, where

\begin{align*} x_{+}& =\frac{\hat{\lambda}(\hat{\lambda}^{2}+4\hat{\mu}(\hat{\lambda}_{0}+\hat{\lambda}_{+}))}{2\hat{\mu}(\hat{\lambda}^{2}+4\hat{\mu}\hat{\lambda}_{+}},\quad x_{-}=\frac{\hat{\lambda}_{0}+\hat{\lambda}_{+}}{\hat{\mu}x_{+}},\\ A_{-}& =\frac{\lambda(\lambda+\alpha)+\hat{\mu}(1-x_{-}(1+\frac{\lambda+\alpha}{\lambda+2\alpha})) +\frac{\lambda_{+}(\lambda+\alpha)}{x_{+}}}{\hat{\mu}(\frac{\lambda+\alpha}{\lambda+2\alpha}x_{-}+x_{+}-1) -\lambda(\lambda+\alpha)-\frac{\lambda_{+}(\lambda+\alpha)}{x_{+}}}. \end{align*}

Moreover, $A_{+}=1$ (due to the asymptotic form of $F(i)$ given below (5.25)), while $A_{-}$ is derived by the boundary condition (5.26). The obtained expansion for region II coincides with the expression (5.21) for $i\gg 1$. Thus, (5.25) is uniformly valid for all $i< j$ with $j\gg 1$. Due to the symmetry of the model, the expansion for $i>j$, $j\gg 1$ is obtained by replacing in the derived expressions $(i,j)$ with $(j,i)$.

5.3. Proof of Lemma 3

The proof is similar to the one presented in Section 5.2, although we need an additional condition to ensure that the approximated joint orbit queue-length distribution is well defined. We start with region I ($i\gg 1,j\gg 1$), where Eqs. (2.1) and (2.2) become

(5.28)\begin{align} R_{m,n}(0)& =\frac{\mu}{\lambda+\alpha_{1}+\alpha_{2}}R_{m,n}(1), \end{align}

(5.29)\begin{align} (\lambda+\alpha_{1}+\alpha_{2})R_{m,n}(1)& =[\lambda_{0}H(1-n)+\lambda_{2}]R_{m-1,n-1}(1)\nonumber\\ & \quad +[\lambda_{0}H(n+1)+\lambda_{1}]R_{m-1,n+1}(1) +\lambda R_{m,n}(0)\nonumber\\& \quad+\alpha_{1}R_{m+1,n-1}(0)+\alpha_{2}R_{m+1,n+1}(0)], \end{align}

which give

(5.30)\begin{align} (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})R_{m,n}(1)& =[\hat{\lambda}_{0}H(1-n)+\hat{\lambda}_{2}]R_{m-1,n-1}(1) \nonumber\\ & \quad +[\hat{\lambda}_{0}H(n+1)+\hat{\lambda}_{2}]R_{m-1,n+1}(1) \nonumber\\ & \quad +\hat{\mu}_{1}R_{m+1,n-1}(1)+\hat{\mu}_{2}R_{m,n}(1), \end{align}

where $\hat {\lambda }=\lambda (\lambda +\alpha _{1}+\alpha _{2})$, $\hat {\lambda }_{l}=\lambda _{l}(\lambda +\alpha _{1}+\alpha _{2})$, $l=0,1,2$, $\hat {\mu }_{l}=\alpha _{l}\mu$, $l=1,2$. Contrary to the symmetric case, we are now seeking for solutions of (5.30) of the form

(5.31)\begin{equation} R_{m,n}(1)=\gamma^{m}\left\{\begin{array}{ll} c_{+}\delta_{+}^{n}, & n\geq 1,\\ c_{0}, & n=0,\\ c_{-}\delta_{-}^{{-}n}, & n\leq{-}1. \end{array}\right. \end{equation}

Substitution of (5.31) to (5.30) yields five equations, corresponding to the cases $n=0$, $n=\pm 1$, $n\geq 2$, $n\leq -2$:

(5.32)\begin{align} (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})c_{0}& =\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{2}}{\gamma} +\hat{\mu}_{1}\gamma\right)c_{-}\delta_{-}+\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{1}}{\gamma} +\hat{\mu}_{2}\gamma\right)c_{+}\delta_{+}, \end{align}

(5.33)\begin{align} (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})c_{+}\delta_{+} & =\left(\frac{\hat{\lambda}_{0}+2\hat{\lambda}_{2}}{2\gamma}+\hat{\mu}_{1}\gamma\right)c_{0} +\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{1}}{\gamma}+\hat{\mu}_{2}\gamma\right)c_{+}\delta_{+}^{2}, \end{align}

(5.34)\begin{align} (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})c_{-}\delta_{-} & =\left(\frac{\hat{\lambda}_{0}+2\hat{\lambda}_{1}}{2\gamma}+\hat{\mu}_{2}\gamma\right)c_{0} +\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{2}}{\gamma}+\hat{\mu}_{1}\gamma\right)c_{-}\delta_{-}^{2}, \end{align}

(5.35)\begin{align} \hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}& =\left(\frac{\hat{\lambda}_{2}}{\gamma}+\hat{\mu}_{1}\gamma\right) \frac{1}{\delta_{+}}+\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{1}}{\gamma}+\hat{\mu}_{2}\gamma\right)\delta_{+}, \end{align}

(5.36)\begin{align} \hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2}& =\left(\frac{\hat{\lambda}_{1}}{\gamma}+\hat{\mu}_{2}\gamma\right) \frac{1}{\delta_{-}}+\left(\frac{\hat{\lambda}_{0}+\hat{\lambda}_{2}}{\gamma}+\hat{\mu}_{1}\gamma\right)\delta_{-}. \end{align}

The system (5.32)–(5.36) is a set of five equations for the five unknowns $\gamma$, $\delta _{+}$, $\delta _{-}$, $c_{+}/c_{0}$, $c_{-}/c_{0}$. Using (5.35) in (5.33) and (5.36) in (5.34), we easily realize that

(5.37)\begin{equation} c_{-}=\frac{\epsilon_{+}}{\epsilon_{-}}c_{+}. \end{equation}

where

(5.38)\begin{equation} \epsilon_{+}=\left(\frac{\hat{\lambda}_{2}}{\gamma}+\hat{\mu}_{1}\gamma\right) \left(\frac{\hat{\lambda}_{0}+2\hat{\lambda}_{1}}{2\gamma}+\hat{\mu}_{2}\gamma\right),\quad \epsilon_{-}=\left(\frac{\hat{\lambda}_{1}}{\gamma}+\hat{\mu}_{2}\gamma\right) \left(\frac{\hat{\lambda}_{0}+2\hat{\lambda}_{2}}{2\gamma}+\hat{\mu}_{1}\gamma\right). \end{equation}

Thus, a natural choice is $c_{-}:=\sqrt {{\epsilon _{+}}/{\epsilon _{-}}}$, and $c_{+}:=\sqrt {{\epsilon _{-}}/{\epsilon _{+}}}$. Note that in the symmetric case, $c_{+}=c_{-}$. Following a procedure similar to the one used in [22], after heavy algebraic manipulations (similar to those given in the proof of Lemma 5) we realize that $\gamma :=\rho ={\hat {\lambda }}/{(\hat {\mu }_{1}+\hat {\mu }_{2})}$, which is crucial for the stability condition, thus, asking $\hat {\lambda }<\hat {\mu }_{1}+\hat {\mu }_{2}$.

For $\gamma =\rho$, (5.35), (5.36) have each one from two roots, namely $\delta _{+}=1$ or $\delta _{+}={(\hat {\lambda }_{2}+\rho ^{2}\hat {\mu }_{1})}/{(\hat {\lambda }_{0}+\hat {\lambda }_{1}+\rho ^{2}\hat {\mu }_{2})}$, and $\delta _{-}=1$ or $\delta _{-}={(\hat {\lambda }_{1}+\rho ^{2}\hat {\mu }_{2})}/{(\hat {\lambda }_{0}+\hat {\lambda }_{2}+\rho ^{2}\hat {\mu }_{1})}$, respectively. Therefore, the feasible candidates are $\delta _{+}={(\hat {\lambda }_{2}+\rho ^{2}\hat {\mu }_{1})}/{(\hat {\lambda }_{0}+\hat {\lambda }_{1}+\rho ^{2}\hat {\mu }_{2})}$, $\delta _{-}={(\hat {\lambda }_{1}+\rho ^{2}\hat {\mu }_{2})}/{(\hat {\lambda }_{0}+\hat {\lambda }_{2}+\rho ^{2}\hat {\mu }_{1})}$. Note, that when $\hat {\lambda }_{0}>|\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|$, both $0<\delta _{+}<1$, and $0<\delta _{-}<1$. Note that this case along with the assumption that $\rho >\max \{\rho _{1},\rho _{2}\}$, where $\rho _{j}=\hat {\lambda }_{j}/\hat {\mu }_{j}$, $j=1,2$, corresponds to the so called strongly pooled case mentioned in [15] Theorem 2 for the standard join the shortest queue model without retrials. Substituting $\gamma$, $\delta _{+}$, $\delta _{-}$ in (5.37) and (5.32), we obtain after some algebra

\begin{align*} c_{0}& =\frac{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\hat{\lambda}(1+\rho)} \sqrt{\frac{\epsilon_{+}}{\epsilon_{-}}}+\frac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}(1+\rho)}\sqrt{\frac{\epsilon_{-}}{\epsilon_{+}}},\\ \epsilon_{-}& =\frac{(\hat{\lambda}_{1}+\hat{\mu}_{2}\rho^{2})(\hat{\lambda}_{0}+2\hat{\lambda}_{2}+2\hat{\mu}_{1}\rho^{2})}{2\rho^{2}},\quad \epsilon_{+}=\frac{(\hat{\lambda}_{2}+\hat{\mu}_{1}\rho^{2})(\hat{\lambda}_{0}+2\hat{\lambda}_{1}+2\hat{\mu}_{2}\rho^{2})}{2\rho^{2}}. \end{align*}

Thus, for region I ($i\gg 1,j\gg 1$), our results are summarized in

(5.39)\begin{equation} \begin{aligned} & p_{i,j}(1)\sim c\rho^{i+j}\left\{\begin{array}{ll} \sqrt{\dfrac{\epsilon_{-}}{\epsilon_{+}}}\left(\dfrac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{j-i}, & j>i,\\ \dfrac{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\hat{\lambda}(1+\rho)} \sqrt{\dfrac{\epsilon_{+}}{\epsilon_{-}}}+\dfrac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}(1+\rho)}\sqrt{\dfrac{\epsilon_{-}}{\epsilon_{+}}}, & i=j,\\ \sqrt{\dfrac{\epsilon_{-}}{\epsilon_{+}}}\left(\dfrac{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{1}}\right)^{i-j}, & j< i. \end{array}\right.\\ & p_{i,j}(0)=\frac{\mu}{\lambda+\alpha_{1}+\alpha_{2}}p_{i,j}(1), \end{aligned} \end{equation}

where $c$ is a multiplicative constant.

Next, we consider region II for $i=0$ or $1$, and $j\gg 1$. We should solve

(5.40)\begin{align} & (\lambda+\alpha_{1}1_{\{i>0\}}+\alpha_{2})p_{i,j}(0)=\mu p_{i,j}(1), \end{align}

(5.41)\begin{align} & (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})p_{i,j}(1)=\hat{\lambda}_{2}p_{i,j-1}(1) +[\hat{\lambda}_{0}+\hat{\lambda}_{1}]p_{i-1,j}(1)+\hat{\mu}_{1}p_{i+1,j}(1)+\hat{\mu}_{2}p_{i,j+1}(1), \end{align}

(5.42)\begin{align} & (\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2})p_{0,j}(1)=\frac{\hat{\mu}_{1}(\lambda+\alpha_{2})}{\lambda+\alpha_{1}+\alpha_{2}}p_{1,j}(1) +\hat{\mu}_{2}p_{0,j+1}(1)+\lambda_{2}(\lambda+\alpha_{2})p_{0.j-1}(1), \end{align}

and

$$p_{i,j}(1)\sim c\sqrt{\frac{\epsilon_{-}}{\epsilon_{+}}} \left(\rho\frac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{j} \left(\rho\frac{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}\right)^{i},\quad i\gg 1.$$

Thus, we seek solutions of (5.41)–(5.42) of the form

$$p_{i,j}(1)\sim c\sqrt{\frac{\epsilon_{-}}{\epsilon_{+}}} \left(\rho\frac{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{j}F(i),$$

and obtain the following set of equations for $F(i)$:

(5.43)\begin{align} (\hat{\lambda}+\hat{\mu}_{1}+\hat{\mu}_{2})F(i)& =(\hat{\lambda}_{0}+\hat{\lambda}_{1})F(i-1) +\left(\frac{\hat{\lambda}_{2}(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})} +\frac{\hat{\mu}_{2}\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right) F(i)\nonumber\\ & \quad +\hat{\mu}_{1}F(i+1), \end{align}

(5.44)\begin{align} (\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2})F(0)& =\hat{\mu}_{1}\frac{\lambda+\alpha_{2}}{\lambda+\alpha_{1}+\alpha_{2}}F(1) +\frac{\hat{\mu}_{2}\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}F(0)\nonumber\\ & \quad +\lambda_{2}(\lambda+\alpha_{2})\frac{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}{\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}F(0), \end{align}

with $F(i)\sim ({\rho (\hat {\lambda }_{0}+\hat {\lambda }_{1}+\rho ^{2}\hat {\mu }_{2})}/{(\hat {\lambda }_{2}+\rho ^{2}\hat {\mu }_{1})})^{i}$, $i\gg 1$. The solution to (5.43), (5.44) is $F(i)=x_{+}^{i}+A_{-}x_{-}^{i}$, where

\begin{align*} x_{+}& =\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}},\quad x_{-}=\frac{(\hat{\lambda}_{0}+\hat{\lambda}_{1})(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\mu}_{1}\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})},\\ A_{-}& =\frac{x_{+}^{2}(\hat{\mu}_{1}\frac{\lambda+\alpha_{2}}{\lambda+\alpha_{1}+\alpha_{2}}+\frac{\lambda_{2}(\lambda+\alpha_{2})}{\rho^{2}})-x_{+}(\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2})+\rho^{2}\hat{\mu}_{2}}{x_{+}(\lambda(\lambda+\alpha_{2})+\hat{\mu}_{2})-[(\lambda_{0}+\lambda_{1})(\lambda+\alpha_{2})+\rho^{2}\hat{\mu}_{2}]-\frac{\lambda_{2}(\lambda+\alpha_{2})}{\rho^{2}}x_{+}^{2}}, \end{align*}

where $A_{-}$ is obtained using (5.44). Moreover, $x_{+}>1$ when $\hat {\lambda }_{2}+\hat {\mu }_{1}\rho ^{2}<\rho (\hat {\lambda }_{0}+\hat {\lambda }_{1}+\hat {\mu }_{2}\rho ^{2})$ (note that this corresponds to the strongly balanced condition mentioned in Theorem 4.1). Furthermore, $x_{-}/x_{+}<1$ when $(\hat {\lambda }_{2}+\hat {\mu }_{1}\rho ^{2}) \sqrt {{(\hat {\lambda }_{0}+\hat {\lambda }_{1})}/{\hat {\mu }_{1}}}<\rho (\hat {\lambda }_{0}+\hat {\lambda }_{1}+\hat {\mu }_{2}\rho ^{2})$. Thus, for $i< j$, $j\gg 1$,

(5.45)\begin{equation} \begin{aligned} & p_{i,j}(1)\sim c\sqrt{\frac{\epsilon_{-}}{\epsilon_{+}}}\left(\frac{\rho(\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{i} \left(\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}\right)^{j} \left[1+A_{-}\left(\frac{x_{-}}{x_{+}}\right)^{i}\right],\\ & p_{i,j}(0)=\frac{\mu}{\lambda+\alpha_{1}1_{\{i>0\}}+\alpha_{2}}p_{i,j}(1). \end{aligned} \end{equation}

An analogous analysis can be performed for region III ($i\gg 1$, $j=0,1$), and results in

(5.46)\begin{equation} \begin{aligned} p_{i,j}(1) & \sim c\sqrt{\frac{\epsilon_{+}}{\epsilon_{-}}}\left(\frac{\rho(\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1}}\right)^{i} \left(\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\right)^{j}\\ & \quad \times \left[1+B_{-}\left(\frac{y_{-}}{y_{+}}\right)^{j}\right],\quad i>j,\ i\gg 1,\\ p_{i,j}(0) & =\frac{\mu}{\lambda+\alpha_{2}1_{\{j>0\}}+\alpha_{1}}p_{i,j}(1), \end{aligned} \end{equation}

where,

\begin{align*} y_{-}& =\frac{(\hat{\lambda}_{0}+\hat{\lambda}_{2})(\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2})}{\hat{\mu}_{2}\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})},\quad y_{+}=\frac{\rho(\hat{\lambda}_{0}+\hat{\lambda}_{2}+\rho^{2}\hat{\mu}_{1})}{\hat{\lambda}_{1}+\rho^{2}\hat{\mu}_{2}}\\ B_{-}& =\frac{y_{+}^{2}(\hat{\mu}_{2}\frac{\lambda+\alpha_{1}}{\lambda+\alpha_{1}+\alpha_{2}}+ \frac{\lambda_{1}(\lambda+\alpha_{1})}{\rho^{2}})-y_{+}(\lambda(\lambda+\alpha_{1})+\hat{\mu}_{1})+\rho^{2}\hat{\mu}_{1}}{y_{+}(\lambda(\lambda+\alpha_{1})+\hat{\mu}_{1}) -[(\lambda_{0}+\lambda_{2})(\lambda+\alpha_{1})+\rho^{2}\hat{\mu}_{1}] -\frac{\lambda_{1}(\lambda+\alpha_{1})}{\rho^{2}}y_{+}^{2}}. \end{align*}

Moreover, $y_{+}>1$ when $\hat {\lambda }_{1}+\hat {\mu }_{2}\rho ^{2}<\rho (\hat {\lambda }_{0}+\hat {\lambda }_{2}+\hat {\mu }_{1}\rho ^{2})$ (note that this corresponds to the strongly balanced condition mentioned in Theorem 4.1), and $y_{-}/y_{+}<1$ when $(\hat {\lambda }_{1}+\hat {\mu }_{2}\rho ^{2}) \sqrt {{(\hat {\lambda }_{0}+\hat {\lambda }_{2})}/{\hat {\mu }_{2}}} <\rho (\hat {\lambda }_{0}+\hat {\lambda }_{2}+\hat {\mu }_{1}\rho ^{2})$.

6. Numerical results

In the following, we numerically validate and compare the theoretical results obtained based on asymptotic analysis in subsection 4.1, with those obtained by the heuristic approach in 4.2. We will see that as $m\to \infty$, i.e., $\min \{N_{1,n},N_{2,n}\}\to \infty$ the expressions derived from the tail asymptotic analysis coincide with the heuristic approximation expressions. Moreover, we notice that even when $m$ takes small values, the difference of the derived expressions is negligible.

Set $\lambda _{0}=0.15$, $\lambda _{1}=0.05$, $\lambda _{2}=0.01$, $\mu =0.44$. Then, for $\alpha _{1}=0.25$, $\alpha _{2}=0.1$ we investigate whether the results obtained by the asymptotic analysis presented in subsection 5.1 agreed with those obtained in subsection 4.2 through our heuristic approach. In this direction, we focus on the absolute difference $|\mathbb {P}(N_{1}=i,N_{2}=j)-\mathbb {P}(\min \{i,j\},j-i)|$, where $\mathbb {P}(N_{1}=i,N_{2}=j)=p_{i,j}(0)+p_{i,j}(1)$ obtained with the aid of Lemma 3 in subsection 4.2, and $\mathbb {P}(\min \{i,j\}=m,j-i=l)=\pi _{m,l}(0)+\pi _{m,l}(1)$ with the aid of Theorem 4.1 in subsection 4.1. Note that under such a setting, $\rho =0.7636>\max \{\rho _{1}=0.2545,\rho _{2}=0.1273\}$, and $\hat {\lambda }_{0}-|\rho ^{2}(\hat {\mu }_{2}-\hat {\mu }_{1})+\hat {\lambda }_{1}-\hat {\lambda }_{2}|=0.0679>0$, and $\rho (\gamma _{2}+\hat {\lambda }_{0})>\gamma _{1}$, $\rho (\gamma _{1}+\hat {\lambda }_{0})>\gamma _{2}$ (i.e., the orbit queues are strongly balanced); see Table 1. Our results show that the stationary approximations through the heuristic approach agreed with the results obtained by the asymptotic analysis.

Table 1. Validation of the asymptotic stationary approximations.

We now focus on the stationary approximations derived through the heuristic approach in subsection 4.2. In Figure 10, we observe the ratio ${Pr(k+1)}/{Pr(k)}$ for increasing values of the total number of jobs in orbits, i.e., $k=N_{1}+N_{2}$, where $Pr(k)=p_{i,j}(1)+p_{i,j}(0)$ with $k=i+j$, and $\lambda _{0}=0.04$, $\lambda _{1}=\lambda _{2}=0.01$, $\mu =0.44$. It is readily seen that $\lim _{k\to \infty } {Pr(k+1)}/{Pr(k)}=\rho ={\hat {\lambda }}/{(\hat {\mu }_{1}+\hat {\mu }_{2})}$. Moreover, we can observe that the more the difference $\alpha _{1}-\alpha _{2}$ get smaller, the faster this ratio tends to $\rho$. In other words, the asymmetry of the retrial rates affects the asymptotic behavior.

Figure 10. The ratio ${Pr(k+1)}/{Pr(k)}$ for increasing values of $k=N_{1}+N_{2}$.

Moreover, we can also observe that the presence of the dedicated traffic also heavily affects the way the ratio ${Pr(k+1)}/{Pr(k)}$ tends to $\rho$. In the following table, we can observe this trend for $\lambda _{1}=0=\lambda _{2}$, $\lambda _{0}=0.06$, $\mu =0.44$, compared with the case where $\lambda _{1}=0.01=\lambda _{2}$, $\lambda _{0}=0.04$, $\mu =0.44$; see Table 2. We observe that in case of no dedicated traffic, the convergence of the ratio ${Pr(k+1)}/{Pr(k)}$ to $\rho$ is really faster compared with the case where we assumed dedicated traffic. Note that in these scenarios, we have assumed that the total arrival rate is fixed and equal to $\lambda =0.06$, so that $\rho =0.1527$ is fixed in both cases. Moreover, we note when $\rho$ increases, the ratio ${Pr(k+1)}/{Pr(k)}$ tends to $\rho$ very fast (from $k\geq 10$ the ratio is very close to $\rho$), as shown in Figure 11, where $\lambda _{0}=0.19$, $\lambda _{1}=\lambda _{1}=0.01$ $\mu =0.44$ (in the case of no dedicated traffic, we assumed $\lambda _{1}=\lambda _{2}=0$, $\lambda _{0}=0.21$).

Figure 11. The ratio ${Pr(k+1)}/{Pr(k)}$ for increasing values of $k=N_{1}+N_{2}$.

Table 2. Effect of dedicated traffic.

7. Conclusion and future work

In this work, we introduced the generalized join the shortest queue system with retrials and two infinite capacity orbit queues. Three independent streams of jobs, namely a smart, and two dedicated streams, flow into a single-server service system, which can hold at most one job. Blocked smart jobs join the shortest orbit queue, while blocked jobs from the dedicated streams join their corresponding orbit queue. We remind that the system is described as a Markov modulated two-dimensional random walk, but this kind of modulation allowed for a completely tractable analysis.

We establish the geometric tail asymptotics along the minimum direction for a fixed value of the difference among orbit queue lengths. Moreover, we apply a heuristic approach to obtain stationary approximations of the number of the joint orbit queue-length distribution, which are accurate when one of the orbit queue lengths is large, and thus agreed with the asymptotic analysis. We have shown that even though the exact solution to the problem may be quite complicated, the asymptotic expansions of $p_{i,j}(k)$ are relatively simple and clearly indicate the dependence of the stationary distribution on the parameters of the model.

We also cope with the stability condition, and based on the stability of the censored chain on the busy states along with extensive simulation experiments, we conjecture that the ergodicity conditions of the censored chain coincide with those of the original model. We postpone the formal proof of this conjecture in a future work. We also leave as a future work the complete investigation of the tail decay rate problem. More precisely, we proved the exact geometric decay under specific conditions given in Theorem 4.1 (which are similar to those given in [15]), but it is still an open problem for the investigation of the decay rates when these conditions collapse.

Finally, simulation experiments indicate that in heavy traffic (i.e., $\rho \to 1$) and when $\lambda _{0}>|\lambda _{1}-\lambda _{2}|$, $\alpha _{1}=\alpha _{2}$, our system exhibits a state-space collapse behavior; see Figure 12 (right). It seems that when the first condition fails, our model does not retain this behavior; see Figure 12 (left), where $\lambda _{1}>\lambda _{0}>\lambda _{2}$, and $\lambda _{0}<|\lambda _{1}-\lambda _{2}|$ (i.e., the impact of dedicated traffic for orbit queue 1 is very crucial). Moreover, it also seems that the condition $\alpha _{1}=\alpha _{2}$ is not crucial (see Figure 9 (left)).

Figure 12. Orbit dynamics for $\rho \to 1$ when $\lambda _{0}<|\lambda _{1}-\lambda _{2}|$ (left), and when $\lambda _{0}>|\lambda _{1}-\lambda _{2}|$ (right).

This means that in heavy-traffic and under these conditions, our load-balancing scheme collapses to a one-dimensional line where all the orbit queue lengths are equal. A similar result was proven for the GJSQ system without retrials in the seminal paper [16]. State-space collapse occurs because the smart traffic flow dominates over the dedicated traffic flows and “forces” the two orbit queues to be equal. Thus, it seems that in the heavy-traffic regime, the sum of the orbit queue lengths can be approximated by the reference system described at the end of subsection 4.1. The reference system behaves as if there is only a single-orbit queue with all the “servers” (i.e., the retrial servers) pooled together as an aggregated “server” (in standard JSQ systems this is called complete resource pooling). This result implies that under specific conditions, GJSOQ is asymptotically optimal, i.e., heavy-traffic delay optimal, since the response time in the pooled single-orbit system is stochastically less than that of a typical load-balancing system, i.e. the reference system seems to serve as a lower bound (in the stochastic sense) on the total orbit queue length of the original model. In a future study, we plan to prove formally this justification.