<figure><img src="image_1.png" /></figure>

HAL
open science

Factorized $A_2$-Leonard pair

Nicolas Crampé, Meri Zaimi

► To cite this version:

Nicolas Crampé, Meri Zaimi. Factorized $A_2$-Leonard pair. The Ramanujan Journal, 2025, 66 (25), 10.1007/s11139-024-01009-y. hal-04595412

HAL Id: hal-04595412

https://hal.science/hal-04595412v1

Submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Factorized $A_2$-Leonard pair

Nicolas Crampé¹*, Meri Zaimi²†,

¹ Institut Denis-Poisson CNRS/UMR 7013 - Université de Tours - Université d'Orléans, Parc de Grandmont, 37200 Tours, France.

¹ Laboratoire d'Annecy-le-Vieux de Physique Théorique LAPTh, Université Savoie Mont Blanc, CNRS, F-74000 Annecy, France.

² Centre de Recherches Mathématiques, Université de Montréal, P.O. Box 6128, Centre-ville Station, Montréal (Québec), H3C 3J7, Canada.

Abstract

The notion of factorized $A_2$-Leonard pair is introduced. It is defined as a rank 2 Leonard pair, with actions in certain bases corresponding to the root system of the Weyl group $A_2$, and with some additional properties. The functions arising as entries of transition matrices are bivariate orthogonal polynomials (of Tratnik type) with bispectral properties. Examples of factorized $A_2$-Leonard pairs are constructed using classical Leonard pairs associated to families of orthogonal polynomials of the ($q$-)Askey scheme. The most general examples are associated to an intricate product of univariate ($q$-)Hahn and dual ($q$-)Hahn polynomials.

**Keywords:** Leonard pairs, bivariate orthogonal polynomials, bispectrality, ($q$-)Askey scheme

**MSC2020 database:** 33C80; 33C45; 16G60

<footnote>*E-mail: crampe1977@gmail.com</footnote>

<footnote>†E-mail: meri.zaimi@umontreal.ca</footnote>
# 1 Introduction

The purpose of this paper is to study a generalized notion of Leonard pair which involves bivariate orthogonal polynomials.

A Leonard pair is formed by two finite-dimensional linear transformations such that there exists a basis where one is diagonal and the other one is irreducible tridiagonal, and vice versa. This notion has been introduced in [42], and it is closely related to the families of orthogonal polynomials in the terminating branch of the (q-)Askey scheme [30]. Indeed, the matrix entries when performing a change of basis for a Leonard pair can be expressed in terms of polynomials satisfying both a three-term recurrence relation and a three-term (q-)difference equation. Leonard pairs are classified according to these families of orthogonal polynomials, the most general being the q-Racah polynomials, which are the discrete analogs of the Askey–Wilson polynomials. The recurrence and q-difference operators associated to both families of orthogonal polynomials satisfy the Askey–Wilson algebra [47]. This algebra, first introduced in [23] to study the symmetry of the 6j-symbol for the quantum group $U_q(sI_2)$, plays an important role in various related contexts such as centralizers of $U_q(sI_2)$ [12, 10], the Skein algebra [5], or symmetries of physical models [20, 29, 33] (for a review see [6]). Leonard pairs are recovered when considering finite-dimensional irreducible representations of the Askey–Wilson algebra. The classification of Leonard pairs in terms of orthogonal polynomials of the (q-)Askey scheme is an algebraic version of Leonard’s theorem [31]. This theorem is particularly useful for classifying association schemes with both P- and Q-polynomial properties in the context of algebraic combinatorics [1]. The connection between such association schemes with polynomial structures and Leonard pairs is obtained when studying the representations of the subconstituent (also called Terwilliger) algebra [39, 40, 41], which is the algebra formed by the adjacency and dual adjacency matrices.

Multivariate generalizations of the orthogonal polynomials of the (q-)Askey scheme have been considered in past works [24, 46, 36, 26, 21, 22, 28]. Algebras associated to them have also been investigated, for instance the higher rank Bannai–Ito, Racah and Askey–Wilson algebras [14, 15, 34, 16, 5, 9, 8]. More recently, multivariate generalizations of the P- and Q-polynomial properties for association schemes have been introduced and studied [3, 2, 13, 4]. However, an analog of Leonard’s theorem for multivariate polynomials is not known at the present time. In order to better understand the structure of higher rank association schemes and of multivariate orthogonal polynomials with bispectral properties, and with the purpose of eventually arriving at some classification of these objects, it is crucial to understand the higher rank generalizations of Leonard pairs.

A generalized notion of Leonard pair has been proposed in [27], based on an algebraic interpretation of the Rahman polynomials. Instead of considering a pair of linear transformations with tridiagonal actions in some eigenbases, one considers a pair of *M*-dimensional subspaces of commuting linear transformations with adjacent actions in some common eigenbases. Such pairs are related to *M*-variate polynomials which satisfy recurrence relations with a number of terms that can be higher than three. However, a systematic and complete study of all pairs satisfying the proposed definition seems difficult. In [3, 13], the bivariate polynomial structures of association schemes are characterized by the types of recurrence relations that the associated polynomials satisfy. Some of the natural examples of association schemes examined in [3] have a type such that one of the recurrence relations has at most three terms. This is the case for the non-binary Johnson schemes. It is shown in [13] that the Terwilliger algebra of these schemes satisfies some interesting properties that can be viewed as a refinement of the definition of higher rank Leonard pairs proposed in [27]. The goal of this paper is to abstract these properties and provide a framework for studying this class of refined higher rank Leonard pairs, that we call factorized Leonard pairs.

The paper will unfold as follows. Section 2 recalls the definition of the classical Leonard pairs introduced in [42], and their description in terms of a parameter array. In Section 3, the factorized Leonard pairs are defined as a refinement of the higher rank Leonard pairs of [27]. We focus in this paper on the case where the adjacency of the actions is of type $A_2$, with reference to the Weyl group root system. Some general properties of these factorized Leonard pairs are also investigated. Section
4 explains how to construct factorized Leonard pairs using classical Leonard pairs. The bivariate polynomials associated to a factorized Leonard pair are discussed in Section 5. Examples of factorized $A_2$-Leonard pairs associated to families of orthogonal polynomials of the Askey scheme (type II) and its q-analog (type I) are provided respectively in Sections 6 and 7. One example (associated to the dual (q-)Hahn and (q-)Hahn polynomials) in each section is described in detail, the others can be obtained through specializations and limits. In Section 8, a higher rank generalized notion of factorized $A_M$-Leonard pair is briefly introduced. Section 9 discusses some outlooks. Finally, in Appendix A, a list of parameter arrays for the classical Leonard pairs is recalled, and useful contiguity relations for the associated polynomials are provided.

## 2 Leonard pairs

We recall in this section the definition of Leonard pairs and their classification [42] (see also [44]).

**Definition 2.1.** Let $V$ denote a vector space over a field $K$ with finite positive dimension. By a Leonard pair on $V$ we mean an ordered pair $(A, A^*)$, where $A: V \to V$ and $A^*: V \to V$ are linear transformations which satisfy the following two properties:

(i) There exists a basis for $V$ with respect to which the matrix representing $A$ is irreducible tridiagonal and the matrix representing $A^*$ is diagonal;

(ii) there exists a basis for $V$ with respect to which the matrix representing $A^*$ is irreducible tridiagonal and the matrix representing $A$ is diagonal.

In this paper, we restrict ourselves to the case where the field $K$ is of characteristic zero.

The classification of the Leonard pairs is well-known and is closely related to the orthogonal polynomials of the (q-)Askey scheme [42, 43] and to the Askey-Wilson algebra [45]. We follow here the notations in [43].

To each Leonard pair $(A, A^*)$ with $\dim(V) = d+1$ ($d$ is called the diameter of the Leonard pair), a set of parameters is associated

$$ \Phi = (\theta_i, \theta_i^*, i = 0, \dots, d; \varphi_j, \phi_j, j = 1, \dots, d), \qquad (2.1) $$

called parameter array. To be precise, there is a one-to-one correspondence between the parameter arrays and the Leonard systems, which are Leonard pairs with orderings of the primitive idempotents. Each Leonard pair is associated to four Leonard systems. We shall always assume that an ordering of the primitive idempotents is fixed for a given Leonard pair such that it forms a Leonard system, but for simplicity we will not mention explicitly the Leonard system.

From the parameter array (2.1), one gets:

* $\{\theta_i | i = 0, \dots, d\}$ are the eigenvalues of $A$;

* $\{\theta_i^* | i = 0, \dots, d\}$ are the eigenvalues of $A^*$;

* in the basis $\{v_0^*, \dots, v_d^*\}$ where $A$ is tridiagonal and $A^*$ is diagonal, they are written as

$$ Av_i^* = b_{i-1}v_{i-1}^* + a_i v_i^* + c_{i+1}v_{i+1}^*, \qquad (2.2) $$

$$ A^*v_i^* = \theta_i^*v_i^*, \qquad (2.3) $$

where, by convention $b_{-1} = c_{d+1} = 0$ and

$$ a_i = a_i(\Phi) = \theta_i + \frac{\varphi_i}{\theta_i^* - \theta_{i-1}^*} + \frac{\varphi_{i+1}}{\theta_i^* - \theta_{i+1}^*}, \qquad (2.4a) $$

$$ b_i = b_i(\Phi) = \varphi_{i+1} \frac{\prod_{h=0}^{i-1} (\theta_i^* - \theta_h^*)}{\prod_{h=0}^{i} (\theta_{i+1}^* - \theta_h^*)}, \qquad c_i = c_i(\Phi) = \phi_i \frac{\prod_{h=0}^{d-i-1} (\theta_i^* - \theta_{d-h}^*)}{\prod_{h=0}^{d-i} (\theta_{i-1}^* - \theta_{d-h}^*)}; \qquad (2.4b) $$
• in the basis {$v_0, \dots, v_d$} where A is diagonal and $A^*$ is tridiagonal, they are written as

$$A v_x = \theta_x v_x, \qquad (2.5)$$

$$A^* v_x = b_{x-1}^* v_{x-1} + a_x^* v_x + c_{x+1}^* v_{x+1}, \qquad (2.6)$$

where, by convention $b_{-1}^* = c_{d+1}^* = 0$,

$$a_i^* = a_i(\Phi^*), \quad b_i^* = b_i(\Phi^*), \quad c_i^* = c_i(\Phi^*), \qquad (2.7)$$

and the transformed parameter array is

$$\Phi^* = (\theta_i^*, \theta_i, i = 0, \dots, d; \varphi_j, \phi_{d-j+1}, j = 1, \dots, d). \qquad (2.8)$$

The classification of Leonard pairs is obtained by providing the complete list of all the possible parameter arrays. We recall the parameter arrays used in this paper in Appendix A.

**Remark 2.2.** Instead of relations (2.5)-(2.6), one can have

$$A v_x = \theta_x v_x, \qquad (2.9)$$

$$A^* v_x = \frac{\alpha_{x-1}}{\alpha_x} b_{x-1}^* v_{x-1} + a_x^* v_x + \frac{\alpha_{x+1}}{\alpha_x} c_{x+1}^* v_{x+1}, \qquad (2.10)$$

for $\alpha_x$ non-vanishing free parameters. They correspond to the normalizations of the vectors $v_x$. Similarly one can modify the normalizations of $v_i^*$ such that relation (2.2) is modified. ▲

For any Leonard pair, the change of basis between both bases is known explicitly and is given in terms of orthogonal polynomials of the (q-)Askey scheme. Explicitly, to each parameter array $\Phi$, we associate the following polynomial of degree $i$ w.r.t. $\lambda$ [43]:

$$u_i(\lambda) = \sum_{k=0}^{i} \prod_{j=0}^{k-1} \frac{(\lambda - \theta_j)(\theta_i^* - \theta_j^*)}{\varphi_{j+1}}. \qquad (2.11)$$

Let us define

$$p_{xi}(\Phi) = k_i(\Phi) u_i(\theta_x), \qquad (2.12)$$

where

$$k_i(\Phi) = \prod_{h=1}^{i} \frac{\varphi_h}{\phi_h} \frac{\prod_{\substack{h=1 \\ h \neq i}}^{d} (\theta_0^* - \theta_h^*)}{\prod_{\substack{h=0 \\ h \neq i}}^{d} (\theta_i^* - \theta_h^*)}. \qquad (2.13)$$

Let us remark that $p_{xi}(\Phi)$ is denoted by $v_i(\theta_x)$ in [43]. These coefficients $p_{xi}(\Phi)$ allow us to connect the two bases:

$$v_i^* = \sum_{x=0}^{d} p_{xi}(\Phi) v_x. \qquad (2.14)$$

These coefficients are given explicitly in Appendix A for different parameter arrays, and satisfy the recurrence and the difference relations:

$$\theta_x p_{xi}(\Phi) = b_{i-1} p_{x,i-1}(\Phi) + a_i p_{xi}(\Phi) + c_{i+1} p_{x,i+1}(\Phi), \qquad (2.15a)$$

$$\theta_i^* p_{xi}(\Phi) = b_x^* p_{x+1,i}(\Phi) + a_x^* p_{xi}(\Phi) + c_x^* p_{x-1,i}(\Phi). \qquad (2.15b)$$
They also satisfy the following orthogonality relation

$$ \sum_{i=0}^{d} p_{xi}(\Phi) p_{iy}(\Phi^{*}) = \delta_{xy} \nu(\Phi), \qquad (2.16) $$

where $\Phi^*$ is given by (2.8) and

$$ \nu(\Phi) = \nu(\Phi^*) = \prod_{h=1}^{d} \frac{(\theta_0 - \theta_h)(\theta_0^* - \theta_h^*)}{\phi_h}. \qquad (2.17) $$

# 3 Definitions and properties

In this section, we provide the definition of a factorized Leonard pair, generalizing the notion of Leonard pairs as well as some of its properties.

## 3.1 Definition of a factorized $A_2$-Leonard pair

Let $N$ be a positive integer and $\mathcal{D}$ be the set constituted of couples of non-negative integers whose sum is at most $N$:

$$ \mathcal{D} = \{(x, y) \mid x, y \in \mathbb{Z}_{\ge 0}, x + y \le N\}. \qquad (3.1) $$

Two couples $(x, y), (x', y') \in \mathcal{D}$ are called $A_2$-adjacent, $(x, y)$ adj $(x', y')$, if

$$ (x - x', y - y') \in A_2, \qquad (3.2) $$

where $A_2$ is the following subset of $\mathbb{Z}^2$:

$$ A_2 = \{(0, 0), (1, -1), (1, 0), (0, 1), (-1, 1), (-1, 0), (0, -1)\}. \qquad (3.3) $$

The notation $A_2$ for the previous set comes from the fact that it corresponds to the root system of the Weyl group $A_2$.

**Definition 3.1.** [27] Let $V$ be a $\mathbb{C}$-vector space. The pair $(H, H^*)$ is a $A_2$-Leonard pair on the domain $\mathcal{D}$ if the following statements are satisfied:

(i) $H$ is a two-dimensional subspace of $End(V)$ whose elements are diagonalizable and mutually commute;

(ii) $H^*$ is a two-dimensional subspace of $End(V)$ whose elements are diagonalizable and mutually commute;

(iii) there exists a bijection $\lambda \mapsto V_\lambda$ from $\mathcal{D}$ to the common eigenspaces of $H$ such that, for all $\lambda \in \mathcal{D}$,

$$ H^*V_\lambda \subseteq \sum_{\substack{\chi \in \mathcal{D} \\ \chi \text{ adj } \lambda}} V_\chi; \qquad (3.4) $$

(iv) there exists a bijection $\lambda \mapsto V_\lambda^*$ from $\mathcal{D}$ to the common eigenspaces $V_\lambda^*$ of $H^*$ such that, for all $\lambda \in \mathcal{D}$,

$$ HV_{\lambda}^{*} \subseteq \sum_{\substack{\chi \in \mathcal{D} \\ \chi \text{ adj } \lambda}} V_{\chi}^{*}; \qquad (3.5) $$

(v) there does not exist a subspace $W$ of $V$ such that $HW \subseteq W$, $H^*W \subseteq W$, $W \neq 0$, $W \neq V$;

(vi) $dim(V_\lambda) = dim(V_\lambda^*) = 1$, for all $\lambda \in \mathcal{D}$;
(vii) there exists a non-degenerate symmetric bilinear form $\langle, \rangle$ on $V$ such that, for $\lambda, \chi \in \mathcal{D}$ and $\lambda \neq \chi$,

$$ \langle V_{\lambda}, V_{\chi} \rangle = 0, \quad \langle V_{\lambda}^{*}, V_{\chi}^{*} \rangle = 0. \qquad (3.6) $$

Let us notice that the dimension of the vector space $V$ is finite and is given more precisely by $\dim(V) = |\mathcal{D}| = \binom{N+2}{2}$.

**Remark 3.2.** The previous definition has been given in the conclusion of [27]. The definition of adjacency used in this paper is slightly modified in comparison to [27], which explains the difference in relations (3.4) and (3.5). ▲

The classification of the usual Leonard pairs has been recalled in Section 2 and is closely related to the orthogonal polynomials of the (q-)Askey scheme. The classification of the $A_2$-Leonard pairs would provide similar results but for bivariate polynomials. For a first step in this direction, we propose to add further conditions to the previous definition such that a classification seems possible. We arrive at the following definition of factorized Leonard pairs.

**Definition 3.3.** Let $(H, H^*)$ be a $A_2$-Leonard pair on the domain $\mathcal{D}$ and $X, Y \in \text{End}(V)$ (resp. $X^*, Y^*$) be a basis of $H$ (resp. $H^*$). The pair $(H, H^*)$ is a factorized $A_2$-Leonard pair if the following statements are satisfied:

(viii) the eigenvalues of $X$ and $Y^*$ satisfy, for any $(x, y), (i, j) \in \mathcal{D}$,

$$ (X - t_x \mathbb{1}) V_{(x,y)} = 0, \quad (Y^* - \theta_j^* \mathbb{1}) V_{(i,j)}^* = 0, \qquad (3.7) $$

where $\mathbb{1}$ is the identity in $\text{End}(V)$ and $t_x \neq t_{x'} \Leftrightarrow x \neq x'$ and $\theta_j^* \neq \theta_{j'}^* \Leftrightarrow j \neq j'$;

(ix) in addition to the commutation relations $[X, Y] = 0$ and $[X^*, Y^*] = 0$, one has

$$ [X, Y^*] = 0; \qquad (3.8) $$

(x) each common eigenspace of $X$ and $Y^*$ has dimension one;

(xi) the following relations hold, for $(i, j), (i + 1, j) \in \mathcal{D}$ and $(x, y), (x, y + 1) \in \mathcal{D}$,

$$ \langle V_{(i+1,j)}^*, X V_{(i,j)}^* \rangle \neq 0, \quad \langle V_{(i,j)}, X V_{(i+1,j)}^* \rangle \neq 0, \qquad (3.9a) $$

$$ \langle V_{(x,y+1)}, Y^* V_{(x,y)} \rangle \neq 0, \quad \langle V_{(x,y)}, Y^* V_{(x,y+1)} \rangle \neq 0. \qquad (3.9b) $$

The ordered 4-tuple $(X, Y; X^*, Y^*)$ is called a basis of the factorized $A_2$-Leonard pair.

**Remark 3.4.** In the previous definition, one could replace the pair $(X, Y^*)$ by $(X^*, Y)$. ▲

In the following, we shall use the shortened name factorized Leonard pair instead of factorized $A_2$-Leonard pair.

## 3.2 Properties of a factorized Leonard pair

In this section, let $(H, H^*)$ be a factorized Leonard pair with basis $(X, Y; X^*, Y^*)$. We provide different properties for this pair, allowing us to study it in detail.

The first property concerns a symmetry of its basis.

**Proposition 3.5.** For $a_1, a_2, a_1^*, a_2^*$ invertible elements of the field $K$ and $b_1, b_2, b_1^*, b_2^*, c, c^* \in K$, the 4-tuple

$$ (a_1 X + b_1 \mathbb{1}, a_2 Y + c X + b_2 \mathbb{1} ; a_1^* X^* + c^* Y^* + b_1^* \mathbb{1}, a_2^* Y^* + b_2^* \mathbb{1}), \qquad (3.10) $$

where $\mathbb{1}$ is the identity in $\text{End}(V)$, is another basis of the factorized Leonard pair $(H, H^*)$.
Proof. It is straightforward to check that all the properties of the definition of a factorized Leonard pair are satisfied by (3.10). □

To study further the factorized Leonard pairs, it is useful to introduce some projectors. Indeed, let us associate to the factorized Leonard pair the following projectors:

$$E_x = \prod_{\substack{x'=0 \\ x' \neq x}}^{N} \frac{X - t_{x'}}{t_x - t_{x'}}, \quad E_j^* = \prod_{\substack{j'=0 \\ j' \neq j}}^{N} \frac{Y^* - \theta_{j'}^*}{\theta_j^* - \theta_{j'}^*}, \qquad (3.11)$$

which are such that $XE_x = t_x E_x$ and $Y^*E_j^* = \theta_j^* E_j^*$. It follows from (3.8) and (3.11) that

$$E_x E_{x'} = \delta_{xx'} E_x, \quad E_j^* E_{j'}^* = \delta_{jj'} E_j, \quad \sum_{x=0}^{N} E_x = \mathbb{1}, \quad \sum_{j=0}^{N} E_j^* = \mathbb{1}, \quad [E_x, E_j^*] = 0. \qquad (3.12)$$

**Lemma 3.6.** For $0 \le x, j \le N$, the dimension of the vector space $E_x E_j^* V$ is given by

$$\dim(E_x E_j^* V) = \begin{cases} 1, & \text{if } x + j \le N, \\ 0, & \text{otherwise.} \end{cases} \qquad (3.13)$$

Proof. $E_x E_j^* V$ is a common vector space of $X$ and $Y^*$ or is zero. Then due to the point (x) of Definition 3.3, one deduces that, for $0 \le x, j \le N$,

$$\dim(E_x E_j^* V) \le 1. \qquad (3.14)$$

Let us remark that $\dim(E_0 V) = N + 1$ (due to the point (iii) of Definition 3.1). It follows from this remark and equation (3.14) that for any $j = 0, 1, \dots, N$,

$$\dim(E_0 E_j^* V) = 1. \qquad (3.15)$$

Since $\dim(E_1 V) = N$ and $\dim(E_N^* V) = 1$ (due to the points (iii), (iv) of Definition 3.1), one deduces using equations (3.14) and (3.15) that for any $j = 0, 1, \dots, N - 1$,

$$\dim(E_1 E_j^* V) = 1, \quad \text{and} \quad \dim(E_1 E_N^* V) = 0. \qquad (3.16)$$

Recursively, one demonstrates the result of the lemma. □

From the previous lemma, one deduces that each common eigenspace $\tilde{V}_{(x,j)}$ ($0 \le x, j \le N$ and $x + j \le N$) of $X$ and $Y^*$ can be written as follows:

$$\tilde{V}_{(x,j)} = E_x E_j^* V = E_j^* E_x V, \qquad (3.17)$$

such that

$$(X - t_x \mathbb{1}) \tilde{V}_{(x,j)} = 0, \quad (Y^* - \theta_j^* \mathbb{1}) \tilde{V}_{(x,j)} = 0. \qquad (3.18)$$

**Lemma 3.7.** For any $j = 0, 1, \dots, N$, $(X, X^*)$ restricted on $E_j^* V$ is a Leonard pair with diameter $N - j$ and, for any $x = 0, 1, \dots, N$, $(Y, Y^*)$ restricted on $E_x V$ is a Leonard pair with diameter $N - x$.

Proof. Let us fix $j \in \{0, 1, \dots, N\}$. On one hand, the subspace $E_j^* V$ can be written as

$$E_j^* V = \sum_{\lambda \in \mathcal{D}} E_j^* V_\lambda = \sum_{i=0}^{N-j} V_{(i,j)}^*. \qquad (3.19)$$
The subspace $V_{(i,j)}^*$ is invariant under $X^*$ for $(i,j) \in \mathcal{D}$. Using condition (iv) of Definition 3.1 and the commutation relation $[X, Y^*] = 0$, one finds for any $i = 0, 1, \dots, N-j$

$$XV_{(i,j)}^* = E_j^*XV_{(i,j)}^* \subseteq E_j^* \sum_{\lambda \text{ adj } (i,j)} V_{\lambda}^* = V_{(i,j)}^* + V_{(i+1,j)}^* + V_{(i-1,j)}^*. \quad (3.20)$$

Recall that the spaces $V_{(i,j)}^*$ have dimension 1 from condition (vi) of Definition 3.1 and denote $V_{i,j}^*$ one non-vanishing vector of this vector space. The set $\{V_{i,j}^*\}_{i=0}^{N-j}$ is a basis of $E_j^*V$. In this basis, $X$ is represented by a tridiagonal matrix and $X^*$ by a diagonal matrix. In addition, the matrix representing $X$ is an irreducible tridiagonal matrix due to relations (3.9a) i.e. there are no vanishing entries in the lower and upper diagonals. One deduces that there does not exist a subspace $W_j$ of $E_j^*V$ such that $XW_j \subseteq W_j$, $X^*W_j \subseteq W_j$, $W_j \neq 0$, $W_j \neq E_j^*V$ (the fact that all the eigenvalues of $X^*$ restricted to $E_j^*V$ are different two by two has also been used).

On the other hand, one can write

$$E_j^*V = \sum_{\lambda \in \mathcal{D}} E_j^* \tilde{V}_{\lambda} = \sum_{x=0}^{N-j} \tilde{V}_{(x,j)}. \quad (3.21)$$

The subspace $\tilde{V}_{(x,j)}$ is invariant under $X$ and one can use condition (iii) of Definition 3.1 together with $[X^*, Y^*] = 0$ to get for any $x = 0, 1, \dots, N-j$

$$X^* \tilde{V}_{(x,j)} = E_j^* X^* E_x V \subseteq E_j^* (E_x V + E_{x+1} V + E_{x-1} V) = \tilde{V}_{(x,j)} + \tilde{V}_{(x+1,j)} + \tilde{V}_{(x-1,j)}. \quad (3.22)$$

Let us denote $\tilde{V}_{x,j}$ a basis vector of the one-dimensional vector space $\tilde{V}_{(x,j)}$. The set $\{\tilde{V}_{x,j}\}_{x=0}^{N-j}$ is a basis of $E_j^*V$. In this basis, $X$ is represented by a diagonal matrix and $X^*$ by an irreducible tridiagonal matrix. It follows that $X$ and $X^*$ restricted on $E_j^*V$ form a Leonard pair.

The proof that $(Y, Y^*)$ restricted on $E_x V$ is a Leonard pair is similar. In this case, $Y$ is diagonal and $Y^*$ is irreducible tridiagonal on $\{V_{(x,y)}\}_{y=0}^{N-x}$, and $Y$ is irreducible tridiagonal and $Y^*$ is diagonal on $\{\tilde{V}_{(x,j)}\}_{j=0}^{N-x}$. $\square$

**Lemma 3.8.** For a factorized Leonard pair, statement (v) of Definition 3.1 concerning the irreducibility is redundant.

*Proof.* The subspace $\tilde{V}_{(x,j)}$ is invariant under $X$ and $Y^*$, and the couples of eigenvalues $(t_x, \theta_j^*)$ are two by two different. Moreover, from the proof of the previous lemma, one gets, for $(x, j), (x+1, j) \in \mathcal{D}$,

$$\langle \tilde{V}_{(x+1,j)}, X^* \tilde{V}_{(x,j)} \rangle \neq 0, \qquad \langle \tilde{V}_{(x,j)}, X^* \tilde{V}_{(x+1,j)} \rangle \neq 0, \qquad (3.23a)$$

and, for $(x, j), (x, j+1) \in \mathcal{D}$,

$$\langle \tilde{V}_{(x,j+1)}, Y \tilde{V}_{(x,j)} \rangle \neq 0, \qquad \langle \tilde{V}_{(x,j)}, Y \tilde{V}_{(x,j+1)} \rangle \neq 0. \qquad (3.23b)$$

Therefore one deduces that $X, Y, X^*$ and $Y^*$ act irreducibly on $V$. $\square$

# 4 Construction of factorized Leonard pairs

In this section, we describe a general method for constructing factorized Leonard pairs.
## 4.1 Parameter arrays associated to a factorized Leonard pair

Let $(X, Y; X^*, Y^*)$ be a basis of a factorized Leonard pair. From Lemma 3.7, $(X, X^*)$ restricted on $E_j^*V$ is a Leonard pair of diameter $N-j$, for any $j=0, \dots, N$. From the result recalled in Section 2, we can associate to each of these Leonard pairs a parameter array given explicitly by

$$ \Phi_j = (t_x, t_x^*(j), x = 0, \dots, N-j; f_x(j), g_x(j), x = 1, \dots, N-j). \quad (4.1a) $$

Similarly $(Y, Y^*)$ restricted on $E_x V$ is a Leonard pair of diameter $N-x$, for any $x=0, \dots, N$ and we can associate to each of these pairs the following parameter array

$$ \Psi_x = (\theta_j(x), \theta_j^*, j = 0, \dots, N-x; \varphi_j(x), \phi_j(x), j = 1, \dots, N-x). \quad (4.1b) $$

Let us emphasize that, for a factorized Leonard pair, the parameters in $\Phi_j$ can depend explicitly on $j$ except the eigenvalues $t_x$ because of statement (viii) of Definition 3.3. Similarly, the dual eigenvalues in $\Psi_x$ do not depend on $x$. We can visualize the different parameter arrays associated to the different Leonard pairs as in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Above each vertex, a couple of commuting operators is written and, below, the vectors of the basis where these operators are diagonal. At each oriented edge is associated the parameter array of the Leonard pair given by the operators above this edge.</figcaption></figure>

Due to Proposition 3.5, the following transformations on the parameter arrays (4.1) associated to a factorized Leonard pair,

$$ \Phi_j = (a_1 t_x + b_1, a_2 t_x^*(j) + c t_x + b_2, x = 0, \dots, N-j; a_1 a_2 f_x(j), a_1 a_2 g_x(j), x = 1, \dots, N-j), \quad (4.2a) $$

$$ \Psi_x = (a_1^* \theta_j(x) + c^* \theta_j^* + b_1^*, a_2^* \theta_j^* + b_2^*, j = 0, \dots, N-x; a_1^* a_2^* \varphi_j(x), a_1^* a_2^* \phi_j(x), j = 1, \dots, N-x), \quad (4.2b) $$

with $a_1, a_2, a_1^*, a_2^*, b_1, b_2, b_1^*, b_2^*, c, c^* \in K$ and $a_1, a_2, a_1^*, a_2^*$ invertible, are also the associated parameter arrays of a factorized Leonard pair.

From the knowledge of the parameter arrays, one can obtain the actions of the elements of the factorized Leonard pair basis $(X, Y; X^*, Y^*)$ on the different bases of the vector space $V$. As before, we denote by $V_{x,y}$ (resp. $\tilde{V}_{x,j}, V_{i,j}^*$) a basis vector of the one-dimensional space $V_{(x,y)}$ (resp. $\tilde{V}_{(x,j)}, V_{(i,j)}^*$). These different vectors are also displayed in Figure 1. In the eigenspace $E_j^*V$, one gets $(Y^* - \theta_j^* \mathbb{I}) E_j^*V = 0$ and thus

$$ Y^* \tilde{V}_{x,j} = \theta_j^* \tilde{V}_{x,j}, \quad Y^* V_{i,j}^* = \theta_j^* V_{i,j}^*. \quad (4.3) $$

By definition of the parameter array $\Phi_j$ one gets:

$$ X \tilde{V}_{x,j} = t_x \tilde{V}_{x,j}, \quad X V_{i,j}^* = b_{i-1}(\Phi_j) V_{i-1,j}^* + a_i(\Phi_j) V_{i,j}^* + c_{i+1}(\Phi_j) V_{i+1,j}^*, \quad (4.4) $$

$$ X^* \tilde{V}_{x,j} = c_{x+1}(\Phi_j^*) \tilde{V}_{x+1,j} + a_x(\Phi_j^*) \tilde{V}_{x,j} + b_{x-1}(\Phi_j^*) \tilde{V}_{x-1,j}, \quad X^* V_{i,j}^* = t_i^*(j) V_{i,j}^*. \quad (4.5) $$

In the eigenspace $E_x V$, one gets $X E_x V = t_x E_x V$ and thus

$$ X V_{x,y} = t_x V_{x,y}, \quad X \tilde{V}_{x,j} = t_x \tilde{V}_{x,j}, \quad (4.6) $$

and by definition of the parameter array $\Psi_x$, one gets:

$$ Y V_{x,y} = \theta_y(x) V_{x,y}, \quad Y \tilde{V}_{x,j} = c_{j+1}(\Psi_x) \tilde{V}_{x,j+1} + a_j(\Psi_x) \tilde{V}_{x,j} + b_{j-1}(\Psi_x) \tilde{V}_{x,j-1}, \quad (4.7) $$

$$ Y^* V_{x,y} = b_{y-1}(\Psi_x^*) V_{x,y-1} + a_y(\Psi_x^*) V_{x,y} + c_{y+1}(\Psi_x^*) V_{x,y+1}, \quad Y^* \tilde{V}_{x,j} = \theta_j^* \tilde{V}_{x,j}. \quad (4.8) $$
**Remark 4.1.** As explained in Remark 2.2, an element of a Leonard pair can always be conjugated by a diagonal matrix. Therefore the actions of $X^*$ and $Y$ on $\tilde{V}_{x,j}$ can be transformed as follows:

$$X^* \tilde{V}_{x,j} = \frac{\alpha_{x+1}^*(j)}{\alpha_x^*(j)} c_{x+1}(\Phi_j^*) \tilde{V}_{x+1,j} + a_x(\Phi_j^*) \tilde{V}_{x,j} + \frac{\alpha_{x-1}^*(j)}{\alpha_x^*(j)} b_{x-1}(\Phi_j^*) \tilde{V}_{x-1,j}, \quad (4.9)$$

$$Y \tilde{V}_{x,j} = \frac{\alpha_{j+1}(x)}{\alpha_j(x)} c_{j+1}(\Psi_x) \tilde{V}_{x,j+1} + a_j(\Psi_x) \tilde{V}_{x,j} + \frac{\alpha_{j-1}(x)}{\alpha_j(x)} b_{j-1}(\Psi_x) \tilde{V}_{x,j-1}. \quad (4.10)$$

These new coefficients could help us to construct a factorized Leonard pair. However, in all the examples found these coefficients are not useful. To simplify the presentation of this article, we choose to present the results without these coefficients. ▲

The entries of the transition matrix associated to the Leonard pair $(X, X^*)$ restricted on $E_j^*V$ are given by $p_{xi}(\Phi_j)$, defined in (2.12). Explicitly, one has

$$V_{i,j}^* = \sum_{x=0}^{N-j} p_{xi}(\Phi_j) \tilde{V}_{x,j}, \quad \tilde{V}_{x,j} = \frac{1}{\nu(\Phi_j)} \sum_{i=0}^{N-j} p_{ix}(\Phi_j^*) V_{i,j}^*. \quad (4.11)$$

In the second relation above, we have used the inverse of the transition matrix and relation (2.16). The entries $p_{xi}(\Phi_j)$ are indicated below the edge associated to the pair $(X, X^*)$ in Figure 1. Similarly, $p_{yj}(\Psi_x)$ is the polynomial associated to the parameter array 4.1b such that

$$\tilde{V}_{x,j} = \sum_{y=0}^{N-x} p_{yj}(\Psi_x) V_{x,y}, \quad V_{x,y} = \frac{1}{\nu(\Psi_x)} \sum_{j=0}^{N-x} p_{jy}(\Psi_x^*) \tilde{V}_{x,j}. \quad (4.12)$$

The previous relations allow us to prove the following proposition.

**Proposition 4.2.** The change of basis from $V_{x,y}$ to $V_{i,j}^*$ is given by

$$V_{i,j}^* = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} p_{xi}(\Phi_j) p_{yj}(\Psi_x) V_{x,y}. \quad (4.13)$$

The fact that this change of basis associated to a factorized Leonard pair is written as a product of two polynomials justifies the name factorized.

It remains to compute the actions of $Y$ on $V_{i,j}^*$, and of $X^*$ on $V_{x,y}$.

## 4.2 Action of Y on $V_{i,j}^*$

Using the change of basis (4.11) and the action (4.7), one gets a first expression for the action of Y on $V_{i,j}^*$

$$Y V_{i,j}^* = \sum_{x=0}^{N-j} p_{xi}(\Phi_j) Y \tilde{V}_{x,j} \quad (4.14)$$

$$= \sum_{x=0}^{N-j} p_{xi}(\Phi_j) \left( c_{j+1}(\Psi_x) \tilde{V}_{x,j+1} + a_j(\Psi_x) \tilde{V}_{x,j} + b_{j-1}(\Psi_x) \tilde{V}_{x,j-1} \right). \quad (4.15)$$

In order to obtain a factorized Leonard pair, this action has to be written as (see statement (iv) of the definition):

$$Y V_{i,j}^* = \sum_{(\epsilon, \epsilon') \in A_2} C_{i,j}^{\epsilon, \epsilon'} V_{i+\epsilon, j+\epsilon'}^* = \sum_{(\epsilon, \epsilon') \in A_2} C_{i,j}^{\epsilon, \epsilon'} \sum_{x=0}^{N-j-\epsilon'} p_{x, i+\epsilon}(\Phi_{j+\epsilon'}) \tilde{V}_{x, j+\epsilon'}, \quad (4.16)$$
where all the coefficients $C_{i,j}$ must be independent of $x$ and $y$. Comparing (4.15) and (4.16), one finds

$$c_{j+1}(\Psi_x)p_{xi}(\Phi_j) = \sum_{\epsilon=-1,0} C_{i,j}^{\epsilon,+} p_{x,i+\epsilon}(\Phi_{j+1}), \quad (4.17a)$$

$$a_j(\Psi_x)p_{xi}(\Phi_j) = \sum_{\epsilon=0,\pm 1} C_{i,j}^{\epsilon,0} p_{x,i+\epsilon}(\Phi_j), \quad (4.17b)$$

$$b_{j-1}(\Psi_x)p_{xi}(\Phi_j) = \sum_{\epsilon=0,\pm 1} C_{i,j}^{\epsilon,-} p_{x,i+\epsilon}(\Phi_{j-1}). \quad (4.17c)$$

To get a factorized Leonard pair, the previous constraints must be solved. These provide strong restrictions on the parameter arrays which are admissible.

## 4.3 Action of $X^*$ on $V_{x,y}$

Similarly as in the previous subsection, the action of $X^*$ on $V_{x,y}$ provides constraints between the parameter arrays. Here one must use the change of basis (4.12).

Demanding that the action of $X^*$ on $V_{x,y}$ takes the following form (see statement (iii) of the definition)

$$X^*V_{x,y} = \sum_{(\epsilon, \epsilon') \in A_2} D_{x,y}^{\epsilon, \epsilon'} V_{x+\epsilon, y+\epsilon'}, \quad (4.18)$$

where all the coefficients $D_{x,y}$ must be independent of $i$ and $j$, leads to the following constraints

$$c_{x+1}(\Phi_j^*) p_{jy}(\Psi_x^*) = \frac{\nu(\Psi_x)}{\nu(\Psi_{x+1})} \sum_{\epsilon'=-1,0} D_{x,y}^{+,\epsilon'} p_{j,y+\epsilon'}(\Psi_{x+1}^*), \quad (4.19a)$$

$$a_x(\Phi_j^*) p_{jy}(\Psi_x^*) = \sum_{\epsilon'=0,\pm 1} D_{x,y}^{0,\epsilon'} p_{j,y+\epsilon'}(\Psi_x^*), \quad (4.19b)$$

$$b_{x-1}(\Phi_j^*) p_{jy}(\Psi_x^*) = \frac{\nu(\Psi_x)}{\nu(\Psi_{x-1})} \sum_{\epsilon'=0,\pm 1} D_{x,y}^{-,\epsilon'} p_{j,y+\epsilon'}(\Psi_{x-1}^*). \quad (4.19c)$$

## 4.4 Resolution of the constraints (4.17) and (4.19)

To solve the constraints (4.17) and (4.19), we pick all the possible choices for the parameter arrays $\Phi_j$ and $\Psi_x$ ($0 \le x, j \le N$), recalled in Appendix A. Then, for each choice, we compute $c_{j+1}(\Psi_x)$, $a_j(\Psi_x)$, $b_{j-1}(\Psi_x)$. We want $a_j(\Psi_x)$ to be an affine transformation of the function depending on $x$ in the L.H.S. of the recurrence relation (2.15a) for $p_{xi}(\Phi_j)$. Namely, we want

$$a_j(\Psi_x) = z_j^0 t_x + u_j, \quad (4.20)$$

where $z_j^0$ and $u_j$ are functions only of $j$ and we recall that $t_x$ are the eigenvalues associated to $\Phi_j$. If $z_j^0$ and $u_j$ exist, we can compute $C_{i,j}^{\epsilon,0}$ and solve the constraints (4.17b):

$$C_{i,j}^{+,0} = z_j^0 b_{i-1}(\Phi_j), \quad (4.21)$$

$$C_{i,j}^{-,0} = z_j^0 a_i(\Phi_j) + u_j, \quad (4.22)$$

$$C_{i,j}^{+,0} = z_j^0 c_{i+1}(\Phi_j). \quad (4.23)$$

For both other constraints (4.17a) and (4.17c), we demand that $c_{j+1}(\Psi_x)$ and $b_{j-1}(\Psi_x)$ be the functions depending on $x$ in the L.H.S. of contiguity equations for $p_{xi}(\Phi_j)$. Such equations have the following forms:

$$\lambda_x^+(\Phi_j) p_{x,i}(\Phi_j) = \gamma_i^{0+}(\Phi_j) p_{x,i}(\Phi_{j+1}) + \gamma_i^{-+}(\Phi_j) p_{x,i-1}(\Phi_{j+1}), \quad (4.24a)$$

$$\lambda_x^-(\Phi_j) p_{x,i}(\Phi_j) = \gamma_i^{0-}(\Phi_j) p_{x,i+1}(\Phi_{j-1}) + \gamma_i^{+-}(\Phi_j) p_{x,i}(\Phi_{j-1}), \quad (4.24b)$$
where $\lambda$ and $\gamma$ are functions to be determined. Let us emphasize that in the R.H.S. of the previous equations there is a shift $j \pm 1$ for the parameter array in the argument of the polynomials. In particular, this shift means that the polynomials on the R.H.S. and L.H.S. of these equations are associated to parameter arrays with different diameters. This type of equations is usually called contiguity relations. In Appendix A, such equations for the different parameter arrays are given. These equations are proved by using the method developed in [9] for the Racah case. Below, we recall this method for the Krawtchouk case. We want now that $c_{j+1}(\Psi_x)$ and $b_{j-1}(\Psi_x)$ correspond to $\lambda_x^+(j)$ and $\lambda_x^-(j)$ as follows

$$c_{j+1}(\Psi_x) = z_j^+ \lambda_x^+(\Phi_j), \qquad (4.25)$$

$$b_{j-1}(\Psi_x) = z_j^- \lambda_x^-(\Phi_j), \qquad (4.26)$$

for some functions $z_j^{\pm}$ (independent of $x$). Then the parameters $C_{i,j}^{\epsilon, \pm}$ can be computed

$$C_{i,j}^{\epsilon, \pm} = z_j^{\pm} \gamma_i^{\epsilon, \pm}(\Phi_j), \qquad (4.27)$$

and the constraints (4.17) are completely solved.

The same procedure is also done to solve the constraints (4.19) and compute $D_{x,y}^{\epsilon, \epsilon'}$.

As just explained, it is important to obtain contiguity relations (4.24) for the considered polynomials. Here, we prove the contiguity recurrence equation (A.41a) for the Krawtchouk polynomial. All the other contiguity equations given in Appendix A are proven similarly.

Relation (A.41a) is proven by induction. For $i = 1$, it is checked by direct computation. For a given $i \ge 1$, we suppose that the equation is valid for any $i' \le i$. From the recurrence relation (2.15a), one gets

$$p_{x,i+1}(\Phi) = \frac{1}{c_{i+1}(\Phi)} \left( (\theta_x - a_i(\Phi))p_{x,i}(\Phi) - b_{i-1}(\Phi)p_{x,i-1}(\Phi) \right). \qquad (4.28)$$

From the induction hypothesis, one may replace $p_{x,i}(\Phi)$ and $p_{x,i-1}(\Phi)$ to get

$$p_{x,i+1}(\Phi) = \frac{1}{c_{i+1}(\Phi)} \left( (\theta_x - a_i(\Phi))(p_{x,i}^{-}(\Phi) + \frac{r}{ss^* - r}p_{x,i-1}^{-}(\Phi)) - b_{i-1}(\Phi)(p_{x,i-1}^{-}(\Phi) + \frac{r}{ss^* - r}p_{x,i-2}^{-}(\Phi)) \right), \qquad (4.29)$$

where we recall that the superscript $-$ of $p_{x,i}^{-}(\Phi)$ means that $d$ is replaced by $d-1$ in its expression. We use again twice the recurrence relation (2.15a) for $\theta_x p_{x,i}^{-}(\Phi)$ and $\theta_x p_{x,i-1}^{-}(\Phi)$ to obtain

$$p_{x,i+1}(\Phi) = \frac{1}{c_{i+1}(\Phi)} \left( b_{i-1}^{-}(\Phi) p_{x,i-1}^{-}(\Phi) + a_i^{-}(\Phi) p_{x,i}^{-}(\Phi) + c_{i+1}^{-}(\Phi) p_{x,i+1}^{-}(\Phi) \right. \\ + \frac{r}{ss^* - r} \left( b_{i-2}^{-}(\Phi) p_{x,i-2}^{-}(\Phi) + a_{i-1}^{-}(\Phi) p_{x,i-1}^{-}(\Phi) + c_i^{-}(\Phi) p_{x,i}^{-}(\Phi) \right) \\ - a_i(\Phi) \left( p_{x,i}^{-}(\Phi) + \frac{r}{ss^* - r} p_{x,i-1}^{-}(\Phi) \right) \\ \left. - b_{i-1}(\Phi) \left( p_{x,i-1}^{-}(\Phi) + \frac{r}{ss^* - r} p_{x,i-2}^{-}(\Phi) \right) \right), \qquad (4.30)$$

where as for $p_{x,i}$, the superscript $-$ of $b$ and $c$ means that $d$ is replaced by $d-1$ in their expression. Then, using the explicit expressions (A.40) for $a_i$, $b_i$ and $c_i$, one proves that relation (A.41a) is also satisfied for $i+1$ which concludes the proof of its validity for any $i$.

## 4.5 Bilinear form

If, for a choice of parameter arrays $\Phi_j$ and $\Psi_x$, the constraints (4.17) and (4.19) are solved, then we want to find sufficient conditions on the parameter arrays for statement (vii) of Definition 3.1 (concerning the bilinear form) to hold.
Let us choose a bilinear form such that

$$ \langle V_{x,y}, V_{u,v} \rangle = \delta_{x,u} \delta_{y,v} F(x, y), \qquad (4.31) $$

for a non-vanishing function $F(x, y)$.

The change of basis (4.13) allows us to compute

$$ \langle V_{i,j}^*, V_{k,l}^* \rangle = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} p_{xi}(\Phi_j) p_{yj}(\Psi_x) \sum_{u=0}^{N-\ell} \sum_{v=0}^{N-u} p_{uk}(\Phi_\ell) p_{v\ell}(\Psi_u) \langle V_{x,y}, V_{u,v} \rangle, \qquad (4.32) $$

$$ = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} p_{xi}(\Phi_j) p_{yj}(\Psi_x) p_{xk}(\Phi_\ell) p_{y\ell}(\Psi_x) F(x, y). \qquad (4.33) $$

The Wilson duality (see for example [43]) states that

$$ p_{ij}(\Phi) = \frac{k_j(\Phi)}{k_i(\Phi^*)} p_{ji}(\Phi^*), \qquad (4.34) $$

where $k_i(\Phi)$ is given in (2.13), which allows us to rewrite (4.33) as follows

$$ \langle V_{i,j}^*, V_{k,l}^* \rangle = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} \frac{k_i(\Phi_j) k_j(\Psi_x)}{k_x(\Phi_j^*) k_y(\Psi_x^*)} p_{ix}(\Phi_j^*) p_{jy}(\Psi_x^*) p_{xk}(\Phi_\ell) p_{y\ell}(\Psi_x) F(x, y). \qquad (4.35) $$

If there exists $G(i, j)$, a function which does not depend on $x$ and $y$, such that

$$ \frac{k_i(\Phi_j) k_j(\Psi_x)}{k_x(\Phi_j^*) k_y(\Psi_x^*)} F(x, y) = \frac{G(i, j)}{\nu(\Psi_x) \nu(\Phi_j)}, \qquad (4.36) $$

then, using two times the orthogonality relation (2.16), one arrives at:

$$ \langle V_{i,j}^*, V_{k,l}^* \rangle = \frac{G(i, j)}{\nu(\Phi_j)} \sum_{x=0}^{N-j} \frac{1}{\nu(\Psi_x)} p_{ix}(\Phi_j^*) p_{xk}(\Phi_\ell) \sum_{y=0}^{N-x} p_{jy}(\Psi_x^*) p_{y\ell}(\Psi_x) \qquad (4.37) $$

$$ = \delta_{i,k} \delta_{j,l} G(i, j). \qquad (4.38) $$

Finally, the statement (vii) of Definition 3.1 of a factorized Leonard pair is satisfied if the constraint (4.36) is satisfied.

To conclude this section, let us summarize the results of this section in the following theorem.

**Theorem 4.3.** The parameters arrays $\Phi_j$ and $\Psi_x$ given by

$$ \Phi_j = (t_x, t_x^*(j), x = 0, \dots, N-j; \quad f_x(j), g_x(j), x = 1, \dots, N-j), \qquad (4.39a) $$

$$ \Psi_x = (\theta_j(x), \theta_j^*, j = 0, \dots, N-x; \quad \varphi_j(x), \phi_j(x), j = 1, \dots, N-x), \qquad (4.39b) $$

are called compatible if there exist $C_{i,j}^{\epsilon, \epsilon'}$ and $D_{x,y}^{\epsilon, \epsilon'}$ such that relations (4.17) and (4.19) are satisfied and the couple of functions $F(x, y)$ and $G(i, j)$ exists such that (4.36) is valid.

A factorized Leonard pair can be associated to a compatible pair of parameter arrays, and the coefficients $C_{i,j}^{\epsilon, \epsilon'}$ and $D_{x,y}^{\epsilon, \epsilon'}$ are uniquely determined. The actions of the elements of the basis $(X, Y; X^*, Y^*)$ of this factorized Leonard pair on the basis $V_{x,y}$ (statement (iii) of Definition 3.1) are given by

$$ X V_{x,y} = t_x V_{x,y}, \quad Y V_{x,y} = \theta_y(x) V_{x,y}, \qquad (4.40a) $$

$$ X^* V_{x,y} = \sum_{(\epsilon, \epsilon') \in A_2} D_{x,y}^{\epsilon, \epsilon'} V_{x+\epsilon, y+\epsilon'}, \qquad (4.40b) $$

$$ Y^* V_{x,y} = b_{y-1}(\Psi_x^*) V_{x,y-1} + a_y(\Psi_x^*) V_{x,y} + c_{y+1}(\Psi_x^*) V_{x,y+1}. \qquad (4.40c) $$
Similarly, on the basis $V_{i,j}^*$ (statement (iv) of Definition 3.1), the actions are given by

$$XV_{i,j}^* = b_{i-1}(\Phi_j)V_{i-1,j}^* + a_i(\Phi_j)V_{i,j}^* + c_{i+1}(\Phi_j)V_{i+1,j}^*, \quad (4.41a)$$

$$YV_{i,j}^* = \sum_{(\epsilon, \epsilon') \in A_2} C_{i,j}^{\epsilon, \epsilon'} V_{i+\epsilon, j+\epsilon'}^*, \quad (4.41b)$$

$$X^*V_{i,j}^* = t_i^*(j)V_{i,j}^*, \quad Y^*V_{i,j}^* = \theta_j^*V_{i,j}^*. \quad (4.41c)$$

The change of basis between both bases is provided by (4.13).

# 5 Polynomials of Tratnik type

In this section, let $\Phi_j$ and $\Psi_x$ be compatible parameter arrays. Before giving a list of possible factorized Leonard pairs, let us discuss the bivariate polynomials

$$T_{i,j}(x, y) = p_{xi}(\Phi_j)p_{yj}(\Psi_x), \quad (5.1)$$

which are the entries of the change of basis between $V_{x,y}$ and $V_{i,j}^*$ given in (4.13). These bivariate functions are called of Tratnik type in reference to [46] where such polynomials have been introduced.

We shall show that these bivariate polynomials are bispectral, i.e. they satisfy two recurrence and two difference relations, and that they are orthogonal.

## 5.1 Difference relations

From the actions of $X^*$ and $Y^*$ on the vectors $V_{i,j}^*$, one deduces the difference relations for $T_{i,j}(x, y)$. Indeed, using the results of Section 4, the action of $X^*$ on $V_{i,j}^*$ can be computed in two ways:

$$X^*V_{i,j}^* = t_i^*(j)V_{i,j}^* = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} T_{i,j}(x, y)t_i^*(j)V_{x,y}, \quad (5.2)$$

and

$$X^*V_{i,j}^* = \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} T_{i,j}(x, y) X^*V_{x,y} \quad (5.3)$$

$$= \sum_{x=0}^{N-j} \sum_{y=0}^{N-x} T_{i,j}(x, y) \sum_{(\epsilon, \epsilon') \in A_2} D_{x,y}^{\epsilon, \epsilon'} V_{x+\epsilon, y+\epsilon'}. \quad (5.4)$$

Comparing the coefficients in front of $V_{x,y}$, one obtains the first difference relation

$$t_i^*(j)T_{i,j}(x, y) = \sum_{(\epsilon, \epsilon') \in A_2} D_{x-\epsilon, y-\epsilon'}^{\epsilon, \epsilon'} T_{i,j}(x - \epsilon, y - \epsilon'). \quad (5.5a)$$

Similarly, the action of $Y^*$ on $V_{i,j}^*$ gives the second difference relation

$$\theta_j^* T_{i,j}(x, y) = b_y(\Psi_x^*) T_{i,j}(x, y+1) + a_y(\Psi_x^*) T_{i,j}(x, y) + c_y(\Psi_x^*) T_{i,j}(x, y-1). \quad (5.5b)$$

## 5.2 Recurrence relations

To obtain the recurrence relations, let us first invert the change of basis (4.13):

$$V_{x,y} = \sum_{j=0}^{N-x} \sum_{i=0}^{N-j} \frac{1}{\nu(\Phi_j)\nu(\Psi_x)} p_{ix}(\Phi_j^*) p_{yj}(\Psi_x^*) V_{i,j}^*. \quad (5.6)$$
Using the Wilson duality (4.34), one gets

$$V_{x,y} = \sum_{j=0}^{N-x} \sum_{i=0}^{N-j} \frac{k_x(\Phi_j^*) k_y(\Psi_x^*)}{\nu(\Psi_x) \nu(\Phi_j) k_i(\Phi_j) k_j(\Psi_x)} T_{i,j}(x,y) V_{i,j}^*. \quad (5.7)$$

Relation (4.36) allows us to simplify the latter relation:

$$V_{x,y} = \sum_{j=0}^{N-x} \sum_{i=0}^{N-j} \frac{F(x,y)}{G(i,j)} T_{i,j}(x,y) V_{i,j}^*. \quad (5.8)$$

One obtains recurrence relations for $T_{i,j}(x,y)$ from the actions of $X$ and $Y$ on the vectors $V_{x,y}$. Indeed, the action of $X$ on $V_{x,y}$ provides the first recurrence relation

$$t_x T_{i,j}(x,y) = \frac{G(i,j)}{G(i+1,j)} b_i(\Phi_j) T_{i+1,j}(x,y) + a_i(\Phi_j) T_{i,j}(x,y) + \frac{G(i,j)}{G(i-1,j)} c_i(\Phi_j) T_{i-1,j}(x,y), \quad (5.9a)$$

and the action of $Y$ on $V_{x,y}$ provides the second one

$$\theta_y(x) T_{i,j}(x,y) = \sum_{(\epsilon, \epsilon') \in A_2} \frac{G(i,j)}{G(i-\epsilon, j-\epsilon')} C_{i-\epsilon, j-\epsilon'}^{\epsilon, \epsilon'} T_{i-\epsilon, j-\epsilon'}(x,y). \quad (5.9b)$$

## 5.3 Degree of the polynomials

Knowing that $p_{0x}(\Phi) = 1$ for any parameter array $\Phi$ (see [43]), one deduces that $T_{0,0}(x,y) = 1$. Then, using the recurrence relations (5.9a) and (5.9b) for $i, j = 0$, one shows that $T_{1,0}(x,y)$ and $T_{0,1}(x,y)$ are polynomials with respect to $t_x$ and $\theta_y(x)$ of total degree 1.

Suppose that, for a given $k$ ($0 < k < N$) and for $i + j < k$, $T_{i,j}(x,y)$ are bivariate polynomials with respect to $t_x$ and $\theta_y(x)$ of total degree $i + j$. Using the recurrence relation (5.9a) for $(i, j) = (k-1, 0)$, $(k-2, 1)$, ..., $(0, k-1)$, one can compute $T_{i,k-i}(x,y)$ for $i = 1, 2, ... k$ in terms of $T_{i,j}(x,y)$ with $i + j < k$ and deduce that they are of total degree $k$. Then, using (5.9b) for $(i, j) = (0, k-1)$, one deduces that $T_{0,k}(x,y)$ is also of total degree $k$ in terms of $t_x$ and $\theta_y(x)$. Therefore, one shows recursively with this method that $T_{i,j}(x,y)$ is a polynomial with respect to $t_x$ and $\theta_y(x)$ of total degree $i + j$, for any couple $(i, j)$.

In fact, it is possible to be a bit more precise. Following the same recursive proof as previously, one can also show that the degree of $T_{i,j}(x,y)$ with respect to $\theta_y(x)$ is exactly $j$.

## 5.4 Orthogonality relation for the Tratnik polynomials

The functions $T_{ij}(x,y)$ are also orthogonal in addition of being bispectral. Indeed, the vectors $V_{x,y}$ can be expressed in terms of the vectors $V_{i,j}^*$ thanks to (5.8), and the latter can in turn be replaced using (4.13). This leads to the expression

$$V_{x,y} = \sum_{x'=0}^{N} \sum_{y'=0}^{N-x'} \sum_{j=0}^{\min(N-x, N-x')} \sum_{i=0}^{N-j} \frac{F(x,y)}{G(i,j)} T_{i,j}(x,y) T_{i,j}(x',y') V_{x',y'}. \quad (5.10)$$

Identifying the coefficient in front of $V_{x,y}$, one gets the orthogonality relation

$$\sum_{j=0}^{\min(N-x, N-x')} \sum_{i=0}^{N-j} \frac{1}{G(i,j)} T_{i,j}(x,y) T_{i,j}(x',y') = \delta_{x,x'} \delta_{y,y'} \frac{1}{F(x,y)}. \quad (5.11)$$

Now, starting by expressing $V_{i,j}^*$ in terms of $V_{x,y}$ thanks to (4.13), and then replacing the latter using (5.8), one gets

$$\sum_{x=0}^{\min(N-j, N-j')} \sum_{y=0}^{N-x} F(x,y) T_{i,j}(x,y) T_{i',j'}(x,y) = \delta_{i,i'} \delta_{j,j'} G(i,j). \quad (5.12)$$
# 6 List of factorized $A_2$-Leonard pairs of type II

In the following sections, we provide couples of compatible parameter arrays, for $j = 0, \dots, N$,

$$ \Phi_j = (t_x, t_x^*(j), x = 0, \dots, N-j; f_x(j), g_x(j), x = 1, \dots, N-j), \quad (6.1a) $$

and, for $x = 0, \dots, N$,

$$ \Psi_x = (\theta_j(x), \theta_j^*, j = 0, \dots, N-x; \varphi_j(x), \phi_j(x), j = 1, \dots, N-x). \quad (6.1b) $$

Then, using Theorem 4.3, this list also gives examples of factorized $A_2$-Leonard pairs. We use the $A_2$-contiguity recurrence equations given in Appendix A to prove that the constraints (4.17) and (4.19) hold. We also verify that (4.36) is satisfied. We shall use the invariance (4.2) to simplify these parameter arrays. The procedure is detailed in the example of Section 6.1.

We restrict in this section to the case where the parameter arrays are of type II. We recall that type II means that the (dual) eigenvalues can be written as $z_2x^2 + z_1x + z_0$, for some parameters $z_i$. Let us also recall that, for one type of eigenvalues, there exist different kinds of polynomials. One names a parameter array according to this associated polynomial. The following conjecture provides a strong restriction on the choice of the parameter arrays.

**Conjecture 6.1.** To get a factorized Leonard pair, it is necessary that $\Phi_j$ for any $j$ (resp. $\Psi_x$ for any $x$) is associated to the same kind of polynomial.

All the examples satisfy this conjecture and we name a pair of compatible parameter arrays $(\Phi_j, \Psi_x)$ with the two names of the associated polynomials.

We do not discuss the irreducibility of each Leonard pair. We refer to [43] where conditions on their parameters are explicitly given such that the associated Leonard pair is irreducible.

## 6.1 Hahn / dual Hahn

Let us consider the case where the parameter array $\Phi_j$ is of Hahn kind and $\Psi_x$ is of dual Hahn kind. The parameter array $\Phi_j$ ($j = 0, 1, \dots, N$) reads as follows, for $x = 0, 1, \dots, N-j$,

$$ t_x = x, \quad t_x^*(j) = m^*x(x+1+s^*+j), \quad (6.2a) $$

and, for $x = 1, \dots, N-j$,

$$ f_x(j) = m^*x(x - N + j - 1)(x + r), \quad g_x(j) = -m^*x(x - N + j - 1)(x + s^* + j - r). \quad (6.2b) $$

The parameter array $\Psi_x$ ($x = 0, 1, \dots, N$) reads as follows, for $j = 0, 1, \dots, N-x$,

$$ \theta_j(x) = \mu j(j + 1 + \sigma + x), \quad \theta_j^* = j, \quad (6.3a) $$

and, for $j = 1, \dots, N-x$,

$$ \varphi_j(x) = \mu j(j - N + x - 1)(j + \rho), \quad \phi_j(x) = \mu j(j - N + x - 1)(j + \rho - \sigma - N - 1). \quad (6.3b) $$

The parameter array $\Phi_j$ is obtained from the one given in (A.30), by taking

$$ \theta_0 \to 0, \quad \theta_0^* \to 0, \quad \sigma \to 1, \quad \mu^* \to m^*, \quad \sigma^* \to s^* + j, \quad \rho \to r, \quad d \to N - j, \quad (6.4) $$

and renaming $x \to i$, $\theta \to t$, $\varphi \to f$, $\phi \to g$. The dependence of $\Phi_j$ on $j$ is therefore present in the diameter $d$ and as a shift in the parameter $\sigma^*$. The parameters $\theta_0$ and $\theta_0^*$ are sent to 0 and $\sigma$ to 1 without loss of generality due to the invariance (4.2). Similarly, the parameter array $\Psi_x$ is obtained from (A.34) by taking

$$ \theta_0 \to 0, \quad \theta_0^* \to 0, \quad \sigma^* \to 1, \quad \sigma \to \sigma + x, \quad d \to N - x, \quad (6.5) $$
and renaming $i \to j$. The dependence of $\Psi_x$ on $x$ is present in the diameter $d$ and as a shift in the parameter $\sigma$.

The recurrence coefficients corresponding to the parameter array $\Phi_j$ are obtained from (A.32):

$$a_i(\Phi_j) = -b_i(\Phi_j) - c_i(\Phi_j), \qquad (6.6a)$$

$$b_i(\Phi_j) = \frac{(i + j - N)(i + 1 + r)(i + 1 + s^* + j)}{(2i + j + s^* + 2)(2i + j + s^* + 1)}, \qquad (6.6b)$$

$$c_i(\Phi_j) = -\frac{i(i + j + s^* - r)(i + 1 + s^* + N)}{(2i + j + s^*)(2i + j + s^* + 1)}, \qquad (6.6c)$$

and those corresponding to the parameter array $\Psi_x$ are obtained from (A.36):

$$a_j(\Psi_x) = -b_j(\Psi_x) - c_j(\Psi_x), \qquad (6.6d)$$

$$b_j(\Psi_x) = \mu(j - N + x)(j + 1 + \rho), \qquad (6.6e)$$

$$c_j(\Psi_x) = \mu j(j - N - 1 + \rho - \sigma). \qquad (6.6f)$$

The coefficient $a_j(\Psi_x)$ can be written as

$$a_j(\Psi_x) = -\mu(j + 1 + \rho)t_x + \mu(N - j)(j + 1 + \rho) + \mu j(N - j + 1 + \sigma - \rho). \qquad (6.7)$$

The constraint (4.17b) can therefore be solved using the recurrence relation (2.15a) of $p_{xi}(\Phi_j)$, giving

$$C_{i,j}^{-,0} = -\mu(j + 1 + \rho)b_{i-1}(\Phi_j), \qquad (6.8a)$$

$$C_{i,j}^{+,0} = -\mu(j + 1 + \rho)c_{i+1}(\Phi_j), \qquad (6.8b)$$

$$C_{i,j}^{0,0} = -\mu(j + 1 + \rho)a_i(\Phi_j) + \mu(N - j)(j + 1 + \rho) + \mu j(N - j + 1 + \sigma - \rho). \qquad (6.8c)$$

The contiguity relation (A.33a) for $p_{xi}(\Phi_j)$ allows one to solve the constraint (4.17a):

$$C_{i,j}^{-,+} = \frac{\mu(j + 1)(j - N + \rho - \sigma)(i + r)(s^* + j + 2)}{(2i + s^* + j)(s^* + j + 1 - r)}, \qquad (6.8d)$$

$$C_{i,j}^{0,+} = \frac{\mu(j + 1)(j - N + \rho - \sigma)(i + 1 + s^* + j - r)(s^* + j + 2)}{(2i + s^* + j + 2)(s^* + j + 1 - r)}. \qquad (6.8e)$$

The constraint (4.17c) is solved using the contiguity relation (A.33b):

$$C_{i,j}^{0,-} = \frac{\mu(j + \rho)(N - j + 1 - i)(i + s^* + j)(r - s^* - j)}{(2i + s^* + j)(s^* + j + 1)}, \qquad (6.8f)$$

$$C_{i,j}^{+, -} = \frac{\mu(j + \rho)(i + 1)(i + s^* + N + 2)(r - s^* - j)}{(2i + s^* + j + 2)(s^* + j + 1)}. \qquad (6.8g)$$

In order to determine the coefficients $D_{x,y}^{\epsilon, \epsilon'}$ which appear in the constraints (4.19), one can observe using (2.8) that $\Phi_j^*$ is the same as $\Psi_x$ with the replacements $j \to x$, $x \to j$, $\mu \to m^*$, $\sigma \to s^*$, $\rho \to r$, and $\Psi_x^*$ is the same as $\Phi_j$ with the replacements $j \to x$, $x \to j$, $m^* \to \mu$, $s^* \to \sigma$, $r \to \rho$. Therefore, comparing the constraints (4.17) and (4.19), one deduces that

$$D_{x,y}^{\epsilon, \epsilon'} = C_{y,x}^{\epsilon', \epsilon} \Big|_{(s^*, r, \mu, \sigma, \rho) \to (\sigma, \rho, m^*, s^*, r)} \times \frac{\nu(\Psi_{x+\epsilon})}{\nu(\Psi_x)}. \qquad (6.9)$$

Using (2.17), one finds

$$\frac{\nu(\Psi_{x+1})}{\nu(\Psi_x)} = \frac{\sigma - \rho + x + 1}{\sigma + x + 2}, \qquad \frac{\nu(\Psi_{x-1})}{\nu(\Psi_x)} = \frac{\sigma + x + 1}{\sigma - \rho + x}. \qquad (6.10)$$
The constraint (4.36) is also satisfied and one gets

$$F(x, y) = \binom{N}{x, y} \frac{(-1)^{N-y}(r+1)_x (\rho+1)_y (\rho-N-\sigma)_{N-x-y}}{(\sigma+x+y+1)_y (\sigma+x+2y+2)_{N-x-y}(r-s^*-N)_x}, \quad (6.11)$$

$$G(i, j) = \binom{N}{i, j} \frac{(r+1)_i [(2+s^*+j)_{N-j}]^2 (\rho+1)_j}{(s^*+j-r+1, s^*+j+1+i)_i (s^*+j+2i+2)_{N-j-i} (\rho-\sigma-N)_j (r-s^*-N)_{N-j}}, \quad (6.12)$$

where $\binom{N}{i,j} = \frac{N!}{(N-i-j)!i!j!}$ and $(x_1, \dots, x_j)_i$ corresponds to the Pochhammer symbols:

$$(x_1, \dots, x_j)_i = \prod_{k=0}^{i-1} (x_1 + k) \dots (x_j + k). \quad (6.13)$$

Therefore, by Theorem 4.3, we conclude that the pair $(\Phi_j, \Psi_x)$ with the parameter arrays (6.2) and (6.3) is associated to a factorized $A_2$-Leonard pair.

The bivariate polynomials (5.1) associated to this factorized Leonard pair are obtained from the (univariate) Hahn polynomials (A.31) and dual Hahn polynomials (A.35). They read explicitly as follows, for $i, j, x, y$ non-negative integers, $i + j \le N$, $x + y \le N$,

$$T_{i,j}(x, y) = \binom{N-j}{i} \binom{N-x}{j} \frac{(-1)^j (r+1)_i (2+s^*+j)_{N-j} (\rho+1)_j}{(s^*+j-r+1, s^*+j+1+i)_i (s^*+j+2i+2)_{N-j-i} (\rho-\sigma-N)_j} \\ \times {}_3F_2 \left( \begin{matrix} -i, i+1+s^*+j, -x \\ r+1, -N+j \end{matrix} \middle| 1 \right) {}_3F_2 \left( \begin{matrix} -j, -y, y+1+\sigma+x \\ \rho+1, -N+x \end{matrix} \middle| 1 \right). \quad (6.14)$$

In the previous equation, ${}_3F_2$ corresponds to the usual hypergeometric function (see e.g. [30]). As explained in Section 5.3, $T_{i,j}(x, y)$, given in the previous relation, are bivariate polynomials of variables $x$ and $\mu y(y+1+\sigma+x)$ of total degree $i+j$. Let us emphasize that the functions obtained here are not the ones studied previously for example in [46, 21]. The latter are not obtained here since one of the four bispectral relations are more complicated and involve nine terms (called $B_2$ in the outlooks of the article).

## 6.2 Krawtchouk / dual Hahn

Let us consider the case Krawtchouk / dual Hahn where $\Phi_j$ is given by, for $j = 0, 1, \dots, N$,

$$t_x = x, \quad t_x^*(j) = x, \quad (6.15a)$$

$$f_x(j) = rx(x - N + j - 1), \quad g_x(j) = (r-1)x(x - N + j - 1), \quad (6.15b)$$

and, $\Psi_x$ is given by (6.3).

The coefficients $C_{i,j}^{\epsilon, \epsilon'}$ can be computed from (A.41a)-(A.41b) and $D_{x,y}^{\epsilon, \epsilon'}$, from (A.33a)-(A.33b). Alternatively, one can take the results of Subsection 6.1, perform the replacements $s^* \to \frac{1}{m^*}$, $r \to \frac{r}{m^*}$, and then take the limit $m^* \to 0$.

## 6.3 Hahn / Krawtchouk

Let us consider the case Hahn / Krawtchouk where $\Phi_j$ is given by (6.2) and $\Psi_x$ is given by, for $x = 0, 1, \dots, N$,

$$\theta_j(x) = j, \quad \theta_j^* = j, \quad (6.16a)$$

$$\varphi_j(x) = \rho j(j - N + x - 1), \quad \phi_j(x) = (\rho - 1)j(j - N + x - 1). \quad (6.16b)$$

The coefficients $C_{i,j}^{\epsilon, \epsilon'}$ can be computed as previously from (A.33a)-(A.33b) and $D_{x,y}^{\epsilon, \epsilon'}$, from (A.41a)-(A.41b). Alternatively, one can take the results of Subsection 6.1, perform the replacements $\sigma \to \frac{1}{\mu}$, $\rho \to \frac{\rho}{\mu}$, and then take the limit $\mu \to 0$.

The functions of Tratnik type for the last two subsections (obtained as limit of the ones given in (6.14)) have been used in [38, 3, 13] to study the non-binary Johnson association scheme.
## 6.4 Krawtchouk / Krawtchouk

Let us consider the case Krawtchouk / Krawtchouk where $\Phi_j$ is given by (6.15) and $\Psi_x$ is given by (6.16).

The coefficients $C_{i,j}^{\epsilon, \epsilon'}$ and $D_{x,y}^{\epsilon, \epsilon'}$ can be computed from (A.41a)-(A.41b). Alternatively, one can take the results of Subsection 6.1, perform the replacements $s^* \to \frac{1}{m^*}$, $r \to \frac{r}{m^*}$, $\sigma \to \frac{1}{\mu}$, $\rho \to \frac{\rho}{\mu}$, and then take the limits $m^* \to 0$, $\mu \to 0$.

The polynomial of Tratnik type in this case can be obtained as a limit of (6.14). This polynomial has already appeared in different contexts and it is a particular case of the Griffiths (or sometimes called Rahman) polynomials [35, 48, 25, 26, 32, 27, 19, 11, 7]. Indeed, up to a normalisation and specifying some parameters, the Griffiths polynomials can be written as a product of two Krawtchouk polynomials and take the form (5.1).

# 7 List of factorized $A_2$-Leonard pairs of type I

In the following, we provide couples of parameter arrays $(\Phi_j, \Psi_x)$ of type I, such that the constraints (4.17), (4.19) and (4.36) can be solved. We recall that type I means that the (dual) eigenvalues can be written as $z^+q^x + z^-q^{-x} + z^0$, for some parameters $z^\pm, z^0$. These provide examples of factorized $A_2$-Leonard pairs of type I.

## 7.1 $q$-Hahn / dual $q$-Hahn

One can consider for $\Phi_j$ the parameter array of the $q$-Hahn kind

$$t_x = q^{-x} - 1, \quad t_x^*(j) = m^* q^{N-j} (q^{-x} - 1) (1 - s^* q^{j+x+1}), \qquad (7.1a)$$

$$f_x(j) = m^* q^{N-j+1-2x} (1-q^x) (1-q^{x-N+j-1}) (1-rq^x), \qquad (7.1b)$$

$$g_x(j) = -m^* q^{N-j+1-x} (1-q^x) (1-q^{x-N+j-1}) (r-s^*q^{x+j}), \qquad (7.1c)$$

and for $\Psi_x$ one takes the parameter array of the dual $q$-Hahn kind

$$\theta_j(x) = \mu q^{N-x} (q^{-j} - 1) (1 - \sigma q^{x+j+1}), \quad \theta_j^* = q^{-j} - 1, \qquad (7.2a)$$

$$\varphi_j(x) = \mu q^{N-x+1-2j} (1-q^j) (1-q^{j-N+x-1}) (1-\rho q^j), \qquad (7.2b)$$

$$\phi_j(x) = \mu q^{2N-2x+2-2j} (1-q^j) (1-q^{j-N+x-1}) (\sigma q^x - \rho q^{j-N+x-1}). \qquad (7.2c)$$

The recurrence coefficients for these parameter arrays are obtained from (A.4) and (A.8):

$$a_i(\Phi_j) = -b_i(\Phi_j) - c_i(\Phi_j), \qquad (7.3a)$$

$$b_i(\Phi_j) = \frac{(1 - s^* q^{i+j+1})(1 - q^{i-N+j})(1 - rq^{i+1})}{(1 - s^* q^{2i+j+2})(1 - s^* q^{2i+j+1})}, \qquad (7.3b)$$

$$c_i(\Phi_j) = \frac{q^{i-N+j}(1 - s^* q^{N+i+1})(1 - q^i)(s^* q^{i+j} - r)}{(1 - s^* q^{2i+j+1})(1 - s^* q^{2i+j})}, \qquad (7.3c)$$

and

$$a_j(\Psi_x) = -b_j(\Psi_x) - c_j(\Psi_x), \qquad (7.3d)$$

$$b_j(\Psi_x) = \mu q^{N-x} (1 - q^{j-N+x})(1 - \rho q^{j+1}), \qquad (7.3e)$$

$$c_j(\Psi_x) = \mu q^{N+1} (1 - q^j)(\sigma - \rho q^{j-N-1}). \qquad (7.3f)$$

One can write

$$a_j(\Psi_x) = -\mu q^N (1 - \rho q^{j+1}) t_x + \mu q^j (1 - q^{N-j})(1 - \rho q^{j+1}) + \mu q^j (1 - q^j)(\rho - \sigma q^{N-j+1}). \qquad (7.3g)$$
Using these recurrence coefficients, the recurrence relation (2.15a) and the contiguity recurrence relations (A.5a)–(A.5b) for the polynomials $p_{xi}(\Phi_j)$, one solves the constraints (4.17):

$$C_{i,j}^{-,0} = -\mu q^N (1 - \rho q^{j+1}) b_{i-1}(\Phi_j), \qquad (7.4a)$$

$$C_{i,j}^{0,0} = -\mu q^N (1 - \rho q^{j+1}) a_i(\Phi_j) + \mu q^j (1 - q^{N-j})(1 - \rho q^{j+1}) + \mu q^j (1 - q^j)(\rho - \sigma q^{N-j+1}), \quad (7.4b)$$

$$C_{i,j}^{+,0} = -\mu q^N (1 - \rho q^{j+1}) c_{i+1}(\Phi_j), \qquad (7.4c)$$

$$C_{i,j}^{-,+} = \frac{\mu q^N (1 - q^{j+1})(\sigma - \rho q^{j-N})(1 - r q^i)(1 - s^* q^{j+2})}{(1 - s^* q^{j+2i})(r - s^* q^{j+1})}, \qquad (7.4d)$$

$$C_{i,j}^{0,+} = \frac{\mu q^{N+i+1}(1 - q^{j+1})(\sigma - \rho q^{j-N})(r - s^* q^{j+i+1})(1 - s^* q^{j+2})}{(1 - s^* q^{j+2i})(r - s^* q^{j+1})}, \qquad (7.4e)$$

$$C_{i,j}^{0,-} = -\frac{\mu q^{i+j}(1 - \rho q^j)(1 - q^{N-j+1-i})(1 - s^* q^{i+j})(r - s^* q^j)}{(1 - s^* q^{2i+j})(1 - s^* q^{j+1})}, \qquad (7.4f)$$

$$C_{i,j}^{+, -} = -\frac{\mu q^j (1 - \rho q^j)(1 - q^{i+1})(1 - s^* q^{i+N+2})(r - s^* q^j)}{(1 - s^* q^{2i+j+2})(1 - s^* q^{j+1})}, \qquad (7.4g)$$

For similar reasons as in Subsection 6.1, the coefficients solving the constraints (4.19) are given by

$$D_{x,y}^{\epsilon, \epsilon'} = C_{y,x}^{\epsilon', \epsilon} \Big|_{(s^*, r, \mu, \sigma, \rho) \to (\sigma, \rho, m^*, s^*, r)} \times \frac{\nu(\Psi_{x+\epsilon})}{\nu(\Psi_x)}, \qquad (7.5)$$

where

$$\frac{\nu(\Psi_{x+1})}{\nu(\Psi_x)} = \frac{\sigma q^x - \rho q^{-1}}{\sigma q^x - q^{-2}}, \qquad \frac{\nu(\Psi_{x-1})}{\nu(\Psi_x)} = \frac{\sigma q^x - q^{-1}}{\sigma q^x - \rho}. \qquad (7.6)$$

The constraint (4.36) is satisfied with the following functions

$$F(x, y) = \begin{bmatrix} N \\ x, y \end{bmatrix}_q \frac{(-1)^N q^{y(x+y-N)} (rq; q)_x (\rho q; q)_y (\rho q^{-N} / \sigma; q)_{N-x-y}}{(qr)^x (q^{N-x+1} s^* / r; q)_x (\sigma q^{x+y+1}; q)_y (q^{-y-N-1} / \sigma; q)_{N-x-y}}, \qquad (7.7)$$

$$G(i, j) = \begin{bmatrix} N \\ i, j \end{bmatrix}_q \frac{q^{-\frac{N}{2}(N+3)} q^{j(j+2)} (rq; q)_i (\rho q; q)_j}{(rq)^i (\sigma q^{N+2})^j (s^* q^{j+1}/r, s^* q^{i+j+1}; q)_i (\rho q^{-N}/\sigma; q)_j} \\ \times \frac{(s^* q^{j+2}; q)_{N-j}^2}{(s^*)^{N-j} (rq^{-N}/s^*; q)_{N-j} (s^* q^{2i+j+2}; q)_{N-i-j}}, \qquad (7.8)$$

where we are using standard notations for the q-Pochhammer symbols and the q-trinomial coefficient:

$$(x; q)_i = \prod_{k=0}^{i-1} (1 - xq^k), \quad (x_1, \dots, x_j; q)_i = \prod_{k=1}^{j} (x_k; q)_i, \qquad (7.9)$$

$$\begin{bmatrix} N \\ i, j \end{bmatrix}_q = \frac{(q; q)_N}{(q; q)_i (q; q)_j (q; q)_{N-i-j}}. \qquad (7.10)$$

By Theorem 4.3, we conclude that the pair $(\Phi_j, \Psi_x)$ with the parameter arrays (7.1) and (7.2) is associated to a factorized $A_2$-Leonard pair.

The bivariate polynomials (5.1) associated to this factorized Leonard pair are obtained from the (univariate) q-Hahn polynomials (A.3) and dual q-Hahn polynomials (A.7), for $i, j, x, y$ non negative integers, $i + j \le N, x + y \le N$:

$$T_{i,j}(x, y) = \frac{(-1)^j q^{\frac{j}{2}(j-3-2N)} (rq; q)_i (s^* q^{j+2}; q)_{N-j} (\rho q; q)_j}{(rq)^i (s^* q^{j+1}/r, s^* q^{i+j+1}; q)_i (s^* q^{2+2i+j}; q)_{N-i-j} \sigma^j (q^{-N} \rho/\sigma; q)_j} \\ \times \begin{bmatrix} N-j \\ i \end{bmatrix}_q \begin{bmatrix} N-x \\ j \end{bmatrix}_q 3\phi_2 \left( \begin{array}{c} q^{-i}, \quad q^{-x}, \quad s^* q^{i+j+1} \\ rq, \quad q^{-N+j} \end{array} \middle| \begin{array}{c} q, \quad q \end{array} \right) 3\phi_2 \left( \begin{array}{c} q^{-j}, \quad q^{-y}, \quad \sigma q^{x+y+1} \\ \rho q, \quad q^{-N+x} \end{array} \middle| \begin{array}{c} q, \quad q \end{array} \right). \qquad (7.11)$$

As explained in Section 5.3, $T_{i,j}(x, y)$, given in the previous relation, are bivariate polynomials of variables $q^{-x} - 1$ and $\mu q^{N-x}(q^{-y} - 1)(1 - \sigma q^{x+y+1})$ of total degree $i + j$.
## 7.2 Limits

As for the factorized Leonard pairs of type II in Section 6, other examples of factorized Leonard pairs of type I can be obtained by taking limits of the case q-Hahn/dual q-Hahn described in Subsection 7.1. We provide here only the limiting procedures leading to couples $(\Phi_j, \Psi_x)$ of different kinds; it is then straightforward to deduce the explicit expressions for the parameter arrays $\Phi_j, \Psi_x$, the coefficients $C_{i,j}^{\epsilon, \epsilon'}$, $D_{x,y}^{\epsilon, \epsilon'}$, the functions $F(x, y)$, $G(i, j)$ and the polynomials $T_{i,j}(x, y)$ using the results of Subsection 7.1.

For the parameter array $\Phi_j$, one can take the following limits starting from the q-Hahn case given by (7.1):

* quantum q-Krawtchouk: $s^* \to \frac{q^{-N-1}}{m^*}$, $r \to \frac{1}{rm^*}$ and then $m^* \to 0$;

* affine q-Krawtchouk: $m^* \to 1$, $s^* \to 0$;

* q-Krawtchouk: $r \to 0$.

Similarly, for the parameter array $\Psi_x$, one can take the following limits starting from the dual q-Hahn case given by (7.2):

* dual quantum q-Krawtchouk: $\sigma \to \frac{q^{-N-1}}{\mu}$, $\rho \to \frac{1}{\rho\mu}$ and then $\mu \to 0$;

* affine q-Krawtchouk: $\mu \to 1$, $\sigma \to 0$;

* dual q-Krawtchouk: $\rho \to 0$.

Any combination of the above limits for $\Phi_j$ and $\Psi_x$ (including the case where no limit is taken for one or both of the parameter arrays) gives a factorized $A_2$-Leonard pair of type I.

The polynomial associated to the case quantum q-Krawtchouk / dual quantum q-Krawtchouk has been studied previously [46, 36, 37, 18]. In these previous papers (see e.g. [18]), it has been written as proportional to, for $n_1, n_2, x_1, x_2$ non negative integers with $n_1 + n_2 \le N$ and $x_1 + x_2 \le N$:

$$ \mathbf{K}_{n_1 n_2}(x_1, x_2) \propto {}_2\phi_1 \left( \begin{array}{c} \tilde{q}^{-n_1}, \tilde{q}^{-x_1} \\ \tilde{q}^{-x_1-x_2} \end{array} \middle| \tilde{q}, \frac{\tilde{q}^{n_1+1}}{\alpha_1^2} \right) {}_2\phi_1 \left( \begin{array}{c} \tilde{q}^{-n_2}, \tilde{q}^{-x_1-x_2+n_1} \\ \tilde{q}^{-N+n_1} \end{array} \middle| \tilde{q}, \frac{\tilde{q}^{n_2+1}}{\alpha_2^2} \right). \quad (7.12) $$

To recover the expression obtained in this paper, one performs the following changes:

$$ x_1 = y, \quad x_2 = N - y - x, \quad n_1 = j, \quad n_2 = N - i - j, \quad \alpha_1^2 = 1/r, \quad \alpha_2^2 = 1/\rho, \quad (7.13) $$

and, using relation (1.4.6) in [17], one rewrites the second factor to get

$$ \mathbf{K}_{j,N-i-j}(y, N-x-y) \propto {}_2\phi_1 \left( \begin{array}{c} \tilde{q}^{-j}, \tilde{q}^{-y} \\ \tilde{q}^{-N+x} \end{array} \middle| \tilde{q}, r\tilde{q}^{j+1} \right) {}_2\phi_1 \left( \begin{array}{c} \tilde{q}^{-i}, \tilde{q}^{-x} \\ \tilde{q}^{-N+j} \end{array} \middle| \tilde{q}, \rho\tilde{q}^{x+1} \right). \quad (7.14) $$

Then, setting $\tilde{q} = q^{-1}$ and using

$$ {}_2\phi_1 \left( \begin{array}{c} a_1, a_2 \\ b \end{array} \middle| q^{-1}, z \right) = {}_2\phi_1 \left( \begin{array}{c} 1/a_1, 1/a_2 \\ 1/b \end{array} \middle| q, \frac{za_1a_2}{qb} \right), \quad (7.15) $$

one obtains

$$ \mathbf{K}_{j,N-i-j}(y, N-x-y) \propto {}_2\phi_1 \left( \begin{array}{c} q^{-j}, q^{-y} \\ q^{-N+x} \end{array} \middle| q, rq^{j+x-N} \right) {}_2\phi_1 \left( \begin{array}{c} q^{-i}, q^{-x} \\ q^{-N+j} \end{array} \middle| q, \rho q^{x+j-N} \right). \quad (7.16) $$

In 7.16, we recognize our expression of (7.11) in the limit quantum q-Krawtchouk / dual quantum q-Krawtchouk.
## 8 Factorized $A_M$-Leonard pair

Up to now, we focused on the $A_2$-Leonard pairs. In this section, we propose the definition of $A_M$-Leonard pair for $M > 2$ which is, as for the case $M = 2$, a refinement of the definition of the rank $M$ Leonard pair (called here $A_M$-Leonard pair) given in [27].

Let $N, M$ be positive integers and let $\mathcal{D}$ be the set constituted of $M$-tuples of non-negative integers whose sum is at most $N$:

$$ \mathcal{D} = \{(x_1, x_2, \dots, x_M) \mid x_1, x_2, \dots, x_M \in \mathbb{Z}_{\ge 0}, x_1 + x_2 + \dots + x_M \le N \}. \quad (8.1) $$

Two $M$-tuples $\mathbf{x} = (x_1, x_2, \dots, x_M)$, $\mathbf{y} = (y_1, y_2, \dots, y_M) \in \mathcal{D}$ are called $A_M$-adjacent, $\mathbf{x}$ adj $\mathbf{y}$, if

$$ (x_1 - y_1, x_2 - y_2, \dots, x_M - y_M) \in A_M, \quad (8.2) $$

where $A_M$ is the following subset of $\mathbb{Z}^M$:

$$ A_M = \{(0, 0, \dots, 0), (1, -1, 0, \dots, 0), (1, 0, \dots, 0), (-1, 0, \dots, 0), \text{permutations} \}. \quad (8.3) $$

The notation $A_M$ for the previous set comes from the fact that it corresponds to the root system of the Weyl group $A_M$.

**Definition 8.1.** [27] Let $V$ be a $\mathbb{C}$-vector space. The pair $(H, H^*)$ is a $A_M$-Leonard pair on the domain $\mathcal{D}$ if the following statements are satisfied:

(i) $H$ is a $M$-dimensional subspace of $\text{End}(V)$ whose elements are diagonalizable and mutually commute;

(ii) $H^*$ is a $M$-dimensional subspace of $\text{End}(V)$ whose elements are diagonalizable and mutually commute;

(iii) there exists a bijection $\mathbf{x} \mapsto V_{\mathbf{x}}^{(M)}$ from $\mathcal{D}$ to the common eigenspaces of $H$ such that, for all $\mathbf{x} \in \mathcal{D}$,

$$ H^*V_{\mathbf{x}}^{(M)} \subseteq \sum_{\substack{\mathbf{y} \in \mathcal{D} \\ \mathbf{y} \text{ adj } \mathbf{x}}} V_{\mathbf{y}}^{(M)}; \quad (8.4) $$

(iv) there exists a bijection $\mathbf{x} \mapsto V_{\mathbf{x}}^{(0)}$ from $\mathcal{D}$ to the common eigenspaces $V_{\mathbf{x}}^{(0)}$ of $H^*$ such that, for all $\mathbf{x} \in \mathcal{D}$,

$$ HV_{\mathbf{x}}^{(0)} \subseteq \sum_{\substack{\mathbf{y} \in \mathcal{D} \\ \mathbf{y} \text{ adj } \mathbf{x}}} V_{\mathbf{y}}^{(0)}; \quad (8.5) $$

(v) there does not exist a subspace $W$ of $V$ such that $HW \subseteq W$, $H^*W \subseteq W$, $W \neq 0$, $W \neq V$;

$$ (\text{vi}) \dim(V_{\mathbf{x}}^{(0)}) = \dim(V_{\mathbf{x}}^{(M)}) = 1, \text{ for all } \mathbf{x} \in \mathcal{D}; $$

(vii) there exists a non-degenerate symmetric bilinear form $\langle \cdot, \cdot \rangle$ on $V$ such that, for $\mathbf{x}, \mathbf{y} \in \mathcal{D}$ and $\mathbf{x} \neq \mathbf{y}$,

$$ \langle V_{\mathbf{x}}^{(0)}, V_{\mathbf{y}}^{(0)} \rangle = 0, \quad \langle V_{\mathbf{x}}^{(M)}, V_{\mathbf{y}}^{(M)} \rangle = 0. \quad (8.6) $$

Let us notice that the dimension of the vector space $V$ is finite and is given more precisely by $\dim(V) = |\mathcal{D}|$.

As previously for the case $A_2$, we want to add some statements to this previous definition in order to obtain a definition that is easier to handle.

**Definition 8.2.** Let $(H, H^*)$ be a $A_M$-Leonard pair on the domain $\mathcal{D}$ and $X_1, X_2, \dots, X_M \in \text{End}(V)$ (resp. $X_1^*, X_2^*, \dots, X_M^*$) be a basis of $H$ (resp. $H^*$). The pair $(H, H^*)$ is a factorized $A_M$-Leonard pair if the following statements are satisfied:
(viii) For $a = 1, \dots, M-1$ and $\mathbf{x} = (x_1, \dots, x_M)$, the eigenvalues of $X_a$ on $V_{\mathbf{x}}^{(M)}$ do not depend on $x_{a+1}, x_{a+2}, \dots, x_M$:

$$ (X_a - \theta^{(a)}(x_1, \dots, x_a) \mathbf{1}) V_{\mathbf{x}}^{(M)} = 0, \qquad (8.7) $$

and $\theta^{(a)}(x_1, \dots, x_{a-1}, x_a) \neq \theta^{(a)}(x_1, \dots, x_{a-1}, x'_a)$ for $x_a \neq x'_a$.

For $a = 2, \dots, M$ and $\mathbf{i} = (i_1, \dots, i_M)$, the eigenvalues of $X_a^*$ on $V_{\mathbf{i}}^{(0)}$ do not depend on $i_1, i_2, \dots, i_{a-1}$:

$$ (X_a^* - \theta^{(a)}(i_a, \dots, i_M) \mathbf{1}) V_{\mathbf{i}}^{(0)} = 0, \qquad (8.8) $$

and $\theta^{(a)}(i_a, i_{a+1}, \dots, i_M) \neq \theta^{(a)}(i'_a, i_{a+1}, \dots, i_M)$ for $i_a \neq i'_a$;

(ix) for $1 \le i < j \le M$, one has

$$ [X_i, X_j^*] = 0; \qquad (8.9) $$

(x) let $1 \le a \le M-1$. There exists a bijection $\mathbf{x} \mapsto V_{\mathbf{x}}^{(a)}$ from $\mathcal{D}$ to the common eigenspaces of $X_1, X_2, \dots, X_a$ and $X_{a+1}^*, \dots, X_M^*$.

(xi) For a given $b \in \{1, 2, \dots, M-1\}$ and for any $a$ and $c$ satisfying $0 < a \le b < c \le M$, the following relations hold, for $(\dots, i_a, \dots), (\dots, i_a + 1, \dots) \in \mathcal{D}$,

$$ \langle V_{(...,i_a+1,...)}^{(a-1)}, X_a V_{(...,i_a,...)}^{(a-1)} \rangle \neq 0, \qquad \langle V_{(...,i_a,...)}^{(a-1)}, X_a V_{(...,i_a+1,...)}^{(a-1)} \rangle \neq 0, \qquad (8.10a) $$

and, for $(\dots, i_c, \dots), (\dots, i_c + 1, \dots) \in \mathcal{D}$,

$$ \langle V_{(...,i_c+1,...)}^{(c)}, X_c^* V_{(...,i_c,...)}^{(c)} \rangle \neq 0, \qquad \langle V_{(...,i_c,...)}^{(c)}, X_c^* V_{(...,i_c+1,...)}^{(c)} \rangle \neq 0. \qquad (8.10b) $$

The ordered tuple $(X_1, X_2, \dots, X_M; X_1^*, X_2^*, \dots, X_M^*)$ is called a basis of the factorized $A_M$-Leonard pair.

We believe that for a generic $M$ it is possible to obtain the same results as those obtained in this paper for the case $M=2$. In particular the polynomials of Tratnik type with $M$ variables may appear naturally. We plan to study the case with a generic $M$ in the future.

# 9 Outlooks

In this paper, we provide the definition of a factorized $A_2$-Leonard pair and start the study of these objects. We show that, in comparison to generic rank 2 Leonard pairs, it is possible to derive interesting properties for factorized Leonard pairs as a consequence of their refined definition, and it is possible to find several classes of examples based on combinations of polynomials of the $(q)$-Askey scheme. Regarding the examples, let us mention the following conjecture.

**Conjecture 9.1.** The list of factorized $A_2$-Leonard pairs given in Sections 6 and 7 is exhaustive.

In other words, we believe that a classification of all the factorized $A_2$-Leonard pairs is obtained. In particular, it seems that there is no factorized $A_2$-Leonard pairs involving the Bannai-Ito polynomials (of type III) or the $(q)$-Racah polynomials.

To prove the above conjecture, we intend to study in detail the algebra associated to a factorized $A_2$-Leonard pair, mimicking the proofs for the usual Leonard pairs. Indeed, as shown in [45], the two
members of a Leonard pair satisfy the Askey-Wilson relations. Using this result, we can show that, for a basis $(X, Y; X^*, Y^*)$ of a factorized $A_2$-Leonard pair, $X$ and $X^*$ also satisfy these relations i.e.

$$X^2X^* - \beta XX^*X + X^*X^2 - \gamma(XX^* + X^*X) - \rho X^* = \gamma^*X^2 + \omega X + \eta I, \qquad (9.1a)$$

$$(X^*)^2 X - \beta X^* X X^* + X(X^*)^2 - \gamma^*(XX^* + X^*X) - \rho^* X = \gamma(X^*)^2 + \omega X^* + \eta^* I, \qquad (9.1b)$$

where the parameters $\beta, \gamma, \gamma^*, \rho, \rho^*, \omega, \eta$ and $\eta^*$ depend on $Y^*$. We recall that $Y^*$ commutes with $X$ and $X^*$, meaning that $Y^*$ is a central element in the algebra generated by $X$ and $X^*$. Similarly, $Y$ and $Y^*$ also satisfy the Askey-Wilson relations:

$$Y^2Y^* - \hat{\beta}YY^*Y + Y^*X^2 - \hat{\gamma}(YY^* + Y^*Y) - \hat{\rho}Y^* = \hat{\gamma}^*Y^2 + \hat{\omega}Y + \hat{\eta}I, \qquad (9.2a)$$

$$(Y^*)^2Y - \hat{\beta}Y^*YY^* + Y(Y^*)^2 - \hat{\gamma}^*(YY^* + Y^*Y) - \hat{\rho}^*Y = \hat{\gamma}(Y^*)^2 + \hat{\omega}Y^* + \hat{\eta}^*I, \qquad (9.2b)$$

where the parameters $\hat{\beta}, \hat{\gamma}, \hat{\gamma}^*, \hat{\rho}, \hat{\rho}^*, \hat{\omega}, \hat{\eta}$ and $\hat{\eta}^*$ depend on $X$ (we recall that $X$ commutes with $Y$ and $Y^*$). It would be necessary to compute the explicit dependence of the parameters in terms of $Y^*$ (or $X$) but also to identify the relations between the four generators $X, Y, X^*$ and $Y^*$. In particular, $[X^*, Y]$ should satisfy some relations with the four previous generators. We hope to come back to this point in the future and pave the way to a proof of the above conjecture.

It would also be interesting to study the cases where the set $A_2$ is replaced by other root systems as it has been mentioned in [27]. The main other example is to replace $A_2$ by $B_2$:

$$B_2 = \{(0,0), (1,-1), (1,0), (0,1), (-1,1), (-1,0), (0,-1), (1,1), (-1,-1)\}. \qquad (9.3)$$

Polynomials of $(q)$-Racah kind should appear in this case. The higher rank factorized $A_M$-Leonard pairs introduced in Section 8 could also be investigated, as well as their analogs for different root systems. In a broader viewpoint, an important question would be to understand if the approach proposed in this paper also permits to study the higher rank Leonard pairs in general and their associated multivariate polynomials.

Another direction for generalizing the results of this paper is to consider other domains. Indeed, in this paper, we focus on the domain given by (3.1) but we believe that, fixing some free parameters to integer values, the domain can be restricted to some other forms. Results of this type have been obtained in [28] for Hahn multivariate polynomials.

Finally, we would like to mention that the bispectrality property of the Griffiths polynomials [24] or the Griffiths polynomials of Racah type [9] leads to examples of rank 2 Leonard pairs which are not factorized Leonard pairs. To our knowledge, these are the only examples of polynomials which are bispectral and not of Tratnik type. To study these polynomials, it would be necessary to generalize the notion of factorized Leonard pair. In Figure 1, we provide a schematic way to understand a factorized Leonard pair. By using the same graphical representation, we believe that the structure associated to the Griffiths type polynomials is represented in Figure 2. This schematic interpretation has been used in [9] as the path in the icosidodecahedron to get the Griffiths polynomials of Racah type (see the case $\partial = 3$ in Figure 3 of [9]). Let us remark that a supplementary operator, called $Z$, is necessary. We intend to start the study of this generalized factorized Leonard pair in the future.

<figure><img src="image_4.png" /><figcaption>Figure 2: Above each vertex, a couple of commuting operators is written and, below, the vectors of the basis where these operators are diagonal. At each oriented edge is associated the parameter array of the Leonard pair given by the operators above this edge.</figcaption></figure>
Acknowledgements:

The authors thank P. Iliev and P. Terwilliger for their valuable comments and LAPTh for its hospital-
ity. They are partially supported by the international research project AAPT of the CNRS. M. Zaimi
holds an Alexander-Graham-Bell scholarship from the Natural Sciences and Engineering Research
Council of Canada (NSERC).

A List of parameter arrays

In [43] all the possible parameter arrays

$$
\Phi = (\theta_i, \theta_i^*, i = 0, \dots, d; \varphi_i, \phi_i, i = 1, \dots, d) \qquad (A.1)
$$

are listed. We recall this list below (except for the cases Racah, q-Racah and Banna-Ito since we
do not find A₂-Leonard pairs based on them) with the associated polynomials pᵣᵢ(Φ) = kᵢ(Φ)uᵢ(θₓ)
(uᵢ(θₓ) is given by the (q-)hypergeometric function and kᵢ(Φ) is the factor in front of it) and the
associated coefficients bᵢ(Φ) and cᵢ(Φ) (we recall that aᵢ(Φ) = -bᵢ(Φ) - cᵢ(Φ) + θ₀). We also provide
for each parameter array A₂-contiguity relations.

A.1 List of parameter arrays of type I

A.1.1 q-Hahn

Parameter array.

$$
\begin{align*}
\theta_i &= \theta_0 + \mu(1 - q^i)q^{-i}, & \theta_i^* &= \theta_0^* + \mu^*(1 - q^i)(1 - \sigma^*q^{i+1})q^{-i}, \\
\varphi_i &= \mu\mu^*q^{1-2i}(1 - q^i)(1 - q^{i-d-1})(1 - \rho q^i), & \phi_i &= -\mu\mu^*q^{1-i}(1 - q^i)(1 - q^{i-d-1})(\rho - \sigma^*q^i). \tag{A.2}
\end{align*}
$$

Polynomial.

$$
p_{xi}(\Phi) = \frac{(q\rho; q)_i (q^2\sigma^*; q)_d}{(q\rho)^i (q\sigma^*/\rho, \sigma^*q^{i+1}; q)_i (\sigma^*q^{2+2i}; q)_{d-i}} \begin{bmatrix} d \\ i \end{bmatrix}_q {}_3\phi_2 \left( \begin{array}{c} q^{-i}, \sigma^*q^{i+1}, q^{-x} \\ \rho q, q^{-d} \end{array} \middle| q, q \right). \quad (A.3)
$$

Recurrence coefficients.

$$
b_i(\Phi) = \frac{\mu(1 - \sigma^* q^{i+1})(1 - q^{i-d})(1 - \rho q^{i+1})}{(1 - \sigma^* q^{2i+2})(1 - \sigma^* q^{2i+1})}, \quad (A.4a)
$$

$$
c_i(\Phi) = \frac{\mu q^{i-d} (1 - \sigma^* q^{d+i+1})(1 - q^i)(\sigma^* q^i - \rho)}{(1 - \sigma^* q^{2i+1})(1 - \sigma^* q^{2i})}. \quad (A.4b)
$$

A₂-contiguity recurrence relations.

$$
\frac{q(\rho - q\sigma^*)}{1 - q^2\sigma^*} p_{xi}(\Phi) = \frac{q^{i+1}(\rho - \sigma^* q^{i+1})}{1 - \sigma^* q^{2i+2}} p_{xi-1}(\Phi) + \frac{1 - \rho q^i}{1 - \sigma^* q^{2i}} p_{xi-1}(\Phi), \quad (A.5a)
$$

$$
(1 - q^{d+1-x}) \frac{1 - q\sigma^*}{\rho - \sigma^*} p_{xi}(\Phi) = \frac{q(1 - q^{i+1})(1 - \sigma^* q^{i+2+d})}{1 - \sigma^* q^{2i+2}} p_{xi+1}^{+}(\Phi) \\
\qquad + \frac{q^{i+1}(1 - q^{d+1-i})(1 - \sigma^* q^i)}{1 - \sigma^* q^{2i}} p_{xi}^{+}(\Phi), \quad (A.5b)
$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon, \sigma^{*} \to \sigma^{*}/q^{\epsilon 1}}$, for $\epsilon = \pm$.
### A.1.2 Dual q-Hahn

**Parameter array.**

$$ 
\begin{aligned}
\theta_i &= \theta_0 + \mu(1 - q^i)(1 - \sigma q^{i+1})q^{-i}, & \theta_i^* &= \theta_0^* + \mu^*(1 - q^i)q^{-i}, \\
\varphi_i &= \mu\mu^* q^{1-2i}(1 - q^i)(1 - q^{i-d-1})(1 - \rho q^i), & \phi_i &= \mu\mu^* q^{d+2-2i}(1 - q^i)(1 - q^{i-d-1})(\sigma - \rho q^{i-d-1}).
\end{aligned} 
\quad (A.6) $$

**Polynomial.**

$$ p_{xi}(\Phi) = (-1)^i q^{\frac{i}{2}(i-3-2d)} \frac{(q\rho; q)_i}{\sigma^i (q^{-d}\rho/\sigma; q)_i} \begin{bmatrix} d \\ i \end{bmatrix}_q 3\phi_2 \left( \begin{array}{c} q^{-i}, q^{-x}, \sigma q^{x+1} \\ \rho q, q^{-d} \end{array} \middle| q, q \right). \quad (A.7) $$

**Recurrence coefficients.**

$$ b_i(\Phi) = \mu(1 - q^{i-d})(1 - \rho q^{i+1}), \quad (A.8a) $$

$$ c_i(\Phi) = q\mu(1 - q^i)(\sigma - \rho q^{i-d-1}). \quad (A.8b) $$

**$A_2$-contiguity recurrence relations.**

$$ q(\rho - \sigma q^d)p_{xi}(\Phi) = q(\rho q^i - \sigma q^d)p_{xi}^{-}(\Phi) + (1 - \rho q^i)p_{xi-1}^{-}(\Phi), \quad (A.9a) $$

$$ \frac{(q^x - q^{d+1})(\sigma q^{d+1} - q^{-1-x})}{\rho - \sigma q^{d+1}} p_{xi}(\Phi) = (q^{i+1} - 1)p_{xi+1}^{+}(\Phi) + (q^{d+1} - q^i)p_{xi}^{+}(\Phi), \quad (A.9b) $$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

### A.1.3 Dual quantum q-Krawtchouk

**Parameter array.**

$$ 
\begin{aligned}
\theta_i &= \theta_0 + \mu(q^i - 1), & \theta_i^* &= \theta_0^* + \mu^*(q^{-i} - 1), \\
\varphi_i &= \mu\mu^*(1 - q^i)(1 - q^{d-i+1})/\rho, & \phi_i &= \mu\mu^*(1 - q^{-i})(1 - q^{d+1-i})(1 - q^i/\rho).
\end{aligned} 
\quad (A.10) $$

We slightly modify the parameters in comparison to [43]: $s \to \mu/q, r \to \mu q^d/\rho, h^* \to \mu^*$. We also change the name in comparison to [43] from quantum q-Krawtchouk to dual quantum q-Krawtchouk in order to match with the usual name of the polynomial used in [30] for example.

**Polynomial.**

$$ p_{xi}(\Phi) = (-1)^i q^{\frac{i(i-1)}{2}} \frac{1}{(\rho q^{-i}; q)_i} \begin{bmatrix} d \\ i \end{bmatrix}_q 2\phi_1 \left( \begin{array}{c} q^{-i}, q^{-x} \\ q^{-d} \end{array} \middle| q, \rho q^{x-d} \right). \quad (A.11) $$

**Recurrence coefficients.**

$$ b_i(\Phi) = \mu(1 - q^{d-i})q^{2i+1}/\rho, \quad (A.12a) $$

$$ c_i(\Phi) = \mu(1 - q^i)(1 - q^i/\rho). \quad (A.12b) $$

**$A_2$-contiguity recurrence relations.**

$$ (\rho - q)p_{xi}(\Phi) = (\rho - q^{i+1})p_{xi}^{-}(\Phi) + q^{i}p_{xi-1}^{-}(\Phi), \quad (A.13a) $$

$$ \frac{\rho q^{d+1}}{\rho - 1} (1 - q^{x-d-1}) p_{xi}(\Phi) = (q^{i+1} - 1) p_{xi+1}^{+}(\Phi) + (q^{d+1} - q^i) p_{xi}^{+}(\Phi), \quad (A.13b) $$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1, \rho \to \rho q^{\epsilon 1}}$, for $\epsilon = \pm$.
### A.1.4 Quantum q-Krawtchouk

We also add in the list the quantum q-Krawtchouk, which is not in the list of [43] (we recall that the parameter array called quantum q-Krawtchouk in [43] is called here dual quantum q-Krawtchouk). It is not in the usual list since it can be obtained by a transformation $q \to q^{-1}$, but it is useful here to also provide the formulas for the quantum q-Krawtchouk.

**Parameter array.**

$$ 
\begin{aligned}
\theta_i &= \theta_0 + \mu(q^{-i} - 1), & \theta_i^* &= \theta_0^* + \mu^*(q^i - 1), \\
\varphi_i &= \mu\mu^*(1 - q^i)(1 - q^{d+1-i})/\rho, & \phi_i &= \mu\mu^*(1 - q^i)(1 - q^{i-d-1})(1 - q^{d+1-i}/\rho).
\end{aligned} 
\quad (A.14) $$

This parameter array is obtained by the transformation 2.8 of the parameter array (A.10).

**Polynomial.**

$$ p_{xi}(\Phi) = (-1)^i q^{\frac{i(i+1)}{2}} \frac{1}{q^{id}(\rho q^{-d}; q)_i} \begin{bmatrix} d \\ i \end{bmatrix}_q 2\phi_1 \left( \begin{array}{c} q^{-x}, \quad q^{-i} \\ q^{-d} \end{array} \middle| \rho, \rho q^{i-d} \right). \quad (A.15) $$

**Recurrence coefficients.**

$$ b_i(\Phi) = -\mu(1 - q^{d-i})q^{-i}/\rho, \quad (A.16a) $$

$$ c_i(\Phi) = \mu(1 - q^{-i})(1 - q^{d+1-i}/\rho). \quad (A.16b) $$

**$A_2$-contiguity recurrence relations.**

$$ (\rho - q^d)p_{xi}(\Phi) = (\rho - q^{d-i})p_{xi-1}^{-}(\Phi) + q^{d+1-i}p_{xi-1}^{-}(\Phi), \quad (A.17a) $$

$$ \frac{\rho}{q^{d+1} - \rho} (1 - q^{d+1-x}) p_{xi}(\Phi) = (1 - q^{-i-1}) q^{d+1} p_{xi+1}^{+}(\Phi) + (q^{d+1-i} - 1) p_{xi}^{+}(\Phi), \quad (A.17b) $$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

### A.1.5 q-Krawtchouk

**Parameter array.**

$$ 
\begin{aligned}
\theta_i &= \theta_0 - \mu(1 - q^{-i}), & \theta_i^* &= \theta_0^* - \mu^*(1 - q^{-i})(1 - \sigma^* q^{i+1}), \\
\varphi_i &= \mu\mu^* q^{1-2i}(1 - q^i)(1 - q^{i-d-1}), & \phi_i &= \mu\mu^*\sigma^* q(1 - q^i)(1 - q^{i-d-1}).
\end{aligned} 
\quad (A.18) $$

**Polynomial.**

$$ p_{xi}(\Phi) = (-1)^i q^{-\frac{i}{2}(3+i)} \frac{(\sigma^* q^2; q)_d}{(\sigma^*)^i (\sigma^* q^{1+i}; q)_i (\sigma^* q^{2+2i}; q)_{d-i}} \begin{bmatrix} d \\ i \end{bmatrix}_q 3\phi_2 \left( \begin{array}{c} q^{-i}, \quad \sigma^* q^{i+1}, \quad q^{-x} \\ 0, \quad q^{-d} \end{array} \middle| q, q \right). \quad (A.19) $$

**Recurrence coefficients.**

$$ b_i(\Phi) = \frac{\mu(1 - q^{i-d})(1 - \sigma^* q^{i+1})}{(1 - \sigma^* q^{2i+1})(1 - \sigma^* q^{2i+2})}, \quad (A.20a) $$

$$ c_i(\Phi) = \frac{\mu\sigma^* q^{2i-d}(1 - q^i)(1 - \sigma^* q^{d+i+1})}{(1 - \sigma^* q^{2i})(1 - \sigma^* q^{2i+1})}. \quad (A.20b) $$
A$_2$-contiguity recurrence relations.

$$p_{xi}(\Phi) = \frac{q^{2i}(1 - q^2\sigma^*)}{1 - \sigma^*q^{2i+2}} p_{xi}^{-}(\Phi) - \frac{(1 - q^2\sigma^*)}{q^2\sigma^*(1 - \sigma^*q^{2i})} p_{xi-1}^{-}(\Phi), \qquad (A.21a)$$

$$(q^{d+1-x} - 1)p_{xi}(\Phi) = \frac{\sigma^*q(1 - q^{i+1})(1 - \sigma^*q^{i+d+2})}{(1 - \sigma^*q^{2i+2})(1 - q\sigma^*)} p_{xi+1}^{+}(\Phi) + \frac{\sigma^*q(q^i - q^{d+1})(1 - \sigma^*q^i)}{(1 - \sigma^*q^{2i})(1 - q\sigma^*)} p_{xi}^{+}(\Phi), \qquad (A.21b)$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1, \sigma^* \to \sigma^*/q^{\epsilon 1}}$, for $\epsilon = \pm$.

### A.1.6 Affine q-Krawtchouk

#### Parameter array.

$$\begin{aligned} \theta_i &= \theta_0 - \mu(1 - q^{-i}), & \theta_i^* &= \theta_0^* - \mu^*(1 - q^{-i}), \\ \varphi_i &= \mu\mu^*q^{1-2i}(1 - q^i)(1 - q^{i-d-1})(1 - \rho q^i), & \phi_i &= -\mu\mu^*\rho q^{1-i}(1 - q^i)(1 - q^{i-d-1}). \end{aligned} \qquad (A.22)$$

#### Polynomial.

$$p_{xi}(\Phi) = \frac{(\rho q; q)_i}{(q\rho)^i} \begin{bmatrix} d \\ i \end{bmatrix}_q 3\phi_2 \left( \begin{array}{c} q^{-i}, \ 0, \ q^{-x} \\ \rho q, \ q^{-d} \end{array} \middle| \begin{array}{c} q, q \end{array} \right). \qquad (A.23)$$

#### Recurrence coefficients.

$$b_i(\Phi) = \mu(1 - q^{i-d})(1 - \rho q^{i+1}), \qquad (A.24a)$$

$$c_i(\Phi) = -\rho\mu q^{i-d}(1 - q^i). \qquad (A.24b)$$

A$_2$-contiguity recurrence relations.

$$p_{xi}(\Phi) = q^i p_{xi}^{-}(\Phi) + q^{-1}(1/\rho - q^i) p_{xi-1}^{-}(\Phi), \qquad (A.25a)$$

$$(1 - q^{d+1-x}) p_{xi}(\Phi) = q\rho(1 - q^{i+1}) p_{xi+1}^{+}(\Phi) + q\rho(q^i - q^{d+1}) p_{xi}^{+}(\Phi), \qquad (A.25b)$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

### A.1.7 Dual q-Krawtchouk

#### Parameter array.

$$\begin{aligned} \theta_i &= \theta_0 + \mu(1 - q^i)(1 - \sigma q^{i+1})q^{-i}, & \theta_i^* &= \theta_0^* + \mu^*(1 - q^i)q^{-i}, \\ \varphi_i &= \mu\mu^*q^{1-2i}(1 - q^i)(1 - q^{i-d-1}), & \phi_i &= \mu\mu^*\sigma q^{d+2-2i}(1 - q^i)(1 - q^{i-d-1}). \end{aligned} \qquad (A.26)$$

#### Polynomial.

$$p_{xi}(\Phi) = (-1)^i \sigma^{-i} q^{\frac{i}{2}(i-2d-3)} \begin{bmatrix} d \\ i \end{bmatrix}_q 3\phi_2 \left( \begin{array}{c} q^{-i}, \ q^{-x}, \ \sigma q^{x+1} \\ 0, \ q^{-d} \end{array} \middle| \begin{array}{c} q, q \end{array} \right). \qquad (A.27)$$

#### Recurrence coefficients.

$$b_i(\Phi) = \mu(1 - q^{i-d}), \qquad (A.28a)$$

$$c_i(\Phi) = \mu\sigma q(1 - q^i). \qquad (A.28b)$$
A$_2$-contiguity recurrence relations.

$$
\sigma p_{xi}(\Phi) = \sigma p_{xi}^{-}(\Phi) - q^{-d-1}p_{xi-1}^{-}(\Phi), \qquad (A.29a)
$$

$$
(q^{-x} - q^{-d-1})(1 - \sigma q^{d+x+2})p_{xi}(\Phi) = \sigma q(1 - q^{i+1})p_{xi+1}^{+}(\Phi) + \sigma q^{i+1}(1 - q^{d-i+1})p_{xi}^{+}(\Phi), \quad (A.29b)
$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

A.2 List of parameter arrays of type II

A.2.1 Hahn

Parameter array.

$$
\begin{align*}
\theta_i &= \theta_0 + \sigma i, & \theta_i^* &= \theta_0^* + \mu^* i (i + 1 + \sigma^*), \\
\varphi_i &= \mu^* \sigma i (i - d - 1) (i + \rho), & \phi_i &= -\mu^* \sigma i (i - d - 1) (i + \sigma^* - \rho). \tag{A.30}
\end{align*}
$$

Polynomial.

$$
p_{xi}(\Phi) = \binom{d}{i} \frac{(\rho+1)_i (2+\sigma^*)_d}{(\sigma^* - \rho + 1, \sigma^* + 1 + i)_i (\sigma^* + 2i + 2)_{d-i}} \ {}_3F_2 \left( \begin{matrix} -i, i+1+\sigma^*, -x \\ \rho+1, -d \end{matrix} \middle| 1 \right). \quad (A.31)
$$

Recurrence coefficients.

$$
b_i(\Phi) = \frac{\sigma(i-d)(i+1+\rho)(i+\sigma^*+1)}{(2i+\sigma^*+1)(2i+\sigma^*+2)}, \qquad (A.32a)
$$

$$
c_i(\Phi) = - \frac{\sigma i (i + d + \sigma^* + 1)(i + \sigma^* - \rho)}{(2i + \sigma^*)(2i + \sigma^* + 1)}. \quad (A.32b)
$$

$A_2$-contiguity recurrence relations.

$$
\frac{\sigma^* + 1 - \rho}{2 + \sigma^*} p_{xi}(\Phi) = \frac{i + \rho}{2i + \sigma^*} p_{x,i-1}^{-}(\Phi) + \frac{i + 1 + \sigma^* - \rho}{2i + \sigma^* + 2} p_{x,i}^{-}(\Phi), \quad (A.33a)
$$

$$
\frac{\sigma^* + 1}{\rho - \sigma^*} (x - d - 1) p_{xi}(\Phi) = \frac{(d + 1 - i)(i + \sigma^*)}{2i + \sigma^*} p_{x,i}^{+}(\Phi) + \frac{(i + 1)(i + \sigma^* + 2 + d)}{2i + \sigma^* + 2} p_{x,i+1}^{+}(\Phi), \quad (A.33b)
$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1, \sigma^* \to \sigma^* - \epsilon 1}$, for $\epsilon = \pm$.

A.2.2 Dual Hahn

Parameter array.

$$
\begin{align*}
\theta_i &= \theta_0 + \mu i (i + 1 + \sigma), & \theta_i^* &= \theta_0^* + \sigma^* i, \\
\varphi_i &= \mu \sigma^* i (i - d - 1) (i + \rho), & \phi_i &= \mu \sigma^* i (i - d - 1) (i + \rho - \sigma - d - 1). \tag{A.34}
\end{align*}
$$

Polynomial.

$$
p_{xi}(\Phi) = (-1)^i \binom{d}{i} \frac{(\rho+1)_i}{(\rho-\sigma-d)_i} {}_3F_2 \left( \begin{matrix} -i, -x, x+1+\sigma \\ \rho+1, -d \end{matrix} \middle| 1 \right). \quad (A.35)
$$

Recurrence coefficients.

$$
b_i(\Phi) = \mu(i-d)(i+1+\rho), \qquad (A.36a)
$$

$$
c_i(\Phi) = \mu i (i + \rho - \sigma - d - 1). \qquad (A.36b)
$$
A$_2$-contiguity recurrence relations.

$$
(d + \sigma - \rho)p_{xi}(\Phi) = -(i + \rho - \sigma - d)p_{\bar{xi}}(\Phi) + (i + \rho)p_{\bar{xi}-1}(\Phi), \quad (A.37a)
$$

$$
\frac{(x - d - 1)(x + d + \sigma + 2)}{d - \rho + \sigma + 1} p_{xi}(\Phi) = (i - d - 1)p_{xi}^{+}(\Phi) - (i + 1)p_{xi+1}^{+}(\Phi), \quad (A.37b)
$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

A.2.3 Krawtchouk

Parameter array.

$$
\begin{align}
\theta_i &= \theta_0 + \sigma i, & \theta_i^* &= \theta_0^* + \sigma^* i, \nonumber \\
\varphi_i &= \rho i (i - d - 1), & \phi_i &= (\rho - \sigma \sigma^*) i (i - d - 1). \tag{A.38}
\end{align}
$$

Polynomial.

$$
p_{xi}(\Phi) = \left( \frac{\rho}{\sigma\sigma^* - \rho} \right)^i \binom{d}{i} {}_2F_1 \left( \begin{array}{c} -i, -x \\ -d \end{array} \middle| \rho^{-1}\sigma\sigma^* \right). \quad (A.39)
$$

Recurrence coefficients.

$$
b_i(\Phi) = \frac{\rho}{\sigma^*} (i-d), \qquad (A.40a)
$$

$$
c_i(\Phi) = \frac{\rho - \sigma\sigma^*}{\sigma^*} i. \qquad (A.40b)
$$

A$_2$-contiguity recurrence relations.

$$
p_{xi}(\Phi) = p_{\bar{xi}}(\Phi) + \frac{\rho}{\sigma\sigma^* - \rho} p_{x,i-1}(\Phi), \qquad (A.41a)
$$

$$
\frac{\sigma\sigma^*}{\rho - \sigma\sigma^*} (x - d - 1) p_{xi}(\Phi) = (d + 1 - i) p_{xi}^{+}(\Phi) + (i + 1) p_{x,i+1}^{+}(\Phi), \quad (A.41b)
$$

where $p_{xi}^{\epsilon}(\Phi) = p_{xi}(\Phi)|_{d \to d+\epsilon 1}$, for $\epsilon = \pm$.

References

[1] E. Bannai, and T. Ito, *Algebraic Combinatorics I: Association Schemes*, Benjamin/Cummings, London, 1984.

[2] E. Bannai, H. Kurihara, D. Zhao, and Y. Zhu, *Multivariate P- and/or Q-polynomial association schemes*, (2023) arXiv:2305.00707.

[3] P.-A. Bernard, N. Crampe, L. Poulain d'Andecy, L. Vinet, and M. Zaimi, *Bivariate P-polynomial association schemes*, Algebraic combinatorics (2024) and arXiv:2212.10824.

[4] P.-A. Bernard, N. Crampe, L. Vinet, M. Zaimi, and X. Zhang, *m-distance-regular graphs and their relation to multivariate P-polynomial association schemes*, (2023) arXiv:2309.16016.

[5] J. Cooke, and A. Lacabanne, *Higher Rank Askey–Wilson Algebras as Skein Algebras*, (2022) arXiv:2205.04414.

[6] N. Crampe, L. Frappat, J. Gaboriaud, L. Poulain d'Andecy, E. Ragoucy, and L. Vinet, *The Askey–Wilson algebra and its avatars*, J. Phys. A 54 (2021) 063001 and arXiv:2009.14815.
[7] N. Crampe, L. Frappat, J. Gaboriaud, E. Ragoucy, L. Vinet, and M. Zaimi, λ-Griffiths polynomials: Bispectrality and biorthogonality, (2023) arXiv:2311.03256.

[8] N. Crampe, L. Frappat, L. Poulain d'Andecy, and E. Ragoucy, The higher-rank Askey-Wilson algebra and its braid group automorphisms, SIGMA 19 (2023) 077 and arXiv:2303.17677.

[9] N. Crampe, L. Frappat, and E. Ragoucy, Representations of the rank two Racah algebra and orthogonal multivariate polynomials, Linear Algebra Its Appl. 664 (2023) 165-215 and arXiv:2206.01031.

[10] N. Crampe, L. Poulain d'Andecy, L. Vinet, and M. Zaimi, Askey-Wilson braid algebra and centralizer of $U_q(s_l_2)$, Ann. Henri Poincaré 24 (2023) 1897-1922 and arXiv:2206.11150.

[11] N. Crampe, W. van de Vijver, and Luc Vinet, Racah problems for the oscillator algebra, the Lie algebra $sl_n$, and multivariate Krawtchouk polynomials, Ann. Henri Poincaré 21 (2020) 3939-3971 and arXiv:1909.12643

[12] N. Crampe, L. Vinet, and M. Zaimi, Temperley-Lieb, Birman-Murakami-Wenzl and Askey-Wilson algebras and other centralizers of $U_q(s_l_2)$, Ann. Henri Poincaré 22 (2021) 3499-3528, arXiv:2008.04905.

[13] N. Crampe, L. Vinet, M. Zaimi, and X. Zhang, A bivariate $Q$-polynomial structure for the non-binary Johnson scheme, J. Comb. Theory A 202 (2024) 105829 and arXiv:2306.01882.

[14] H. De Bie, H. De Clercq, and W. van de Vijver, The higher rank $q$-deformed Bannai-Ito and Askey-Wilson algebra, Commun. Math. Phys. 374 (2020) 277, doi:10.1007/s00220-019-03562-w and arXiv:1805.06642.

[15] H. De Bie, V.X. Genest, W. van de Vijver, and L. Vinet, A higher rank Racah algebra and the $Z_n^2$ Laplace-Dunkl operator, J. Phys. A: Math. Theor. 51 (2018) 025203, doi:10.1088/1751-8121/aa9756 and arXiv:1610.02638.

[16] H. De Clercq, Higher Rank Relations for the Askey-Wilson and $q$-Bannai-Ito Algebra, SIGMA 15 (2019) 099, doi:10.3842/SIGMA.2019.099 and arXiv:1908.11654.

[17] G. Gasper, and M. Rahman, Basic Hypergeometric Series, Second Edition (Corrected and Extended). Cambridge University Press, 2004.

[18] V.X. Genest, S. Post, and L. Vinet, An algebraic interpretation of the multivariate $q$-Krawtchouk polynomials, Ramanujan J. 43 (2017) 415-445 and arXiv:1508.07770.

[19] V.X. Genest, L. Vinet, and A. Zhedanov, The multivariate Krawtchouk polynomials as matrix elements of the rotation group representations on oscillator states, J. Phys. A: Math. Theor. 46 (2013) 505203 and arXiv:1306.4256.

[20] V.X. Genest, L. Vinet, and A. Zhedanov, Superintegrability in two dimensions and the Racah-Wilson algebra, Lett. Math. Phys. 104 (2014) 931-952 and arXiv:1307.5539.

[21] J. Geronimo, and P. Iliev, Bispectrality of multivariable Racah-Wilson polynomials, Constr. Approx. 31 (2010) 417-457 and arXiv:0705.1469.

[22] J.S. Geronimo, and P. Iliev, Multivariable Askey-Wilson function and bispectrality, Ramanujan J. 24 (2011) 273-287, doi:10.1007/s11139-010-9244-3.

[23] Y.I. Granovskii, and A.S. Zhedanov, Hidden symmetry of the Racah and Clebsch-Gordan problems for the quantum algebra $sl_q(2)$, J. Group Theoret. Methods Phys. 1 (1993) 161-171 and arXiv:hep-th/9304138.
[24] R. Griffiths, *Orthogonal polynomials on the multinomial distribution*, Austral. J. Statist. 13 (1971) 27-35; Corrigenda Austral. J. Statist. 14 (1972) 270.

[25] F.A. Grünbaum, *The Rahman polynomials are bispectral*, SIGMA 3 (2007) 65-75 and arXiv:0705.0468.

[26] M.R. Hoare, and M. Rahman, *A probabilistic origin for a new class of bivariate polynomials*, SIGMA 4 (2008) 89-106.

[27] P. Iliev, and P. Terwilliger, *The Rahman polynomials and the Lie algebra $sl_3(C)$*, Trans. Amer. Math. Soc. 364 (2012) 4225-4238 and arXiv:1006.5062.

[28] P. Iliev, and Y. Xu, *Hahn polynomials for hypergeometric distribution*, Adv. in Appl. Math. 139 (2022) 102364 and arXiv:2012.12168.

[29] E.G. Kalnins, J.M. Kress, and Jr.W. Miller, *Second-order superintegrable systems in conformally flat spaces. I. Two-dimensional classical structure theory*, J. Math. Phys. 46 (2005) 053509.

[30] R. Koekoek, P.A. Lesky, and R.F. Swarttouw, *Hypergeometric orthogonal polynomials and their q-analogues*, Springer Monographs in Mathematics, Springer Berlin, Heidelberg, 2010.

[31] D. Leonard, *Orthogonal polynomials, duality, and association schemes*, SIAM J. Math. Anal. 13 (1982) 656-663.

[32] H. Mizukawa, and H. Tanaka, $(n+1, m+1)$-hypergeometric functions associated to character algebras, *Proceedings of the American Mathematical Society* 132 (2004) 2613-2618.

[33] S. Post, *Models of quadratic algebras generated by superintegrable systems in 2D*, SIGMA 7 (2011) 036 and arXiv:1104.0734.

[34] S. Post, and A. Walter, *A higher rank extension of the Askey-Wilson Algebra*, arXiv:1705.01860.

[35] A.V. Rosenblyum, *Spectral analysis of generators of representations of the group $U(3)$*, Theor. Math. Phys. 73 (1987) 1352-1356.

[36] H. Rosengren, *Multivariable orthogonal polynomials and coupling coefficients for discrete series representations*, SIAM J. Math. Anal. 30 (1998) 233-272.

[37] F. Scarabotti, *Multidimensional Hahn polynomials, intertwining functions on the symmetric group and Clebsch-Gordon coefficients*, Methods Appl. Anal. 14 (2007) 355-386 and arXiv:0805.0670.

[38] H. Tarnanen, M.J. Aaltonen, and J.-M. Goethals, *On the Nonbinary Johnson Scheme*, Eur. J. Comb. 6 (1985) 279-285.

[39] P. Terwilliger, *The Subconstituent Algebra of an Association Scheme, (PartI)*, J. Algebr. Comb. 1 (1992) 363-388.

[40] P. Terwilliger, *The Subconstituent Algebra of an Association Scheme, (PartII)*, J. Algebr. Comb. 2 (1993) 73-103.

[41] P. Terwilliger, *The Subconstituent Algebra of an Association Scheme, (PartIII)*, J. Algebr. Comb. 2 (1993) 177-210.

[42] P. Terwilliger, *Two linear transformations each tridiagonal with respect to an eigenbasis of the other*, Linear Algebra Its Appl. 330 (2001) 149-203, doi:10.1007/s10623-004-4862-7 and arXiv:math/0406555.
[43] P. Terwilliger, *Two linear transformations each tridiagonal with respect to an eigenbasis of the other; an algebraic approach to the Askey scheme of orthogonal polynomials*, (2004) [arXiv:math/0408390].

[44] P. Terwilliger, *Notes on the Leonard system classification*, (2004) [arXiv:math/0305356].

[45] P. Terwilliger, and R. Vidunas, *Leonard pairs and the Askey–Wilson relations*, *J. Algebra Its Appl.* **03** (2004) 411–426 and [arXiv:math/0305356].

[46] M.V. Tratnik, *Some multivariable orthogonal polynomials of the Askey tableau-discrete families*, *J. Math. Phys.* **32** (1991) 2337–2342.

[47] A.S. Zhedanov, “*Hidden symmetry*” of the Askey–Wilson polynomials, *Theor. Math. Phys.* **89** (1991) 1146–1157, doi:10.1007/BF01015906.

[48] A. Zhedanov, *9j-symbols of the oscillator algebra and Krawtchouk polynomials in two variables*, *J. Phys. A* **30** (1997) 8337–8357.