# The *q*-Racah polynomials from scalar products of Bethe states II

Pascal Baseilhac, Rodrigo A. Pimenta

► To cite this version:

Pascal Baseilhac, Rodrigo A. Pimenta. The *q*-Racah polynomials from scalar products of Bethe states II. Journal of Physics A: Mathematical and Theoretical, 2025, 58 (12), pp.125205. 10.1088/1751-8121/adbfe7. hal-05394690

HAL Id: hal-05394690

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

Submitted on 3 Dec 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
THE *q*-RACAH POLYNOMIALS FROM SCALAR PRODUCTS OF BETHE STATES II

PASCAL BASEILHAC* AND RODRIGO A. PIMENTA**

**ABSTRACT.** The theory of Leonard triples is applied to the derivation of normalized scalar products of on-shell and off-shell Bethe states generated from a Leonard pair. The scalar products take the form of linear combinations of *q*-Racah polynomials with coefficients depending on the off-shell parameters. Upon specializations, explicit solutions for the corresponding Belliard-Slavnov linear systems are obtained. It implies the existence of a determinant formula in terms of inhomogeneous Bethe roots for the *q*-Racah polynomials. Also, a set of relations that determines solutions (Bethe roots) of the corresponding Bethe equations of inhomogeneous type in terms of solutions of Bethe equations of homogeneous type is obtained.

MSC: 33D45; 81R50; 81U15.

**Keywords:** Askey-Wilson algebra; Leonard pairs; Orthogonal polynomials; Bethe ansatz

# 1. INTRODUCTION

In a recent series of papers [BaP19, BerCSV20, BerCCV22, BP22], a correspondence between the representation theory of the Askey-Wilson algebra [T87, Zh91] or Racah algebra [LL65, GZ88] and the framework of the algebraic Bethe ansatz [STF79, Sk88] has been investigated. More precisely, this correspondence aims to relate two subjects that have been independently developed in the context of representation theory and mathematical physics, respectively. On one hand, irreducible finite dimensional representations of the Askey-Wilson or Racah algebras have been extensively studied in recent years. They are classified using the theory of Leonard pairs [TV03, T04]. In this context, it is shown that the entries of the transition matrices relating different eigenbases of elements of a Leonard pair correspond to discrete orthogonal polynomials that belong to the Askey-scheme [T03]. On the other hand, the framework of the algebraic Bethe ansatz has a long history (for a review and references therein, see e.g. [Sl22]) and provides an efficient approach for deriving exact results in a large class of quantum integrable systems. In order to solve quantum integrable models with $U(1)$ symmetry breaking boundary fields, the original approach is upgraded to the *modified* algebraic Bethe ansatz approach developed in [BelC13, Bel15, C15, ABGP15, BelP15]. In both approaches, the spectral problem for the Hamiltonian of the system (eigenvalues and eigenvectors) is formulated in terms of Bethe roots satisfying a system of transcendental equations known as the Bethe equations whereas the eigenvectors are Bethe states. Depending on the system under consideration, the Bethe equations are either of homogeneous type or of inhomogeneous type [WYCS15]. In particular, some simple examples of integrable systems<fnref n="1" /> and corresponding Bethe equations can be studied using the theory of Leonard pair [BaP19]. Thus, establishing a Leonard pair/algebraic Bethe ansatz correspondence aims to get new insights in each subject. Also, having in mind that the theory of Leonard pairs generalizes to the theory of tridiagonal pairs [T99], a Leonard pair/algebraic Bethe ansatz correspondence can be viewed as a toy model for a tridiagonal pair/algebraic Bethe ansatz correspondence. The latter should find applications in the theory of special functions generalizing the Askey-scheme and quantum integrable systems with non-trivial boundary conditions and finite dimensional quantum space. More details are given in the last section.

The starting point of the Leonard pair/algebraic Bethe ansatz correspondence is the observation [BaP19] that the universal exchange relations which typically arise in the algebraic Bethe ansatz analysis of various examples of quantum integrable spin chains with boundaries are satisfied by so-called dynamical operators generated from the Askey-Wilson algebra. A similar observation holds for the Racah algebra [BerCCV22]. This opens the possibility to reconsider the main objects of the theory of Leonard pair (eigenbases, transition coefficients, discrete orthogonal polynomials of the Askey-scheme,...) from the perspective of the algebraic Bethe ansatz and its modified version (on-shell Bethe states and related Bethe equations, scalar products of Bethe states, Baxter Q-polynomial,...) and reciprocally.

<footnote n="1">¹The *q*-analog of the quantum Euler top and the '1-site' open XXZ spin chain [BaP19, Subsection 5.1]; A 3-sites Heisenberg chain with inhomogeneous couplings [BaP19, Subsection 5.2].</footnote>
For instance, in [BaP19] the spectral problem for certain combinations of elements of a Leonard pair of q-Racah type - known as the Heun-Askey-Wilson operator - has been solved using the theory of Leonard pairs combined with the framework of the modified algebraic Bethe ansatz [BaP19]. For the limiting case $q = 1$ associated with the Racah algebra and related Leonard pairs, the analog analysis is done in [BerCSV20, BerCCV22]. In [BP22], different types of eigenbases and dual eigenbases of Leonard pairs of q-Racah type are constructed using the framework of the modified algebraic Bethe ansatz: it is shown that eigenbases for a Leonard pair can be built from Bethe states associated with either homogeneous or inhomogeneous Bethe equations. This lead to an interpretation of the q-Racah polynomials<fnref n="2" /> as certain ratios of scalar products of various types of on-shell Bethe states.

The present paper is the second part of [BP22]. Here, we investigate further the Leonard pair/algebraic Bethe ansatz correspondence motivated by the following problems.

(a) In [BP22], it was shown that certain ratios of scalar products of on-shell versus on-shell Bethe states are identical to q-Racah polynomials. A first goal of this paper is to extend the analysis to the computation of scalar products of on-shell versus off-shell Bethe states. The strategy is to use the theory of Leonard pairs [T01, T03, T04, TV03] and Leonard triples [Curt07] (see also [H11]).

(b) For some class of scalar products of on-shell versus off-shell Bethe states, compact determinant formulas involving Bethe roots can be derived by solving certain systems of linear equations constructed by Belliard and Slavnov [BelS19]. A second goal of this paper is to identify the Belliard-Slavnov system of linear equations satisfied by some specialization of the scalar products computed in (a). This aims to derive a new presentation of q-Racah polynomials using determinants of matrices associated with inhomogenous Bethe roots.

(c) In the modified algebraic Bethe ansatz framework, inhomogeneous Bethe equations arise. Contrary to the standard algebraic Bethe ansatz where Bethe roots are classified using an integer $N$, there is no known criteria to classify corresponding Bethe roots. The third goal of this paper is to identify some criteria for the class of inhomogenous Bethe equations under consideration.

Note that (a), (b) and (c) can be viewed as a toy model for studying scalar products of Bethe states in the open XXZ spin chain from the perspective of the theory of tridiagonal pairs, see more details in the Concluding Remarks section.

The text is organized as follows. In Section 2, some known results on the Leonard pair of q-Racah type/algebraic Bethe ansatz correspondence are reviewed based on [BP22]. The exchange relations associated with a Leonard pair $(A, A^*)$ and corresponding dynamical operators that generate off-shell Bethe states are briefly recalled. Then, eigenbases for $(A, A^*)$ built from on-shell Bethe states associated with either homogeneous or inhomogeneous Bethe equations are given. In Section 3, the goal is to compute explicitly normalized scalar products of on-shell Bethe states versus off-shell Bethe states using the Leonard pairs data, in order to extend the results of [BP22]. To this end, the Leonard triple of q-Racah type $(A, A^*, A^\circ)$ is introduced. Eigenbases for the triple $(A, A^*, A^\circ)$ and related transition matrices are given explicitly. Using the triple, the normalized scalar products are computed, see Theorem 1. Upon specializations of previous results, a relation characterizing the Bethe roots of inhomogenous type associated with the eigenvalue $\theta_N^*, N = 0, \dots, 2s$ in terms of Bethe roots of homogeneous type associated with the same eigenvalue $\theta_N^*$ is obtained, see Proposition 3.4. In Section 4, the Belliard-Slavnov system of linear equations satisfied by certain specializations of the normalized scalar products of Theorem 1 are derived, see Proposition 4.1. As shown in Proposition 4.2, a determinant formula exist for such normalized scalar products. Upon specialization to on-shell Bethe states, it implies the existence of a determinant formula for the q-Racah polynomials, see Corollary 4.1. A summary of the results and some perspectives are given in the last section. Some formulas and explicit examples are reported in Appendix A and B, respectively.

**Notations:** The parameter $q$ is assumed not to be a root of unity. We write $[X, Y]_q = qXY - q^{-1}YX$ and $[n]_q = \frac{q^n - q^{-n}}{q - q^{-1}}$. We use the standard $q$-shifted factorials:

$$ (1.1) \qquad (a; q)_n = \prod_{k=0}^{n-1} (1 - a q^k), \qquad (a_1, a_2, \dots, a_k; q)_n = \prod_{j=1}^{k} (a_j; q), $$

<footnote n="2">²All other families of orthogonal polynomials of the scheme can be reached from the q-Racah polynomials by various limit transitions (either in the scalar parameters entering into the definition of the polynomials, or $q \to 1$). See [KS96, K10] for details.</footnote>
and define

$$ (1.2) \qquad b(x) = x - x^{-1}. $$

## 2. ALGEBRAIC BETHE ANSATZ AND LEONARD PAIRS OF q-RACAH TYPE

In this section, the correspondence between the theory of Leonard pairs of q-Racah type and the algebraic Bethe ansatz formalism is reviewed (see [BaP19, BP22] and references therein). Off-shell Bethe states are built using the image of dynamical operators satisfying the exchange relations associated with the Askey-Wilson algebra, and eigenbases for each operator of the Leonard pair ($A, A^*$) are given in the form of on-shell Bethe states associated either with Bethe equations of homogeneous type or Bethe equations of inhomogeneous type.

### 2.1. Dynamical operators and the Askey-Wilson algebra.
In the context of the quantum inverse scattering method, the entries {$\mathcal{A}(u), \mathcal{B}(u), \mathcal{C}(u), \mathcal{D}(u)$} with $u \in \mathbb{C}^*$ of the so-called (fundamental spin-$\frac{1}{2}$) monodromy matrix or double-row monodromy matrix are the key ingredients in the analysis of transfer matrix' spectral problem, scalar products of Bethe states, correlation functions and form factors. For a review, see [Sl22]. Importantly, according to the integrable model under consideration those entries satisfy a set of exchange relations that are equivalent to either a Yang-Baxter equation or a reflection equation.

In order to implement a solution based on the algebraic Bethe ansatz approach to some of the fundamental examples of quantum integrable models (see e.g. [TF79, CLSW03]) the analysis however requires to introduce certain linear combinations of those entries, which are called dynamical operators. Let us introduce the fundamental dynamical operators

$$ (2.1) \qquad \{\mathcal{A}^{\epsilon}(u, m), \mathcal{B}^{\epsilon}(u, m), \mathcal{C}^{\epsilon}(u, m), \mathcal{D}^{\epsilon}(u, m) | m \in \mathbb{N}^{*}, \epsilon = \pm 1\} $$

subject to the set of universal exchange relations<fnref n="3" />:

$$ (2.2) \qquad \mathcal{B}^{\epsilon}(u, m + 2)\mathcal{B}^{\epsilon}(v, m) = \mathcal{B}^{\epsilon}(v, m + 2)\mathcal{B}^{\epsilon}(u, m), $$

$$ (2.3) \qquad \mathcal{A}^{\epsilon}(u, m + 2)\mathcal{B}^{\epsilon}(v, m) = f(u, v)\mathcal{B}^{\epsilon}(v, m)\mathcal{A}^{\epsilon}(u, m) \\ \qquad \qquad +g(u, v, m)\mathcal{B}^{\epsilon}(u, m)\mathcal{A}^{\epsilon}(v, m) + w(u, v, m)\mathcal{B}^{\epsilon}(u, m)\mathcal{D}^{\epsilon}(v, m), $$

$$ (2.4) \qquad \mathcal{D}^{\epsilon}(u, m + 2)\mathcal{B}^{\epsilon}(v, m) = h(u, v)\mathcal{B}^{\epsilon}(v, m)\mathcal{D}^{\epsilon}(u, m), \\ \qquad \qquad +k(u, v, m)\mathcal{B}^{\epsilon}(u, m)\mathcal{D}^{\epsilon}(v, m) + n(u, v, m)\mathcal{B}^{\epsilon}(u, m)\mathcal{A}^{\epsilon}(v, m), $$

$$ (2.5) \qquad \mathcal{C}^{\epsilon}(u, m + 2)\mathcal{B}^{\epsilon}(v, m) = \mathcal{B}^{\epsilon}(v, m - 2)\mathcal{C}^{\epsilon}(u, m) \\ \qquad \qquad +q(u, v, m)\mathcal{A}^{\epsilon}(v, m)\mathcal{D}^{\epsilon}(u, m) + r(u, v, m)\mathcal{A}^{\epsilon}(u, m)\mathcal{D}^{\epsilon}(v, m) \\ \qquad \qquad +s(u, v, m)\mathcal{A}^{\epsilon}(u, m)\mathcal{A}^{\epsilon}(v, m) + x(u, v, m)\mathcal{A}^{\epsilon}(v, m)\mathcal{A}^{\epsilon}(u, m) \\ \qquad \qquad +y(u, v, m)\mathcal{D}^{\epsilon}(u, m)\mathcal{A}^{\epsilon}(v, m) + z(u, v, m)\mathcal{D}^{\epsilon}(u, m)\mathcal{D}^{\epsilon}(v, m) $$

and

$$ (2.6) \qquad \mathcal{C}^{\epsilon}(u, m - 2)\mathcal{C}^{\epsilon}(v, m) = \mathcal{C}^{\epsilon}(v, m - 2)\mathcal{C}^{\epsilon}(u, m), $$

$$ (2.7) \qquad \mathcal{C}^{\epsilon}(v, m + 2)\mathcal{A}^{\epsilon}(u, m + 2) = f(u, v)\mathcal{A}^{\epsilon}(u, m)\mathcal{C}^{\epsilon}(v, m + 2) \\ \qquad \qquad +g(u, v, m)\mathcal{A}^{\epsilon}(v, m)\mathcal{C}^{\epsilon}(u, m + 2) + w(u, v, m)\mathcal{D}^{\epsilon}(v, m)\mathcal{C}^{\epsilon}(u, m + 2), $$

$$ (2.8) \qquad \mathcal{C}^{\epsilon}(v, m + 2)\mathcal{D}^{\epsilon}(u, m + 2) = h(u, v)\mathcal{D}^{\epsilon}(u, m)\mathcal{C}^{\epsilon}(v, m + 2) \\ \qquad \qquad +k(u, v, m)\mathcal{D}^{\epsilon}(v, m)\mathcal{C}^{\epsilon}(u, m + 2) + n(u, v, m)\mathcal{A}^{\epsilon}(v, m)\mathcal{C}^{\epsilon}(u, m + 2), $$

where the explicit expressions for the coefficients $f(u, v), g(u, v, m), \dots$ are reported in Appendix A.

The above set of exchange relations is universal in the following sense. It follows from a reflection equation - up to conventions, see [CLSW03, eq. (12)] - which provides a presentation for the alternating central extension of the q-Onsager algebra $\mathcal{A}_q$ [BS09, T21]. By construction, the corresponding dynamical operators are power series of the inverse spectral parameter $u^{-1}$ (or $u$ depending on the convention chosen) with coefficients in $\mathcal{A}_q$. For several examples of quantum integrable models such as the open XXZ spin-$j$ chain [D07], the trigonometric solid-on-solid model [FK10], the alternating spin chain [CYSW14], the relativistic Toda chain [ZCYSW17], the lattice sine-Gordon model [MNP17] or the three-sites interacting model constructed in

<footnote n="3">³Those relations have been first introduced in the analysis of the open XXZ spin-1/2 chain with generic integrable boundary conditions [CLSW03]. In this case, the dynamical operators specialize to operators acting on finite dimensional tensor product representations of $U_q(\overline{sl}_2)$.</footnote>
[BaP19, Section 5.2], the fundamental exchange relations take the form (2.2)-(2.8) where the dynamical opera-
tors {$\mathcal{A}^\epsilon(u, m)$, $\mathcal{B}^\epsilon(u, m)$, $\mathcal{C}^\epsilon(u, m)$, $\mathcal{D}^\epsilon(u, m)$} specialize to certain Laurent polynomials in $u$ with coefficients in
a certain quotient of $\mathcal{A}_q$.

In continuity with the previous works [BaP19, BP22], in this paper we consider the quotient of $\mathcal{A}_q$ known
as the Askey-Wilson algebra. The Askey-Wilson algebra (AW) is generated by $A, A^*$ subject to the relations
[T87, Zh91]

$$
(2.9) \qquad [A, [A, A^*]]_{q-1} = \rho A^* + \omega A + \eta \mathcal{I},
$$

$$
(2.10) \qquad [A^*, [A^*, A]]_q|_{q-1} = \rho A + \omega A^* + \eta^* \mathcal{I},
$$

where the structure constants are {$\rho, \omega, \eta, \eta^*$} $\in \mathbb{C}^*$ and the identity element is denoted $\mathcal{I}$. In this case, the
explicit expressions of the dynamical operators that solve the exchange relations (2.2)-(2.8) are given in [BaP19,
Appendix 4]. As Laurent polynomials in $u$ with coefficients in AW, the dynamical operators depend on the free
scalars $\alpha, \beta$ (see Appendix A) and $\chi \in \mathbb{C}^*$ (called ‘gauge parameters’). For instance, one has:

$$
\begin{equation} \tag{2.11}
\begin{split}
\mathcal{B}^+(u,m) = \frac{\beta b(u^2)}{\alpha q^{2m+2} - \beta} & \left( \frac{\chi q^m}{\beta \rho u} [A^*, A]_q - \frac{\beta q^{-m}}{u \chi} [A, A^*]_q + \frac{(q^2+1)}{qu} A - \left(\frac{1}{qu^3} + qu\right) A^* \right. \\
& \left. \qquad + \frac{q^{-m-3}}{\beta \rho \chi b(q^2)} \left( (q^2+1)u^{-1} \left(\omega(\chi^2 q^{2m+2} - \beta^2 q^2 \rho) + \beta(q^4-1)\eta^* \chi q^m \right) - q^2 \rho (q^2 u + u^{-3}) (\beta^2 \rho - \chi^2 q^{2m}) \right) \right),
\end{split}
\end{equation}
\begin{equation*}
\begin{split}
\mathcal{B}^-(u,m) = \frac{\beta b(u^2) q^{2m+1}}{\alpha q^{-2} - \beta q^{2m}} & \left( \frac{u \chi q^{-m-1}}{\beta \rho} [A^*, A]_q - \frac{\beta u q^{m-1}}{\chi} [A, A^*]_q + \frac{(q^2 u^4 + 1)}{q^2 u} A - \left(\frac{q^2+1}{q^2} u \right) A^* \right. \\
& \left. \qquad - \frac{q^{-m-4}}{\beta \rho \chi b(q^2)} \left( \rho q^2 (q^2 u^3 + u^{-1}) (\beta^2 \rho q^{2m} - \chi^2) + (q^2+1)u (\omega(\beta^2 \rho q^{2m+2} - q^2 \chi^2) + \beta(q^4-1)\eta \chi q^m) \right) \right).
\end{split}
\end{equation*}
$$

2.2. Leonard pairs and off-shell/on-shell Bethe states. For most of the known examples of quantum
integrable models associated with the exchange relations (2.2)-(2.8), the dynamical operators act on a finite
dimensional vector space. Following [BaP19, BP22], here we consider irreducible finite dimensional representa-
tions of AW which relate to the concept of Leonard pair of q-Racah type. Eigenbases for a Leonard pair (A, A*)
and their duals are first recalled. Then, following [BP22] off-shell Bethe states associated with the Leonard pair
and the correspondence between on-shell Bethe states and eigenvectors of the Leonard pair is given.

2.2.1. *The Leonard pair (A, A*) and related eigenbases.* Let $\mathcal{V}$ be a vector space of positive finite dimension $\dim(\mathcal{V}) = 2s + 1$, where $s$ is an integer or half-integer. Define $\pi : \text{AW} \rightarrow \text{End}(\mathcal{V})$. Assume $\pi(\text{A}), \pi(\text{A}^*)$ are diagonalizable on $\mathcal{V}$, each multiplicity-free and $\mathcal{V}$ is irreducible. Then, $\pi(\text{A}), \pi(\text{A}^*)$ is a Leonard pair [TV03, Definition 1.1]. For further convenience, we denote:

$$
(2.12) \qquad \pi(\mathcal{A}) = \mathcal{A}, \quad \pi(\mathcal{A}^*) = \mathcal{A}^*.
$$

Let {$\theta_M$}$_{M=0}^{2s}$ denote the eigenvalue sequence associated with $A$ (resp. {$\theta_N^*$}$_{N=0}^{2s}$ the eigenvalue sequence associated with $A^*$). By definition of a Leonard pair [T01], there exists an eigenbasis with vectors {$|\theta_M\rangle$}$_{M=0}^{2s}$ (resp. an eigenbasis with vectors {$|\theta_N^*\rangle$}$_{N=0}^{2s}$) such that

$$
(2.13) \qquad A|\theta_M\rangle = \theta_M|\theta_M\rangle, \quad A^*|\theta_M\rangle = A_{M+1,M}^{(*,*)}|\theta_{M+1}\rangle + A_{M,M}^{(*,*)}|\theta_M\rangle + A_{M-1,M}^{(*,*)}|\theta_{M-1}\rangle,
$$

$$
(2.14) \qquad A^* |\theta_N^*\rangle = \theta_N^* |\theta_N^*\rangle, \quad A |\theta_N^*\rangle = A_{N+1,N}^{(*)} |\theta_{N+1}^*\rangle + A_{N,N}^{(*)} |\theta_N^*\rangle + A_{N-1,N}^{(*)} |\theta_{N-1}^*\rangle.
$$

Here $A_{-1,0}^{(*,*)} = A_{2s+1,2s}^{(*,*)} = A_{-1,0}^{(\cdot,*)} = A_{2s+1,2s}^{(\cdot,*)} = 0$ and all other coefficients are non-zero, the representation being irreducible. Given a Leonard pair, the transition matrix $P^{\{\cdot,\cdot\}}$ from the basis $\{|\theta_M\rangle\}_{M=0}^{2s}$ to the basis $\{|\theta_N^*\rangle\}_{N=0}^{2s}$ (resp. the inverse transition matrix $P^{\{\cdot,\cdot\}}^{-1}$ from the basis $\{|\theta_N^*\rangle\}_{N=0}^{2s}$ to the basis $\{|\theta_M\rangle\}_{M=0}^{2s}$) is known explicitly in terms of $q$-Racah polynomials [T04]:

$$
(2.15) \qquad |\theta_N^*\rangle = \sum_{M=0}^{2s} P_{MN}^{\{\cdot,\cdot\}} |\theta_M\rangle \quad \text{and} \quad |\theta_M\rangle = \sum_{N=0}^{2s} (P^{\{\cdot,\cdot\}})^{-1}_{NM} |\theta_N^*\rangle
$$

where the coefficients of the transition matrices are given by:

$$
(2.16) \qquad P_{MN}^{\{\cdot,\cdot\}} = k_N R_M(\theta_N^*) \quad \text{and} \quad (P^{\{\cdot,\cdot\}})^{-1}_{NM} = \nu_0^{-1} k_M^* R_M(\theta_N^*) .
$$
More details on the explicit expressions for the coefficients {$A_{M\pm1,M}^{(\pm,\cdot)}, A_{M,M}^{(\pm,\cdot)}$} and {$A_{N\pm1,N}^{(\pm,\cdot)}, A_{N,N}^{(\pm,\cdot)}$}, the coefficients $\nu_0, k_N, k_M^*$ and the q-Racah polynomials $R_M(\theta_N^*)$ are given in Subsection 3.1.

Let $\tilde{\mathcal{V}}$ be the vector space dual of $\mathcal{V}$, i.e. the vector space $\tilde{\mathcal{V}}$ of all linear functionals from $\mathcal{V}$ to $\mathbb{C}$. Define the family of covectors {$\langle\theta_M|\rangle_{M=0}^{2s} \in \tilde{\mathcal{V}}$ (resp. {$\langle\theta_M^*|\rangle_{N=0}^{2s} \in \tilde{\mathcal{V}}$} associated with the eigenvalue sequence {$\theta_M\rangle_{M=0}^{2s}$ (resp. {$\theta_N^*\rangle_{N=0}^{2s}$})). One has:

$$ (2.17) \qquad \langle \theta_M | A \rangle = \langle \theta_M | \theta_M \rangle, \quad \langle \theta_M | A^* \rangle = \langle \theta_{M+1} | \tilde{A}_{M+1,M}^{(\pm,\cdot)} + \langle \theta_M | \tilde{A}_{M,M}^{(\pm,\cdot)} + \langle \theta_{M-1} | \tilde{A}_{M-1,M}^{(\pm,\cdot)} \rangle, $$

$$ (2.18) \qquad \langle \theta_N^* | A^* \rangle = \langle \theta_N^* | \theta_N^* \rangle, \quad \langle \theta_N^* | A \rangle = \langle \theta_{N+1}^* | \tilde{A}_{N+1,N}^{(\pm,\cdot)} + \langle \theta_N^* | \tilde{A}_{N,N}^{(\pm,\cdot)} + \langle \theta_{N-1} | \tilde{A}_{N-1,N}^{(\pm,\cdot)} \rangle. $$

Introduce the scalar product $\langle .| . \rangle : \tilde{\mathcal{V}} \times \mathcal{V} \to \mathbb{C}$. From the orthogonality relations

$$ (2.19) \qquad \langle \theta_{M'} | \theta_M \rangle = \delta_{MM'} \xi_M, \quad \langle \theta_{N'}^* | \theta_N^* \rangle = \delta_{NN'} \xi_N^*, $$

where $\delta_{ij}$ denotes the Kronecker symbol, one gets:

$$ (2.20) \qquad \tilde{A}_{M',M}^{(\pm,\cdot)} = A_{M,M'}^{(\pm,\cdot)} \frac{\xi_M}{\xi_{M'}} \quad \text{for } M' = M+1, M, M-1, $$

$$ (2.21) \qquad \tilde{A}_{N',N}^{(\pm,\cdot)} = A_{N,N'}^{(\pm,\cdot)} \frac{\xi_N}{\xi_{N'}} \quad \text{for } N' = N+1, N, N-1. $$

The Leonard pair $(A, A^*)$ associated with the AW relations (2.9), (2.10) is known to be of q-Racah type: the eigenvalue sequences are of the form [T99, Theorem 4.4 (case I)]:

$$ (2.22) \qquad \theta_M = b q^{2M} + c q^{-2M}, \quad \theta_N^* = b^* q^{2N} + c^* q^{-2N}, $$

where $b, c, b^*, c^* \in \mathbb{C}^*$ are such that multiplicities in the spectra do not occur. More details are given in Subsection 3.1.

### 2.2.2. Off-shell and on-shell Bethe states.
In the algebraic Bethe ansatz setting, the starting point of the construction of Bethe states associated with the Leonard pair $(A, A^*)$ is the identification of so-called reference states. Let $m_0$ be an integer. By [BaP19, Propositions 3.1, 3.2], if the parameter $\alpha$ is such that:

$$ (2.23) \qquad (q^2 - q^{-2}) \chi^{-1} \alpha c^* q^{m_0} = 1 \quad (\text{resp. } (q^2 - q^{-2}) \chi^{-1} \alpha b q^{-m_0} = -1) $$

then

$$ (2.24) \qquad \pi(\mathcal{C}^+(u, m_0))|\theta_0^*\rangle = 0 \quad (\text{resp. } \pi(\mathcal{C}^-(u, m_0))|\theta_0\rangle = 0). $$

A correspondence between the references states $|\Omega^{\pm}\rangle$ and the fundamental eigenvectors for Leonard pairs follows [BP22, Definition 3.1]:

$$ (2.25) \qquad |\theta_0\rangle = |\Omega^-angle, \quad |\theta_0^*\rangle = |\Omega^+\rangle. $$

More generally, different types of Bethe states are built from successive actions of the dynamical operators $\mathcal{B}^{\pm}(u, m)$ on each reference state $|\Omega^{\pm}\rangle$. Namely, consider the string of dynamical operators

$$ (2.26) \qquad B^{\epsilon}(\bar{u}, m, M) = \mathcal{B}^{\epsilon}(u_1, m + 2(M-1)) \cdots \mathcal{B}^{\epsilon}(u_M, m) $$

where we denote $\bar{u}$ the set of variables $\bar{u} = \{u_1, u_2, \dots, u_M\}$.

**Definition 2.1.** The Bethe states generated from the Leonard pair $(A, A^*)$ are $(m_0 \in \mathbb{N})$:

$$ (2.27) \qquad |\Psi_M^-(\bar{u}, m_0)\rangle = \pi(B^-(\bar{u}, m_0, M))|\Omega^-\rangle \quad \text{for } (q^2 - q^{-2})\chi^{-1}\alpha b q^{-m_0} = -1 \quad \text{and} \quad \beta = 0, $$

$$ (2.28) \qquad |\Psi_+^M(\bar{w}, m_0)\rangle = \pi(B^+(\bar{w}, m_0, M))|\Omega^+\rangle \quad \text{for } (q^2 - q^{-2})\chi^{-1}\alpha c^* q^{m_0} = 1 \quad \text{and} \quad \beta = 0. $$

As usual, if the set of variables $\bar{u}$ or $\bar{w}$ is arbitrary, the Bethe states (2.27), (2.28) are called ‘off-shell’. If the set of variables $\bar{u}$ or $\bar{w}$ satisfies certain Bethe ansatz equations, the Bethe states are called ‘on-shell’.

**Remark 1.** The Bethe states of Definition 2.1 associated with $\beta = 0$ are suitable for studying integrable models with an Hamiltonian of the form $A$ (or $A^*$) or $\kappa A + \kappa^* A^*$. See e.g. [BaP19, Section 5.2]. To handle an Hamiltonian of the so-called Heun-Askey-Wilson form $\kappa A + \kappa^* A^* + \kappa_+[A, A^*]_q + \kappa_-[A^*, A]_q$, the case $\beta \neq 0$ has to be considered instead.
A correspondence between on-shell Bethe states (associated with either homogeneous Bethe equations or inhomogeneous Bethe equations) and eigenvectors of Leonard pairs is derived in [BP22]. We refer the reader to this work for more details. Recall the functions $f(u, v)$ and $h(u, v)$ in Appendix A and introduce the notations $(\epsilon = \pm 1, \zeta \in \mathbb{C}^*)$:

$$ (2.29) \qquad \Lambda_1^\epsilon(u) = \frac{q^{-2s-1}}{u^\epsilon} (q^{2s+1}u\zeta^{-1} - u^{-1}\zeta) (q^{2s+1}u\zeta - u^{-1}\zeta^{-1}) \\ \qquad \times (u\mathbf{c}^*q^{-2s} + u^{-1}\mathbf{b}q^{2s}) \left( u \left( \frac{\mathbf{c}}{\mathbf{c}^*} \right)^{\frac{1-\epsilon}{2}} + u^{-1} \left( \frac{\mathbf{b}^*}{\mathbf{b}} \right)^{\frac{1+\epsilon}{2}} \right), $$

$$ (2.30) \qquad \Lambda_2^\epsilon(u) = \frac{(u^2 - u^{-2})q^{-2s-1}}{u^\epsilon(qu^2 - q^{-1}u^{-2})} (q^{2s-1}u^{-1}\zeta - u\zeta^{-1}) (q^{2s-1}u^{-1}\zeta^{-1} - u\zeta) \\ \qquad \times (q^2ubq^{2s} + u^{-1}\mathbf{c}^*q^{-2s}) \left( q^2u \left( \frac{\mathbf{b}^*}{\mathbf{b}} \right)^{\frac{1+\epsilon}{2}} + u^{-1} \left( \frac{\mathbf{c}}{\mathbf{c}^*} \right)^{\frac{1-\epsilon}{2}} \right). $$

Note that $\Lambda_i^\epsilon(u)$, $i = 1, 2$, are the eigenvalues associated with the dynamical operators $\mathcal{A}^\epsilon(u, m)$, $\mathcal{D}^\epsilon(u, m)$ acting on the reference state $|\Omega^\epsilon\rangle$ [BaP19, Lemma 3.3], respectively. As will be detailed in Subsections 3.1, 3.2, the parameter $\zeta$ finds a natural interpretation with respect to the parameter sequences associated with a Leonard triple.

* On-shell Bethe states of homogeneous type:

Define the set of functions:

$$ (2.31) \qquad E_{\pm}^{M}(u_i, \bar{u}_i) = - \frac{b(u_i^2)}{b(qu_i^2)} \prod_{j=1, j \neq i}^{M} f(u_i, u_j) \Lambda_1^{\pm}(u_i) + \prod_{j=1, j \neq i}^{M} h(u_i, u_j) \Lambda_2^{\pm}(u_i), $$

for $i = 1, \dots, M$, where we denote $\bar{u}_i = \bar{u} \setminus u_i$. For each integer $M$ (resp. $N$) with $0 \le M, N \le 2s$, assume there exists at least one set of non trivial admissible$^4$ Bethe roots $S_-^{M(h)} = \{u_1, \dots, u_M\}$ (resp. $S_+^{*N(h)} = \{w_1, \dots, w_N\}$) such that

$$ (2.32) \qquad E_-^M(u_i, \bar{u}_i) = 0 \quad \text{for} \quad \bar{u} = S_-^{M(h)}, \quad (\text{resp.} \quad E_+^N(u_i, \bar{u}_i) = 0 \quad \text{for} \quad \bar{u} = S_+^{*N(h)}). $$

**Definition 2.2.** The equations $E_-^M(u_i, \bar{u}_i) = 0$ (resp. $E_+^N(u_i, \bar{u}_i) = 0$) for $i = 1, \dots, M$ (resp. $i = 1, \dots, N$) are called the Bethe ansatz equations of homogeneous type associated with the set of Bethe roots $\bar{u} = S_-^{M(h)}$ (resp. $\bar{u} = S_+^{*N(h)}$).

As shown in [BP22, Lemma 3.5] the eigenvectors of $A$ (resp. $A^*$) can be written as on-shell Bethe states of homogeneous type:

$$ (2.33) \qquad |\theta_M\rangle = \mathcal{N}_M(\bar{u})|\Psi_-^M(\bar{u}, m_0)\rangle \quad \text{for} \quad \bar{u} = S_-^{M(h)}, $$

$$ (2.34) \qquad |\theta_N^*\rangle = \mathcal{N}_N^*(\bar{w})|\Psi_+^N(\bar{w}, m_0)\rangle \quad \text{for} \quad \bar{w} = S_+^{*N(h)} $$

with normalization factors $\mathcal{N}_0(\cdot) = \mathcal{N}_0^*(\cdot) = 1$ and

$$ (2.35) \qquad \mathcal{N}_M(\bar{u}) = \prod_{k=1}^{M} (qu_k b(u_k^2) A_{k,k-1}^{(*,*)})^{-1}, \qquad \mathcal{N}_N^*(\bar{w}) = \prod_{k=1}^{N} (-q^{-1}w_k^{-1}b(w_k^2)A_{k,k-1}^{(*,*)})^{-1}, $$

where (2.13), (2.14).

<footnote>$^4$Note that the set of equations $E_\pm^M(u_i, \bar{u}_i) = 0$ for $i = 1, \dots, M$ contain trivial solutions where $u_i^2 = u_i^{-2}$ or $u_i = 0$ for $i = 1, \dots, M$, or solutions such that $U_i = U_j$ where $U_i = (qu_i^2 + q^{-1}u_i^{-2})/(q + q^{-1})$. Such solutions are not admissible. For further informations, see [BaP19, BP22] and [Sl22, Subsection 2.3].</footnote>
• On-shell Bethe states of inhomogeneous type:

Define the functions:

$$
(2.36) \quad E_{\pm}(u_i, \bar{u}_i) = \frac{b(u_i^2)}{b(qu_i^2)} u_i^{\pm 1} \prod_{j=1, j \neq i}^{2s} f(u_i, u_j) \Lambda_1^{\pm}(u_i) - (q^2 u_i^3)^{\mp 1} \prod_{j=1, j \neq i}^{2s} h(u_i, u_j) \Lambda_2^{\pm}(u_i) \\
\qquad + \nu_{\pm} \frac{u_i^{\mp 2} b(u_i^2)}{b(q)} \frac{\prod_{k=0}^{2s} b(q^{1/2+k-s} \zeta u_i) b(q^{1/2+k-s} \zeta^{-1} u_i)}{\prod_{j=1, j \neq i}^{2s} b(u_i u_j^{-1}) b(q u_i u_j)},
$$

where

$$
(2.37) \qquad \nu_+ = q^{-1-4s}c^*, \quad \nu_- = q^{1+4s}b,
$$

for $i = 1, \dots, 2s$. For each integer $M$ (or $N$) with $0 \le M, N \le 2s$, there exists at least one set of non trivial admissible Bethe roots $S_+^{M(i)} = \{u_1, \dots, u_{2s}\}$ (resp. $S_-^{*N(i)} = \{w_1, \dots, w_{2s}\}$) such that

$$
(2.38) \quad E_+(u_j, \bar{u}_j) = 0 \quad \text{for} \quad \bar{u} = S_+^{M(i)} \quad \text{(resp. } E_-(w_j, \bar{w}_j) = 0 \quad \text{for} \quad \bar{w} = S_-^{*N(i)}),
$$

and

$$
(2.39) \quad \theta_M = q^{-4s} \left( c^* (\zeta^2 + \zeta^{-2}) [2s]_q + q^{2s} (bq^{2s} + cq^{-2s}) - q c^* \sum_{j=1}^{2s} (qu_j^2 + q^{-1}u_j^{-2}) \right) \quad \text{for } \bar{u} = S_+^{M(i)}
$$

$$
(2.40) \quad (\text{resp. } \theta_N^* = q^{4s} (b(\zeta^2 + \zeta^{-2})[2s]_q + q^{-2s}(b^*q^{2s} + c^*q^{-2s}) - q^{-1}b \sum_{j=1}^{2s} (qw_j^2 + q^{-1}w_j^{-2})) \\
\qquad \text{for } \bar{w} = S_-^{*N(i)}).
$$

**Definition 2.3.** The equations $E_+(u_i, \bar{u}_i) = 0$ (resp. $E_-(u_i, \bar{u}_i) = 0$) for $i = 1, \dots, 2s$ are called the Bethe ansatz equations of inhomogeneous type associated with the set of Bethe roots $\bar{u} = S_+^{M(i)}$ (resp. $\bar{u} = S_-^{*N(i)}$).

As shown in [BP22, Lemma 3.7] the eigenvectors of $A$ (resp. $A^*$) can be written as on-shell Bethe states of
inhomogeneous type:

$$
(2.41) \qquad |\theta_M\rangle = \mathcal{N}_M^{(i)}(\bar{u}') |\Psi_+^{2s}(\bar{u}', m_0)\rangle \quad \text{for } \bar{u}' = S_+^{M(i)},
$$

$$
(2.42) \qquad |\theta_N^*\rangle = \mathcal{N}_N^{*(i)}(\bar{w}') |\Psi_-^{2s}(\bar{w}', m_0)\rangle \quad \text{for } \bar{w}' = S_-^{*N(i)}
$$

with

$$
(2.43) \qquad \mathcal{N}_M^{(i)}(\bar{u}') = \mathcal{N}_{2s}^*(\bar{u}') (P^{\{\cdot,\cdot\}})^{-1}_{2s,M}, \qquad \mathcal{N}_N^{*(i)}(\bar{w}') = \mathcal{N}_{2s}^*(\bar{w}') P_{2s,N}^{\{\cdot,\cdot\}},
$$

where (2.35) and (2.16) are used.

3. SCALAR PRODUCTS OF BETHE STATES FROM THE LEONARD TRIPLE

In this section, the theory of Leonard triples that extends the theory of Leonard pairs is applied to the
derivation of certain scalar products of off-shell Bethe states (2.27), (2.28) with eigenvectors $|\theta_M\rangle$. To this
end, the Leonard triple $(A, A^*, A^\diamond)$ is built from the Leonard pair $(A, A^*)$, and the transition matrices relating
the respective eigenbases associated with the triple are given. Using the properties of $A^\diamond$, two different types
of normalized scalar products of Bethe states are computed. Specializing those scalar products, we obtain an
explicit relation characterizing Bethe roots of inhomogenous type in terms of Bethe roots of homogeneous type.

3.1. **The Leonard triple (A, A*, A^①).** The concept of Leonard triple is a natural generalization of the concept of Leonard pair. According to [Curt07, Definition 1.2], a Leonard triple on V is a 3-tuple of maps in End(V) such that for each map, there exists a basis of V with respect to which the matrix representing that map is diagonal and the matrices representing the other two maps are irreducible tridiagonal.

A Leonard triple of $q$-Racah type $(A, A^*, A^\diamond)$ is constructed as follows. Recall the Leonard pair of $q$-Racah type $(A, A^*)$ with spectra (2.22) introduced in Subsection 2.2.1. Denote by $\mathbb{I}$ the $(2s+1) \times (2s+1)$ identity matrix. Introduce the notations:

$$
(3.1) \qquad \theta_M^a = b^a q^{2M} + c^a q^{-2M}, \quad M = 0, 1, \dots, 2s, \quad a \in \{\cdot, *, \diamond\}
$$
with $b^a, c^a \in \mathbb{C}^*$ and

$$ (3.2) \quad \omega^{\{a,b,c\}} = -(q-q^{-1})^2 \left( \theta_s^a \theta_s^b - \frac{1}{r_0} (q^{2s+1} + q^{-2s-1}) \theta_s^c \right) \quad \text{with} \quad r_0^{-2} = bc = b^*c^* = b^\diamond c^\diamond. $$

Define

$$ (3.3) \qquad A^{\diamond} = \frac{r_0}{(q^2 - q^{-2})} [A^*, A]_q + \frac{r_0 \omega^{\{\cdot,*,\diamond\}}}{(q - q^{-1})(q^2 - q^{-2})} \mathbb{I}. $$

By [TV03, Theorem 1.5] the Leonard pair $(A, A^*)$ satisfies the image by $\pi$ of the Askey-Wilson relations (2.9)-(2.10) for some structure constants. The proof of the following Lemma is straightforward, thus we skip it.

**Lemma 3.1.** Assume that the structure constants in (2.9)-(2.10) take the form:

$$ (3.4) \qquad \rho = - \frac{(q^2 - q^{-2})^2}{r_0^2}, \quad \omega = \omega^{\{\cdot,*,\diamond\}}, \quad \eta = \frac{(q + q^{-1})}{r_0} \omega^{\{\cdot,\diamond,\cdot\}}, \quad \eta^* = \frac{(q + q^{-1})}{r_0} \omega^{\{*,\diamond,\cdot\}}. $$

The triple $(A, A^*, A^\diamond)$ satisfy the Askey-Wilson relations:

$$ (3.5) \qquad [A^a, [A^a, A^b]]_q^{q-1} = - \frac{(q^2 - q^{-2})^2}{r_0^2} A^b + \omega^{\{a,b,c\}} A^a + \frac{(q + q^{-1})}{r_0} \omega^{\{a,c,b\}} \mathbb{I} \quad \text{for any } a \neq b \in \{\cdot,*,\diamond\}. $$

**Remark 2.** The Askey-Wilson relations (3.5) for $(a, b, c) \in \{((., \cdot, \diamond), (\cdot, \cdot, \diamond)), (\diamond, \cdot, \cdot), (\cdot, \cdot, \diamond)\}$ together with the relation (3.3) can be alternatively put into the $\mathbb{Z}_3$-symmetric form:

$$ (3.6) \qquad \frac{r_0}{(q^2 - q^{-2})} [A^b, A^a]_q - A^c + \frac{r_0 \omega^{\{a,b,c\}}}{(q - q^{-1})(q^2 - q^{-2})} = 0 \quad \text{for } (a,b,c) \in \{((., \cdot, \diamond), (\cdot, \cdot, \diamond), (\diamond, \cdot, \cdot)), (\diamond, \cdot, \cdot)\}. $$

Note that the Askey-Wilson relations considered in [H11, eqs. (27)-(29)] are obtained from the following specialization of (3.6):

$$ (3.7) \qquad \begin{aligned} A|_{q \to q^{-1}} &\to A, & A^*|_{q \to q^{-1}} &\to A^*, & A^\diamond|_{q \to q^{-1}} &\to A^\epsilon, & 2s \to d, & r_0 \to 1, \\ b &\to a^{-1}q^{2s}, & c &\to aq^{-2s}, & b^* &\to b^{-1}q^{2s}, & c^* &\to bq^{-2s}, & b^\diamond &\to -c^{-1}q^{2s}, & c^\diamond &\to -cq^{-2s}. \end{aligned} $$

**Remark 3.** The relationship between the algebraic Bethe ansatz parameter $\zeta$ introduced in Subsection 2.2 (see also [BaP19, BP22]) and the parameters $b^\diamond, c^\diamond$ in (3.1) is determined through a comparison between the structure constants [BP22, eqs. (2.12)-(2.14)] and the structure constants in (3.5) for $(a, b, c) = (\cdot, \cdot, \diamond)$. One finds:

$$ (3.8) \qquad \zeta^{-2} = b^{\diamond} r_0 q^{2s} \quad \text{and} \quad \zeta^2 = c^{\diamond} r_0 q^{-2s}. $$

**3.1.1.** Eigenbases and respective actions of the Leonard triple. We now turn to the explicit construction of three eigenbases $\{|\theta_M^a\rangle\}_{M=0}^{2s}$, $\{|\theta_M^b\rangle\}_{M=0}^{2s}$, $\{|\theta_M^c\rangle\}_{M=0}^{2s}$ and respective actions associated with the triple $(A^a, A^b, A^c)$ for the sets $(a, b, c) \in \{(\cdot, \cdot, \diamond), (\cdot, \diamond, \cdot), (\diamond, \cdot, \cdot)\}$.

**Lemma 3.2.** Assume $(A^a, A^b)$ is a Leonard pair of $q$-Racah type with spectra (3.1). There exists eigenbases $\{|\theta_M^a\rangle\}_{M=0}^{2s}$ and $\{|\theta_M^b\rangle\}_{M=0}^{2s}$ such that:

$$ (3.9) \qquad A^a |\theta_M^a\rangle = \theta_M^a |\theta_M^a\rangle, \quad A^b |\theta_M^a\rangle = A_{M+1,M}^{(b,a)} |\theta_{M+1}^a\rangle + A_{M,M}^{(b,a)} |\theta_M^a\rangle + A_{M-1,M}^{(b,a)} |\theta_{M-1}^a\rangle, $$

$$ (3.10) \qquad A^b |\theta_M^b\rangle = \theta_M^b |\theta_M^b\rangle, \quad A^a |\theta_M^b\rangle = A_{M+1,M}^{(a,b)} |\theta_{M+1}^b\rangle + A_{M,M}^{(a,b)} |\theta_M^b\rangle + A_{M-1,M}^{(a,b)} |\theta_{M-1}^b\rangle $$

with

$$ (3.11) \qquad A_{M,M-1}^{(b,a)} = \frac{q^{2-4s} \frac{(1-q^{2M})(c^a - b^a q^{2M+4s})(b^b b^c r_0 q^{4s-1} + b^a q^{2M-2})(c^a c^c r_0 q^{-1} + c^b q^{2M-2})}{(c^a - b^a q^{4M-2})(c^a - b^a q^{4M})}}{(c^a - b^a q^{4M-2})(c^a - b^a q^{4M})}, $$

$$ (3.12) \qquad A_{M-1,M}^{(b,a)} = \frac{(1 - q^{2M-4s-2})(c^a - b^a q^{2M-2})(c^a + b^b b^c r_0 q^{2M+4s-1})(c^b + b^a c^c r_0 q^{2M-1})}{(c^a - b^a q^{4M-4})(c^a - b^a q^{4M-2})}, $$

$$ (3.13) \qquad A_{M,M}^{(b,a)} = \theta_0^b - A_{M,M+1}^{(b,a)} - A_{M,M-1}^{(b,a)}, $$

and

$$ (3.14) \qquad A_{M,N}^{(a,b)} = A_{M,N}^{(b,a)} |b^a \leftrightarrow b^b, c^a \leftrightarrow c^b, b^c q^{4s} \leftrightarrow c^c|. $$
Proof. We adapt the results of [T04] to our conventions: in [T04, eqs. (134)-(137)], firstly we substitute $q \to q^2$, $\theta_0 \to h(1+sq^2)$, $\theta_0^* \to h^*(1+s^*q^2)$ and replace

$$ (3.15) \qquad h \to c^a, \quad h^* \to c^b, \quad d \to 2s, \quad s \to \frac{b^a}{c^a}q^{-2}, \quad s^* \to \frac{b^b}{c^b}q^{-2}, $$

$$ (3.16) \qquad r_1 \to -\frac{b^b b^c}{c^a} r_0 q^{4s-1}, \quad r_2 \to -\frac{b^a c^c}{c^b} r_0 q^{-1}. $$

Then, according to the above substitutions, expressions for the entries of $A^a$ (resp. $A^b$) in the eigenbasis of $A^b$ (resp. $A^a$) are computed using [T04, bases $d^*0*0d$ and $d00*d^*$ page 23]. $\square$

We now compute the action of $A^c$ as defined in (3.6) on the eigenbases $\{|\theta_M^a\rangle\}_{M=0}^{2s}$ and $\{|\theta_M^b\rangle\}_{M=0}^{2s}$ such that (3.9), (3.10) hold. Below, the proof of the first lemma is given. The second lemma is proven along the same line.

**Lemma 3.3.** Assume that $(A^c, A^a)$ is a Leonard pair of q-Racah type with spectra (3.1). Then

$$ (3.17) \qquad A^c |\theta_M^a\rangle = A_{M+1,M}^{(c,a)} \frac{f_{M+1}}{f_M} |\theta_{M+1}^a\rangle + A_{M,M}^{(c,a)} |\theta_M^a\rangle + A_{M-1,M}^{(c,a)} \frac{f_{M-1}}{f_M} |\theta_{M-1}^a\rangle, $$

where

$$ (3.18) \qquad f_M = f_0 \prod_{k=0}^{M-1} f(k), \quad f(k) = \frac{c^a q^{-2k-2} + r_0 q^{-1} c^b b^c}{b^a q^{2k} + r_0 q^{-1} c^b b^c} $$

for some scalar $f_0 \in \mathbb{C}^*$. One sets $f_{-1} = 0$.

*Proof.* On one hand, using (3.6) we find:

$$ (3.19) \qquad A^c |\theta_M^a\rangle = \hat{A}_{M+1,M}^{(c,a)} |\theta_{M+1}^a\rangle + \hat{A}_{M,M}^{(c,a)} |\theta_M^a\rangle + \hat{A}_{M-1,M}^{(c,a)} |\theta_{M-1}^a\rangle, $$

where

$$ (3.20) \qquad \hat{A}_{M+1,M}^{(c,a)} = q^{-2M-1} r_0 c^a A_{M+1,M}^{(b,a)}, \quad \hat{A}_{M-1,M}^{(c,a)} = q^{2M-1} r_0 b^a A_{M-1,M}^{(b,a)}, $$

$$ (3.21) \qquad \hat{A}_{M,M}^{(c,a)} = \frac{r_0 \omega^{\{a,b\}}}{(q-q^{-1})(q^2-q^{-2})} + \frac{r_0}{q+q^{-1}} \theta_M^a A_{M,M}^{(b,a)}. $$

On the other hand, since $(A^c, A^a)$ is a Leonard pair, there exists a basis of $A^a$ with eigenvectors $|\theta_M^a\rangle = f_M|\theta_M^a\rangle$ such that:

$$ (3.22) \qquad A^a |\bar{\theta}_M^a\rangle = \theta_M^a |\bar{\theta}_M^a\rangle, \quad A^c |\bar{\theta}_M^a\rangle = A_{M+1,M}^{(c,a)} |\bar{\theta}_M^a\rangle + A_{M,M}^{(c,a)} |\bar{\theta}_M^a\rangle + A_{M-1,M}^{(c,a)} |\bar{\theta}_M^a\rangle. $$

Comparing the second equation in (3.22) to (3.17), after simplifications one gets the desired coefficients. In particular, the identity relating the diagonal coefficients in (3.17) and (3.22) leads to the non-trivial relation

$$ (3.23) \qquad \frac{r_0 \omega^{\{a,b,c\}}}{(q-q^{-1})(q^2-q^{-2})} + \frac{r_0}{q+q^{-1}} \theta_M^a \left( \theta_0^b - A_{M,M+1}^{(b,a)} - A_{M,M-1}^{(b,a)} \right) = \theta_0^c - A_{M,M+1}^{(c,a)} - A_{M,M-1}^{(c,a)}, $$

which is readily checked. $\square$

**Lemma 3.4.** Assume that $(A^b, A^c)$ is a Leonard pair of q-Racah type with spectra (3.1). Then

$$ (3.24) \qquad A^c |\theta_M^b\rangle = A_{M+1,M}^{(c,b)} \frac{g_{M+1}}{g_M} |\theta_{M+1}^b\rangle + A_{M,M}^{(c,b)} |\theta_M^b\rangle + A_{M-1,M}^{(c,b)} \frac{g_{M-1}}{g_M} |\theta_{M-1}^b\rangle, $$

where

$$ (3.25) \qquad g_M = g_0 \prod_{k=0}^{M-1} g(k), \quad g(k) = q^{4k} \frac{b^b c^b q^{-2k} + r_0 q c^a b^c}{c^b b^b q^{2k} + r_0 q^{-1} c^a b^c}. $$

for some scalar $g_0 \in \mathbb{C}^*$. One sets $g_{-1} = 0$.

We now introduce the eigenbasis $\{|\theta_M^c\rangle\}_{M=0}^{2s}$ of $A^c$ and give the action of $A^b$ and $A^a$ on its vectors. Adapting the results in [T04, eqs. (134)-(137)] to our conventions, one gets the following lemma.
**Lemma 3.5.** Assume that $(A^b, A^c)$ is a Leonard pair of $q$-Racah type. There exists an eigenbasis $\{|\theta_M^c\rangle\}_{M=0}^{2s}$ of $A^c$ such that:

$$ (3.26) \qquad A^b |\theta_M^c\rangle = A_{M+1,M}^{(b,c)} |\theta_{M+1}^c\rangle + A_{M,M}^{(b,c)} |\theta_M^c\rangle + A_{M-1,M}^{(b,c)} |\theta_{M-1}^c\rangle, $$

where $A_{M,N}^{(b,c)}$ is given by (3.14) for $(a, b, c) \to (b, c, a)$.

The previous lemma fixes the normalizations of the eigenvectors $|\theta_M^c\rangle$. For the proof of the following lemma, we use (3.6) with the substitution $(a, b, c) \to (b, c, a)$ and (3.26), and follow the same reasoning as for Lemmas 3.3, 3.4.

**Lemma 3.6.** Assume that $(A^c, A^a)$ is a Leonard pair of $q$-Racah type with spectra (3.1). Then

$$ (3.27) \qquad A^a |\theta_M^c\rangle = A_{M+1,M}^{(a,c)} \frac{h_{M+1}}{h_M} |\theta_{M+1}^c\rangle + A_{M,M}^{(a,c)} |\theta_M^c\rangle + A_{M-1,M}^{(a,c)} \frac{h_{M-1}}{h_M} |\theta_{M-1}^c\rangle, $$

where

$$ (3.28) \qquad h_M = h_0 \prod_{k=0}^{M-1} \mathfrak{h}(k), \quad \mathfrak{h}(k) = q^{4k} \frac{\mathfrak{b}^c \mathfrak{c}^c q^{-2k} + r_0 \mathfrak{q} \mathfrak{c}^b \mathfrak{b}^a}{\mathfrak{c}^c \mathfrak{b}^c q^{2k} + r_0 \mathfrak{q}^{-1} \mathfrak{c}^b \mathfrak{b}^a} $$

for some scalar $h_0 \in \mathbb{C}^*$. One sets $h_{-1} = 0$.

**Remark 4.** The structure constants in (3.5) are recovered by applying [TV03, Theorem 4.5 and 5.3] for the spectra (3.1) and diagonal entries (3.13).

Up to now, we have assumed that $(A^a, A^b)$, $(A^b, A^c)$ and $(A^c, A^a)$ are Leonard pairs. In order to show that $(A^a, A^b, A^c)$ is a Leonard triple, it remains to identify the conditions such that those assumptions hold.

**Proposition 3.1.** If the following conditions

(i) $b^a(c^a)^{-1} \neq q^{-2M}$, $b^b(c^b)^{-1} \neq q^{-2M}$, $b^c(c^c)^{-1} \neq q^{-2M}$, for $1 \le M \le 4s-1$,

(ii) $c^b \neq -r_0 q^{-2M+1} c^a c^c$, $b^a \neq -r_0 q^{-2M+1+4s} b^b b^c$, $c^a \neq -r_0 q^{2M+4s-1} b^b b^c$, $c^b \neq -r_0 q^{2M-1} b^a c^c$, for $1 \le M \le 2s$,

are satisfied, then $(A^a, A^b, A^c)$ is a Leonard triple of $q$-Racah type.

*Proof.* Firstly, we identify the conditions (i) that ensure the spectra (3.1) have no multiplicities. From (3.1), one has:

$$ \theta_M^a - \theta_N^a = (q^{M-N} - q^{-M+N})(b^a q^{M+N} - c^a q^{-M-N}). $$

Recall $q$ is not a root of unity. Observe that this latter equation is vanishing only if the conditions (i) are not satisfied. Now assume (i). Secondly, we identify the conditions (ii) such that for each map (e.g. $A^a$), with respect to its eigenbasis (e.g. $\{|\theta_M^a\rangle\}_{M=0}^{2s}$) the matrices representing the other two maps (e.g. $A^b, A^c$) are irreducible tridiagonal. In other words, we identify the relations between the parameters such that off-diagonal upper and lower entries entering in (3.9), (3.10), (3.24), (3.26) and (3.17), (3.27) are not vanishing. Assuming (i), (ii) are satisfied, we conclude $(A^a, A^b, A^c)$ is a Leonard triple. $\square$

**Remark 5.** The conditions (i), (ii) specialize to [H11, Definition 7.1 (RQRAC3), (RQRAC4)] using (3.7).

3.1.2. Transition matrices. From the theory of Leonard pairs, it is known that the transition matrix relating the two different eigenbases can be expressed in terms of orthogonal polynomials of the Askey-scheme [T04]. In previous subsection, three different eigenbases associated with a Leonard triple $(A^a, A^b, A^c)$ have been introduced. Below, the transition matrices relating the three eigenbases are given explicitly in terms of $q$-Racah polynomials. For convenience, we denote:

$$ (3.29) \qquad P_{MN}^{\{a,b\}} = k_N^{\{a,b,c\}} R_M^{\{a,c\}}(\theta_N^b) \quad \text{and} \quad (P^{\{a,b\}})^{-1}_{NM} = \nu_0^{\{a,b,c\}}^{-1} k_M^{\{b,a,c\}} R_M^{\{a,c\}}(\theta_N^b) $$

and

$$ (3.30) \qquad R_M^{\{a,c\}}(\theta_N^b) = 4\phi_3 \left[ \frac{q^{-2M}, \frac{b^a}{c^a}q^{2M}, q^{-2N}, \frac{b^b}{c^b}q^{2N}}{\frac{b^a c^c}{c^b}r_0q, -\frac{b^b b^c}{c^a}r_0q^{4s+1}, q^{-4s}; q^2, q^2} \right], $$
$$k_N^{\{a,b,c\}} = \frac{\left(-\frac{b^b}{c^a} r_0 q^{4s+1}, -\frac{b^a}{c^b} r_0 q, \frac{b^b}{c^c} q^{-4s}; q^2\right)_N}{(q^2, -\frac{b^b}{c^a} r_0 q, -\frac{c^a}{c^b} r_0 q^{1-4s}, \frac{b^b}{c^b} q^{4s+2}; q^2)_N} \frac{(1-\frac{b^b}{c^b} q^{4N})}{(\frac{b^a}{c^a})^N (1-\frac{b^b}{c^b})}, \quad k_M^{\{b,a,c\}} = k_M |_{b^a \leftrightarrow b^b, c^a \leftrightarrow c^b, b^b q^{4s} \leftrightarrow c^c}$$

$$\nu_0^{\{a,b,c\}} = \frac{\left(\frac{b^a}{c^a} q^2, \frac{b^b}{c^b} q^2; q^2\right)_{2s}}{\left(-\frac{b^b}{c^b} r_0 q^{4s+1}\right)_{2s} \left(-\frac{b^a}{c^a} r_0 q^{1-4s}, -\frac{c^a}{c^b} r_0 q^{1-4s}; q^2\right)_{2s}}.$$

**Remark 6.** The correspondence between the variables entering in the standard definition of the q-Racah polynomial $R_n(\mu(x); \alpha, \beta, \gamma, \delta|q)$ where $\mu(x) = q^{-x} + \gamma\delta q^{x+1}$ in [KS96, eq. (3.2.1)], and the Bethe ansatz parameters $b^a, c^a, s, a \in \{*, \diamond\}$, entering in the Bethe ansatz equations (2.32) and (2.38) with (2.39), (2.40) is given by:

$$\begin{align*}
& q \rightarrow q^2, \quad n \rightarrow M, \quad x \rightarrow N, \\
& \alpha \rightarrow -\frac{b^a c^c}{c^b} r_0 q^{-1}, \quad \beta \rightarrow -\frac{b^c c^b}{c^a} r_0 q^{-1}, \quad \gamma \rightarrow q^{-4s-2}, \quad \delta \rightarrow \frac{b^b}{c^b} q^{4s}.
\end{align*}$$

**Proposition 3.2.** The transition matrix relating the basis $\{|\theta_N^b\rangle\}_{N=0}^{2s}$ and $\{|\theta_M^a\rangle\}_{M=0}^{2s}$ is such that:

$$(3.31) \qquad |\theta_N^b\rangle = \sum_{M=0}^{2s} P_{MN}^{\{a,b\}} |\theta_M^a\rangle \quad \text{and} \quad |\theta_M^a\rangle = \sum_{N=0}^{2s} (P^{\{a,b\}})^{-1}_{NM} |\theta_N^b\rangle;$$

The transition matrix relating the basis $\{|\theta_N^c\rangle\}_{N=0}^{2s}$ and $\{|\theta_M^b\rangle\}_{M=0}^{2s}$ is such that:

$$(3.32) \qquad |\theta_N^c\rangle = \sum_{M=0}^{2s} g_M P_{MN}^{\{b,c\}} |\theta_M^b\rangle \quad \text{and} \quad |\theta_M^b\rangle = \sum_{N=0}^{2s} g_M^{-1} (P^{\{b,c\}})^{-1}_{NM} |\theta_N^c\rangle;$$

The transition matrix relating the basis $\{|\theta_N^a\rangle\}_{N=0}^{2s}$ and $\{|\theta_M^c\rangle\}_{M=0}^{2s}$ is such that:

$$(3.33) \qquad |\theta_N^a\rangle = \sum_{M=0}^{2s} f_N^{-1} h_M P_{MN}^{\{c,a\}} |\theta_M^c\rangle \quad \text{and} \quad |\theta_M^c\rangle = \sum_{N=0}^{2s} f_N h_M^{-1} (P^{\{c,a\}})^{-1}_{NM} |\theta_N^a\rangle,$$

with (3.18), (3.25), (3.28) and (3.29).

**Proof.** Firstly, consider (3.31). With respect to the coefficients (3.9), (3.10), the explicit expressions for the transition matrix with entries $P_{MN}^{\{a,b\}}$ - as well as the inverse transfer matrix' entries - are obtained from the results in [T04, Section 16] adapted to our conventions (3.15), (3.16). Starting from (3.24), (3.26) and (3.17), (3.27), eqs. (3.32) and (3.33) are derived similarly. $\square$

According to Lemmas 3.3, 3.4 and 3.6, the three eigenbases associated with the Leonard triple $(A^a, A^b, A^c)$ are defined up to overall scalars $f_0, g_0, h_0$. Their precise relationship follow from consistency of the eigenbases' construction, as we now show.

**Proposition 3.3.** The scalars $f_0, g_0, h_0$ are such that:

$$(3.34) \qquad f_0 h_0^{-1} g_0^{-1} = (-1)^{2s} q^{2s(2s-1)} \frac{\left(q^2 r_0^2 b^2; q^2\right)_{2s} \left(q^2 r_0^2 (b^{\diamond})^2; q^2\right)_{2s} \left(q^2 r_0^2 (b^{*})^2; q^2\right)_{2s}}{\left(-\frac{q r_0 b b^{\diamond}}{b^*}; q^2\right)_{2s} \left(-\frac{q r_0 b b^{*}}{b^{\diamond}}; q^2\right)_{2s} \left(-\frac{q r_0 b^{\diamond} b^{*}}{b}; q^2\right)_{2s}}.$$

**Proof.** Use $|\theta_M^b\rangle$ given by (3.31) in (3.32), and compare it to (3.33). One finds:

$$(3.35) \qquad f_0 h_0^{-1} g_0^{-1} = \frac{\prod_{k=0}^{N-1} \mathfrak{h}(k) \sum_{M=0}^{2s} P_{M'M}^{\{a,b\}} P_{MN}^{\{b,c\}} \prod_{k=0}^{M-1} \mathfrak{g}(k)}{(P^{\{c,a\}})^{-1}_{M'N} \prod_{k=0}^{M'-1} \mathfrak{f}(k)}$$

for all $M' = 0, \dots, 2s$ and $N = 0, \dots, 2s$. The l.h.s. being independent of $M', N$, we set $M' = N = 0$ in the r.h.s. From (3.29), one has $P_{0M}^{\{a,b\}} = k_M^{\{a,b,c\}}$, $P_{M0}^{\{b,c\}} = 1$ and $(P^{\{c,a\}})^{-1}_{00} = \nu_0^{\{c,a,b\}}^{-1}$. It follows:

$$(3.36) \qquad f_0 h_0^{-1} g_0^{-1} = \nu_0^{\{c,a,b\}} \sum_{M=0}^{2s} k_M^{\{a,b,c\}} \prod_{k=0}^{M-1} \mathfrak{g}(k).$$

Then, using the relation

$$(3.37) \qquad (a; q^{-2})_M = (-a)^M q^{-M(M-1)} (a^{-1}; q^2)_M$$
and $r_0^{-2} = b^a c^a = b^b c^b = b^c c^c$, we obtain:

$$ (3.38) \qquad \prod_{k=0}^{M-1} g(k) = \frac{\left(\frac{-qr_0 b^a b^c}{b^a}; q^2\right)_M}{\left(-\frac{qr_0 b^a c^c}{c^b}; q^2\right)_M} q^{M(M-1)} \left(\frac{q^{-1} r_0 b^c c^a}{b^b}\right)^{-M}. $$

Moreover, $\nu_0^{\{c,a,b\}}$ simplifies to:

$$ (3.39) \qquad \nu_0^{\{c,a,b\}} = \frac{(q^2 r_0^2 (b^a)^2, q^2 r_0^2 (b^c)^2; q^2)_{2s}}{\left(-\frac{qr_0 b^a b^b}{b^c}, -\frac{qr_0 b^b b^c}{b^a}; q^2\right)_{2s}} \frac{1}{\left(-\frac{qr_0 b^a b^c}{b^b}\right)^{2s}}. $$

Using $k_M^{\{a,b,c\}}$ below (3.30) and (3.38), one gets:

$$ (3.40) \qquad \sum_{M=0}^{2s} k_M^{\{a,b,c\}} \prod_{k=0}^{M-1} g(k) = \sum_{M=0}^{2s} \frac{\left(\frac{-qr_0 b^a b^c}{c^b} q^{4s}, \frac{b^b}{c^b}, q^{-4s}; q^2\right)_M}{\left(-\frac{b^b}{qr_0 b^a b^c} q^{-4s+2}, \frac{b^b}{c^b} q^{4s+2}, q^2; q^2\right)_M} \frac{(1 - \frac{b^b}{c^b} q^{4M})}{\left(\frac{q^{-1} r_0 b^a b^c}{b^b}\right)^M (1 - \frac{b^b}{c^b})} q^{M(M-1)}. $$

In order to show (3.34), it remains to prove that the r.h.s. of (3.34) equals the r.h.s. of (3.36). This is equivalent to the equality

$$ (3.41) \qquad \sum_{M=0}^{2s} \frac{\left(\frac{-qr_0 b^a b^c}{c^b} q^{4s}, \frac{b^b}{c^b}, q^{-4s}; q^2\right)_M}{\left(-\frac{b^b}{qr_0 b^a b^c} q^{-4s+2}, \frac{b^b}{c^b} q^{4s+2}, q^2; q^2\right)_M} \frac{(1 - \frac{b^b}{c^b} q^{4M})}{\left(1 - \frac{b^b}{c^b}\right) \left(\frac{q^{-1} r_0 b^a b^c}{b^b}\right)^M} q^{M(M-1)} = \frac{(q^2 r_0^2 (b^a)^2; q^2)_{2s}}{\left(-\frac{q^{1-4s} b^b}{r_0 b^a b^c}; q^2\right)_{2s}}. $$

This latter relation is actually a limiting case and specialization of the orthogonality relation satisfied by the $q$-Racah polynomials, that is obtained as follows. Set $n = m = 0$ and substitute $q \to q^2$ in [KS96, eq. (3.2.2)] and consider its limit $\beta \to \infty$. Then, substitute

$$ (3.42) \qquad N \to 2s, \quad x \to M, \quad \alpha \to q^{-4s-2}, \quad \gamma \to -\frac{b^a b^c}{c^b} r_0 q^{4s-1}, \quad \delta \to -\frac{b^b}{r_0 b^a b^c} q^{-4s-1}, $$

in the resulting equation, we find (3.41). $\square$

In Appendix B, an explicit example of a Leonard triple is given for $s = 1/2$.

**3.2. Computation of scalar products.** The representation theory of Leonard triples of $q$-Racah type is now applied to the computation of scalar products of any eigenvector of $A$ with certain off-shell Bethe states. Recall the definition of the off-shell Bethe states (2.27), (2.28) and the relation between the parameter $\zeta$ and the parameters $b^{\diamond}, c^{\diamond}$ given by (3.8). For convenience, let us introduce the notation:

$$ (3.43) \qquad U_i = \frac{q u_i^2 + q^{-1} u_i^{-2}}{q + q^{-1}}. $$

**Lemma 3.7.** For $\beta = 0$, one has:

$$ (3.44) \qquad \pi(B^\epsilon(\bar{u}, m, M)) = \mathcal{G}_M^\epsilon(\bar{u}) \prod_{i=1}^M \left( U_i - \frac{r_0}{q + q^{-1}} A^\diamond \right) $$

with

$$ (3.45) \qquad \mathcal{G}_M^\epsilon(\bar{u}) = \left( \epsilon(q + q^{-1}) q^{-\epsilon(M+1)} b^{\frac{1-\epsilon}{2}} c^{*\frac{1+\epsilon}{2}} \right)^M \left( \prod_{i=1}^M b(u_i^2) u_i^{-\epsilon} \right) $$

*Proof.* For $\beta = 0$, (2.11) simplifies to

$$ (3.46) \qquad \pi(\mathcal{B}^\epsilon(u, m)) = b(u^2)u^{-\epsilon} \frac{q^{-\epsilon(m+2)}}{\alpha} \left( A^\diamond - \frac{(q+q^{-1})U}{r_0} \right) $$

with the identification (3.3). Then, the image of the string of dynamical operators $\pi(B^\epsilon(\bar{u}, m, M))$ with (2.26) is reduced, and depends on $\alpha$ and $m$. Now, consider the constraint on $\alpha$ according to the definition of the reference state. Choosing $\chi = -(q^2 - q^{-2})/r_0$ in (2.27) (resp. (2.28)) one has $\alpha = q^m/r_0b$ for $\epsilon = -$ (resp. $\alpha = -q^{-m}/r_0c^*$ for $\epsilon = +$). In each case $\epsilon = \pm$, one replaces the corresponding value of $\alpha$ in the reduced expression of $\pi(B^\epsilon(\bar{u}, m, M))$. After simplifications, one gets (3.44) $\square$

The off-shell Bethe states (2.27), (2.28) are now expressed in the eigenbasis of $A$.
**Lemma 3.8.**

$$ (3.47) \quad |\Psi_{-}^{M''}(\bar{u}, m)\rangle = \mathcal{G}_{M''}^{-}(\bar{u}) \sum_{N=0}^{2s} \sum_{M'=0}^{2s} \frac{f_N}{f_0} (P^{\{\diamond, \cdot\}})^{-1}_{NM'} \prod_{i=1}^{M''} \left(U_i - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) |\theta_N\rangle $$

$$ (3.48) \quad |\Psi_{+}^{M''}(\bar{u}, m)\rangle = \mathcal{G}_{M''}^{+}(\bar{u}) \sum_{N=0}^{2s} \sum_{M'=0}^{2s} \frac{f_N}{\nu_0^{\{*, \diamond, \cdot\}} g_0 h_{M'}} (P^{\{\diamond, \cdot\}})^{-1}_{NM'} \prod_{i=1}^{M''} \left(U_i - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) |\theta_N\rangle . $$

*Proof.* To show (3.47), we proceed as follows. Consider the l.h.s. of (3.47) defined by (2.27) with (2.25). Using the first relation in (3.33) for $N = 0, M \to M'$ and $(c, a, b) \to (\diamond, \cdot, *)$, the reference state $|\theta_0\rangle$ decomposes into the eigenvectors $|\theta_{M'}^{\diamond}\rangle$ of $A^{\diamond}$. Note that $P_{M'0}^{\{\diamond, \cdot\}} = 1$. Apply (3.44) to this combination. Then, rewrite the resulting expression in terms of the eigenvectors of $A$ using the second relation in (3.33). It reduces to the r.h.s. of (3.47). The proof of (3.48) is done similarly, observing that $(P^{\{*, \diamond\}})^{-1}_{M'0} = \nu_0^{\{*, \diamond, \cdot\}}^{-1}$ is independent of $M'$. $\square$

We are now ready to give the main result of this paper, the normalized scalar products of on-shell versus off-shell Bethe states, using the properties of the Leonard triple and related eigenbases. Recall (3.45), (3.34), (3.18), (3.28), (3.1), (3.29).

**Theorem 1.**

$$ (3.49) \quad \frac{\langle \theta_M | \Psi_{-}^{M''}(\bar{u}, m) \rangle}{\langle \theta_M | \theta_M \rangle} = \mathcal{G}_{M''}^{-}(\bar{u}) \frac{f_M}{f_0} \sum_{M'=0}^{2s} (P^{\{\diamond, \cdot\}})^{-1}_{MM'} \prod_{i=1}^{M''} \left( U_i - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right), $$

$$ (3.50) \quad \frac{\langle \theta_M | \Psi_{+}^{M''}(\bar{u}, m) \rangle}{\langle \theta_M | \theta_M \rangle} = \mathcal{G}_{M''}^{+}(\bar{u}) \frac{f_M}{g_0 \nu_0^{\{*, \diamond, \cdot\}}} \sum_{M'=0}^{2s} \frac{1}{h_{M'}} (P^{\{\diamond, \cdot\}})^{-1}_{MM'} \prod_{i=1}^{M''} \left( U_i - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right). $$

*Proof.* To show (3.49), we apply $\langle \theta_M |$ to (3.47) and use (2.19). We show (3.50) similarly. $\square$

In view of the next section, let us display the results for one of the simplest case, $s = 1/2$ and $M''(= 2s) = 1$. For the two examples below, we denote:

$$ (3.51) \quad X_M^{\pm, off}(u_1) = \frac{\langle \theta_M | \Psi_\pm^1(u_1, m) \rangle}{\langle \theta_M | \theta_M \rangle}, \quad M = 0, 1. $$

Recall the notation (1.2) and the identification (3.8).

**Example 1.** Consider $s = 1/2$ and $\epsilon = -$. From (3.49) we obtain:

$$ (3.52) \quad X_0^{-, off}(u_1) = \frac{q^3 r_0^2 b^2}{b(q) (q^2 r_0^2 b^2; q^2)_1} u_1 b(u_1^2) Y_-(u_1|u_1), $$

$$ (3.53) \quad X_1^{-, off}(u_1) = - \frac{qb(q)b\left(-\frac{qr_0b^\diamond b^*}{b^*}; q^2\right)_1 \left(-\frac{qr_0b^\diamond b^*}{b}; q^2\right)_1}{r_0 b^\diamond (q^2 r_0^2 b^2; q^2)_1} u_1 b(u_1^2), $$

with

$$ (3.54) \quad Y_-(x|x) = \frac{b(q)}{q^3 r_0^3 b^* b^\diamond} \left( \frac{qr_0 b^\diamond b^* (q^2 r_0^2 b^2 - 1)}{x^2} + q^3 x^2 r_0 b^\diamond b^* (q^2 r_0^2 b^2 - 1) + b^* (1 - q^4 r_0^2 b^2) \right. \\ \left. - (q^2 - 1) qr_0 b^\diamond (q^2 r_0^2 (b^*)^2 + 1) + q^2 r_0^2 (b^\diamond)^2 b^* (1 - q^4 r_0^2 b^2) \right). $$

**Example 2.** Consider $s = 1/2$ and $\epsilon = +$. From (3.50) we obtain:

$$ (3.55) \quad X_0^{+, off}(u_1) = \frac{1}{q^4 r_0^3 b^* b^\diamond} u_1^{-1} b(u_1^2) Z_1(u_1), $$

$$ (3.56) \quad X_1^{+, off}(u_1) = \frac{1}{q^4 r_0^3 b^* b^\diamond} u_1^{-1} b(u_1^2) Z_2(u_1), $$
where $Z_i(x)$ are Laurent polynomials in $x$ given by,

$$ (3.57) \qquad Z_1(x) = q^3 r_0 b^{\diamond} x^2 + \frac{qr_0 b^{\diamond}}{x^2} - 1 + q^2 r_0^2 b^{\diamond} (q(q^2 - 1) r_0 b b^* - b^{\diamond}), $$

$$ (3.58) \qquad Z_2(x) = q^3 r_0 b b^{\diamond} x^2 + \frac{qr_0 b b^{\diamond}}{x^2} + q(q^2 - 1) r_0 b^{\diamond} b^* - b(q^2 r_0^2 (b^{\diamond})^2 + 1). $$

**Remark 7.** As an example of application of Theorem 1 to an elementary integrable system, on one hand consider the Hamiltonian $H(\lambda)$ obtained from specializing [BaP19, eq. (5.2)] for $\kappa_{\pm} = 0$ and identifying $\lambda = \frac{\kappa^*}{\kappa}$:

$$ H(\lambda) = S_1^x S_2^x + S_1^y S_2^y + \Delta S_1^z S_2^z + \lambda(S_2^x S_3^x + S_2^y S_3^y + \Delta S_2^z S_3^z) \\ + \frac{1}{q-q^{-1}}(S_1^z + (\lambda-1)S_2^z - \lambda S_3^z), \qquad \Delta = \frac{q+q^{-1}}{2}, $$

where

$$ (3.59) \qquad S^x = \frac{1}{2} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad S^y = \frac{1}{2} \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad S^z = \frac{1}{2} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. $$

This Hamiltonian is the image in End($\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$) of $A + \lambda A^*$ via the specialization of the map [BaP19, eq. (2.14), (2.15)]. We refer the reader to [BaP19, Subsection 5.2.1] where the solution (spectrum and eigenstates) is discussed in details. Using the framework of Leonard pairs, the Bethe eigenstates of $H(\lambda)$ are given by (2.27) (or (2.28)) for $M = 2s$ where $\bar{u}$ satisfy the inhomogeneous Bethe equations [BaP19, eq. (3.75)] for $\epsilon = -$ (or $\epsilon = +$). On the other hand, consider the Hamiltonian $H' = H(0)$. The eigenstates are given by (2.33). Thus, the ratios in Theorem 1 for $\bar{u}$ satisfying [BaP19, eq. (3.75)] determine the overlap between the eigenstates of $H(\lambda)$ and $H'$. In the context of quantum phase transitions where $\lambda$ is viewed as the driving parameter (see e.g. [Gu08, Section III]), those ratios determine the fidelity when considering the ground states of each Hamiltonian.

**3.3. Relating inhomogeneous and homogeneous Bethe roots.** In Subsection 2.2.2, we have seen that the eigenvector $|\theta_N^*\rangle$ of $A^*$ can be either written as the on-shell Bethe state (2.34) associated with the Bethe equations of homogeneous type $E_+^N(u_j, \bar{u}_j) = 0$ with (2.31), $j = 0, ..., N$, or as the on-shell Bethe state (2.42) associated with the Bethe equations of inhomogeneous type $E_-^N(u_j, \bar{u}_j) = 0$, $j = 0, ..., 2s$ with (2.36), (2.40). In this section, a set of equations relating the inhomogenous and homogeneous Bethe roots associated with the eigenvalue $\theta_N^*$, $N = 0, ..., 2s$, is obtained.

Let us denote the symmetric Bethe roots for the homogeneous case and inhomogeneous case as follows:

$$ (3.60) \qquad U_j^{(i)} \quad \text{for} \quad \bar{u} = S_-^{*N(i)}, \quad U_j^{(h)} \quad \text{for} \quad \bar{u} = S_+^{*N(h)} \quad \text{with} \quad (3.43). $$

We also introduce the notation $\bar{U}_j = \{U_1, U_2, ..., U_{j-1}, U_{i+1}, ..., U_N\}$. From [BaP19, Proposition 3.5], each system of Bethe equations is equivalent to a system of polynomial equations:

$$ (3.61) \qquad E_+^N(u_j, \bar{u}_j) = 0 \quad \Leftrightarrow \quad P_+^N(U_j^{(h)}, \bar{U}_j^{(h)}) = 0 \quad \text{for} \quad j = 1, ..., N, $$

$$ (3.62) \qquad E_-(u_j, \bar{u}_j) = 0 \quad \Leftrightarrow \quad P_-^{2s}(U_j^{(i)}, \bar{U}_j^{(i)}) = 0 \quad \text{for} \quad j = 1, ..., 2s. $$

respectively. Here, $P_+^N(U_j, \bar{U}_j)$ (resp. $P_-^{2s}(U_j, \bar{U}_j)$) has maximal degree $N$ (resp. degree $2s+1$) in the variables $U_1, ..., U_N$ (resp. $U_1, ..., U_{2s}$). For more details, see [BaP19, Subsection 3.5.1]. Following [BaP19, Conjectures 2,3] based on numerical supporting evidences:

**Hypothesis 1.** The system of polynomial equations (3.61) admits a unique admissible solution $\bar{U}^{(h)} = \{U_1^{(h)}, U_2^{(h)}, ..., U_N^{(h)}\}$ up to permutation.

**Hypothesis 2.** The system of polynomial equations (3.62) admits $2s+1$ distinct admissible solution $\bar{U}^{(i)} = \{U_1^{(i)}, U_2^{(i)}, ..., U_{2s}^{(i)}\}$ up to permutation.

In order to relate the two systems of Bethe roots, we specialize Theorem 1. The following lemmas give two relations between $q$-Racah polynomials associated with different parameter sequences, that involve either homogeneous Bethe roots or inhomogenous Bethe roots.
**Lemma 3.9.** Assume Hypothesis 1. The q-Racah polynomial $R_M^{\{\diamond, *\}}(\theta_N^*)$ admits the following decomposition of homogeneous type in terms of $R_{M'}^{\{\diamond, *\}}(\theta_M)$:

$$ (3.63) \quad R_M^{\{\diamond, *\}}(\theta_N^*) = \frac{f_M \sum_{M'=0}^{2s} \prod_{j=1}^N \left( U_j^{(h)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) \frac{1}{h_{M'}} k_{M'}^{\{\diamond, *\}} R_{M'}^{\{\diamond, *\}}(\theta_M)}{f_0 \sum_{M'=0}^{2s} \prod_{j=1}^N \left( U_j^{(h)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) \frac{1}{h_{M'}} k_{M'}^{\{\diamond, *\}}} \quad \text{with } \bar{u} = S_+^{*N(h)}. $$

*Proof.* Firstly, we compute the normalized scalar product $\langle \theta_M | \theta_N^* \rangle / \langle \theta_M | \theta_M \rangle$. To do that, we multiply both sides of the equality in (3.50) by $\mathcal{N}_N^*(\bar{u})$ with (2.34) and specialize the result to $\bar{u} = S_+^{*N(h)}$. After simplifications using (3.29), it follows:

$$ \frac{\langle \theta_M | \theta_N^* \rangle}{\langle \theta_M | \theta_M \rangle} = \mathcal{N}_N^*(\bar{u}) \mathcal{G}_N^+(\bar{u}) \nu_0^{\{*, \diamond, *\}}^{-1} \frac{f_M}{g_0} \sum_{M'=0}^{2s} \prod_{j=1}^N \left( U_j^{(h)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) \frac{1}{h_{M'}} k_{M'}^{\{\diamond, *, *\}} R_{M'}^{\{\diamond, *, *\}}(\theta_M) \quad \text{for } \bar{u} = S_+^{*N(h)}. $$

Combining above expression with the fact that the q-Racah polynomials are given by (see [T03, Theorem 14.6 and 15.6]):

$$ (3.64) \quad R_M^{\{\diamond, *\}}(\theta_N^*) = \frac{\langle \theta_M | \theta_N^* \rangle}{\langle \theta_0 | \theta_N^* \rangle} \frac{\langle \theta_0 | \theta_0 \rangle}{\langle \theta_M | \theta_M \rangle}, $$

we get (3.63).

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

Similar arguments are used to derive the following lemma, thus we skip the proof.

**Lemma 3.10.** Assume Hypothesis 2. The q-Racah polynomial $R_M^{\{\diamond, *\}}(\theta_N^*)$ admits the following decomposition of inhomogeneous type in terms of $R_{M'}^{\{\diamond, *\}}(\theta_M)$:

$$ (3.65) \quad R_M^{\{\diamond, *\}}(\theta_N^*) = \frac{f_M \sum_{M'=0}^{2s} \prod_{j=1}^2 \left( U_j^{(i)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) k_{M'}^{\{\diamond, *, *\}} R_{M'}^{\{\diamond, *, *\}}(\theta_M)}{f_0 \sum_{M'=0}^{2s} \prod_{j=1}^2 \left( U_j^{(i)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) k_{M'}^{\{\diamond, *, *\}}} \quad \text{with } \bar{u} = S_-^{*N(i)}. $$

Combining above lemmas, it follows:

**Proposition 3.4.** Assume Hypothesis 1 and 2. Let $U_j^{(h)}$, $j = 1, ..., N$ be the symmetric Bethe roots associated with the eigenvalue $\theta_N^*$ given in (2.22) and homogeneous Bethe equations $E_+^N(u_j, \bar{u}_j) = 0$ with (2.31). The inhomogenous symmetric Bethe roots $U_j^{(i)}$, $j = 0, ..., 2s$ associated with $\theta_N^*$ satisfy the relations:

$$ (3.66) \quad \frac{\prod_{j=1}^{2s} \left( U_j^{(i)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right)}{\sum_{M'=0}^{2s} \prod_{j=1}^{2s} \left( U_j^{(i)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) k_{M'}^{\{\diamond, *, *\}}} = \frac{\frac{1}{h_{M'}} \prod_{j=1}^{N} \left( U_j^{(h)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right)}{\sum_{M'=0}^{2s} \frac{1}{h_{M'}} \prod_{j=1}^{N} \left( U_j^{(h)} - \frac{r_0}{q+q^{-1}} \theta_{M'}^{\diamond} \right) k_{M'}^{\{\diamond, *, *\}}}, \quad \text{for } M' = 0, 1, ..., 2s. $$

*Proof.* The difference between the r.h.s. of (3.65) and the r.h.s. of (3.63) is vanishing. Using the fact that the q-Racah polynomials $R_{M'}^{\{\diamond, *\}}(\theta_M)$ are linearly independent, one gets (3.66).

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

**Example 3.** Consider $s = 1/2$. From (3.66), we obtain$^5$:

$$ (3.67) \quad \text{For } N=0: \quad U_1^{(i)} = \left( \frac{r_0}{q+q^{-1}} \right) \left( \frac{h_0 \theta_0^{\diamond} - h_1 \theta_1^{\diamond}}{h_0 - h_1} \right), $$

$$ (3.68) \quad \text{For } N=1: \quad U_1^{(i)} = \left( \frac{r_0}{q+q^{-1}} \right) \left( \frac{\left( \frac{r_0}{q+q^{-1}} \right) (h_1 - h_0) \theta_0^{\diamond} \theta_1^{\diamond} + (h_0 \theta_0^{\diamond} - h_1 \theta_1^{\diamond}) U_1^{(h)}}{\left( \frac{r_0}{q+q^{-1}} \right) (h_1 \theta_0^{\diamond} - h_0 \theta_1^{\diamond}) + (h_0 - h_1) U_1^{(h)}} \right). $$

From (3.61) with (2.31), we find a unique admissible solution of (2.31,3.61) given by

$$ (3.69) \quad U_1^{(h)} = \frac{\mathbf{b}(q^2 \mathbf{b}^* - q^{-2} \mathbf{c}^*)(1 + q^2 r_0^2 (\mathbf{b}^*)^2) + (q - q^{-1}) r_0^{-1} \mathbf{b}^* (1 + q^2 r_0^2 \mathbf{b}^2)}{(q + q^{-1}) r_0 \mathbf{b}^* (q^2 \mathbf{b}^* - \mathbf{c}^*)}, $$

in agreement with Hypothesis 1. Independently, it is checked that $E_-(u_i, \bar{u}_i) = 0$ with (2.36) admits only two admissible solutions in agreement with Hypothesis 2, and that they coincide with (3.67) and (3.68) with (3.69).

<footnote>$^5$Note that for $N = 0$ there are no homogeneous Bethe roots in (3.66).</footnote>
For higher values of $s = 1, 3/2$, we proceed as follows. Firstly, admissible homogeneous symmetric Bethe roots $U_j^{(h)}$, $j = 1, ..., N$, are determined numerically from (3.61), and Hypothesis 1 is checked. Replacing the numerical values in (3.66), for each $N = 0, ..., 2s$, the inhomogeneous symmetric Bethe roots $U_j^{(i)}$, $j = 1, ..., 2s$, are determined numerically. Independently, the admissible solutions of $E_-(u_i, \bar{u}_i) = 0$ with (2.36) are computed, and Hypothesis 2 is checked. Those admissible solutions are found to agree with the solutions previously computed from (3.66).

**Example 4.** Consider $s = 1$ and the numerical values $q = 3, r_0 = 1, b = 5, b^* = 7$ and $b^\diamond = 1/2$. In Table 1 the numerical solutions of the system (3.66) for $N = 0, 1, 2$ are given.

<table><caption>TABLE 1. Homogeneous and inhomogeneous Bethe roots that solve (3.61), (3.66) for the scalars $q = 3, r_0 = 1, b = 5, b^* = 7$ and $b^\diamond = 1/2$.</caption><thead><tr><th></th><th>{U<sub>1</sub><sup>(h)</sup>, U<sub>2</sub><sup>(h)</sup>}</th><th>{U<sub>1</sub><sup>(i)</sup>, U<sub>2</sub><sup>(i)</sup>}</th></tr></thead><tbody><tr><td>N = 0</td><td>-</td><td>{3.50405, 11.9071}</td></tr><tr><td>N = 1</td><td>{18.5087, -}</td><td>{3.208, 12.0789}</td></tr><tr><td>N = 2</td><td>{2.72742, 12.0749}</td><td>{2.01305, 12.1539}</td></tr></tbody></table>

### 4. SCALAR PRODUCTS FOR $M'' = 2s$ AND BELLIARD-SLAVNOV LINEAR SYSTEM

In this section, we consider the normalized scalar products of on-shell and off-shell Bethe states (3.49), (3.50) at the special value $M'' = 2s$ from a different perspective. It is shown that those scalar products satisfy systems of homogeneous linear equations of the form introduced and studied in [BelS19]. It implies the existence of a determinant formula for the solutions of Theorem 1 at $M'' = 2s$, and for $q$-Racah polynomials when those solutions are specialized to Bethe roots of inhomogeneous type.

**4.1. Belliard-Slavnov linear system associated with a Leonard pair.** We apply the method recently proposed in [BelS19]. Below, we denote:

$$ (4.1) \qquad \bar{\mathcal{Y}} = \{y_1, y_2, \dots, y_{2s}, y_{2s+1}\}, \quad \bar{\mathcal{Y}}_j = \bar{\mathcal{Y}} \setminus y_j \quad \text{and} \quad \bar{\mathcal{Y}}_{j,k} = \bar{\mathcal{Y}}_j \setminus y_k. $$

Introduce the convenient notations:

$$ (4.2) \qquad g(u, \bar{u}) = \prod_{u_i \in \bar{u}} g(u, u_i) \quad \text{with} \quad g(u, u_i) = (b(u/u_i)b(quu_i))^{-1} $$

and

$$ (4.3) \qquad Y_{\epsilon}(u|\bar{u}) = \frac{u^{2\epsilon+1} \Lambda_{1}^{\epsilon}(u)}{b(u^2)b(qu^2)} \prod_{u_i \in \bar{u}} b(u/qu_i)b(uu_i) + \frac{q^{-\epsilon-1} u^{-1} \Lambda_{2}^{\epsilon}(u)}{b(u^2)b(q^2u^2)} \prod_{u_i \in \bar{u}} b(qu/u_i)b(q^2uu_i) \\ - \delta_{\epsilon,+} \nu_{+} \prod_{j=0}^{2s} b(q^{1/2+j-s}u\zeta)b(q^{1/2+j-s}u\zeta^{-1}), $$

where $\epsilon = \pm$, $\delta_{+,+} = 1$, $\delta_{-,+} = 0$ and $\nu_+$ is given in (2.37).

Given a Leonard pair $(A, A^*)$, different systems of linear equations can be derived, either starting from $A$ or from $A^*$. In order to relate with the results of Theorem 1, here we focus on a system derived from studying the operator $A$ only. As a preliminary, the action of $A$ on two different types of off-shell Bethe states is computed. Recall the structure constants $\eta, \eta^*$ given by (3.4) with (3.2) and the identification (3.8).

**Lemma 4.1.**

$$ (4.4) \qquad A|\Psi_{\epsilon}^{2s}(\bar{\mathcal{Y}}_j, m)\rangle = \sum_{k=1}^{2s+1} L_{jk}^{\epsilon}|\Psi_{\epsilon}^{2s}(\bar{\mathcal{Y}}_k, m)\rangle, \quad j = 1, \dots, 2s+1, $$

with

$$ (4.5) \qquad L_{jk}^{\epsilon} = \left(\frac{y_j}{y_k}\right)^{\epsilon} \frac{b(y_k^2)}{b(y_j^2)} g(y_k, \bar{\mathcal{Y}}_k) Y_{\epsilon}(y_k | \bar{\mathcal{Y}}_j) + \frac{(q+q^{-1})^2}{\rho b(y_j^2) b(q^2 y_j^2)} \left(\eta^* + \eta \frac{q y_j^2 + q^{-1} y_j^{-2}}{q + q^{-1}}\right) \delta_{jk}. $$
Proof. Firstly, we consider the action of A on off-shell Bethe states. It is obtained by specializing the equations below [BaP19, eq. (3.76)] for $\kappa = 1, \kappa^* = 0$. One gets the off-shell equation:

$$ (4.6) \qquad \left( A - \frac{(q+q^{-1})^2}{\rho b(u^2)b(q^2u^2)} \left( \eta^* + \eta \frac{qu^2+q^{-1}u^{-2}}{q+q^{-1}} \right) \right) |\Psi_{\epsilon}^{2s}(\bar{u}, m)\rangle \\ = \Lambda_{\epsilon}^{2s}(u, \bar{u}) |\Psi_{\epsilon}^{2s}(\bar{u}, m)\rangle + \sum_{k=1}^{2s} \frac{(q-q^{-1})u^{\epsilon}u_{k}^{(\epsilon+1)/2}\bar{E}_{\epsilon}(u_{k}, \bar{u}_{k})}{(u^2-u^{-2})b(uu_{k}^{-1})b(quu_{k})} |\Psi_{\epsilon}^{2s}(\{u, \bar{u}_{k}\}, m)\rangle $$

where

$$ \Lambda_{\epsilon}^{2s}(u, \bar{u}) = \frac{u^{2\epsilon+1}\Lambda_{1}^{\epsilon}(u)}{b(u^2)b(qu^2)} \prod_{j=1}^{2s} f(u, u_j) + \frac{q^{-\epsilon-1}u^{-1}\Lambda_{2}^{\epsilon}(u)}{b(u^2)b(q^2u^2)} \prod_{j=1}^{2s} h(u, u_j) - \delta_{\epsilon,+\nu_{+}} \frac{\prod_{j=0}^{2s} b(q^{1/2+j-s}\zeta u) b(q^{1/2+j-s}\zeta^{-1}u)}{\prod_{i=1}^{2s} b(uu_i^{-1}) b(quu_i)} $$

and $\bar{E}_-(u_k, \bar{u}_k) = -E_-(u_k, \bar{u}_k)$ with (2.31), $\bar{E}_+(u_k, \bar{u}_k) = E_+(u_k, \bar{u}_k)$, with (2.36).

Secondly using the expressions (4.2), (4.3), one rewrites:

$$ A|\Psi_{\epsilon}^{2s}(\bar{u}, m)\rangle = \left( g(u, \bar{u})Y_{\epsilon}(u|\bar{u}) + \frac{(q+q^{-1})^2}{\rho b(u^2)b(q^2u^2)} \left( \eta^* + \eta \frac{qu^2+q^{-1}u^{-2}}{q+q^{-1}} \right) \right) |\Psi_{\epsilon}^{2s}(\bar{u}, m)\rangle \\ - \sum_{u_k \in \bar{u}} \left( \frac{u}{u_k} \right)^{\epsilon} \frac{b(u_k^2)}{b(u^2)} g(u, u_k) g(u_k, \bar{u}_k) Y(u_k|\bar{u}) |\Psi_{\epsilon}^{2s}(\{u, \bar{u}_k\}, m)\rangle. $$

Set $u = y_j, \bar{u} = \bar{Y}_j$. It follows:

$$ A|\Psi_{\epsilon}^{2s}(\bar{Y}_j, m)\rangle = \left( g(y_j, \bar{Y}_j)Y_{\epsilon}(y_j|\bar{Y}_j) + \frac{(q+q^{-1})^2}{\rho b(y_j^2)b(q^2y_j^2)} \left( \eta^* + \eta \frac{qy_j^2+q^{-1}y_j^{-2}}{q+q^{-1}} \right) \right) |\Psi_{\epsilon}^{2s}(\bar{Y}_j, m)\rangle \\ - \sum_{k=1}^{2s} \left( \frac{y_j}{y_k} \right)^{\epsilon} \frac{b(y_k^2)}{b(y_j^2)} \underbrace{g(y_j, y_k)g(y_k, \bar{Y}_{j,k})}_{=g(y_k, \bar{Y}_k)} Y_{\epsilon}(y_k|\bar{Y}_j) |\Psi_{\epsilon}^{2s}(\{y_j, \bar{Y}_{j,k}\}, m)\rangle, $$

which gives (4.5).

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

Then, independently of Theorem 1 we have:

**Proposition 4.1.** For any $M = 0, 1, \dots, 2s$, the normalized scalar product

$$ (4.7) \qquad X_{M}^{\epsilon, \text{off}}(\bar{u}) = \frac{\langle \theta_{M} | \Psi_{\epsilon}^{2s}(\bar{u}, m) \rangle}{\langle \theta_{M} | \theta_{M} \rangle} $$

solves the Belliard-Slavnov linear system

$$ (4.8) \qquad \sum_{k=1}^{2s+1} \mathcal{M}_{jk}^{\epsilon,M} X_{M}^{\epsilon,\text{off}}(\bar{Y}_k) = 0 \quad \text{with} \quad \mathcal{M}_{jk}^{\epsilon,M} = L_{jk}^{\epsilon} - \delta_{jk}\theta_M \quad \text{for } j=1,2,\dots,2s+1 $$

with (4.5) and $\bar{Y}$ are distinct arbitrary complex numbers.

Proof. From (4.4), one has:

$$ (4.9) \qquad (A - \theta_M \mathbb{I}) |\Psi_{\epsilon}^{2s}(\bar{Y}_j, m)\rangle = \sum_{k=1}^{2s+1} (L_{jk}^{\epsilon} - \theta_M \delta_{jk}) |\Psi_{\epsilon}^{2s}(\bar{Y}_k, m)\rangle \\ = \sum_{k=1}^{2s+1} \mathcal{M}_{jk}^{\epsilon,M} |\Psi_{\epsilon}^{2s}(\bar{Y}_k, m)\rangle. $$

Acting from the left with $\langle \theta_M |$ and using (2.17), (2.19), one gets (4.8).

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

Independently of Theorem 1, we may now ask about the existence of non-trivial solutions of the system of homogenous linear equations (4.8).

**Lemma 4.2.** Let $\mathcal{M}^{\epsilon,M}$, $\epsilon = \pm$, $M = 0, \dots, 2s$ denote the matrices with entries $\mathcal{M}_{jk}^{\epsilon,M}$, $j, k = 1, \dots, 2s+1$, given in (4.8).

$$ (4.10) \qquad \mathrm{rank}(\mathcal{M}^{\epsilon,M}) = 2s. $$
Proof. The vectors $|\Psi_{\epsilon}^{2s}(\bar{Y}_k, m)\rangle$, $k = 1, ..., 2s+1$, given by (3.47), (3.48) with $M'' = 2s$ are linearly independent. By (4.9), the matrix with entries $M_{jk}^{\epsilon,M}$ represents the linear map $(\pi(\mathbf{A}) - \theta_M I)$ with respect to the $2s+1$ vectors $|\Psi_{\epsilon}^{2s}(\bar{Y}_k, m)\rangle$. By linear algebra, $\mathrm{rank}(\pi(\mathbf{A}) - \theta_M I) = \mathrm{dim}\mathcal{V} - \mathrm{dim}(\mathrm{kernel}(\pi(\mathbf{A}) - \theta_M I))$. By definition, we have $\mathrm{dim}\mathcal{V} = 2s + 1$ and one knows that the dimension of the each eigenspace of $\pi(\mathbf{A})$ is 1. So, $\mathrm{rank}(M^{\epsilon,M}) = \mathrm{rank}(\pi(\mathbf{A}) - \theta_M I) = 2s$. $\square$

The above lemma implies that the linear system (4.8) admits non-trivial solutions. Results of the previous section provide explicit solutions to (4.8). The following is an immediate consequence of Theorem 1 and Proposition 4.1.

**Theorem 2.** The normalized scalar product of on-shell and off-shell Bethe states given by (3.49) (resp. (3.50)) at $M'' = 2s$ solves the linear system (4.8) at $\epsilon = -$ (resp. $\epsilon = +$).

For instance, recall Examples 1, 2. Using the expressions (3.52), (3.53), (3.55), (3.56), it is readily checked that the linear system (4.7)

$$ \begin{aligned} & M_{11}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(y_2) + M_{12}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(y_1) = 0, \\ & M_{21}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(y_2) + M_{22}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(y_1) = 0, \end{aligned} $$

is satisfied for $M = 0, 1, \epsilon = \pm$. For higher values of $s = 1, 3/2$, using Mathematica it can be independently checked that the expressions (3.49) and (3.50) for $M'' = 2s$ solve the linear system (4.7) for $\epsilon = -$ and $\epsilon = +$, respectively.

**4.2. Determinant formula and q-Racah polynomials.** The fact that the normalized scalar products (4.7) solve the systems of homogeneous linear equations (4.8) for $\epsilon = \pm$ implies the existence of a determinant formula for (3.49), (3.50) at $M'' = 2s$. Let $T$ be a $2s+1 \times 2s+1$ matrix, and denote by $T_{[k,l]}$ the submatrix of $T$ obtained by deleting the $k$-th row and $l$-column.

**Proposition 4.2.** The normalized scalar product (4.7) admits a determinant formula given by:

$$ (4.11) \qquad X_M^{\epsilon,\mathrm{off}}(\bar{Y}_k) = \psi_M^{\epsilon}(-1)^k \det \tilde{M}_{[2s+1,k]}^{\epsilon,M}, $$

where $\tilde{M}^{\epsilon,M} = W^{\epsilon,M}M^{\epsilon,M}$ for some non-degenerate $2s+1 \times 2s+1$ matrix $W^{\epsilon,M}$ such that:

$$ (i) \quad (W^{\epsilon,M}M^{\epsilon,M})_{2s+1,k} = 0 \quad \text{for all } k = 1, \dots, 2s+1, $$

$$ (ii) \quad \det \tilde{M}_{[2s+1,k]}^{\epsilon,M} \quad \text{is independent of } y_k, $$

and $\psi_M^{\epsilon}$ some function of $M, \epsilon$.

**Proof.** By Lemma 4.2, given $M, \epsilon$ fixed the solutions of the system of homogeneous linear equations (4.8) can be written in terms of minors of a certain matrix denoted $\tilde{M}^{\epsilon,M}$ - up to an overall factor $\psi_M^{\epsilon}$ - using Cramer's rule. The matrix $\tilde{M}^{\epsilon,M}$ is determined as follows. Firstly, we are looking for a $2s+1 \times 2s+1$ matrix $W^{\epsilon,M}$ such that conditions (i) hold, which turns the system (4.8) into the new system:

$$ (4.12) \qquad \sum_{k=1}^{2s} \tilde{M}_{jk}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(\bar{Y}_k) = - \tilde{M}_{j2s+1}^{\epsilon,M} X_M^{\epsilon,\mathrm{off}}(\bar{Y}_{2s+1}) \quad \text{for } j = 1, 2, \dots, 2s. $$

Looking for non-trivial solutions using Cramer's rule, we assume $W^{\epsilon,M}$ is non-degenerate such that $\det \tilde{M}_{[2s+1,k]}^{\epsilon,M} (= \det(W^{\epsilon,M}M^{\epsilon,M})_{[2s+1,k]}) \neq 0$. Moreover, because the l.h.s. of (4.11) is independent of $y_k$ we are looking for a matrix $W^{\epsilon,M}$ such that conditions (ii) hold. Applying Cramer's rule, the solutions of (4.12) are written in terms of $X_M^{\epsilon,\mathrm{off}}(\bar{Y}_{2s+1})$:

$$ (4.13) \qquad X_M^{\epsilon,\mathrm{off}}(\bar{Y}_k) = (-1)^{k+2s+1} \frac{\det \tilde{M}_{[2s+1,k]}^{\epsilon,M}}{\det \tilde{M}_{[2s+1,2s+1]}^{\epsilon,M}} X_M^{\epsilon,\mathrm{off}}(\bar{Y}_{2s+1}), $$

which leads to (4.11) for some overall factor $\psi_M^{\epsilon}$ that is independent of the variables $y_j, j = 1, ..., 2s+1$. Finally, this factor is determined by comparing (3.49), (3.50) for $M'' = 2s$ at the limit $\bar{u} \to \infty$ to the r.h.s. of (4.11) with (4.7) at the limit $\bar{Y}_k \to \infty$. $\square$
By specializing the parameters $\tilde{Y}_k$ in (4.11) to certain Bethe roots, the off-shell Bethe states entering in the ratio (4.7) reduce to on-shell Bethe states that are eigenvectors of the operator $A$. This leads to a determinant formula for $q$-Racah polynomials. Recall the parameters' correspondence (3.8).

**Corollary 4.1.** The $q$-Racah polynomial $R_M^{\{\cdot,\cdot\}}(\theta_N^*)$ with (3.30) admits a determinant formula given by:

$$ (4.14) \qquad R_M^{\{\cdot,\cdot\}}(\theta_N^*) = \frac{\bar{\psi}_M \det \tilde{\mathcal{M}}_{[2s+1,2s+1]}^{-,M}}{\psi_0 \det \tilde{\mathcal{M}}_{[2s+1,2s+1]}^{-,0}} \quad \text{for } \bar{\mathcal{Y}}_{2s+1} = S_-^{*N(i)}. $$

*Proof.* Multiply (4.7) for $\epsilon = -$ by the second normalization factor in (2.43) with $\bar{u} = \bar{\mathcal{Y}}_{2s+1}$. Let us specialize $\bar{\mathcal{Y}}_{2s+1} = S_-^{*N(i)}$ in the resulting ratio. Using (2.42), it follows:

$$ (4.15) \qquad \frac{\langle \theta_M | \theta_N^* \rangle}{\langle \theta_M | \theta_M \rangle} = \mathcal{N}_N^*(\bar{\mathcal{Y}}_{2s+1}) \frac{\langle \theta_M | \Psi_{-}^{2s} (\bar{\mathcal{Y}}_{2s+1}, m) \rangle}{\langle \theta_M | \theta_M \rangle} \quad \text{for } \bar{\mathcal{Y}}_{2s+1} = S_-^{*N(i)}. $$

Then, we use (4.11) with (4.7) to express the l.h.s. of (4.15) in terms of a determinant. Finally, we use (3.64) to get (4.14). $\square$

**4.3. Examples.** To illustrate Proposition 4.2 and Corollary 4.1, we consider some examples below.

* The case $s = 1/2$, $\epsilon = -$: Consider $M = 0$. Firstly, one shows the identity

$$ g(y_1, y_2)Y_-(y_1|y_2) + \frac{(q+q^{-1})^2}{\rho b(y_1^2)b(q^2y_1^2)} \left( \eta^* + \eta \frac{qy_1^2 + q^{-1}y_1^{-2}}{q + q^{-1}} \right) - \theta_0 = g(y_1, y_2)Y_-(y_1|y_1), $$

from which we find the entries of $\tilde{\mathcal{M}}^{-,0}$ are given by

$$ \tilde{\mathcal{M}}_{11}^{-,0} = \frac{g(y_1, y_2)Y_-(y_1|y_1)}{y_2b(y_2^2)} \left( y_2b(y_2^2)\mathcal{W}_{11}^{-,0}(y_1, y_2) + y_1b(y_1^2)\mathcal{W}_{12}^{-,0}(y_1, y_2) \right), $$

$$ \tilde{\mathcal{M}}_{12}^{-,0} = \frac{g(y_2, y_1)Y_-(y_2|y_2)}{y_1b(y_1^2)} \left( y_2b(y_2^2)\mathcal{W}_{11}^{-,0}(y_1, y_2) + y_1b(y_1^2)\mathcal{W}_{12}^{-,0}(y_1, y_2) \right), $$

$$ \tilde{\mathcal{M}}_{21}^{-,0} = \frac{g(y_1, y_2)Y_-(y_1|y_1)}{y_2b(y_2^2)} \left( y_2b(y_2^2)\mathcal{W}_{21}^{-,0}(y_1, y_2) + y_1b(y_1^2)\mathcal{W}_{22}^{-,0}(y_1, y_2) \right), $$

$$ \tilde{\mathcal{M}}_{22}^{-,0} = \frac{g(y_2, y_1)Y_-(y_2|y_2)}{y_1b(y_1^2)} \left( y_2b(y_2^2)\mathcal{W}_{21}^{-,0}(y_1, y_2) + y_1b(y_1^2)\mathcal{W}_{22}^{-,0}(y_1, y_2) \right). $$

Then, conditions (i) imply

$$ (4.16) \qquad \mathcal{W}_{21}^{-,0}(y_1, y_2) = -\frac{y_1b(y_1^2)}{y_2b(y_2^2)}\mathcal{W}_{22}^{-,0}(y_1, y_2) $$

with $\mathcal{W}_{22}^{-,0}(y_1, y_2)$ arbitrary. To satisfy conditions (ii), one imposes $\frac{\partial \tilde{\mathcal{M}}_{12}^{-,0}}{\partial y_1} = \frac{\partial \tilde{\mathcal{M}}_{11}^{-,0}}{\partial y_2} = 0$. The associated differential equations, omitted here for brevity, are solved by

$$ (4.17) \qquad \mathcal{W}_{12}^{-,0}(y_1, y_2) = \frac{y_2b(y_2^2)}{y_1b(y_1^2)} \left( \frac{qy_1b(y_1^2)}{g(y_1, y_2)} C^{-,0} - \mathcal{W}_{11}^{-,0}(y_1, y_2) \right), $$

where $C^{-,0}$ is a free scalar and $\mathcal{W}_{11}^{-,0}(y_1, y_2)$ remains undetermined. Using (4.16),(4.17), one finds that the nonzero entries of $\tilde{\mathcal{M}}^{-,0}$ are given by

$$ (4.18) \qquad \tilde{\mathcal{M}}_{11}^{-,0} = qy_1b(y_1^2)Y_-(y_1|y_1)C^{-,0}, \quad \tilde{\mathcal{M}}_{12}^{-,0} = -qy_2b(y_2^2)Y_-(y_2|y_2)C^{-,0}. $$

Comparing (3.52) to (4.11) where $\det \tilde{\mathcal{M}}_{[2,1]}^{-,0} = \tilde{\mathcal{M}}_{12}^{-,0}$ we finally identify $\psi_0^-$:

$$ (4.19) \qquad \psi_0^- = \frac{q^2 r_0^2 b^2}{b(q) (q^2 r_0^2 b^2; q^2) {}_1C^{-,0}}. $$

Now, consider $M = 1$. Firstly, one shows the identity

$$ g(y_1, y_2)Y_-(y_1|y_2) + \frac{(q+q^{-1})^2}{\rho b(y_1^2)b(q^2y_1^2)} \left( \eta^* + \eta \frac{qy_1^2 + q^{-1}y_1^{-2}}{q + q^{-1}} \right) - \theta_1 = g(y_1, y_2)Y_-(y_2|y_2), $$
and from conditions (ii) it follows

$$ (4.28) \qquad W_{12}^{+0}(y_1, y_2) = \frac{y_2^{-1}b(y_2)}{Z_2(y_1)} \left( \frac{q}{g(y_1, y_2)} C^{+0} - \frac{Z_2(y_2)}{y_1^{-1}b(y_1^2)} W_{11}^{+0}(y_1, y_2) \right), $$

where $C^{+,0}$ is a free scalar and $W_{22}^{+,0}(y_1, y_2)$, $W_{11}^{+,0}(y_1, y_2)$ are arbitrary. Then, one finds the following nonzero entries of $\tilde{\mathcal{M}}^{+,0}$:

$$ (4.29) \qquad \tilde{\mathcal{M}}_{11}^{+,0} = - \frac{1}{q^4 r_0^4 \mathbf{b}\mathbf{b}^* (\mathbf{b}^\diamond)^2} y_1^{-1} b(y_1^2) Z_1(y_1) C^{+,0}, \qquad \tilde{\mathcal{M}}_{12}^{+,0} = \frac{1}{q^4 r_0^4 \mathbf{b}\mathbf{b}^* (\mathbf{b}^\diamond)^2} y_2^{-1} b(y_2^2) Z_1(y_2) C^{+,0}. $$

Comparing $\det \tilde{\mathcal{M}}_{[2,1]}^{+,0} = \tilde{\mathcal{M}}_{12}^{+,0}$ to (3.55) we finally identify $\psi_0^+$:

$$ (4.30) \qquad \psi_0^+ = - \frac{r_0 \mathbf{b}\mathbf{b}^\diamond}{C^{+,0}}. $$

Now, consider $M=1$. Firstly, one shows

$$ (4.31) \qquad g(y_1, y_2) Y_+(y_1|y_2) + \frac{(q+q^{-1})^2}{\rho b(y_1^2)b(q^2 y_1^2)} \left( \eta^* + \eta \frac{q y_1^2 + q^{-1} y_1^{-2}}{q + q^{-1}} \right) - \theta_1 = - \frac{g(y_1, y_2)}{q^5 r_0^4 \mathbf{b}\mathbf{b}^* (\mathbf{b}^\diamond)^2} Z_1(y_2) Z_2(y_1) $$

and proceed as before. From conditions (i) we find

$$ (4.32) \qquad W_{21}^{+,1}(y_1, y_2) = - \frac{y_1^{-1} b(y_1^2) Z_1(y_1)}{y_2^{-1} b(y_2^2) Z_1(y_2)} W_{22}^{+,1}(y_1, y_2) $$

and from conditions (ii) it follows

$$ (4.33) \qquad W_{12}^{+,1}(y_1, y_2) = \frac{y_2^{-1}b(y_2^2)}{Z_1(y_1)} \left( \frac{q}{g(y_1, y_2)} C^{+,1} - \frac{Z_1(y_2)}{y_1^{-1}b(y_1^2)} W_{11}^{+,1}(y_1, y_2) \right) $$

where $C^{+,1}$ is a free scalar and $W_{22}^{+,1}(y_1, y_2)$, $W_{11}^{+,1}(y_1, y_2)$ are arbitrary. Then, one finds the following nonzero entries of $\tilde{\mathcal{M}}^{+,1}$:

$$ (4.34) \qquad \tilde{\mathcal{M}}_{11}^{+,1} = - \frac{1}{q^4 r_0^4 \mathbf{b}\mathbf{b}^* (\mathbf{b}^\diamond)^2} y_1^{-1} b(y_1^2) Z_2(y_1) C^{+,1}, \qquad \tilde{\mathcal{M}}_{12}^{+,1} = \frac{1}{q^4 r_0^4 \mathbf{b}\mathbf{b}^* (\mathbf{b}^\diamond)^2} y_2^{-1} b(y_2^2) Z_2(y_2) C^{+,1}. $$

Comparing $\det \tilde{\mathcal{M}}_{[2,1]}^{+,1} = \tilde{\mathcal{M}}_{12}^{+,1}$ to (3.56) we finally identify $\psi_1^+$:

$$ (4.35) \qquad \psi_1^+ = - \frac{r_0 \mathbf{b}^\diamond}{C^{+,1}}. $$

## 5. CONCLUDING REMARKS

In this paper, the Leonard pair/algebraic Bethe ansatz correspondence initiated in [BaP19, BP22] is further studied. The main results are explicit expressions for normalized scalar products of on-shell versus off-shell Bethe states in terms of q-Racah polynomials, see Theorem 1. This extends the results of [BP22] obtained for the on-shell versus on-shell Bethe states. Upon specializations, the normalized scalar products here computed are shown to satisfy a Belliard-Slavnov system of linear equations, see Theorem 2. This implies the existence of a determinant formula for q-Racah polynomials, which involves Bethe roots satisfying a Bethe equation of inhomogeneous type; See Corollary 4.1. Also, a set of relations that determine the solutions of Bethe equations of inhomogeneous type in terms of solutions of Bethe equations of homogeneous type is obtained, see Proposition 3.4.

The results obtained show that the finite dimensional representation theory of the Askey-Wilson algebra (i.e. the theory of Leonard pairs and triples) gives: (i) a way to compute normalized scalar products of on-shell and off-shell Bethe states that solve a class of Belliard-Slavnov systems of linear equations; (ii) a way to relate a class of inhomogenous Bethe equations to homogenous Bethe equations by specialization of scalar products. Reciprocally, the approach followed leads to (iii) a new interpretation of discrete (q-Racah) orthogonal polynomials of the Askey-scheme as determinants of matrices associated with Bethe roots of inhomogeneous type.

Some perspectives are now presented. Firstly, it is clear that the Leonard pair/algebraic Bethe ansatz correspondence for the Racah case ($q=1$) can be completed along the same lines starting from the material in [BerCSV20, BerCCV22]. Secondly, it is expected that the Leonard pair/algebraic Bethe ansatz correspondence
generalizes to a tridiagonal pair/algebraic Bethe ansatz correspondence. Indeed, it is known that the Askey-Wilson algebra with generators $A, A^*$ is a quotient of the $q$-Onsager algebra $O_q$ with generators $W_0, W_1$. From the point of view of the algebraic Bethe ansatz, it is possible to show that the exchange relations (2.2)-(2.8) are solved by dynamical operators expressed in terms of $W_0, W_1$. Now, recall that irreducible finite dimensional representations of $O_q$ are classified according to the theory of tridiagonal pairs of $q$-Racah type [IT09] that generalizes the theory of Leonard pairs. Similarly to the Leonard pair/algebraic Bethe ansatz correspondence, a tridiagonal pair/algebraic Bethe ansatz correspondence is thus expected. From the point of view of representation theory and special functions, this correspondence should lead to some interesting results. For instance, for a tridiagonal pair of $q$-Racah type $(\pi(W_0), \pi(W_1))$ with $\pi: O_q \to \text{End}(V)$, a characterization of the entries of the transition matrices relating certain eigenbases of $\pi(W_0)$ to certain eigenbases of $\pi(W_1)$ in terms of special functions is an open problem - see however [BVZ16] for some examples, and [CGT25] for the case of tridiagonal pairs of type II. Within a tridiagonal pair/algebraic Bethe ansatz correspondence, we expect the special functions of interest - generalizing the $q$-Racah polynomials - can be computed in terms of normalized scalar products of on-shell Bethe states. From the point of view of quantum integrable models, a tridiagonal pair/algebraic Bethe ansatz correspondence should also give new insights for the open XXZ spin chain with integrable boundary conditions or certain partition functions of the six-vertex model generalizing [BelPS24]. In particular, for the case of special right (or left) boundary conditions for the open XXZ spin chain, the spectral problem for the Hamiltonian should admit a solution either associated with Bethe equations of homogeneous type, or associated with Bethe equations of inhomogenous type. In both cases, the on-shell Bethe states diagonalizing the Hamiltonian should find an interpretation in terms of eigenvectors of $\pi(W_0)$ (or $\pi(W_1)$). Related scalar products should also find an interpretation in terms of special functions generalizing the $q$-Racah polynomials.

Some of these problems will be discussed elsewhere.

**Acknowledgments:** We thank S. Belliard for discussions and N. Crampé for comments on the manuscript. P.B. is supported by C.N.R.S. R.A.P. thanks the Institut Denis-Poisson for hospitality.

# APPENDIX A. COEFFICIENTS OF THE EXCHANGE RELATIONS

Let $\alpha, \beta \in \mathbb{C}^*$ or $\alpha = 0, \beta \in \mathbb{C}^*$ or $\alpha \in \mathbb{C}^*, \beta = 0$. The coefficients of the exchange relations (2.2)-(2.8) are given by

$$f(u, v) = \frac{b(qv/u)b(uv)}{b(v/u)b(quv)}, \quad h(u, v) = \frac{b(q^2uv)b(qu/v)}{b(quv)b(u/v)},$$

$$g(u, v, m) = \frac{\gamma(u/v, m+1) b(q)b(v^2)}{\gamma(1, m+1) b(quv^2) b(\frac{u}{v})}, \quad w(u, v, m) = -\frac{\gamma(uv, m) b(q)}{\gamma(1, m+1) b(quv)},$$

$$k(u, v, m) = \frac{\gamma(v/u, m+1) b(q)b(q^2u^2)}{\gamma(1, m+1) b(qu^2) b(\frac{u}{v})}, \quad n(u, v, m) = \frac{\gamma(1/(uv), m+2) b(q)b(v^2) b(q^2u^2)}{\gamma(1, m+1) b(qu^2) b(qv^2) b(quv)},$$

$$q(u, v, m) = \frac{\gamma(u/v, m) b(q)b(uv)}{\gamma(1, m+1)b(u/v)b(quv)}, \quad r(u, v, m) = \frac{b(q)b(u^2)\gamma(1, m)\gamma(v/u, m+1)}{\gamma(1, m+1)^2 b(qu^2) b(v/u)},$$

$$s(u, v, m) = \frac{b(q)^2 b(u^2) \gamma(v^{-2}, m+1) \gamma(v/u, m+1)}{\gamma(1, m+1)^2 b(qu^2) b(qv^2) b(\frac{u}{v})}, \quad x(u, v, m) = \frac{b(q)b(u^2)b(qu/v)\gamma(1/(uv), m+1)}{\gamma(1, m+1)b(qu^2)b(u/v)b(quv)},$$

$$y(u, v, m) = -\frac{b(q)^2 \gamma(v^{-2}, m+1) \gamma(uv, m)}{\gamma(1, m+1)^2 b(qv^2) b(quv)}, \quad z(u, v, m) = -\frac{b(q)\gamma(1, m)\gamma(uv, m)}{\gamma(1, m+1)^2 b(quv)}.$$

where

$$(A.1) \qquad \gamma^{\epsilon}(u, m) = \alpha^{\frac{1-\epsilon}{2}} \beta^{\frac{\epsilon+1}{2}} q^{-m} u - \alpha^{\frac{\epsilon+1}{2}} \beta^{\frac{1-\epsilon}{2}} q^m u^{-1}.$$

# APPENDIX B. EXAMPLE OF A LEONARD TRIPLE AND EIGENBASES FOR $s = 1/2$

B.1. A Leonard triple from a spin-$s$ representation of $U_q(sl_2)$. Recall the quantum algebra $U_q(sl_2)$ with generators $S_{\pm}, s_3$ and defining relations:

$$(B.1) \qquad [s_3, S_{\pm}] = \pm S_{\pm}, \quad [S_{+}, S_{-}] = \frac{q^{2s_3} - q^{-2s_3}}{q - q^{-1}}.$$
Now let $(\pi, \mathcal{V})$ denote the irreducible (spin-$s$) representation of $U_q(sl_2)$ of dimension $\dim(\mathcal{V}) = 2s+1$ on which the generators $\{q^{\pm s_3}, S_{\pm}\}$ act as:

$$ \pi(q^{\pm s_3}) = \sum_{k=1}^{2s+1} q^{\pm(s+1-k)} E_{kk}, \quad \pi(S_{-}) = \sum_{k=1}^{2s} \sqrt{[k]_q[2s+1-k]_q} E_{k+1k}, \quad \pi(S_{+}) = \sum_{k=1}^{2s} \sqrt{[k]_q[2s+1-k]_q} E_{kk+1}, $$

where the matrix $E_{ij}$ has a unit at the intersection of the $i$-th row and the $j$-th column, and all other entries are zero. The following gives an example of Leonard triple $A, A^*, A^\circ$ of q-Racah type satisfying the Askey-Wilson relations (3.5):

$$ A = -(q-q^{-1})r_0^{-1/2}(b^{\diamond})^{1/2}q^{s+1/2}\pi(S_+q^{s_3}) + (q-q^{-1})r_0^{-1/2}(c^{\diamond})^{1/2}q^{-s-1/2}\pi(S_-q^{s_3}) + \theta_s\pi(q^{2s_3}), $$

$$ A^* = -(q-q^{-1})r_0^{-1/2}(c^{\diamond})^{1/2}q^{-s-1/2}\pi(S_+q^{-s_3}) + (q-q^{-1})r_0^{-1/2}(b^{\diamond})^{1/2}q^{s+1/2}\pi(S_-q^{-s_3}) + \theta_s^*\pi(q^{-2s_3}). $$

The expression for $A^\circ$ is computed using (3.3). One gets$^6$:

$$ \begin{aligned} A^{\circ} = (q^{2s-1}b^{\diamond} - q^{1-2s}c^{\diamond}) \pi \left( \frac{q^{2s_3} - q^{-2s_3}}{q + q^{-1}} \right) &+ r_0^{-1} (q - q^{-1})^2 \pi(S_-^2) - \frac{(q - q^{-1})^2}{q + q^{-1}} \theta_s^{\diamond} \pi(S_- S_+) \\ &+ (q - q^{-1}) r_0^{1/2} \pi(S_-) \left( q^{-s+1/2}(c^{\diamond})^{1/2} \theta_s^* \pi(q^{-s_3}) + q^{s-1/2}(b^{\diamond})^{1/2} \theta_s \pi(q^{s_3}) \right) \\ &+ \frac{r_0}{q + q^{-1}} \left( \theta_s \theta_s^* + \frac{\omega^{\{\cdot,*,\diamond\}}}{(q - q^{-1})^2} \right) \mathbb{I}. \end{aligned} $$

**Remark 8.** With respect to the conventions chosen in [BaP19], we have the following identifications:

$$ b \rightarrow \frac{1}{2} e^{-\mu} q^{\frac{1}{2}(\nu+\nu')-2s}, \quad b^* \rightarrow \frac{1}{2} e^{\mu'} q^{\frac{1}{2}(\nu+\nu')-2s}, \quad b^{\diamond} \rightarrow \frac{1}{2} v^2 q^{\frac{1}{2}(\nu+\nu')-2s}, $$

$$ c \rightarrow \frac{1}{2} e^{\mu} q^{\frac{1}{2}(\nu+\nu')+2s}, \quad c^* \rightarrow \frac{1}{2} e^{-\mu'} q^{\frac{1}{2}(\nu+\nu')+2s}, \quad c^{\diamond} \rightarrow \frac{1}{2} v^{-2} q^{\frac{1}{2}(\nu+\nu')+2s}, $$

$$ (B.2) \qquad r_0 \rightarrow 2q^{-\frac{1}{2}(\nu+\nu')}, \quad v \rightarrow \zeta^{-1}, \quad \chi \rightarrow -\frac{1}{2}q^{\nu}(q^2 - q^{-2}), $$

where $\{\mu, \mu', \nu, \nu', v\}$ are complex parameters.

B.2. **Eigenvectors and their duals for $s = 1/2$.** In this subsection, as a simple example, the eigenvectors of the Leonard triple $(A, A^*, A^\circ)$ and respective normalizations for the case $s = 1/2$ are given.

From (B.2), the eigenvectors of $A$ are

$$ (B.3) \qquad |\theta_0\rangle = n_0 \left( \frac{\frac{q^{-1/2}r_0^{-1/2}(b^{\diamond})^{1/2}}{b}}{1} \right), \quad |\theta_1\rangle = n_1 \left( \frac{q^{3/2}r_0^{3/2}(b^{\diamond})^{1/2}b}{1} \right), $$

where $n_{0,1}$ are some scalars. The action (2.13) fixes the ratio

$$ (B.4) \qquad \frac{n_1}{n_0} = - \frac{b^{\diamond} + qr_0 b b^*}{qr_0 b (b^* + qr_0 b b^{\diamond})}. $$

For the dual eigenvectors, one has

$$ (B.5) \qquad \langle \theta_0 | = m_0 \left( 1, -q^{3/2} r_0^{3/2} (b^{\diamond})^{1/2} b \right), \quad \langle \theta_1 | = m_1 \left( 1, -\frac{q^{-1/2} r_0^{-1/2} (b^{\diamond})^{1/2}}{b} \right), $$

where $m_{0,1}$ are some scalars. The dual action (2.17) fixes

$$ (B.6) \qquad \frac{m_1}{m_0} = \left( \frac{\xi_1}{\xi_0} \right) \frac{qr_0 b(b^* + qr_0 b b^*)}{(b^* + qr_0 b b^*)}. $$

Choosing

$$ (B.7) \qquad m_0 = \frac{q^{1/2} r_0^{1/2} b \xi_0}{(b^{\diamond})^{1/2} (1 - q^2 r_0^2 b^2) n_0}, $$

one has - see (2.19):

$$ (B.8) \qquad \langle \theta_0 | \theta_0 \rangle = \xi_0, \quad \langle \theta_1 | \theta_1 \rangle = \xi_1. $$

<footnote>$^6$For simplicity, we choose the branch cut $\sqrt{1/(r_0^2 b^\diamond)} = 1/(r_0 \sqrt{b^\diamond})$.</footnote>
From (B.2), the eigenvectors of $A^*$ are

$$ (B.9) \qquad |\theta_0^*\rangle = n_0^* \left( -\frac{q^{1/2} r_0^{1/2} b^*}{(b^\diamond)^{1/2}} \middle| \frac{1}{1} \right), \quad |\theta_1^*\rangle = n_1^* \left( -\frac{1}{q^{3/2} r_0^{3/2} (b^\diamond)^{1/2} b^*} \middle| \frac{1}{1} \right), $$

where $n_{0,1}^*$ are some scalars. The action (2.14) fixes

$$ (B.10) \qquad \frac{n_1^*}{n_0^*} = - \frac{qr_0 b^* (b^\diamond + qr_0 b^* b^\diamond)}{(b + qr_0 b^* b^\diamond)}. $$

For the dual vectors, one has

$$ (B.11) \qquad \langle \theta_0^* | = m_0^* \left( 1, \frac{1}{q^{3/2} r_0^{3/2} (b^\diamond)^{1/2} b^*} \right), \quad \langle \theta_1^* | = m_1^* \left( 1, \frac{q^{1/2} r_0^{1/2} b^*}{(b^\diamond)^{1/2}} \right), $$

where $m_{0,1}^*$ are some scalars. The dual action (2.18) fixes

$$ (B.12) \qquad \frac{m_1^*}{m_0^*} = \left( \frac{\xi_1^*}{\xi_0^*} \right) \frac{(b + qr_0 b^* b^\diamond)}{qr_0 b^* (b^\diamond + qr_0 b^* b^*)}. $$

Choosing

$$ (B.13) \qquad m_0^* = \frac{q^{3/2} r_0^{3/2} (b^\diamond)^{1/2} b^* \xi_0^*}{(1 - q^2 r_0^2 (b^*)^2) n_0^*}, $$

one has

$$ (B.14) \qquad \langle \theta_0^* | \theta_0^* \rangle = \xi_0^*, \quad \langle \theta_1^* | \theta_1^* \rangle = \xi_1^*. $$

From (B.2), the eigenvectors of $A^\diamond$ are

$$ (B.15) \qquad |\theta_0^\diamond\rangle = n_0^\diamond \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad |\theta_1^\diamond\rangle = n_1^\diamond \left( \frac{q^{1/2} r_0^{1/2} b^* (q^2 r_0^2 (b^\diamond)^2 - 1)}{(b^\diamond)^{1/2} (q^2 r_0^2 b^* (b^\diamond + qr_0 b^\diamond) + b + qr_0 b^\diamond b^*)} \middle| \frac{1}{1} \right), $$

where $n_{0,1}^\diamond$ are some scalars. The action (3.26) for $(b, c) = (*, \diamond)$ fixes the ratio

$$ (B.16) \qquad \frac{n_1^\diamond}{n_0^\diamond} = - \frac{q^2 r_0^2 b^* (b^* + qr_0 b^* b^\diamond) + b + qr_0 b^* b^*}{qr_0 (b^* + qr_0 b^* b^\diamond) (b^\diamond + qr_0 b^* b^*)}. $$

The transitions (3.31) and (3.32) for $(a, b, c) = (\cdot, *, \diamond)$ determine the ratios

$$ (B.17) \qquad \frac{n_0^*}{n_0^\diamond} = \frac{1}{g_0(1 - q^2 r_0^2 (b^*)^2)}, \quad \frac{n_0}{n_0^\diamond} = - \frac{qr_0 b (b^* + qr_0 b b^\diamond)}{g_0 b^\diamond (q^2 r_0^2 b^2 - 1) (q^2 r_0^2 (b^*)^2 - 1)}. $$

while (3.33) gives (3.34) for $s = 1/2$:

$$ (B.18) \qquad f_0 h_0^{-1} g_0^{-1} = - \frac{(1 - q^2 r_0^2 b^2)(1 - q^2 r_0^2 (b^\diamond)^2)(1 - q^2 r_0^2 (b^*)^2)}{\left( \frac{qr_0 b b^\diamond}{b^*} + 1 \right) \left( \frac{qr_0 b b^*}{b^\diamond} + 1 \right) \left( \frac{qr_0 b^\diamond b^*}{b} + 1 \right)}. $$

### B.3. Examples of scalar products. For the normalized scalar products

$$ (B.19) \qquad \frac{\langle \theta_M | \theta_N^* \rangle}{\langle \theta_M | \theta_M \rangle}, $$

we obtain:

$$ (B.20) \qquad \frac{\langle \theta_0 | \theta_0^* \rangle}{\langle \theta_0 | \theta_0 \rangle} = 1, \quad \frac{\langle \theta_0 | \theta_1^* \rangle}{\langle \theta_0 | \theta_0 \rangle} = - \frac{(b^\diamond + qr_0 b^* b^*) (1 + q^3 r_0^3 b^* b^\diamond)}{qr_0 (b + qr_0 b^* b^\diamond) (b^* + qr_0 b^* b^\diamond)}, \quad \frac{\langle \theta_1 | \theta_0^* \rangle}{\langle \theta_1 | \theta_1 \rangle} = 1, \quad \frac{\langle \theta_1 | \theta_1^* \rangle}{\langle \theta_1 | \theta_1 \rangle} = -1. $$

From those expressions, the $q$-Racah polynomials can be computed (see [T03, Theorem 14.6 and 15.6]):

$$ (B.21) \qquad R_M^{(\cdot, \diamond)}(\theta_N^*) = \frac{\langle \theta_M | \theta_N^* \rangle}{\langle \theta_0 | \theta_N^* \rangle} \frac{\langle \theta_0 | \theta_0 \rangle}{\langle \theta_M | \theta_M \rangle} = \frac{\langle \theta_N^* | \theta_M \rangle}{\langle \theta_0^* | \theta_0^* \rangle} \frac{\langle \theta_0^* | \theta_N^* \rangle}{\langle \theta_N^* | \theta_N^* \rangle}. $$
Using (B.20), it follows $R_0^{\{\cdot,\circ\}}(\theta_0^*) = R_0^{\{\cdot,\circ\}}(\theta_1^*) = R_1^{\{\cdot,\circ\}}(\theta_0^*) = 1$ and

$$ (B.22) \qquad \begin{aligned} R_1^{\{\cdot,\circ\}}(\theta_1^*) &= \frac{\langle\theta_1|\theta_1^*\rangle \langle\theta_0|\theta_0\rangle}{\langle\theta_0|\theta_1^*\rangle \langle\theta_1|\theta_1\rangle} = \frac{\langle\theta_1^*|\theta_1\rangle \langle\theta_0^*|\theta_0^*\rangle}{\langle\theta_0^*|\theta_1\rangle \langle\theta_1^*|\theta_1^*\rangle} = \frac{qr_0(\mathbf{b} + qr_0\mathbf{b}^*\mathbf{b}^*)(\mathbf{b}^* + qr_0\mathbf{b}\mathbf{b}^*)}{(\mathbf{b}^* + qr_0\mathbf{b}\mathbf{b}^*)(1 + q^3r_0^3\mathbf{b}\mathbf{b}^*\mathbf{b}^*)} \\ &= 1 - \frac{(q^{-2}, \frac{\mathbf{b}}{c}q^2; q^2)_1 q^2}{(-\frac{\mathbf{b}}{c}^*r_0q, -\frac{\mathbf{b}}{c}^*r_0q^3, q^{-2}; q^2)_1(q^2; q^2)_1 c^*} (\theta_1^* - \mathbf{b}^* - \mathbf{c}^*), \end{aligned} $$

where the second line matches with (3.30) for $M = N = 1$ and $s = 1/2$.

## REFERENCES

[ABGP15] J. Avan, S. Belliard, N. Grosjean and R.A. Pimenta, Modified algebraic Bethe ansatz for XXZ chain on the segment - III - Proof, Nucl. Phys. B 899 (2015) 229, arXiv:1506.02147.

[BaP19] P. Baseilhac and R.A. Pimenta, Diagonalization of the Heun-Askey-Wilson operator, Leonard pairs and the algebraic Bethe ansatz, Nucl. Phys. B 949 (2019) 114824, arXiv:1909.02464.

[BP22] P. Baseilhac and R.A. Pimenta, The q-Racah orthogonal polynomials from scalar products of Bethe states, J. Phys. A 56 (2023) 385202, arXiv:2211.14727.

[BS09] P. Baseilhac and K. Shigechi, A new current algebra and the reflection equation, Lett. Math. Phys. 92 (2010), 47-65, arXiv:0906.1482.

[BVZ16] P. Baseilhac, L. Vinet and A. Zhedanov, The q-Onsager algebra and multivariable q-special functions, J. Phys. A: Math. Theor. 50 (2017) 395201, arXiv:1611.09250.

[Bel15] S. Belliard, Modified algebraic Bethe ansatz for XXZ chain on the segment - I: Triangular cases, Nucl. Phys. B 892 (2015) 1, arXiv:1408.4840.

[BelS19] S. Belliard and N. Slavnov, Why scalar products in the algebraic Bethe ansatz have determinant representation, J. High Energ. Phys. 10 (2019) 103, arXiv:1908.00032.

[BelC13] S. Belliard and N. Crampé, Heisenberg XXX Model with General Boundaries: Eigenvectors from Algebraic Bethe Ansatz, SIGMA 9 (2013), 072, arXiv:1309.6165.

[BelP15] S. Belliard and R.A. Pimenta, Modified algebraic Bethe ansatz for XXZ chain on the segment - II - general cases, Nucl. Phys. B 894 (2015) 527, arXiv:1412.7511.

[BelPS21] S. Belliard, R.A. Pimenta and N.A. Slavnov, Scalar product for the XXZ spin chain with general integrable boundaries, J. Phys. A: Math. Theor. 54 (2021) 344001, arXiv:2103.12501.

[BelPS24] S. Belliard, R.A. Pimenta and N.A. Slavnov, Modified rational six vertex model on the rectangular lattice, SciPost Phys. 16 (2024) 009, arXiv:2310.05850.

[BerCCV22] P-A. Bernard, G. Carcone, N. Crampé and L. Vinet, Bethe ansatz diagonalization of the Heun-Racah operator, Lett. Math. Phys. 113 (2023) 8, arXiv:2209.09213.

[BerCSV20] P-A. Bernard, N. Crampé, D. Shaaban Kabakibo and L. Vinet, Heun operator of Lie type and the modified algebraic Bethe ansatz, J. Math. Phys. 62 (2021) 083501, arXiv:2011.11659.

[C15] N. Crampé, Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries, J. Phys. A: Math. Theor. 48 (2015) 08FT01, arXiv:1411.7954.

[CGT25] N. Crampé, J. Gaboriaud and S. Tsujimoto, Change of basis for the tridiagonal pairs of type II, arXiv:2503.01231.

[CLSW03] J. Cao, H.-Q. Lin, K. Shi and Y. Wang, Exact solution of XXZ spin chain with unparallel boundary fields, Nucl. Phys. B663 (2003) 487.

[CYSW14] J. Cao, W-L. Yang, K. Shi and Y. Wang, Exact solution of the XXZ alternating spin chain with generic non-diagonal boundaries, Ann. of Phys. 361 (2015) 91, arXiv:1409.3646.

[Curt07] B. Curtin, Modular Leonard triples, Linear Algebra and its Applications. 424 (2007) 510-539.

[D07] A. Doikou, A note on the boundary spin s XXZ chain, Phys. Lett. A, 366 (2007) 556-562, arXiv:hep-th/0612268.

[FK10] G. Filali and N. Kitanine, Partition function of the trigonometric SOS model with reflecting end, J. Stat. Mech. (2010) L06001, arXiv:1004.1015.

[GZ88] Y.A. Granovskii and A. Zhedanov, Nature of the symmetry group of the 6j-symbol, Journal of Experimental and Theoretical Physics 94 (1988) 49-54.

[Gu08] S-J. Gu, Fidelity approach to quantum phase transitions, Int. J. Mod. Phys. B 24 (2010) 4371-4458, arXiv:0811.3127.

[H11] H-W. Huang, The classification of Leonard triples of QRacah type, Linear Alg. Appl. 436 (2012) 1442-1472, arXiv:1108.0458.

[IT09] T. Ito and P. Terwilliger, The augmented tridiagonal algebra, Kyushu Journal of Mathematics 64 (2010) No. 1 81-144, arXiv:0904.2889v1.

[KS96] R. Koekoek and R. Swarttuw, The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue, arXiv:math.ca/9602214v1.

[K10] T.H. Koornwinder, On the limit from q-Racah polynomials to big q-Jacobi polynomials, SIGMA 7 (2011) 040, 8 pages, arXiv:1011.5585v4.

[LL65] J.-M. Lévy-Leblond J.-M. and M. Lévy-Nahas, Symmetrical coupling of three angular momenta, J. Math. Phys. 6 (1965), 1372-1380.

[Sk88] E. K. Sklyanin, Boundary conditions for integrable quantum systems, J. Phys. A 21 (1988), 2375-2389.

[MNP17] J. M. Maillet, G. Niccoli and B. Pezelier, Transfer matrix spectrum for cyclic representations of the 6-vertex reflection algebra I, SciPost Phys. 2 (2017) 009, arXiv:1802.08853.
[STF79] E. K. Sklyanin, L. A. Takhtadzhyan and L. D. Faddeev, *The Quantum Inverse Problem Method. I*, Theor. Math. Phys. 40 (1979), 688.

[Sl22] N. Slavnov, *Algebraic Bethe ansatz and correlation functions: an advanced course*, World Scientific Publishing (2022).

[T87] P. Terwilliger, *Leonard pairs and dual polynomial sequences*, Preprint 1987.

[T99] P. Terwilliger, *Two relations that generalize the q-Serre relations and the Dolan-Grady relations*, Proceedings of the Nagoya 1999 International workshop on physics and combinatorics. Editors A. N. Kirillov, A. Tsuchiya, H. Umemura. 377-398, arXiv:math/0307016.

[T01] P. Terwilliger, *Two linear transformations each tridiagonal with respect to an eigenbasis of the other*, Linear Alg. and Appl. 330 (2001) 149-203, arXiv:math/0406555.

[T03] P. Terwilliger, *Leonard pairs and the q-Racah polynomials*, Linear Alg. and Appl. 387 (2004) 235-276, arXiv:math/0306301.

[T04] P. Terwilliger, *Leonard pairs from 24 points of view*, Rocky Mountain J. Math. 32 (2002), 827-888, arXiv:math/0406577.

[T21] P. Terwilliger, *The alternating central extension of the q-Onsager algebra*, Comm. Math. Phys. 387 (2021) 1771-1819; arXiv:2103.03028.

[TF79] L.A. Takhtadzhan and L.D. Faddeev, *The quantum method of the inverse problem and the Heisenberg XYZ model*, Russian Math. Surveys 34 5 (1979), 11-68.

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

[WYCS15] Y. Wang, W.-L. Yang, J. Cao and K. Shi, *Off-Diagonal Bethe Ansatz for Exactly Solvable Models*, Springer-Verlag Berlin Heidelberg, (2015).

[ZCYSW17] X. Zhang, J. Cao, W.-L. Yang, K. Shi and Y. Wang, *Exact solution of the relativistic quantum Toda chain*, J. Phys. A: Math. Theor. 50 (2017) 124003, arXiv:1609.07385 [math-ph].

[Zh91] A.S. Zhedanov, "*Hidden symmetry*" of the Askey-Wilson polynomials, Theoret. and Math. Phys. 89 (1991) 1146-1157.

* INSTITUT DENIS-POISSON CNRS/UMR 7013 - UNIVERSITÉ DE TOURS - UNIVERSITÉ D'ORLÉANS PARC DE GRAMMONT, 37200 TOURS, FRANCE

Email address: pascal.baseilhac@univ-tours.fr

** INSTITUTO DE FÍSICA - UNIVERSIDADE FEDERAL DO RIO GRANDE DO SUL, PORTO ALEGRE, 91501-970, BRAZIL

Email address: rodrigo.pimenta@ufrgs.br