Solution of Baxter equation for the $q$-Toda and Toda$_2$ chains by NLIE

We construct a basis of solutions of the scalar $\boldsymbol{ \texttt{t} }- \boldsymbol{ \texttt{Q} }$ equation describing the spectrum of the $q$-Toda and Toda$_2$ chains by using auxiliary non-linear integral equations. Our construction allows us to provide quantisation conditions for the spectra of these models in the form of thermodynamic Bethe Ansatz-like equations.


Introduction
The Ruijsenaars-Schneider q-Toda chain introduced in [17], appears as a natural q deformation of the quantum mechanical Toda chain. As show in [2], the q-Toda chain can be embedded in a larger family of models which can be interpreted as stemming from the quantisation of the classical Toda chain endowed with the second Hamilton structure. The quantisation of the spectrum of the q-Toda chain was derived in [11,14] in the form of a scalar t − Q equations by implementing the separation of variables transform for the chain. This form of the quantisation conditions is however non-efficient for any practical study of the spectrum and needs to be recast in a more convenient way. Due to the connection of the q-Toda chain, in particular in the two-particle sector, with certain operators arising as quantisation of mirror curves, there has been recently a renewed interest in describing the spectrum of the chain [9,10,12,18]. The works [12,18] describe the quantisation conditions for the spectrum by solving the finite difference scalar t − Q equation in non-explicit terms, in the spirit of the early quantisation conditions for the quantum mechanical Toda chain [5,7,8]. The approach of [9,12] takes its roots in the work [15] which proposed a heuristic scheme allowing one to provide quantisation conditions of many quantum integrable models by means of TBA-like equations which involve auxiliary solutions to non-linear integral equations. The conjectured quantisation conditions were established, for the quantum mechanical Toda chain, in [13]. The recent works [9,12] proposed a conjectural form of quantisation conditions appropriate for describing the spectrum of the q-Toda chain. These results were checked numerically in the case of a low number of particles.
In the recent work [1], the authors have obtained a characterisation of the spectrum of the Ruijsenaars-Schneider q-Toda and the Toda 2 chain by means of a scalar t − Q equation. In the present paper, we show that this scalar t − Q equation can be solved in terms of a one parameter family of Fredholm determinants, or generalisations thereof, acting on ℓ 2 (N) generalising the construction obtained in [5,7,8] for the Toda chain. We then show that this equation may as well be solved within the non-linear integral equation approach developed, for the Toda chain, in [13]. This leads to quantisation conditions for these two models of interest in the form of thermodynamic Bethe Ansatz-like equations. So far, we were not able to compare our form of the quantisation conditions with the ones obtained in [9,12]. This may be due to the fact that the works [9,12] consider the |q| = 1 case while we focus on the |q| < 1 regime.
The paper is organised as follows. Section 2 introduces the two models of interest: the q-Toda and the Toda 2 chain. The Hamiltonians of these models are given in Subsection 2.1. Subsection 2.2 recalls the main result of [1], which are the quantisation conditions for the spectrum of these models in terms of a scalar Baxter t−Q equation. This is the starting point of the analysis that is carried out in this work. Section   Here, ω 1 , ω 2 are two auxiliary complex parameters which satisfy the constraint ℑ(ω 1 /ω 2 ) > 0 so that |q| < 1. Throughout this work, we shall always consider ω 1 , ω 2 to be generic. (xn−x n−1 ) . κ 1 , κ 2 appearing above are some coupling constants. By means of the quantum inverse scattering method, these Hamiltonians can be embedded into a larger family of commuting operators H = H • 0 , H • 1 , . . . , H • N , with • = q − Toda or Toda 2 . Furthermore, elements of H commute with their respective duals which are obtained by the substitution (ω 1 , ω 2 ) ↔ (ω 2 , ω 1 ). These a priori provide one with an independent set of operators so that, in order to have a well-posed spectral problem, one should consider the joint diagonalisation of the family H and its dual family H. Here and in the following, dual quantities will be denoted with a tilde, viz. for the Toda 2 chain. One defines analogously the dual quantities.

The Baxter equation quantising the spectrum
Irrespectively of the considered model, t(λ), t(λ) both commute with the total momentum operator Thus one may immediately reduce the dimensionality of the spectral problem by projecting onto the subspace h p 0 , h = ⊕ h p 0 dp 0 , (2.5) associated with the Fourier mode p 0 of the total momentum operator P tot . Upon projection, the operator t(λ) restricts to the operator t(λ; p 0 ) on the reduced space h p 0 . The reduced operator is expected to have a pure point-wise spectrum although this property has not yet been established. Taken the hyperbolic polynomial form of t(λ), resp. t(λ), the Eigenvalues of t(λ; p 0 ), resp. t(λ; p 0 ), can be parameterised by N complex roots τ = (τ 1 , . . . , τ N ), resp. τ = (τ 1 , . . . ,τ N ) and take the form: Any collection of roots τ ,τ parameterising the Eigenvalues of t(λ; p 0 ), resp. t(λ; p 0 ), are subject to the constraint for the q-Toda chain and for the Toda 2 chain. These constraints translate the fact that the Eigenvalues are associated with Eigenfunctions belonging to h p 0 .
We are now in position to formulate precisely the spectral problem associated with the q-Toda and Toda 2 chains. The latter consists in finding a collection of roots τ ,τ , satisfying to the constraints (2.8) or (2.9) so that there exists a joint, self-dual, entire solution q to the set of self-dual Baxter equations subordinate to the associated hyperbolic polynomials t τ and tτ where the polynomials are defined in (2.6) or (2.7), depending on the model of interest. The parameters appearing in the above equations take the form : The main difference between the q-Toda and Toda 2 chains is that, in the former model, κ depends explicitly on the zero mode p 0 and the transfer matrix Eigenvalue polynomial only grows in the direction ℜ λ/ω 2 → −∞. For further applications, it will be important to study the analyticity properties in p 0 of the solution q to the t − q equations governing the spectrum of the Toda 2 chain. We leave this to a subsequent publication.

Solution of the Baxter equation.
We shall construct solutions (t τ ,tτ , q) to the spectral problem described in Subsection 2.2 by following the below procedure: i) we first construct two fundamental, meromorphic and linearly independent, solutions q ± to the set of self-dual Baxter equations (2.10)-(2.11), this for any value of the roots τ and τ subject to the constraints (2.8), (2.9). Their expression involves certain determinants of infinite matrices, in the spirit of Gutzwiller [7,8].
ii) We show that these determinants can be re-expressed in terms of certain solutions to non-linear integral equations and that this induces a re-parametrisation of the solutions q ± in terms of a new set of variables δ,δ which can be constructed from τ andτ occurring in i). Here, we follow the strategy devised in [13] for the Toda chain model.
iii) We show that producing a linear combination solving the set of self-dual Baxter equations q(λ) = P + (λ)q + (λ) + P − (λ)q − (λ), with P ± being elliptic functions on the lattice (iω 1 , iω 2 ), and such that q is entire can only be possible if δ =δ and δ satisfy a set of Bethe equations.
iv) We show that the characterisation of the solutions to the spectral problem in terms of solutions to non-linear integral equation is complete, in that with every such solution one can associate two functions q ± which constitute a set of fundamental solutions to a set of self-dual Baxter equations (2.10)-(2.11) associated with some polynomials t τ ,tτ .
The fact that |q| < 1, |q| > 1 play an important role in the analysis developed below. In particular, taking the |q|, |q| → 1 limits, even on the level of the final formulae does not appear evident despite the fact that this limit is formally smooth on the level of the scalar t − Q equation. The extension of the techniques developed in this work so as to deal with the |q| = 1 case demands a separate investigation and will not be considered here.

The fundamental system of solutions
From now on, it will appear convenient to introduce the parameter and two constants ζ,ζ such that Note that, for the q-Toda chain, owing to the form of the constraints on τ ,τ (2.8), one may take ζ =ζ = 0, and it is this choice that will be made in the following. Finally, in the case of the Toda 2 chain, we shall restrict ourselves to the regime of parameters max a=1,2 e 2p 0 ρ ωa < 1 . 3.1.1 The infinite determinants K ± , K ± By taking suitable infinite size limits of finite size determinants, we construct in Subsection A four meromorphic functions K ± and their dual K ± on C. The functions K ± • are iω 2 periodic; • have simple poles at {τ k ∓ iω 1 N + iω 2 Z}; • given λ = ixω 1 + iyω 2 , they have the asymptotic expansion • solve the second order, finite-difference, equations: The hyperbolic polynomials t τ are given by eq. (2.6) or eq. (2.7), depending on the model of interest. We do stress that the parameters τ defining the polynomials t τ (λ) are taken arbitrary here, it particular they don't necessarily correspond to a parametrisation of an Eigenvalue of the transfer matrix t(λ).
It seems also worthy to stress that, in the Toda 2 chain case, the finite difference equations satisfied by K ± are not the dual ones of those satisfied by K ± .
It is easy to see that H, resp. H, have simple poles on the lattice τ k + iω 1 Z + iω 2 Z . Thus, there exists constants h, h and N roots δ = (δ 1 , . . . , δ N ) along with their dualsδ = (δ 1 , . . . ,δ N ) satisfying, independently of the model, to the constraints such that one has the zero/pole factorisation Here, given any set of N -parameters λ = (λ 1 , . . . , λ N ), we introduced the shorthand notation where θ(λ) is closely related to the Jacobi θ 1 function and θ to its dual. See Appendix D.1, eqn. (D.2) in particular, for more details. δ,δ are representatives of the lattice of zeroes of the quasi-elliptic functions H, H. In fact, by using the infinite determinant representation for the functions K ± , K ± , one may interpret H and H as being double-sided infinite determinants of Hill type. This interpretation is rather direct for the q-Toda chain but demands some work in the Toda 2 chain.
Finally, when the modular parameters (ω 1 , ω 2 ) and coupling constants κ 1 , κ 2 satisfy the reality condition ω 1 = ω 2 and κ a ∈ R, the two Hill determinants are related as what is a direct consequence of the definitions (3.13), (3.14) and the conjugation properties for K ± (3.12).

The solutions q ±
We are now in position to present the explicit form of the fundamental system of solutions q ± to the set of dual Baxter equations (2.10)-(2.11). The formulae describing q ± involve p-infinite product (z; p) for |p| < 1 and it is convenient to introduce the compact notations (w(λ); p) δ = N k=1 (w(λ − δ k ); p) for any function w . We refer to Appendix D.1 for more details on these functions. Finally, recall that to any collection of variables τ , τ one can associate, through the Hill determinant construction, two constants h, h and a collection of variables δ, δ. Then, we introduce the functions (3.21) The building blocks of Q ± are constructed with the help of the functions K ± and K ± as: (3.23) The expression for the functions f (±) p 0 (λ) depends on whether one considers the q-Toda or the Toda 2 chain. In the q-Toda case, they take the form while, for the Toda 2 chain, they read p 0 solve the finite difference equations (C.12) or (C.11), depending on the model. By using the finite difference equations satisfied by K ± , K ± and elementary transformation properties under iω 1 , iω 2 shifts for the q 2 , q −2 products and the θ, θ functions, one readily checks that q ± defined above solve eqns. (2.10)-(2.11).

Some heuristics leading to the construction of q ±
We now briefly discuss the role played by the various building blocks of the functions q ± . The two solutions to the set of dual Baxter equations (2.10)-(2.11) can be factorised as The prefactors involving θ-function correspond to a particular choice of a quasi-elliptic function. multiplying a given solution. They should be thought of as corresponding to a particularly convenient normalisation of solutions. The functions q (0) (3.26) correspond to the solutions to the "asymptotic" t − q equations, where the second (+) or the first (−) term in the set of self-dual t − q equations (2.10)-(2.11) has been dropped. Such equations describe the spectrum of appropriate open Toda chains, and take the form (3.27) as well as The functions K ± , resp. K ± , correct the behaviour of the full solution for the iω 1 , resp. dual iω 2 , regime so that one may hope to obtain a solution to the full self-dual t − q equations (2.10)-(2.11).

The auxiliary functions Y δ and Yδ
For further purpose, we single out two strips in the complex plane: and its dual B = z ∈ C : z = ixω 1 + iyω 2 (x, y) ∈] − 1/2; 1/2[×R . (3.30) They are depicted in Fig. 1. Further, we observe that the Hill determinants H(λ) and H(λ) can be decomposed as and where we have set It is easy to see that the functions λ → K ± λ ∓ iω 1 /2 /u ± λ − iω 1 /2 have only simple poles and that these form the lattice δ k ∓ iω 1 /2 ∓ iNω 1 + iZω 2 . Furthermore, owing to the asymptotic expansions of K ± , one has that for λ = ixω 1 + iyω 2 ∈ B and where 1 Toda 2 = 1 for the Toda 2 model and 1 Toda 2 = 0 for the q-Toda chain. Thus, there exists some x 0 ≥ 0 such that these functions have no zeroes in the region corresponding to ±x > x 0 .
Thus, one may introduce a curve C δ;ρ , starting at ix ′ 0 ω 1 −iω 2 /2 and ending at ix ′ 0 ω 1 +iω 2 /2, for some x ′ 0 , that divides the strip B in two disjoint domains B ± such that ±it ∈ B ± when t > 0 and large enough. The curve is oriented in such a way that B + is to its left, see Fig. 1. The curve C δ;ρ is chosen in such a way that: • the poles and zeroes of • the poles and zeroes of This property ensures that there exist holomorphic and iω 2 periodic determinations of the logarithms 35) which, owing to (3.34), are bounded functions in B ± . In fact one can even choose the determination so that l + (ixω 1 + iyω 2 ) → 0 for x → +∞, and uniformly in y ∈ R.
Further, we introduce two functions on C δ;ρ : It is easy to see that V δ admits a continuous determination of its logarithm on C δ;ρ given by and that one has the explicit representation .
The latter is a consequence of the identity Finally, the functional relation between the Hill determinant and K ± entails that Y δ and V δ are related as .

(3.41)
Analogously, one introduces the dual quantities. We only insist on the differences occurring in the dual setting. There, one has the asymptotic expansions for λ = ixω 1 + iyω 2 ∈ B. One then defines the curve Cδ ;ρ analogously to C δ;ρ , although one should pay attention to the fact that the dual curve has to be oriented so that the domain B + is to its left, see Fig. 1. One likewise defines dual, iω 1 periodic, holomorphic determinations of the logarithms 43) which, owing to (3.42), are bounded functions in B ± . In fact, one can even choose the determination so that l − (ixω 1 + iyω 2 ) → 0 for y → −∞, and uniformly in x ∈ R.
Define two kernels K, K as Then, the functions Y δ introduced in (3.36) and its dual counterpart Yδ satisfy to the below non-linear integral equations: and its dual The logarithms in the integrand are to be understood as being given by the determinations of ln V δ (λ) and ln Vδ (λ) that have been provided above. The non-linear integral equation property can be established as follows. Let λ ∈ C δ;ρ and set (3.47) By using the expression (3.38) and the fact that λ → l ± (λ) are analytic in B ± , iω 2 periodic and bounded when λ → ∞, λ ∈ B ± , one may recast Ψ in the form Here ∂B ± stands for the canonically oriented boundary of B ± . Indeed, the contribution of the lateral boundaries does not contribute to each of the integrals due to the iω 2 periodicity of the integrand while the boundaries at infinity produce vanishing contributions since, for is meromorphic on B and, for λ ∈ C δ;ρ , has simple poles at τ = λ ± iω 1 ∈ B ± with residues ±1/(2iπ), one concludes that hence the claim. The dual case is dealt with in much the same way.

Integral representations for K ± , K ±
To start with, one introduces two auxiliary Here, V δ is as defined through (3.37) and ln V δ is the appropriate determination for its logarithm on C δ;ρ given in (3.38). Similarly, one introduces the dual functions on These functions are closely related to K ± and K ± . Indeed, it holds for λ ∈ C δ;ρ that and, similarly for the dual analogues, given λ ∈ Cδ ;ρ , We remind that the functions u ± , resp. u ± , have been defined in (3.32), resp. (3.33). Note that (3.54)-(3.55) hold in fact everywhere by meromorphic continuation. These factorisations can be established by residue calculations analogous to those used in the proof that Y δ and Yδ solve the non-linear integral equations (3.45) and (3.46). We only specify that the +1 term in the integrands are chosen in such a way that (3.56) This allows one to ensure that when transforming C δ;ρ , resp. Cδ ;ρ , into ±∂B ± , resp. ±∂ B ± , the boundary at infinity of −∂B − , resp. ∂ B + does not produce additional contributions stemming from the fact that l − (τ ), resp. l + (τ ), does not approach zero when τ → ∞ with τ ∈ B − , resp. with τ ∈ B + .
Various properties of the functions v ↑/↓ which follow from their definition (3.50)-(3.51) and the fact that Y δ solve (3.45) are obtained in Subsection B.3. In particular, it is shown there that λ → v ↑ (λ) has simple poles at δ k − iN * ω 1 + iZω 2 with non-vanishing residues. Likewise, λ → v ↓ (λ − iω 1 ) has simple poles at δ k + iN * ω 1 + iZω 2 with non-vanishing residues. Although the pole property is rather clear on the level of the representations (3.54), the fact that these are genuine poles, viz. that the residues are non-vanishing, is, however, not directly apparent.
It is interesting to observe the form of the finite-difference equations satisfied by the v ↑/↓ and their duals. By using their relation to K ± , K ± and the finite difference equations satisfied by the latter one gets that The dual equations follow by the duality transformation. One should note that, in these equations, the parameters τ and δ are not independent but related through (3.16), viz. the δ a 's are the representatives of the lattice of zeroes of the Hill determinant that is built up from the parameters τ a . Observe also that these equation allow one to express the polynomial t τ (λ) solely in terms of v ↑/↓ and of the polynomials t δ . By using the Wronskian relation satisfied by v ↑/↓ (B.37), one gets that An analogous equation holds for the dual case.

New representation for the functions ψ ± , ψ ±
The results of the previous analysis allow one to recast the building blocks ψ ± , ψ ± in terms of quantities solely depending on the parameters δ,δ and hence drop the connection with the parameters τ ,τ . These functions are expressed as From the preceding discussions, it follows that: • ψ ± is entire • ψ ± does not vanish on the lattice δ k + iZω 1 + iZω 2 as follows from the properties of v ↑/↓ discussed earlier on; • Dual properties hold for ψ ± ; The parametrisation of q ± 's building blocks in terms of the functions v ↑/↓ , v ↑/↓ allows one to prove that q ± are two linearly independent solutions of the set of dual Baxter t − q equations. The linear independence is a consequence of q ± not sharing a common set of zeroes. Indeed, let Z, resp. Z, be the set of zeroes of V δ , resp. Vδ, and let Z ± = Z ∩ B ± and Z ± = Z ∩ B ± . The results obtained in Subsection B.3 of the appendix ensure that q ± has zeroes at Z ± ± iω 1 /2 + iZω 2 ∪ Z ± ± iω 2 /2 + iZω 1 . These sets are obviously disconnected. It remains to show that Z − is non-empty. This is a consequence of the fact that e Y δ (λ) has poles at δ k ± iω 1 /2 ± iN * ω 1 + iZω 2 .

Quantisation conditions in the form of Bethe equations for the δ k 's
For the purpose of stating the result, it is convenient to introduce two functions Let q be an entire solution of the set of self-dual Baxter equations (2.10)-(2.11) associated with the polynomials t τ and tτ . Let δ andδ be the representatives of the zeroes of the associated Hill determinant (3.13) and its dual (3.14), such that the coordinates δ k are pairwise distinct modulo the lattice iω 1 Z + iω 2 Z Then it holds that δ =δ and, up to an overall inessential constant, q can be expressed in terms of the elementary solutions q ± (3.20) as (3.62) The set of N + 1 parameters (δ, ξ) satisfies to the constraint as well as to the Bethe equations which take the form (3.65) These equations are expressed in terms of the quantum dilogarithm defined in (D.15) . Prior to giving a proof of this statement, some remarks are in order.
i) The Bethe equations (3.64), (3.65) are manifestly modular invariant. Furthermore, if the reality condition ω 1 = ω 2 and κ a ∈ R holds, then the Bethe equations are compatible with real roots δ l and ξe Ω 2 p 0 being of modulus one. Indeed then, due to (D. 16), it holds Also in this situation, one can readily check that C δ;ρ = Cδ ;ρ , where the complex conjugation also takes care of the change of orientation. Further, assuming that δ ∈ R N , one obtains that exp Y δ (λ) = exp Y δ (λ) as follows from the explicit expression (3.39) for these quantities, the conjugation properties of the zeroes of the transfer matrices Eigenvalues τ =τ and those of the Hill determinants (3.18) and of K ± (3.12). From there, one deduces that one may choose branches appropriately so that ln[V δ (λ) = ln[ V δ (λ), and thus, by taking complex conjugation, one gets that I δ (λ) = − I δ (λ). This property entails that, in such a case, one has that I δ (λ) + I δ (λ) is purely imaginary for real λ.
ii) The expressions for q ± simplify in the case when δ =δ and read where S is the double sine function. It follows immediately from this representation that, when δ =δ, q ± are manifestly modular invariant.
iii) The statement of the quantisation conditions when the coordinates of δ are not pairwise distinct is more complicated and we chose to omit it here since we believe this situation to be very non-generic. Still, we would like to stress that in the non-pairwise distinct case, one cannot a priori discard the possibility that δ =δ.
The proof of the above characterisation of the spectrum goes as follows. Since q ± is a fundamental system of solutions to the set of self-dual Baxter equations (2.10)-(2.11), by virtue of the structure of solutions to the set of self-dual t − q equations which is obtained in Subsection C.2 of the appendix, for any solution to these equations, the entire solution q in particular, the exists elliptic functions P ± (λ) such that As q is entire, the self-dual Wronskian of q, c.f. (C.4), is also an entire function. The functional equations (C.5) it satisfies then determines it completely to be W[q] = C q κ −iλ for some normalisation constant C q . On the other hand, the self-dual Wronskian can be computed explicitly by means of (C.14) leading to where Cδ is a constant defined through (C. 15)-(C. 16). Both information entail that Suppose that δ =δ and let δ k ∈ δ\δ. Then, P + (λ)P − (λ) has simple zeroes at δ k +iω 1 Z+iω 2 Z. Each zero being simple, only one of the two elliptic functions may vanish on the lattice Thence, should such a situation occur, one would have Res q(λ)dλ, λ = δ k = 0, (3.72) hence contradicting that q is entire. Accordingly, it holds that δ =δ so that: If the elliptic functions are non-constant, then they necessarily have zeroes and poles. Suppose that z is a zero of P + . Then P + vanishes on the lattice z + iZω 1 + iZω 2 . Thus, P − has poles on this lattice. However, it follows from the characterisation of v ↓ (λ − iω 1 ), resp.ṽ ↓ (λ − iω 2 ), given in Subsection B.3 that, for a bounded |λ|, does not vanish provided that n and m are large enough. Thence, it holds that the function q − P − has a pole at z − inω 1 − imω 2 for some n, m large enough. Since q + P + is regular 1 at z − inω 1 − imω 2 , this contradicts that q is entire. Thus, P + and P − are both constant, and up to an overall inessential constant one may take P + = 1 and P − = −ξ. Thus q is expressed by the below linear combination where we used that The denominator in (3.74) has simple zeroes on the lattice δ k + iω 1 Z + iω 2 Z, k ∈ [[1; N ]]. Since δ =δ, Q ± do not vanish on this lattice. This ensures that the parameters δ characterising an entire solution q have to satisfy to the constraints (3.76) 1 If z ∈ δ + iZω1 + iZω2, then the pole at z of q+ is compensated by the zero of P+.
for k = 1, . . . , N and any (n, m) ∈ Z 2 . By using eqns. (C.8)-(C.10) which ensure thať In other words, all constraints (3.76) will be satisfied, provided that they holds for n = m = 0 and k = 1, . . . , N . Thus, q defined as q = q + − ξq − will be entire, provided that it holds A straightforward calculation shows that Here, 1 Toda 2 = 1 in the case of the Toda 2 model and 1 Toda 2 = 0 for the q-Toda chain.
To conclude, we remind that the additional constraint (3.63) on the roots δ a is already present from the very start of the analysis, c.f. (3.15).

Equivalence of T Q equations and non-linear integral equations
In this section we show that the description of the spectrum through non-linear integral equations is fully equivalent to the original, dual t − q equation based description ??. Namely, we show that starting from a given set of solutions Y δ , Yδ to the appropriately defined nonlinear integral equations (3.45), (3.46), one can introduce two functions q ± . These functions enjoy certain property that allow one to produce two combinations of these functions giving rise to two hyperbolic polynomials t τ andtτ whose roots depend on the sets of original roots δ,δ arising in the non-linear integral equations for Y δ , Yδ. These polynomials are such that the q ± solve the set of dual t − q equations subordinate to the polynomials t τ andtτ . In the discussions to come, we will only consider the case of generic, in particular pairwise distinct sets of roots δ,δ. However, all properties continue on to the non-generic case, upon evident modifications due to multiplicities, by taking limits. Next, we assume that one is given two solutions Y δ and Yδ to the non-linear integral equations (3.45), (3.46). Since, these non-linear integral equations are the starting point of the analysis, some care is needed in discussing the contours C δ;ρ , Cδ ;ρ . These are defined as follows. C δ;ρ is a contour in B starting at ixω 1 −iω 2 /2 and ending at ixω 1 +iω 2 /2 for some well chosen x ∈ R. This entails that C δ;ρ cuts B, in two disjoint domains B ± such that ±itω 1 ∈ B ± for t > 0 and large. The contour C δ;ρ is oriented so that the domain B + is to its left.
Thus, the non-linear integral equations are joint equations for the unknown functions Y δ and Yδ and the contours C δ;ρ and Cδ ;ρ .
Several comments are in order i) It is shown in Appendix B.1 that, provided |ρ| is small enough, there exist a unique solution Y δ , resp. Yδ, to (3.45), resp. (3.46), and that, for |ρ| small enough, one can take the contours C δ;ρ , resp. Cδ ;ρ , to be ρ independent.
iii) Appendix B.2 discusses the analytic properties of Y δ and Yδ, solving According to Subsection B.3, λ → v ↑ (λ) extends to a meromorphic function on C. It is non-vanishing in B + and has simple poles at δ k − iN * ω 1 + iZω 2 with non-vanishing residues. Likewise, λ → v ↓ (λ − iω 1 ) extends to a meromorphic function on C. It is non-vanishing in B − and has simple poles at δ k + iN * ω 1 + iZω 2 with non-vanishing residues.

The functions of main interest
Finally, let (4.4) The building blocks of Q ± are constructed with the help of the functions v ↑/↓ andṽ ↓/↑ : It will be of use for the following, to establish a few overall properties of the functions Q ± and its building blocks which follow from properties of solutions to the non-linear integral equation for Y δ and Yδ: • ψ ± is entire and does not vanish on the lattice δ k + iZω 1 + iZω 2 as follows from the properties of v ↑/↓ discussed earlier on; • dual properties hold for ψ ± ; • this ensures that Q ± are both entire functions.

t τ andtτ are trigonometric polynomials.
Let q ± be as given by eq. (4.3) with δ,δ satisfying to (4.1). In addition, in the Toda 2 case, assume that ρ ω 1 e 2ω 1 p 0 < 1 and ρ ω 2 e 2ω 2 p 0 < 1 . Then, there exists N roots τ = {τ a } N 1 and N dual rootsτ = {τ a } N 1 such that The roots {τ a } and {τ a } satisfy to (4.10) In the Toda 2 case, we are only able to establish the hyperbolic polynomiality of the quantities appearing in the rhs of (4.8)-(4.9) under the hypothesis (4.7). We have no idea whether the upper bound in this hypothesis represents some genuine limitation of the method or whether it may be omitted without altering the conclusion of the above statement by a suitable modification of the form of the non-linear integral equation.
In order to establish the statement, we denote by T δ (λ) the rhs of (4.8). Then a direct calculation shows that andW can be computed in a closed form; the explicit expression can be found in (C.8). Also, we have introduced Note that the functions ψ + and ψ − could have been simplified between the numerator and denominator in (4.11) by using their iω 1 periodicity. The functions Q (ω 1 ) ± are entire functions what entails that T δ is a meromorphic function whose poles correspond to the zeroes of This expression in a consequence of (C.8) and the iω 1 periodicity of ψ + and ψ − . Thus, the potential poles of T δ (λ) can only be located at the zeroes of θ δ (λ), viz. belong to the lattice δ k + iω 1 Z + iω 2 Z. It follows from the properties justified at the end of Subsection 4.1 that the functions Q (ω 1 ) ± do not vanish on this lattice. Writing out explicitly thatW ω 1 [Q . (4.17) In particular, with ℓ = 2 this entails thatǓ ω 1 [Q + , Q − ](δ k + i(n − 1)ω 1 + imω 2 ) = 0. Since δ k + i(n − 1)ω 1 + imω 2 is a simple zero of the denominator in (4.11), this entails that T δ is entire.
One can further simplify the expression for T δ and recast it as From there it becomes clear that T δ is iω 2 periodic and, owing to the asymptotic behaviour of v ↑/↓ at ∞ which is established in Sub-section B.3 of the appendix, we infer that (4.19) Here, ǫ > 0 is arbitrary but small enough, C ǫ is some ǫ-dependent constant and D z,ǫ = s ∈ C : |s − z| ≤ ǫ . (4.20) Since T δ (λ) is entire, it admits the integral representation Thus, immediate bounds based on this representation and (4.19) leads to Thus, upon replacing C ǫ ֒→ C ǫ /ǫ, the bounds (4.19) hold, in fact, on C. By the hyperbolic variant of Liouville's theorem, we do get that there exist N roots τ k such that T δ = t τ . Furthermore, upon looking at the leading e    (1)) and recalling the explicit expression for ρ in the Toda 2 model. The reasonings in the dual case are quite similar, with the sole difference that, if T δ (λ) denotes the rhs of (4.9) one gets andW The ω 2 reduced Wronskian appearing in (4.25), can be expressed aš The rest of the reasonings in the dual case are formally the same.

The Baxter equation associated with q ±
Let q ± be as given by (4.3) and let the hyperbolic polynomials t τ andtτ be as given by (4.8)-(4.9). Then q ± provide one with the two linearly independent solutions of the self-dual t − q equations (2.10)-(2.11) associated with these polynomials. The proof of the statement goes as follows. The fact that q ± are linearly independent is a consequence of q ± not sharing a common set of zeroes. Indeed, let Z, resp. Z, be the set of zeroes of V δ , resp. Vδ, and let Z ± = Z ∩ B ± and Z ± = Z ∩ B ± . The results obtained in Subsection B.3 of the appendix ensure that q ± has zeroes at Z ± ± iω 1 /2 + iZω 2 ∪ Z ± ± iω 2 /2 + iZω 1 . These sets are obviously disconnected. It remains to show that Z − is nonempty. This is a consequence of e Y δ (λ) = 1 + o(1) when λ → ∞, λ ∈ B + , and of the fact that e Y δ (λ) has poles at δ k ± iω 1 /2 ± iN * ω 1 + iZω 2 .
Next, introduce the iω 1 -Wronskian It holds, The last two terms in each bracket cancel out. It is discussed in Subsection C.1 of the appendix that W ω 1 [q + , q − ](λ) satisfies to the finite difference equation . Thus, the handlings can be simplified further, leading to (4.31) The closed formula (C.6)-(C.8) for W ω 1 [q + , q − ] ensures that this Wronskian is a meromorphic function on C that, furthermore, is not identically zero. It can thus vanish at most on a locally finite set. Thus, it holds that q + solves the iω 1 t − q equation everywhere, with the exception of an at most locally finite set and away from its pole. By continuity, it satisfies this equation on C \ {δ k + iZω 1 + iZω 2 }. The reasoning is much the same regarding to the dual, iω 2 t − q equation, as well as for q − . We leave the details to the reader.

Spectrum of the transfer matrices
In this last section, we establish how to reconstruct the spectrum of the transfer matrix and its dual from the knowledge of the parameters δ describing a given joint Eigenvector of these matrices. One could, of course, use the relation (3.59). However, the latter seems rather impractical from the point of view of concrete applications. The reconstruction we propose below is rather direct and relies on objects that are natural from the point of view of the non-linear integral equation approach.
Let δ be a collection of pairwise distinct parameters solving the Bethe equations (3.64) or (3.65). Assume that the contours C δ;ρ and C δ;ρ are such that the roots τ andτ associated with the Eigenvalues t τ of t(λ) andtτ of t(λ) are such that Then, given V δ , resp. Vδ, as defined in (4.2), it holds where k = 1, . . . , N − 1.
Prior to discussing the proof of this representation, some remarks are in order.
i) The integral representation reconstructs the Newton polynomials iii) The hypotheses (5.1) may appear as a posteriori data which are hard to verify a priori. However, one should note that, taken the definition of H, for |ρ| small enough, the sets δ will be close to the poles of H, viz. the τ . Hence, the very definition of the contour C δ;ρ ensures that hypotheses (5.1) does hold. For greater values of |ρ|, the validity of (5.1) is less clear, although, in practice, one may always reach the large values of |ρ| by first starting from |ρ| small enough, and then deforming the parameter. That would then allow one to find, for the value of ρ of interest, p k , p k ∈ Z such that The reconstruction (5.2)-(5.3) would then still hold, provided one makes the substitution in the corresponding formulae.
We now establish the result. We only focus on the Toda 2 model, and in that case solely establish eqn. (5.2). All the other cases are dealt with in a similar fashion, up to small, evident, modifications. Let Upon using that ln V δ is smooth on the integration contour and that it is iω 2 periodic there, one may integrate by parts and, upon using that the α k 's are also iω 2 periodic one gets that the boundary terms cancel out, so that: The functions v ↑/↓ , resp. v ↑/↓ , built from the given solution to the non-linear integral equation define polynomials t τ andtτ . One builds then, just as in Section 3, two elementary solutions q ± . From these one builds V δ as in (3.37). Thus, one has the explicit representation Observe that one has the below bounds in the λ → ∞ regime for some constant C k > 0, as well as Thus, the decay at infinity and the iω 2 periodicity of the integral allow one to split the integral into three pieces The first two parts can be computed, at this stage, by residues and each yields 0 since the integrand is holomorphic. Relatively to the third integral, one may close the contours by virtue of iω 2 -periodicity. ln ′ H µ − i ω 1 2 has simple poles at µ ∈ {i ω 1 2 + δ k + iZω 1 + iZω 2 } with 2iπ residue +1 and simple poles at µ ∈ {i ω 1 2 + τ k + iZω 1 + iZω 2 } with 2iπ residue −1. By construction of the curve C δ;ρ , one has that δ k ± iω 1 /2 ∈ B ± . Since δ k + 3iω 1 /2 is already at distance greater than 1 along the iω 1 axis and in units of iω 1 , solely δ k + iω 1 /2 lies inside the integration contour of the third integral. Likewise, by virtue of hypothesis (5.1), solely the poles at τ k + iω 1 /2 are located inside of the integration contour. Hence, what entails the claim.

Conclusion
In this work, we constructed the explicit solution to the scalar t − Q equations describing the spectrum of the q-Toda and Toda 2 chains. The solution q was expressed through Fredholm determinants and we showed that it is equivalently expressed in terms of solutions to nonlinear integral equations. Ultimately, we obtained a set of quantisation conditions, in the form of Thermodynamic Bethe Ansatz like equations.
The proof of these properties goes as follows. We first assume that convergence holds and establish periodicity and the second order finite difference equation. We treat the case of K + as all others are dealt similarly.
We first prove the periodicity. One obviously has K (n) + (λ + iω 2 ) = K (n) + (λ) Hence, taking the limit, one gets K + (λ+iω 2 ) = K + (λ). The recurrence equations (3.6) follows by developing in respect to the first column: and then taking the limit on the level of this representation. Finally, the presence/absence of poles at some point z ∈ C follows upon taking the n → +∞ limit on the level of Res(K (n) + (λ)dλ, λ = z). Hence, it remains to establish convergence and asymptotics. Here, we treat the case of K + and K − as both determinants demands slightly different techniques to deal with.
To start with, we observe that given an n × n matrix with entries M One can represent K (n) and δ ab is the Kronecker symbol. Further, take ǫ > 0 and small enough and define where D z,ǫ is the open disk of radius ǫ centred at z. Then, uniformly in λ ∈ D (+) τ ;ǫ of the form λ = ixω 1 + iyω 2 , with x ≥ x 0 for some x 0 ∈ R, it holds, for some x 0 and ǫ-dependent constant C , that 1 t τ (λ + ikω 1 ) ≤ C e γπ ω 2 N λ q kγ with γ = 1 q − Toda γ = 2 Toda 2 . (A.14) Thence, for some constant C ′ , one can bound C[δK τ . Since the limit issues from a uniform convergence on compacts of a sequence on holomorphic functions, one gets that K + is holomorphic on D (±) τ . Finally, the asymptotics follow from the estimates on C[δK (n) + ] given above and the bound (A.11). In the case of the q-Toda chain, the case of K (n) − is dealt with by very analogous handlings. However, the Toda 2 case demands more work. In that case, one can represent K  − ab = b a δ a,b−1 + c a δ a,b+1 and the coefficients take the form for some constant C and uniformly in λ ∈ D (−) τ ;ǫ , with λ = ixω 1 + iyω 2 with x ≤ x 0 and x 0 fixed. The inverse of T can be computed explicitly, see e.g. [4] where the inversion of general tri-diagonal matrices is performed. Set α such that e −α = b. Observe that the hypothesis |b| < 1 implies that ℜ(α) > 0. Then, one has Then, it is straightforward to obtain, that, uniformly in k, ℓ ∈ [[1; n]], one has Thus, since one has what allows one to conclude regarding to convergence to K − as well as to its asymptotics along B − . It is useful to remark that, in the q-Toda chain case, all quantities K ± , K ± correspond to one parameter λ family of Fredholm determinants of id + trace class operators on ℓ 2 (N). They can thus be immediately defined as the corresponding infinite determinants, see e.g. [6]. However, this is no longer the case for the Toda 2 chain as the approximants of K − and K + do not have a standard form. Of course, upon multiplication by T −1 one still recovers an expression for the form id + rank n approximation of a trace class operator, but the structure of the obtained matrix makes the second order finite difference equation structure much less apparent.

B Non-linear integral equations: solvability and properties
Recall the definitions of the strips B given in (3.29) and B given in (3.30).
Furthermore, given a collection of points δ k ∈ B andδ k ∈ B one introduces a nonintersecting contour • Γ δ in B starting at ixω 1 − iω 2 /2 and ending at ixω 1 + iω 2 /2 for some x ∈ R. The contour Γ δ cuts B in two disjoint domains B ± such that ±itω 1 ∈ B ± for t > 0 and large. The contour is oriented in such a way that B + is on its + side, where the + side is to the left of the contour along its orientation. Finally, the contour Γ δ is chosen such that the points δ k ± i ω 1 2 ± iNω 1 are all in B ± .
• Γδ in B starting at ixω 2 + iω 1 /2 and ending at ixω 2 + iω 1 /2 for somex ∈ R. The contour Γδ cuts B in two disjoint domains B ± such that ±itω 2 ∈ B ± for t > 0 and large. The contour is oriented in such a way that B + is on its + side, where the + side is to the left of the contour along its orientation. Finally, the contour Γδ is chosen such that the B.1 Unique solvability at small ρ of the NLIEs of interest Fix R > 0 and introduce two Banach spaces Then, there exists ρ 0 > 0 and small enough such that, for any |ρ| < ρ 0 , i) for any f ∈ S R and f ∈ S R , the functions have no zeroes on Γ δ , resp. Γ δ , and have continuous principal value logarithms on these curves; ii) the non-linear integral equation has a unique solution in S R . Similarly, the dual non-linear integral equation has a unique solution in S R . We remind that K has been introduced in eqn. (3.44), while its dual is defined as K = K |ω 1 ↔ω 2 .
Here, we give the proof in the case of equation (B.3). The dual case can be dealt with in much the same way. The first statement is clear provided that |ρ| is small enough. Thence, the operator is well-defined on S R . Also it stabilises this set provided that |ρ| is small enough (its maximal magnitude depends on R). Indeed, by using that and upon taking |ρ| small enough, direct bounds lead to where We In the following, given a collection of variables δ along with their duals δ, we assume that we are given solutions Y δ and Yδ to the two non-linear integral equations appearing in the body of the paper, eqns. (3.45) and (3.46), along with their associated contours C δ;ρ and Cδ ;ρ .
Also, below, we will use the shorthand notations V δ and Vδ introduced in (4.2). Finally, we assume δ,δ to be generic and, in particular, pairwise distinct modulo the lattice iω 1 Z + iω 2 Z.
We stress that each of these simple poles is genuine meaning that it has non-vanishing residues.
The functions V δ , resp. Vδ are meromorphic on C: • V δ has simple poles at δ k ± i ω 1 2 ± iNω 1 + iZω 2 • Vδ has simple poles atδ k ± i ω 2 2 ± iNω 2 + iZω 1 and theses are their sole poles in C. Furthermore, given any ǫ > 0 small enough, Y δ has the asymptotic behaviour • q-Toda where D z,ǫ is the disk of radius ǫ and centred at z.
and e Y δ (λ) = 1 Analogous asymptotics holds for the dual case.
Thus, from the already gathered data, one gets that e Y δ (λ) can be meromorphically continued to B (2) as − \ B (1) .

(B.17)
Let Z (1) denote the zeroes of the meromorphic function V δ in B (1) . Thus, e Y δ (λ) , as a meromorphic function on B (2) , has zeroes in Z (1) ± ± iω 1 , with Z (1) ± = Z (1) ∩ B ± and simple poles in δ k ± i3ω 1 /2. Since the residues of V δ do not vanish at δ k ± iω 1 /2, c.f. eqn. (B.14), the same property holds for the residues of e Y δ (λ) at δ k ± i3ω 1 /2 since e J δ does not vanish on B (2) \ B (1) . The above also entails that V δ can be meromorphically continued to B (2) . Its only poles are simple and are located at δ k ± iω 1 (1 + p)/2 for p = 0, 1 and all have non-vanishing residues. Furthermore, denote by Z (2) the zeroes of V δ in B (2) . In its turn this information allows to meromorphically continue Y δ to B (3) . By straightforward induction, the claim follows. The situation is similar for the dual case. It thus follows that .

(B.18)
It remains to establish the statement relative to the asymptotics. The q-Toda and Toda 2 chains demand a separate treatment.
Upon making the change of unknown function and setting Let π n (γ) = n p=n 0 p (γ) . π n (γ) is bounded in n and γ ∈ g, and is uniformly away from 0. Upon changing the unknown function ϕ n (γ) = π n (γ)ψ n (γ) one infers that ψ n (γ) = ψ n 0 −1 (γ) + n k=n 0 If q 2 χ −1 > 1, then This entails the claim for the range of parameters stated in the claim.

B.3 The auxiliary functions v ↑/↓
Let Y δ , C δ;ρ and their duals be as introduced in Subsection B.2 above and let δ,δ satisfy to the constraint (4.1). Recall the definition (4.2) of V δ and Vδ. Let v ↑ (λ), v ↓ (λ − iω 1 ), resp. v ↑ (λ), v ↓ (λ − iω 2 ), be defined on C δ;ρ , resp. Cδ ;ρ , by means of the integral representations given in the core of the paper, eqns. (3.50), (3.51) and eqns. (3.52), (3.53) for their duals. Then λ → v ↑ (λ) has no zeroes and no poles in B + . It extends into a meromorphic function on C with • simple poles at δ k − iN * ω 1 + iZω 2 ; • zeroes at Z − − i ω 1 2 . Similarly, λ → v ↓ (λ − iω 1 ) has no zeroes and poles in B − and extends into a meromorphic function on C with • simple poles at δ k + iN * ω 1 + iZω 2 ; and satisfy to the Wronskian relation It is convenient to introduce the compact notation Finally, for any ǫ > 0 small enough, v ↓ has the λ → ∞ asymptotic behaviour The asymptotic λ → ∞ behaviour of v ↑ depends, however, on the model and takes the form Here 1 Toda 2 = 1 in the case of the Toda 2 model and 1 Toda 2 = 0 in the case of the q-Toda chain. Dual properties hold for the dual functions, namely for any ǫ > 0 small enough, v ↑ has the λ → ∞ asymptotic behaviour The asymptotic λ → ∞ behaviour of v ↑ depends, however, on the model and takes the form The proof of this statement follows from the results obtained in Subsection B.2 and the analytic continuation formulae v ↑ (λ) = exp −

C General properties of solutions to the Baxter equation
In this section of the appendix, we discuss various overall properties satisfied by solutions q to the system of dual Baxter equations (2.10)-(2.11). We shall assume that ω 1 and ω 2 are generic and that they satisfy to the condition ℑ(ω 1 /ω 2 ) > 0 so that |q| < 1.

(C.3)
Cδ is a constant depending on the model.
(C. 16) Here, again, the results follow from a direct computation and the use of Wronskian relations eq. (B.37).

C.2 General form of the solutions
We have now introduced enough notations so as to characterise the form of any solution to the set of dual Baxter equations (2.10)-(2.11).