1 Introduction

In last few years, there has been an increasing activity in the study of the spectral representation of quasi-birth-and-death (QBD) processes, extending the pioneering work of Karlin and McGregor [13,14,15] in the 1950s (see also the recent monograph [4]). These processes are a natural extension of the so-called birth–death chains, where the state space, instead of \(\mathbb {N}_0\), is given by pairs of the form (n, k), where \(n\in \mathbb {N}_0\) is usually called the level, while \(1\le k\le r_n\) is referred to as the phase (which may depend on the different levels). For a general setup, see [19]. The transition probability matrix (discrete-time) or the infinitesimal operator matrix (continuous-time) of the QBD process is then block tridiagonal. If \(r_n=1\) for all \(n\in \mathbb {N}_0\), then we go back to classical birth–death chains. If \(r_n=N\) for all \(n\in \mathbb {N}_0\), where N is a positive integer, then all blocks in the Jacobi matrix have the same dimension \(N\times N\). In this case, the spectral analysis can be performed using matrix-valued orthogonal polynomials (see [3, 7] for the discrete-time case and [2] for the continuous-time case). Many examples have been analyzed in this direction using spectral methods in the last few years (see [1, 3, 7, 8, 10,11,12]).

A natural source of examples of more complicated QBD processes comes from the theory of multivariate orthogonal polynomials (of dimension d), where now the number of phases is given by \(r_n=\left( {\begin{array}{c}n+d-1\\ n\end{array}}\right) \). In [6], we performed the spectral analysis in the general setting of this situation as well as obtained results about recurrence and the invariant measure of these processes in terms of the spectral measure supported on some region \(\Omega \subset \mathbb {R}^d\). We also applied our results to several examples of bivariate orthogonal polynomials (\(d=2\)), namely product orthogonal polynomials, orthogonal polynomials on a parabolic domain and orthogonal polynomials on the triangle. The aim of this paper is to continue our previous work but now we will focus on the so-called first family of bivariate Jacobi–Koornwinder polynomials (see [5, Section 2.7]), first introduced by Koornwinder [16] (see also the review paper [17] where they are called Class VI). These polynomials are supported in the so-called swallow tail region (see Fig. 1) and they are eigenfunctions of two independent differential operators of orders two and four. Some properties such as a Rodrigues-type expression or an expansion in terms of James-type zonal polynomials can be found in [18, 20]. They are considered a highly non-trivial generalization of the Jacobi polynomials. Yuan Xu proved some cubature rules for specific values of the parameters [21, 22].

The paper is organized as follows. In Sect. 2, we introduce the Jacobi–Koornwinder polynomials we will be working with. Then, we normalize the polynomials in such a way that they are equal to 1 at one of the corners of the swallow tail region [specifically at the point (1, 1)]. With this family of polynomials, we derive the coefficients of the two three-term recurrence relations (one for each variable) in terms of the coefficients of the three-term recurrence relation of the classical Jacobi polynomials on [0, 1]. In Sect. 3, we will study under what conditions we may provide a probabilistic interpretation of the linear convex combination of the two Jacobi matrices associated with the three-term recurrence relations. Under these conditions, we compute the Karlin–McGregor formula, the invariant measure and study recurrence of the family of discrete-time QBD processes. Finally, in Sect. 4, we give an urn model associated with one of the QBD processes introduced in Sect. 3, for the special case of \(\beta =\alpha \).

2 Bivariate Jacobi–Koornwinder Polynomials

In [5, Section 2.7] and [17], the Jacobi–Koornwinder polynomials are constructed in terms of the Jacobi weight function supported on \([-1,1]\). The swallow tail region \(\Omega \) is then contained in the bounded rectangle \([-2,2]\times [-1,1]\). For convenience, we consider a change of variables (\(u\mapsto 4u-2, v\mapsto 2v-1\)) such that the swallow tail region \(\Omega \) is contained inside the unit square \([0,1]\times [0,1]\). Then, the region \(\Omega \) is given by (see Fig. 1)

$$\begin{aligned} \Omega =\left\{ (u,v){:}\,2u+v-1>0, 1-2u+v>0, 2u^2-2u-v+1>0 \right\} . \end{aligned}$$
Fig. 1
figure 1

Swallow tail region \(\Omega \) where the weight function is defined

Full size image

The weight function acting on this region will be given by

$$\begin{aligned} W_{\alpha ,\beta ,\gamma }(u,v)=\frac{1}{C}(1-2u+v)^{\alpha } (2u+v-1)^{\beta } (2u^2-2u-v+1)^{\gamma }, \end{aligned}$$
(2.1)

where C is the normalizing constant

$$\begin{aligned} C=\frac{2^{\alpha +\beta -\gamma +2}\Gamma (\alpha +1)\Gamma (\beta +1)\Gamma (\gamma +1)\Gamma (2\alpha +2\gamma +2)\Gamma (2\beta +2\gamma +2)\Gamma (\alpha +\beta +\gamma +3)}{\Gamma (\alpha +\gamma +1)\Gamma (\beta +\gamma +1)\Gamma (\alpha +\beta +2\gamma +3)\Gamma (2\alpha +2\beta +2\gamma +5)}, \end{aligned}$$

such that \(\int _{\Omega }W_{\alpha ,\beta ,\gamma }(u,v)dudv=1\). To ensure integrability, we need to have \(\alpha ,\beta ,\gamma >-1\), \(\alpha +\gamma +3/2>0\) and \(\beta +\gamma +3/2>0\). This normalized constant was computed in [20, Lema 6.1].

As it was pointed out in [5, Proposition 2.7.3] the monic Jacobi–Koornwinder polynomials \(P_{n,k}^{\alpha ,\beta ,\gamma }(u,v)\) satisfy the following second-order partial differential equation (after the change of variables):

$$\begin{aligned} \mathcal {D}_{\alpha ,\beta ,\gamma }P_{n,k}^{\alpha ,\beta ,\gamma }(u,v)=-\lambda _{n,k}^{\alpha ,\beta ,\gamma }P_{n,k}^{\alpha ,\beta ,\gamma }(u,v), \end{aligned}$$

where

$$\begin{aligned} \mathcal {D}_{\alpha ,\beta ,\gamma }&=\left[ u(1-u)-(1-v)/4\right] \partial _{uu}+(1-v)(2u-1)\partial _{uv}+\left[ (2u-1)^2+v(1-2v)\right] \partial _{vv}\\ {}&\quad +\left[ \beta +\gamma +3/2-(\alpha +\beta +2\gamma +3)u\right] \partial _u+2\left[ (\beta -\alpha )u-(\alpha +\beta +\gamma +5/2)v+\alpha +1\right] \partial _v,\\ \lambda _{n,k}^{\alpha ,\beta ,\gamma }&=n(n+\alpha +\beta +2\gamma +2)+k(k+\alpha +\beta +1). \end{aligned}$$

With this partial differential equation, it is possible to generate all the monic Jacobi–Koornwinder polynomials for any values of \(\alpha , \beta , \gamma \). For the special cases of \(\gamma =\pm 1/2\), it is possible to write them in terms of classical Jacobi polynomials (see [5, Proposition 2.7.2]). Another way to compute the Jacobi–Koornwinder polynomials is using the Rodrigues-type formula found in [20, Section 5].

Let us now introduce a new set of polynomials \(Q_{n,k}^{\alpha ,\beta ,\gamma }(u,v)\) normalized in such a way that \(Q_{n,k}^{\alpha ,\beta ,\gamma }(1,1)=1\) for all \(n\ge 0\) and \(0\le k\le n\). \(Q_{n,k}^{\alpha ,\beta ,\gamma }(u,v)\) can also be defined as

$$\begin{aligned} Q_{n,k}^{\alpha ,\beta ,\gamma }(u,v)=\sigma _{n,k}^{-1}P_{n,k}^{\alpha ,\beta ,\gamma }(u,v), \end{aligned}$$

where \(\sigma _{n,k}=P_{n,k}^{\alpha ,\beta ,\gamma }(1,1)\) is given by

$$\begin{aligned} \sigma _{n,k}=\frac{2^{2k-n+1}(2\gamma +2)_{n-k-1}(\alpha +1)_k(\alpha +\gamma +3/2)_n}{(\gamma +3/2)_{n-k-1}(k+\alpha +\beta +1)_k(n+k+\alpha +\beta +2\gamma +2)_{n-k}(n+\alpha +\beta +\gamma +3/2)_k}. \nonumber \\ \end{aligned}$$
(2.2)

Here, we are using the standard notation for the Pochhammer symbol \((a)_0=1, (a)_n=a(a+1)\ldots (a+n-1), n\ge 1\). The vector polynomials \({\mathbb {Q}}_n=(Q_{n,0}^{\alpha ,\beta ,\gamma },Q_{n,1}^{\alpha ,\beta ,\gamma },\dots , Q_{n,n}^{\alpha ,\beta ,\gamma }), n\ge 0,\) satisfy the three-term recurrence relations

$$\begin{aligned} \begin{aligned} u\, {{\mathbb {Q}}}_n(u,v)&= A_{n,1} {{\mathbb {Q}}}_{n+1}(u,v)+ B_{n,1} {{\mathbb {Q}}}_n(u,v) + C_{n,1} {{\mathbb {Q}}}_{n-1}(u,v),\\ v\, {{\mathbb {Q}}}_n(u,v)&= A_{n,2} {{\mathbb {Q}}}_{n+1}(u,v)+ B_{n,2} {{\mathbb {Q}}}_n(u,v) + C_{n,2} {{\mathbb {Q}}}_{n-1}(u,v). \end{aligned} \end{aligned}$$
(2.3)

It is possible to compute explicitly the coefficients \(A_{n,i}, B_{n,i},C_{n,i}, i=1,2,\) using [20, Section 9] and (2.2). For that let us introduce the following notation:

$$\begin{aligned} \begin{aligned} a_n&=\frac{(n+\alpha +1)(n+\alpha +\beta +1)}{(2n+\alpha +\beta +1)(2n+\alpha +\beta +2)},\\ b_n&=\frac{(n+\alpha +1)(n+1)}{(2n+\alpha +\beta +1)(2n+\alpha +\beta +2)}+\frac{(n+\beta )(n+\alpha +\beta )}{(2n+\alpha +\beta )(2n+\alpha +\beta +1)},\\ c_n&=\frac{n(n+\beta )}{(2n+\alpha +\beta )(2n+\alpha +\beta +1)},\\ \delta _{n,k}&=(n-k)(n+k+\alpha +\beta +1). \end{aligned} \end{aligned}$$
(2.4)

Observe that \(a_n,b_n,c_n\) are the coefficients of the three-term recurrence relation satisfied by the classical Jacobi polynomials \(Q_n^{(\beta ,\alpha )}(x)\) on [0, 1] normalized by \(Q_n^{(\beta ,\alpha )}(1)=1\) (see [9, Section 5] for instance). In particular, we always have that \(a_n,c_n>0, b_n\ge 0,\) and \(a_n+b_n+c_n=1\), i.e., they are probabilities. The norms of these Jacobi polynomials (which will be used later) are given by

$$\begin{aligned} \Vert Q_n^{(\beta ,\alpha )}\Vert _w^2=\frac{\Gamma (\alpha +1)\Gamma (\alpha +\beta +2)\Gamma (n+1)\Gamma (n+\beta +1)}{\Gamma (\beta +1)\Gamma (n+\alpha +1)\Gamma (n+\alpha +\beta +1)(2n+\alpha +\beta +1)}, \end{aligned}$$
(2.5)

where w is the normalized Jacobi weight (see (5.2) of [9]).

On one side, the matrices \(A_{n,1}\), \(B_{n,1}\) and \(C_{n,1}\) in (2.3) are of the form

$$\begin{aligned} \begin{array}{c} A_{n,1}=\left[ \begin{array}{ccccccc} a_{n,0} &{}\quad &{}\quad &{}\quad &{}\quad 0 \\ &{}\quad a_{n,1} &{}\quad &{}\quad &{}\quad \vdots \\ &{}\quad &{}\quad \ddots &{}\quad &{}\quad \\ &{}\quad &{}\quad &{}\quad a_{n,n} &{}\quad 0 \end{array} \right] , \quad C_{n,1}=\left[ \begin{array}{cccccccc} c_{n,0} &{}\quad &{}\quad &{}\quad \\ &{}\quad c_{n,1} &{}\quad &{}\quad \\ &{}\quad &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad &{}\quad c_{n,n-1} \\ 0 &{}\quad &{}\quad \dots &{}\quad 0 \end{array} \right] ,\\ B_{n,1}=\left[ \begin{array}{ccccccc} b_{n,0} &{}\quad e_{n,0}&{}\quad &{}\quad &{}\quad \\ d_{n,1} &{}\quad b_{n,1}&{}\quad e_{n,1}&{}\quad &{}\quad \\ &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad d_{n,n-1} &{}\quad b_{n,n-1}&{}\quad e_{n,n-1} \\ &{}\quad &{}\quad &{}\quad d_{n,n} &{}\quad b_{n,n} \end{array} \right] , \end{array} \end{aligned}$$
(2.6)

where the entries of \(A_{n,1}\), \(B_{n,1}\) and \(C_{n,1}\) [see (2.4)] are given by

$$\begin{aligned} \begin{aligned} a_{n,k}&=\frac{1}{2}a_{n+\gamma +\frac{1}{2}}\frac{\delta _{n+2\gamma +1,k}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 0, \quad k=0,1,\ldots ,n, \\ c_{n,k}&=\frac{1}{2}c_{n+\gamma +\frac{1}{2}}\frac{\delta _{n,k}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1, \quad k=0,1,\ldots ,n-1,\\ e_{n,k}&=\frac{1}{2}a_{k}\frac{\delta _{n+\gamma +\frac{1}{2},k+\gamma +\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1, \quad k=0,1,\ldots ,n-1,\\ d_{n,k}&=\frac{1}{2}c_{k}\frac{\delta _{n+\gamma +\frac{1}{2},k-\gamma -\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1, \quad k=1,2,\ldots ,n,\\ b_{n,k}&=1-a_{n,k}-c_{n,k}-d_{n,k}-e_{n,k},\quad n\ge 0, \quad k=0,1,\ldots ,n. \end{aligned} \end{aligned}$$
(2.7)

Remark 2.1

The coefficient \(b_{n,k}\) can also be written as

$$\begin{aligned} b_{n,k}=\frac{1}{2}(b_{n+\gamma +1/2}+b_k)+\frac{1-4\gamma ^2}{4(\beta ^2-\alpha ^2)}(2b_{n+\gamma +1/2}-1)(2b_{k}-1). \end{aligned}$$
(2.8)

Also, from (2.4) it is possible to see that

$$\begin{aligned} \frac{\delta _{n+2\gamma +1,k}}{\delta _{n+\gamma +\frac{1}{2},k}}+\frac{\delta _{n,k}}{\delta _{n+\gamma +\frac{1}{2},k}}+\frac{\delta _{n+\gamma +\frac{1}{2},k-\gamma -\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}}+\frac{\delta _{n+\gamma +\frac{1}{2},k+\gamma +\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}}=4. \end{aligned}$$

Remark 2.2

Observe, from (2.4) and recalling that \(\alpha ,\beta ,\gamma >-1\), \(\alpha +\gamma +3/2>0\) and \(\beta +\gamma +3/2>0\), that the coefficients \(a_{n,k}, c_{n,k}, e_{n,k}, d_{n,k}\) in (2.7) are all positive. The coefficient \(a_{n,k}\) may have a different sign if \(n=k=0\) or \(k=n\). But a direct computation looking at the factors of \(a_{n,k}\) shows that \(a_{n,n}>0,n\ge 0\). The same applies to \(e_{n,k}\) for \(k=0\) and \(d_{n,k}\) for \(k=n\). As for the coefficient \(b_{n,k}\) in (2.7) [see also (2.8)] it is also possible to see that \(b_{n,k}\ge 0\). Indeed, let us call \(T_{n,k}^{\alpha ,\beta ,\gamma }\) the right-hand term in (2.8). A direct computation using (2.4) gives that

$$\begin{aligned} T_{n,k}^{\alpha ,\beta ,\gamma }=\frac{(1-4\gamma ^2)(\beta ^2-\alpha ^2)}{4(2k+\alpha +\beta )(2k+\alpha +\beta +2)(2n+\alpha +\beta +2\gamma +1)(2n+\alpha +\beta +2\gamma +3)}. \end{aligned}$$

If \(\alpha ^2=\beta ^2\) then \(b_{n,k}=1/2\) and there is nothing to prove. For \(n=0\) and \(k=0\), we have

$$\begin{aligned} b_{0,0}=\frac{2\beta +2\gamma +3}{2(\alpha +\beta +2\gamma +3)}>0. \end{aligned}$$

Therefore, it is enough to prove that \(b_{n,k}\ge 0\) for \(n\ge 1\) and \(k=0,1,\ldots ,n\). On one side, if \(T_{n,k}^{\alpha ,\beta ,\gamma }\ge 0\), then we immediately have that \(b_{n,k}\ge 0\) since \(b_n\ge 0\). On the contrary, for the values of \(\alpha ,\beta ,\gamma \) such that \(T_{n,k}^{\alpha ,\beta ,\gamma }<0\), we always have that

$$\begin{aligned} \frac{\partial }{\partial n} T_{n,k}^{\alpha ,\beta ,\gamma }=-\frac{(1-4\gamma ^2)(\beta ^2-\alpha ^2)(2n+\alpha +\beta +2\gamma +2)}{4(2k+\alpha +\beta )(2k+\alpha +\beta +2)(2n+\alpha +\beta +2\gamma +1)^2(2n+\alpha +\beta +2\gamma +3)^2}. \end{aligned}$$

It is easy to see that \(T_{n,k}^{\alpha ,\beta ,\gamma }<0\) if and only if either \(1<4\gamma ^2\) and \(\beta ^2>\alpha ^2\) or \(1>4\gamma ^2\) and \(\beta ^2<\alpha ^2\). For these values of \(\alpha ,\beta ,\gamma \), we have that \(T_{n,k}^{\alpha ,\beta ,\gamma }\), as a function of \(n\ge 1\), is always increasing and converging to 0. Hence, the minimum negative value of \(T_{n,k}^{\alpha ,\beta ,\gamma }\) is attained at \(n=1\). Therefore, it is enough to check that \(b_{1,1}\ge 0\) and \(b_{1,0}\ge 0\). On one hand, \(b_{1,1}\) can be written as

$$\begin{aligned} b_{1,1}=\frac{(\beta +1)(2\beta +2\gamma +3)+(\alpha +1)(2\beta +2\gamma +7)}{2(\alpha +\beta +2)(\alpha +\beta +2\gamma +5)}>0. \end{aligned}$$

On the other hand, \(b_{1,0}\) can be written as

$$\begin{aligned} b_{1,0}=\frac{(2\beta +2\gamma +5)(\alpha +\beta +2)^2+(2\gamma +1)(\alpha +\beta +2)(2\beta +2\gamma +3)+4(\beta +1)(2\gamma +1)}{2(\alpha +\beta +2)(\alpha +\beta +2\gamma +3)(\alpha +\beta +2\gamma +5)}. \end{aligned}$$

The previous expression is always nonnegative if \(\gamma \ge -1/2\). If \(-1<\gamma <-1/2\), then we must have that \(\alpha ,\beta >-1/2\). In this situation, we can rewrite the numerator of \(b_{1,0}\) as the following expression:

$$\begin{aligned}&(2\beta +2\gamma +5)(\alpha +1/2)^2+(\alpha +1/2)\left[ 4(\beta +1/2)^2+(8\gamma +14)(\beta +1/2)+4\gamma ^2+10\gamma +10\right] \\ {}&\quad +2(\beta +1/2)^3+(6\gamma +10)(\beta +1/2)^2+2(\gamma +1)(2\gamma +9)(\beta +1/2)+4(\gamma +1)(\gamma +2), \end{aligned}$$

which is always positive. Finally, we need to see that \(b_{n,k}, n\ge 1,\) does not have any minimum for some value of n or k. A direct computation shows that

$$\begin{aligned} \frac{\partial }{\partial n}b_{n,k}=\frac{(\alpha ^2-\beta ^2)((2k+\alpha +\beta +1)^2-4\gamma ^2)(2n+\alpha +\beta +2\gamma +2)}{(2k+\alpha +\beta )(2k+\alpha +\beta +2)(2n+\alpha +\beta +2\gamma +1)^2(2n+\alpha +\beta +2\gamma +3)^2}. \end{aligned}$$

Hence the sign of \(\frac{\partial }{\partial n}b_{n,k}\) is always the same and \(b_{n,k}\) is always strictly increasing or decreasing to the value of \(\frac{1}{2}b_k\ge 0\). For the possible values of nonnegative integer values of k such that \((2k+\alpha +\beta +1)^2=4\gamma ^2\), i.e., \(k_0^{\pm }=-\frac{1}{2}\left( \alpha +\beta \pm 2\gamma +1\right) \), we always have that \(\frac{\partial }{\partial n}b_{n,k_0^{\pm }}=0\). A direct computation shows that

$$\begin{aligned} b_{n,k_0^{\pm }}=\frac{1}{2}+\frac{\alpha ^2-\beta ^2}{4(1-4\gamma ^2)}. \end{aligned}$$

Observe that \(b_{n,k_0^{\pm }}\) is a constant positive value under the conditions on the parameters \(\alpha ,\beta ,\gamma \) such that \(T_{n,k}^{\alpha ,\beta ,\gamma }<0\). Therefore, \(b_{n,k}\ge 0\) for \(n\ge 0\) and \(k=0,1,\ldots ,n\). In particular, that means that the block tridiagonal Jacobi matrix constructed from the coefficients in (2.6) [see \(J_1\) in (3.1) below] is a stochastic matrix.

On the other side, the matrices \(A_{n,2}\), \(B_{n,2}\) and \(C_{n,2}\) are tridiagonal matrices of the form

$$\begin{aligned} \begin{array}{c} A_{n,2}=\left[ \begin{array}{ccccccc} a_{n,0}^{(2)} &{}\quad a_{n,0}^{(3)}&{}\quad &{}\quad &{}\quad \\ a_{n,1}^{(1)} &{}\quad a_{n,1}^{(2)}&{}\quad a_{n,1}^{(3)} &{}\quad &{}\quad \\ &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad a_{n,n}^{(1)} &{}\quad a_{n,n}^{(2)} &{}\quad a_{n,n}^{(3)} \end{array} \right] , \\ \quad B_{n,2}=\left[ \begin{array}{ccccccc} b_{n,0}^{(2)} &{}\quad b_{n,0}^{(3)}&{}\quad &{}\quad &{}\quad \\ b_{n,1}^{(1)} &{}\quad b_{n,1}^{(2)}&{}\quad b_{n,1}^{(3)} &{}\quad &{}\quad \\ &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad b_{n,n-1}^{(1)} &{}\quad b_{n,n-1}^{(2)} &{}\quad b_{n,n-1}^{(3)} \\ &{}\quad &{}\quad &{}\quad b_{n,n}^{(1)} &{}\quad b_{n,n}^{(2)} \end{array} \right] , \\ C_{n,2}=\left[ \begin{array}{ccccccc} c_{n,0}^{(2)} &{}\quad c_{n,0}^{(3)}&{}\quad &{}\quad &{}\quad \\ c_{n,1}^{(1)} &{}\quad c_{n,1}^{(2)}&{}\quad c_{n,1}^{(3)} &{}\quad &{}\quad \\ &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad c_{n,n-2}^{(1)} &{}\quad c_{n,n-2}^{(2)} &{}\quad c_{n,n-2}^{(3)} \\ &{}\quad &{}\quad &{}\quad c_{n,n-1}^{(1)} &{}\quad c_{n,n-1}^{(2)} \\ &{}\quad &{}\quad &{}\quad &{}\quad c_{n,n}^{(1)} \end{array} \right] , \end{array} \end{aligned}$$
(2.9)

where the entries of \(A_{n,2}\), \(B_{n,2}\) and \(C_{n,2}\) [see again (2.4)] are given by

$$\begin{aligned} \begin{aligned} a_{n,k}^{(1)}&=2c_ka_{n+\gamma +\frac{1}{2}}\frac{\delta _{n+\gamma +\frac{1}{2},k-\gamma -\frac{1}{2}}\delta _{n+\gamma +1,k-\gamma -1}}{\delta _{n+\frac{1}{2}(\gamma +\frac{1}{2}),k-\frac{1}{2}(\gamma +\frac{1}{2})}\delta _{n+\frac{1}{2}(\gamma +\frac{3}{2}),k-\frac{1}{2}(\gamma +\frac{3}{2})}},\quad n\ge 1,\quad k=1,\ldots ,n,\\ a_{n,k}^{(2)}&=(2b_k-1)a_{n+\gamma +\frac{1}{2}}\frac{\delta _{n+2\gamma +1,k}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 0, \quad k=0,1,\ldots ,n,\\ a_{n,k}^{(3)}&=2a_ka_{n+\gamma +\frac{1}{2}}\frac{\delta _{n+\gamma +\frac{1}{2},k+\gamma +\frac{1}{2}}\delta _{n+\gamma +1,k+\gamma +1}}{\delta _{n+\frac{1}{2}(\gamma +\frac{1}{2}),k+\frac{1}{2}(\gamma +\frac{1}{2})}\delta _{n+\frac{1}{2}(\gamma +\frac{3}{2}),k+\frac{1}{2}(\gamma +\frac{3}{2})}},\quad n\ge 0,\quad k=0,1,\ldots ,n, \\ b_{n,k}^{(1)}&=(2b_{n+\gamma +\frac{1}{2}}-1)c_k\frac{\delta _{n+\gamma +\frac{1}{2},k-\gamma -\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1, \quad k=1,\ldots ,n, \\ b_{n,k}^{(2)}&=1-b_{n,k}^{(1)}-b_{n,k}^{(3)}-a_{n,k}^{(1)}-a_{n,k}^{(2)}-a_{n,k}^{(3)}-c_{n,k}^{(1)}-c_{n,k}^{(2)}-c_{n,k}^{(3)},\quad n\ge 0, \quad k=0,1,\ldots ,n, \\ b_{n,k}^{(3)}&=(2b_{n+\gamma +\frac{1}{2}}-1)a_k\frac{\delta _{n+\gamma +\frac{1}{2},k+\gamma +\frac{1}{2}}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1, \quad k=0,1,\ldots ,n-1, \\ c_{n,k}^{(1)}&=2c_kc_{n+\gamma +\frac{1}{2}}\frac{\delta _{n,k}\delta _{n-\frac{1}{2},k-\frac{1}{2}}}{\delta _{n+\frac{1}{2}(\gamma -\frac{1}{2}),k+\frac{1}{2}(\gamma -\frac{1}{2})}\delta _{n+\frac{1}{2}(\gamma +\frac{1}{2}),k+\frac{1}{2}(\gamma +\frac{1}{2})}},\quad n\ge 1, \quad k=1,\ldots ,n, \\ c_{n,k}^{(2)}&=(2b_k-1)c_{n+\gamma +\frac{1}{2}}\frac{\delta _{n,k}}{\delta _{n+\gamma +\frac{1}{2},k}},\quad n\ge 1,\quad k=0,1,\ldots ,n-1,\\ c_{n,k}^{(3)}&=2a_kc_{n+\gamma +\frac{1}{2}}\frac{\delta _{n,k}\delta _{n-\frac{1}{2},k+\frac{1}{2}}}{\delta _{n+\frac{1}{2}(\gamma -\frac{1}{2}),k-\frac{1}{2}(\gamma -\frac{1}{2})}\delta _{n+\frac{1}{2}(\gamma +\frac{1}{2}),k-\frac{1}{2}(\gamma +\frac{1}{2})}},\quad n\ge 2,\quad k=0,1,\ldots ,n-2. \end{aligned} \nonumber \\ \end{aligned}$$
(2.10)

Remark 2.3

Observe again, from (2.4) and recalling that \(\alpha ,\beta ,\gamma >-1\), \(\alpha +\gamma +3/2>0\) and \(\beta +\gamma +3/2>0\), that the coefficients \(a_{n,k}^{(1)}, a_{n,k}^{(3)}, c_{n,k}^{(1)}, c_{n,k}^{(3)}\) in (2.10) are all positive. The coefficient \(a_{n,k}^{(3)}\) may have a different sign if \(k=n=0\) or \(k=n\). But a direct computation looking at the factors of \(a_{n,n}^{(3)}\) shows that \(a_{n,n}^{(3)}>0,n\ge 0\). The same applies to \(a_{n,k}^{(1)}\) for \(k=n\). Unlike the coefficient \(b_{n,k}\) in (2.8), we were not able to find a simplified expression of \(b_{n,k}^{(2)}\).

The coefficients for \(A_{n,i}, B_{n,i},C_{n,i}, i=1,2,\) can be significantly simplified for the values of \(\gamma =\pm 1/2\). For \(\gamma =-1/2\), we get

$$\begin{aligned} a_{n,k}=\frac{1}{2}a_n,\quad c_{n,k}=\frac{1}{2}c_n,\quad e_{n,k}=\frac{1}{2}a_k,\quad d_{n,k}=\frac{1}{2}c_k,\quad b_{n,k}=\frac{1}{2}(b_n+b_k), \end{aligned}$$

and

$$\begin{aligned} a_{n,k}^{(1)}&=2a_nc_k,\quad a_{n,k}^{(2)}=a_n(2b_k-1),\quad a_{n,k}^{(3)}=2a_na_k,\\ b_{n,k}^{(1)}&=(2b_n-1)c_k,\quad b_{n,k}^{(2)}=\frac{1}{2}\left( 1+(2b_n-1)(2b_k-1)\right) ,\quad b_{n,k}^{(3)}=(2b_n-1)a_k,\\ c_{n,k}^{(1)}&=2c_nc_k,\quad c_{n,k}^{(2)}=c_n(2b_k-1),\quad c_{n,k}^{(3)}=2c_na_k, \end{aligned}$$

while for \(\gamma =1/2\), we obtain

$$\begin{aligned} a_{n,k}&=\frac{1}{2}a_{n+1}\frac{\delta _{n+2,k}}{\delta _{n+1,k}},\quad c_{n,k}=\frac{1}{2}c_{n+1}\frac{\delta _{n,k}}{\delta _{n+1,k}},\quad e_{n,k}=\frac{1}{2}a_k\frac{\delta _{n+1,k+1}}{\delta _{n+1,k}},\\ d_{n,k}&=\frac{1}{2}c_k\frac{\delta _{n+1,k-1}}{\delta _{n+1,k}},\quad b_{n,k}=\frac{1}{2}(b_{n+1}+b_k), \end{aligned}$$

and

$$\begin{aligned} a_{n,k}^{(1)}&=2a_{n+1}c_k\frac{\delta _{n+\frac{3}{2},k-\frac{3}{2}}}{\delta _{n+\frac{1}{2},k-\frac{1}{2}}},\quad a_{n,k}^{(2)}=a_{n+1}(2b_k-1)\frac{\delta _{n+2,k}}{\delta _{n+1,k}},\quad a_{n,k}^{(3)}=2a_{n+1}a_k\frac{\delta _{n+\frac{3}{2},k+\frac{3}{2}}}{\delta _{n+\frac{1}{2},k+\frac{1}{2}}},\\ b_{n,k}^{(1)}&=(2b_{n+1}-1)c_k\frac{\delta _{n+1,k-1}}{\delta _{n+1,k}},\quad b_{n,k}^{(2)}=\frac{1}{2}\left( 1+(2b_{n+1}-1)(2b_k-1)\right) ,\\ {}&b_{n,k}^{(3)}=(2b_{n+1}-1)a_k\frac{\delta _{n+1,k+1}}{\delta _{n+1,k}},\\ c_{n,k}^{(1)}&=2c_{n+1}c_k\frac{\delta _{n-\frac{1}{2},k-\frac{1}{2}}}{\delta _{n+\frac{1}{2},k+\frac{1}{2}}},\quad c_{n,k}^{(2)}=c_{n+1}(2b_k-1)\frac{\delta _{n,k}}{\delta _{n+1,k}},\quad c_{n,k}^{(3)}=2c_{n+1}a_k\frac{\delta _{n-\frac{1}{2},k+\frac{1}{2}}}{\delta _{n+\frac{1}{2},k-\frac{1}{2}}}. \end{aligned}$$

Remark 2.4

The normalization of the polynomials \(Q_{n,k}^{\alpha ,\beta ,\gamma }(u,v)\) such that \(Q_{n,k}^{\alpha ,\beta ,\gamma }(1,1)=1\) will guarantee us that the sum of all rows of the corresponding Jacobi matrices \(J_1\) and \(J_2\) [see (3.1) below] is exactly 1. This does not mean that both \(J_1\) and \(J_2\) are stochastic matrices or have some probabilistic interpretation, something that we will discuss in the next section. We could have used another “corner” of the region \(\Omega \) (see Fig. 1) like (0, 1) or (1/2, 0). On one side, it turns out that normalization at the point (1/2, 0) will not provide us Jacobi matrices with probabilistic interpretation. On the other side, normalization at the point (0, 1) is somehow “symmetric” to the normalization at the point (1, 1). Indeed, we have that \(P_{n,k}^{\alpha ,\beta ,\gamma }(0,1)=(-1)^{n+k}\sigma _{n,k}\), where \(\sigma _{n,k}\) is given by (2.2), and the corresponding new vector polynomials \({\widetilde{{\mathbb {Q}}}}_n\) satisfy the three-term recurrence relations

$$\begin{aligned} \begin{aligned} -u\, \widetilde{{\mathbb {Q}}}_n(u,v)&= {\widetilde{A}}_{n,1} \widetilde{{\mathbb {Q}}}_{n+1}(u,v)+ {\widetilde{B}}_{n,1} \widetilde{{\mathbb {Q}}}_n(u,v) + {\widetilde{C}}_{n,1} \widetilde{{\mathbb {Q}}}_{n-1}(u,v), \\ v\,\widetilde{{\mathbb {Q}}}_n(u,v)&= {\widetilde{A}}_{n,2} \widetilde{{\mathbb {Q}}}_{n+1}(u,v)+ {\widetilde{B}}_{n,2} \widetilde{{\mathbb {Q}}}_n(u,v) + {\widetilde{C}}_{n,2} \widetilde{{\mathbb {Q}}}_{n-1}(u,v), \end{aligned} \end{aligned}$$

where the coefficients \({\widetilde{A}}_{n,i}, {\widetilde{B}}_{n,i}, {\widetilde{C}}_{n,i}, i=1,2,\) are exactly the same as the coefficients \(A_{n,i}, B_{n,i},C_{n,i}, i=1,2,\) but changing \(\alpha \) by \(\beta \) and \(\beta \) by \(\alpha \), except for \({\widetilde{B}}_{n,1}\) where we have \({\widetilde{B}}_{n,1}=B_{n,1}-I\) (changing \(\alpha \) by \(\beta \) again and vice versa). In this case, we have that the sum of the rows of the Jacobi matrix \(J_1\) is 0, while the sum of the rows of the Jacobi matrix \(J_2\) is 1. For more comments about the choice of normalizing corners the reader can consult [6, Section 6].

3 QBD Processes Associated with Jacobi–Koornwinder Bivariate Polynomials

In this section, we will study under what conditions we may provide a probabilistic interpretation of the coefficients of the three-term recurrence relations (2.6) and (2.9). From the recurrence relations (2.3), we can define the following two block tridiagonal Jacobi matrices:

$$\begin{aligned} J_{1}=\left( \begin{array}{cccccc} B_{0,1} &{}\quad A_{0,1} &{}\quad &{}\quad &{}\bigcirc \\ C_{1,1} &{}B_{1,1} &{}\quad A_{1,1} &{}\quad &{}\quad \\ &{}C_{2,1} &{}\quad B_{2,1} &{}A_{2,1} &{}\quad \\ \bigcirc &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots \end{array}\right) ,\quad J_{2}=\left( \begin{array}{cccccc} B_{0,2} &{}\quad A_{0,2} &{}\quad &{}\quad &{}\bigcirc \\ C_{1,2} &{}B_{1,2} &{}\quad A_{1,2} &{}\quad &{}\quad \\ &{}C_{2,2} &{}\quad B_{2,2} &{}A_{2,2} &{}\quad \\ \bigcirc &{}\quad &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots \end{array}\right) .\nonumber \\ \end{aligned}$$
(3.1)

By construction of the vector polynomials \(\mathbb {Q}_n\) [see (2.3)], we always have that \(J_i\varvec{e}=\varvec{e}, i=1,2\), where \(\varvec{e}=\left( 1, 1, 1, \ldots \right) ^T\) is the semi-infinite vector with all components equal to 1. We now consider the linear convex combination of \(J_1\) and \(J_2\) in the following way:

$$\begin{aligned} \varvec{P}=(1-\tau )J_1+\tau J_2,\quad 0\le \tau \le 1. \end{aligned}$$
(3.2)

We would like to see under what conditions we get a probabilistic interpretation of \(\varvec{P}\). In particular, we will see when \(\varvec{P}\) is a stochastic matrix. We immediately have that \(\varvec{P}\varvec{e}=\varvec{e}\) but now we need all entries of \(\varvec{P}\) to be positive (except possibly for the main block diagonal, where we only need to be nonnegative). Therefore, looking at the nonzero entries of \(\varvec{P}\), we need to have

$$\begin{aligned} \begin{aligned} \tau a_{n,k}^{(1)}&> {} 0,&\quad n\ge 1,&\quad k=1,\ldots ,n, \\ (1-\tau )a_{n,k}+\tau a_{n,k}^{(2)}&> {} 0,&\quad n\ge 0,&\quad k=0,1,\ldots ,n, \\ \tau a_{n,k}^{(3)}&> {} 0,&\quad n\ge 0,&\quad k=0,1,\ldots ,n, \\ (1-\tau )d_{n,k}+\tau b_{n,k}^{(1)}&\ge 0,&\quad n\ge 1,&\quad k=1,\ldots ,n, \\ (1-\tau )b_{n,k}+\tau b_{n,k}^{(2)}&\ge 0,&\quad n\ge 0,&\quad k=0,1,\ldots ,n,\\ (1-\tau )e_{n,k}+\tau b_{n,k}^{(3)}&\ge 0,&\quad n\ge 1,&\quad k=0,1,\ldots ,n-1,\\ \tau c_{n,k}^{(1)}&> {} 0,&\quad n\ge 1,&\quad k=1,\ldots ,n, \\ (1-\tau )c_{n,k}+\tau c_{n,k}^{(2)}&> {} 0,&\quad n\ge 1,&\quad k=0,1,\ldots ,n-1, \\ \tau c_{n,k}^{(3)}&> {} 0,&\quad n\ge 2,&\quad k=0,1,\ldots ,n-2. \end{aligned} \end{aligned}$$
(3.3)

First, from Remark 2.2, we have that the case \(\tau =0\) (\(\varvec{P}=J_1\)) gives rise to a stochastic matrix, so it is enough to study the inequalities in (3.3) for \(0<\tau \le 1\). Also, from Remark 2.3, we have that the first, third, seventh and ninth inequalities in (3.3) hold for any \(0<\tau \le 1\). From the definition of the coefficients (2.7) and (2.10), we have the following properties:

$$\begin{aligned} a_{n,k}-a_{n,k}^{(2)}&=a_{n,k}(3-4b_k),\quad c_{n,k}-c_{n,k}^{(2)}=c_{n,k}(3-4b_k),\\ d_{n,k}-b_{n,k}^{(1)}&=d_{n,k}(3-4b_{n+\gamma +1/2}),\quad e_{n,k}-b_{n,k}^{(3)}=e_{n,k}(3-4b_{n+\gamma +1/2}), \end{aligned}$$

where \(b_n\) is defined by (2.4). These properties can be used in the second, fourth, sixth and eighth inequalities in (3.3) in the following way:

$$\begin{aligned} \begin{aligned} (1-\tau )a_{n,k}+\tau a_{n,k}^{(2)}&=a_{n,k}(1-\tau (3-4b_k))>0,\quad n\ge 0,\quad k=0,1,\ldots ,n,\\ (1-\tau )d_{n,k}+\tau b_{n,k}^{(1)}&=d_{n,k}(1-\tau (3-4b_{n+\gamma +1/2}))>0,\quad n\ge 1,\quad k=1,\ldots ,n,\\ (1-\tau )e_{n,k}+\tau b_{n,k}^{(3)}&=e_{n,k}(1-\tau (3-4b_{n+\gamma +1/2}))>0,\quad n\ge 1,\quad k=0,1,\ldots ,n-1,\\ (1-\tau )c_{n,k}+\tau c_{n,k}^{(2)}&=c_{n,k}(1-\tau (3-4b_k))>0,\quad n\ge 1,\quad k=0,1,\ldots ,n-1. \end{aligned} \end{aligned}$$
(3.4)

Since the coefficients \(a_{n,k}, c_{n,k}, e_{n,k}, d_{n,k}\) are all positive (see Remark 2.2), all the analysis is reduced to study the inequality \((1-\tau (3-4b_k))>0\) for \(k\in \mathbb {N}_0\) (the rest of cases are just simple translations of the previous inequality (see below)).

Let us start with the first and fourth inequalities in (3.4), i.e., \(1-\tau (3-4b_k)>0\) for \(k\in \mathbb {N}_0\). If there is some \(k\in \mathbb {N}_0\) such that \(3-4b_k\le 0\), then the inequality \(1-\tau (3-4b_k)>0\) always holds, so we can take any \(\tau \) such that \(0<\tau \le 1\). Otherwise, if \(3-4b_k>0\), then \(1-\tau (3-4b_k)>0\) holds if and only if

$$\begin{aligned} 0<\tau \le \inf _{\{k\in \mathbb {N}_0:3-4b_k>0\}}\left( \frac{1}{3-4b_k}\right) . \end{aligned}$$
(3.5)

Therefore, let us call

$$\begin{aligned} F_k^{\alpha ,\beta }=\frac{1}{3-4b_k}=\frac{(2k+\alpha +\beta )(2k+\alpha +\beta +2)}{4k(k+\alpha +\beta +1)+(\alpha +\beta )(3\alpha -\beta +2)}=\frac{N_k^{\alpha ,\beta }}{D_k^{\alpha ,\beta }}.\nonumber \\ \end{aligned}$$
(3.6)

We have the following properties:

  • If \(\alpha ^2=\beta ^2\) then \(F_k^{\alpha ,\beta }=1\).

  • \(\lim _{k\rightarrow \infty }F_k^{\alpha ,\beta }=1\).

  • \(F_k^{\alpha ,\beta }, k\ge 0,\) only vanishes at \(k_0=-\frac{1}{2}(\alpha +\beta )\in (-\infty ,1)\) and \(D_{k_0}^{\alpha ,\beta }=2(\alpha ^2-\beta ^2)\).

  • The values of \(F_k^{\alpha ,\beta }\) at \(k=0\) and \(k=1\) are given by

    $$\begin{aligned} F_0^{\alpha ,\beta }=\frac{\alpha +\beta +2}{3\alpha -\beta +2},\quad F_1^{\alpha ,\beta }=\frac{(\alpha +\beta +2)(\alpha +\beta +4)}{4(\alpha +\beta +2)+(\alpha +\beta )(3\alpha -\beta +2)}. \end{aligned}$$
    (3.7)

  • The derivative of \(F_k^{\alpha ,\beta }\) with respect to k is given by

    $$\begin{aligned} \frac{\partial }{\partial k}F_k^{\alpha ,\beta }=\frac{8(\alpha ^2- \beta ^2)(2k+\alpha +\beta +1)}{\left( D_k^{\alpha ,\beta }\right) ^2}. \end{aligned}$$
    (3.8)

The behavior of \(F_k^{\alpha ,\beta }, k\ge 0,\) depends strongly on the behavior of the denominator \(D_k^{\alpha ,\beta }\). In the following lemma, we study the properties of \(D_k^{\alpha ,\beta }\).

Lemma 3.1

As a function of k, \(D_k^{\alpha ,\beta }\) is a parabola with vertex located at \(k_0=-\frac{1}{2}(\alpha +\beta +1)\in (-\infty ,1/2)\) and \(D_{k_0}^{\alpha ,\beta }=2\alpha ^2-2\beta ^2-1\). If \(1+2\beta ^2-2\alpha ^2\ge 0\), then the real zeros of \(D_k^{\alpha ,\beta }\) are given by

$$\begin{aligned} x_0^{\pm }=-\frac{1}{2}\left( \alpha +\beta +1\mp \sqrt{1+2\beta ^2-2\alpha ^2}\right) . \end{aligned}$$
(3.9)

The location of the zeros \(x_0^{\pm }\) depends on the values of \(\alpha ,\beta \). For \(x_0^{+}\), we have

  • If \(\alpha +\beta +1>0\) then \(x_0^{+}>0\) if and only if \((\alpha +\beta )(3\alpha -\beta +2)<0\) and \(x_0^{+}<0\) otherwise.

  • If \(\alpha +\beta +1<0\) then \(x_0^{+}>0\).

For \(x_0^{-}\), we have

  • \(x_0^{-}<1/2\).

  • If \(\alpha +\beta +1>0\) then \(x_0^{-}<0\).

  • If \(\alpha +\beta +1<0\) then \(x_0^{-}>0\) if and only if \((\alpha +\beta )(3\alpha -\beta +2)>0\) and \(x_0^{-}<0\) otherwise.

Proof

It is a standard exercise of calculus, so we omit the details. \(\square \)

From the previous analysis, we see that the graph of \(F_k^{\alpha ,\beta }\) depends strongly on the values of the parameters \(\alpha \) and \(\beta \). In particular, the monotonicity [see (3.8)] depends on the sign of \(\alpha ^2-\beta ^2\) and \(\alpha +\beta +1\), the existence of the zeros of the parabola \(D_k^{\alpha ,\beta }\) depends on the sign of the hyperbola \(1+2\beta ^2-2\alpha ^2\) [see (3.9)] and their location also on the sign of \(3\alpha -\beta +2\) and \(\alpha +\beta \) (see the end of Lemma 3.1). Therefore, taking in mind all these curves, let us consider a subdivision of the two-dimensional region \(\alpha ,\beta >-1\) given by Fig. 2.

Fig. 2
figure 2

Subdivision of the two-dimensional region \(\alpha ,\beta >-1\) where \(F_k^{\alpha ,\beta }\) may have different graphs

Full size image

Next, we will study the graph of \(F_k^{\alpha ,\beta }\) at each of these regions and determine the upper bound of \(\tau \) in (3.5). In the following analysis, we omit the cases where \(\alpha ^2=\beta ^2\) and hence \(F_k^{\alpha ,\beta }=1\) [see the first property after (3.6)].

  • \(A_1=\{\alpha <\beta ,\beta \le 3\alpha +2,\alpha >-\beta \}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(k\ge 0\) to the limit value of 1, \(F_0^{\alpha ,\beta }>1\) and \(x_0^+\le 0\) (see Lemma 3.1). Therefore, the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is attached at \(k\rightarrow \infty \) and \(0<\tau \le 1\) (see the first graph of Fig. 3). If \(\beta =3\alpha +2\), then \(F_{0^+}^{\alpha ,\beta }\rightarrow +\infty \) and we still have \(0<\tau \le 1\).

  • \(A_2=\{\alpha>-1,\beta>3\alpha +2,\alpha >-\beta \}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(k\ge 0\) to the limit value of 1, \(F_0^{\alpha ,\beta }<0\), \(x_0^-<0\) and \(x_0^+>0\) (see Lemma 3.1). The graph of \(F_k^{\alpha ,\beta }\) has a pole at \(x_0^+\), but since it is strictly decreasing and \(F_0^{\alpha ,\beta }<0\), the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is again attached at \(k\rightarrow \infty \), and therefore, \(0<\tau \le 1\) (see the second graph of Fig. 3).

  • \(B_1=\{\beta >-1,1+2\beta ^2\le 2\alpha ^2\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly increasing for \(k\ge 0\) to the limit value of 1, \(0<F_0^{\alpha ,\beta }<1\) and \(D_k^{\alpha ,\beta }>0\) (there are no poles). Hence, the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is attached at \(k=0\), and therefore, \(0<\tau \le F_0^{\alpha ,\beta }\) (see the third graph of Fig. 3). If \(1+2\beta ^2=2\alpha ^2\), then we have the same graph since \(x_0^{\pm }=-\frac{1}{2}(\alpha +\beta +1)<0\) and we still have \(0<\tau \le F_0^{\alpha ,\beta }\).

  • \(B_2=\{\alpha>\beta ,\alpha>-\beta ,1+2\beta ^2>2\alpha ^2\}\). In this region we have that \(F_k^{\alpha ,\beta }\) is strictly increasing for \(k\ge 0\) to the limit value of 1, \(0<F_0^{\alpha ,\beta }<1\) and \(x_0^+<0\) (see Lemma 3.1). The graph of \(F_k^{\alpha ,\beta }\) for \(k\ge 0\) is exactly the same as in region \(B_1\) (see the third graph of Fig. 3). Hence, the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is attached at \(k=0\) and we have \(0<\tau \le F_0^{\alpha ,\beta }\).

  • \(B_3=\{\alpha>\beta ,\alpha <-\beta ,\alpha +\beta>-1,\beta >-1\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(k\ge 0\) to the limit value of 1, \(0<F_0^{\alpha ,\beta }<1\), \(x_0^-<0\) and \(x_0^+>0\) (see Lemma 3.1). The graph of \(F_k^{\alpha ,\beta }\) has a pole at \(x_0^+\), but since it is strictly decreasing and \(0<F_0^{\alpha ,\beta }<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is attached at \(k=0\), and therefore, \(0<\tau \le F_0^{\alpha ,\beta }\) (the graph is similar to the second one in Fig. 3 but with \(0<F_0^{\alpha ,\beta }<1\)).

  • \(B_4=\{\alpha>\beta ,\alpha +\beta \le -1,\beta >-1\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly increasing for \(0\le k< -\frac{1}{2}(\alpha +\beta +1)\) to a local maximum attached at \(k=-\frac{1}{2}(\alpha +\beta +1)\) and then it is strictly decreasing to the limit value of 1, \(0<F_0^{\alpha ,\beta }<1\), \(x_0^-<0\) and \(x_0^+>0\) (see Lemma 3.1). If \(x_0^+<1\), then the positive minimum of \(F_k^{\alpha ,\beta }\) [see (3.5)] is attached at \(k=0\), and therefore, \(0<\tau \le F_0^{\alpha ,\beta }\). If \(x_0^+>1\), we have the same bound since \(F_k^{\alpha ,\beta }<0\) for \(1\le k<x_0^+\) (the graph is similar to the second one in Fig. 3 but with \(0<F_0^{\alpha ,\beta }<1\) and a local maximum at \(k=-\frac{1}{2}(\alpha +\beta +1)\ge 0\)).

  • \(C_1=\{\alpha <-\beta ,1+2\beta ^2\ge 2\alpha ^2,\beta \ge 3\alpha +2\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly increasing for \(k\ge 0\) to the limit value of 1, \(F_0^{\alpha ,\beta }<0\) and \(x_0^+\le 0\) (see Lemma 3.1). The graph of \(F_k^{\alpha ,\beta }\) is similar to the one in region \(B_1\) but \(F_0^{\alpha ,\beta }<0\) (see the third graph in Fig. 3). Since \(F_k^{\alpha ,\beta }\) vanishes at \(-\frac{1}{2}(\alpha +\beta )<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) is now attached at \(k=1\) (recall that \(k\in \mathbb {N}_0\) for the upper bound in (3.5)), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\).

  • \(C_2=\{1+2\beta ^2<2\alpha ^2,\alpha +\beta \ge -1,\alpha >-1\}\). As in region \(C_1\), \(F_k^{\alpha ,\beta }\) is strictly increasing for \(k\ge 0\) to the limit value of 1 and \(F_0^{\alpha ,\beta }<0\) (there are no poles). Since \(F_k^{\alpha ,\beta }\) vanishes at \(-\frac{1}{2}(\alpha +\beta )<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\) (see the third graph in Fig. 3 but with \(F_0^{\alpha ,\beta }<0\)).

  • \(C_3=\{1+2\beta ^2\le 2\alpha ^2,\alpha +\beta <-1,\alpha >-1\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(0\le k< -\frac{1}{2}(\alpha +\beta +1)\) to a local minimum attached at \(k=-\frac{1}{2}(\alpha +\beta +1)\) and then it is strictly increasing to the limit value of 1, \(F_0^{\alpha ,\beta }<0\) (there are no poles if \(1+2\beta ^2<2\alpha ^2\)). Since \(F_k^{\alpha ,\beta }\) vanishes at \(-\frac{1}{2}(\alpha +\beta )<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\) (see the fourth graph in Fig. 3). If \(1+2\beta ^2=2\alpha ^2\), then the local minimum becomes a pole at \(x_0^{\pm }=-\frac{1}{2}(\alpha +\beta +1)\). The positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\).

  • \(C_4=\{1+2\beta ^2>2\alpha ^2,\beta>3\alpha +2,\alpha >-1\}\). This is the only region where both \(x_0^{\pm }>0\). In fact, we have \(x_0 ^-<-\frac{1}{2}(\alpha +\beta +1)<x_0^+<-\frac{1}{2}(\alpha +\beta )<1\). The function \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(0\le k< -\frac{1}{2}(\alpha +\beta +1)\) to a local minimum attached at \(k=-\frac{1}{2}(\alpha +\beta +1)\) (with a pole \(x_0^-\) in the middle) and then it is strictly increasing to the limit value of 1 (with a pole \(x_0^+\) in the middle). Since \(F_0^{\alpha ,\beta }<0\) the positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\) (see the fifth graph in Fig. 3).

  • \(C_5=\{\alpha <\beta ,\beta \le 3\alpha +2,\alpha +\beta \le -1\}\). In this region, we have that \(F_k^{\alpha ,\beta }\) is strictly decreasing for \(0\le k< -\frac{1}{2}(\alpha +\beta +1)\) to a local minimum attached at \(k=-\frac{1}{2}(\alpha +\beta +1)\) and then it is strictly increasing to the limit value of 1 with a pole \(x_0^+<-\frac{1}{2}(\alpha +\beta )<1\) in the middle. Since \(F_0^{\alpha ,\beta }>1\) and \(F_k^{\alpha ,\beta }\) vanishes at \(-\frac{1}{2}(\alpha +\beta )<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\), and therefore, \(0<\tau \le F_1^{\alpha ,\beta }\) (see the sixth graph in Fig. 3). If \(\beta =3\alpha +2\), then \(F_{0^+}^{\alpha ,\beta }\rightarrow +\infty \) and if \(\alpha +\beta =-1\) then \(k=0\) is the local minimum. In both cases, we still have \(0<\tau \le F_1^{\alpha ,\beta }\).

  • \(C_6=\{\alpha<\beta ,\beta<3\alpha +2,\alpha +\beta >-1,\alpha <-\beta \}\). The situation is exactly the same as in region \(C_5\) but now \(F_k^{\alpha ,\beta }\) is strictly increasing for \(k\ge 0\) to the limit value of 1 with a pole \(x_0^+<-\frac{1}{2}(\alpha +\beta )<1\) in the middle. Since \(F_0^{\alpha ,\beta }>1\) and \(F_k^{\alpha ,\beta }\) vanishes at \(-\frac{1}{2}(\alpha +\beta )<1\), the positive minimum of \(F_k^{\alpha ,\beta }\) is again attached at \(k=1\) and therefore \(0<\tau \le F_1^{\alpha ,\beta }\) (see the sixth graph in Fig. 3 but without a local minimum).

Fig. 3
figure 3

The graph of \(F_k^{\alpha ,\beta }, k\ge 0\) for different values of \(\alpha ,\beta \)

Full size image

The previous analysis only applies for the first and fourth inequalities in (3.4), i.e., \(1-\tau (3-4b_k)>0\) for \(k\in \mathbb {N}_0\). For the second and third inequalities in (3.4) (both are the same), we need to have \(1-\tau (3-4b_{n+\gamma +1/2})>0\) for \(n\in \mathbb {N}\), i.e., \(n\ge 1\). Now, we need to study the behavior of \(F_{n+\gamma +1/2}^{\alpha ,\beta }\) for \(n\ge 1\) with a possible dependence on the parameter \(\gamma \). But observe that the graph of \(F_{n+\gamma +1/2}^{\alpha ,\beta }\) is just the graph of \(F_{k}^{\alpha ,\beta }\) displaced to the left the quantity of \(\gamma +1/2\) if \(\gamma +1/2\ge 0\) or to the right the quantity of \(-(\gamma +1/2)\) if \(\gamma +1/2<0\). Since \(\gamma >-1\), we can only displace the graph of \(F_{k}^{\alpha ,\beta }\) at most 1/2 to the right, but any quantity to the left. Going through case by case in the graphs of \(F_{n+\gamma +1/2}^{\alpha ,\beta }\) for every region in Fig. 2, we immediately have that the positive minimum of \(F_{n+\gamma +1/2}^{\alpha ,\beta }\) never improves the upper bound of \(\tau \) in (3.5) except when \(\gamma +1/2<0\) and \(\alpha ,\beta \) are chosen inside any of the regions \(C_i, i=1,\ldots ,6\). Only under these conditions, we have that the positive minimum is attached at \(n=1\), and therefore, \(0<\tau \le F_{\gamma +3/2}^{\alpha ,\beta }\). This is the only situation where the upper bound of \(\tau \) in (3.5) may depend on \(\gamma \).

Finally, we still have to check one more inequality, the one given in the fifth place in (3.3), i.e., \((1-\tau )b_{n,k}+\tau b_{n,k}^{(2)}\ge 0\) for \(n\ge 0\) and \(k=0,1,\ldots ,n\). If \(b_{n,k}\le b_{n,k}^{(2)}\), then we can always take any \(\tau \) such that \(0<\tau \le 1\). But if \(b_{n,k}>b_{n,k}^{(2)}\), then the inequality holds if and only if

$$\begin{aligned} 0<\tau \le \inf _{\{n\in \mathbb {N}_0, 0\le k\le n:b_{n,k}>b_{n,k}^{(2)}\}}\left( \frac{b_{n,k}}{b_{n,k}-b_{n,k}^{(2)}}\right) . \end{aligned}$$
(3.10)

The expression \(b_{n,k}(b_{n,k}-b_{n,k}^{(2)})^{-1}\) is actually intractable. We have not been able to find a suitable expression so that we can analyze \(b_{n,k}(b_{n,k}-b_{n,k}^{(2)})^{-1}\) and compare with \(F_{k}^{\alpha ,\beta }\). Therefore, we have nothing left to do but performing extensive computations on Maple and/or Mathematica in order to get some insight on the values of \(b_{n,k}(b_{n,k}-b_{n,k}^{(2)})^{-1}\) and compare them with \(F_{k}^{\alpha ,\beta }\) or \(F_{n+\gamma +1/2}^{\alpha ,\beta }\). After these computations, checking different values of \(\alpha ,\beta ,\gamma \) in each of the regions of Fig. 2 and for different values of \(n\ge 0\) and \(0\le k\le n\), we have that the positive minimum of \(b_{n,k}(b_{n,k}-b_{n,k}^{(2)})^{-1}\) never improves the upper bound of \(\tau \) in (3.5).

In summary, from Fig. 2, let us call the regions (see now Fig. 4)

$$\begin{aligned} {\textbf {A}}=\bigcup _{i=1}^2 A_i,\quad {\textbf {B}}=\bigcup _{i=1}^4 B_i,\quad {\textbf {C}}=\bigcup _{i=1}^6 C_i. \end{aligned}$$
Fig. 4
figure 4

The regions A, B and C where \(\tau \) may take different values in order to have a stochastic matrix \(\varvec{P}\) (courtesy of C. Juarez)

Full size image

Then, we have the following:

  • In the region A\(=\left\{ \alpha>-1,\beta>\alpha ,\beta >-\alpha \right\} \), we can choose any

    $$\begin{aligned} 0\le \tau \le 1, \end{aligned}$$

    and the matrix \(\varvec{P}\) in (3.2) will always be stochastic.

  • In the region B\(=\left\{ \beta>-1,\alpha >\beta \right\} \), we can choose any

    $$\begin{aligned} 0\le \tau \le F_0^{\alpha ,\beta }, \end{aligned}$$

    where \(F_0^{\alpha ,\beta }\) can be found in (3.7), and the matrix \(\varvec{P}\) in (3.2) will always be stochastic.

  • In the region C\(=\left\{ \alpha>-1,\beta >\alpha ,\beta <-\alpha \right\} \), we can choose any

    $$\begin{aligned} 0\le \tau \le \min \{F_1^{\alpha ,\beta },F_{\gamma +3/2}^{\alpha ,\beta }\}, \end{aligned}$$

    where \(F_1^{\alpha ,\beta }\) can be found in (3.7) and \(F_{\gamma +3/2}^{\alpha ,\beta }\) can be computed from (3.6), and the matrix \(\varvec{P}\) in (3.2) will always be stochastic. A straightforward computation shows that \(F_{\gamma +3/2}^{\alpha ,\beta }\le F_1^{\alpha ,\beta }\) if and only if \(\gamma \le -1/2\).

Therefore, for all values of \(\tau \) in the ranges described above for the regions A, B and C, we have a family of discrete-time QBD processes \(\{Z_t: t = 0, 1, \ldots \}\) with transition probability matrix \(\varvec{P}=(1-\tau )J_1+\tau J_2\). Thus the Karlin–McGregor representation formula (see formula (2.13) of [6]) for the (i, j) block entry of the matrix \(\varvec{P}\) is given by

$$\begin{aligned} \varvec{P}_{i,j}^n=\left( \int _\Omega [(1-\tau )u+\tau v]^n\mathbb {Q}_i(u,v)\mathbb {Q}_j^T(u,v)W_{\alpha ,\beta ,\gamma }(u,v)dudv\right) \Pi _j, \end{aligned}$$
(3.11)

where \(\mathbb {Q}_n,n\ge 0,\) are the vector polynomials satisfying (2.3) and \(W_{\alpha ,\beta ,\gamma }\) is the normalized weight function (2.1). The matrices \(\Pi _j, j\ge 0,\) are the inverses of the norms of the corresponding vector polynomials \(\mathbb {Q}_j, j\ge 0\). Using [6, Lemma 2.1], we have one way of giving an explicit expression of these norms (another way could be using [20, Section 6]). Indeed, it is possible to see that a generalized inverse \(G_n=(G_{n,1},G_{n,2})\) of \(C_n^T=(C_{n,1},C_{n,2})^T\) is given by

$$\begin{aligned} G_n=\left[ \begin{array}{ccccc|ccc} 1/c_{n,0} &{} 0&{}\ldots &{}0&{}0&{}0&{}\ldots &{}0 \\ 0 &{} 1/c_{n,1}&{}\ldots &{}0&{}0&{}0&{}\ldots &{}0 \\ \vdots &{} \vdots &{}\ddots &{}\vdots &{}\vdots &{}\vdots &{}&{}\vdots \\ 0 &{}0&{}\ldots &{}0&{}1/c_{n,n-1}&{}0&{}\ldots &{}0\\ 0 &{}0&{}\ldots &{}-\frac{c_{n,n-2}^{(3)}}{c_{n,n}^{(1)}c_{n,n-2}}&{}-\frac{c_{n,n-1}^{(2)}}{c_{n,n}^{(1)}c_{n,n-1}}&{}0&{}\ldots &{}1/c_{n,n}^{(1)} \end{array}\right] . \end{aligned}$$

Since the representation of \(\Pi _n\) is independent of the choice of the generalized inverse \(G_n\) (see [6, Lemma 2.1]), we have after some straightforward computations, that \(\Pi _n\) is a diagonal matrix of the form \(\Pi _n=\text{ diag }\left( \Pi _{n,0},\ldots ,\Pi _{n,n}\right) \) where

$$\begin{aligned}&\Pi _{n,k} \nonumber \\ =&\frac{(\alpha +1)_k(\alpha +\gamma +3/2)_n(\alpha +\beta +\gamma +5/2)_{n-1}(\alpha +\beta +2\gamma +3)_{n+k-1}(2\gamma +2)_{n-k-1}(\alpha +\beta +2)_{k-1}}{(\beta +1)_k(\beta +\gamma +3/2)_n(\alpha +\beta +2)_{n+k}(\gamma +3/2)_n}\nonumber \\ {}&\quad \times \frac{(2n+\alpha +\beta +2\gamma +2)(n+k+\alpha +\beta +\gamma +3/2)(2n-2k+2\gamma +1)(2k+\alpha +\beta +1)}{k!(n-k)!}. \end{aligned}$$

Another way of writing these norms in terms of the norms of the Jacobi polynomials (2.5) is

$$\begin{aligned} \Pi _{n,k}&=\Vert Q_{n+\gamma +1/2}^{(\beta ,\alpha )}\Vert ^{-2}\Vert Q_{k}^{(\beta ,\alpha )}\Vert ^{-2}\Vert Q_{\gamma +1/2}^{(\beta ,\alpha )}\Vert ^{2}\\&\quad \times \frac{(\alpha +\beta +2\gamma +2)_{n+k}(2\gamma +2)_{n-k-1}(n+k+\alpha +\beta +\gamma +3/2)(2n-2k+2\gamma +1)}{(\alpha +\beta +2)_{n+k}(n-k)!(\alpha +\beta +\gamma +3/2)}. \end{aligned}$$

In particular, we have that the family of polynomials \(\mathbb {Q}_n,n\ge 0,\) is mutually orthogonal. Therefore, another way to write the Karlin–McGregor (3.11) formula is entry by entry

$$\begin{aligned}&\left( \varvec{P}_{i,j}^n\right) _{i',j'}=\mathbb {P}\left[ Z_t=(j,j')\,|\,Z_0=(i,i')) \right] =\frac{\Pi _{j,j'}}{C\sigma _{i,i'}\sigma _{j,j'}}\sum _{k=0}^n\left( {\begin{array}{c}n\\ k\end{array}}\right) (1-\tau )^{n-k}\tau ^k\\ {}&\quad \times \left( \int _\Omega u^{n-k}v^{k}P_{i,i'}^{\alpha ,\beta ,\gamma }(u,v)P_{j,j'}^{\alpha ,\beta ,\gamma }(u,v)(1-2u+v)^\alpha (2u+v-1)^\beta (2u^2-2u-v+1)^\gamma dudv \right) . \end{aligned}$$

According to [6, Theorem 2.5], we can construct an invariant measure \(\varvec{\pi }\) for the QBD process which is given by

$$\begin{aligned} \varvec{\pi }&=\left( \Pi _0;\left( \Pi _1\varvec{e}_2\right) ^T;\left( \Pi _3\varvec{e}_3\right) ^T;\ldots \right) \\ {}&=\left( 1; \frac{(2\alpha +2\gamma +3)(\alpha +\beta +2\gamma +4)(2\alpha +2\beta +2\gamma +5)}{(\alpha +\beta +2)(2\beta +2\gamma +3)},\right. \\ {}&\quad \left. \qquad \frac{(\alpha +1)(2\alpha +2\gamma +3)(2\alpha +2\beta +2\gamma +7)(\alpha +\beta +2\gamma +3)_2}{(\beta +1)(2\beta +2\gamma +3)(2\gamma +3)(\alpha +\beta +2)};\ldots \right) . \end{aligned}$$

Here, \(\varvec{e}_N\) denotes the N-dimensional vector with all components equal to 1. Finally, it is also possible to study recurrence of the family of discrete-time QBD processes using (2.21) of [6]. The process is recurrent if and only if

$$\begin{aligned} \int _\Omega \frac{W_{\alpha ,\beta ,\gamma }(u,v)}{1-\tau v-(1-\tau )u}dudv=\infty . \end{aligned}$$

After some computations it turns out that, in the range of the values of \(\tau \) for which \(\varvec{P}\) is stochastic, this integral is divergent, and therefore, (null) recurrent, if and only if \(\alpha +\gamma \le -1\). Otherwise the QBD process is transient. The QBD process can never be positive recurrent since the spectral measure is absolutely continuous and does not have any jumps (see the end of Section 2 of [6] for more details).

Remark 3.2

Instead of (3.2), we could have considered the situation where \(\varvec{P}=\tau _1J_1+\tau _2J_2\) and \(\varvec{P}\varvec{e}=\varvec{0}\), in which case we would have had a continuous-time QBD process. Since \(J_i\varvec{e}=\varvec{e}, i=1,2\) then we need \(\tau _2=-\tau _1=-\tau \), and therefore, \(\varvec{P}=\tau (J_2-J_1)\). All off-diagonal entries of \(\varvec{P}\) must be nonnegative while the entries of the main diagonal must be nonpositive. A closer look to these conditions entry by entry shows that it is never possible to have a continuous-time QBD process in this context.

Remark 3.3

Going back to Remark 2.4, we could have studied under what conditions we may provide a probabilistic interpretation of a linear combination of \(J_1\) and \(J_2\) of the form \(\varvec{P}=\tau _1J_1+\tau _2J_2\) for the case where we normalize the polynomials at the point (0, 1). For that there are at least two possibilities, either a continuous or a discrete-time QBD process. If we want to have a continuous-time QBD process, then we need \(\varvec{P}\varvec{e} = \varvec{0}\) and nonnegative off-diagonal entries. But this is possible if and only if \(\tau _2=0\) and \(\tau _1>0\), i.e., a positive scalar multiple of \(J_1\). If we want to have a discrete-time QBD process, then we need \(\varvec{P}\varvec{e} = \varvec{e}\) and nonnegative (scalar) entries. This is possible if and only if \(\tau _2=1\) and the parameter \(\tau _1=\tau \) is chosen in such a way that all entries of \(\varvec{P}\) are nonnegative. The entries of \(\varvec{P}=\tau J_1+J_2\) are nonnegative if and only if

$$\begin{aligned} a_{n,k}^{(1)}>0&, \quad \tau a_{n,k}+a_{n,k}^{(2)}>0, \quad a_{n,k}^{(3)}>0,\\ \tau d_{n,k}+b_{n,k}^{(1)}\ge 0&, \quad \tau (b_{n,k}-1)+b_{n,k}^{(2)}\ge 0, \quad \tau e_{n,k}+b_{n,k}^{(3)}\ge 0,\\ c_{n,k}^{(1)}>0&, \quad \tau c_{n,k}+c_{n,k}^{(2)}>0, \quad c_{n,k}^{(3)}>0. \end{aligned}$$

Now, from the definition [see (2.7) and (2.10)], we have

$$\begin{aligned} \frac{a_{n,k}^{(2)}}{a_{n,k}}=\frac{c_{n,k}^{(2)}}{c_{n,k}}=2(2b_k-1),\quad \frac{b_{n,k}^{(1)}}{d_{n,k}}=\frac{b_{n,k}^{(3)}}{e_{n,k}}=2(2b_{n+\gamma +1/2}-1). \end{aligned}$$

Therefore, as before, the lower bounds for \(\tau \) (depending also on \(\alpha ,\beta ,\gamma \)) will depend on the behavior of \(2(1-2b_{k})\) and \(2(1-2b_{n+\gamma +1/2})\). Additionally, the condition \(\tau (b_{n,k}-1)+b_{n,k}^{(2)}\ge 0\) is equivalent to \(\tau \le b_{n,k}^{(2)}/(1-b_{n,k})\), meaning that will also have upper bounds for \(\tau \).

4 An Urn Model for the Jacobi–Koornwinder Bivariate Polynomials

In this section, we will give an urn model associated with one of the QBD models introduced in the previous section. For simplicity, we will study the case of the discrete-time QBD process (3.2) with \(\tau =0\) (therefore, \(\varvec{P}=J_1\)) and \(\beta =\alpha \). In this section, we will assume that \(\alpha \) and \(\gamma \) are nonnegative integers. Consider \(\{Z_t{:}\,t=0,1,\ldots \}\) the discrete-time QBD process on the state space \(\{(n,k){:}\,0\le k\le n,n\in \mathbb {N}_0\}\) whose one-step transition probability matrix is given by the coefficients \(A_{n,1}, B_{n,1}, C_{n,1}\) in (2.6)–(2.7) [see also (2.4)]. At every time step \(t = 0, 1, 2,\ldots \) the state (n, k) will represent the number of n blue balls inside the kth urn \(\hbox {A}_k, k = 0,1,\ldots ,n\). Observe that the number of urns available goes with the number of blue balls at every time step. All the urns we use sit in a bath consisting of an infinite number of blue and red balls.

Since \(\beta =\alpha \) the coefficients in (2.7) are simplified and given explicitly by

$$\begin{aligned} \begin{aligned} a_{n,k}&=\frac{(2n+4\alpha +2\gamma +3)(n-k+2\gamma +1)(n+k+2\alpha +2\gamma + 2)}{4(n+\alpha +\gamma +1)(2n-2k+2\gamma +1)(2n+2k+4\alpha +2\gamma +3)},\quad k=0,1,\ldots ,n, \\ c_{n,k}&=\frac{(2n+2\gamma +1)(n-k)(n+k+2\alpha +1)}{4(n+\alpha +\gamma +1)(2n-2k+2\gamma +1)(2n+2k+4\alpha +2\gamma +3)},\quad k=0,1,\ldots ,n-1,\\ e_{n,k}&=\frac{(k+2\alpha +1)(n-k)(n+k+2\alpha +2\gamma +2)}{(2k+2\alpha +1)(2n-2k+2\gamma +1)(2n+2k+4\alpha +2\gamma +3)},\quad k=0,1,\ldots ,n-1,\\ d_{n,k}&=\frac{k(n-k+2\gamma +1)(n+k+ 2\alpha +1)}{(2k+2\alpha +1)(2n-2k+2\gamma +1)(2n+2k+4\alpha +2\gamma +3)},\quad k=1,2,\ldots ,n-1,\\ b_{n,k}&=\frac{1}{2},\quad k=0,1,\ldots ,n. \end{aligned}\nonumber \\ \end{aligned}$$
(4.1)

In Fig. 5, we can see a diagram of all possible transitions of this discrete-time QBD process.

Fig. 5
figure 5

Diagram of all possible transitions of the discrete-time QBD process corresponding with \(J_1\) for the bivariate Jacobi–Koornwinder polynomials

Full size image

At time \(t=0\) the initial state is \(Z_0=(n,k)\). The urn model will be divided in two steps. First, we consider two auxiliary urns \(\hbox {U}_1\) and \(\hbox {U}_2\). In urn \(\hbox {U}_1\), we put \(n+k+2\alpha +2\gamma +2\) blue balls and \(n+k+2\alpha +1\) red balls from the bath, and in urn \(\hbox {U}_2\), we put \(n-k+2\gamma +1\) blue balls and \(n-k\) red balls also from the bath. Then, we draw independently one ball from urn \(\hbox {U}_1\) and urn \(\hbox {U}_2\) at random with the uniform distribution. We have four possibilities:

  1. (1)

    Both balls from \(\hbox {U}_1\) and \(\hbox {U}_2\) are blue with probability

    $$\begin{aligned} \frac{n+k+2\alpha +2\gamma +2}{2n+2k+4\alpha +2\gamma +3}\times \frac{n-k+2\gamma +1}{2n-2k+2\gamma +1}. \end{aligned}$$

    Observe that this number is included in the coefficient \(a_{n,k}\) in (4.1). Then, we remove all the balls in urn \(\hbox {A}_k\) and put \(2n+4\alpha +2\gamma +3\) blue balls and \(2n+2\gamma +1\) red balls in urn \(\hbox {A}_k\). Draw one ball from \(\hbox {A}_k\) at random with the uniform distribution. If we get a blue ball, then we remove all balls in urn \(\hbox {A}_k\) and add \(n+1\) blue balls to the urn \(\hbox {A}_k\) and start over. Therefore, we have

    $$\begin{aligned} \mathbb {P}\left[ Z_{1}=(n+1,k)\, |\, Z_0=(n,k)\right] =a_{n,k}. \end{aligned}$$
  2. (2)

    Both balls from \(\hbox {U}_1\) and \(\hbox {U}_2\) are red with probability

    $$\begin{aligned} \frac{n+k+2\alpha +1}{2n+2k+4\alpha +2\gamma +3}\times \frac{n-k}{2n-2k+2\gamma +1}. \end{aligned}$$

    Observe that this number is included in the coefficient \(c_{n,k}\) in (4.1). Then, we remove all the balls in urn \(\hbox {A}_k\) and put \(2n+4\alpha +2\gamma +3\) blue balls and \(2n+2\gamma +1\) red balls in urn \(\hbox {A}_k\). Draw one ball from \(\hbox {A}_k\) at random with the uniform distribution. If we get a red ball, then we remove all balls in urn \(\hbox {A}_k\) and add \(n-1\) blue balls to the urn \(\hbox {A}_k\) and start over. Therefore, we have

    $$\begin{aligned} \mathbb {P}\left[ Z_{1}=(n-1,k)\, |\, Z_0=(n,k)\right] =c_{n,k}. \end{aligned}$$
  3. (3)

    The ball from \(\hbox {U}_1\) is blue and the ball from \(\hbox {U}_2\) is red with probability

    $$\begin{aligned} \frac{n+k+2\alpha +2\gamma +2}{2n+2k+4\alpha +2\gamma +3}\times \frac{n-k}{2n-2k+2\gamma +1}. \end{aligned}$$

    Observe that this number is included in the coefficient \(e_{n,k}\) in (4.1). Then, we remove all the balls in urn \(\hbox {A}_k\) and put \(k+2\alpha +1\) blue balls and k red balls in urn \(\hbox {A}_k\). Draw one ball from \(\hbox {A}_k\) at random with the uniform distribution. If we get a blue ball, then we remove all balls in urn \(\hbox {A}_k\) and add n blue balls to the urn \(\hbox {A}_{k+1}\) and start over. Therefore, we have

    $$\begin{aligned} \mathbb {P}\left[ Z_{1}=(n,k+1)\, |\, Z_0=(n,k)\right] =e_{n,k}. \end{aligned}$$
  4. (4)

    The ball from \(\hbox {U}_1\) is red and the ball from \(\hbox {U}_2\) is blue with probability

    $$\begin{aligned} \frac{n+k+2\alpha +1}{2n+2k+4\alpha +2\gamma +3}\times \frac{n-k+2\gamma +1}{2n-2k+2\gamma +1}. \end{aligned}$$

    Observe that this number is included in the coefficient \(d_{n,k}\) in (4.1). Then, we remove all the balls in urn \(\hbox {A}_k\) and put \(k+2\alpha +1\) blue balls and k red balls in urn \(\hbox {A}_k\). Draw one ball from \(\hbox {A}_k\) at random with the uniform distribution. If we get a red ball, then we remove all balls in urn \(\hbox {A}_k\) and add n blue balls to the urn \(\hbox {A}_{k-1}\) and start over. Therefore, we have

    $$\begin{aligned} \mathbb {P}\left[ Z_{1}=(n,k-1)\, |\, Z_0=(n,k)\right] =d_{n,k}. \end{aligned}$$

In each of the previous four possibilities there is a complementary probability. In cases (1) and (3), we may have a red ball in the second step, while in cases (2) and (4), we may have a blue ball in the second step. In all these four possibilities, we remove all balls in urn \(\hbox {A}_k\) and add n blue balls to the urn \(\hbox {A}_{k}\) and start over. The addition of these four probabilities gives 1/2. Therefore, we have

$$\begin{aligned} \mathbb {P}\left[ Z_{1}=(n,k)\, |\, Z_0=(n,k)\right] =b_{n,k}=\frac{1}{2}. \end{aligned}$$

Remark 4.1

If \(\beta \ne \alpha \) the probabilities in (2.7) will have an extra factor, so we will have to add an extra step to the previous urn model. However, the urn model is not as clear as the previous one.

Remark 4.2

It would be possible to consider an urn model taking \(\tau =1\) in (3.2) (therefore, \(\varvec{P}=J_2\)). But in this case, the coefficients \(A_{n,2}, B_{n,2}, C_{n,2}\) in (2.9)–(2.10) are way more complicated than the case we studied here. The diagram of all possible transitions will look like Figure 3 of [6]. In [6], an urn model was proposed for the orthogonal polynomials on the triangle as a consequence of finding a simple stochastic LU factorization of \(\varvec{P}\). Although it may be possible to find a LU factorization of \(\varvec{P}\) in this situation, each of the factors are not as simple as the original one, so this method is no longer convenient to find an urn model.