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 [Reference Miyazawa30], in which the method developed in [Reference Miyazawa29] 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 [Reference Li, Miyazawa and Zhao26], 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})$.

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 [Reference Miyazawa30] Theorem 1.5 by using a different framework based on [Reference Miyazawa29], 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 [Reference Li, Miyazawa and Zhao26] 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 [Reference Li, Miyazawa and Zhao26]. 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 [Reference Li, Miyazawa and Zhao26] 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.

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*}

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*}

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 [Reference Foley and McDonald15]. 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.1. Proof of Theorem 4.1

The proof follows the lines in [Reference Li, Miyazawa and Zhao26], 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 [Reference Li, Miyazawa and Zhao26]:

Proposition 1 (Theorem 2.1 in [Reference Li, Miyazawa and Zhao26])

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}

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}

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}

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 [Reference Foley and McDonald15]).

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 [Reference Foley and McDonald15], 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 [Reference Foley and McDonald15] Proposition 1, following [Reference Meyn and Tweedie28] Theorem 14.3.7, so further details are omitted.

Following [Reference Li, Miyazawa and Zhao26] 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 [Reference Knessl, Matkowsky, Schuss and Tier22], 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 [Reference Foley and McDonald15] 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})$.