QBD Processes Associated with Jacobi–Koornwinder Bivariate Polynomials and Urn Models

We study a family of quasi-birth-and-death (QBD) processes associated with the so-called first family of Jacobi–Koornwinder bivariate polynomials. These polynomials are orthogonal on a bounded region typically known as the swallow tail. We will explicitly compute the coefficients of the three-term recurrence relations generated by these QBD polynomials and study the conditions under we can produce families of discrete-time QBD processes. Finally, we show an urn model associated with one special case of these QBD processes.


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 S. Karlin and J. 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 N 0 , is given by pairs of the form (n, k), where n ∈ N 0 is usually called the level, while 1 ≤ k ≤ 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 ∈ N 0 then we go back to classical birth-death chains.If r n = N for all n ∈ N 0 , where N is a positive integer, then all blocks in the Jacobi matrix have the same dimension N × N .In this case, the spectral analysis can be performed by 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 by 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 = n+d−1 n .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 Ω ⊂ 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 T. Koornwinder in [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 Figure 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 [20] and [18].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 Section 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 Section 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 Section 4, we give an urn model associated with one of the QBD processes introduced in Section 3, for the special case of β = α.

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 Then the region Ω is given by (see Figure 1) The weight function acting on this region will be given by where C is the normalizing constant To ensure integrability we need to have α, β, γ > −1, α+γ +3/2 > 0 and β + γ + 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 (u, v) satisfy the following second-order partial differential equation (after the change of variables): With this partial differential equation it is possible to generate all the monic Jacobi-Koornwinder polynomials for any values of α, β, γ.For the special cases of γ = ±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 by using the Rodrigues-type formula found in [20, Section 5].

Let us now introduce a new set of polynomials
n,k (u, v) can also be defined as where σ n,k = P α,β,γ n,k (1, 1) is given by Here we are using the standard notation for the Pochhammer symbol (a) 0 = 1, (a) n = a(a + 1) 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: (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 (1) = 1 (see [9,Section 5] for instance).
In particular we always have that a n , c n > 0, b n ≥ 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 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 where the entries of A n,1 , B n,1 and C n,1 (see (2.4)) are given by Remark 2.1.The coefficient b n,k can also be written as Also, from (2.4) it is possible to see that On the other side, the matrices A n,2 , B n,2 and C n,2 are tridiagonal matrices of the form n,1 a n,n a n,1 c ) where the entries of A n,2 , B n,2 and C n,2 (see again (2.4)) are given by , k = 0, 1, . . ., n − 2. (2.9) The coefficients for A n,i , B n,i , C n,i , i = 1, 2, can be significantly simplified for the values of γ = ±1/2.For γ = −1/2 we get and n,k = 2c n a k , while for γ = 1/2 we obtain and n,k = a n+1 (2b k − 1) Remark 2.2.The normalization of the polynomials Q α,β,γ n,k (u, v) such that Q α,β,γ n,k (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 Ω (see Figure 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 (0, 1) = (−1) n+k σ n,k , where σ n,k is given by (2.2), and the corresponding new vector polynomials Q n satisfy the three-term recurrence relations where the coefficients A n,i , B n,i , 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 α by β and β by α, except for B n,1 where we have B n,1 = B n,1 − I (changing α by β again and viceversa).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].

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.8).From the recurrence relations (2.3) we can define the following two block tridiagonal Jacobi matrices By construction of the vector polynomials Q n (see (2.3)) we always have that J i e = e, i = 1, 2, where e = (1, 1, 1, . ..) 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 We would like to see under what conditions we get a probabilistic interpretation of P .In particular we will see when P is a stochastic matrix.We immediately have that P e = e but now we need all entries of 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 P , we need to have τ a From the definition of the coefficients (2.7) and (2.9) we have the following properties , where b n is defined by (2.4).First, from (2.7), we have that a n,k , c n,k > 0 and d n,k , e n,k ≥ 0 as long as γ > −1 for α, β ≥ −1/2 and γ + 3/2 ≥ max{−α, −β} in any other case.Under these conditions on γ we have, for instance for the case n,k > 0, that the parameter τ must be chosen so that

Similar considerations can be made for the inequalities
n,k > 0. That means that the behavior of (3 − 4b n ) −1 as n ≥ 0 will be the main ingredient in order to find these upper bounds for the parameter τ .In particular, we have which behavior, for different values of α and β, will be the key to analyze the inequalities given above.Therefore let us define the following γ-dependent constant .
• In the region C= {α > −1, β > α, β < −α}, it turns out that we can choose any 0 ≤ τ ≤ min{C 1/2 , C γ+1 }, and the matrix P in (3.2) will always be stochastic.If α ≥ −1/2 then this is possible for any γ > −1, while for −1 < α < −1/2 we need to have γ + α + 3/2 > 0. This is the only case where the upper bound may depend on γ.A straightforward computation shows that Therefore, for all values of τ 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, . ..} with transition probability matrix P = (1 − τ )J 1 + τ J 2 .Thus the Karlin-McGregor representation formula (see formula (2.13) of [6]) for the (i, j) block entry of the matrix P is given by where Q n , n ≥ 0, are the vector polynomials satisfying (2.3) and W α,β,γ is the normalized weight function (2.1).The matrices Π j , j ≥ 0, are the inverses of the norms of the corresponding vector polynomials Q j , j ≥ 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 T is given by .
Since the representation of Π n is independent of the choice of the generalized inverse G n (see [6, Lemma 2.1]) we have, after some straightforward computations, that Π n is a diagonal matrix of the form Another way of writing these norms in terms of the norms of the Jacobi polynomials (2.5) is In particular we have that the family of polynomials Q n , n ≥ 0, is mutually orthogonal.Therefore, another way to write the Karlin-McGregor (3.4) formula is entry by entry According to [6, Theorem 2.5] we can construct an invariant measure π for the QBD process which is given by Here 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 After some computations it turns out that, in the range of the values of τ for which P is stochastic, this integral is divergent, and therefore (null) recurrent, if and only if −3/2 < α + γ ≤ −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.1.Instead of (3.2), we could have considered the situation where P = τ 1 J 1 + τ 2 J 2 and P e = 0, in which case we would have had a continuous-time QBD process.Since J i e = e, i = 1, 2 then we need τ 2 = −τ 1 = −τ and therefore P = τ (J 2 − J 1 ).All off-diagonal entries of 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.2.Going back to Remark 2.2 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 P = τ 1 J 1 + τ 2 J 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 continuoustime QBD process then we need P e = 0 and nonnegative off-diagonal entries.But this is possible if and only if τ 2 = 0 and τ 1 > 0, i.e. a positive scalar multiple of J 1 .If we want to have a discrete-time QBD process then we need P e = e and nonnegative (scalar) entries.This is possible if and only if τ 2 = 1 and the parameter τ 1 = τ is chosen in such a way that all entries of P are nonnegative.The entries of P = τ J 1 + J 2 are nonnegative if and only if n,k > 0. Now, from the definition (see (2.7) and (2.9)), we have Therefore, as before, the lower bounds for τ (depending also on α, β, γ) will probably depend on the behavior of the constant value n,k /(1 − b n,k ), meaning that will also have upper bounds for τ .We leave the details to the reader.

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 τ = 0 (therefore P = J 1 ) and β = α.In this section we will assume that α and γ are nonnegative integers.Consider {Z t : t = 0, 1, . ..} the discrete-time QBD process on the state space {(n, k) : 0 ≤ k ≤ n, n ∈ 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, . . . the state (n, k) will represent the number of n blue balls inside the k-th urn A k , k = 0, 1, . . ., 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 β = α the coefficients in (2.7) are simplified and given explicitly by In Figure 3 we can see a diagram of all possible transitions of this discrete-time QBD process.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 U 1 and U 2 .In urn U 1 we put n + k + 2α + 2γ + 2 blue balls and n + k + 2α + 1 red balls from the bath, and in urn U 2 we put n − k + 2γ + 1 blue balls and n − k red balls also from the bath.Then we draw independently one ball from urn U 1 and urn U 2 at random with the uniform distribution.We have four possibilities: (1) Both balls from U 1 and U 2 are blue with probability Observe that this number is included in the coefficient a n,k in (4.1).Then we remove all the balls in urn A k and put 2n + 4α + 2γ + 3 blue balls and 2n + 2γ + 1 red balls in urn A k .Draw one ball from A k at random with the uniform distribution.If we get a blue ball then we remove all balls in urn A k and add n + 1 blue balls to the urn A k and start over.Therefore we have (2) Both balls from U 1 and U 2 are red with probability Observe that this number is included in the coefficient c n,k in (4.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 A k and add n blue balls to the urn A k and start over.The addition of these four probabilities gives 1/2.Therefore we have Remark 4.1.If β = α 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 τ = 1 in (3.2) (therefore P = J 2 ).But in this case the coefficients A n,2 , B n,2 , C n,2 in (2.8)-(2.9)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 P .Although it may be possible to find a LU factorization of 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.

Figure 1 .
Figure 1.Swallow tail region Ω where the weight function is defined.

Figure 2 .
Figure 2. The regions A, B and C where τ may take different values in order to have a stochastic matrix P (courtesy of C. Juarez).

2 d 3 , 3 Figure 3 .
Figure 3. Diagram of all possible transitions of the discrete-time QBD process corresponding with J 1 for the bivariate Jacobi-Koornwinder polynomials.

( 3 )( 4 ) 1 .
1).Then we remove all the balls in urn A k and put 2n + 4α + 2γ + 3 blue balls and 2n + 2γ + 1 red balls in urn A k .Draw one ball from A k at random with the uniform distribution.If we get a red ball then we remove all balls in urn A k and add n − 1 blue balls to the urn A k and start over.Therefore we haveP [Z 1 = (n − 1, k) | Z 0 = (n, k)] = c n,k .The ball from U 1 is blue and the ball from U 2 is red with probabilityn + k + 2α + 2γ + 2 2n + 2k + 4α + 2γ + 3 × n − k 2n − 2k + 2γ + 1 .Observe that this number is included in the coefficient e n,k in (4.1).Then we remove all the balls in urn A k and put k + 2α + 1 blue balls and k red balls in urn A k .Draw one ball from A k at random with the uniform distribution.If we get a blue ball then we remove all balls in urn A k and add n blue balls to the urn A k+1 and start over.Therefore we haveP [Z 1 = (n, k + 1) | Z 0 = (n, k)] = e n,k.The ball from U 1 is red and the ball from U 2 is blue with probabilityn + k + 2α + 1 2n + 2k + 4α + 2γ + 3 × n − k + 2γ + 1 2n − 2k + 2γ +Observe that this number is included in the coefficient d n,k in (4.1).Then we remove all the balls in urn A k and put k + 2α + 1 blue balls and k red balls in urn A k .Draw one ball from A k at random with the uniform distribution.If we get a red ball then we remove all balls in urn A k and add n blue balls to the urn A k−1 and start over.Therefore we have P [Z 1 = (n, k − 1) | Z 0 = (n, k)] = d n,k .