Diagonalization of the Heun-Askey-Wilson operator, Leonard pairs and the algebraic Bethe ansatz

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: 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]: giving a second presentation [Z91]. Alternatively, for the pair (H, iP) one also gets Askey-Wilson relations, with different structure constants: 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|θ * n = √ 2n|θ * n−1 , 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: n |θ * 0 = 0 and |θ * n = (X + [X, H]) n |θ * 0 . (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 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: Alternatively, it is also known that the spectral problem (1.6) can be formulated within the analytical Bethe ansatz framework, see e.g. [Gr17, Section 2.1]. Namely, consider the elementary Baxter T-Q relation: Given n fixed, (1.9) implies that the zeroes {x i |i = 1, ..., n} satisfy the set of Bethe ansatz equations .., 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) . 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 problem 2 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 [BTVZ18] 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 . 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 [WZ95,DE19]. Define the Heun-Askey-Wilson element [BTVZ18,CVZ19]: where κ, κ * , κ ± and 3 χ = 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 generality 4 . 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 [Za95,B04] associated with non-diagonal reflection matrices [DeG93,GZ94], 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, see 5 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 [CYSW13], 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 states 6 , which in particular allows the computation of scalar products between on-shell/off-shell Bethe states, see [BP15a,BP15b,BS19a,BS19b] 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 [T03,TV03] 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 problem 7 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 3 The parameter χ is introduced for further convenience. 4 For a certain image of (1.16) in Uq(sl 2 ), see however [WZ95]. 5 For a more complete list of references and an account of historic details on this subject, we refer the reader to the Introduction of [Be15]. 6 The on-shell Bethe states can also be retrieved from the off-diagonal Bethe ansatz [ZLCYSW15]. 7 Note that for the generic case, the eigenvalue expression can be obtained from the results in [CYSW15].
degree M follows from the 'transition' coefficient (see Lemma 4.2): 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): For the generic case κ, κ * , κ ± = 0 and a different realization of the Heun-Askey-Wilson operator, a similar picture holds as briefly mentioned.
(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 [WZ95,Tu16] 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 [WZ95,B04]. Using the transfer matrix formalism and the theory of Leonard pairs [T03,TV03], 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 mentioned. 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 The identity element is denoted I. We will use the standard q-shifted factorials (also called q-Pochhammer functions) [KS96], q-number and binomial, respectively: We also use the notation For the elementary symmetric polynomials in the variables {x i |i = 1, ..., n}, we use the notation: 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 [BTVZ18] follows from the transfer matrix associated with the solution of the reflection equation constructed in [Za95,B04]. 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 where u is the so-called spectral parameter. Then R(u) satisfies the quantum Yang-Baxter equation in the space Consider the reflection equation (also called the boundary quantum Yang-Baxter equation) introduced in the context of the boundary quantum inverse scattering theory (see [Ch84], [S88] for details). For the U q ( sl 2 ) R−matrix (2.1), the reflection equation reads: 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 [BK07]). Here, we consider the simplest solution of this type. We refer the reader to [Za95,B04] for details.
Note that the scalar parameter χ does not enter in the structure constants of the Askey-Wilson relations but it is introduced for further convenience.
Let us mention that this example can be derived using the 'dressing' procedure [S88] that generates solutions of the reflection equation. For details, see [Za95].
Example 2. [GZ93b] (see also [H16]) The Askey-Wilson algebra is embedded into The following gives an example of elements A, A * satisfying the Askey-Wilson relations (1.14), (1.15): where the structure constants are given by: 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 [DeG93,GZ94]: 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): 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)/(q − q −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: where the explicit form of the scalar function γ 0 (u) is omitted for simplicity and It follows that Γ satisfies 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 [BB17] 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: Then, it is easy to show that the couple (A, A ⋄ ) satisfies the Askey-Wilson relations: 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 [Cu07,H12] 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: Then, by direct calculation one finds that (Ā * ,Ā) satisfy 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 [CLSW03,D07] and its modified version [Be15,BP15,ABGP15] combined with the theory of Leonard pairs [T03,TV03]. 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.
where h is an integer such that 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 onV, it is not possible to construct a state that is annihilated by the off-diagonal entries. Indeed, given the representationV 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 [CLSW03] 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): such that the off-diagonal entries ofK(u|m) admit a reference state.
The gauge transformation is built as follows. For more details, we refer the reader to [CLSW03] and [BP15] for the notations used here. Let ǫ = ±1, α, β be generic complex parameters and m be an integer. Introduce the covariant vectors and the contravariant vectors Applying the gauge transformation to K(u), one finds: 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: 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: 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: and 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: 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: 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.
Lemma 3.2. Let m 0 be an integer. If the parameters α, β are such that: 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:π 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 For the choice of gauge parameter α fixed by (3.21), on one hand using (3.2) one finds: On the other hand, using (3.2) one finds: Combining both expressions, the off-diagonal contribution in |θ * 1 vanishes. It follows: 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: and define the parameter χ as Adapting the results of [T02a, T02b] the structure constants are given by: Then, the eigenvalues of the dynamical operators take the factorized form: To summarize, given a Leonard pair there are at least 9 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.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: 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: , (3.53) where, for further convenience, we have introduced: The spectral problem is now solved. Choose the gauge parameter α according to Lemma 3.1 or Lemma 3.2. For convenience, the parametrization (3.32)-(3.36) is used. Recall the notations (3.44), (3.45).
One has:π where the setū satisfies the Bethe equations: Proof. Recall the notation (1.21). Using the (off-shell) relations (3.46), (3.47) and the actions (3.31), one finds: with (3.37) and , 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 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: 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.
One has:π where the setū satisfies the Bethe equations: 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: 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: 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: According to the gauge transformation chosen, combining the expressions for A, A * it follows: 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: A crucial ingredient for the solution of the spectral problem of (3.66) is the following lemma 10 . 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: where we denote We now turn to the solution of the spectral problem.
where the setū satisfies the (inhomogeneous) Bethe equations: Proof. For convenience, define the element: Explicitly, using (3.69) for β = 0 it gives: where the notation 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: 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: 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: 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: 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π where Λ ǫ 1 (u) and Λ ǫ 2 (u) are given by (3.37) and 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: ≡ 0 determines the coefficient c M + in terms of the Leonard pair's data and the gauge parameters α, β: δ gπ (B ǫ (u, m 0 + 4s))|Ψ 2s g,ǫ (ū, m 0 ) = (3.92) 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 [ABGP15] 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).
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,ū) + 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: 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 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}.
Recall that the Bethe ansatz equations follow from requiring 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 fulfilled. 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 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 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 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 spaceV, as expected. We formulate the following conjecture: 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 }.

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: In terms of the elements A, A * , B, from the Askey-Wilson relations (1.14), (1.15) one obtains: 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): 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 }: 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: 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.
Replacing these expressions in (4.1), the proof of the following proposition is immediate.
According to the results of Section 3, the roots of the Q-polynomial are determined by (3.107) and the spectrum is given by 4.1.2. Special case κ = 0 and κ * = κ ± = 0. As mentioned 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: .
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: 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.

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 → C[z, z −1 ] such that: Note that π(A * ) is a specialization of the Askey-Wilson second-order q-difference operator, see e.g. [KS96] 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: 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: For convenience, we denote respectively by Q M (Z),Q 2s (Z) andQ 2s (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.
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). where P −1 ≡ 0 and From the three-term recurrence relation (4.22), one extracts the leading term of degree M in the variable x: 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: Expanding the Askey-Wilson polynomials in the variable Z = (z + z −1 )/(q + q −1 ) using (4.20) (see [KS96] for definitions), for the Q-polynomial (4.24) we obtain: Example 8. With the identification (4.25): 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 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: 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: 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 finds 12 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 [CGO16] 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 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 [Tu16] (see also [WZ95]). It is a quantum version of the Zhukovsky-Volterra gyrostat of classical mechanics [Ba09,LOZ06], and arises in spin systems with anisotropy [ZU87]. 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: 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 ).

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 [T03,TV03]. 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. [Ta17,KST18,Ta19].
• 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 [BK05,BK07], 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 [BB17]. 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 combination ( * * ) is a polynomial inÃ ⋄ of degree 2s + 1 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 =Ã ⋄ , one finds: We now consider ( * * ) on the finite dimensional representationV. 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 onV. 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.
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 2 i ) to get the final result.