Pascal Baseilhac
email: pascal.baseilhac@idpoisson.fr

Rodrigo A Pimenta
email: pimenta@ifsc.usp.br

DIAGONALIZATION OF THE HEUN-ASKEY-WILSON OPERATOR, LEONARD PAIRS AND THE ALGEBRAIC BETHE ANSATZ

Keywords: 81R50, 81R10, 81U15 Askey-Wilson algebra, Leonard pairs, Reflection equation, Bethe ansatz, Integrable systems

An operator of Heun-Askey-Wilson type is diagonalized within the framework of the algebraic Bethe ansatz using the theory of Leonard pairs. For different specializations and the generic case, the corresponding eigenstates are constructed in the form of Bethe states, whose Bethe roots satisfy Bethe ansatz equations of homogeneous or inhomogenous type. For each set of Bethe equations, an alternative presentation is given in terms of 'symmetrized' Bethe roots. Also, two families of on-shell Bethe states are shown to generate two explicit bases on which a Leonard pair acts in a tridiagonal fashion. In a second part, the (in)homogeneous Baxter T-Q relations are derived. Certain realizations of the Heun-Askey-Wilson operator as second q-difference operators are introduced. Acting on the Q-polynomials, they produce the T-Q relations. For a special case, the Q-polynomial is identified with the Askey-Wilson polynomial, which allows one to obtain the solution of the associated Bethe ansatz equations. The analysis presented can be viewed as a toy model for studying integrable models generated from the Askey-Wilson algebra and its generalizations. As examples, the q-analog of the quantum Euler top and various types of three-sites Heisenberg spin chains in a magnetic field with inhomogeneous couplings, three-body and boundary interactions are solved. Numerical examples are given. The results also apply to the time-band limiting problem in signal processing.

Introduction

Although not pointed out in standard textbooks of quantum mechanics, the quantum harmonic oscillator is among the basic examples of quantum integrable systems generated from the so-called Askey-Wilson algebra introduced in [Z91] (see (1.14), (1.15) below). Introducing the Heisenberg algebra with generators a, a † and defining relations [a, a † ] = 1, for an appropriate change of variables the Hamiltonian, position and momentum operators read, respectively:

H = a † a + 1 2 , X = 1 √ 2 (a † + a) , P = i √ 2 (a † -a) . (1.1)
The triplet (H, X, iP) generates a specialization of the Askey-Wilson algebra. Considering the first presentation of the Askey-Wilson algebra with three generators given in [Z91], one routinely finds: H, X = -iP, X, iP = -1 and iP, H = X. It follows that the pair (H, X) satisfies the so-called Askey-Wilson relations [Z91]:

X, X, H = -1 , H, H, X = X (1.2)
giving a second presentation [Z91]. Alternatively, for the pair (H, iP) one also gets Askey-Wilson relations, with different structure constants:

H, H, iP = iP , iP, iP, H = 1 . (1.3)
In the literature, the analysis of the spectral problem for the Hamiltonian H is usually based on two different representations 1 of the Heisenberg algebra. There exists a basis {|θ * n , n ∈ Z + } such that the generators a, a † and the so-called number operator N = a † a act respectively as a|θ (1.4) 1 We assume the reader is familiar with the bases {|x } and {|n } of standard textbooks in quantum mechanics. For convenience, below we introduce x| = e -x 2 /2 θx| and |n = (2 n n!) -1/2 |θ * n . Note that a rigorous mathematical definition of the basis {|x } requires the framework of rigged Hilbert spaces, see [M01, M05] and references therein.

On the other hand, there exists a basis {|θ x , x ∈ R} to which one associates the linear functional denoted θ x | such that the generators a, a † and N act as

θ x |a = 1 √ 2 d dx θ x | , θ x |a † = 1 √ 2 2x - d dx θ x | , θ x |N = 1 2 2x - d dx d dx θ x | . (1.5)
Then, the transition coefficients θ x |θ * n solve a bispectral problem which reads as a second-order differential equation and a three term recurrence relation given by:

θ x |H|θ * n = (n + 1 2 ) θ x |θ * n ⇔ 2x - d dx d dx H n (x) = 2nH n (x) with H n (x) = θ x |θ * n , (1.6) θ x |X|θ * n = θ x |θ * n-1 + 2n θ x |θ * n+1 ⇔ 2xH n (x) = H n-1 (x) + 2nH n-1 (x) . (1.7)
This system coincides with the well-known defining relations of the orthogonal Hermite polynomials H n (x) [START_REF] Koekoek | The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue[END_REF]eqs. (1.13.3), (1.13.4)] located at the bottom of the Askey-scheme:

H n (x) = [n/2] k=0 (-1) k n! k!(n -2k!) (2x) n-2k e.g. H 0 (x) = 1 , H 1 (x) = 2x , H 2 (x) = 4x 2 -2, ... (1.8)
Alternatively, it is also known that the spectral problem (1.6) can be formulated within the analytical Bethe ansatz framework. Namely, consider the elementary Baxter T-Q relation:

2xQ ′ n (x) -Q ′′ n (x) = Λ n Q n (x) where Q n (x) = n i=1 (x -x i ) . (1.9)
Given n fixed, (1.9) implies that the zeroes {x i |i = 1, ..., n} satisfy the set of Bethe ansatz equations

Q ′′ n (x i ) Q ′ n (x i ) = 2x i ⇔ n j=1,j =i 1 x i -x j =
x i for all i = 1, ..., n, (1.10) and the spectrum is given by Λ n = 2n. Thus, by comparison of (1.9) with (1.6) it follows that the Q-polynomial in (1.9) can be interpreted as the the transition coefficients θ x |θ * n : Q n (x) = 2 -n θ x |θ * n = 2 -n H n (x) . (1.11) Furthermore, using (1.1) and (1.5) together with (1.4), the Q-polynomial can be written in terms of the pair (H, X) as:

Q n (x) = θ x | (X + [X, H]) n |θ * 0 ⇔ Q n (x) = 2 -n 2x - d dx n 1 . (1.12)
More generally, given the triplet (H, X, iP) satisfying the Askey-Wilson relations (1.2) and (1.3) it is natural to consider the combination (recall that iP = [X, H])

I ho = κX + κ * H + κ + [X, H] (1.13)
and the corresponding spectral problem. From an algebraic perspective, the pair (X, I ho ) (or alternatively (H, I ho )) generates the simplest specialization of the family of Heun algebras of Lie type recently introduced in [CVZ19, Definition 2.1]. For generic parameters κ, κ * , κ + , it is well-known that the spectral problem2 for I ho can still be solved in terms of Hermite polynomials (see e.g. [CDL]).

The possibility of generalizing the above picture to the Askey-Wilson algebra with generic structure constants [Z91] and to the Heun-Askey-Wilson algebra recently introduced in [START_REF] Baseilhac | The Heun-Askey-Wilson algebra and the Heun operator of Askey-Wilson type[END_REF] is a natural problem to be considered, for various reasons. Some motivations are briefly described in Section 6. Let ρ, ω, η, η * be generic scalars. The Askey-Wilson algebra is generated by A, A * subject to the relations [Z91] A, A, A * q q -1 = ρ A * + ω A + ηI , (1.14) A * , A * , A q q -1 = ρ A + ω A * + η * I .

(1.15)

In the literature, these relations are called the Askey-Wilson relations. It is known that the Askey-Wilson algebra gives an algebraic framework for all the orthogonal polynomials of the Askey-scheme [Z91]. Also, it is known that the zeroes of the Askey-Wilson polynomials satisfy a class of Bethe ansatz equations [START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF][START_REF] Van Diejen | Solutions of convex Bethe Ansatz equations and the zeros of (basic) hypergeometric orthogonal polynomials[END_REF]. Define the Heun-Askey-Wilson element [START_REF] Baseilhac | The Heun-Askey-Wilson algebra and the Heun operator of Askey-Wilson type[END_REF][START_REF] Crampé | Heun algebras of Lie type[END_REF]:

I(κ, κ * , κ + , κ -) = κ A + κ * A * + κ + χ -1 [A, A * ] q + κ -χ [A * , A] q (1.16)
where κ, κ * , κ ± and3 χ = 0 are scalars. To our knowledge, the diagonalization of a Heun-Askey-Wilson operator associated with (1.16) has not been considered yet in full generality4 . For finite dimensional representations, as we will show it turns out that this problem can be solved within the boundary quantum inverse scattering method [S88]. In fact, the Heun-Askey-Wilson operator commutes with a given transfer matrix [START_REF] Zabrodin | Quantum transfer matrices for discrete and continuous quasi-exactly solvable problems[END_REF][START_REF] Baseilhac | An integrable structure related with tridiagonal algebras[END_REF] associated with non-diagonal reflection matrices [START_REF] Vega | Boundary K-matrices for the XXX,XXZ and XYZ spin chains[END_REF][START_REF] Ghoshal | Boundary S-matrix and boundary states in two-dimensional integrable field theory[END_REF], see (2.21) below. Such kind of reflection matrices contain arbitrary constant parameters (boundary couplings) which in general break the U (1) symmetry of the transfer matrix, and it prevents the application of the standard Bethe ansatz techniques to study its diagonalization, unless that restrictions on the boundary couplings are imposed, see5 e.g. [CLSW03, N04, YZ07, D07, FNR07, BCR12, PL13]. An important progress for the unrestricted cases was achieved by the introduction of the offdiagonal Bethe ansatz [START_REF] Cao | Off-diagonal Bethe ansatz and exact solution of a topological spin ring[END_REF], a method that proposes an inhomogeneous Baxter T-Q equation as solution of the spectral problem for integrable models without U (1) symmetry [WYCS]. Beyond the computation of the spectrum, a modification of the algebraic Bethe ansatz was developed in [BC13, Be15, C14, BP15, ABGP15] providing the construction of the associated off-shell Bethe states6 , which in particular allows the computation of scalar products between on-shell/off-shell Bethe states, see [START_REF] Belliard | Slavnov and Gaudin-Korepin Formulas for Models without U (1) Symmetry: the Twisted XXX Chain[END_REF][START_REF] Belliard | Slavnov and Gaudin-Korepin Formulas for Models without U (1) Symmetry: the XXX chain on the segment[END_REF][START_REF] Belliard | Scalar Products in Twisted XXX Spin Chain[END_REF][START_REF] Belliard | Why scalar products in the algebraic Bethe ansatz have determinant representation[END_REF] and references therein. The main feature of the modified algebraic Bethe ansatz are the off-shell relations satisfied by the 'creation operators', see e.g. Lemma 3.5 and Conjecture 1. In the algebraic Bethe ansatz perspective, these off-shell relations are the origin of the inhomogeneous term in the Baxter T-Q equation.

In the present paper, for any irreducible finite dimensional representation of the Askey-Wilson algebra and q = 1, the construction of eigenvectors of (1.16) within the (modified) algebraic Bethe ansatz framework and the analog of the construction (1.9), (1.10), (1.11) is considered. In this case, the theory of Leonard pairs developed in [START_REF] Terwilliger | Leonard pairs and the q-Racah polynomials[END_REF][START_REF] Terwilliger | Leonard pairs and the Askey-Wilson relations[END_REF] offers the proper mathematical setting. The main results are the following:

(i) The Heun-Askey-Wilson operator associated with (1.16) is diagonalized on any irreducible finite dimensional representation (π, V) of dimension 2s + 1 (s is an integer or half-integer), using the theory of Leonard pairs combined with the algebraic Bethe ansatz technique. Three cases of (1.16) indexed by 'a' are considered in details: the special cases (i.e. a = sp) κ = 0, κ * = κ ± = 0 or κ * = 0, κ = κ ± = 0; the diagonal case (i.e. a = d) κ, κ * = 0, κ ± = 0; the generic case (i.e. a = g) κ, κ * , κ ± = 0. In all cases, the solution of the spectral problem7 takes the form: π (I(κ, κ * , κ + , κ -)) |Ψ M a,ǫ (ū) = Λ M a,ǫ |Ψ M a,ǫ (ū) , (1.17) where the eigenstates |Ψ M a,ǫ (ū) are given in the form of on-shell Bethe states (with M = 0, 1, ..., 2s for a = sp and M = 2s for a = d, g) obtained from the so-called 'reference state' with index ǫ ∈ {±} through successive actions of a 'creation' operators {B ǫ (u i , m i )}, thus generalizing (1.12). According to the choice of parameters and reference state, the Bethe roots ū = {u 1 , ..., u M } satisfy homogeneous or inhomogeneous Bethe ansatz equations. See Propositions 3.1, 3.2, 3.3, 3.4. For the case a = d, the inhomogeneous term is related with the characteristic polynomial of an operator denoted A ⋄ that forms a Leonard triple with A, A * . Importantly, this provides a proof of the off-shell relations in Lemma 3.5 only based on the representation theory of the Askey-Wilson algebra.

(ii) An infinite dimensional representation (π, V) of the Askey-Wilson algebra is introduced. In this representation, the operators π(A), π(A * ) act as second-order q-difference operators whereas the Heun-Askey-Wilson operator π(I(κ, κ * , κ + , κ -)) acts in general as a fourth-order q-difference operator. For the special case κ = κ ± = 0, the spectral problem for the Heun-Askey-Wilson operator q-difference operator (reduced to the one for π(A * )) produces an homogeneous Baxter T-Q relation generalizing (1.9). In this case, the Baxter Q-polynomial of degree M follows from the 'transition' coefficient (see Lemma 4.2):

Q M (Z) = N M (ū) -1 z|Ψ M sp,+ (ū) = M i=1
Z -qu 2 i + q -1 u -2 i q + q -1 where Z = z + z -1 q + q -1 , (1.18) thus generalizing (1.11). Furthermore, it solves a bispectral problem for the Askey-Wilson polynomial, as expected. For the diagonal case κ, κ * = 0, κ ± = 0, the action of the q-difference operator π(I(κ, κ * , 0, 0)) on the Baxter Q-polynomial produces an inhomogeneous Baxter T-Q relation generalizing (1.9). In this case, the Baxter Q-polynomial of degree 2s is also expressed in terms of the 'transition' coefficient (see Lemma 4.2):

Q 2s (Z) = N 2s (ū) -1 z|Ψ 2s d,+ (ū) . (1.19)
For the generic case κ, κ * , κ ± = 0 and a different realization of the Heun-Askey-Wilson operator, a similar picture holds as briefly mentionned.

(iii) For each system of homogeneous or inhomogeneous Bethe ansatz equations derived in Section 3, an alternative presentation in terms of the 'symmetrized' Bethe roots (3.103) is obtained. In particular, this presentation gives a simpler characterization of the Bethe roots than the one based on the standard Bethe equations. See Proposition 3.6.

(iv) An algebraic Bethe ansatz solution for the q-analog of the quantum Euler top [START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF][START_REF] Turbiner | The Heun operator as a Hamiltonian[END_REF] and various examples of three-sites Heisenberg spin chains with three-body terms, inhomogeneous couplings and boundary interactions. See Section 5 where numerical examples are given.

The paper is organized as follows. In Section 2, we review the derivation of the Heun-Askey-Wilson element (1.16) based on the connection between the reflection equation algebra and the Askey-Wilson algebra [START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF][START_REF] Baseilhac | An integrable structure related with tridiagonal algebras[END_REF]. Using the transfer matrix formalism and the theory of Leonard pairs [START_REF] Terwilliger | Leonard pairs and the q-Racah polynomials[END_REF][START_REF] Terwilliger | Leonard pairs and the Askey-Wilson relations[END_REF], the diagonalization of the Heun-Askey-Wilson operator associated with (1.16) is considered within the framework of the algebraic Bethe ansatz in Section 3. Namely, for various choices of parameters {κ, κ * , κ + , κ -}, the eigenstates of the Heun-Askey-Wilson operator are given in the form of Bethe states and associated Bethe equations, and in terms of Leonard pairs' data. In each case, an alternative presentation of the Bethe ansatz equations is given in terms of a system of polynomial equations in the 'symmetrized' variables (3.103). Also, for any off-shell or on-shell Bethe state an expansion formula is given using the Poincaré-Birkhoff-Witt basis of the Askey-Wilson algebra. It follows that certain families of Bethe states generate explicit bases for Leonard pairs. In Section 4, four different types of Baxter T-Q relations are derived, either homogeneous or inhomogeneous. The corresponding system of Bethe equations produces the Bethe equations derived in Section 3 and the solutions are expressed as symmetric polynomials, identified as Q-polynomials. Independently, a second-order q-difference operator realization of the Heun-Askey-Wilson element (1.16) is given for the special and diagonal cases. It is shown that its action on the Q-polynomial produces the T-Q relations. For the generic case, a different realization has to be considered, as briefly mentionned. For the special case κ = κ ± = 0 the Q-polynomial is expressed in terms of the orthogonal Askey-Wilson polynomials. An interpretation of the Q-polynomial as transition coefficients arises naturally. In Section 5, we apply the previous results to the diagonalization of various examples of integrable models: a q-analog of the quantum Euler top and Heisenberg spin chains with three-sites, that can be viewed as deformations of the three-sites U q (sl 2 )-invariant XXZ spin chain. In Section 6, some perspectives are briefly described. In Appendices A, B, formulas are collected. In Appendix C, the proof of Lemma 3.5 is given. Appendix D is devoted to the proof of Proposition 3.5. In Appendix E, different realizations of the Askey-Wilson algebra in terms of first or second-order q-difference operators are given. In Appendix F, realizations of the dynamical operators in terms of the q-difference operators are displayed.

Notations:

The parameter q is assumed not to be a root of unity and is different than 1. We write [X, Y ] q = qXYq -1 Y X, the commutator [X, Y ] = [X, Y ] q=1 . The identity element is denoted I. We will use the standard q-shifted factorials (also called q-Pochhammer functions) [START_REF] Koekoek | The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue[END_REF], q-number and binomial, respectively:

(a; q) n = n-1 k=0 (1 -aq k ), [n] q = q n -q -n q -q -1 , [0] q = 1 , m n = m! n!(m -n)! . (1.20)
We also use the notation

b(x) = x -x -1 . (1.21)
For the elementary symmetric polynomials in the variables {x i |i = 1, ..., n}, we use the notation:

e k (x 1 , x 2 , ..., x n ) = 1≤j1<j2<•••<j k ≤n x j1 x j2 • • • x j k (1.22)
2. The Heun-Askey-Wilson element from the reflection equation In this section, we show how the Heun-Askey-Wilson element (1.16) introduced in [START_REF] Baseilhac | The Heun-Askey-Wilson algebra and the Heun operator of Askey-Wilson type[END_REF] follows from the transfer matrix associated with the solution of the reflection equation constructed in [START_REF] Zabrodin | Quantum transfer matrices for discrete and continuous quasi-exactly solvable problems[END_REF][START_REF] Baseilhac | An integrable structure related with tridiagonal algebras[END_REF]. Most of the material given in this section is taken from the works [S88, Za95, B04, BK07], so we skip the details.

Let the operator-valued function R : C * → End(C 2 ⊗ C 2 ) be the intertwining operator (quantum R-matrix) between the tensor product of two fundamental representations associated with the quantum algebra U q ( sl 2 ). The element R(u) depends on the deformation parameter q and is defined by

R(u) =     uq -u -1 q -1 0 0 0 0 u -u -1 q -q -1 0 0 q -q -1 u -u -1 0 0 0 0 uq -u -1 q -1     , (2.1)
where u is the so-called the spectral parameter. Then R(u) satisfies the quantum Yang-Baxter equation in the space

C 2 ⊗ C 2 ⊗ C 2 . Using the standard notation R ij (u) ∈ End((C 2 ) i ⊗ (C 2 ) j ), it reads R 12 (u/v)R 13 (u)R 23 (v) = R 23 (v)R 13 (u)R 12 (u/v) ∀u, v. (2.2)
Consider the reflection equation (also called the boundary quantum Yang-Baxter equation) introduced in the context of the boundary quantum inverse scattering theory (see [START_REF] Cherednik | Factorizing particles on the half-line and root systems[END_REF], [S88] for details). For the U q ( sl 2 ) R-matrix (2.1), the reflection equation reads:

R(u/v) (K(u) ⊗ II) R(uv) (II ⊗ K(v)) = (II ⊗ K(v)) R(uv) (K(u) ⊗ II) R(u/v) , (2.3)
where K(u) is a 2 × 2 matrix. We are now interested in a certain class of solutions of the reflection equation (2.3) for which the entries of the matrix K(u) are assumed to be Laurent polynomials in the spectral parameter whose coefficients are elements in certain homomorphic images of the q-Onsager algebra (see [START_REF] Baseilhac | A deformed analogue of Onsager's symmetry in the XXZ open spin chain[END_REF]). Here, we consider the simplest solution of this type. We refer the reader to [START_REF] Zabrodin | Quantum transfer matrices for discrete and continuous quasi-exactly solvable problems[END_REF][START_REF] Baseilhac | An integrable structure related with tridiagonal algebras[END_REF] for details.

Proposition 2.1. Assume {χ, ρ, ω, η, η * } are generic scalars in C * . Let K 11 (u) ≡ A(u), K 12 (u) ≡ B(u), K 21 (u) ≡ C(u), K 22 (u) ≡ D(u) be the elements of the 2 × 2 square matrix K(u) such that:

A(u) = (u 2 -u -2 ) qu A -q -1 u -1 A * -(q + q -1 )ρ -1 ηu + η * u -1 , (2.4) D(u) = (u 2 -u -2 ) qu A * -q -1 u -1 A -(q + q -1 )ρ -1 η * u + ηu -1 , (2.5) B(u) = χ(u 2 -u -2 ) ρ -1 [A * , A] q + ω q -q -1 + qu 2 + q -1 u -2 q 2 -q -2 , (2.6) C(u) = ρχ -1 (u 2 -u -2 ) ρ -1 [A, A * ] q + ω q -q -1 + qu 2 + q -1 u -2 q 2 -q -2 . (2.7)
Then K(u) satisfies the reflection equation (2.3) provided A, A * generate the Askey-Wilson algebra with defining relations (1.14),(1.15).

Note that the scalar parameter χ does not enter in the structure constants of the Askey-Wilson relations but it is introduced for further convenience.

Two explicit examples are now given. Let {q ±s3 , S ± } denote the generators of the quantum algebra U q (sl 2 ) with defining relations:

[s 3 , S ± ] = ±S ± , [S + , S -] = q 2s3 -q -2s3 q -q -1 . (2.8)
The Casimir element C = (qq -1 ) 2 S -S + + q 2s3+1 + q -2s3-1 (2.9) is central. Introduce the coproduct ∆ : U q (sl 2 ) → U q (sl 2 ) ⊗ U q (sl 2 ) such that:

∆(S + ) = S + ⊗ II + q 2s3 ⊗ S + ∆(S -) = S -⊗ q -2s3 + II ⊗ S -, ∆(q s3 ) = q s3 ⊗ q s3 . (2.10) Example 1. [GZ93a] Define the parameters k + = - 1 2 q ν q -q -1 , k -= 1 2 q ν ′ q -q -1 , ǫ + = cosh(µ)q 1 2 (ν+ν ′ ) , ǫ -= cosh (µ ′ ) q 1 2 (ν+ν ′ ) , (2.11)
where µ, µ ′ , ν, ν ′ , v are generic scalars. The following gives an example of elements A, A * satisfying the Askey-Wilson relations (1.14), (1.15):

A → k + vq 1/2 S + q s3 + k -v -1 q -1/2 S -q s3 + ǫ + q 2s3 , (2.12) A * → k + v -1 q -1/2 S + q -s3 + k -vq 1/2 S -q -s3 + ǫ -q -2s3 , (2.13)
where the structure constants ρ, ω, η, η * are given by (3.33)-(3.36) with the substitution q 2s+1 + q -2s-1 → C.

Let us mention that this example can be derived using the 'dressing' procedure [S88] that generates solutions of the reflection equation. For details, see [START_REF] Zabrodin | Quantum transfer matrices for discrete and continuous quasi-exactly solvable problems[END_REF].

Another example follows from [START_REF] Granovskii | Hidden symmetry of the Racah and Clebsch-Gordan problems for the quantum algebra slq(2)[END_REF] (see also [H16]). It will connect the Heun-Askey-Wilson element (1.16) to Hamiltonians of 3-sites spin chains, see Section 5.

Example 2. [START_REF] Granovskii | Hidden symmetry of the Racah and Clebsch-Gordan problems for the quantum algebra slq(2)[END_REF] (see also [H16]) The Askey-Wilson algebra is embedded into

U q (sl 2 ) ⊗ U q (sl 2 ) ⊗ U q (sl 2 ). Define c 1 = C ⊗ II ⊗ II, c 2 = II ⊗ C ⊗ II, c 3 = II ⊗ II ⊗ C and c 4 = (id × ∆) • ∆(C)
. The following gives an example of elements A, A * satisfying the Askey-Wilson relations (1.14), (1.15):

A → ∆(C) ⊗ II , (2.14) A * → II ⊗ ∆(C) , (2.15)
where the structure constants are given by:

ρ = -(q 2 -q -2 ) 2 , (2.16) ω = -(q -q -1 ) 2 (c 1 c 3 + c 2 c 4 ) , (2.17) η * = (q 2 -q -2 ) 2 (q + q -1 ) (c 1 c 2 + c 3 c 4 ) , (2.18) η = (q 2 -q -2 ) 2 (q + q -1 ) (c 2 c 3 + c 1 c 4 ) . (2.19)
Given a solution of the reflection equation, it is known that a generating function for mutually commuting quantities is provided by the so-called transfer matrix [S88]. In the present case, consider the most general scalar solution of the so-called "dual" reflection equation given by [START_REF] Vega | Boundary K-matrices for the XXX,XXZ and XYZ spin chains[END_REF][START_REF] Ghoshal | Boundary S-matrix and boundary states in two-dimensional integrable field theory[END_REF]:

K + (u) = quκ + q -1 u -1 κ * κ + (q 2 u 2 -q -2 u -2 ) κ -ρ(q 2 u 2 -q -2 u -2 ) quκ * + q -1 u -1 κ , (2.20)
where κ ± , κ, κ * are generic scalars in C. The transfer matrix reads t(u) = tr (K + (u)K(u)), where the trace is taken over the two-dimensional auxiliary space. Expanded in the spectral parameter u, for K(u) with (2.4)-(2.7) the transfer matrix produces the Heun-Askey-Wilson element (1.16):

t(u) = (q 2 u 2 -q -2 u -2 )(u 2 -u -2 ) κ A + κ * A * + κ + χ -1 [A, A * ] q + κ -χ [A * , A] q + F 0 (u) , (2.21)
where F 0 (u) is a scalar function 8 . It follows that the spectral problem for the Heun-Askey-Wilson operator associated with (1.16) and the transfer matrix (2.21) are identical. In the next section, we will use this connection to diagonalize the Heun-Askey-Wilson operator using the algebraic Bethe ansatz.

Let us mention that the known central element (the so-called Casimir element) of the Askey-Wilson algebra can be extracted from the so-called Sklyanin's quantum determinant Γ(u) given in [S88]. For the reflection equation algebra, recall that the quantum determinant is defined as: Γ(u) = tr P - 12 (K(u) ⊗ II) R(u 2 q)(II ⊗ K(uq)) , (2.22) 8 The expression of F 0 (u) is not needed in further analysis, so we omit its explicit expression.

where P - 12 = (1 -P )/2 with P = R(1)/(qq -1 ). As shown in [S88], Γ(u), (K(u)) ij = 0 for any i, j. Inserting K(u) with (2.4)-(2.7) into (2.22), one gets:

Γ(u) = (u 2 q 2 -u -2 q -2 ) 2(q -q -1 ) ∆(u) - 2ρ (q -q -1 ) , with ∆(u) = - 2(q -q -1 ) ρ u 2 -u -2 q 2 u 2 -q -2 u -2 (Γ + γ 0 (u))
where the explicit form of the scalar function γ 0 (u) is omitted for simplicity and

Γ = q 2 -q -2 q -1 AA * [A * , A] q -q -2 [A * , A] 2 q + q -2 ρA 2 + ρq 2 A * 2 + ω(AA * + A * A) (2.23) + η * (1 + q -2 )A + q 2 η(1 + q -2 )A * . It follows that Γ satisfies [Γ, A] = [Γ, A * ] = 0 (2.24)
i.e. Γ is central in the Askey-Wilson algebra. Note that this result provides an alternative derivation of the central element given in [Z91, eq. (1.3)]. Also, let us mention that Γ is a restriction of the expression computed in [START_REF] Baseilhac | An attractive basis for the q-Onsager algebra[END_REF] for the q-Onsager algebra. Indeed, the Askey-Wilson algebra can be viewed as a certain quotient of the q-Onsager algebra by the relations (1.14), (1.15).

To conclude this section, let us make few additional comments. As mentioned in the introduction, for the quantum harmonic oscillator each couple (H, X) and (H, iP) satisfies the simplified Askey-Wilson relations (1.2), (1.3). It is thus natural to ask for such triplet in the context of the Askey-Wilson algebra. Let A, A * satisfy the Askey-Wilson relations (1.14), (1.15). Define the element:

A ⋄ = i √ ρ [A * , A] q + iω (q -q -1 ) √ ρ . (2.25)
Then, it is easy to show that the couple (A, A ⋄ ) satisfies the Askey-Wilson relations:

A, A, A ⋄ q q -1 = ρ A ⋄ + i(q -q -1 ) √ ρ η A - i √ ρ (q -q -1 )
ωI , (2.26)

A ⋄ , A ⋄ , A q q -1 = ρ A + i(q -q -1 ) √ ρ η A ⋄ + η * I . (2.27)
Similar relations are easily derived for the couple (A * , A ⋄ ). In the literature, the triplet (A, A * , A ⋄ ) is related with the concept of Leonard triple. We refer the reader to [START_REF] Curtin | Modular Leonard triples[END_REF][START_REF] Huang | The classification of Leonard triples of Q-Racah type[END_REF] for details. Some properties of the element (2.25) will be used in Appendix C.

Let us also mention that the transformation from the pair (H, X) to the pair (H,iP) relating the Askey-Wilson relations (1.2), (1.3) generalizes to the q-deformed case as follows. Given (A * , A) satisfying (1.14), (1.15) with structure constants ρ, ω, η, η * , define a new pair of elements ( Ā * , Ā) by:

Ā * = A * , Ā = i √ ρ [A, A * ] q + iω (q -q -1 ) √ ρ . (2.28)
Then, by direct calculation one finds that ( Ā * , Ā) satisfy

Ā, Ā, Ā * q q -1 = ρ Ā * + ω Ā + ηI , (2.29) Ā * , Ā * , Ā q q -1 = ρ Ā + ω Ā * + η * I (2.30) with the structure constants ω = i(q -q -1 )η * √ ρ , η * = - iω √ ρ q -q -1 , η = η . (2.31)
In Section 4 and Appendix E we will give various examples of operators acting on an infinite dimensional vector space which provide realizations of the Askey-Wilson algebra (1.14), (1.15). In Appendix E, the invertible transformation (2.28) will be used.

Diagonalization of the Heun-Askey-Wilson operator via the algebraic Bethe ansatz

In this section, using (2.21) the Heun-Askey-Wilson operator associated with (1.16) is diagonalized on any irreducible finite dimensional representation of the Askey-Wilson algebra (1.14)-(1.15) for different choices of parameters κ, κ * , κ ± using the framework of the algebraic Bethe ansatz [START_REF] Cao | Exact solution of XXZ spin chain with unparallel boundary fields[END_REF][START_REF] Doikou | A Note on the boundary spin s XXZ chain[END_REF] and its modified version [Be15, BP15, ABGP15] combined with the theory of Leonard pairs [START_REF] Terwilliger | Leonard pairs and the q-Racah polynomials[END_REF][START_REF] Terwilliger | Leonard pairs and the Askey-Wilson relations[END_REF]. The eigenstates are obtained in the form of Bethe states, and the spectrum of the Heun-Askey-Wilson operator is given in terms of solutions (the so-called Bethe roots) of Bethe equations. For each choice of the parameters, it is shown that the corresponding Bethe ansatz equations admit an alternative presentation in terms of the 'symmetrized' variables U i given by (3.103). Also, up to an overall factor the Bethe states are expressed in the Poincaré-Birkhoff-Witt basis of the Askey-Wilson algebra with polynomial coefficients in U i . Two eigenbases for Leonard pairs in terms of Bethe states are derived in this framework.

3.1. Preliminaries.

3.1.1. Leonard pairs. Let (π, V) denote a finite dimensional representation on which the elements of the Askey-Wilson algebra A, A * act as π(A), π(A * ). Assume π(A), π(A * ) are diagonalizable on V, the spectra are multiplicityfree and V is irreducible. Then, by [TV03, Theorem 6.2] the operators π(A), π(A * ) form a Leonard pair, see [TV03, Definition 1.1]. Define dim( V) = 2s + 1 with s an integer or half-integer. Given the eigenvalue sequence θ 0 , θ 1 , ..., θ 2s associated with π(A) (resp. the eigenvalue sequence θ * 0 , θ (i) in the eigenbasis of π(A), then π(A * ) acts as a tridiagonal matrix;

(ii) in the eigenbasis of π(A * ), then π(A) acts as a tridiagonal matrix. According to the properties (i)-(ii) the action of the Leonard pair π(A), π(A * ) on the eigenvectors takes the form:

π(A)|θ M = θ M |θ M , π(A * )|θ M = a M,M+1 |θ M+1 + a M,M |θ M + a M,M-1 |θ M-1 , (3.1) π(A * )|θ * M = θ * M |θ * M , π(A)|θ * M = a * M,M+1 |θ * M+1 + a * M,M |θ * M + a * M,M-1 |θ * M-1 , (3.2)
where the coefficients a 0,-1 = a 2s,2s+1 = a * 0,-1 = a * 2s,2s+1 = 0, and the set of coefficients {a M,M±1 }, {a M,M } and {a * M,M±1 }, {a * M,M } are determined according to the representation chosen. Furthermore, the eigenvalue sequences are such that (see [T99, Theorem 4.4 (case I)]):

θ M = bq 2M + cq -2M , θ * M = b * q 2M + c * q -2M , (3.3)
where b, c, b * , c * are generic scalars. Due to the existence of the Askey-Wilson relations (2.26), (2.27), on the vector space V it is clear that (A, A ⋄ ) and (A * , A ⋄ ) both give examples of Leonard pairs. Let |θ ⋄ 0 , ..., |θ ⋄ 2s denote the eigenvectors of π(A ⋄ ) with respective eigenvalues θ ⋄ 0 , ..., θ ⋄ 2s . Adapting the results of [H12], one finds that the eigenvalues of A ⋄ take the form:

θ ⋄ M = b ⋄ q 2M + c ⋄ q -2M . (3.4)
where b ⋄ , c ⋄ are generic scalars. In this parametrization, note that the structure constant ρ is given by [T99, Lemma 4.5]:

ρ = -bc(q 2 -q -2 ) 2 = -b * c * (q 2 -q -2 ) 2 = -b ⋄ c ⋄ (q 2 -q -2 ) 2 . (3.5)
On the finite dimensional vector space V, in addition to the Askey-Wilson relations (1.14), (1.15) by the Cayley-Hamilton theorem one has:

2s M=0 (A -θ M ) = 0 , 2s M=0 (A * -θ * M ) = 0 , 2s M=0 (A ⋄ -θ ⋄ M ) = 0 . (3.6) Example 3. [GZ93a]
With respect to Example 1, let (π, V (s)) denote the irreducible (spin-s) representation of U q (sl 2 ) of dimension dim(V (s)) = 2s + 1 on which the generators {q ±s3 , S ± } act as:

π(q ±s3 ) = 2s+1 k=1 q ±(s+1-k) E kk , π(S -) = 2s+1 k=1 [k] q [2s + 1 -k] q E k+1k , π(S + ) = 2s+1 k=1 [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. For generic scalar parameters k ± , ǫ ± , v and q in (2.12), (2.13), then π(A), π(A * ) is a Leonard pair. There exists two bases on which π(A), π(A * ) act as (3.1), (3.2) and the spectra take the form (3.3) with (3.32). See e.g. [START_REF] Baseilhac | The q-Onsager algebra and multivariable q-special functions[END_REF] for details.

Example 4. [START_REF] Granovskii | Hidden symmetry of the Racah and Clebsch-Gordan problems for the quantum algebra slq(2)[END_REF][START_REF] Huang | An embedding of the universal Askey-Wilson algebra into Uq(sl 2 ) ⊗ Uq(sl 2 ) ⊗ Uq(sl 2 )[END_REF] With respect to Example 2, consider the tensor product of three irreducible representations of U q (sl 2 ) denoted V (j 1 ) ⊗ V (j 2 ) ⊗ V (j 3 ). Following [H16], define Σ the set consisting of all pairs (ℓ, k) of integers such that

0 ≤ ℓ ≤ min{2j 1 + 2j 2 , 2j 2 + 2j 3 , 2j 1 + 2j 3 , j 1 + j 2 + j 3 } , 0 ≤ k ≤ 2j 1 + 2j 2 + 2j 3 -2ℓ . (3.7)
Define V k (ℓ) as the vector subspace of V (j 1 ) ⊗ V (j 2 ) ⊗ V (j 3 ) spanned by the simultaneous eigenvectors of (id×∆)•∆(q 2s3 ) and c 4 = (id×∆)•∆(C) associated with eigenvalues q 2j1+2j2+2j3-2(k+ℓ) and q 2j1+2j2+2j3-2ℓ+1 + q -2j1-2j2-2j3+2ℓ-1 respectively. Then

V (j 1 ) ⊗ V (j 2 ) ⊗ V (j 3 ) = (ℓ,k)∈Σ V k (ℓ) . (3.8) The irreducible finite dimensional representation (π, V k (ℓ)) is such that dim(V k (ℓ)) = min(2j 1 , ℓ) + min(2j 2 , ℓ) + min(2j 3 , ℓ) -2ℓ + 1 . (3.9)
For all (ℓ, k) ∈ Σ, π(A), π(A * ) is a Leonard pair. There exists two bases on which π(A), π(A * ) act as (3.1), (3.2) and the spectra take the form (3.3) with (3.32) with the identification

b = q -2 min(2j3,ℓ)-2j1-2j2-1+2ℓ , c = b -1 and M = h -ℓ + min(2j 3 , ℓ) , (3.10) b * = q -2 min(2j1,ℓ)-2j2-2j3-1+2ℓ , c * = b * -1 and M = h -ℓ + min(2j 1 , ℓ)
where h is an integer such that

ℓ -min(2j 3 , ℓ) ≤ h ≤ min(2j 1 , ℓ) + min(2j 2 , ℓ) -ℓ. (3.11)
See e.g. [H16, Theorem 5.4, Lemma 5.6] for details.

3.1.2. Gauge transformation and reference states. Our aim is to diagonalize the Heun-Askey-Wilson operator associated with (1.16) on a finite dimensional vector space within the framework of the algebraic Bethe ansatz and using the theory of Leonard pairs. In this context, the first step is the construction of the so-called reference (or vacuum) state. Consider the K-matrix with non-diagonal entries of the form (2.6), (2.7). As A, A * act as a Leonard pair on V, it is not possible to construct a state that is annihilated by the off-diagonal entries. Indeed, given the representation V on which A, A * act as π(A), π(A * ), according to (i), (ii) there is no state |Ω such that π(C(u))|Ω = 0 (and similarly for π(B(u))) . To circumvent this problem, an idea [START_REF] Cao | Exact solution of XXZ spin chain with unparallel boundary fields[END_REF] is to apply a gauge transformation (parameterized by an integer m) to the K-matrices K(u) with (2.4)-(2.7) and K + (u) given by (2.20):

t(u) = tr K + (u)K(u) = tr K+ (u|m)K(u|m) , (3.12)
such that the off-diagonal entries of K(u|m) admit a reference state.

The gauge transformation is built as follows. For more details, we refer the reader to [START_REF] Cao | Exact solution of XXZ spin chain with unparallel boundary fields[END_REF] and [START_REF] Belliard | Modified algebraic Bethe ansatz for XXZ chain on the segment -II -general cases[END_REF] for the notations used here. Let ǫ = ±1, α, β be generic complex parameters and m be an integer. Introduce the covariant vectors

|X ǫ (u, m) = αq ǫm u ǫ 1 , |Y ǫ (u, m) = βq -ǫm u ǫ 1 (3.13)
and the contravariant vectors

Xǫ (u, m)| = -ǫ q -ǫ u -ǫ γ ǫ (1, m -1) -1 αq ǫm u ǫ , Ỹ ǫ (u, m)| = -ǫ q -ǫ u -ǫ γ ǫ (1, m + 1) 1 -βq -ǫm u ǫ (3.14) where γ ǫ (u, m) = α 1-ǫ 2 β ǫ+1 2 q -m u -α ǫ+1 2 β 1-ǫ 2 q m u -1 . (3.15)
Applying the gauge transformation to K(u), one finds:

K(u|m) = A ǫ (u, m) B ǫ (u, m) C ǫ (u, m) D ǫ (u, m) , with A ǫ (u, m) = Ỹ ǫ (u, m -2)|K(u)|X ǫ (u -1 , m) , (3.16) B ǫ (u, m) = Ỹ ǫ (u, m)|K(u)|Y ǫ (u -1 , m) , (3.17) C ǫ (u, m) = Xǫ (u, m)|K(u)|X ǫ (u -1 , m) , (3.18) D ǫ (u, m) = γ ǫ (1, m + 1) γ ǫ (1, m) Xǫ (u, m + 2)|K(u)|Y ǫ (u -1 , m) - (q -q -1 )γ ǫ (u -2 , m + 1) (qu 2 -q -1 u -2 )γ ǫ (1, m) A ǫ (u, m). (3.19)
In the literature, the elements above are called the 'dynamical operators'. Their explicit expressions in terms of the elements A, A * are reported in Appendix A. In terms of the dynamical operators, the transfer matrix (3.12) reads:

t(u) = (q 2 u 2 -q -2 u -2 ) a(u, m)A ǫ (u, m) + d(u, m)D ǫ (u, m) + b(u, m)B ǫ (u, m) + c(u, m)C ǫ (u, m) (3.20) where a(u, m) = αu 2ǫ κu + κ * u -1 + u ǫ (u 2 -u -2 )q -(m+1)ǫ (κ + -αβκ -ρ) -βq -(2m+2)ǫ κ * u + κu -1 (α -βq -ǫ(2m+2) )(qu 2 -q -1 u -2 ) , d(u, m) = -βu 2ǫ q -(2m+1)ǫ qκu + q -1 κ * u -1 -u ǫ q -(m+1)ǫ q 2 u 2 -q -2 u -2 (κ + -αβκ -ρ) + αq -ǫ qκ * u + q -1 κu -1 (α -βq -ǫ(2m+2) )(q 2 u 2 -q -2 u -2 ) , b(u, m) = u ǫ α 2 κ -ρq (m+2)ǫ -αǫq ǫ κ ǫ+1 2 κ * 1-ǫ 2 -κ + q -mǫ α -βq -2mǫ , c(u, m) = u ǫ -β 2 κ -ρq (2-3m)ǫ + βǫq (1-2m)ǫ κ ǫ+1 2 κ * 1-ǫ 2 + κ + q -mǫ α -βq -2mǫ .
According to a certain gauge transformation parametrized by the scalars α, β and ǫ = ±1, a reference state and its dual can be identified using the theory of Leonard pairs. Recall that the Leonard pair is formed by the operators π(A), π(A * ).

Lemma 3.1. Let m 0 be an integer. If the parameters α, β are such that:

(q 2 -q -2 )χ -1 αc * q m0 = 1 (resp. (q 2 -q -2 )χ -1 βc * q -m0 = 1) (3.21) then |Ω + ≡ |θ * 0 satisfies π(C + (u, m 0 ))|Ω + = 0 (resp. π(B + (u, m 0 ))|Ω + = 0 ). (3.22)
Proof. We show the first relation in (3.22). By (3.18), C + (u, m 0 ) is given by (A.3). By (3.2), one has:

π αq m0 uχ [A, A * ] q - χq -m0 αρu [A * , A] q - q 2 + 1 qu A + 1 qu 3 + qu A * + g 0 (u)I |θ * 0 (3.23) = u -1 const 1 (u)|θ * 1 + u -1 const 0 (u)|θ * 0
where we denote

g 0 (u) = q -m0-3 αρχ(q 2 -q -2 ) u -1 ρq 3 (qu 2 + q -1 u -2 ) α 2 ρq 2m0 -χ 2 +u -1 q 2 + 1 (α 2 ρq 2m0 -χ 2 )q 2 ω -αη * χq m0 (q 4 -1)
and

const 1 (u) = αq m0 χ (qθ * 0 -q -1 θ * 1 ) - χq -m0 αρ (qθ * 1 -q -1 θ * 0 ) -q + q -1 a * 0,1 , (3.24) const 0 (u) = αq m0 χ - χq -m0 αρ (q -q -1 )θ * 0 -(q + q -1 ) a * 0,0 + (qu 2 + q -1 u -2 )θ * 0 + ug 0 (u). (3.25)
Requiring π(C + (u, m 0 ))|Ω + = 0, it implies the conditions const 1 (u) = 0 and const 0 (u) = 0 to be satisfied. First, let us consider const 1 (u) = 0. As a * 0,1 = 0 in (3.2), using (3.3) for M = 0, 1, the equation const 1 (u) = 0 gives a quadratic equation in α, which admits two solutions:

α = χq -m0 (q 2 -q -2 )c * or α ′ = χq -m0+2 (q 2 -q -2 )c * . (3.26)
We now turn to const 0 (u) = 0. It is found that only the first solution α is such that the coefficients of u 2 and u -2 in const 0 (u) are vanishing. For this choice of α, it follows:

const 0 (u) = 0 ⇔ (q -q -1 ) 2 θ * 0 2 + ρ a * 0,0 + θ * 0 ω + η * = 0 . (3.27)
To show that the coefficient a * 0,0 satisfies the second equation, we use the defining relations of the Askey-Wilson algebra. Indeed, by (1.15) one has: Similarly, the following lemma is shown. The proof follows the same steps as previously, so we omit it.

π A * , A * , A q q -1 -ρ A -ω A * -η * I |θ * 0 = 0 . ( 3
Lemma 3.2. Let m 0 be an integer. If the parameters α, β are such that:

(q 2 -q -2 )χ -1 αbq -m0 = -1 (resp. (q 2 -q -2 )χ -1 βbq m0 = -1) (3.29) then |Ω -≡ |θ 0 satisfies π(C -(u, m 0 ))|Ω -= 0 (resp. π(B -(u, m 0 ))|Ω -= 0 ). (3.30)
For a proper choice of the gauge parameters α, β, the action of the 'diagonal' dynamical operators A ± (u, m 0 ) and D ± (u, m 0 ) on the reference states |Ω ± is easily computed. Lemma 3.3. Let α be fixed by (3.21) for |Ω + or fixed by (3.29) for |Ω -. Then, the dynamical operators are such that:

π(A ± (u, m 0 ))|Ω ± = Λ ± 1 (u)|Ω ± and π(D ± (u, m 0 ))|Ω ± = Λ ± 2 (u)|Ω ± , (3.31)
where Λ ± 1 (u), Λ ± 2 (u) are ratios of Laurent polynomials in the variable u. Proof. We show π(A + (u, m 0 ))|Ω + = Λ + 1 (u)|Ω + . By (3.16), A + (u, m 0 ) is given by (A.1). The action of A + (u, m 0 ) on the reference state

|Ω + ≡ |θ * 0 reads π(A + (u, m 0 ))|Ω + = u 2 -u -2 u(αq 2m0 -q 2 β) π χq m0 ρ [A * , A] q - αβq m0+2 χ [A, A * ] q + β + αq 2m0 qA - αq 2m0 + βq 4 u 4 qu 2 A * + f 0 (u)I |θ * 0
where we have denoted

f 0 (u) = - 1 ρq 3 u 2 χ (q 2 -q -2 ) (u 2 -u -2 ) ρq m0+2 (αβq 2 ρ -χ 2 )(u 6 q 2 -u -2 ) +(q 2 + 1)q m0 αη * χq m0 (q 4 -1) -q 2 ω(αβq 2 ρ -χ 2 ) -(q 2 + 1)q 2 u 4 βη * χ(q 4 -1) -q m0 ω(αβq 2 ρ -χ 2 ) +(q 2 -1)u 2 (αηχq 2m0 -βηχq 2 )(1 + q 2 ) 2 -ρq m0+2 (αβq 2 ρ -χ 2 ) .
For the choice of gauge parameter α fixed by (3.21), on one hand using (3.2) one finds:

π χq m0 ρ [A * , A] q - αβq m0+2 χ [A, A * ] q |θ * 0 = χq m0 + βb * q 2 (q 2 -q -2 ) ρ (q -q -1 )θ * 0 a * 0,0 |θ * 0 + χq m0+1 b * (q 2 -q -2 ) ρ -βq a * 0,1 |θ * 1 .
On the other hand, using (3.2) one finds:

π β + αq 2m0 qA |θ * 0 = βq -χq m0+1 b * (q 2 -q -2 ) ρ a * 0,0 |θ * 0 - χq m0+1 b * (q 2 -q -2 ) ρ -βq a * 0,1 |θ * 1 .
Combining both expressions, the off-diagonal contribution in |θ * 1 vanishes. It follows:

π(A + (u, m 0 ))|Ω + = (u 2 -u -2 )(q 2 -q -2 )c * u(q m0 χ -c * β(q 4 -1)) χq m0 + βb * q 2 (q 2 -q -2 ) ρ (q -q -1 )θ * 0 + βq -χq m0+1 b * (q 2 -q -2 ) ρ a * 0,0 - χq m0 c * (q 2 -q -2 ) q -2 u -2 + βq 2 u 2 qθ * 0 + f 0 (u) |Ω + .
Inserting the expression of a * 0,0 according to (3.27) and using (3.3), Λ + 1 (u) is obtained as a ratio of Laurent polynomials in u, depending on b * , c * and the structure constant ρ, ω, η * . The other polynomial expressions Λ - 1 (u), Λ ± 2 (u) are obtained following the same steps, so we omit the details. Without loss of generality, for a convenient choice of parametrization the eigenvalues Λ ± 1 (u), Λ ± 2 (u) can be written in a factorized form. The proof of the following result is straightforward. Let µ, µ ′ , ν, ν ′ , v be generic scalars. Introduce the parametrization:

b = 1 2 q 1 2 (ν+ν ′ ) e -µ q -2s , c = 1 2 q 1 2 (ν+ν ′ ) e µ q 2s , (3.32) b * = 1 2 q 1 2 (ν+ν ′ ) e µ ′ q -2s , c * = 1 2 q 1 2 (ν+ν ′ ) e -µ ′ q 2s , b ⋄ = 1 2 q 1 2 (ν+ν ′ ) v 2 q -2s , c ⋄ = 1 2 q 1 2 (ν+ν ′ ) v -2 q 2s ,
and define the parameter χ as

χ = - 1 2 q ν (q 2 -q -2 ) .
Adapting the results of [T02a, T02b] the structure constants are given by:

ρ = - 1 4 q ν+ν ′ (q 2 -q -2 ) 2 , (3.33) ω = 1 4 q ν+ν ′ (q -q -1 ) 2 (v 2 + v -2 )(q 2s+1 + q -2s-1 ) -4 cosh(µ) cosh(µ ′ ) , (3.34) η = - 1 4 q 3 2 (ν+ν ′ ) (q 2 -q -2 ) 2 (q + q -1 ) cosh(µ)(v 2 + v -2 ) -cosh(µ ′ )(q 2s+1 + q -2s-1 ) , (3.35) η * = - 1 4 q 3 2 (ν+ν ′ ) (q 2 -q -2 ) 2 (q + q -1 ) cosh(µ ′ )(v 2 + v -2 ) -cosh(µ)(q 2s+1 + q -2s-1 ) . (3.36)
Then, the eigenvalues of the dynamical operators take the factorized form:

Λ ǫ 1 (u) = 1 2u ǫ e ǫ(µ-µ ′ )/2-(µ+µ ′ )/2 q (ν+ν ′ )/2-(2s+1) q 2s+1 uv -1 -u -1 v q 2s+1 uv -u -1 v -1 (3.37) × ue 1 2 (1-ǫ)(µ-µ ′ ) + u -1 e 1 2 (ǫ+1)(µ ′ +µ) ue 1 2 (1-ǫ)(µ ′ +µ) + u -1 e 1 2 (ǫ+1)(µ ′ -µ) , Λ ǫ 2 (u) = (u 2 -u -2 ) 2u ǫ (qu 2 -q -1 u -2 ) e ǫ(µ-µ ′ )/2-(µ+µ ′ )/2 q (ν+ν ′ )/2-(2s+1) q 2s-1 u -1 v -uv -1 q 2s-1 u -1 v -1 -uv × q 2 ue 1 2 (ǫ+1)(µ ′ -µ) + u -1 e 1 2 (1-ǫ)(µ ′ +µ) q 2 ue 1 2 (ǫ+1)(µ ′ +µ) + u -1 e 1 2 (1-ǫ)(µ-µ ′ ) .
To summarize, given a Leonard pair there are at least9 two possible choices of gauge transformation for which a reference state can be identified. In each case, the reference state is simply given by the fundamental eigenvector of either π(A) or π(A * ).

3.1.3. Eigenvalue of the Casimir Γ. On the vector space generated by successive actions of the dynamical operators {A ǫ (u, m), B ǫ (u, m), C ǫ (u, m), D ǫ (u, m)} on a reference state (the so-called Bethe states), the eigenvalue of the central element Γ is computed as follows. The dynamical operators being combinations of A, A * , they commute with Γ. So, it is sufficient to compute the action of Γ on a reference state, for instance |Ω + . Recall |Ω + ≡ |θ * 0 . Applying (2.23) and using (3.2), one gets:

Γ|Ω + = Γ 0 |Ω + + a * 01 Γ 1 |θ * 1 + a * 01 a * 12 Γ 2 |θ * 2 , (3.38)
where Γ 0 ( ≡ Γ 0 (a * 00 , a * 01 a * 10 )), Γ 1 (≡ Γ 1 (a * 00 , a * 11 )) and Γ 2 are functions of the spectral data θ * 0 , θ * 1 , θ * 2 and the structure constants ρ, ω, η, η * . Using (3.3), (3.5), one finds Γ 2 = 0. Applying (1.14), (1.15), to |θ * 0 , one extracts the expressions of a * 00 , a * 11 as well as the product a * 01 a * 10 in terms of the ρ, ω, η, η * . Inserted in Γ 1 (a * 00 , a * 11 ), one finds Γ 1 (a * 00 , a * 11 ) = 0. Similarly, inserted in the eigenvalue Γ 0 ≡ Γ 0 (a * 00 , a * 01 a * 10 ), using the parametrization (3.32) and (3.33)-(3.36) one gets finally:

Γ|Ω + = Γ 0 |Ω + (3.39) with Γ 0 = q -q -1 2 q 2(ν+ν ′ -2s) q 4s 4 q 2 + q -2 -(q -2s-1 + q 2s+1 )(v 2 + v -2 ) cosh(µ) cosh(µ ′ ) + 1 2 (cosh(2µ) + cosh(2µ ′ )) -cosh(2µ) cosh(2µ ′ ) - 1 16 q 4s -1 q 2 q 4s -q -2 (v 4 + v -4 ) + 1 16
1 + q 4 q 8s + -q -4 + 4q -2 + 2 + 4q 2q 4 q 4s + q -4 + 1 .

3.1.4. Commutation relations between the dynamical operators. In the following analysis, we will need the commutation relations between the entries of K(u|m). These relations are derived from the reflection equation (2.3). By straightforward calculations, in particular one has:

B ǫ (u, m + 2)B ǫ (v, m) = B ǫ (v, m + 2)B ǫ (u, m), (3.40) A ǫ (u, m + 2)B ǫ (v, m) = f (u, v)B ǫ (v, m)A ǫ (u, m) (3.41) +g(u, v, m)B ǫ (u, m)A ǫ (v, m) + w(u, v, m)B ǫ (u, m)D ǫ (v, m), D ǫ (u, m + 2)B ǫ (v, m) = h(u, v)B ǫ (v, m)D ǫ (u, m), (3.42) +k(u, v, m)B ǫ (u, m)D ǫ (v, m) + n(u, v, m)B ǫ (u, m)A ǫ (v, m), C (u, m + 2)B(v, m) = B(v, m -2)C (u, m) (3.43) +q(u, v, m)A (v, m)D(u, m) + r(u, v, m)A (u, m)D(v, m) +s(u, v, m)A (u, m)A (v, m) + x(u, v, m)A (v, m)A (u, m) +y(u, v, m)D(u, m)A (v, m) + z(u, v, m)D(u, m)D(v, m),
where the coefficients are collected in appendix B. Using these commutation relations, the action of the entries on products of off-diagonal operators can be derived. Let ū denote the set of variables ū = {u 1 , u 2 , . . . , u M }. Define the following strings of length M of operators B ǫ (u i , m):

B ǫ (ū, m, M ) = B ǫ (u 1 , m + 2(M -1)) • • • B ǫ (u M , m) , (3.44) B ǫ ({u, ūi }, m, M ) = B ǫ (u 1 , m + 2(M -1)) • • • B ǫ (u, m + 2(M -i)) . . . B ǫ (u M , m) . (3.45)
Using the dynamical commutation relations (3.40), (3.41), (3.42), one shows that the action of the diagonal dynamical operators {A ǫ (u, m), D ǫ (u, m)} on the string (3.44) is given by

A ǫ (u, m + 2M )B ǫ (ū, m, M ) = M i=1 f (u, u i )B ǫ (ū, m, M )A ǫ (u, m) (3.46) + M i=1 g(u, u i , m + 2(M -1)) M j=1,j =i f (u i , u j )B ǫ ({u, ūi }, m, M )A ǫ (u i , m) + M i=1 w(u, u i , m + 2(M -1)) M j=1,j =i h(u i , u j )B ǫ ({u, ūi }, m, M )D ǫ (u i , m) and D ǫ (u, m + 2M )B ǫ (ū, m, M ) = M i=1 h(u, u i )B ǫ (ū, m, M )D ǫ (u, m) (3.47) + M i=1 k(u, u i , m + 2(M -1)) M j=1,j =i h(u i , u j )B ǫ ({u, ūi }, m, M )D ǫ (u i , m) + M i=1 n(u, u i , m + 2(M -1)) M j=1,j =i f (u i , u j )B ǫ ({u, ūi }, m, M )A ǫ (u i , m) .
In the next subsections, we will use extensively the above relations for various choices of gauge parameters ǫ, α, β in order to solve the spectral problem for the Heun-Askey-Wilson operator associated with (1.16).

3.1.5. Action of π(B ǫ (u, m + 4s)) on V. According to the choice of the gauge parameter α, recall that the action of the dynamical operators on each reference state |Ω ± is given in Lemmas 3.1, 3.2, 3.3. The vector space V being finite dimensional, we now study the action of the dynamical operator B + (u, m 0 + 4s) (resp. B -(u, m 0 + 4s)) on the eigenvector |θ * 2s (resp. |θ 2s ). Recall (3.1), (3.2).

Lemma 3.4. Assume (q 2q -2 )χ -1 βb * q -m0 = 1 (resp. (q 2q -2 )χ -1 βcq m0 = -1). Then

π(B + (u, m 0 + 4s))|θ * 2s = 0 (resp. π(B -(u, m 0 + 4s))|θ 2s = 0) . (3.48)
Proof. We show the first relation, the second one is derived similarly. Using (A.2), acting on the eigenstate |θ * j of A * with eigenvalue θ * j one finds:

π(B + (u, m))|θ * j = β(u 2 -u -2 ) (αq 2m+2 -β)u B jj+1 (m)|θ * j+1 + B jj (U, m)|θ * j + B jj-1 (m)|θ * j-1 (3.49) where B 0-1 = B 2s2s+1 = 0, B jj+1 (m) = χq m βρ (qθ * j+1 -q -1 θ * j ) - βq -m χ (qθ * j -q -1 θ * j+1 ) + (q + q -1 ) a * jj+1 , B jj (U, m) = χq m βρ - βq -m χ (q -q -1 )θ * j + (q + q -1 ) a * jj -U (q + q -1 )θ * j + (χ 2 q m -β 2 ρq -m ) βχ(q -q -1 ) U + ω ρ + η * (q + q -1 ) ρ , B jj-1 (m) = χq m βρ (qθ * j-1 -q -1 θ * j ) - βq -m χ (qθ * j -q -1 θ * j-1 ) + (q + q -1 ) a * jj-1 ,
and we have denoted U = (qu 2 + q -1 u -2 )/(q + q -1 ). From the Askey-Wilson relations, one gets:

a * jj = - ωθ * j + η * ρ + (q -q -1 ) 2 (θ * j ) 2 .
(3.50)

For j = 2s it follows:

π(B + (u, m))|θ * 2s = β(u 2 -u -2 ) (αq 2m+2 -β)u B 2s2s (U, m)|θ * 2s + B 2s2s-1 (m)|θ * 2s-1 .
In particular, for m = m 0 + 4s and β such that (q 2q -2 )χ -1 βb * q -m0 = 1 one finds:

B 2s2s (U, m 0 + 4s) = 0 , B 2s2s-1 (m 0 + 4s) = 0 (3.51)
which implies the first relation in (3.48).

3.2. Diagonalization of the 'special' case κ * = κ ± = 0 or κ = κ ± = 0. We start by considering two simple specializations of the Heun-Askey-Wilson operator (1.16), namely:

I(κ, 0, 0, 0) = κ A or I(0, κ * , 0, 0) = κ * A * . (3.52)
In these cases, the spectral problem for the specialization of the Heun-Askey-Wilson operator acting on an irreducible finite dimensional vector space reduces to the spectral problem for either π(A) or π(A * ). Below, the corresponding eigenstates are written as Bethe states, and the eigenvalues are derived. As a byproduct, the two families of Bethe states provide two explicit bases for the Leonard pair π(A), π(A * ), see subsection 3.5.3 .

In the framework of the algebraic Bethe ansatz, we first need to write the two elements A and A * in terms of the dynamical operators (3.16)-(3.19) according to the reference state on which these operators act, either |Ω -or |Ω + . From (2.4), (2.5) and using (3.16), (3.19), according to the choice of gauge transformation one gets for instance:

A = A -(u, m) + q u η(u) + q -1 u -1 η(u -1 ) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) , (3.53) A * = A + (u, m) + q u η(u -1 ) + q -1 u -1 η(u) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) with η(u) = (q + q -1 )ρ -1 ηu + η * u -1 , (3.54)
where, for further convenience, we have introduced:

A -(u, m) = u -1 (u 2 -u -2 ) 1 (qu 2 -q -1 u -2 ) A -(u, m) + 1 (q 2 u 2 -q -2 u -2 ) D -(u, m) , (3.55) A + (u, m) = u (u 2 -u -2 ) 1 (qu 2 -q -1 u -2 ) A + (u, m) + 1 (q 2 u 2 -q -2 u -2 ) D + (u, m) . (3.56)
The spectral problem is now solved. Choose the gauge parameter α according to Lemma 3.1 or Lemma 3.

π (I(κ, 0, 0, 0)) |Ψ M sp,-(ū, m 0 ) = κ 2 q 1 2 (ν+ν ′ ) e -µ q -2s+2M + e µ q 2s-2M |Ψ M sp,-(ū, m 0 ) (3.58)
where the set ū satisfies the Bethe equations:

M j=1,j =i b(u i /(qu j ))b(u i u j ) b(qu i /u j )b(q 2 u i u j ) = qe µ ′ u i + q -1 e µ u -1 i qe -µ u i + q -1 e µ ′ u -1 i b q 1 2 -s vu i b q 1 2 -s v -1 u i e µ ′ u i + e -µ u -1 i e µ u i + e µ ′ u -1 i b q s+ 1 2 vu i b q s+ 1 2 v -1 u i for i = 1, . . . , M .
Proof. Recall the notation (1.21). Using the (off-shell) relations (3.46), (3.47) and the actions (3.31), one finds:

π(A -(u, m 0 + 2M ))|Ψ M sp,-(ū, m 0 ) = Λ M -(u, ū)|Ψ M sp,-(ū, m 0 ) (3.59) - u -1 b(q) b(u 2 ) M i=1 E M sp (u i , ūi ) b(uu -1 i )b(quu i ) |Ψ M sp,-({u, ūi }, m 0 )
where we denote

|Ψ M sp,-({u, ūi }, m 0 ) = π(B ǫ ({u, ūi }, m 0 , M ))|Ω -, Λ M -(u, ū) = u -1 (u 2 -u -2 )   1 (qu 2 -q -1 u -2 ) M j=1 f (u, u j )Λ - 1 (u) + 1 (q 2 u 2 -q -2 u -2 ) M j=1 h(u, u j )Λ - 2 (u)   (3.60)
with (3.37) and

E M sp (u i , ūi ) = - b(u 2 i ) b(qu 2 i ) M j=1,j =i f (u i , u j )Λ - 1 (u i ) + M j=1,j =i h(u i , u j )Λ - 2 (u i ) (3.61) for i = 1, . . . , M . Requiring E M
sp (u i , ūi ) = 0, one gets the so-called Bethe equations given below (3.58) using (3.37). To determine the eigenvalues in the r.h.s of (3.58), we proceed as follows. From (3.60), observe that

Λ M -(u, ū) + q u η(u) + q -1 u -1 η(u -1 ) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) (3.62)
is a meromorphic function in the variable u. To be equal to a constant (i.e. independent of u), we need to study the singular part of this function. Using (3.37) and the expressions in Appendix B, one finds that the singular points are located at:

u ∈ {±u j , ±q -1 u -1 j , j = 1, ..., M } . (3.63)
The sum of all the residues at these points vanishes. Extracting the constant part of (3.62), one obtains (3.58).

Similarly, the eigenstates of the Heun-Askey-Wilson operator for the special case κ = κ + = κ -= 0 are derived. The proof being analog to the previous case, we skip the details.

Proposition 3.2. Define |Ψ M sp,+ (ū, m 0 ) = π(B + (ū, m 0 , M ))|Ω + . (3.64)
One has:

π (I(0, κ * , 0, 0)) |Ψ M sp,+ (ū, m 0 ) = κ * 2 q 1 2 (ν+ν ′ ) e -µ ′ q 2s-2M + e µ ′ q -2s+2M |Ψ M sp,+ (ū, m 0 ) (3.65)
where the set ū satisfies the Bethe equations:

M j=1,j =i b(u i /(qu j ))b(u i u j ) b(qu i /u j )b(q 2 u i u j ) = qe -µ u i + q -1 e -µ ′ u -1 i qe µ ′ u i + q -1 e -µ u -1 i b q 1 2 -s vu i b q 1 2 -s v -1 u i e -µ u i + e µ ′ u -1 i e -µ ′ u i + e -µ u -1 i b q s+ 1 2 vu i b q s+ 1 2 v -1 u i for i = 1, . . . , M .
Let us point out that the spectra of the two specializations (3.52) do not depend on the Bethe roots. Furthermore, as expected the structure of the eigenvalues matches with (3.3) in agreement with [T99, Theorem 4.4 (case I)].

3.3. Diagonalization of the 'diagonal' case κ = 0, κ * = 0 and κ ± = 0. This case is associated with the diagonal form of the K + (u) matrix (2.20). The Heun-Askey-Wilson operator that we will diagonalize below using the algebraic Bethe ansatz is given by:

I(κ, κ * , 0, 0) = κ A + κ * A * . (3.66)
Similarly to the special case discussed in the previous subsection, the Heun-Askey-Wilson operator associated with (3.66) is first expressed in terms of the dynamical operators (3.16)-(3.19). If we choose the reference state to be |Ω -, according to the corresponding gauge transformation the element A is given by (3.53) whereas A * is given by:

A * = Ã-(u, m) + q u η(u -1 ) + q -1 u -1 η(u) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) (3.67) with Ã-(u, m) = u -1 (u 2 -u -2 ) γ -q -1 u -2 , m (qu 2 -q -1 u -2 )γ -(1, m + 1) A -(u, m) + γ -qu 2 , m (q 2 u 2 -q -2 u -2 )γ -(1, m + 1) D -(u, m) + αq -m-1 γ -(1, m) B -(u, m) - βq m-1 γ -(1, m) C -(u, m)
where the notation (3.15) is used. On the other hand, if the reference state is |Ω + , the element A * is given by (3.54) whereas A is now given by:

A = Ã+ (u, m) + q u η(u) + q -1 u -1 η(u -1 ) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) , (3.68) where Ã+ (u, m) = u (u 2 -u -2 ) γ + q -1 u -2 , m (qu 2 -q -1 u -2 )γ + (1, m + 1) A + (u, m) + γ + qu 2 , m (q 2 u 2 -q -2 u -2 )γ + (1, m + 1) D + (u, m) + αq m+1 γ + (1, m) B + (u, m) - βq 1-m γ + (1, m) C + (u, m) .
According to the gauge transformation chosen, combining the expressions for A, A * it follows:

κ A + κ * A * = u ǫ αu ǫ κu + κ * u -1 -βq -(2m+2)ǫ u -ǫ κ * u + κu -1 (u 2 -u -2 ) (qu 2 -q -1 u -2 ) α -βq -(2m+2)ǫ A ǫ (u, m) (3.69) + q -ǫ u ǫ αu -ǫ qκ * u + q -1 κu -1 -βq -2mǫ u ǫ qκu + q -1 κ * u -1 (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) α -βq -(2m+2)ǫ D ǫ (u, m) - ǫ αq ǫ κ ǫ+1 2 κ * 1-ǫ 2 u ǫ (u 2 -u -2 ) (α -βq -2mǫ ) B ǫ (u, m) + ǫ βκ ǫ+1 2 κ * 1-ǫ 2 q -(2m-1)ǫ u ǫ (u 2 -u -2 ) (α -βq -2mǫ ) C ǫ (u, m) + (q + q -1 ) 2 ρ (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) ηκ * + η * κ + (ηκ + η * κ * ) qu 2 + q -1 u -2 q + q -1 .
Recall that (3.12) holds for any choice of gauge transformation, i.e. any choice of parameters ǫ, α and β. Having fixed the gauge parameter α according to the choice of reference state |Ω + or |Ω -, to simplify the analysis of the spectral problem for (3.66), without loss of generality we choose to fix the gauge parameter β in order to eliminate the term C ǫ (u, m) in (3.69). To this end, we set:

β = 0 .
A crucial ingredient for the solution of the spectral problem of (3.66) is the following lemma10 . For simplicity, the proof is reported Appendix C. Recall the notation (3.44), (3.45). For a generic set of parameters ū = {u 1 , u 2 , ..., u 2s }, define the Bethe vector:

|Ψ 2s d,ǫ (ū, m 0 ) = π(B ǫ (ū, m 0 , 2s))|Ω ǫ , |Ψ 2s
d,ǫ ({u, ūi }, m 0 ) = π(B ǫ ({u, ūi }, m 0 , 2s))|Ω ǫ . Lemma 3.5. For M = 2s and generic {u, u i }, one has:

π(B ǫ (u, m 0 + 4s))|Ψ 2s d,ǫ (ū, m 0 ) = (3.70) δ d u -ǫ b(u 2 ) 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) 2s i=1 b(uu -1 i )b(q -1 u -1 u -1 i ) |Ψ 2s d,ǫ (ū, m 0 ) -δ d 2s i=1 u -ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) b(uu -1 i )b(q -1 u -1 u -1 i ) 2s j=1,j =i b(u i u -1 j )b(q -1 u -1 i u -1 j ) |Ψ 2s d,ǫ ({u, ūi }, m 0 )
where we denote

δ d = - ǫ(-1) 2s+1 2 e -µ(1-ǫ)/2-µ ′ (1+ǫ)/2 q (ν+ν ′ )/2-ǫ(2s+2) . (3.71)
We now turn to the solution of the spectral problem.

Proposition 3.3. For ǫ = ±1, one has:

π (I(κ, κ * , 0, 0)) |Ψ 2s d,ǫ (ū, m 0 ) = Λ 2s d,ǫ |Ψ 2s d,ǫ (ū, m 0 ) (3.72) with Λ 2s d,+ = κ * θ * 2s + κe µ-µ ′ b (v 2 + v -2 )[2s] q + 2e µ ′ cosh(µ) -q 2s j=1 (qu 2 j + q -1 u -2 j ) , (3.73) Λ 2s d,-= κθ 2s + κ * e µ ′ -µ c * (v 2 + v -2 )[2s] q + 2e µ cosh(µ ′ ) -q -1 2s j=1 (qu 2 j + q -1 u -2 j ) , (3.74)
where the set ū satisfies the (inhomogeneous) Bethe equations:

b(u 2 i ) b(qu 2 i ) (κ u i + κ * u -1 i ) 2s j=1,j =i f (u i , u j )Λ ǫ 1 (u i ) -q -ǫ u -2ǫ i (qκ * u i + q -1 κu -1 i ) 2s j=1,j =i h(u i , u j )Λ ǫ 2 (u i ) (3.75) +(-1) 2s ǫ(q -q -1 ) -1 q ǫ κ (1+ǫ)/2 κ * (1-ǫ)/2 δ d u -2ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) 2s j=1,j =i b(u i u -1 j )b(qu i u j ) = 0
for i = 1, . . . , 2s.

Proof. For convenience, define the element:

W d (u, m) = κ A + κ * A * - (q + q -1 ) 2 ρ (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) ηκ * + η * κ + (ηκ + η * κ * ) qu 2 + q -1 u -2 q + q -1 .
Explicitly, using (3.69) for β = 0 it gives:

W d (u, m) = u ǫ ∆ d (u) (u 2 -u -2 ) (qu 2 -q -1 u -2 ) A ǫ (u, m) + u ǫ ∆ d (q -1 u -1 ) (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) D ǫ (u, m) - ǫ q ǫ κ ǫ+1 2 κ * 1-ǫ 2 u ǫ u 2 -u -2 B ǫ (u, m) ,
where the notation

∆ d (u) = u ǫ κu + κ * u -1 (3.76)
is introduced. From the multiple actions (3.46), (3.47), Lemma 3.3 and Lemma 3.5, it follows that the action of W d (u, m 0 + 4s) on the Bethe vector |Ψ 2s d,ǫ (ū, m 0 ) is given by:

W d (u, m 0 + 4s)|Ψ 2s d,ǫ (ū, m 0 ) = Λ 2s d,ǫ (u, ū)|Ψ 2s d,ǫ (ū, m 0 ) - 2s i=1 (q + q -1 )u ǫ E d (u i , ūi ) (u 2 -u -2 )b(uu -1 i )b(quu i ) |Ψ 2s d,ǫ ({u, ūi }, m 0 )
where

Λ 2s d,ǫ (u, ū) = u ǫ ∆ d (u) (u 2 -u -2 ) (qu 2 -q -1 u -2 ) 2s j=1 f (u, u j )Λ ǫ 1 (u) + u ǫ ∆ d (q -1 u -1 ) (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) 2s j=1 h(u, u j )Λ ǫ 2 (u) +(-1) 2s+1 ǫq ǫ κ (1+ǫ)/2 κ * (1-ǫ)/2 δ d 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) 2s i=1 b(uu -1 i )b(quu i ) (3.77) and E d (u i , ūi ) = b(u 2 i ) b(qu 2 i ) ∆ d (u i ) 2s j=1,j =i f (u i , u j )Λ ǫ 1 (u i ) + ∆ d (q -1 u -1 i ) 2s j=1,j =i h(u i , u j )Λ ǫ 2 (u i ) -(-1) 2s ǫ(q -q -1 ) -1 q ǫ κ (1+ǫ)/2 κ * (1-ǫ)/2 δ d u -ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) 2s j=1,j =i b(u i u -1 j )b(qu i u j ) . (3.78) Requiring E d (u i , ūi ) = 0 for i = 1, . . . ,
Λ 2s d,ǫ (u, ū) + (q + q -1 ) 2 ρ (u 2 -u -2 ) (q 2 u 2 -q -2 u -2 ) ηκ * + η * κ + (ηκ + η * κ * ) qu 2 + q -1 u -2 q + q -1 (3.79)
is a meromorphic function in the variable u. To be equal to a constant, we study its singular part. Using (3.37) and the expressions in Appendix B, one finds that the singular points are located at:

u ∈ {±u j , ±q -1 u -1 j , j = 1, ..., M } . (3.80)
The sum of all the residues at these points vanishes. Extracting the constant part of (3.79), one obtains (3.73), (3.74).

3.4. Diagonalization of the generic case κ = 0, κ * = 0, κ ± = 0. We now consider the most general case, associated with the non-diagonal matrix K + (u) given by (2.20). Similarly to the special and diagonal cases, the first step is to express the Heun-Askey-Wilson operator (1.16) in terms of the dynamical operators (3.16)-(3.19). According to the gauge transformation parametrized by α, β applied in (3.12), recall the expressions (3.53), (3.54), (3.67), (3.68) for the elements A and A * . In addition, from (2.4)-(2.7) and (3.16)-(3.19) one gets the following expressions for [A, A * ] q and [A * , A] q in terms of the dynamical operators:

[A * , A] q = - αβρχ -1 q -ǫ(m+1) u ǫ α -q -2ǫ(m+1) β 1 qu 2 -q -1 u -2 A ǫ (u, m) - 1 u 2 -u -2 D ǫ (u, m) (3.81) + ρχ -1 u ǫ (α -q -2ǫm β)(u 2 -u -2 ) α 2 q ǫ(m+2) B ǫ (u, m) -β 2 q ǫ(-3m+2) C ǫ (u, m) -ρ qu 2 + q -1 u -2 q 2 -q -2 + ω q -q -1 , [A, A * ] q = χq -ǫ(m+1) u ǫ α -q -2ǫ(m+1) β 1 qu 2 -q -1 u -2 A ǫ (u, m) - 1 u 2 -u -2 D ǫ (u, m) - χe -mǫ u ǫ (α -q -2ǫm β)(u 2 -u -2 ) (B ǫ (u, m) -C ǫ (u, m)) -ρ qu 2 + q -1 u -2 q 2 -q -2 + ω q -q -1 .
For generic parameters α, β and integer m, in terms of the dynamical operators recall that the transfer matrix is given by (3.20). In order to simplify the analysis of the spectral problem, as a second step let us fix the gauge parameters α, β and m such that the coefficients of B ǫ (u, m) and C ǫ (u, m) in (3.20) vanish. For convenience and without loss of generality, let us choose the following parametrization:

κ + = q ϕ 2 , κ -= - q ϕ′ 2 , κ = ρ 1/2 q (ϕ+ϕ ′ )/2 cosh(ξ) , κ * = ρ 1/2 q (ϕ+ϕ ′ )/2 cosh(ξ ′ ) . (3.82)
Then, for the choice of gauge parameters α = -ǫρ -1/2 q (ϕ-ϕ ′ )/2-ǫ(1+m0+2M) e -ξ(1+ǫ)/2+ξ ′ (1-ǫ)/2 , (3.83) β = -ǫρ -1/2 q (ϕ-ϕ ′ )/2-ǫ(1-m0-2M) e ξ(1+ǫ)/2-ξ ′ (1-ǫ)/2 , from the expression of the transfer matrix (3.20) at m = m 0 + 2M and (2.21), for any choice of ǫ = ±1 one finds the following expression for the Heun-Askey-Wilson element (1.16):

I(κ, κ * , κ + , κ -) = 1 (u 2 -u -2 ) u ǫ ∆ g (u) (qu 2 -q -1 u -2 ) A ǫ (u, m 0 + 2M ) + u ǫ ∆ g (q -1 u -1 ) (q 2 u 2 -q -2 u -2 ) D ǫ (u, m 0 + 2M ) (3.84) + q (ϕ+ϕ ′ )/2 (q + q -1 ) 2 ρ 1/2 (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) η cosh(ξ ′ ) + η * cosh(ξ) + (η cosh(ξ) + η * cosh(ξ ′ )) qu 2 + q -1 u -2 q + q -1 - (q ϕ χ -1 -q ϕ ′ χ) 2 ρ qu 2 + q -1 u -2 q 2 -q -2 + ω q -q -1
where we have denoted

∆ g (u) = 1 2 e -ξ q (ϕ+ϕ ′ )/2 ρ 1/2 (e -ξ ′ u + e ξ u -1 )(e ξ(1-ǫ)/2+ξ ′ (1+ǫ)/2 u + e ξ(1+ǫ)/2+ξ ′ (1-ǫ)/2 u -1 ) . (3.85)
Below, we will consider the action of the dynamical operators on the reference state |Ω + or |Ω -. For the choice of the gauge parameters α, β (3.83), the action of the dynamical operators A ǫ (u, m 0 ) and D ǫ (u, m 0 ) is now considered 11 . Lemma 3.6. For the choice of gauge parameters (3.83), the action of the operators A ǫ (u, m 0 ) and D ǫ (u, m 0 ) on |Ω ǫ is given by

π(A ǫ (u, m 0 ))|Ω ǫ = Λ ǫ 1 (u)|Ω ǫ + c M ǫ π(B ǫ (u, m 0 -2))|Ω ǫ , (3.86) π(D ǫ (u, m 0 ))|Ω ǫ = Λ ǫ 2 (u)|Ω ǫ - (q 2 u 2 -q -2 u -2 ) (qu 2 -q -1 u -2 ) c M ǫ π(B ǫ (u, m 0 -2))|Ω ǫ , (3.87)
11 Similar relations were derived in [START_REF] Avan | Modified algebraic Bethe ansatz for XXZ chain on the segment -III -Proof[END_REF] where Λ ǫ 1 (u) and Λ ǫ 2 (u) are given by (3.37) and

c M ǫ = e -ξ(1+ǫ)+ξ ′ (1-ǫ) q -2ǫ(2M+1) q ϕ/2 -e (µ ′ +ξ)(1+ǫ)/2+(µ-ξ ′ )(1-ǫ)/2 q (ϕ ′ +ν-ν ′ )/2+ǫ(2M-2s+1) ρ 1/2 q ϕ/2 -e (µ ′ -ξ)(1+ǫ)/2+(µ+ξ ′ )(1-ǫ)/2 q (ϕ ′ +ν-ν ′ )/2-ǫ(2M+2s+1) ρ 1/2 . (3.88)
Proof. We show (3.86) for ǫ = +1. Recall |Ω + ≡ |θ * 0 . Using the explicit expressions (A.1) and (A.2) for m → m -2 together with (3.2) it follows:

π(A ǫ (u, m 0 )) -c M + π(B ǫ (u, m 0 -2)) |Ω + = u -1 (u 2 -u -2 ) (αq 2m0 -βq 2 ) const ′ 1 (c M + )a * 01 |θ * 1 + eigen(u, c M + )|θ * 0 where const ′ 1 (c M + ) = (1 -c M + ) χq m0 ρ (qθ * 1 -q -1 θ * 0 ) - (αβq m0+2 -c M + β 2 q -m0+4 ) χ (qθ * 0 -q -1 θ * 1 )
+ βq + αq 2m0+1c M + β(q 3 + q) . Requiring const ′ 1 (c M + ) ≡ 0 determines the coefficient c M + in terms of the Leonard pair's data and the gauge parameters α, β:

c M + = (χ 2 + αβρ)q m0+2 θ * 1 -(χ 2 + αβρq 4 )q m0 θ * 0 + (αq 2m0+2 + βq 2 )ρχ (χ 2 q m0+2 + β 2 ρq -m0+4 )θ * 1 -(χ 2 q m0 + β 2 ρq -m0+6 )θ * 0 + (q 4 + q 2 )βρχ . (3.89)
Inserting (3.3) for M = 0, 1 in the above expression, for the choice of gauge parameters (3.83) and using the parametrization (3.32) one gets (3.88). Then, inserting (3.88) in eigen(u, c M + ), one finds eigen 0 (u, c M + ) = Λ + 1 (u). This completes the proof of (3.86) for ǫ = +1. The other relations (3.86) for ǫ = -1 and (3.87) are shown similarly, so we skip the details.

As for the case of 'diagonal' parameters studied in the previous section, a crucial ingredient for the solution of the spectral problem of (3.84) is the conjecture below (see also similar relations in [BP15, ABGP15]) which is a generalization of Lemma 3.5. For a generic set of parameters ū = {u 1 , u 2 , ..., u 2s }, define the Bethe vector:

δg π(B ǫ (u, m 0 + 4s))|Ψ 2s g,ǫ (ū, m 0 ) = (3.92) δ g u -ǫ b(u 2 ) 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) 2s i=1 b(uu -1 i )b(q -1 u -1 u -1 i ) |Ψ 2s g,ǫ (ū, m 0 ) -δ g 2s i=1 u -ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) b(uu -1 i )b(q -1 u -1 u -1 i ) 2s j=1,j =i b(u i u -1 j )b(q -1 u -1 i u -1 j ) |Ψ 2s g,ǫ ({u, ūi }, m 0 )
where

δ g = (-1) 2s q ν ′ 4 (q ϕ/2 -e (µ-ξ ′ )(1-ǫ)/2-(µ ′ +ξ)(1+ǫ)/2 q (ϕ ′ +ν-ν ′ )/2-2s-1 ρ 1/2 ) (3.93) ×(q ϕ/2 -e -(µ-ξ ′ )(1-ǫ)/2+(µ ′ +ξ)(1+ǫ)/2 q (ϕ ′ +ν-ν ′ )/2+2s+1 ρ 1/2 ) , δg = e -3ξ(1+ǫ)/2+ξ ′ (1-ǫ)/2 q (ϕ+ϕ ′ )/2-(1+2ǫ)(4s+1) (1 -e ξ ′ (1-ǫ)+ξ(1+ǫ) q 2(4s+1) ) × (q ϕ/2 -e (µ-ξ ′ )(1-ǫ)/2+(µ ′ +ξ)(1+ǫ)/2 q (ϕ ′ +ν-ν ′ )/2+ǫ(2s+1) ρ 1/2 )ρ 1/2 2(q ϕ/2 -e (µ+ξ ′ )(1-ǫ)/2+(µ ′ -ξ)(1+ǫ)/2 q (ϕ ′ +ν-ν ′ )/2-ǫ(6s+1) ρ 1/2 ) .
This conjecture has been checked with Mathematica for small values of s = 1/2, 1, 3/2. Note that the case s = 1/2 has been proved in [START_REF] Avan | Modified algebraic Bethe ansatz for XXZ chain on the segment -III -Proof[END_REF] using the separation of variables (SoV) basis, and the method can be generalized for arbitrary s. For generic s, it might be interesting to give a proof by analogy with Appendix C using the theory of Leonard pairs. This might be studied elsewhere.

We now turn to the solution of the spectral problem. Recall the parametrization (3.82).

Proposition 3.4. For ǫ = ±1, one has:

π (I(κ, κ * , κ + , κ -)) |Ψ 2s g,ǫ (ū, m 0 ) = Λ 2s g,ǫ |Ψ 2s g,ǫ (ū, m 0 ) (3.94) with Λ 2s g,+ = κ * θ * 2s + κ * θ * 2s | µ ′ →ξ cosh(µ) cosh(ξ ′ ) + - κ * 2 cosh(ξ ′ ) θ * 3s+1/2 | µ ′ →µ ′ +ξ + (-1) 2s+1 δ g [2s + 1] q (v 2 + v -2 )
-ω (χ -1 κ + + χκ -) (qq -1 ) + κ * q ν+ν ′ 4 cosh(ξ ′ ) (qq -1 )(q 2s e µ ′ +ξq -2s e -µ ′ -ξ ) -(-1) 2s+1 δ g 2s j=1 (qu 2 j + q -1 u -2 j ) ,

Λ 2s g,-= κθ 2s + κθ 2s | µ→-ξ ′ cosh(µ ′ ) cosh(ξ) + - κ 2 cosh(ξ) θ 3s+1/2 | µ→µ-ξ ′ + (-1) 2s+1 δ g [2s + 1] q (v 2 + v -2 ) -ω (χ -1 κ + + χκ -) (q -q -1 ) + κ 2 cosh(ξ) (q -q -1 )θ 2s | µ→µ-ξ -(-1) 2s+1 δ g 2s j=1 (qu 2 j + q -1 u -2 j ) ,
where the set ū satisfies the (inhomogeneous) Bethe equations:

b(u 2 i ) b(qu 2 i ) 1 2 e -ξ q (ϕ+ϕ ′ )/2 ρ 1/2 (e -ξ ′ u i + e ξ u -1 i )(e ξ(1-ǫ)/2+ξ ′ (1+ǫ)/2 u i + e ξ(1+ǫ)/2+ξ ′ (1-ǫ)/2 u -1 i ) 2s j=1,j =i f (u i , u j )Λ ǫ 1 (u i ) - 1 2 e -ξ q (ϕ+ϕ ′ )/2 ρ 1/2 (qe ξ u i + q -1 e -ξ ′ u -1 i )(qe ξ(1+ǫ)/2+ξ ′ (1-ǫ)/2 u i + q -1 e ξ(1-ǫ)/2+ξ ′ (1+ǫ)/2 u -1 i ) 2s j=1,j =i h(u i , u j )Λ ǫ 2 (u i ) -(-1) 2s δ g (q -q -1 ) -1 u -ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) 2s j=1,j =i b(u i u -1 j )b(qu i u j ) = 0 for i = 1, . . . , 2s.
Proof. For convenience, define the element:

W g (u, m) = I(κ, κ * , κ + , κ -) + (q ϕ χ -1 -q ϕ ′ χ) 2 ρ qu 2 + q -1 u -2 q 2 -q -2 + ω q -q -1 - ρ -1/2 q (ϕ+ϕ ′ )/2 (q + q -1 ) 2 (u 2 -u -2 )(q 2 u 2 -q -2 u -2 )
η cosh(ξ ′ ) + η * cosh(ξ) + (η cosh(ξ) + η * cosh(ξ ′ )) qu 2 + q -1 u -2 q + q -1 .

Explicitly, using (3.84) evaluated at the point m = m 0 + 2M we have:

W g (u, m 0 + 2M ) = 1 (u 2 -u -2 ) u ǫ ∆ g (u) (qu 2 -q -1 u -2 ) A ǫ (u, m 0 + 2M ) + u ǫ ∆ g (q -1 u -1 ) (q 2 u 2 -q -2 u -2 ) D ǫ (u, m 0 + 2M ) . (3.95)
From the multiple actions (3.46), (3.47), Lemma 3.6 and Conjecture 1, it follows that the action of W g (u, m 0 +4s) on the Bethe vector |Ψ 2s g,ǫ (ū, m 0 ) is given by:

W g (u, m 0 + 4s)|Ψ 2s g,ǫ (ū, m 0 ) = Λ 2s g,ǫ (u, ū)|Ψ 2s g,ǫ (ū, m 0 ) - 2s i=1 (q + q -1 )u ǫ E g (u i , ūi ) (u 2 -u -2 )b(uu -1 i )b(quu i ) |Ψ 2s g,ǫ ({u, ūi }, m 0 ) (3.96) where Λ 2s g,ǫ (u, ū) = u ǫ ∆ g (u) (u 2 -u -2 )(qu 2 -q -1 u -2 ) 2s j=1 f (u, u j )Λ ǫ 1 (u) + u ǫ ∆ g (q -1 u -1 ) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) 2s j=1 h(u, u j )Λ ǫ 2 (u) +(-1) 2s δ g 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) 2s i=1 b(uu -1 i )b(quu i ) , (3.97) and E g (u i , ūi ) = - b(u 2 i ) b(qu 2 i ) ∆ g (u i ) 2s j=1,j =i f (u i , u j )Λ ǫ 1 (u i ) + ∆ g (q -1 u -1 i ) 2s j=1,j =i h(u i , u j )Λ ǫ 2 (u i ) (3.98) +(-1) 2s δ g (q -q -1 ) -1 u -ǫ i b(u 2 i ) 2s k=0 b(q 1/2+k-s vu i )b(q 1/2+k-s v -1 u i ) 2s j=1,j =i b(u i u -1 j )b(qu i u j )
.

Requiring E g (u i , ūi ) = 0 for i = 1, . . . , 2s, one gets the Bethe equations below (3.94). To determine the eigenvalues Λ 2s g,ǫ in (3.94), observe that Λ 2s g,ǫ (u, ū) +

q (ϕ+ϕ ′ )/2 (q + q -1 ) 2 ρ 1/2 (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) η cosh(ξ ′ ) + η * cosh(ξ) + (η cosh(ξ) + η * cosh(ξ ′ )) qu 2 + q -1 u -2 q + q -1 - (q ϕ χ -1 -q ϕ ′ χ) 2 ρ qu 2 + q -1 u -2 q 2 -q -2 + ω q -q -1 (3.99)
is a meromorphic function in the variable u. To be equal to a constant, we study its singular part. The singular points are located at:

u ∈ {±u j , ±q -1 u -1 j , j = 1, ..., M } . (3.100)
The sum of all the residues at these points vanishes. Extracting the constant part of (3.99), one obtains the eigenvalues Λ 2s g,± below (3.94). 3.5. Bethe ansatz equations and Bethe states revisited. In this subsection, it is shown that each system of Bethe ansatz equations previously derived can be rewritten in terms of the polynomials (3.104) in the 'symmetrized' Bethe roots (3.103), see Proposition 3.5. Thus, solving the Bethe ansatz equations is reduced to finding the solutions of a system of polynomial equations instead of Laurent polynomial equations, see Proposition 3.6. Also, an expansion formula for the Bethe states in the Poincaré-Birkhoff-Witt basis of the Askey-Wilson algebra is given, exhibiting the explicit dependence in the variables {U i |i = 1, ..., M }, see Corollary 3.1. This leads to the construction of two different eigenbases for the Leonard pairs, see Proposition 3.7.

3.5.1. An alternative presentation for the Bethe ansatz equations. For the three cases studied in the previous subsections, observe that each system of Bethe ansatz equations in Propositions 3.1, 3.2, 3.3 and 3.4 enjoys the symmetries

u i ←→ ±q -1 u -1 i , u i ←→ -u i (3.101) u j ←→ ±q -1 u -1 j , u j ←→ -u j for j = i . (3.102)
This suggests that each system of Bethe equations admits an alternative presentation as a system of equations written solely in terms of the 'symmetrized' variables

U i = qu 2 i + q -1 u -2 i q + q -1 with i = 1, ..., M, (3.103)
as we now show. Recall the notation Ūi = {U 1 , U 2 , ..., U i-1 , U i+1 , ..., U M } with (3.103). Define the polynomial:

P M a (U i , Ūi ) = M-1 k=0 (-1) k (q + q -1 ) M-1 2 M-1-k e k ( Ūi ) M-1-k l=0 M -1 -k l (q -q -1 ) l U M-1-k-l i (q 2 + q -2 ) 1+k+l-M g 2[ l 2 ] 0 (U i )g (p[l]) a,ǫ (U i ) + ∆a H(U i ) (3.
=        0 for a = sp, -2(-1) 2s ǫ q -(ν+ν ′ )/2+ǫ (q-q -1 ) κ * 1-ǫ 2 κ 1+ǫ 2 δ d for a = d, 2 (-1) 2s q -(ν+ν ′ )/2
(q-q -1 ) δ g for a = g .

Note that P M a (U i , Ūi ) is of maximal degree M in the variable U i for the case a = sp, and of maximal degree 2s + 1 for a ∈ {d, g}.

Example 5. For the special case a = sp and ǫ = +, for M = 1, 2, 3, the polynomials are given by:

P 1 sp (U 1 ) = (q + q -1 )s 1-2s U 1 -r 0 -s 1 (v 2 + v -2 ) , P 2 sp (U 1 , U 2 ) = (q + q -1 ) 2 s 3-2s U 2 1 -s 1-2s U 1 U 2 -(q + q -1 ) (r 2 + s 3 (v 2 + v -2 ))U 1 -(r 0 + s 1 (v 2 + v -2 ))U 2 +(q 2 -q -2 )(q 2s -q -2s )s 1 + (q + q -1 )r 2s+1 (v 2 + v -2 )) , P 3 sp (U 1 , {U 2 , U 3 }) = (q + q -1 ) 3 (s 5-2s U 3 1 -s 3-2s U 2 1 (U 2 + U 3 ) + s 1-2s U 1 U 2 U 3
) +U 1 q -2s-3 (qq -1 )(q + q -1 ) 2 (e µ ′ (q 4s+2 + q 4s+6 -2q 6 ) + e -µ ′ (-2q 4s + q 4 + 1))

+(q + q -1 )r 2s+4 (v 2 + v -2 ) +(U 2 + U 3 )(q + q -1 )(s 1 (q -2s -q 2s )(q 2 -q -2 ) -r 2s+2 (v 2 + v -2 )) -U 2 U 3 (q + q -1 ) 2 (r 0 + s 1 (v 2 + v -2 )) +(U 1 U 2 + U 1 U 3 )(q + q -1 ) 2 (r 2 + s 3 (v 2 + v -2 )) -U 2 1 (q + q -1 ) 2 (r 4 + s 5 (v 2 + v -2 )) -(q 2 -q -2 )((q 2s -q -2s )r 2s+2 -s 1 (q 2 -q -2 )(v 2 + v -2 ))
where r j = 2 cosh(µ)(q 2s-jq -2s+j ) , s j = e µ ′ q je -µ ′ q -j .

Example 6. For the diagonal case a = d and ǫ = +, for M = 1 (s = 1/2) and M = 2 (s = 1), the polynomials are given by:

P 1 d (U 1 ) = - (q + q -1 ) 2 (q -q -1 ) e -µ ′ κU 2 1 +2(q + q -1 )(sinh(µ ′ )κ * + (-cosh(µ) + e -µ ′ q -1 (q -q -1 ) -1 (v 2 + v -2 ))κ)U 1 +(-2 cosh(µ)(q -q -1 ) -(qe µ ′ -q -1 e -µ ′ )(v 2 + v -2 ))κ * -(e µ ′ (q -q -1 ) 2 -2q -1 (q -q -1 ) cosh(µ)(v 2 + v -2 ) + q -2 e -µ ′ (v 2 + v -2 ) 2 )(q -q -1 ) -1 κ , P 2 d (U 1 , U 2 ) = - (q + q -1 ) 3 (q -q -1 ) q -1 e -µ ′ κ(U 3 1 + (1 -q 2 )U 2 1 U 2 )
+U 2 ((-2q -1 (q + q -1 ) cosh(µ)(v 2 + v -2 ) + q -1 (qq -1 )(q + q -1 ) 2 (e µ ′ qe -µ ′ q -1 ))κ +(2(qq -1 )(q + q -1 ) 2 cosh(µ) + (q + q -1 )(e µ ′ qe -µ ′ q -1 )(v 2 + v -2 ))κ * ) +U 1 U 2 (q + q -1 ) 2 ((2q cosh(µ)q -2 e -µ ′ (v 2 + v -2 ))κ -(e µ ′ q -1e -µ ′ q)κ * )

+U 2 1 (q + q -1 ) 2 (q -2 ( e -µ ′ (q + 2q -1 ) (qq -1 ) (v 2 + v -2 ) -2q cosh(µ))κ + (e µ ′ qe -µ ′ q -1 )κ * ) +U 1 ((2 cosh(µ)q -3 (q + q -1 )(v 2 + v -2 ) -e -µ ′ 1 + 2q -2 + 2q -4 + q -6 (qq -1 ) (v 4 + v -4 ) -e -µ ′ (q + q -1 )(q 3 + 2q -3 + 3q -5 ) (qq -1 ) )κ -(q + q -1 )(e

µ ′ q 3 -e -µ ′ q -3 )(v 2 + v -2 )κ * ) +(q 2 -q -2 )((q 2 -q -2 )(e µ ′ q -e -µ ′ q -1 ) -2 cosh(µ)(v 2 + v -2 ))κ * +(-2q -1 (q 2 -q -2 ) 2 cosh(µ) + e -µ ′ q -3 (q -q -1 ) -1 (v 6 + v -6 ) -e µ ′ (q 2 -q -2 )(v 2 + v -2 ) +e -µ ′ (2q -q -1 + q -5 + q -7 )(q -q -1 ) -1 (v 2 + v -2 ))κ .
The proof of the following proposition is given in Appendix D. Recall (3.61), (3.78) or (3.98).

Proposition 3.5.

E a (u i , ūi ) = u -ǫ i b(u 2 i )q (ν+ν ′ )/2 2 M j =i b(u i /u j )b(qu i u j ) P M a (U i , Ūi ) (3.105)
with M = 0, 1, ..., 2s for a = sp, M = 2s for a ∈ {d, g} .

Recall that the Bethe ansatz equations follow from requiring E a (u i , ūi ) = 0 (3.106) for all i = 1, ..., M . Let {u 1 , u 2 , ..., u M } denote the Bethe roots that satisfy the Bethe ansatz equations in Propositions 3.1, 3.2, 3.3 and 3.4. According to (3.105) and (3.106), it follows: Proposition 3.6. The Bethe roots {u 1 , u 2 , ..., u M } are determined from (3.103) where {U 1 , U 2 , ..., U M } are solutions of the system of polynomial equations:

P M a (U i , Ūi ) = 0 for i = 1, ..., M. (3.107)
In general, there may be solutions such that U i = U j for i = j. However, these solutions are not compatible with (3.105). In the following, a solution is called admissible if the condition U i = U j for any i = j is fullfiled. We now study the subset of admissible solutions Ū = {U 1 , U 2 , ..., U M } of (3.107) for the special a = sp, diagonal a = d and generic a = g cases.

Admissible solutions for a = sp: For M = 1, it is clear that P 1 sp (U 1 ) has a unique solution. For M = 2, we have to solve the set of equations

P 2 sp (U 1 , U 2 ) = 0 , P 2 sp (U 2 , U 1 ) = 0. (3.108)
Numerically, we find that this system admits 4 = 2 2 solutions. Two solutions are such that U 1 = U 2 , which are discarded (not admissible). The other two solutions are distinct and related by permutation. So, up to permutation, the admissible solution of (3.108) associated with (3.105) is unique. For M = 3, we have to solve the set of equations

P 3 sp (U 1 , {U 2 , U 3 }) = 0 , P 3 sp (U 2 , {U 1 , U 3 }) = 0 , P 3 sp (U 3 , {U 1 , U 2 }) = 0. (3.109)
Numerically, we find 27 = 3 3 solutions. Among them, only 6 = 3! are admissible. Furthermore, they are all related by a permutation. So, up to permutation, again we conclude that the admissible solution of (3.109) is unique. For M = 4, numerically we find 256 = 4 4 solutions to (3.107). Among these, only 24 = 4! are admissible and related by permutation. Again, this manifests the fact that, up to permutation, the admissible solution of the polynomial equations (3.107) is unique.

More generally, by Bezout's theorem [START_REF] Griffiths | Principles of Algebraic Geometry[END_REF]p. 670] the total number of solutions of the system (3.107) for a = sp is M M . Previous numerical analysis suggests that the number of admissible solutions is M !, all related by permutation. Although we have no proof at the moment for M generic, we formulate the following conjecture:

Conjecture 2. The system of polynomial equations (3.107) for a = sp admits a unique admissible solution Ū = {U 1 , U 2 , ..., U M } up to permutation.

Admissible solutions for a = d and a = g: For M = 1 (s = 1/2), P 1 a (U 1 ) is a polynomial of degree 2 with two distinct solutions. For M = 2 (s = 1), we have to solve the set of equations

P 2 a (U 1 , U 2 ) = 0 , P 2 a (U 2 , U 1 ) = 0 (3.110)
where the two polynomials are of total degree 3. Numerically, we find 9 = 3 2 solutions. Among those, only 6 solutions are admissible. Up to permutation, we have 3 distinct solutions. For M = 3 (s = 3/2), we have to solve the set of equations

P 3 a (U 1 , {U 2 , U 3 }) = 0 , P 3 a (U 2 , {U 1 , U 3 }) = 0 , P 3 a (U 3 , {U 1 , U 2 }) = 0 (3.111)
where the three polynomials are of total degree 3. Numerically, we find 64 = 4 3 solutions. Among those, only 4 are admissible and distinct.

For generic M = 2s, by Bezout's theorem the total number of solutions of the system (3.107) for a = d or a = g is (2s + 1) 2s . Previous numerical analysis suggests that the number of admissible solutions is 2s + 1 which matches with the dimension of the vector space V, as expected. We formulate the following conjecture: Conjecture 3. The system of polynomial equations (3.107) for a = d, g with M = 2s admits 2s + 1 distinct admissible solutions.

For all three cases a = {sp, d, g}, we wish to observe that the numerical analysis of (3.107) is simpler compared with the usual analysis of the Bethe ansatz equations in terms of the original Bethe roots {u i }.

3.5.2. Bethe states and the PBW basis of Askey-Wilson algebra. In the framework of the algebraic Bethe ansatz applied to the K-matrix given in Proposition 2.1, by construction any Bethe state is a polynomial in the elements A, A * of the Askey-Wilson algebra acting on a certain reference state, |Ω + or |Ω -. In general, this polynomial is not written in terms of linearly independent monomials in A, A * . Furthermore, according to the previous subsection, it is expected that any Bethe state can be written, up to an overall factor, in terms of the 'symmetrized' variables (3.103). Below, an expansion formula for any (off-shell or on-shell) Bethe state is given in the linear basis of the Askey-Wilson algebra, and the dependency in the variables U i is exhibited.

As a preliminary, recall that a linear basis for the Askey-Wilson algebra is known, see e.g. [T11, Theorem 4.1]. For convenience, we introduce the element B such that:

B = A, A * q (⇒ A * A = q 2 AA * -qB) . (3.112)
In terms of the elements A, A * , B, from the Askey-Wilson relations (1.14), (1.15) one obtains:

A, B q -1 = ρ A * + ω A + ηI ⇒ BA = q -2 AB -ρq -1 A * -ωq -1 A -ηq -1 , (3.113) A * , B q = -ρ A -ω A * -η * I ⇒ BA * = q 2 A * B + ρqA + ωqA * + η * q .
(3.114) Thus, using the ordering relations (3.112),(3.113),(3.114), any polynomial in the elements A, A * , B can be written in terms of the linearly independent monomials

A i A * j B k for any i, j, k ≥ 0 (3.115)
that form a Poincaré-Birkhoff-Witt basis (see definition in [D93, p. 299]) of the Askey-Wilson algebra. We refer the reader to [T11] for more details.

Consider the dynamical operators B ± (u, m). Recall their expressions in terms of the elements A, A * of the Askey-Wilson algebra given in Appendix A. Define

U = qu 2 + q -1 u -2 q + q -1 . (3.116)
According to the ordering prescription (3.115), in the PBW basis of the Askey-Wilson algebra (3.115) we obtain (the notation (1.21) is used):

B + (u, m) = βb(u 2 )(q+q -1 ) u(αq 2m+2 -β) χq m+1 b(q) βρ AA * -(χ 2 q m+2 +β 2 ρq -m ) (q+q -1 )βρχ B -U A * + A + h + 0 (U, m)I , (3.117) B -(u, m) = βub(u 2 )(q+q -1 ) (αq -2m-2 -β) χq -m+1 b(q) βρ AA * -(χ 2 q -m+2 +β 2 ρq m ) (q+q -1 )βρχ B + U A -A * + h - 0 (U, m)I (3.118)
where

h ǫ 0 (U, m) = (χ 2 q ǫm -β 2 ρq -ǫm ) βχ(q 2 -q -2 ) U + ω ρ + η * ρ for ǫ = +1 , -η ρ for ǫ = -1
.

For an arbitrary product of the dynamical operators B ± (u, m), it is easy to extract the general structure in terms of the ordered monomials (3.115). For instance, consider a product of two operators B ǫ (u 1 , m ′ ), B ǫ (u 2 , m). In its ordered form, it reads as a polynomial of total degree four in A, A * and B with the prescription that the power of each element is at most two. More generally, by induction it follows:

Lemma 3.7. For any integer M ≥ 1 and any set ū = {u 1 , u 2 , ..., u M }:

B ǫ (u 1 , m + 2(M -1)) • • • B ǫ (u M , m) = β α M (q + q -1 ) M q 2M(M-m-2) ( β α q -2m-2 ; q 4 ) M M l=1 (u 2 l -u -2 l ) u ǫ l i, j, k ≤ M i + j + k ≤ 2M ζ [M] i,j,k ( Ū )A i A * j B k , where ζ [M] i,j,k ( Ū ) are polynomials of total degree M in the variables U i = qu 2 i +q -1 u -2 i q+q -1
.

Remark 1. Note that both expressions (3.117), (3.118) are regular in the parameter β. Thus, for the special case β = 0 the product of B ǫ (u i , m i ) in Lemma 3.7 is well-defined.

Applying the ordering prescription (3.115) to any product of dynamical operators B ǫ (u, m), the polynomials ζ M i,j,k (ū) can be derived recursively.

Example 7. For ǫ = +, the non-vanishing coefficients ζ

[M] i,j,k (ū) are given by:

For M = 1: ζ [1] 110 = ζ 110 (m) = χq m+1 (q -q -1 ) βρ , ζ [1] 001 = ζ 001 (m) = - (χ 2 q m+2 + β 2 ρq -m ) (q + q -1 )βρχ , ζ [1] 100 = ζ 100 = 1, ζ [1] 010 (U 1 ) = ζ 010 (U 1 ) = -U 1 , ζ [1] 000 (U 1 ) = ζ 000 (U 1 , m) = h + 0 (U 1 , m).
For M = 2:

ζ [2] 220 = q 2 ζ 110 (m -2)ζ 110 (m), ζ [2] 111 = ζ 110 (m -2)ζ 001 (m) + ζ 110 (m)ζ 001 (m -2) -q 3 ζ 110 (m -2)ζ 110 (m) ζ [2] 101 = ζ 100 ζ 001 (m) + q -2 ζ 001 (m -2) -qζ 110 (m -2) , ζ [2] 011 = q 2 ζ 001 (m -2)ζ 010 (U 2 ) + ζ 001 (m)ζ 010 (U 1 ) -q 3 ζ 110 (m)ζ 010 (U 1 ), ζ [2] 120 = q 2 ζ 110 (m)ζ 010 (U 1 ) + ζ 110 (m -2)ζ 010 (U 2 ), ζ [2] 210 = ζ 100 (q 2 ζ 110 (m -2) + ζ 110 (m)), ζ [2] 110 = ζ 110 (m -2)ζ 000 (U 2 , m) + ζ 110 (m)ζ 000 (U 1 , m -2) +ζ 100 (q 2 ζ 010 (U 1 ) + ζ 010 (U 2 )) -q 2 ωζ 110 (m -2)ζ 110 (m) ζ [2] 200 = ρζ 110 (m)(q -1 ζ 001 (m -2) -q 2 ζ 110 (m -2)) + ζ 2 100 , ζ [2] 020 = ζ 010 (U 1 )ζ 010 (U 2 ) -ρq -1 ζ 001 (m -2)ζ 110 (m), ζ [2] 002 = ζ 001 (m -2)ζ 001 (m), ζ [2] 100 = ρqζ 001 (m -2)ζ 010 (U 2 ) -ρq 2 ζ 110 (m)ζ 010 (U 1 ) -η * q 2 ζ 110 (m -2)ζ 110 (m) -ωq -1 ζ 001 (m -2)ζ 100 + η * q -1 ζ 001 (m -2)ζ 110 (m) + (ζ 000 (U 1 , m -2) + ζ 000 (U 2 , m))ζ 100 , ζ [2] 010 = -ρq -1 ζ 001 (m -2)ζ 100 -ηζ 001 (m -2)ζ 110 (m) + ωqζ 001 (m -2)ζ 010 (U 2 ) -ωq 2 ζ 010 (U 1 )ζ 110 (m) + ζ 000 (U 1 , m -2)ζ 010 (U 2 ) + ζ 000 (U 2 , m)ζ 010 (U 1 ), ζ [2] 001 = -qζ 010 (U 1 )ζ 100 + ζ 000 (U 1 , m -2)ζ 001 (m) + ζ 000 (U 2 , m)ζ 001 (m -2), ζ [2] 000 = ζ 000 (U 1 , m -2)ζ 000 (U 2 , m) -ηq -1 ζ 001 (m -2)ζ 100 -η * q 2 ζ 010 (U 1 )ζ 110 (m) -qζ 010 (U 2 ) ζ 001 (m -2) .
Let V be the irreducible finite dimensional vector space on which the Leonard pair π(A), π(A * ) act. Denote the (off-shell or on-shell) Bethe state:

|Ψ ǫ (ū, m) = π(B ǫ (u 1 , m + 2(M -1)) • • • B ǫ (u M , m))|Ω ǫ . (3.119)
For M ≤ 2s, there are no relations besides the Askey-Wilson relations (1.14), (1.15), whereas for M ≥ dim( V) = 2s + 1, additional relations occur from the characteristic polynomials of π(A) and π(A * ), see (3.6). By Lemma 3.7, it follows: Corollary 3.1. For 0 ≤ M ≤ 2s, any Bethe state admits the expansion formula:

|Ψ M ǫ (ū, m) = β α M (q + q -1 ) M q 2M(M-m-2) ( β α q -2m-2 ; q 4 ) M M l=1 (u 2 l -u -2 l ) u ǫ l i, j, k ≤ M i + j + k ≤ 2M ζ M i,j,k ( Ū )π A i A * j B k |Ω ǫ . (3.120)
3.5.3. Eigenbases for Leonard pairs. The results of previous subsections suggest that the two families of onshell Bethe states (3.57) and (3.64) for M = 0, 1, 2, ..., 2s, provide two explicit examples of bases on which the Leonard pair π(A), π(A * ) act as (3.1), (3.2). Recall that solving the system of Bethe equations in Proposition 3.1 or 3.2 reduces to solving (3.107) for a = sp. Proposition 3.7. For each M = 0, 1, ..., 2s, assume the system of polynomial equations (3.107) for a = sp admits at least one admissible solution U = {U 1 , ..., U M } with (3.103). Given the Leonard pair π(A), π(A * ) defined in subsection 3.1.1, the set of Bethe states (3.57) (resp. (3.64)) with M = 0, 1, ..., 2s, forms a basis of the vector space V. We identify:

|θ M ≡ M l=1 u ǫ l (u 2 l -u -2 l ) |Ψ M sp,-(ū, m 0 ) for α according to Lemma 3.2 , |θ * M ≡ M l=1 u ǫ l (u 2 l -u -2 l ) |Ψ M sp,+ (ū, m 0 )
for α according to Lemma 3.1 .

Proof. By definition, recall that dim( V) = 2s + 1, see subsection 3.1.1. For q not a root of unity by (3.3) the (2s + 1)-eigenvalues θ M and θ * M are multiplicity-free, with M = 0, 1, ..., 2s. By Propositions 3.1 and 3.2, the Bethe states (3.57) and (3.64) are the eigenvectors of π(A) and π(A * ), respectively. Given M fixed and according to previous comments, each eigenspace has dimension one. Also, in view of the form (3.120) for different values of M the eigenvectors in the set {|Ψ M sp,ǫ (ū, m 0 ) } are linearly independent (no additional relations besides the Askey-Wilson relations).

Note that this is consistent with the numerical analysis of the polynomial equations (3.107) for a = sp: up to permutation, (3.107) admits a unique admissible solution {U 1 , ..., U M }.

Finally, let us mention that the numerical analysis for s = 1/2, 1, 3/2 done in subsection 3.5.1 suggests that it should be possible to generalize Proposition 3.7 to the cases a = d or a = g. Indeed, in these cases and M = 1, 2, 3, it was found that the total number of distinct admissible solutions is exactly 2s + 1 = dim( V). However, for generic values of M = 2s the proof that the spectrum in (3.72) or (3.94) is multiplicity-free -a key ingredient in Proposition 3.7 -is missing.

Baxter T-Q relations and the Heun-Askey-Wilson q-difference operator

In this section, homogeneous and inhomogeneous Baxter T-Q relations are deduced from the results of the previous section. Independently, for the special and diagonal cases we construct a q-difference operator realization of the Heun-Askey-Wilson element (1.16) acting on an infinite dimensional representation (π, V). For each choice of parameters, π(I(κ, κ * , 0, 0)) gives a specialization of the Heun-Askey-Wilson q-difference operator introduced in [BTVZ18, Proposition 5]. Its action on the Q-polynomial produces the T-Q relations. We also briefly comment on the similar result for the generic case. For completeness, the action of the dynamical operators on the unit is given. These results suggest an interpretation of the Q-polynomials as transition matrix coefficients.

(In)homogeneous Baxter T-Q relations.

4.1.1. Special case κ * = 0 and κ = κ ± = 0. For the special case κ * = κ ± = 0, κ = 0 and the reference state |Ω -, a functional relation for the spectrum of π(A) is given by (3.60) with (3.59), (3.53). For the special case κ = κ ± = 0, κ * = 0 and the reference state |Ω + , the functional relation for the spectrum of π(κ * A * ) -denoted below Λ * M sp,+ -that occurs in the proof of Proposition 3.2 is given by:

Λ * M sp,+ = κ * u (u 2 -u -2 )   Λ + 1 (u) (qu 2 -q -1 u -2 ) M j=1 f (u, u j ) + Λ + 2 (u) (q 2 u 2 -q -2 u -2 ) M j=1 h(u, u j )   (4.1) +κ * q u η(u -1 ) + q -1 u -1 η(u) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) .
Introduce the q-difference operator T ± such that T ± (f (u 2 )) = f (q ±2 u 2 ). By elementary computations, one finds:

M j=1 h(u, u j ) = T + Q M (U ) Q M (U ) and M j=1 f (u, u j ) = T -Q M (U ) Q M (U ) with Q M (U ) = M j=1 (U -U j ) (4.2)
where the notations (3.116), (3.103) are used. Furthermore

q u η(u -1 ) + q -1 u -1 η(u) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) = (q + q -1 ) 2 (η + η * U ) ρ(u 2 -u -2 )(q 2 u 2 -q -2 u -2 )
.

Replacing these expressions in (4.1), the proof of the following proposition is immediate.

Proposition 4.1. The eigenvalues Λ * M sp,+ of the Heun-Askey-Wilson operator π(I(0, κ * , 0, 0)) are given by the homogeneous Baxter T-Q relation

(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) Λ * M sp,+ Q M (U ) = κ * uΛ + 2 (u)T + Q M (U ) + κ * uΛ + 1 (u) (q 2 u 2 -q -2 u -2 ) (qu 2 -q -1 u -2 ) T -Q M (U ) +κ * (q + q -1 ) 2 ρ (η + η * U )Q M (U )
with (3.33), (3.35), (3.36).

According to the results of Section 3, the roots of the Q-polynomial are determined by (3.107) and the spectrum is given by

Λ * M sp,+ = κ * 2 q 1 2 (ν+ν ′ ) e -µ ′ q 2s-2M + e µ ′ q -2s+2M .
4.1.2. Special case κ = 0 and κ * = κ ± = 0. As mentionned previously, eq. (3.60) together with (3.59), (3.53) determine the spectrum of π(κA). Following the same arguments as for the special case κ * = 0 just described, an homogenous Baxter T-Q relation characterizing the spectrum of π(κA) can be written. However, it is interesting to observe that specializing the results of the diagonal case allows to exhibit an inhomogeneous T-Q relation characterizing the spectrum of π(κA) too. Consider (3.77) for ǫ = + and κ * = 0. It reads:

Λ 2s d,+ (u, ū)| κ * =0 = κu (u 2 -u -2 )   u 2 Λ + 1 (u) (qu 2 -q -1 u -2 ) 2s j=1 f (u, u j ) + q -2 u -2 Λ + 2 (u) (q 2 u 2 -q -2 u -2 ) 2s j=1 h(u, u j )   (4.3) +(-1) 2s+1 qκδ d 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) 2s i=1 b(uu -1 i )b(quu i )
.

Using (D.5),

2s i=1 b(u/u i )b(quu i ) = (q + q -1 ) 2s Q 2s (U )
and

q u η(u) + q -1 u -1 η(u -1 ) (u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) = (q + q -1 ) 2 (η * + ηU ) ρ(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) , (4.4) it follows:
Proposition 4.2. The eigenvalues Λ 2s sp,+ of the Heun-Askey-Wilson operator π(I(κ, 0, 0, 0)) are given by the inhomogeneous Baxter T-Q relation

(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) Λ M sp,+ Q 2s (U ) = κq -2 u -1 Λ + 2 (u)T + Q 2s (U ) + κu 3 Λ + 1 (u) (q 2 u 2 -q -2 u -2 ) (qu 2 -q -1 u -2 ) T -Q 2s (U ) +κ (q + q -1 ) 2 ρ (η * + ηU )Q 2s (U ) +κqδ d (-1) 2s+1 (q + q -1 ) 2s-2 (U 2 -1)H(U )
with (3.33), (3.35), (3.36), (3.71) and (D.5).

According to the results of Section 3, the roots of the Q-polynomial are now determined by the specialization of (3.107) for a = d, ǫ = + and κ * = 0. Also, the spectrum follows from Λ 2s d,+ | κ * =0 . Contrary to (3.65), it depends on the Bethe roots satisfying the inhomogeneous Bethe ansatz equation (3.107): {0.4 -0.646529i, 0.4 + 0.646529i} {0.3 + 0.676757i, 0.3 -0.676757i} {-0.4 + 0.859069i, -0.4 -0.859069i} 3/2 2.16667 (1) 3.04167 (0) 6.16667 (2) 24.0417 (3) {0.249828, 0.645086 + 1.62544i, 0.645086 -1.62544i} {0.201668, 0.529166 -1.67841i, 0.529166 + 1.67841i} {0.0420546, 0.108973 + 1.83288i, 0.108973 -1.83288i} {-2.42474 + 2.13318i, -2.42474 -2.13318i, -0.610518} 2 2.08333 (1) 3.33333 (2) 5.52083 (0) 12.0833 (3) 48.0208 (4) {0.299681 + 0.596429i, 0.299681 -0.596429i, 1.22532 -3.42503i, 1.22532 + 3.42503i} {0.214101 + 0.612971i, 0.214101 -0.612971i, 0.910899 -3.58263i, 0.910899 + 3.58263i} {0.0723844 -0.630621i, 0.0723844 + 0.630621i, 0.352616 + 3.81015i, 0.352616 -3.81015i} {-1.3877 -4.22879i, -1.3877 + 4.22879i, -0.287303 + 0.617224i, -0.287303 -0.617224i} {-11.8331 -3.15869i, -11.8331 + 3.15869i, -1.84146, -0.842389}

Λ 2s sp,+ = κ 2 q 1 2 (ν+ν ′ ) e -µ ′ q -2s (v 2 + v -2 )[2s] q + 2e µ ′ cosh(µ) -(1 + q 2 ) 2s i=1 U i . (4.
Table 1. Spectrum of π(κA) for κ = ν = ν ′ = v = 1, q = 2, µ = log(1/3), µ ′ = log(1/5).
Combined with the fact that the spectrum of π(κA) is always of the form (3.3) with (3.32), it implies that for each set of symmetrized Bethe roots satisfying (3.107), there exists an integer M ∈ {0, 1, ..., 2s} such that the following equality holds:

e -µ q -2s+2M + e µ q 2s-2M = e -µ ′ q -2s (v 2 + v -2 )[2s] q + 2e µ ′ cosh(µ) -(1 + q 2 ) 2s i=1 U i . (4.6)
This equality has been checked numerically for small values of s = 1/2, 1, 3/2, 2. Given s fixed, for each set of symmetrized Bethe roots the eigenvalues Λ 2s sp,+ are displayed in Table 1. The corresponding value of M is given in parenthesis.

(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) Λ 2s d,+ Q 2s (U ) = (4.7) = u∆ d (q -1 u -1 )Λ + 2 (u)T + Q 2s (U ) + u∆ d (u)Λ + 1 (u) (q 2 u 2 -q -2 u -2 ) (qu 2 -q -1 u -2 ) T -Q 2s (U ) + (q + q -1 ) 2 ρ (κη * + κ * η + (κη + κ * η * )U )Q 2s (U ) + κqδ d (-1) 2s+1 (U 2 -1) (q + q -1 ) 2s-2 H(U )
with (3.33), (3.35), (3.36), (3.71), (3.76) and (D.5).

According to the results of Section 3, the roots of the Q-polynomial are determined by (3.107) for a = d, ǫ = + and the spectrum of π(I(κ, κ * , 0, 0)) is given by (3.73). 4.1.4. Generic case κ, κ * , κ ± = 0. For the generic case it is also straightforward to derive the corresponding inhomogeneous Baxter T-Q relation from (3.97):

Proposition 4.4. The eigenvalues Λ 2s g,+ of the Heun-Askey-Wilson operator π(I(κ, κ * , κ + , κ -)) are given by the inhomogeneous Baxter T-Q relation

(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) Λ 2s g,+ Q 2s (U ) = (4.8) = u∆ g (q -1 u -1 )Λ + 2 (u)T + Q 2s (U ) + u∆ g (u)Λ + 1 (u) (q 2 u 2 -q -2 u -2 ) (qu 2 -q -1 u -2 ) T -Q 2s (U ) +ζ g (U )Q 2s (U ) + δ g (-1) 2s (U 2 -1) (q + q -1 ) 2s-2 H(U )
where

ζ g (U ) = (q + q -1 ) 2 κη * + κ * η + (κη + κ * η * )U ρ - (κ + χ -1 + κ -χ) (q -q -1 ) (U 2 -1) ρU + ω (4.9)
with (3.33), (3.35), (3.36), (3.93), (3.85), (3.82) and (D.5).

Here the roots of the Q-polynomial are determined by (3.107) for a = g, ǫ = + and the spectrum of π(I(κ, κ * , κ + , κ -)) is given below (3.94). 4.2. Heun-Askey-Wilson q-difference operator and the Q-polynomial. The starting point is a realization of the Askey-Wilson algebra (1.14), (1.15) in terms of q-difference operators. Define the elementary q-difference operators T ± such that T ± (f (z)) = f (q ±2 z). In Appendix E, examples are given. Applying the invertible transformation (2.28) to the realization (E.1), (E.2) of the Askey-Wilson algebra given in [T99, Section 5], we obtain the linear transformation denoted π: AW → [z, z -1 ] such that:

π(A) = q -1 z -1 φ(z)(T + -1) + q -1 zφ(z -1 )(T --1) (4.10) + 1 2 q ν+ν ′ 2 e -µ ′ q 2s 2e µ ′ cosh(µ) -(v 2 + v -2 )q -2s-1 + q -1 (z + z -1 ) , π(A * ) = φ(z)(T + -1) + φ(z -1 )(T --1) + 1 2 q (ν+ν ′ )/2 (e µ ′ q -2s + e -µ ′ q 2s ) (4.11) where φ(z) = 1 2 q ν+ν ′ 2 e -µ ′ q 2s (1 + qe -µ+µ ′ z)(1 + qe µ+µ ′ z)(1 -q -2s v 2 z)(1 -q -2s v -2 z) (1 -z 2 )(1 -q 2 z 2 ) . (4.12)
Note that π(A * ) is a specialization of the Askey-Wilson second-order q-difference operator, see e.g. [START_REF] Koekoek | The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue[END_REF] for details. By straightforward calculations, one finds that the corresponding structure constants in (1.14), (1.15) are given by (3.33)-(3.36). Using (4.10), (4.11), the realization of the q-commutators in A, A * is also computed. For instance, one finds: (4.13) π([A * , A] q ) = -q ν+ν ′ (qq -1 ) 4 (q + q -1 )(z + z -1 ) -(q 2s+1 + q -2s-1 )(v 2 + v -2 ) + 4 cosh(µ) cosh(µ ′ ) .

For simplicity, the expression of π([A, A * ] q ) is reported at the end of Appendix E. In general, the image of the element (1.16) by the map π reads as a fourth order q-difference operator that generalizes the Askey-Wilson operator (4.11). However, for the special and diagonal cases this expression reduces to a second-order q-difference operator. In these cases, it is easy to see that it is a specialization of the Heun-Askey-Wilson operator introduced in [BTVZ18, Proposition 5]. Below, we consider the action of π(I(κ, κ * , 0, 0)) on the Q-polynomial Q M (Z) with (4.2) where we denote:

Z = (z + z -1 ) (q + q -1 ) . (4.14)
For convenience, we denote respectively by Q M (Z), Q2s (Z) and Q2s (Z) the Q-polynomials associated with the symmetrized Bethe roots (3.103) satisfying (3.107) for a = sp, for a = d and κ * = 0, for a = d, respectively. Lemma 4.1. According to the parameters κ, κ * , the action of the Heun-Askey-Wilson q-difference operator π(I(κ, κ * , 0, 0)) on the corresponding Q-polynomial is such that:

π(I(0, κ * , 0, 0))Q M (Z) = Λ M sp,+ Q M (Z) , (4.15)
π(I(κ, 0, 0, 0)) Q2s (Z) = Λ * 2s sp,+ Q2s (Z)κqδ d (-1) 2s+1 (q + q -1 ) 2s H(Z) , (4.16) π(I(κ, κ * , 0, 0)) Q2s (Z) = Λ 2s d,+ Q2s (Z)κqδ d (-1) 2s+1 (q + q -1 ) 2s H(Z) . (4.17)

Proof. Firstly, we show (4.15). Consider the Baxter T-Q relation (4.1). Recall (3.37). Observe that:

Λ + 1 (u) = u -1 (u 2 -u -2 )(qu 2 -q -1 u -2 )φ(q -1 u -2 ) , (4.18) Λ + 2 (u) = u -1 (u 2 -u -2 )(q 2 u 2 -q -2 u -2 )φ(qu 2 ) (4.19)
with (4.12). Also, one easily show that:

(q + q -1 ) 2 (η + η * U ) ρ(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) = 1 2 q 1 2 (ν+ν ′ ) e -µ ′ q 2s + e µ ′ q -2s -φ(qu 2 ) -φ(q -1 u -2 ) .
Inserting these expressions in (4.1), one obtains (4.15) using (4.11) through the identification z = qu 2 . Secondly, we show (4.17). Inserting (4.18), (4.19) and

(q + q -1 ) 2 (η * + ηU ) ρ(u 2 -u -2 )(q 2 u 2 -q -2 u -2 ) = 1 2 q 1 2 (ν+ν ′ ) e -µ ′ q 2s (1 + q -2 )U + 2e µ ′ cosh(µ) -(v 2 + v -2 )q -2s-1 -q -2 u -2 φ(qu 2 ) -u 2 φ(q -1 u -2 )
into (4.7), one obtains (4.17) using (4.10), (4.11) through the identification z = qu 2 . The proof of (4.16) is the specialization κ * = 0 of the proof of (4.17).

Although not reported here, the connection between the Baxter T-Q relation for the generic case (4.8) and the Heun-Askey-Wilson q-difference operator introduced in [BTVZ18] can be established using a different realization of the Askey-Wilson algebra. In this case, the action of the Heun-Askey-Wilson q-difference operator π(I(κ, κ * , κ + , κ -)) on the Baxter Q-polynomial also takes a form similar to (4.17). 4.3. The Q-polynomial for the special case and the Askey-Wilson polynomials. It is well-known that the eigenfunctions of the second-order q-difference operator (4.11) are given by the Askey-Wilson orthogonal polynomials [KS96, eq. (3.1.6)]. Recall that the Askey-Wilson polynomials P M (x), M = 0, 1, 2, ... are defined by: P M (x; a, b, c, d) = 4 φ 3 q -2M , abcdq 2M-2 , az, az -1 ab, ac, ad ; q 2 , q 2 with x = z + z -1 . (4.20) Denoting the Askey-Wilson second-order q-difference by D according to [KS96, eq. (3.1.6)], the Askey-Wilson polynomials solve the bispectral problem:

DP M (x) = (q -2M + abcdq 2M-2 )P M (x), (4.21) xP M (x) = b M P M+1 (x) + a M P M (x) + c M P M-1 (x) (4.22)
where P -1 ≡ 0 and

b M = (1 -abq 2M )(1 -acq 2M )(1 -adq 2M )(1 -abcdq 2M-2 ) a(1 -abcdq 4M-2 )(1 -abcdq 4M ) , c M = a(1 -q 2M )(1 -bcq 2M-2 )(1 -bdq 2M-2 )(1 -cdq 2M-2 ) (1 -abcdq 4M-4 )(1 -abcdq 4M-2 ) , a M = a + a -1 -b M -c M .
From the three-term recurrence relation (4.22), one extracts the leading term of degree M in the variable x: P M (x) = (ab; q 2 ) M (ac; q 2 ) M (ad; q 2 ) M (abcdq -2 ; q 2 ) M a M (abcdq -2 ; q 4 ) M (abcd; q 4 ) M -1

x M + • • • (4.23)
For generic parameters q, a, b, c, d, the spectrum in (4.21) is non-degenerate, of the form (3.3). Based on this observation, a comparison of (4.21) and (4.23) with (4.15) immediately yields to: Proposition 4.5. For the special case κ = κ ± = 0, the Q-polynomial (4.2) of Proposition 4.1 is given by

Q M (Z) =
(ab; q 2 ) M (ac; q 2 ) M (ad; q 2 ) M (abcdq -2 ; q 2 ) M (q + q -1 ) M a M (abcdq -2 ; q 4 ) M (abcd; q 4 ) M P M (z + z -1 ; a, b, c, d) (4.24)

with a = -qe -µ+µ ′ , b = -qe µ+µ ′ , c = q -2s v 2 , d = q -2s v -2 . (4.25)
Expanding the Askey-Wilson polynomials in the variable Z = (z + z -1 )/(q + q -1 ) using (4.20) (see [START_REF] Koekoek | The Askey-scheme of hypergeometric orthogonal polynomials and its q-analogue[END_REF] for definitions), for the Q-polynomial (4.24) we obtain:

Q M (Z) = M l=0 (-1) l Q l,M Z M-l (4.26) with Q l,M =
(-1) M (ab; q 2 ) M (ac; q 2 ) M (ad; q 2 ) M (abcdq -2 ; q 2 ) M (q + q -1 ) l a M (abcdq -2 ; q 4 ) M (abcd; q 4 ) M (4.27) × M k=M-l (q -2M ; q 2 ) k (abcdq 2M-2 ; q 2 ) k q k(k+1) a k (ab; q 2 ) k (ac; q 2 ) k (ad; q 2 ) k (q 2 ; q 2 ) k e l+k-M (x 1 , x 2 , ..., x k )

where

x l = aq 2l-2 + a -1 q 2-2l .
Example 8. With the identification (4.25):

Q 1 (Z) = Z - (a + b + c + d -abc -abd -acd -bcd) (q + q -1 ) (1 -abcd) , Q 2 (Z) = Z 2 - q(a + b + c + d -q 2 (abc + abd + acd + bcd)) (1 -abcdq 4 ) Z +     q 2 -1 1 -a 2 b 2 c 2 d 2 q 4 + (ab + ac + ad + bd + bc + cd) q 2 + 1 + a 2 + b 2 + c 2 + d 2 q 2 + a 2 b 2 c 2 + a 2 b 2 d 2 + a 2 c 2 d 2 + b 2 c 2 d 2 q 4 + abcd (ab + ac + ad + bd + bc + cd) q 4 + q 6 -q 4 + q 2 + 1 abc 2 + ab 2 c + abd 2 + ab 2 d + acd 2 + ac 2 d + a 2 cd + a 2 bd + a 2 bc + b 2 cd + bc 2 d + bcd 2 -abcd q 2 + 1 3 + abc (a + b + c) + abd (a + b + d) + acd (a + c + d) + bcd (b + c + d)     (q + q -1 ) 2 (1 -abcdq 2 )(1 -abcdq 4 ) .
According to Proposition 4.5, the roots of the Q-polynomial can be computed for large values of M , a regime in which usually they are difficult to access by solving directly the Bethe equations of Proposition 3.2 (or similarly (3.107) for a = sp).

To conclude this subsection, let us remark that the identification of the roots of the Q-polynomial (4.26) with the admissible Bethe roots satisfying Proposition 3.2 (or (3.107) for a = sp) imply the existence of certain relations that are now described for completeness. For instance, given M fixed a comparison between Q M (Z) (see (4.2)) and (4.26) leads to a system of equations for the symmetrized Bethe roots (3.103). Equating the coefficients of both polynomials, one finds that the symmetrized Bethe roots (3.103) satisfy the following set of relations:

e l (U 1 , U 2 , ..., U M ) = Q l,M for l = 0, 1, • • • , M . (4.28)
For small values of M = 1, 2, 3, these relations are equivalent to the following polynomial equations.

Example 9.

For M = 1:

U 1 = Q 1,1 , For M = 2: U 2 = Q 1,2 -U 1 , U 2 1 -U 1 Q 1,2 + Q 2,2 = 0 , For M = 3: U 3 = Q 1,3 -U 1 -U 2 , U 2 2 + U 2 (U 1 -Q 1,3 ) + U 2 1 -U 1 Q 1,3 + Q 2,3 = 0 , U 3 1 -U 2 1 Q 13 + U 1 Q 2,3 -Q 3,3 = 0 .
Let us make some comments about the polynomial equations (4.28) compared with the polynomial equations (3.107). For small values of M , using (4.25) and (4.27) we have computed numerically the solutions of (4.28). It is found that they correspond to the subset of adminissible solutions of (3.107). A direct proof of this fact for generic values of M remains to be done.

As a second consequence of the relation between the Q-polynomial and the Askey-Wilson polynomials, an explicit relation between sets of Bethe roots associated with different values of M can be extracted from the recurrence relation (4.22). For convenience, introduce the notation

Q l,M ≡ e l (U [M] 1 , U [M] 2 , ..., U [M] M ) (4.29)
where the indices are added in order to distinguish between different root systems {U

[M] 1 , U [M] 2 , ..., U [M]
M } that solve (3.107). From (4.22) and (4.24) it follows:

Q 1,M+1 = Q 1,M + a M (q + q -1 ) , M > 0 Q l+2,M+1 = Q l+2,M + a M (q + q -1 ) Q l+1,M - c M b M-1 (q + q -1 ) 2 Q l,M-1 for l = 0, 1, • • • , M -2, M > 1 Q M+1,M+1 = a M (q + q -1 ) Q M,M - c M b M-1 (q + q -1 ) 2 Q M-1,M-1 .
4.4. The Q-polynomial as a transition coefficient. For the dynamical operators (A.2), consider the choice of gauge parameter β = 0. Using (4.13), after straightforward simplifications one gets:

π(B + (u, m)) = χb(u 2 ) α(q -q -1 )q 2+m u U - z + z -1 q + q -1 . (4.30)
Considering in particular the product of dynamical operators (3.44) entering in the definition of the Bethe states of Propositions 3.1, 3.2 and 3.3, the following lemma is easily shown. Lemma 4.2. Assume the gauge parameter α satisfies (3.21). For a = sp (M = 0, 1, ..., 2s) or a = d (M = 2s), the Baxter Q-polynomial (4.2) is given by

Q M (Z) = N M (ū) -1 z|Ψ M a,+ (ū) with N M (ū) = (-1) M (q + q -1 ) M 2 M q ν+ν ′ 2 e -µ ′ q 2s-M-1 M M i=1 b(u 2 i ) u i
for (3.107) and the following notation:

z|Ψ M a,+ (ū) = π(B + (u 1 , m 0 + 2(M -1)) • • • B + (u M , m 0 ))| β=0 1 . (4.31)
Here we have denoted z|Ψ 0 a,+ (ū) = z|Ω + = Q 0 (Z) ≡ 1 . For completeness, the action of the other dynamical operators on 1 is computed in a straightforward manner using T ± 1 = 1. Using the expressions reported in Appendix F, for which α satisfies (3.21) and β = 0, for any integer m 0 one finds12

π(A + (u, m 0 ))1 = u -1 b(u 2 )b(qu 2 )φ(q -1 u -2 ) , π(D + (u, m 0 ))1 = u -1 b(u 2 )b(q 2 u 2 )φ(qu 2 ) , π(C + (u, m 0 ))1 = 0 .
Suppose a rigorous mathematical definition of a basis {|z } within the framework of rigged Hilbert spaces (see e.g. [M01]) is given for the Askey-Wilson algebra, extending the case of the quantum harmonic oscillator and Hermite functions [START_REF] Celeghini | Applications of rigged Hilbert spaces in quantum mechanics and signal processing[END_REF] to the realm of Askey-Wilson orthogonal polynomials [KS96, p. 50]. To our knowledge, this problem has not been considered yet in the literature. If solved, then the notation z|Ψ M sp,+ (ū) would find a natural interpretation as a transition coefficient connecting the continuous basis {|z } and the discrete basis {|θ * M } given in Proposition 3.7.

Applications

In this Section, we apply the algebraic Bethe ansatz solution for the Heun-Askey-Wilson operator of Section 3 to the diagonalization of the q-analog of the quantum Euler top as well as to various examples of Hamiltonians of 3-sites Heisenberg spin chains in a magnetic field, inhomogeneous couplings, three-body terms and boundary interactions.

5.1. Algebraic Bethe ansatz solution for the q-analog of the quantum Euler top. The relation between the Hamiltonian of a quantum Euler top in a magnetic field built from sl 2 (R) and the Heun operator has been recently studied in [START_REF] Turbiner | The Heun operator as a Hamiltonian[END_REF] (see also [START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF]). It is a quantum version of the Zhukovsky-Volterra gyrostat of classical mechanics [START_REF] Basak | Explicit solution of the Zhukovski-Volterra gyrostat[END_REF][START_REF] Levin | Painlevé VI, rigid tops and reflection equation[END_REF], and arises in spin systems with anisotropy [START_REF] Zaslavskii | Periodic effective potentials for spin systems and new exact solutions of the onedimensional Schrödinger equation for the energy bands[END_REF]. Considering the realization of the Askey-Wilson algebra given in Example 1, it is straightforward to derive a q-deformed analog of the Euler top [Tu16, eq. ( 2)] generated from U q (sl 2 ) starting from (1.16). Replacing (2.12), (2.13) in (1.16) one gets the bilinear expression considered in [WZ95, Section 3.2]. It gives a Hamiltonian of the form: π(I(κ, κ * , κ + , κ -)) = t +-S + S -+ t 00 q 2s3 + t ′ 00 q -2s3 + t ++ S 2 + + t --S 2 -+t 0+ S + q s3 + t 0-S -q s3 + t ′ 0+ S + q -s3 + t ′ 0-S -q -s3 + I 0 , where the coupling constants t +-, t 00 , t ′ 00 , t ±± , t 0± , t ′ 0± are expressed in terms of the parameters k ± , ǫ ± , v, κ, κ * , κ ± , χ introduced in Example 1 and I 0 central in U q (sl 2 ).

For an irreducible finite dimensional representation of U q (sl 2 ) of dimension 2s + 1 on which this Hamiltonian acts, we consider the parametrization (2.11) and apply the results of Section 3. For generic parameters κ, κ * , κ ± , the spectrum and Bethe eigenstates are given by Proposition 3.4. In the Tables 2, 3 and 4, we give the numerical results for s = 1/2, s = 1, s = 3/2, respectively.

κ = 10i(1+ √ 5) √ 3 , κ * = 20i √ 3 , κ + = √ 3 2 , κ -= -3 2 χ = -40 3 √ 3 294
= 3, ν = ν ′ = 1, µ = 0.2, µ ′ = 0.3, v = 1.1.
κ = 10i(1+ √ 5) √ 3 , κ * = 20i √ 3 , κ + = √ 3 2 , κ -= -3 2 χ = -40 3
= 3, ν = ν ′ = 1, µ = 0.2, µ ′ = 0.3, v = 1.
κ = 10i(1+ √ 5) √ 3 , κ * = 20i √ 3 , κ + = √ 3 2 , κ -= -3 2 χ = -40 3 √ 3 -20905.2 +
= 3, ν = ν ′ = 1, µ = 0.2, µ ′ = 0.3, v = 1.1.
5.2.1. The Heisenberg chain for j 1 = j 2 = j 3 = 1/2. Denote S a i , a = x, y, z with i = 1, 2, 3 as the operators acting on the representation V (j i ) such that

S x = 1 2 0 1 1 0 , S y = 1 2 0 -i i 0 , S z = 1 2 1 0 0 -1 (5.1) Using π(q ±s3 ) = ± q -1 √ q S z + q + 1 2 √ q , π(S ± ) = S x ± iS y ,
it follows that the Heun-Askey-Wilson operator is given by π 6. Numerical results for the parameters q = 2, ν = ν ′ = 1.

(I(κ, κ * , κ + , κ -)) = 2(q -q -1 ) 2 a i,j J a ij S a i S a j + a =b i,j K ab ij S a i S b j + {a,b,c} i,j,k L abc ijk S a i S b j S c k + i b z i S z i + J 0 , ( 5

Perspectives

Besides the generalization to the Askey-Scheme of the basic quantum harmonic oscillator construction described in the Introduction, there are four main motivations for the present paper.

• The Askey-Wilson algebra with generators A, A * provides the algebraic framework for all orthogonal polynomials of the Askey-scheme [Z91]: the bispectral problem with respect to A, A * produces the well-known recurrence and second-order q-difference, difference or differential equations satisfied by these polynomials. If A, A * act as a Leonard pair (irreducible finite dimensional representation of the Askey-Wilson algebra are considered), the entries of the transition matrix between the two respective eigenbasis of the Leonard pair are given in terms of the orthogonal polynomials [START_REF] Terwilliger | Leonard pairs and the q-Racah polynomials[END_REF][START_REF] Terwilliger | Leonard pairs and the Askey-Wilson relations[END_REF]. By analogy, the Heun-Askey-Wilson algebra [BTVZ18, Definition 2.1] with generators A, I is a generalization of the Askey-Wilson algebra. So, it should provide the algebraic framework for a class of special functions beyond the Askey-scheme. Thus, investigating the spectral problem with respect to A, I is an important issue in this direction. Related works in these directions are e.g. [START_REF] Takemura | On q-deformations of the Heun equation[END_REF][START_REF] Kojima | Ultradiscrete limit of the spectral polynomial of the q-Heun equation[END_REF][START_REF] Takemura | Heun's differential equation and its q-deformation[END_REF].

• The q-Onsager algebra with generators W 0 , W 1 introduced in [T99] (see also [B04]) is known to be a homomorphic pre-image of the Askey-Wilson algebra with W 0 → A, W 1 → A * . From that point of view, the theory of tridiagonal pairs developed by Terwilliger et al. generalizes the theory of Leonard pairs [T99]. In [START_REF] Baseilhac | A new (in)finite dimensional algebra for quantum integrable models[END_REF][START_REF] Baseilhac | A deformed analogue of Onsager's symmetry in the XXZ open spin chain[END_REF], elements denoted {I 2k+1 , k ∈ Z + } that generate a commutative subalgebra of the q-Onsager algebra have been constructed. In general, they read as polynomials in W 0 , W 1 of maximal degree 2k + 2 [START_REF] Baseilhac | An attractive basis for the q-Onsager algebra[END_REF]. In particular, the image of the element I 1 in the Askey-Wilson algebra is the Heun-Askey-Wilson element (1.16). Thus, the analysis presented here can be viewed as a warm up for the diagonalization of the mutually commuting elements {I 2k+1 , k ∈ Z + } within the algebraic Bethe ansatz. In quantum integrable systems, it is important to stress that the elements {I 2k+1 , k ∈ Z + } are the basic building quantities for mutually commuting quantities, for instance the Hamiltonian of the open XXZ spin chain with generic integrable boundary conditions.

• The Askey-Wilson algebra admits an embedding into U q (sl 2 ) ⊗ U q (sl 2 ) ⊗ U q (sl 2 ) [START_REF] Granovskii | Hidden symmetry of the Racah and Clebsch-Gordan problems for the quantum algebra slq(2)[END_REF] (see also [H16]), where the generators map as A → ∆(C) ⊗ 1, A * → 1 ⊗ ∆(C) with C, ∆, respectively the Casimir element and coproduct of U q (sl 2 ). In the recent literature [PW17, DDV18], a generalization of the Askey-Wilson algebra indexed by N has been introduced. For N = 3, it reduces to the Askey-Wilson algebra. Importantly, an embedding of this algebra in terms of 'intermediate' Casimir elements of (U q (sl 2 )) ⊗N has been given. Thus, the spectral problem solved in the present paper can be viewed as a toy model for studying generalizations of (1.16) to N > 3. In the context of quantum integrable spin chains, this approach differs from the usual one based on Sklyanin's construction [S88]. Indeed, for N = 3 considering an irreducible finite dimensional representation for each U q (sl 2 ) labelled by j 1 , j 2 , j 3 , for j 1 = j 2 = j 3 = j, the Heun-Askey-Wilson operator associated with (1.16) gives the Hamiltonian of a three-sites spin-j chain in a magnetic field with inhomogeneous couplings, three-body and boundary interactions. Thus, generalizations of the Askey-Wilson algebra should naturally generate Heisenberg spin chains with possible inhomogeneous couplings and long range interactions.

• In the context of signal treatment, the optimal reconstruction of a signal from limited observational data is a central problem. In particular, the diagonalization of Heun type operators play a crucial role in the so-called band-time limiting problem [P87, [START_REF] Grünbaum | Time-band limiting and the bispectral problem[END_REF][START_REF] Gr Ünbaum | Algebraic Heun operator and band-time limiting[END_REF]. For a recent review see [START_REF] Bergeron | Signal processing, orthogonal polynomials, and Heun equations[END_REF] and references therein. For the q-deformed case, it is straigthforward to derive the q-analogs of the conditions [GVZ17, eqs. (5.31), (5.32), (5.33)] such that π(I(κ, κ * , κ + , κ -) commutes with projectors on certain eigenspaces of π(A), π(A * ). Thus, the results here presented apply in this context.

Other perspectives may be also briefly mentionned. For instance, for the realization (4.10), (4.11), the Heun-Askey-Wilson operator gives a fourth-order q-difference operator. Possible connections with q-Krall polynomials [VYZ01] could be explored. Another interesting direction is to study the relation with the separation of variables framework, see e.g. [START_REF] Faldella | Complete spectrum and scalar products for open spin-1/2 XXZ quantum chains with non-diagonal boundary terms[END_REF] and references therein. In particular, the eigenvectors |θ ⋄ M introduced in subsection 3.1.1 diagonalize the off-diagonal entry π(B(u)) given by (2.6). In addition, it should be possible to study the spectrum of the Heun-Askey-Wilson q-difference operator using an homogeneous Baxter T-Q relation with nonpolynomial solution [START_REF] Lazarescu | Bethe Ansatz and Q-operator for the open ASEP[END_REF]. Also, for q a root of unity and cyclic representations of the Askey-Wilson algebra, the results here presented may be extended in light of [START_REF] Baseilhac | Cyclic tridiagonal pairs, higher order Onsager algebras and orthogonal polynomials[END_REF][START_REF] Huang | Finite-dimensional irreducible modules of the universal Askey-Wilson algebra at roots of unity[END_REF]. A + (u, m) = b(u 2 ) u(αq 2mq 2 β) χq m ρ [A * , A] q -αβq m+2 χ [A, A * ] q + (β + αq 2m )qA -(αq 2m + βq 4 u 4 ) qu 2 A * -1 ρq 3 u 2 χb(q 2 )b(u 2 ) ρq m+2 (αβq 2 ρχ 2 )(u 6 q 2u -2 ) +(q 2 + 1)q m (αη * χq m (q 4 -1)q 2 ω(αβq 2 ρχ 2 )) -(q 2 + 1)q 2 u 4 (βη * χ(q 4 -1)q m ω(αβq 2 ρχ 2 ))

+(q 2 -1)u 2 (ηχ(αq 2mβq 2 )(1 + q 2 ) 2ρq m+2 (αβq 2 ρχ 2 )) ,

A -(u, m) = b(u 2 ) αq 2 -βq 2m uχq m+2 ρ [A * , A] q - αβuq m χ [A,
A * ] q + (βq 2m + αq 4 u 4 ) qu Aqu(α + βq 2m )A * + 1 ρq 3 u 2 χb(q 2 )b(u 2 ) ρq m+2 (q 2 u 7u -1 )(q 2 χ 2αβρ) -(q 2 -1)u 3 ((q 2 + 1) 2 η * χ(αq 2βq 2m ) +ρq m+2 (q 2 χ 2αβρ)) -(q 2 + 1)q 2 u 5 (ωq m (αβρq 2 χ 2 ) + α(q 4 -1)ηχ) +(q 2 + 1)uq m (β(q 4 -1)ηχq m + ω(αβq 2 ρq 4 χ 2 )) ,

B + (u, m) = βb(u 2 ) αq 2m+2 -β χq m βρu [A * , A] q - βq -m uχ [A, A * ] q + (q 2 + 1) qu A -( 1 qu 3 + qu)A * (A.2) +
q -m-3 βρχb(q 2 ) (q 2 + 1)u -1 (ω(χ 2 q 2m+2β 2 q 2 ρ) + β(q 4 -1)η * χq m )q 2 ρ(q 2 u + u -3 )(β 2 ρχ 2 q 2m ) ,

B -(u, m) = βb(u 2 )q 2m+1 αq -2βq 2m uχq -m-1 βρ [A * , A] q -βuq m-1 χ [A, A * ] q + (q 2 u 4 + 1) q 2 u A -(q 2 + 1)u q 2 A * -q -m-4 βρχb(q 2 ) ρq 2 (q 2 u 3 + u -1 )(β 2 ρq 2mχ 2 ) + (q 2 + 1)u(ω(β 2 ρq 2m+2q 2 χ 2 ) + β(q 4 -1)ηχq m ) ,

C + (u, m) = αb(u 2 ) α -βq 2-2m
αq m uχ [A, A * ] q -χq -m αρu [A * , A] q -(q 2 + 1) qu A + ( 1 qu 3 + qu)A * (A.3) + q -m-3 αρχb(q 2 ) ρq 2 (q 2 u + u -3 )(α 2 ρq 2mχ 2 ) + (q 2 + 1)u -1 (ω(α 2 ρq 2m+2q 2 χ 2 )α(q 4 -1)η * χq m ) ,

C -(u, m) = αqb(u 2 ) αq 2 -βq 2m αuq 1-m χ [A, A * ] q - uχq m+1 αρ [A *
, A] q -(q 2 u 4 + 1) u A + u(q 2 + 1)A * + q -m αρχ(q 4 -1) ρq 2 (q 2 u 3 + u -1 )(α 2 ρχ 2 q 2m ) + (q 2 + 1)u(ω(α 2 q 2 ρχ 2 q 2m+2 ) + α(q 4 -1)ηχq m ) , D + (u, m) = b(u 2 )b(q 2 u 2 ) ub(qu 2 )(αq 2mβq 2 ) αβq m+2 χ [A, A * ] q -χq m ρ [A * , A] qq(β + αq 2m )A + q(β + αu 4 q 2m ) u 2 A * (A.4) + 1 ρq 3 u 2 b(q 2 )b(q 2 u 2 ) (q 2 -1)u 2 χ -1 (ρq m+2 (αβq 2 ρχ 2 ) + ηχ(β(q 3 + q) 2α(q 2 + 1) 2 q 2m ))

+ρχ -1 q m (q 6 u 6u -2 )(αβq 2 ρχ 2 ) + (q 2 + 1)χ -1 (ωq m (χ 2αβq 2 ρ) + β(q 4 -1)η * χ) -(q 2 + 1)u 4 q m+2 χ -1 (α(q 4 -1)η * χq m + q 2 ω(χ 2αβq 2 ρ)) ,

One has b(qu 2 i ) N = g (-1) N/2-k N k ((q + q -1 )U i ) N -2k +(N -2k)

(N -2k)/2 l=1 (-1) l l N -2kl -1 l -1 ((q + q -1 )U i ) N -2k-2l + N N/2   .

Collecting terms of same degree in U i , the claim follows.

Let us introduce ∆ d (u) and ∆ g (u) respectively given by (3.76) and (3.85), and define ∆ sp (u) = 1. Define F ǫ (u i ) = e -(µ-µ ′ )/2+ǫ(µ-µ ′ )/2 q -(2s+1) (q 2s+1 u i v -1u -1 i v)(q 2s+1 u i vu -1 i v -1 ) ×(e (µ-µ ′ )(1-ǫ)/2 u i + e (µ+µ ′ )(1+ǫ)/2 u -1 i )(e (µ+µ ′ )(1-ǫ)/2 u i + e (µ ′ -µ)(1+ǫ)/2 u -1 i ).

Lemma D.3. There exist polynomials {g

- ∆ a (u i ) b(qu 2 i ) F ǫ (u i ) + ∆ a (q -1 u -1 i ) b(qu 2 i ) F ǫ (q -1 u -1 i ) = g (even)
a,ǫ b(qu j /u i )b(u i u j )F ǫ (u i ) + ∆ a (q -1 u -1 i ) b(qu 2 i ) M j=1,j =i b(qu i /u j )b(q 2 u i u j )F ǫ (q -1 u -1 i ) = M-1 k=0 (-1) k (q + q -1 ) M-1 2 M-1-k e k ( Ūi )

(U i ) , (D.2) ∆ a (u i )F ǫ (u i ) + ∆ a (q -1 u -1 i )F ǫ (q -1 u -1 i ) = g (odd
M-1-k l=0 M -1 -k l (q -q -1 ) l U M-1-k-l i (q 2 + q -2 ) 1+k+l-M g 2[ l 2 ] 0 (U i )g (p[l]) a,ǫ (U i ) ,
where p[l] = even (resp. odd) for l even (resp. odd). (D.4) Also, we will need: Lemma D.4. For any integer or half-integer s, define: H(U i ) = (q + q -1 ) 2s+1 2s+1 k=0 (-1) k e k (X 0 , X 1 , ..., X 2s )U 2s+1-k i with X k = q 2k-2s v 2 + q -2k+2s v -2 q + q -1 .

One has: 2s k=0 b(q 1/2+k-s vu)b(q 1/2+k-s v -1 u) = H(U ) . (D.5) Proof. Observe that the product (D.5) is left invariant under the transformation (3.101). Thus, it is a polynomial in U i . In terms of the elementrary symmetric polynomials (1.22), one routinely gets (D.5). constants {ρ, ω, η, η * } given by (E.3). In particular, one finds that π 1 (A), π 1 (A * ) and π(A), π(A * ) given by (4.10), (4.11) are related through the transformation (2.28). For completeness, in addition to (4.10), (4.11), (4.13), one finds: (E.6) π([A, A * ] q ) = b(q 2 )φ(z)φ(q 2 z)q -2 z -1 T 2 + + b(q 2 )φ(z -1 )φ(q 2 z -1 )q -2 zT 2 -+ b(q)q (ν+ν ′ )/2 2 (2 cosh(µ)q 2se -µ ′ q -1 (v 2 + v -2 )) + b(q)q (ν+ν ′ )/2 2 (e -µ ′ q 2s-2 (q + q -1 ) + e µ ′ q -2s-1 + e -µ ′ q 2s-1 )z -1 -b(q 2 )q -2 z -1 φ(q 2 z)b(q)q -1 z -1 φ(z)b(qz)q -1 φ(z -1 ) + b(z)φ(q -2 z -1 ) φ(z)T + + b(q)q (ν+ν ′ )/2 2 (2 cosh(µ)q 2se -µ ′ q -1 (v 2 + v -2 )) + b(q)q (ν+ν ′ )/2 2 (e -µ ′ q 2s-2 (q + q -1 ) + e µ ′ q -2s-1 + e -µ ′ q 2s-1 )z -b(q 2 )q -2 zφ(q 2 z -1 )b(q)q -1 zφ(z -1 )b(qz -1 )q -1 φ(z) + b(z -1 )φ(q -2 z) φ(z -1 )T - + b(q)q ν+ν ′ 4 (e µ ′ q -2s + e -µ ′ q 2s )(2 cosh(µ)q 2s + e -µ ′ q 2s-1 (z + z -1 )e -µ ′ q -1 (v 2 + v -2 )) -b(q)q (ν+ν ′ )/2 2 ((2 cosh(µ)q 2se -µ ′ q -1 (v 2 + v -2 )) + (e µ ′ q -2s-1 + 2e -µ ′ q 2s-1 )z -1 + e -µ ′ q 2s-1 z)φ(z) -b(q)q (ν+ν ′ )/2 2 ((2 cosh(µ)q 2se -µ ′ q -1 (v 2 + v -2 )) + (e µ ′ q -2s-1 + 2e -µ ′ q 2s-1 )z + e -µ ′ q 2s-1 z -1 )φ(z -1 ) +b(q)q -1 z -1 φ(z) 2 + b(q)q -1 zφ(z -1 ) 2 + b(q)q -1 (z + z -1 )φ(z)φ(z -1 ) +b(z)(φ(z -1 )φ(q -2 z)φ(z)φ(q -2 z -1 )) .

Appendix F. Dynamical operators and the q-difference realization

Fixing the gauge parameters β = 0 and α by (3.21), using (4.10), (4.11), (4.13) and (E.6), we obtain the following expressions for the dynamical operators (A.1), (A.3) and (A.4) : π(A + (u, m 0 )) = -b(u 2 )q -1 u -3 z -1 φ(z)(zqu 2 )T +b(u 2 )q -1 u -3 zφ(z -1 )(z -1qu 2 )T - + q (ν+ν ′ )/2 (q + q -1 ) b(qz)b(qz -1 ) (U -z + z -1 q + q -1 ) -1 2 e -µ ′ b(u 2 )q 2s u -1 (z 2 + z -2 ) +((q -2s-1 + q 2s+1 )(cosh(µ)u -3 + cosh(µ ′ )u -1 )u -3 (v -2 + v 2 )(u 2 cosh(µ) + cosh(µ ′ )))(z + z -1 ) -(q + q -1 )(v -2 + v 2 )(cosh(µ)u -3 + cosh(µ ′ )u -1 ) + cosh(µ)(q -2s-1 + q 2s+1 )(q + q -1 )u -1 +u -3 ( 1 2 e -µ ′ ((q -2sq 2s )q -2 + q -2s + q 2s ) + 1 2 e µ ′ (q -2s-2 + q -2s + q 2s + q 2s+2 ))

+ 1 2 e -µ ′ q 2s (q 2 + q -2 )u , π(D + (u, m 0 )) = -b(q 2 u 2 )b(u 2 ) b(qu 2 ) φ(z)(u -1 z -1qu)T + -b(q 2 u 2 )b(u 2 ) b(qu 2 ) φ(z -1 )(u -1 zqu)T - + q (ν+ν ′ )/2 (q + q -1 )b(u 2 ) b(qz)b(qz -1 )b(qu 2 ) (U -z + z -1 q + q -1 ) 1 2 e -µ ′ b(q 2 u 2 )q 2s u -1 (z 2 + z -2 ) +((q -2s-1 + q 2s+1 )(cosh(µ)q 2 + cosh(µ ′ )u -1 )u -1 (v -2 + v 2 )(cosh(µ) + q 2 u 2 cosh(µ ′ )))(z + z -1 ) -(q + q -1 )(v -2 + v 2 )(cosh(µ)q 2 u + cosh(µ ′ )u -1 ) + cosh(µ ′ )(q -2s-1 + q 2s+1 )(q + q -1 )q 2 u +u -1 ( 1 2 e -µ ((q -2sq 2s )q -2 + q -2s + q 2s ) + 1 2 e µ (q -2s-2 + q -2s + q 2s + q 2s+2 ))

-1 2 q 2s (q 2 + q -2 )(e -µ ′ b(q 2 u 2 )e -µ )u -1 , + (e -2µ ′ q 4s + 1)z (q + q -1 )uq + e 2µ ′ q -4s + 2e -2µ ′ q 4sq 2 + 1q -2 (q + q -1 )uzq + q -1 u -3 φ(z) b(u 2 ) .

  , a † |θ * n = 1 √ 2 |θ * n+1 and N |θ * n = n|θ * n where |θ * 0 denotes the lowest weight (or 'reference') state such that a|θ * 0 = 0. In terms of the operators (H, X), one has: (X -[X, H]) n |θ * 0 = 0 and |θ * n = (X + [X, H]) n |θ * 0 .

  .28) Using (3.1), (3.2), one finds that the Askey-Wilson relation (1.15) implies (3.5) and (3.27). Thus, const 0 (u) = 0. The second relation in (3.22) is shown similarly.

  2. For convenience, the parametrization (3.32)-(3.36) is used. Recall the notations (3.44), (3.45). Proposition 3.1. Define |Ψ M sp,-(ū, m 0 ) = π(B -(ū, m 0 , M ))|Ω -. (3.57) One has:

  |Ψ 2s g,ǫ (ū, m 0 ) = π(B ǫ (ū, m 0 , 2s))|Ω ǫ , (3.90) |Ψ 2s g,ǫ ({u, ūi }, m 0 ) = π(B ǫ ({u, ūi }, m 0 , 2s))|Ω ǫ . (3.91) Conjecture 1. For M = 2s and generic {u, u i }, one has:

  104) with (D.1), (D.2), (D.3), (D.4), (D.5) and ∆a

4.1. 3 .

 3 Diagonal case κ, κ * = 0 and κ ± = 0. Following a similar analysis, for the diagonal case it is straightforward to derive the inhomogeneous Baxter T-Q relation from (3.77). It yields to: Proposition 4.3. The eigenvalues Λ 2s d,+ of the Heun-Askey-Wilson operator π(I(κ, κ * , 0, 0)) are given by the inhomogeneous Baxter T-Q relation

Appendix A .

 . Dynamical operators in the Askey-Wilson presentation The dynamical operators are polynomials of maximum degree two in the elements A, A * of the Askey-Wilson algebra. From (3.16)-(3.19), one finds: (A.1)

  i ). Proof. Using [GR43, eqs. 1.320, 1.331], one has:

  (U i )|p ∈ {even, odd}, a ∈ {sp, d, g}, ǫ = ±} of the form g (p) a,ǫ (U i ) = ,[k] are scalars such that

  ) a,ǫ (U i ) . (D.3) Proof. Routine. As a consequence of Lemmas D.1, D.2, D.3 and Remark 2, it follows: Corollary D.1.(-1) M ∆ a (u i ) b(qu 2 i ) M j=1,j =i

  2s, one gets the Bethe equations (3.75).

	To determine the eigenvalues
	Λ 2s d,± in (3.72), observe that

Table 2 .

 2 Numerical results for the parameters q

	.909 -337.018 i -1006.7 + 746.697 i	294.909 -337.018 i -1006.7 + 746.697 i	{0.454184 + 0.509566 i} {1.36299 -0.412627 i}

Table 3 .

 3 Numerical results for the parameters q

			-2394.67 + 986.732 i	-2394.67 + 986.732 i	{2.98826 -0.846233 i, -0.155658 + 1.20672 i}
	√	3	-6079.21 + 1505.54 i -1543.12 -1249.58 i	-6079.21 + 1505.54 i -1543.12 -1249.58 i	{4.06015 + 0.244047 i, 1.69724 -0.997537 i} {2.43438 + 1.09148 i, 0.117738 + 1.26215 i}

Table 4 .

 4 Numerical results for the parameters q

	2621.26 i	-20905.2 + 2621.26 i	{1.04709 -1.16555 i, 8.06887 -1.51037 i
			9.48236 + 1.52911 i}
	-15139.6 + 3210.92 i	-15139.6 + 3210.92 i	{-0.255185 + 3.99139 i, 3.64641 -1.51464 i
			11.0049 -0.370951 i}
	-9508.34 -2492.84 i	-9508.34 -2492.84 i	{0.658762 + 0.810534 i, 4.03758 + 7.37438 i,
			9.56462 + 1.26855 i}
	-6679.41 + 393.291 i	-6679.41 + 393.291 i	{2.80458 + 9.49603 i, 0.718949 + 0.525899 i,
			7.03427 -0.69735 i}

Table

	Spin-1 chain	Direct diagonalization (degeneracy) 128.125 (7)	Diagonalization via ABA Λ 2s a,+ (s) 128.125 (0)	Bethe roots {U 1 , . . . , U 2s } -
	κ = 1 , κ * = 1/3 , κ ± = 0	106.128 (5) 54.4972 (5) 100.033 (3) 45.8128 (3) 24.7794 (3) 32.5 (1)	106.128 (1/2) 54.4972(1/2) 100.033 (1) 45.8128 (1) 24.7794 (1) 32.5 (0)	{-3.82841} {-3.39815} {-5.84346, -1.26701} {-4.28657,-1.01656} {-3.59681, -1.00521} -
		230.13 (7)	230.13 (0)	-
	κ = -5 4 √ κ + = 1 2 , κ * = -9 4 √ 2 , 8 , κ -= -1 16 4 χ = -15	33.5628 + 33.3326 i (5) 33.5628 -33.3326 i (5) -2.26847 + 22.221 i (3) -2.26847 -22.221 i (3) -13.7437 (3)	33.5628 + 33.3326 i (1/2) 33.5628 -33.3326 i (1/2) -2.26847 + 22.221 i (1) -2.26847 -22.221 i (1) -13.7437 (1)	{ -3.51786 -5.42098 i} {-3.51786 + 5.42098 i} {-8.41833 -8.60097 i, -0.252667 -0.78235 i} {-8.41833 + 8.60097 i, -0.252667 + 0.78235 i} {-1.91267 + 5.49269 i, -1.91267 -5.49269 i}
		-0.200512 (1)	-0.200512 (0)	-
	.2)

However, adding a term such as {X, H} to (1.13) where {., .} is the anticommutator leads to a more complicated spectral problem.

The parameter χ is introduced for further convenience.

For a certain image of (1.16) in Uq(sl 2 ), see however[START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF].

For a more complete list of references and an account of historic details on this subject, we refer the reader to the Introduction of[START_REF] Belliard | Modified algebraic Bethe ansatz for XXZ chain on the segment -I: Triangular cases[END_REF].

The on-shell Bethe states can also be retrieved from the off-diagonal Bethe ansatz[START_REF] Zhang | Retrieve the Bethe states of quantum integrable models solved via off-diagonal Bethe Ansatz[END_REF].

Note that for the generic case, the eigenvalue expression can be obtained from the results in[START_REF] Cao | Exact solution of the alternating XXZ spin chain with generic non-diagonal boundaries[END_REF].

Note that if the constraints π(B ± (u, m 0 ))|Ω ± = 0 are satisfied then the excited states would be created by successive actions of the dynamical operators π(C ± (u, m)) on |Ω ± .

Similar relations have been obtained in[START_REF] Crampé | Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries[END_REF][START_REF] Belliard | Modified algebraic Bethe ansatz for XXZ chain on the segment -I: Triangular cases[END_REF].

Note that after acting on 1 with the expressions in Appendix F the dependence on z vanishes and the results drastically simplify.

Acknowledgments: We thank S. Belliard, O. Brodier, N. Crampé and P. Terwilliger for discussions and comments. We thank H-W. Huang for explanations on [H16]. P.B. thanks H. Saleur for pointing out reference [F19], C. Fewster for pointing out reference [M05] and comments on rigged Hilbert spaces, R. de la Madrid for comments on the construction of an '|x ′ basis for the Askey-scheme using rigged Hilbert space, G. Roux for comments on potential applications of the Hamiltonians of Section 5, and L. Vinet and A. Zhedanov for joint collaborations around the Heun operator and comments on applications of the results here presented to the time-band limiting problem in signal processing. R.P. thanks the Institut Denis Poisson for hospitality. R.P. is supported by the São Paulo Research Foundation (FAPESP) and by the Coordination for the Improvement of Higher Education Personnel (CAPES), grants # 2017/02987-8 and #88881.171877/2018-01. P.B. is supported by C.N.R.S.

with the coupling constants, J x 12 = J y 12 = κ + 3q 4 -2q 2 + 2q -2 -3q -4 4 (κ + χ -1 + κ -χ), J x 23 = J y 23 = κ * + 3q 4 -2q 2 + 2q -2 -3q -4 4 (κ + χ -1 + κ -χ), J z 12 = ∆κ + qq -1 2 (q 2 + q -2 ) 2 (κ + χ -1 + κ -χ), J z 23 = ∆κ * + qq -1 2 (q 2 + q -2 ) 2 (κ + χ -1 + κ -χ), J x 13 = J y 13 = J z 13 = (qq -1 ) 3 2 (κ + χ -1 + κ -χ),

b z 1 = qq -1 4 (κ + (q 3 + q -3 )(qq -1 )(κ + χ -1 + κ -χ)), b z 3 = -qq -1 4 (κ * + (q 3 + q -3 )(qq -1 )(κ + χ -1 + κ -χ)), J 0 = (q + q -1 + 3q 3 + 3q -3 ) (κ + κ * ) 4

+(12q -12q -1 -5q 3 + 5q -3 + q 5q -5 + 2q 7 -2q -7 ) (κ + χ -1 + κ -χ) 4 .

For the diagonal case κ ± = 0, the three-body terms vanish and the Heun-Askey-Wilson operator simplifies to π(I(κ, κ * , 0, 0)) = (qq -1 ) 2 2 κ(S x 1 S x 2 + S y 1 S y 2 + ∆S z 1 S z 2 ) + κ * (S x 2 S x 3 + S y 2 S y 3 + ∆S z 2 S z 3 ) -qq -1 2 (κ * S z 3 -κS z 1 + (κκ * )S z 2 ) + (q + q -1 ) 2 (κ + κ * ) 3 2 (q 2 + q -2 ) -1 .

Following Example 4, we have set j 1 = j 2 = j 3 = 1/2. Then the set Σ = (ℓ, k) such that:

In the algebraic Bethe ansatz results of Proposition 3.2 and 3.3, we fix:

In Table 5, for numerical values of the scalar parameters κ, κ * , κ ± the eigenvalues of the Hamiltonian π(I(κ, κ * , κ + , κ -)) obtained by direct diagonalization and the ones obtained from the algebraic Bethe ansatz (ABA) are displayed. 5.2.2. The Heisenberg chain for j 1 = j 2 = j 3 = 1. Denote S a i , a = x, y, z with i = 1, 2, 3 as the operators acting on the representation V (j i ) such that

Spin-

-0.200512 (4) -6.25895 + 3.32745 i (2) -6.25895 -3.32745 i (2) -0.200512 (0) -6.25895 + 3.32745 i (1/2) -6.25895 -3.32745 i (1/2) -{-0.793147 -1.40509 i} {-0.793147 + 1.40509 i} Table 5. Numerical results for the parameters q = 2, ν = ν ′ = 1.

Note that the term (S z ) 2 will lead to a proliferation of higher-order terms in the basis {S x i , S y i , S z i }, and for that reason we only write the expression of the Heun-Askey-Wilson operator for the diagonal case It is given by

Following Example 4, we have set j 1 = j 2 = j 3 = 1, and the set Σ = (ℓ, k) such that:

In the algebraic Bethe ansatz results of Proposition 3.2 and 3.3, we fix:

In Table 6, for numerical values of the scalar parameters κ, κ * , κ ± the eigenvalues of the Hamiltonian π(I(κ, κ * , κ + , κ -)) obtained by direct diagonalization and the ones obtained from the algebraic Bethe ansatz (ABA) are displayed.

+ρq m (q 6 u 7u -1 )(q 2 χ 2αβρ) + (q 2 + 1)u(ωq m (αβρq 2 χ 2 ) + α(q 4 -1)ηχ)

-(q 2 + 1)u 5 q m+2 (β(q 4 -1)ηχq m + ω(αβq 2 ρq 4 χ 2 )) .

Appendix B. Coefficients of the commutation relations

The coefficients of the dynamical commutation relations (3.41), (3.42), (3.43) are given by

Appendix C. Proof of Lemma 3.5

In this section, we show Lemma 3.5 for the choice ǫ = +. The proof for ǫ =is done along the same line, so we omit the details. As a preliminary, consider (A.2) for β = 0. It reads:

where the notation Ã⋄

with (2.25) is introduced for convenience. Recall the notation (3.44), (3.45). Successively, one has:

Also, one has:

Consider now the combination:

In this expression, the gauge parameter α is given by (3.21) with (3.32), and we use (3.71) for ǫ = +. Using the previous expressions, after straightforward simplifications of ( * ) and factorizing out the common overall factor, one gets the new combination:

.

The combination ( * * ) is a polynomial in Ã⋄ of degree 2s with coefficients that are meromorphic functions of U, U i . It has poles located at U = U i and U i = U j , and the residues at these points vanishes. Setting U = C, one finds:

We now consider ( * * ) on the finite dimensional representation V. Recall (C.2) with (2.25). Using (3.4) with (3.32), observe that X k is the spectrum of Ã⋄ . Thus, ( * * ) is proportional to the characteristic polynomial of A ⋄ which is vanishing on V. This concludes the proof of Lemma 3.5 for ǫ = +.

Appendix D. Proof of Proposition 3.5

To prepare the proof of Proposition 3.5, we first derive some intermediate results. Recall the notation Ūi = {U 1 , ..., U i-1 , U i+1 , ..., U M }, (3.103), (1.21) and (1.22).

Proof. Observe that:

Expand this expression in terms of the elementary symmetric polynomials in the variables U j , j = i. Insert:

Then, expand in U i and b(qu

Lemma D.2. For any even integer N , define:

We are now ready to show Proposition 3.5. First, we extract from E a (u i , ūi ) the part that is left invariant under the transformations (3.101), (3.102). Introduce the rational functions in u i , i = 1, ..., M :

It is easy to show that: Appendix E. second-order q-difference operators and the Askey-Wilson algebra Below, we display two examples of second-order q-difference operator realizations of the Askey-Wilson algebra. The realization π 1 is taken from [T99], whereas π 2 is taken from [START_REF] Baseilhac | The q-Onsager algebra and multivariable q-special functions[END_REF]. Using the invertible non-linear transformation (2.28) that relates different Askey-Wilson algebras, in particular it is shown that the realization π in (4.10), (4.11) follows from π 1 . In each case, the structure constants are given.

Adapting the notations of [T99] for our purpose, let us consider the following linear transformation denoted

where T ± is such that T ± (f (z)) = f (q ±2 z) and (4.12). Inserting (E.1), (E.2) in (1.14), (1.15), the structure constants of the Askey-Wilson algebra are given by: ρ = -1 4 q ν+ν ′ (q 2q -2 ) 2 , (E.3) ω = 1 2 (qq -1 ) 2 q ν+ν ′ cosh(µ)(q 2s+1 + q -2s-1 )cosh(µ

The structure constants in this case are given by (3.33)-(3.36). This realization of the Askey-Wilson algebra has been used in [START_REF] Baseilhac | The q-Onsager algebra and multivariable q-special functions[END_REF], based on [START_REF] Wiegmann | Algebraization of difference eigenvalue equations related to Uq(sl 2 )[END_REF].

Using the invertible non-linear transformation (2.28) that relates two different Askey-Wilson algebras, other realizations can be derived. Let A, A * satisfy the Askey-Wilson with structure constants (3.33)-(3.36). Then, using -q ν+ν ′ -4 (q 4 -1)

2 → i q 4 -1 q 1 2 (ν+ν ′ -4) one finds that Ā, Ā * satisfy (2.29), (2.30) with the structure

) + e -2µ ′ q 2s-2q -2s ) q 2q -2 -2e µ ′ φ(q 2 z -1 )q 1 2 (-ν-ν ′ )-2s-2 + e 2µ ′ q -4s-1 + e -2µ ′ q 4s-1q + q -1 q + q -1 z + e -2µ ′ q 4s z -1 b(u 2 )u -1 φ(z -1 )T - + -2e µ ′ (-q -2s + e -2µ ′ q 2s-2 + (q -2se -2µ ′ q 2s+2 )z 2 )φ(q -2 z -1 )φ(z)q (ν+ν ′ )/2 (q 2q -2 )uz + 2e µ ′ ((e -2µ ′ q 2s+2q -2s )z -1 + (q -2se -2µ ′ q 2s-2 )z)φ(z -1 )φ(q -2 z)q (ν+ν ′ )/2 (q 2q -2 )u e -3µ ′ (e 2µ ′ q 2 -1)(e µ ′q 2s+1 )(q 2s+1 + e µ ′ )(v 2 + v -2 )q ν+ν ′ 2 -3 2(q + q -1 )u e -2µ ′ (q 2sq -2s )(e 2µ ′q 4s+2 ) cosh(µ)q ν+ν ′ 2 -2 (q + q -1 )u + 2e -µ ′ (q 4s + e 2µ ′ )zφ(z -1 ) 2 q -(ν+ν ′ )/2-2s-1 (q + q -1 )u + 2e -µ ′ (q 4s + e 2µ ′ )φ(z) 2 q -(ν+ν ′ )/2-2s-1 (q + q -1 )uz + 2e -µ ′ (q 4s + e 2µ ′ )(z + z -1 )φ(z -1 )φ(z)q -(ν+ν ′ )/2-2s-1 (q + q -1 )u + e -3µ ′ (-e 2µ ′ (q 4 + 1)q 4s + q 8s+2 + e 4µ ′ q 2 )(z + z -1 )q ν+ν ′ 2 -2s-3

2(q + q -1 )u

--e -2µ ′ (q 4s + e 2µ ′ )(v 2 + v -2 )q -2s-1 (q + q -1 )u + uq + e -µ-µ ′ (1 + e 2µ )(q 4s + e 2µ ′ ) (q + q -1 )u + (e 2µ ′ q -4s + 2e -2µ ′ q 4sq 2 + 1q -2 )z (q + q -1 )uq + e -2µ ′ q 4s + 1 (q + q -1 )uzq + q -1 u -3 φ(z -1 )

--e -2µ ′ (q 4s + e 2µ ′ )(v 2 + v -2 )q -2s-1 (q + q -1 )u + uq + e -µ-µ ′ (1 + e 2µ )(q 4s + e 2µ ′ ) (q + q -1 )u