A Relativistic Conical Function and its Whittaker Limits

In previous work we introduced and studied a function $R(a_{+},a_{-},{\bf c};v,\hat{v})$ that is a generalization of the hypergeometric function ${}_2F_1$ and the Askey-Wilson polynomials. When the coupling vector ${\bf c}\in{\mathbb C}^4$ is specialized to $(b,0,0,0)$, $b\in{\mathbb C}$, we obtain a function ${\mathcal R}(a_{+},a_{-},b;v,2\hat{v})$ that generalizes the conical function specialization of ${}_2F_1$ and the $q$-Gegenbauer polynomials. The function ${\mathcal R}$ is the joint eigenfunction of four analytic difference operators associated with the relativistic Calogero-Moser system of $A_1$ type, whereas the function $R$ corresponds to $BC_1$, and is the joint eigenfunction of four hyperbolic Askey-Wilson type difference operators. We show that the ${\mathcal R}$-function admits five novel integral representations that involve only four hyperbolic gamma functions and plane waves. Taking their nonrelativistic limit, we arrive at four representations of the conical function. We also show that a limit procedure leads to two commuting relativistic Toda Hamiltonians and two commuting dual Toda Hamiltonians, and that a similarity transform of the function ${\mathcal R}$ converges to a joint eigenfunction of the latter four difference operators.

1 Introduction This article may be viewed as a continuation of our previous work on a 'relativistic' generalization R of the Gauss hypergeometric function 2 F 1 , introduced in [1]. The latter paper and two later parts in a series [2,3] will be referred to as I, II and III in the sequel. The definition of the R-function in I is in terms of a contour integral, whose integrand involves eight hyperbolic gamma functions. (We review this in Section 2, cf. (2.1)-(2.5).) In recent years, van de Bult [4] tied in the R-function with the notion of modular double of the quantum group U q (sl(2, C)), as defined by Faddeev [5]. As a spin-off, he obtained a new representation of the R-function. Also, van de Bult, Rains and Stokman [6] have shown (among other things) that the 8-variable R-function R(a + , a − , c; v,v), a + , a − , v,v ∈ C, a + /a − / ∈ (−∞, 0], c ∈ C 4 , (1.1) can be obtained as a limit of Spiridonov's 9-variable hyperbolic hypergeometric function [7]. Their novel viewpoint leads to a third representation for the R-function. (See Proposition 4.20 and Theorem 4.21 in [6] for the latter two representations.) In this paper we are concerned with a 5-variable specialization of the R-function, defined by R(a + , a − , b; x, y) ≡ R(a + , a − , (b, 0, 0, 0); x, y/2). (1.2) Suitable discretizations of this function give rise to the q-Gegenbauer polynomials, whereas discretizations of the R-function yield the Askey-Wilson polynomials, cf. I; moreover, the nonrelativistic limit of the R-function yields the conical function specialization of 2 F 1 . Hence it may be viewed as corresponding to the Lie algebra A 1 , whereas the R-function can be tied in with BC 1 .
The key new result of this paper concerning R consists of the integral representation Here and throughout the paper we use parameters α ≡ 2π/a + a − , a ≡ (a + + a − )/2, (1.4) G(a + , a − ; z) is the hyperbolic gamma function (cf. Appendix A), and the dependence on a + , a − is suppressed. (We shall often do this when no confusion can arise.) Furthermore, in (1.3) we choose at first (1.5) By contrast to the previous three integral representations following from I, [4] and [6], the integrand in (1.3) involves only four hyperbolic gamma functions. We also obtain several closely related representations that involve in addition plane waves, cf. (3.47)-(3.51). As will transpire in Section 3, upon using the first one (1.3) of these novel representations (which we dub 'minimal' representations) to introduce the R-function, it is possible to rederive in a more transparent and self-contained way a great many features that also follow upon specialization of the Rfunction theory, developed not only in I, II and III, but also in our later papers [8] and [9]. Moreover, special cases and limits of the R-function are far more easily obtained from the minimal representations than from the original integral representation of I or from the alternative representations following from [4] and [6]. (The integrands of these earlier representations involve at least eight hyperbolic gamma function factors.) A survey of the results of I-III and [8] can be found in [10], but the definition (1.2) of the A 1 -analog of the (BC 1 ) R-function dates back to the more recent paper [9]. In Section 2, we review in particular the pertinent results from [9]. However, we have occasion to add a lot more information that follows by specializing previous findings concerning the R-function and related functions to their A 1 counterparts. This includes the asymptotic behavior and Hilbert space properties obtained in II and III, resp., which are adapted to the A 1 setting in Subsection 2.2, and the parameter shifts obtained in [8], which we focus on in Subsection 2.4. Moreover, in (2.27) we detail the connection of the renormalized function R r (a + , a − , b; x, y) ≡ G(ib − ia) G(2ib − ia) R(a + , a − , b; x, y), (1.6) to the A 1 type functions M (ma + + na − ; x, y), m, n ∈ Z, which featured in our previous papers [11] and [12]. We present the proof of (2.27) in Subsection 2.3, together with various corollaries. Altogether, Section 2 invokes a considerable amount of information from our previous work. We have attempted to sketch this in such a way that the reader need only consult the pertinent papers for quite technical aspects (in case of doubt and/or inclination, of course). Even so, it is probably advisable to skim through Section 2 at first reading, referring back to it when the need arises.
By contrast, Section 3 (combined with Appendices A and C) is largely self-contained. Its starting point is a hyperbolic functional identity that first arose as a specialization of elliptic functional identities expressing the relation of certain Hilbert-Schmidt integral kernels to the elliptic BC 1 relativistic Calogero-Moser difference operators introduced by van Diejen [13]. We need not invoke these identities (which can be found in [14], cf. also [15]), since the relevant hyperbolic version is quite easily proved directly. The key point is that the hyperbolic identities can be rewritten in terms of two pairs of hyperbolic A 1 -type relativistic Calogero-Moser difference operators A ± (b j ; x), j = 1, 2, with distinct couplings b 1 , b 2 . The difference operators are given by More specifically, the function B(b; x, y) is the Fourier transform of the hyperbolic kernel function, which is a product of four hyperbolic gamma functions. When we write the integrand of the integral defining B as a product of two factors that involve only two hyperbolic gamma functions, we can use the Plancherel relation and the explicit Fourier transform formula for factors of this type (derived in Appendix C) to obtain two new integral representations for B. In particular, this leads to a function C(b; x, y) given by However, we need a further study of the C-function (1.11) to arrive at a proof of this relation to the function R r , as defined originally by (1.6) and (1.2). Indeed, as already alluded to below (1.5), we can use (1.11) as a starting point to derive many features that C and R r have in common.
In particular, the general analysis in Appendix B of I can be applied to the integral on the r.h.s. of (1.11), which yields a complete elucidation of the behavior of C(b; x, y) under meromorphic continuation. Moreover, via the A∆Es (analytic difference equations) (1.10) and the manifest invariance of C under interchanging x and y, it follows that C(b; x, y) is a joint eigenfunction of the four A∆Os (analytic difference operators) with eigenvalues 2c + (y), 2c − (y), 2c + (x), 2c − (x), (1.14) resp. This is also the case for R r (b; x, y) and, moreover, the equality (1.12) can be shown for the special case y = ib by a further application of Appendix C. The general case then follows by a uniqueness argument already used in Subsection 2.3. We reconsider the special b-values b mn ≡ ma + + na − , m, n ∈ Z, (1.15) in Subsection 4.1, inasmuch as they satisfy b mn ∈ (0, 2a). Indeed, the new Fourier transform representations in Section 3 are only well defined for b ∈ (0, 2a), but they can be explicitly evaluated by a residue calculation when b is of this form. The key point is that the G-ratios in the integrand can then be written in terms of the hyperbolic cosines c ± (z) by using the G-A∆Es (A.2). In principle, this yields again the functions M (b mn ; x, y) from [11], but we have not tried to push through a direct equality proof (as opposed to appealing to uniqueness). Subsection 4.2 deals with the nonrelativistic limit. Specializing the results of I yields the hypergeometric function in terms of which the conical function can be expressed (cf. Chapter 14 in [16]). The five minimal representations (3.47)-(3.51) of the R-function lead to four representations (4.45)-(4.48) of the limit function. Rewriting them in terms of the conical function, three of these can be found in the literature (by looking rather hard). This is reassuring, since just as in I we were not able to get rigorous control on the nonrelativistic limits.
In order to describe the results of Section 5, we begin by recalling that in our paper [17] we arrived at relativistic nonperiodic Toda N -particle systems by taking a limit of the relativistic hyperbolic Calogero-Moser N -particle systems. In this limit the self-duality of the latter is not preserved, inasmuch as the dual commuting Hamiltonians have a very different character from the defining Hamiltonian and its commuting family. Specialized to the present context, this limit can be used to obtain a joint eigenfunction of two Toda Hamiltonians H T ± (η; x) and two dual Toda HamiltoniansĤ T ± (η; y), with the real parameter η playing the role of a coupling constant.
The limit transition proceeds in two stages. The first step is to set b = a − iγ, γ ∈ R.
(1. 16) At the classical level the analogous b-choice still yields real-valued Hamiltonians with a welldefined self-dual action-angle map and scattering theory [18]. Correspondingly, the four reduced N = 2 quantum Hamiltonians at issue here are still formally self-adjoint for this b-choice. (They are similarity transforms of the A∆Os (1.13) with a weight function factor.) Moreover, restricting attention to their joint eigenfunction remains real-valued, although this reality property is no longer manifest: It hinges on a symmetry property under taking b to 2a − b, which translates into evenness in the parameter γ. The next step is to substitute 18) and take Λ to ∞. In this limit the Hamiltonians and their joint eigenfunction converge, whereas the dual Hamiltonians must be multiplied by a factor e δ (−Λ) to obtain a finite limit. This can be understood from their Λ-dependent eigenvalues 2c δ (x + Λ) following from the x-shift (1.18), cf. (1.14). Indeed, after multiplication by this renormalizing factor the eigenvalues have the finite limits e δ (x), δ = +, −. The five representations of the R-function give rise to four representations of the relativistic Toda eigenfunction F T (η; x, y), namely (5.25), (5.26), (5.32) and (5.33). Suitably paired off, however, these different formulas express real-valuedness with (1.17) in effect. Taking this into account, we wind up with two essentially different representations that are intertwined via the Plancherel formula for the Fourier transform. The key formula involved here is derived in Corollary C.2.
The results just delineated can be found in Subsection 5.1. In Subsection 5.2 we first study the asymptotic behavior of F T (η; x, y) for x → ±∞ and y → ∞. We then clarify the analyticity properties of F T (η; x, y) by introducing a similarity transform H(x − η, y). Using the four representations (5.54)-(5.57) of the latter, we show that the function H(x, y) is holomorphic for (x, y) ∈ C 2 . Subsection 5.3 deals with the joint eigenfunction properties of F T (η; x, y) and its similarity transforms. Formally, these follow from those of the R-function. However, the Toda limit is not easy to control analytically, and the direct derivation of the eigenvalue equations is not too hard and quite illuminating.
Our results in Section 5 have some overlap with earlier results by Kharchev, Lebedev and Semenov-Tian-Shansky [19], who obtained functions closely related to F T (η; x, y) from the viewpoint of harmonic analysis for Faddeev's modular double of a quantum group [5]. The nonrelativistic nonperiodic Toda eigenfunctions are widely known as Whittaker functions, and meanwhile it has become customary to call eigenfunctions for q-Toda Hamiltonians Whittaker functions as well. In particular, q-Whittaker functions were introduced by Olshanetsky and Rogov for rank 1 (their work can be traced from [20]) and by Etingof for arbitrary rank [21], and these functions have been further studied in various later papers (see e.g. [22] and references given there). We would like to stress that these functions are quite different from the ones at issue here and in [19]. The crux is that the former are only well defined for q not on the unit circle, whereas here and in [19] the eigenfunctions have a symmetric dependence on two generically distinct q's, given by This state of affairs is closely related to the different character of the trigonometric gamma function (more widely known as the q-gamma function, with the restriction |q| = 1 being indispensable) and the hyperbolic gamma function (which depends on parameters a + and a − in the right half plane). In Section 6 we study the nonrelativistic limit of the representations of the relativistic eigenfunction, arriving at two distinct representations for the nonperiodic Toda eigenfunction that have been known for a long time. Just as for the relativistic case, its property of being also an eigenfunction for a dual Hamiltonian seems not to have been observed before. (These duality features are the quantum counterparts of duality features of the pertinent action-angle maps, first pointed out in [17].) To control one of the two pertinent limits, a novel limit transition for the hyperbolic gamma function is needed, whose proof is relegated to Appendix B.
2 The R-function as a special case of the R-function

The functions R and R r
The R-function (1.1) is defined as a contour integral over a variable z, with the z-dependence of the integrand encoded in a product of eight hyperbolic gamma function factors. (See Appendix A for a review of the relevant features of the hyperbolic gamma function.) Specifically, with suitable restrictions on the eight variables, the R-function is given by Here we havê , and K is given by with new parameters Also, recall a and α are defined by (1.4).
We do not need the definition of the contour C for general variable choices (this is discussed in I and Section 4 of the survey [10]); instead we presently define C for the cases at issue. For the special c-choice in (1.2) we can use the duplication formula (A.10) to obtain Using also the reflection equation (A.6) we deduce that R is given by For the variable choice that is most relevant for Hilbert space purposes, namely, a + , a − , b, x, y > 0, (2.9) the contour C may be chosen equal to the real line in the z-plane, indented downwards near 0 so as to avoid a pole of K(b; z). From (A.17)-(A.18) it follows that the poles of K(b; z) are located on the imaginary axis at (2.10) Thus they are above the contour, whereas the remaining z-poles of the integrand at are below C. From the above representation it is immediate that R is symmetric under the interchange of the parameters a + and a − : R(a + , a − , b; x, y) = R(a − , a + , b; x, y). (2.12) It is not at all clear, however, that R is also symmetric under the interchange of the positions x and y: This self-duality feature follows in particular from a second relation between R and R, namely, (This is equation (4.8) in [9].) Indeed, this second c-choice yields the same function K(b; z) as the first one, so that substitution of (2.1) (with the same contour C) now yields (2.8) with x and y interchanged on the r.h.s. There are two more c-choices that lead from R to R, namely, (b, 0, b, 0) and (b, b, 0, 0). Specifically, from equations (4.6) and (4.7) in [9] we have From the definition of the R-function we then obtain alternative integral representations for the R-function from which the self-duality property (2.13) is manifest. (Indeed, since we have c 0 =ĉ 0 = b for these two choices, the integrand is invariant under the interchange of x and y.) On the other hand, the modular invariance property (2.12) is not at all clear, since the integral representations involve the hyperbolic gamma function with a − replaced by 2a − . Using (A.11), they can be re-expressed in terms of the modular invariant function G(a + , a − ; z). However, the resulting integrand is then still not modular invariant. Since it seems not to simplify and does not look illuminating, we do not detail it any further. The analyticity properties of the R-function are known in great detail from Theorem 2.2 in I, cf. also Section 4 in the survey [10]. Combining this theorem with the definition (1.2) of R and its self-duality property (2.13), we deduce in particular that R extends from the intervals (2.9) to a meromorphic function in b, x and y, whose poles in x and y can only occur at the locations We proceed to list further consequences of the R-function theory for R. Two features that are clear from each of the above integral representations are evenness and scale invariance (given scale invariance of G): the eigenvalue A∆E (analytic difference equation) for A + (b; y) entails In view of (2.20), it follows from this that R n is of the form R n (x) = P n (c + (x)), (2.24) where P n (z) is a polynomial in z of degree n and parity (−) n . The relation of these polynomials to the q-Gegenbauer polynomials and to the Askey-Wilson polynomials associated with the four relevant c-choices is detailed at the end of Section 4 of [9]. The renormalized R-function R r given by (1.6) is the counterpart of the renormalized Rfunction R r obtained from (2.1) by omitting the product j G(is j ) in K, cf. (2.4 The renormalizing factor in the function R r ensures that it has no poles that are independent of x and y, cf. Theorem 2.2 in I. More precisely, R r (a + , a − , b; x, y) extends to a function that is meromorphic in the domain and whose poles can only occur at the locations (2.17).
It is not obvious, but true that for the special b-choices b mn (1.15) we have an equality where M (b mn ; x, y) is the function defined at the end of Section III in our paper [11]. Therefore, R r (b; x, y) is the continuous (indeed, real-analytic) interpolation to arbitrary b ∈ R of the function given by equation (3.74) in [11], which is defined only for the b-values b mn . (Note the latter are dense in R when the ratio a + /a − is irrational.) It will not cause surprise that in the free case we have M (b 00 ; x, y) = exp(iαxy/2). (2.28) It would take us too far afield, however, to detail all of the functions M (b mn ; x, y), m, n ∈ Z. For our purposes it is enough to specify their general structure: They are 'elementary' in the sense that they can be written M (b mn ; x, y) = exp(iαxy/2)R mn (e + (x), e − (x), e + (y), e − (y)), (2.29) where R mn is a rational function of its four arguments, cf. Section III in [11]. In Subsection 2.4 we deduce this structure in another way (namely, by exploiting parameter shifts). Moreover, for the case where m and n are not both positive or both non-positive, this structure can be understood from the novel Fourier transform representations (3.48)-(3.51), cf. Subsection 4.1.
We postpone the proof of the equality assertion (2.27) to Subsection 2.
3. An ingredient of this proof is the asymptotic behavior of R r (b; x, y) as x goes to ∞, and this is most easily obtained as a corollary of the asymptotics of a closely related function E(b; x, y), defined by (2.38).

The functions E and F
The function E(b; x, y) can be viewed as a specialization of the function denoted E(γ; v,v) in II and [10]. The relation between γ and c reads γ(a + , a − , c) = (c 0 − a, c 1 − a − /2, c 2 − a + /2, c 3 ). (2.30) In particular, the 'free' case c = 0 yields The switch from c to γ is crucial for uncovering further symmetries: The function E(γ; v,v) is invariant under D 4 transformations on γ (i.e., permutations and even sign changes), cf. II. (In [6] this D 4 symmetry has been reobtained in a quite different way.) It is defined by Here, the generalized (BC 1 ) Harish-Chandra c-function is given by 34) and the constant by Denoting the γ-vectors corresponding to the two one-parameter families by γ (1) , γ (2) , it is easy to check that a straightforward calculation (using the duplication formula (A.10)) yields where we have introduced the constant and generalized (A 1 ) Harish-Chandra c-function Recalling (2.14) and using (2.37), it readily follows that we also have E(b; x, y) = E γ (2) ; x/2, y . (2.42) It involves more work to obtain the relations between E and E corresponding to (2.15) and (2.16). Setting (2.43) these are given by (2.44) (These formulas amount to special cases of the doubling identity for the E-function obtained in Section 6 of [9].) The relation (2.39) between E and R r yields a similarity transformation turning the A∆Os From this it is easy to verify that these A∆Os are formally self-adjoint operators on L 2 (R) (by contrast to the A∆Os (1.13)), and that they are invariant under the transformation It is not obvious, but true that we also have This symmetry property can be derived from (2.38) and the D 4 invariance of the E-function: We have corresponding to the factor G(x + ia)G(y + ia), and b-dependent poles at The main disadvantage of the function E(b; x, y) compared to the R r -function is that it is not even in x and y, since the c-functions in (2.39) are not even. Instead, it satisfies On the other hand, E inherits all other important properties of R r , and is the simplest function to use for Hilbert space purposes. In particular, it has the 'unitary asymptotics' cf. Theorem 1.2 in II. Here, the u-function encodes the scattering associated with the A∆Os A ± (b; x), reinterpreted as commuting self-adjoint operators on the Hilbert space L 2 ((0, ∞), dx). More precisely, using corresponding results on the E-function from II and III, it follows that the generalized Fourier transform The self-adjointness of the operatorsÂ ± (b) on H associated to the A∆Os A ± (b; x) for b in this interval can then be easily understood from the unitarity of F: they are the pullbacks to H under F of the self-adjoint operators of multiplication by 2c ± (y) onĤ.
As already mentioned, these statements follow from II and III by specialization, but it may help to look first at Section 9 in the survey [10]. Starting from the representation (2.38), the vector γ (1) belongs to the polytope P given by equation (9.2) in [10], provided b ∈ (0, 2a). Therefore, the transform associated with E(γ (1) ; v,v) (defined by equation (9.4)) is an isometry. As a consequence, the transform (2.56) is an isometry. (The normalization factor in (2.56) differs from that in equation (9.4) in [10] to accommodate the scale factor 1/2 in the y-dependence of E in (2.38).) Isometry of the inverse transform is then clear from the self-duality of E(b; x, y). Next, for b = 0 we have the identity E(0; x, y) = R r (0; x, y) = 2 cos(αxy/2), (2.58) cf. (2.39)-(2.41) and equation (7.33) in [10] (also, note that for b = 0 we have γ (1) = γ f , cf. (2.49) and (2.31)). In view of the symmetry (2.48), it follows that F amounts to the cosine transform for b = 0 and b = 2a, so these transforms are unitary as well. More generally, we obtain a family of unitary operators which is strongly continuous on the parameter set Π u and satisfies with φ(b) given by (2.40). Hence the u-function (2.53) has asymptotics Also, the reflection equation (A.6) and the complex conjugation relation (A.9) entail Thus, if we set (with the square root phase factors reducing to 1 for b = 0), then F has asymptotics and if we replace E by F in the above unitary transform (2.56), we retain unitarity. Introducing the weight function we have and we can also write F in terms of R r as F(b; x, y) = w(b; x) 1/2 w(b; y) 1/2 R r (b; x, y), b, x, y > 0, (2.68) with the positive square roots understood. Hence, F is a joint eigenfunction of the A∆Os H ± (b; x) and H ± (b; y) with eigenvalues (1.14), where For b ∈ (0, 2a) the reduced weight function w r (b; x) is positive for all real x, and since it is also even, its positive square root for x > 0 has a real-analytic extension to an even positive function on all of R. By contrast, it is clear from (2.70) that w(b; x) 1/2 , x > 0, extends to an odd real-analytic function on R. As a consequence, one can also view the transform associated with F(b; x, y), b ∈ (0, 2a), as a unitary transform from the odd subspace of L 2 (R, dy) onto the odd subspace of L 2 (R, dx). This is the viewpoint taken in [12], where we studied this transform (among other ones) for the special b-values N a + with N ∈ N * . As shown there, for b > 2a unitarity and self-adjointness generically break down in a way that can be understood in great detail.
To be sure, the precise connection between the above functions F and R r and the functions F r and E r from [12] is not clear at face value. But the latter are derived from the functions M ((N + 1)a + ; x, y) of [11], as specified below equation (1.42) in [12], so this connection is encoded in the identities (2.27) for the special cases (m, n) = (N + 1, 0), N ∈ N.

The identities (2.27) and their consequences
We proceed to prove the general identities (2.27). Our reasoning involves in particular a comparison of the behavior for x → ∞ of the functions on the l.h.s. and r.h.s. For R r (b; x, y) this asymptotics easily follows upon combining (2.39), (2.54) and (2.61): Next, we consider the functions M (b mn ; ±x, y). To begin with, they are eigenfunctions of the four A∆Os (1.13) (where b = b mn ) with eigenvalues (1.14), cf. Theorem II.3 in [11]. Their 'elementary' form (2.29) follows from equations (3.65)-(3.68) in [11]. The function K N + ,N − (a + , a − ; x, y) occurring in these formulas is specified in equation (3.2), with S N δ given by equation (2.21). In turn, the coefficients in equation (2.21) are defined via equations (2.2)-(2.5) in [11]. (See also Subsection 4.1 for more information on these special cases.) It is straightforward to obtain the asymptotics for Re (x) → ∞ from these explicit formulas. Specifically, this yields The decay rate ρ is the minimum of the two numbers 2π/a ± , and the implied constant can be chosen uniform for Im (x) varying over R. Comparing (2.71) and (2.72), it follows that the functions on the l.h.s. and r.h.s. of (2.27) have the same asymptotics for x → ∞. It therefore suffices to prove that for fixed a + , a − , y > 0 and m, n ∈ Z, they must be proportional as functions of x. Moreover, we may as well assume a + /a − is irrational, since equality for this case entails equality for all a + , a − > 0. (Indeed, the functions M (b mn ; ±x, y) are manifestly real-analytic in a + and a − for a + , a − > 0, and this real-analyticity property is also valid for R r (b; x, y), cf. I.) The key consequence of the irrationality assumption is that the vector space of meromorphic joint solutions f (x) to the A∆Es is two-dimensional. To explain why this is so, we first note that the functions M (b mn ; ±x, y) are independent solutions to (2.73), their independence already being clear from their general form (2.29). Moreover, it follows from their uniform asymptotics (2.72) that there exists a positive number Λ, depending on the fixed variables a + , a − , m, n and y, but not on x, such that in the half plane Re (x) > Λ both functions are zero-free, and satisfy We are now in the position to invoke a result from Section 1 in [23], to the effect that the above suffices for any joint meromorphic solution f (x) of (2.73) to be a linear combination of the two functions M (b mn ; ±x, y). Since R r (b mn ; x, y) is an even meromorphic joint solution, the functions on the l.h.s. and r.h.s. of (2.27) are proportional, so their equality now follows. In particular, for the free case m = n = 0 we recover the identity (2.58) from (2.27)-(2.28). Moreover, taking y = ib mn in (2.27), we can invoke (2.25) to deduce the corollary Using the G-A∆Es (A.2), the r.h.s. can be rewritten in terms of sine-functions. For the special case m = N + 1, n = 0, the resulting identity amounts to equation (2.78) in [11], cf. also Subsection 4.1.
We would like to add in passing that it is very plausible that (2.74) is not necessary for twodimensionality. Indeed, denoting by P c the field of meromorphic functions with period c ∈ C * , any third independent joint meromorphic solution would have to be both of the form and of the form Since the intersection of the fields P ia + and P ia − reduces to the constants when a + /a − is irrational, we expect (but are unable to prove) that this simultaneous representation should lead to a contradiction without appealing to (2.74). Now that we have proved (2.27), it follows that the function R r (a + , a − , b; x, y), which is real-analytic on the parameter set is the continuous interpolation of the functions on the r.h.s. of (2.27), which are only defined for the dense subset of 'elementary' parameters The natural question whether another linear combination of M (b mn ; ±x, y) that is independent from the even one admits a continuous interpolation as well remains open. In this connection we should point out that our reasoning at the end of Section 3 of [24] renders this extremely unlikely, but is not conclusive. Indeed, we cannot rule out that the sequence of functions Q − given by equation (3.15) in [24], with N + ∈ N, gives rise to an infinity of distinct limits L − , corresponding to distinct subsequences. (This oversight is of no consequence for the later sections in [24].) For the same reason, the analogous assertion about the R-function, made at the end of [10], has not been completely proved. Before turning to parameter shifts, we derive a non-obvious reality feature of R r from the relations (2.27), namely, The point is that from the explicit formulas for the functions M it is apparent that we have

Parameter shifts
We continue to summarize results concerning parameter shifts from [8], inasmuch as they apply to the present A 1 context. In Section 1 of [8] we introduced the up-shifts (2.83) and the down-shifts where δ, δ ′ = +, −. Clearly, the up-shifts S − (x) commute, and the down-shifts S (d) x) commute as well. The shifts are also related by where δ = +, −. It is a matter of straightforward calculations to verify the formulas (2.83) and (2.85)-(2.87). Starting from the joint eigenfunctions exp(±iαxy/2) of A ± (0; x) with eigenvalues 2c ± (y), one can now obtain joint eigenfunctions with the same eigenvalues for A ± (b mn ; x) by acting with the shifts on the plane waves. By construction, these joint eigenfunctions are of the elementary form (2.29). Choosing a + /a − irrational, it follows from two-dimensionality of the joint eigenspace that these eigenfunctions are (generally y-dependent) multiples of M (b mn ; ±x; y).
A more telling action of the shifts is encoded in Indeed, these relations hold for arbitrary b. For b = b mn , it then follows by using (2.27) that they also hold for the summands M (b mn ; ±x, y). (This is because their plane wave factors are independent, cf. (2.29).) The equations (2.89) and (2.88) follow from a suitable specialization of equations (3.11) and (3.13) in [8]. But in the present A 1 case we can also derive them quite easily by using the elementary joint eigenfunctions M (b mn ; x, y) with a + /a − / ∈ Q. Indeed, once we have shown that (2.89), (2.88) hold for y > 0, b = b mn , m, n ∈ Z, and with R r replaced by M , it is easy to deduce (2.89), (2.88) from (2.27) and interpolation. (Note that the four shifts commute with parity.) Their validity for these special cases can be readily verified: One need only show that the functions on the l.h.s. and r.h.s. have the same x → ∞ asymptotics, and using (2.72) this causes little difficulty.
Next, we obtain the counterparts of the A δ (b; x)-and R r -shifts for the A∆Os A δ (b; x) and their joint eigenfunction E(b; x, y). (For the BC 1 setting we did this in Section 8 of [9]; as in previous cases, it is in fact simpler and more illuminating to obtain the relevant formulas by direct means, instead of by specialization.) They are given by A moment's thought shows that this entails the validity of (2.85)-(2.87) with S, A replaced by S, A. Also, using the definition (2.41) of the c-function and the G-A∆Es (A.2), we obtain the explicit formulas Notice that they imply Finally, a straightforward calculation yields the counterparts of (2.88) and (2.89): To conclude this section, we point out that the eight shifts acting on x have duals acting on y given by the formulas (2.82), (2.84) and (2.90)-(2.93) with x → y. By self-duality, their respective actions on R r (b; x, y) and E(b; x, y) follow from the above by interchanging x and y.

Five minimal representations of the R-function
We begin this section by focusing on the kernel function We have established that this function satisfies three independent kernel identities. We expect that these might be useful in other contexts than the present one. Indeed, here we only need the special case (3.14) of the first of the identities. We collect the three identities in the following proposition.
Proof . To prove (3.2), we divide l.h.s. and r.h.s. by K(b; x − ia −δ , v) and use the A∆Es (A.2) to write the result as Both sides are 2ia δ -periodic functions of x with equal limits The residues at the (generically simple) poles x = 0, x = ia δ in the period strip clearly cancel. By Liouville's theorem, it remains to check that the residues at the poles x = ±v − ib/2 ± ia δ cancel as well, and this is a routine calculation. Next, we divide (3.3) by K(b; x, v) and use (A.2) to obtain Both sides are 2ia δ -periodic functions of x with equal limits The residues at x = 0 and x = ia δ manifestly cancel. It is a straightforward calculation to verify that the residues at the remaining poles x = ±ia δ ± ia −δ /2 ± v + ib cancel, too. Hence (3.3) follows. Finally, to prove (3.4), we divide both sides by K(b; x − ia −δ , v) and use (A.2) to get as the counterpart of (3.5): Both sides are 2ia δ -periodic functions of x with equal limits As before, residue cancellation at x = 0 and x = ia δ is immediate, whereas the verification that the residues at the remaining poles x = ±ia δ ± v − ib/2 cancel as well involves a bit more work.
From (1.7) we see that the identity (3.2) can be rewritten as At first sight, one might think that the identities (3.3) and (3.4) can also be rewritten by using a rescaled version of the two commuting A 1 difference operators A ± (b; z). The two difference operators corresponding to (3.3) do commute. Even so, one can only rescale one of the operators such that it takes the A 1 form (1.7), but not both at once. For the purpose of studying the A 1 operators, then, we can only make use of (3.2). More specifically, our starting point is the special case d = 0: In order to exploit this identity, we introduce the Fourier transform The integral is well defined, since the b-restriction ensures that the v-poles of the integrand at These features also imply that B extends from the real x-axis to a function that is holomorphic in the strip Im Next, we temporarily assume b ∈ (0, a s /2), a s ≡ min(a + , a − ). (3.18) Then the action of the shifts in the A∆Os A ± (b; x) given by (1.7) is well defined on B(b; x, y), provided we restrict x to a strip |Im x| < a s /2. Moreover, we may take the shifts under the integral sign in (3.15) and use the kernel identity (3.14) to obtain |Im x| < a s /2.
Upon shifting contours R → R ± ia −δ , no poles are met, and so we obtain The integrands of both terms are now equal, and the contours can be shifted back to R without changing the value of the integrals. Hence we deduce the eigenvalue equations Reverting to our previous assumption b ∈ (0, 2a), we proceed to obtain two different representations of B(b; x, y). To this end we use the Plancherel relation with the Fourier transform defined bŷ we now define We can calculate the Fourier transforms of these four functions by using Proposition C.1. Doing so, we use the Plancherel relation (3.22) and then replace q by z + y/2 to obtain the two representations announced above: Next, we compare (3.28) to (3.24), deducing that it can be rewritten as Also, defining a new function C(b; x, y) by (1.11), we see that (3.27) amounts to The function C(a + , a − , b; x, y) is of central importance for what follows. Writing it as where we have introduced we infer that its behavior under analytic continuation in its 5 variables is immediate from the general analysis in Appendix B of I (with N specialized to 2). We proceed to summarize the salient information.
To this end, we need the function the product function P (a + , a − , b; x, y) extends from (0, ∞) 2 × (0, a + + a − ) × R 2 to a function that is holomorphic in the domain Hence C is meromorphic in D(a + , a − , b), with poles occurring solely at the zeros of the E-product, cf. (A.21); moreover, the maximal multiplicity of a pole at z = z 0 , with z = x, y, is given by the zero multiplicity at z = z 0 of the pertinent E-factor. The corresponding meromorphy properties of B are now clear from its relation (3.30) to C: It continues meromorphically to the domain D(a + , a − , 2a − b). From (3.29) we then deduce that B(a + , a − , b; x, y) has a meromorphic extension to the larger domain D + (2.26). Using (3.30) again, it now follows that C has a meromorphic extension to D + as well.
We proceed to obtain further information on the function C. First, we list features that are immediate from its definition (1.11) and properties of the G-function, cf. Appendix A: x, y ∈ R, (real-valuedness). (3.40) Clearly, the relations (3.36)-(3.39) are well defined and hold true on D + (2.26). Also, the property (3.40) can be rendered manifest by substituting the integral representation (A.5) in the four G-functions and combining factors to obtain Second, combining (3.37) with (3.29) and (3.30), we deduce Third, we can use Proposition C.1 once more to obtain from (1.11) the explicit result Last but not least, we claim that we have the joint eigenvalue equations To prove this claim, we first note that the eigenvalue equations (3.21) continue meromorphically to D + . Next, we observe that the G-A∆Es (A.2) imply which is equivalent to the first two A∆Es in (3.44). The last two are then clear from the self-duality relation (3.37). All of the properties of C just derived also hold true for the function G(ib − ia)R r defined by (1.6) and (1.2), cf. Section 2. By using solely the eigenvalue properties (3.44), the evenness properties (3.38), and the normalization (3.43), we can now show that these two functions coincide, as announced in the Introduction, cf. (1.12). Specifically, applying the uniqueness argument in Subsection 2.3 to C in its dependence on x, we obtain the equality (1.12) up to a proportionality factor p(a + , a − , b, y). Repeating this argument for the y-dependence, we see that the proportionality factor can only depend on the parameters a + , a − , b. From the normalization relation (3.43) it then follows that p = 1, thus proving (1.12).
Let us now collect the resulting minimal representations of the R-function. From (1.12) and (1.6) we obtain Next, combining (3.24) and (3.30), we deduce and using (3.29) we infer Finally, using the self-duality property of R, we obtain from (3.48) and (3.49) the representations Taking stock of the above developments, we note that we might have started from the first minimal representation (3.47) to define the R-function. Then many of its properties follow quite easily. On the other hand, it seems not feasible to give a direct proof of its crucial joint eigenfunction property. With hindsight, however, this can be shown by first obtaining the second representation (3.48) (say) via Proposition C.1, and then using the identity (3.14) to arrive at (3.21). From this the joint eigenfunction property (3.44) follows as before.
Another important property of R is its asymptotic behavior for x → ∞. Like other features addressed in this section, this can already be gleaned from Section 2, via the specialization of the more general asymptotics of the function E(γ; v,v) obtained in Theorem 1.2 of II. However, provided we restrict b to the interval (0, 2a), it is quite easy to obtain the x → ∞ asymptotics directly from the new representations of R in terms of a Fourier transform.
To detail this, let us first note that we need only consider the function E(b; x, y), which we can now view as being defined via (2.39)-(2.41). (Indeed, there is no difficulty in obtaining the asymptotics of the c-function; in this connection, compare (2.61), (2.54) and (2.71).) Using (3.49), we deduce that for b ∈ (0, 2a) the E-function has the representation Letting y ∈ (0, ∞), we can shift the contour up by a − b/2 + ǫ, where ǫ > 0 is small enough so that only the simple poles at z = ±y/2 − ib/2 + ia, with the multiplier given by Now from (A.13) we see that M (b; x) converges to 1 for x → ∞. To recover the asymptotics (2.54), therefore, it is enough to show that the r.h.s. of (3.52) with z replaced by z + ia − ib/2 + iǫ vanishes for x → ∞.
To prove this, we write the shifted contour integral as Now one need only use (A.13) to verify that the integrand is bounded by a multiple of exp(−αb|z|/2), which implies the function (3.56) does converge to 0 for x → ∞.
We stress that this short argument only yields (2.54) under the restriction b ∈ (0, 2a). In particular, by contrast to the previous contour integral representation used in II, one must cope with an inevitable contour pinching when one tries to use (3.52) to go beyond this b-interval.
Another issue is that stronger asymptotic estimates than just obtained are necessary to recover the Hilbert space transform features for the E-function sketched in Subsection 2.3, cf. (2.55)-(2.60). It is beyond our scope to study this further, but we would like to repeat that the b-interval [0, 2a] cannot be enlarged without losing the critical unitarity and self-adjointness properties [12].
At face value, the new representations (3.47)-(3.51) seem to hold promise for a direct proof of the shift properties of R r , cf. (2.88)-(2.89). Even so, we were unable to push this through. To date, therefore, the only reasoning yielding the properties for general b is to first derive them for special b-values, as sketched in Subsection 2.4. We now turn to a study of the R r -function for these special values.

Specializations and nonrelativistic limit 4.1 Elementary special cases
As already mentioned, the functions R N (a + , a − ; x, y) ≡ R r (a + , a − , (N + 1)a + ; x, y), N ∈ N, (4.1) have been extensively studied before. They were first obtained more than twenty years ago [25], and then reconsidered from an algebraic and function-theoretic viewpoint in [11] and from a representation-theoretic viewpoint in a paper by van Diejen and Kirillov [26]. The corresponding Hilbert space transforms were studied in great detail in [12]. Our first goal in this section is to demonstrate how the elementary character of these functions can be directly understood from the Fourier transform representations (3.48)-(3.51). Indeed, thus far the relation encoded in (2.27) has only been shown by appealing to a uniqueness argument. The crux is that for the choices b = (N + 1)a ± one can use the G-A∆Es (A.2) to obtain integrals that can be explicitly evaluated by a residue calculation.
Specifically, let us start from (3.49) to obtain first Now we recall that (3.49) is valid for b ∈ (0, 2a), which implies we have N a + < a − in (4.2). Taking y > 0 from now on, it follows that the integrand has 2N + 2 simple poles in the strip Im z ∈ (0, a − ). The product yields a function that is ia − -periodic in z. Thus, denoting the integral by I N , we have The residues are easily calculated, and hence we obtain Also, the prefactor can be calculated by using (A.2) once more, combined with (A.12). Introducing we get For N = 0 the product of (4.4) and (4.6) yields More generally, the product can be written as where we have set Introducing the phase factor it is not hard to see from (4.9) that K N is of the form Thus far, our conclusions about R N and K N were only based on an explicit evaluation of the representation (3.49). However, a lot more information follows upon using the features of R r (b; x, y). In particular, the function K N (x, y)/P N (x)P N (y) is a joint eigenfunction of the four A∆Os (1.13) (where b = (N + 1)a + ) with eigenvalues (1.14), since this holds true for R N (x, y). Also, the self-duality and evenness of R N imply 13) and this entails that the coefficients in (4.12) have the symmetry properties (4.14) Of course, this can be easily checked for small N , but for arbitrary integers (4.14) is not at all obvious from (4.9). One more feature of the coefficients is that they are Laurent polynomials in q with integer coefficients. Like the symmetries (4.14), it is not a routine matter to show this from (4.9). The crux is, however, that the above functions K N coincide with those of [11], by virtue of the uniqueness argument used in Subsection 2.3. The coefficients were studied in detail in Section II of [11], and there the interested reader can find explicit formulas for the coefficients as Laurent polynomials in q. See also the paper by van Diejen and Kirillov [26], where yet different representations of the functions K N (x, y) were obtained.
With a little more effort, the elementary character of R r for the more general can also be understood from (3.49). Indeed, from (A.2) it follows by a straightforward calculation that we have an identity (4.16) Using this identity several times (together with (A.12) for the factor G(ib +− − ia)), we deduce from (3.49) and (1.6) the representation The denominator of the integrand has no zero for z ∈ R unless M is odd and N is even. The zeros of the corresponding factor s − (z ± y/2) are then matched by the zeros of the factor s + (z ± y/2) of the numerator. After a suitable contour shift, we can expand the numerator product into exponentials, yielding a sum of convergent integrals (recall we require b +− ∈ (0, 2a)). When M is even or N is odd, we can do the same without a contour shift. Each of the integrals is then basically of the same form as previously evaluated for the case M = 0. (More in detail, the same 2N + 2 poles arise in the period strip for the c − -product.) From these observations the general structure anticipated in (2.27) and (2.29) readily follows, provided n > 0 and m ≤ 0, or m > 0 and n ≤ 0. In Section III of [11] we studied the functions (4.17) in considerable detail, but it is beyond our scope to derive the explicit form used there directly from their representation (4.17). We do add that it seems plausible that the factorization exhibited in equations (3.3)-(3.4) of [11] can be understood by a more refined analysis of the above sum of contour integrals. In any case, we repeat that equality of the pertinent functions follows from the uniqueness argument explained in Subsection 2.3.

The nonrelativistic limit
We begin this subsection with a remark addressed to physicist readers, who may care about dimension issues. In our paper [11], which we had occasion to cite several times in the previous subsection, the variable y of the present paper was denoted by p. This is a widely used notation for the spectral variable in nonrelativistic quantum mechanics, where p is viewed as a momentum. In our relativistic setting, however, it is far more natural to view the scale parameters a + and a − as having dimension [position], and then the 'geometric' and 'spectral' variables x and y both have dimension [position] as well. (To be more specific, one of the scale parameters can be viewed as an interaction range, and the other one as the Compton wave length /mc of the relativistic particles under consideration.) This goes to explain our change from p to y.
Of course, from a mathematical viewpoint such notational issues and the notion of dimension may be ignored. When taking the nonrelativistic limit, however, these physical considerations point the way. We need to let the speed of light c go to ∞, so one of the scale parameters should go to 0. In particular, we cannot retain modular invariance. Accordingly, we first set a + = 2π/µ, a − = β, µ, , β > 0. Here, we view β as 1/mc, with m the particle mass, and we trade a + for a parameter µ with dimension [position] −1 to avoid a great many factors π. Next, the spectral variable y (dual position) is replaced by the momentum variable p = µy/β. With these changes in effect, it is easy to verify the expansion where so that the eigenvalue of A becomes p 2 /4, while the eigenvalues of the monodromy operator M and dual A∆OÂ remain 2 cosh(πp/ µ) and 2 cosh(µx/2), resp. Likewise, the similarity-transformed A∆Os A ± (2.45) and Hamiltonians H ± (2.69) yield where the Harish-Chandra c-function is given bŷ , (4.34) and the operatorsÂ andĤ by similarity withŵ nr (g/ ; p/ µ) 1/2 , wherê w nr (λ; k) ≡ 1/ĉ nr (λ; ±k). Next, we study the nonrelativistic limit of the R-function. To this end we first use the definition (1.1), the scaling property (2.19) and the limit lim t→0 R(π, t, tc; v, tu) established and discussed in I. We write the limit lim β→0 R(2π/µ, β, gβ; x, βp/µ) =: ψ nr (g/ ; µx/2, p/ µ), (4.39) in terms of the dimensionless quantities which is a well-known quadratic transformation, cf. e.g. [28, p. 125].
Using other familiar features of the hypergeometric function, it is not difficult to verify that the operators A, M andÂ do have the expected eigenvalues p 2 /4, 2 cosh(p/ µ) and 2 cosh(µx/2) on the limit function ψ nr (g/ ; µx/2, p/ µ). More specifically, for A this amounts to the ODE satisfied by 2 F 1 and forÂ this involves the contiguous relations. The M -eigenfunction property follows by using the known analytic continuation of 2 F 1 (a, b, c; w) across the logarithmic branch cut w ∈ [1, ∞).
The above specialization of the hypergeometric function basically yields the so-called conical (or Mehler) function. To be specific, the latter can be defined by . We now consider the nonrelativistic limit of the minimal representations of R(b; x, y) derived in Section 3. We were not able to obtain a sensible limit for the second one, given by (3.48). For the remaining four, however, the limit can be handled in a sense to be explained shortly. For expository reasons we first list the resulting representations for ψ nr (λ; r, k): (4.48) We proceed to discuss these formulas. First, we note that they are derived under the assumption λ, r, k ∈ (0, ∞), (4.49) and that this implies that the integrals in (4.45), (4.47) and (4.48) are absolutely convergent. The integral in (4.46), however, is only absolutely convergent for λ > 1/2; For λ ∈ (0, 1/2] it should be viewed as a Fourier transform in the sense of tempered distributions. Second, we compare these formulas to results in [16], where a host of representations for 2 F 1 and its conical function specialization are listed. Formula (4.48) can be readily found there: It can be obtained from equation (14.12.4), which can be written (4.50) (This involves the duplication formula of the gamma function, equation (5.5) in [16].) As they stand, the three remaining representations do not occur in [16]. However, as was pointed out by a referee, they can also be tied in with results in the vast literature connected to the hypergeometric function. To begin with, formula (4.47) can be derived (with some effort) by combining equations (15.8.14) and (15.6.7) in [16]. It seems that the formulas (4.45) and (4.46) cannot be obtained by using only [16] or some other standard reference book. Even so, they agree with known results. Indeed, (4.45) amounts to equation (2.3) in the paper [29], whereas (4.46) can be derived by combining several sources. Specifically, the integral in (4.46) can be viewed as a special case of the Meyer G-function, cf. Section 16.17 in [16] and p. 144, (2) of [30]. After contour deformation, a residue calculation leads to a formula involving a linear combination of two 2 F 1 's with gamma function coefficients, cf. equations (16.17.2) and (16.17.3) in [16] or (7) in [30]. Finally, it follows from equation (3.2(27)) in [31] that the latter formula yields the conical function as represented by (4.50). Third, none of the four representations (4.45)-(4.48) has been obtained with complete rigor. The difficulty is to obtain uniform tail bounds on the pertinent integrands that allow an application of the dominated convergence theorem. (In fact, to date a similar snag has not yet been obviated for the limit transition (4.37) either.) In this connection we should add that we were unable to verify directly that each of the four formulas gives rise to a joint eigenfunction of the operators A, M andÂ with eigenvalues p 2 /4, 2 cosh(p/ µ) and 2 cosh(µx/2). On the other hand, a direct check of the joint eigenfunction properties of the five relativistic representations (3.47)-(3.51) seems not feasible either.
We continue by sketching the derivation of the four formulas (4.45)-(4.48). First, we observe that any factor of the form with t not depending on a − , can be treated via (A.25). Indeed, a scaling by a + yields In particular, this yields not only the asymptotics of the numerical prefactors, viz., but also that of the y-dependent prefactor in (3.50)-(3.51): To handle the x-dependent prefactor in (3.48)-(3.49), however, (A.25) is of no help and (A.26) seems not to apply either. But in fact we can use the G-A∆Es (A.2) to first write and then (A.26) can be invoked to deduce the dominant asymptotics For the first and second plane wave we now switch to a new variable t given by so that they turn into exp(itr) and exp(itk), resp. Likewise, in (3.47) we change z to t/µ to get dimensionless G-function arguments. With these variable changes in place, we proceed to look at the asymptotic behavior of the G-ratios in (3.47)-(3.51). For the first case this is immediate from (A. 26), and this easily yields (4.45) when the pointwise limit is interchanged with the integration. (As alluded to above, an L 1 -bound uniform for β near 0 is needed to make the interchange rigorous. As well as in the next cases, no such bound is available for now.) For the second case (3.48), it seems not possible to get from the pointwise behavior of the integrand (with x > 0) as β goes to 0 a factor exp(−πx/β ) that takes care of the diverging factor exp(πx/β ) coming from the prefactor (4.57). By contrast, for the third case we may and will make a shift with c ∈ R chosen so as to stay away from poles while shifting and taking β to 0. This cancels the diverging factor and results in (4.46) via (A.25). Finally, an application of (A.25) and (A.26) leads to the limits (4.47) and (4.48), resp. To conclude this subsection, we point out that in view of (2.39) and (1.2) the nonrelativistic limit of E(b; x, y) is given by 2w nr (λ; r) 1/2 c nr (λ; k) ψ nr (λ; r, k) =: E nr (λ; r, k), with the scattering function The latter is normalized so that it equals 1 for λ = 1, just as u(b; z) (2.53) is normalized to equal 1 for b = a ± . In this connection we would like to add that from (4.7) and (2.39)-(2.41) one readily deduces E(a ± ; x, y) = 2i sin(πxy/a + a − ). Accordingly, the 'free' theory with which the scattering is compared is given by the sine transform (and not by the cosine transform, which arises for b = λ = 0, cf. (2.58)).

Taking the relativistic Toda limit
Throughout this section, we require that the parameters a + and a − be positive. It is also convenient to require (x, y) ∈ (0, ∞) 2 , (5.1) until further notice. In keeping with our outline in the Introduction, we begin by considering the b-values (1.16). In this case the w-function (2.66) reads and hence is no longer real-valued for real z and γ = 0. By contrast, the u-function (2.53) is given by so it is still unitary for real z; moreover, it is even in γ. The Hamiltonians (2.69) can be written so they remain formally self-adjoint for real z; they are also even in γ. Next, we consider the five representations of the joint eigenfunction F(a − iγ; x, y) of the four Hamiltonians H ± (a − iγ; x) and H ± (a − iγ; y). Combining (2.68) and (1.6) with (3.47)-(3.51), these are given by (The square roots are positive for γ = 0.) As they stand, none of these representations yields a manifestly real-valued function for γ = 0. However, comparing (5.7) and (5.8), we see that these formulae are related by a complex conjugation (take z → −z in one of them to check this). Likewise, (5.6) and (5.9) are related by a complex conjugation. Since the five formulae yield the same function F(a − iγ; x, y), this function is in fact real-valued. This reality feature can be tied to the b → 2a − b symmetry of the E-function, cf. (2.48). Indeed, the latter invariance implies that E(a−iγ; x, y) is even in γ. Now since the u-function (5.3) and the phase φ(a − iγ) (given by (2.40)) are manifestly even in γ, it follows that F(a − iγ; x, y) is even in γ, cf. (2.64). Comparing once again (5.7) with (5.8), and (5.6) with (5.9), we see that these formulae are also related by flipping the sign of γ, in accord with evenness. Substituting we are now prepared to study the Toda limit Λ → ∞. First, for r ∈ R we have Thus we obtain relativistic nonperiodic Toda Hamiltonians given by (The square roots are positive for x → ∞.) Clearly, these are formally self-adjoint on L 2 (R, dx).
In this connection we point out that in view of the diverging x-shift, we may and will from now on allow x to vary over R in the Toda quantities, whereas we continue to require that y be positive. Next, we note lim Λ→∞ e δ (−2Λ)c δ (y + ira −δ /2 + η + Λ)c δ (y + ira −δ /2 − η − Λ) = e δ (2η)/4. (5.13)
To obtain the Toda limit of the joint eigenfunction F(a − iγ; x, y) involves a greater effort. The key tool is the asymptotics (A.13) of the hyperbolic gamma function. This enables us to show that the limit exists for each of the five representations (5.5)-(5.9). The details now follow.
To start with, the asymptotic behavior for Λ → ∞ of the five prefactors can be assembled from the three formulae

17)
G(±(x + Λ) + ia) ∼ exp(παa(x + Λ)). (5.18) Next, consider the integrand of (5.5) with the substitutions (5.10). Two of the four G-functions are invariant, and the remaining two yield If we now combine the Λ-dependent terms coming from the prefactor in (5.5), then we see that they cancel the Λ-dependent term in (5.19). The product of the remaining terms is readily verified to be given by Finally, since we integrate z in (5.20) over R, we may shift z by y/2, yielding the limit function Turning to the integrand of (5.6), the substitution (5.10) again leaves two of the four Gfunctions unchanged, as well as the plane wave factor. The remaining G-product has asymptotics Combining this with the asymptotics of the prefactor following from (5.16)-(5.18), the Λdependent terms drop out. Taking z → −z in the resulting limit function yields Comparing this representation to (5.22), we see that each of the factors on the right-hand side is matched by its complex-conjugate. Thus, the equality of (5.22) and (5.24) is in keeping with the real-valuedness of F T (η; x, y) following from its being the limit of a real-valued function. Proceeding in the same way for (5.9), we obtain as its limit again (5.22). Of course this should be expected, since the factors on the right-hand side of (5.9) and (5.6) are related by complex conjugation. On the other hand, the equality of the limits of (5.5) and (5.9) yields a nontrivial check of the substantial limit calculations.
At face value, the representations (5.7) and (5.8) seem not to give rise to a sensible Toda limit. In fact, however, they do, but it is expedient to postpone the details. First, we rewrite the representations (5.22) and (5.24) in a more telling form, by bringing in the G-cousins G L and G R , cf. (A.27)-(A.28). Indeed, a straightforward calculation yields the equivalent representations (5.25) and we can use (C.45) with s = a/4 to compute the Fourier transforms. This yieldŝ The integral in (5.25) is therefore equal to The new representation thus obtained can be somewhat simplified by reverting to the Gfunction, and by shifting the contour down by a/4 (recall y > 0). In this way we obtain from (5.25) an alternative representation (5.32) Proceeding in the same way for (5.26), we arrive at a fourth representation, namely Comparing it to (5.32), we deduce once again real-valuedness of the function F T (η; x, y) on R 2 × (0, ∞).
Having these two new representations at hand, we can see with hindsight that they can also be obtained from (5.8) and (5.7), respectively. Indeed, when we let z → z + ia/2 − η/2 − Λ/2 (5.34) in the integrand of (5.8), then we obtain the two G-functions featuring in (5.32), times two Λ-dependent ones. If we now proceed in the same way as before, using (5.16)-(5.18) to handle the asymptotics of the prefactor, then we arrive once more at (5.32), yielding a check on the rather extensive calculations. Likewise, (5.7) gives rise to (5.33).

Asymptotic and analytic properties of F T (η; x, y)
With the various representations of the function F T (η; x, y) at our disposal, several salient features can be readily derived. First, it is remarkably easy to show from (5.22) that it has exponential decay for x → −∞ (as might be expected from the exponential divergence of the 'potential' factors in the Hamiltonians H T ± (η; x) (5.12)). To be specific, we have a bound Inspecting (5.22), it is clear that we need only show that the integral yields a function that is O(1) for x → −∞. To this end we point out that from (A.13) we have estimates and that no poles arise for real v. Hence the function v → G(v − ia/2) is bounded on R. If we now take z → z + (x − η)/2 in (5.22), then we obtain a factor G(z − ia/2) times a factor that is bounded for x, y, z, η ∈ R. Thus we can invoke the bound (5.36) on the first factor to deduce that the integral is in fact bounded for x, y, η varying over R, completing the proof of (5.35). It is also not hard to obtain the x → ∞ asymptotics. To this end we start from the representation (5.32) and follow the reasoning below (3.52). Thus we shift the contour down by ǫ, where ǫ > 0 is small enough so that only the simple poles at z = ±y/2 are passed (recall our standing assumption y > 0). The residues of the integral then follow from (A. 19), yielding a residue sum Using the G-asymptotics (A.13), it is easily verified that the remainder integral vanishes for x → ∞, so that we deduce F T (η; x, y) ∼ u T (η; y) 1/2 exp(iαxy/2) + u T (η; −y) 1/2 exp(−iαxy/2), x → ∞. The corresponding weight function is given by w T (y) ≡ 1/c T (η; ±y) = G(±y + ia) = 4s + (y)s − (y), (5.42) where we used the G-A∆Es (A.2) in the last step.
The dual counterparts of these formulae are not obvious. To begin with, we have been unable to establish the large-y asymptotics of F T (η; x, y). We conjecture, however, that this is given by Even when this can be shown, it is not clear whether the function G R (x − η) can be viewed as an S-matrix for the dual dynamics. Indeed, the dual scattering theory seems quite unusual, just as at the classical level [17]. Moreover, like the w-function w(a − iη − iΛ; x + Λ), the u-function u(a − iη − iΛ; x + Λ) has no limit for Λ → ∞.
On the other hand, the similarity transforms can also be obtained as the limits cf. (2.45) and (5.11). Note that they have holomorphic coefficients, whereas the dual A∆Oŝ (which are the counterparts of A δ (b; y)), have meromorphic coefficients. Surprisingly, the A∆Os A ± (a − iη − iΛ; y) have no limit, whereas they do have obvious Toda counterparts, namelŷ Even though w(a − iη − iΛ; x + Λ) has no limit either, there exists a function Hence, ignoring phases and quadratic exponentials, this u-function encodes the conjectured large-y asymptotics (5.43). Finally, the similarity-transformed A∆Os Thus far, we have kept x real and y positive in the function F T (η; x, y). We proceed to study its analyticity features. To this end, consider the function The weight functions w T and w T (given by (5.42) and (5.48)) are well understood from an analytic viewpoint, so we need only clarify the character of H(x, y). The first point to note is that for each of the four integral representations (5.22), (5.24), (5.32) and (5.33) of F T (η; x, y) the weight function factors on the right-hand side of (5.53) ensure that the prefactors of the z-integrals become entire functions of x and y. Indeed, this is clear from the corresponding representations We have singled out the function H(x, y), because it extends from the real x-axis and the positive y-axis (where it takes real values) to a holomorphic (i.e., entire) function in x and y. Taking this assertion for granted, the analytic character of F T (η; x, y) can be read off from (5.53). We continue by proving the holomorphy claim.
Consider first the function defined by the integral in (5.54). The integrand is a meromorphic function I(z), whose asymptotics for Re z → ±∞ readily follows from the G-asymptotics (A.13). Specifically, we get I(z) = O(exp(−αRe z[−Im x/2 + a/2 + Im z + Im y]), Re z → ∞, (5.58) Therefore, exponential decay for Re z → ∞ can be achieved by taking and for Re z → −∞ by taking Im z < −Im x/2 + a/2 − Im y. Since the two G-functions do not depend on y, this already implies that H(x, y) extends to a holomorphic function of y. Indeed, when we continue y off the positive axis, we need only move the contour R up on the right and down on the left (whenever need be) so as to retain exponential decay.
For the x-continuation there is also no problem coming from the tail ends of the contour, but we need to avoid that the contour gets pinched between the upward and downward pole sequences (cf. (A.16)-(A.17)), Turning to the representation (5.55), we can argue in the same way to conclude that exponential decay for Re z → ∞ can be achieved by taking Im z < Im x/2 + a/2 − Im y, (5.64) and for Re z → −∞ by taking Im z > −Im x/2 − a/2 − Im y.
Here we get poles for as x is continued off the real axis. Thus we should require so as to avoid contour pinching. It therefore follows that H(x, y) extends to a holomorphic function outside the half line (5.67). Combining these two conclusions, we deduce that H(x, y) extends to a holomorphic function on C 2 , as asserted. It also follows that the contour integral in (5.54) extends to a meromorphic function of x and y, with poles only at x = −ia − z kl . Likewise, the contour integral in (5.55) yields a meromorphic function with poles only at x = ia + z kl .
To conclude this account of analyticity features, we point out that the holomorphy of H(x, y) can also be derived in a somewhat different way from the two representations (5.56)-(5.57). For these cases the two downward pole sequences at z = ∓y/2 − z kl can always be avoided by moving the contour up. Now, however, another type of restriction arises from the requirement of exponential decay on the contour tails. For the integrand in (5.56) we need Im z < a/2 − Im x/2 on the right tail, but to obtain exponential decay on the left tail we must require Im x > −a. Thus we can only deduce holomorphy for Im x > −a. Likewise, for (5.57) we need Im x < a and Im z < a/2 + Im x/2 on the left tail, so we can only infer holomorphy for Im x < a.
Even so, from these two findings we can again conclude holomorphy on C 2 . Moreover, it follows that the contour integrals in (5.56) and (5.57) give rise to meromorphic functions of x and y with poles only at x = −ia − z kl and x = ia + z kl , respectively.

Joint eigenfunction properties
To complete this section, we verify that the joint eigenfunction properties have survived the Λ → ∞ limit. Due to the simpler analyticity properties of the pertinent contour integrals (compared to the hyperbolic case), this is rather straightforward. First, to show that F T (η; x, y) is an eigenfunction of the Hamiltonians H T ± (η; x) (5.12) with eigenvalues 2c ± (y), we need only show that the A∆Os A T ± (0; x) (5.44) have the latter eigenvalues on the function G R (x) −1/2 F T (0; x, y). To this end we invoke the representation (5.26). It follows from our analysis of the contour integral in (5.55) that the contour integral in (5.26) with η = 0 defines a function M(x, y) that extends to a meromorphic function of x and y with poles occurring solely for x = ia + z kl . We may write this function for x / ∈ i[a, ∞) as where we have introduced The crux is now that we have a kernel identity (This identity can be readily checked by dividing first by K T (x, z − ia −δ /2) and then using the G L -A∆Es (A.30).) Thus we obtain Shifting C up and down by a −δ /2, no poles are met, and on the resulting contours C + and C − there is still exponential decay at the left and right. Hence we get e δ (−y) Since the two integrands are now equal and the integrals yield the same value M(x, y), the eigenvalue property (5.70) follows. We point out that the kernel identity (5.71) plays a role similar to the kernel identity (3.14) of the hyperbolic case. In the Toda case, however, we can ensure that the z-poles stay at an arbitrary distance from the real axis by choosing Im x appropriately, by contrast to the v-poles in (3.15), cf. (3.16). Moreover, in the Toda case the decay properties on the horizontal tails of the contour depend on Im z, whereas the choice of Im z is irrelevant in the hyperbolic case (since the 'z 2 -terms' drop out in the |Re z| → ∞ asymptotics).
Next, we note that to prove that F T (η; x, y) is a joint eigenfunction of the dual Hamilto-niansĤ T ± (η; y) (5.14) with eigenvalues e ± (x), we need only show that the A∆OsÂ T ± (0; y) (5.46) have these eigenvalues on the function w T (y) −1/2 F T (0; x, y).
To this end we use the representation (5.32). Accordingly we introducê Recalling our analysis of the contour integral in (5.56), we see that we should require first of all Im x > −a in (5.74). Then the integral is well defined when we choose the horizontal tails ofĈ equal to R (say) on the left and having Im z < a/2−Im x/2 on the right, while the middle part is above the pole sequences at z = ±y/2 − z kl . Furthermore, the functionM(x, y) is holomorphic in x and y for Im x > −a and extends to a meromorphic function with poles at x = −ia − z kl . In view of these features, we need only provê for y varying over a square Re y ∈ (−a, a), Im y ∈ (−a, a) (say), while keeping x real. To do so, we choose the middle part of the contourĈ equal to 2ia + (−2a, 2a), and connect this part to (−∞, −3a) and (3a, ∞) in the obvious way. Then the y-shifts can be taken under the integral sign without z-poles hitting the contour. With these analytic preliminaries in place, the key algebraic point is once more a kernel identity, namely, Shifting the contour up by a −δ /2, no poles are met, and so the joint eigenvalue equations (5.76) follow.
We conclude this subsection with some remarks. It follows from the eigenvalue features just proved that the holomorphic function H(x, y) is a joint eigenfunction of the four A∆Os A T ± (x) (5.52) andÂ T ± (0; y) (5.46) with eigenvalues 2c ± (y) and e ± (x), respectively. The coefficients of the former A∆Os are entire, whereas the coefficients of the latter are meromorphic. The coefficients of the A∆Os A T ± (0; x) (5.44) are entire as well, but their joint eigenfunction G R (x) −1/2 F T (0; x, y) is meromorphic in x, with poles for x = ia + z kl .
This shows by example that the meromorphic vs. entire character of the coefficients of the type of commuting A∆O pairs at issue in this paper is compatible both with entire and with meromorphic joint eigenfunctions. In this connection, it should be noted that when the ratio a + /a − is irrational, then the only multipliers that do no destroy the joint eigenfunction property are the constants.
Another consequence worth pointing out consists of the relations Indeed, these follow from the eigenvalue equations The nonrelativistic Toda case In this section we obtain the nonrelativistic counterparts of the quantities in Section 5, along the lines laid out in Subsection 4.2 for the hyperbolic case. Thus we switch to the parameters (4.18) and momentum variable (4.19), whereas the Toda analog of (4.20) is the substitution η = 2 µ ln(βµg), g > 0. (6.1) For the operators H T + (η; x) (5.12) and A T + (η; x) (5.44), these substitutions entail the expansion First, we have Comparing this to the limit (A.33), we see that it applies to the β → 0 limit of (6.12), with the parameter λ equal to 2. Therefore, the limit of G R (x − η) equals 1, a circumstance that also explains why we get coinciding limits for H T δ and A T δ in (6.2)-(6.4), cf. (5.44).
Next, from (6.8) we see that the above substitutions imply lim β→0 G(±y + ia) a + a − = µp π sinh(πp/ µ). (6.13) Moreover, αy 2 clearly vanishes for β → 0, so it now follows that the limits of the prefactors of the four representations are all equal to the square root of the right-hand side of (6.13).
Turning to the integrand in (5.25), the plane wave exp(iαzy) becomes exp(izp/ ). Furthermore, the substitutions on the two G R -functions yield Comparing once more to the limit (A.33), we see that it now applies to the β → 0 limit of (6.14) with λ equal to 1. Therefore, a short calculation gives the limit Taking z → t/µ in the integral, we wind up with the limit function It is readily verified that (5.26) also leads to the representation (6.17) for the limit function. Proceeding with the above substitutions for (5.32), we see that we can only get convergence of the exponentials for β → 0 when we first replace the integration variable z by βw (say). Hence we get a factor β up front, and the plane wave exp(−iαzx) becomes exp(−iµwx); moreover, the quadratic exponential converges to 1 for β → 0. It therefore remains to consider the ratio β exp(2iw ln(βµg))/G 2π µ , β; − βw ± βp 2µ Scaling the G-functions by µ/2π, we can invoke (A.25) to deduce that (6.18) has β → 0 limit (2πµ) −1 exp(2iw ln(g/ ))Γ(−iw ± ip/2 µ). (6.19) Putting the pieces together, we now get the limit function (6.16) represented as The representation (5.33) also leads to (6.20). Thus we obtain two different representations for the limit function (6.16). The resulting identity  (6.17) it is obvious thatŵ nr (k) −1/2 F T nr (λ; r, k) extends from the positive k-axis to an entire function of k. Neither from (6.17) nor from (6.20) entireness in r is manifest, but this is well known (it can be inferred, e.g., from ODE theory). As already mentioned in the Introduction, the eigenfunction property for the dual operators seems not to occur in the standard sources. It is most easily checked forÂ (6.6) by proceeding as in the relativistic case. Here it follows from the kernel identity (This identity is the counterpart of (5.77).) The large-r asymptotics of F T nr (λ; r, k) is well known. It is given by F T nr (λ; r, k) ∼û nr (λ; k) 1/2 e irk +û nr (λ; −k) 1/2 e −irk , r → ∞, (6.24) whereû nr (λ; k) =ĉ nr (λ; k)/ĉ nr (λ; −k), (6.25) and readily verified from (6.20). It seems much harder to obtain the large-k asymptotics from the above representations. Assuming plane-wave behavior, a consideration of the dual operators leads to the expectation where φ ∈ [−π, π). Indeed, just like for its relativistic counterpart (5.43), this seems the simplest behavior that is consistent with the eigenvalues and dual eigenvalues. To be sure, for neither case it is a priori clear that the asymptotic behavior must involve plane waves. At any rate, a result pertinent to (6.26) can be found in the literature: An asymptotic expansion for K ip (x) with p > x > 0 occurs in [32, Section 7.13.2, formula (19)]. The dominant asymptotics does give rise to (6.26) with φ = −π/4, but an unsettling O(x −1 ) error term is present. (If dependence on x is included, then one would rather expect an O(x) error term. Indeed, x → 0 corresponds to v → ∞, a limit for which the Toda potential is exponentially vanishing.)

A The hyperbolic gamma function
The hyperbolic gamma function was introduced and studied in [33] as a so-called minimal solution of a special first order analytic difference equation. It is basically the same as Kurokawa's double sine [34], Faddeev's quantum dilogarithm [35], and Woronowicz's quantum exponential function [36]. (The precise connections between these functions are spelled out in Appendix A of our paper [37].) In this appendix we review features of the hyperbolic gamma function G(a + , a − ; z) that are used in the present paper; if need be, see [33] for proofs.
Unless specified otherwise, we choose a + , a − > 0, (A.1) and suppress the dependence of G on a + , a − . To begin with, G(z) can be defined as the unique minimal solution of one of the two analytic difference equations that has modulus 1 for real z and satisfies G(0) = 1 (recall (1.9) for the notation used here); remarkably, this entails that the other one is then satisfied as well. It is meromorphic in z, and for z in the strip no poles and zeros occur. Hence we have with the function g(z) being holomorphic in S. Explicitly, g(z) has the integral representation g(a + , a − ; z) = where the decay rate r can be any positive number satisfying r < α min(a + , a − ), (A.14) and where Defining z kl ≡ ika + + ila − , k, l ∈ N ≡ {0, 1, 2, . . .}, (A. 16) the hyperbolic gamma function has its poles at 17) and its zeros at The pole at −ia is simple and has residue lim z→−ia In view of these features, G(z) can be written as a ratio of entire functions, where E(a + , a − ; z) has its zeros at The function E(a + , a − ; z) we have occasion to employ is defined in Appendix A of I; it is closely related to Barnes' double gamma function [38]. We need two more of its properties. First, from equations (A.41) and (A.43) in I we have where z belongs to the strip S (A.3). Second, we need the A∆Es it satisfies, namely, cf. equations (A.46)-(A.47) in I. We also state two zero step size limits of the hyperbolic gamma function, which we need for taking nonrelativistic limits. The first one yields the relation to the Euler gamma function: lim κ↓0 G(1, κ; κz + i/2) exp iz ln(2πκ) − ln(2π)/2 = 1/Γ(iz + 1/2). (A.25) For the second one we need to require that z stay away from cuts given by ±i[a + /2, ∞). Then we have uniformly on compact subsets of the cut plane. For the relativistic Toda setting it is expedient to employ two slightly different hyperbolic gamma functions defined by These functions are the unique minimal solutions of the analytic difference equations The properties of the functions G R and G L just stated are easy to infer from the corresponding properties of the hyperbolic gamma function. (In Appendix A of [37] we already introduced functions S R and S L that differ from G R and G L by the shift z → z − ia.) Finally, we have occasion to use the limits which hold uniformly for z varying over arbitrary compact subsets of C. To our knowledge, these limits have not been obtained before. We present their proof in the next appendix.
B Proof of (A.33) Since we have we need only show (A.33) for g R . Our proof actually yields a stronger result, which may be useful in other contexts. To state this result, we fix and choose δ satisfying Then we shall show where the implied constant can be chosen uniformly for z varying over compact subsets of C.
A key ingredient of our proof is the comparison function we used in Subsection III A of [33] to obtain the G-asymptotics (A.13). Specifically, we first focus on the difference where I(a + , a − ; y) ≡ 1 2y a + a − sinh(a + y) sinh(a − y) − A 2 sinh 2 (Ay) .
The A-choice (B.6) is the unique one guaranteeing that I(y) has no pole at y = 0. Specifically, we easily calculate where c(a + , a − ) is a polynomial in a + and a − of degree 4. Since we let a − go to 0, we may and will assume from now on Next, we shift the y-contour up by On account of (B.10), this ensures that only the simple pole at y = iπ/a + is passed. The residue at this pole is readily calculated, yielding the representation d(z) = a + a − 2 sin(πa − /a + ) e + (−2z) + ρ(z), z ∈ S, (B.12) where ρ is the remainder integral ρ(z) ≡ 1 2i exp(−2rz) R duI(u + ir) exp(2iuz), r = πδ/a + , z ∈ S. (B.13) We are now prepared to replace z by z + λs(a + , a − ), |Im z| ≤ a + /2, (B.14) so that we get, using (B.11), d(z + λs) = a + a − a λ a + w sinh(a + w) − a 2 + /2 sinh 2 (a + w/ with the implied constant uniform on compact subsets of the strip |Im z| ≤ a + /2. Next, we claim that we have to verify this estimate.) As a consequence, we obtain A 2 a + a − g R (A, A; z + λs) uniformly for z in C-compacts. Thus (a stronger version of) our claim follows. This concludes the proof of (B.4), and so (A.33) follows as an obvious corollary.
Note that the proportionality constant can now be checked by taking the cosine transform of (C.35). We proceed by deriving a corollary of the above result, which we need in Subsection 5.1. First, we fix ν in the strip Im ν ∈ (0, a), so that G(x − ν) has exponential decay for x → ±∞, and choose µ real, so that G(x − µ) is a phase for real x. Consider now the integral α 2π In virtue of the G-asymptotics, the phase factor in square brackets converges to exp(iαx 2 /4) for µ → −∞. By dominated convergence, the integral therefore converges to α 2π 1/2 R dx exp(iαxy)G(x − ν) exp(iαx 2 /4), y ∈ R, Im ν ∈ (0, a).
Thus we have recovered the Fourier transform (A.21) in [37].