Lax matrices for a 1-parameter subfamily of van Diejen--Toda chains

In this paper, we construct Lax matrices for certain relativistic open Toda chains endowed with a one-sided 1-parameter boundary interaction. Built upon the Lax representation of the dynamics, an algebraic solution algorithm is also exhibited. To our best knowledge, this particular 1-parameter subfamily of van Diejen--Toda chains has not been analyzed in earlier literature.


Introduction
As is known, the relativistic story of the Toda systems began with the seminal work of Ruijsenaars [R90], and has received lots of attention from the outset [BR88,BR89,Su90,Su91,Su96,Su97,Su03,Su18]. For a comprehensive review of the early developments we recommend [KMZ] and references therein. From our perspective we must pay particular attention to the work [Su90], since this is the first paper that introduced non-trivial integrable deformations of the relativistic Toda systems. The deformed systems of type-I defined in [Su90] are characterized by 4 parameters, the type-II lattices contain 2 parameters, whereas the type-III deformations contain no couplings. The next major step in this direction is due to van Diejen, who derived integrable deformations of the relativistic Toda chains with 9 coupling constants [D94], including Suris' deformed systems as special cases. Since the appearance of these integrable deformations we have been witnessing a continuous development at the quantum level [D95,KT,E,Se,Cher,DE]. However, at the classical level the theory of the deformed relativistic Toda lattices still appears to be in its infancy. The only reasonable explanation of this situation is the lack of Lax representation of the dynamics for the most generic deformations. Indeed, to our best knowledge, so far only two papers have addressed the construction of Lax matrices for certain special variants of the deformed models. Essentially by a folding procedure, Ruijsenaars [R90] derived C-type and BC-type models from the translational invariant open lattices associated with the A-type root systems. Though Ruijsenaars' folded models naturally inherit Lax pairs from the translational invariant lattices, they contain no (essential) coupling parameters. The most satisfactory results in this direction can be found in [Su90], where spectral parameter dependent 2 × 2 Lax matrices are provided for Suris' type-I/II/III deformations. Furthermore, the underlying r-matrix structure is also exhibited. Still, there is a huge gap between the 4-parameter family of type-I systems and the 9-parameter family of the van Diejen-Toda chains.
Although the full picture with the maximal number of coupling constants is still out of reach, in this paper we get one step closer to the solution of this long-standing open problem by providing a detailed analysis for a particular 1-parameter subfamily of van Diejen-Toda chains, that has not been studied earlier. In order to describe these chains, take an arbitrary integer n ≥ 3 and consider the phase space P = { ζ = (ξ, η) | ξ 1 , . . . , ξ n , η 1 , . . . , η n ∈ R } = R 2n (1.1) equipped with its standard smooth manifold structure. Also, we introduce a global coordinate system on it by the family of functions q a (ζ) = ξ a and θ a (ζ) = η a (ζ = (ξ, η) ∈ P, 1 ≤ a ≤ n). (1.2) As customary in the theory of the aforementioned relativistic integrable many-body systems, the coordinates q a and θ a are called the particle positions and the particle rapidities, respectively. Regarding P as a model of the cotangent bundle of R n , it naturally carries the symplectic form As concerns the inter-particle interaction of our integrable chains, it proves handy to introduce the single variable smooth function f µ : R → (0, ∞), r → f µ (r) = » 1 + µ 2 e −r , (1.5) where µ is a real parameter. In passing we remark that for µ = 0 we have f µ (r) > 1.
(1.6) Now, let β ∈ (0, ∞) and κ ∈ R be arbitrary constants and consider the Hamiltonian + cosh(βθ n )f β (q n−1 − q n )f β (2q n )f κβ (2q n ) + κβ 2 e −2qn + 1 2 κβ 4 e −(q n−1 +qn) . (1.7) Since we shall keep fixed the 'inverse speed of light' β throughout the paper, it should not cause confusion that we suppress the dependence of the Hamiltonian on this parameter. However, parameter κ does play a decisive role in our discussion. In this respect our first trivial observation is that for κ < 0 the Hamiltonian H κ is not bounded from below. As numerical experiments show, in this case even the completeness of the Hamiltonian flow is questionable, and so the techniques we wish to present in this paper would require serious modifications. Therefore, from now on we assume κ ≥ 0. (1.8) Recalling (1.6), it is plain that this condition entails the lower estimate H κ > n. It is also clear that the classical mechanical system (P, ω, H κ ) belongs to the class of deformed relativistic Toda systems introduced by van Diejen [D94]. In order to clarify this connection, let us call to mind the Hamiltonian given in equation (37) of [D94], which depends on nine (complex) parameters: g, g 0 , g ′ 0 , g 1 , g ′ 1 , k 0 , k ′ 0 , k 1 , k ′ 1 . Now, one can easily verify that by setting g = β, g 0 = g ′ 0 = g 1 = g ′ 1 = 0, (1.9) van Diejen's multi-parametric Hamiltonian reduces to R * H κ , where R is the canonical transformation R : P → P, (ξ, η) → (−ξ, −η). (1.11) The appearance of the reflection map R is due to the exponent −r in the definition of function f µ (1.5). Indeed, the majority of the papers, including [D94], apply the plus sign convention in the exponent. In any case, we see that H κ (1.7) describes an open relativistic Toda chain endowed with a one-sided 1-parameter boundary interaction.
Having identified our models in the family of the van Diejen-Toda chains, we are now in a position to make a comparison with the deformed models appearing in [R90,Su90]. Looking at equation (6.24) in [R90], from (1.7) it is straightforward to see that in the special case κ = 0 the composition R * H 0 reproduces Ruijsenaars' C n -type Hamiltonian. On the other hand, the comparison with the deformed models defined in [Su90] is a bit more subtle. Nevertheless, the type-III lattices can be ruled out immediately, since they have boundary interactions on both sides. By inspecting the type-I/II deformations, one may recognize immediately that only the type-I lattices bear a resemblance to our models H κ (1.7), but exact identification can be made only in the special case κ = 0. Indeed, H 0 corresponds to Suris' C n -type lattice characterized by the parameter set (a 1 , b 1 , a n , b n ) = (0, 0, 0, 1) as given below equation (6) in [Su90]. Incidentally, this particular lattice coincides with Ruijsenaars' C n -type model. To sum up, it is safe to say that Lax matrices related to our chains H κ (1.7) have appeared in the literature only for κ = 0, and so the real novelty of our paper is the construction of Lax representation to the κ > 0 case. In passing we mention that at a more appropriate place of the paper, in subsection 2.2, we perform a detailed comparison with Suris' deformed models.
Though the dependence of H κ (1.7) on κ is clearly visible, we still have to convince ourselves that different values of the parameter lead to 'essentially different' Hamiltonians. Being cautious is not without reasons: it has been observed even at the level of the nonrelativistic Toda systems that by shifting the particle positions, one can introduce 'fake' coupling constants into the models. To illustrate this phenomenon in the context of our systems, for any χ = (χ 1 , . . . , χ n ) ∈ (0, ∞) n define the canonical transformation T χ : P → P, (ξ 1 , . . . , ξ n , η 1 , . . . , η n ) → (ξ 1 + ln(χ 1 ), . . . , ξ n + ln(χ n ), η 1 , . . . , η n ). (1.12) Notice that by composing H κ (1.7) with T χ we obtain (1.13) leading to the proliferation of the coupling constants in a trivial manner. Working backward, this observation can also be utilized to investigate whether parameter κ can be eliminated: given κ 1 , κ 2 ∈ [0, ∞), we say that Hamiltonians H κ 1 and H κ 2 are equivalent up to change of couplings, in notation H κ 1 ∼ H κ 2 , if there is an n-tuple χ ∈ (0, ∞) n such that H κ 2 = T * χ H κ 1 . It is evident that ∼ is an equivalence relation on the set of our distinguished Hamiltonians (1.7). Exploiting the explicit formula (1.13) one can also verify that H κ 1 ∼ H κ 2 if and only if κ 1 = κ 2 or κ 1 κ 2 = 1. (1.14) This characterization of the equivalence may seem a bit strange at first. However, if κ > 0, then with the special n-tuple κ = (κ, . . . , κ) from (1.13) we do get The content of the above line can be interpreted by saying that the small coupling and the strong coupling regimes are in duality. Utilizing this duality, we could safely assume that 0 ≤ κ ≤ 1. Moreover, if κ 1 , κ 2 ∈ [0, 1] and κ 1 = κ 1 , then by (1.14) we have H κ 1 ≁ H κ 2 , implying that we cannot get rid of parameter κ from H κ (1.7) by shifting the particle coordinates. In this sense, κ turns out to be an essential parameter. In the light of the above discussion, one may wonder why we incorporate the inverse speed of light into the couplings of the Hamiltonian H κ (1.7). Indeed, with the special n-tuple χ β = (β 2n−1 , . . . , β 3 , β), from (1.13) one sees immediately that (1.16) and so β disappears from the expressions governing the particle interactions. In spite of this apparent simplification, we keep β for two reasons. First, it makes easier to control the transition from the relativistic models to their non-relativistic counterparts. Specifically, by expanding H κ (1.7) in β, one finds can be identified with the Hamiltonian of the C n -type non-relativistic open Toda lattice. Note that with the special n-tuple χ κ = (1 + κ, . . . , 1 + κ) the shift defined in (1.12) yields κ is purely artificial. In this sense the one-sided 1-parameter boundary interaction given in H κ (1.7) has no non-trivial footprint on the non-relativistic level. Admittedly, the second reason for keeping β has its roots in wishful thinking. In [Su90] Suris observed a remarkable relationship between certain discrete time generalized Toda lattices and Ruijsenaars' relativistic Toda lattices. An essential ingredient of this picture is the correspondence between the discrete time step-size and the inverse speed of light. Thus, by emphasizing the role of β in our paper, our hope is that it may facilitate to find a discrete time generalized Toda lattice interpretation of the deformed relativistic Toda model H κ (1.7).
Having described the integrable systems of our interest, now we wish to briefly outline the content of the paper. Since H κ (1.7) is not a textbook Hamiltonian, in Section 2 we investigate the dynamics, with particular emphasis on the completeness of the Hamiltonian flow. Also, motivated by the work of Suris [Su90], in this section we write down the Hamiltonian equations of motion in a distinguished set of canonical variables. Incidentally, in this Darboux system the relationship between a particular instance of Suris' generalized Toda lattice of type-I and the deformed relativistic Toda model H κ (1.7) with κ = 0 also becomes transparent, as discussed below Proposition 2.6. The ultimate goal of Section 3 is to construct Lax matrices for the dynamics generated by H κ (1.7). First, as formulated in Theorem 3.1, we set up a so-called 'Lax triad' for the dynamics, which can be seen as a weaker form of the Lax equation. Built upon this intermediate step, an honest Lax representation of the dynamics also emerges, as summarized in Theorem 3.3. The members of the proposed Lax pair (L, A) are defined in (3.117) and (3.43), respectively. It is worth mentioning that in the most interesting (new) cases, that is for κ > 0, the Lax matrix L has lower bandwidth 2, whereas A is pentadiagonal. By exploring further the relationships among H κ , L and A, in Section 4 we provide an algebraic solution algorithm for the Hamiltonian dynamics generated by H κ . As can be seen in Theorem 4.2, and in the subsequent discussion, the time evolution of the particle coordinates can be recovered from the LDU factorization of certain exponential matrix flow constructed with the aid of the Lax matrix L. Finally, in Section 5 we discuss some open problems related to the dynamical system (1.7).

Analyzing the dynamics
In the first half of this section we address the issue of completeness of the Hamiltonian flow generated by H κ (1.7). Our analysis hinges on a time reversal argument, which is a well-known technique to the experts of the area. Nevertheless, since completeness plays a crucial role in our investigations, we present this material, too, in a concise manner.
In the second half of the section we introduce a distinguished set of canonical variables, that will pave the way to the construction of Lax matrices. Let us also note that starting from this section we shall keep the non-negative parameter κ fixed, and shall apply the shorthand notation H = H κ .
Finally, a further piece of notation: with any strictly positive integer m ∈ N we associate the finite subset (2.1) 2.1. Completeness of the flow. Recalling the phase space P (1.1), consider the smooth map T : P → P, (ξ, η) → (ξ, −η).

(2.2)
Since T is an involution, it is automatically invertible with inverse T −1 = T . Remembering the standard coordinates (1.2), it is evident that T * q a = q a and T * θ a = −θ a (a ∈ N n ). (2.3) To put it simple, T reverses the rapidities. Since the cosh function is even, the Hamiltonian H (1.7) is invariant under the reversal of the rapidities; that is, T * H = H. Notice also that for the pullback of the symplectic form (1.3) by T we get meaning that T is actually an anti-symplectomorphism. To proceed, with the aid of the Poisson bracket (1.4) we also introduce the Hamiltonian vector field X H ∈ X(P ) corresponding to our distinguished Hamiltonian function H (1.7). Namely, for its action on the family of smooth functions we employ the usual convention Taking an arbitrary point ζ ∈ P , let be the (unique) maximally defined integral curve of X H satisfying the initial condition That is, the curve is determined by the differential equatioṅ (2.8) Here and below, the dot refers to the differentiation with respect to time. Note also that the endpoints of the maximal domain are appropriate (unique) constants obeying Now, making use of T (2.2), for any ζ ∈ P consider the well-defined smooth curve (2.10) Since H is invariant under the anti-Poisson map T , it is plain that for the action of the tangent vectorċ ζ (t) ∈ T c ζ (t) P on any smooth function f ∈ C ∞ (P ) we can writė (2.11) That is, c ζ in an integral curve of X H satisfying the initial condition Thus, simply by comparing c ζ with the maximal integral curve from the maximality of γ T (ζ) we infer that and also c ζ (t) = γ T (ζ) (t) (t ∈ (−b ζ , −a ζ )).
(2.15) Notice that, on account of (2.14), for all ζ ∈ P we have (2.16) However, since T is an involution, utilizing the above inequalities we can also write (2.17) Now, keeping in mind (2.10) and (2.15), the comparison of (2.16) and (2.17) leads to the following result immediately.
As the Lemma suggests, T (2.2) is sometimes called the time reversal map. In any case, the message of the above result is crystal clear: if one wishes to analyze the Hamiltonian flow generated by a T -invariant Hamiltonian function, any property of the backward flow can be inferred effortlessly from the study of the forward flow. This observation is of completely general nature, since it applies to any Hamiltonian invariant under the reversal of rapidities. Of course, in order to derive sharper results on the flow, one has to face the peculiarities of the Hamiltonian at hand. With this in mind, we shall need estimates on the time evolution of the coordinate functions (1.2).
Starting with the particle positions q a (a ∈ N n ), from (1.4) and (2.5) we find Thus, giving a glance at H (1.7), it is straightforward that whereas for 2 ≤ a ≤ n − 1 we have Proof. According to Lemma 2.1, it suffices to show that the Hamiltonian flow generated by H is forward-complete. In other words, we must show that for each ζ ∈ P we have b ζ = ∞.
Suppose to the contrary that there is a point ζ ∈ P such that b ζ < ∞. Since H is a first integral of X H , on account of (2.24) it is obvious that ∀a ∈ N n and ∀t ∈ [0, b ζ ) we can write (2.28) That is, upon introducing the closed cube As concerns the time evolution of θ a (a ∈ N n ), from (2.26) it is plain that ∀t ∈ [0, b ζ ) we have |θ a (γ ζ (t))| < β −1 ln(2H(γ ζ (t))) = β −1 ln(2H(ζ)). (2.31) That is, with the aid of the closed cube the content of (2.31) can be rephrased as Combining this observation with (2.30), we conclude On the other hand, since b ζ < ∞, the theory of ordinary differential equations guarantees that γ ζ 'escapes' from the compact subset B ζ × C ζ ⊆ P (see e.g. [Si,Theorem 3.7]). More precisely, one can find a small positive number δ ∈ (0, b ζ ) such that contradicting (2.34). As a consequence, the assumption b ζ < ∞ must be rejected, and so the proof is complete.
2.2. New set of canonical variables. It was observed during the early developments of the relativistic Toda chains that many calculations could be made simpler in variables different from the original coordinates provided by the particle positions and rapidities [BR88,BR89]. From Suris' paper [Su90] it is also clear that these special change of coordinates are instrumental in finding the link between certain discrete time generalized Toda lattices and the relativistic Toda chains. Therefore, taking the lead of Section 5 in [Su90], we find it convenient to introduce a family of smooth functions as follows. Let whereas for 2 ≤ a ≤ n − 1 we define (2.38) Built upon the above auxiliary objects, we also define (2.39) and the map S : P → P, ζ → (q 1 (ζ), . . . , q n (ζ), p 1 (ζ), . . . , p n (ζ)). (2.40) As it turns out, these new functions (2.39) provide a convenient framework to analyze the Hamiltonian system (1.7). As an added bonus, the comparison with [Su90] also becomes more transparent.
We initiate the study of the functions (2.39) with an elementary observation: by working out their first order partial derivatives, there is no difficulty in verifying that and so the following is immediate.
Proposition 2.3. The differential 1-form n c=1 s c dq c is closed. Incidentally, since P (1.1) is topologically trivial, the closed differential 1-form featuring the above Proposition is automatically exact. In other words, there is a smooth function S ∈ C ∞ (P ), unique up to a constant, such that n c=1 s c dq c = dS. (2.42) However, from our perspective it is more important that the symplectic form (1.3) is also exact, and it can be written as (2.43) Combining this observation with Proposition 2.3, one can easily verify following result.
Proof. Note that the auxiliary functions s a are independent of the rapidities. Thus, simply by inspecting the defining formula of S (2.40), we see that it is an invertible map with inverse Since both S and S −1 are smooth, S is a diffeomorphism. Next, pulling back the coordinate functions (1.2) by S, it is also evident that S * q a = q a and S * θ a = p a (a ∈ N n ). (2.45) Consequently, taking into account (2.43), Proposition 2.3 entails which clearly shows that S is a symplectic map.
Besides the above Darboux system, in the following we shall use extensively a distinguished family of non-canonical variables, too. Namely, for each a ∈ N n we define (2.47) as well as Though it should not cause confusion that the dependence of the above functions on the parameters β and κ is not indicated, we wish to stress that these parameters are still part of our analysis. Finally, it proves also handy to introduce the shorthand notation ∆ = 1 + κ 2 w n . (2.49) Lemma 2.6. In terms of the variables (2.47) and (2.48), the Hamiltonian function (1.7) acquires the form Proof. Utilizing (1.5), (2.48) and (2.49), we see that (2.52) Recalling (2.36) and (2.37), from (2.47) we obtain Now, simply by plugging the above expressions into H (1.7), the desired formula (2.50) follows immediately.
As we promised in the Introduction, at this point we are prepared to make a thorough comparison with the deformed systems defined in [Su90]. Giving a glance at the definitions (2.47) and (2.48), from (2.50) we find immediately that (2.55) Looking at the corresponding Hamiltonian J 1 at the end of Section 3 in [Su90], it is clear that only the type-I deformations share similarity with (2.55). Indeed, by identifying the discrete step-size with the inverse speed of light, in our notations the Hamiltonian J 1 of the type-I deformations translates into where J − = a 1 β 2 e q 1 (1 + e βp 1 ) + b 1 β 2 e 2q 1 +βp 1 , (2.57) J + = a n β 2 e −qn (1 + e −βpn ) + b n β 2 e −2qn−βpn . (2.58) In the above formulas the coupling parameters a 1 , b 1 , a n , and b n control the strength of the interaction at the boundaries. Now, by comparing (2.55) with (2.56), we see that they coincide if and only if κ = 0 and (a 1 , b 1 , a n , b n ) = (0, 0, 0, 1). Consequently, the Lax matrices introduced in [Su90] are applicable for our models (1.7) only in the special case κ = 0, which case was analyzed in [R90], too. This observation fully justifies our effort in the κ > 0 cases. Apart from clarifying the connection with Suris' work [Su90], our hope is that formula (2.55) may prove to be the first step toward the discrete time integrable system interpretation of (1.7). Turning back to our analysis, our next goal is to describe the dynamics using the global canonical coordinates given in Corollary 2.5. In particular, we need their derivatives along the Hamiltonian vector field X H (2.5). However, thanks their canonicity, for all a ∈ N n we can write (2.59) That is, our task amounts to calculating the first order partial derivatives of the Hamiltonian H (2.50) with respect to the new coordinate functions.
Lemma 2.7. In terms of the functions defined in (2.47) and (2.48), the derivative of the particle position q a (a ∈ N n ) along the Hamiltonian vector field X H has the form (2.60) As concerns the functions p a (a ∈ N n ) defined in (2.39), we have whereas for 2 ≤ a ≤ n − 1 we can write (2.63) Proof. Giving a glance at (2.47), we see Turning to the functions defined in (2.48), it is also plain that ∀a ∈ N n we have whilst ∂w n ∂q a = −2δ n,a w n and ∂w n ∂p a = 0.
(2.66) Thus, simply by exploiting the explicit expression given in (2.50), the formulas featuring the Lemma can be worked out without any difficulty.
Now an important remark is in order. Similarly to the translational invariant cases [BR88,BR89,Su90], from the above two Lemmas it is transparent that the functions introduced in (2.47), (2.48) and (2.49) may greatly facilitate our calculations. Indeed, the Hamiltonian (2.50) in Lemma 2.6, as well as the derivatives featuring Lemma 2.7 are polynomial expressions of x a , x −1 a , w a and ∆. For this reason, in the rest of the paper we shall express almost all objects only in terms of these functions. With this in mind, we find it advantageous to record here the derivatives of (2.47) and (2.48) along the Hamiltonian vector field X H . As concerns x a (a ∈ N n ), it is plain that ( 2.67) Thus, thanks to the explicit expressions (2.61), (2.62) and (2.63), the control over the derivatives X H [x a ] is complete. Turning to w a (a ∈ N n ), notice that (2.68) Therefore, utilizing (2.60), we can obtain formulas for these derivatives as well. Since we shall encounter them quite often, below we make them explicit. First, clearly we have ( 2.69) Next, for 2 ≤ a ≤ n − 2 we can write and finally (2.72)

Lax matrices
Mainly for fixing our conventions, in the beginning of this section we gather some rudimentary facts from the theory of matrices. Subsequently, we associate a Lax triad and a Lax equation to the dynamics governed by the Hamiltonian (1.7).
3.1. Algebraic preliminaries. Taking an arbitrary positive integer m ∈ N, let 0 m and 1 m denote the m × m zero matrix and the m × m unit matrix, respectively. Also, for all a, b ∈ N m consider the m × m elementary matrix E a,b defined with the entries As is known, they obey the commutation relations Next, in the real matrix algebra R m×m , for each k ∈ Z we introduce a distinguished subspace As a vector space decomposition, we clearly have Regarding its compatibility with the matrix multiplication, from the above commutation relations (3.2) it arises that is actually a Z-grading. Sometimes it is called the standard grading, or the principal gradation of the matrix algebra R m×m . Note that, on account of (3.4), each m × m matrix X can be uniquely decomposed as It is plain that X 0 is the diagonal part of X, and in general X k may have nontrivial entries only on its |k|th superdiagonal or subdiagonal, depending on the sign of k. In association with the Z-grading (3.4) we find it convenient to introduce the shorthand notations for the lower triangular and the upper triangular part of X, respectively. Relatedly, for the strictly lower triangular and the strictly upper triangular part of X we shall write As a further piece of preparation, let us recall that a square matrix is called a unit (lower or upper) triangular matrix, if it is (lower or upper) triangular with ones on the main diagonal. We conclude this subsection with an important remark. In the rest of this section we shall employ the above notations only for n × n matrices. However, starting from Section 4, we shall apply the notations (3.6), (3.7) and (3.8) exclusively for (2n) × (2n) matrices.
3.2. Lax triad. As we have already stipulated, till the end of this section the notations introduced in subsection 3.1 are in effect only for quadratic matrices of size n. In particular, E a,b stands for the n × n elementary matrix featuring 1 in the intersection of the ath row with the bth column and 0 everywhere else. We shall also frequently encounter the n × n reversal matrix, which is the distinguished permutation matrix (3.9) Since R 2 = 1 n , the matrix R is invertible, and It is also evident that ∀a, b ∈ N n we have and so (3.12) Incidentally, by virtue of the above relations it is plain that ∀k ∈ Z we have (3.13) In particular, the R-conjugate of a lower triangular matrix is upper triangular, and vice versa.
Having all the necessary algebraic objects at our disposal, utilizing the functions defined in (2.47) and (2.49) we now build up the n × n invertible diagonal matrix (3.14) With the aid of its derivative along the Hamiltonian vector field X H (2.5), we also define It is plain that whereas for 2 ≤ a ≤ n − 1 we have (3.18) However, the last diagonal entry of D requires a more careful treatment. For this reason, we find it convenient to introduce the shorthand notations (3.20) Now, looking back to (2.72), we see that and so from the explicit form of ∆ (2.49) we deduce that Therefore, giving a glance at (2.67) and Lemma 2.7, for the remaining diagonal entry of D (3.15) we can cook up the formula (3.23) To proceed, let us also introduce the n × n strictly lower triangular matrix (3.24) By acting on it with X H we get Note that each entry of the above matrix can be made fairly explicit by exploiting the equations highlighted in (2.69), (2.70) and (2.71). So far we have encountered only with quadratic matrices of size n. However, in our theory the (2n) × (2n) matrices partitioned into four n × n blocks play a more prominent role. Accordingly, by a (2n) × (2n) block matrix we shall always mean a (2n) × (2n) matrix partitioned into four blocks of equal size. Thus, it proves handy to introduce the shorthand notation N = 2n.
( 3.26) What is more important, based upon the n × n matrices D (3.14) and W (3.24), we now construct the N × N block matrices that allows us to introduce one of the most important objects in our study by the invertible lower bidiagonal matrix (3.28) As will transpire, this matrix do play a distinguished role in the construction of a Lax pair for the deformed relativistic Toda model (1.7). To this end, we shall need two more matrix valued functions. First, we introduce n auxiliary smooth functions by setting whereas for 2 ≤ a ≤ n − 1 we define (3.31) Now, utilizing these functions, let us consider the N × N block matrix where the diagonal blocks are defined by the n × n tridiagonal matrices whilst the off-diagonal blocks are given by the n × n spare matrices Second, in order to save space, we introduce the shorthand notations whereas for 2 ≤ a ≤ n − 1 we define (3.39) Making use of the above smooth functions, let us also consider the N × N block matrix ( 3.40) where the pertinent blocks are built upon the n × n matrices Notice that the diagonal block (3.41) is tridiagonal. Therefore, by inspecting the offdiagonal blocks (3.35), (3.36) and (3.42), too, it is clear that in the most interesting cases characterized by the inequality κ > 0, both A (3.32) and B (3.40) are pentadiagonal matrices. Finally, to be able to formulate one of the main results of the paper, from D (3.27) and A (3.32) we construct the N × N matrix which is also pentadiagonal for κ > 0. Proof. The verification of (3.44) turns out to be quite laborious, based upon the careful study of the N × N matrix Remembering (3.28) and (3.43), Leibniz rule yields ( 3.48) where the diagonal blocks are built up form the n × n matrices whereas the off-diagonal blocks are given by Recall also that D (3.15) is diagonal, W (3.24) has non-trivial entries only in its first subdiagonal, whereas matrices A ++ (3.33) and B ++ (3.41) are tridiagonal. Therefore, invoking the decomposition (3.6), we can write 55) meanwhile the remaining two terms are given by the more complicated expressions By inspecting the simpler looking formulas displayed in (3.55), from (3.33), (3.41) and (3.24) we find It is also clear that A ++ 1 = B ++ 1 , so we deduce immediately that Ξ ++ −2 = Ξ ++ 1 = 0 n . (3.59) Turning to the much more involved formula given in (3.56), notice that for all a ∈ N n−2 we can write whereas for 2 ≤ a ≤ n − 1 we also have As a consequence, we can write Ξ ++ 0 = 0 n . So, in the graded decomposition (3.54) each homogeneous term vanishes, whence Ξ ++ = 0 n .
Having grasped the idea of the proof, for the remaining three blocks we shall highlight only the most essential steps. Turning to (3.50), for the last term on the right hand side we can write E n,n B +− R = B +− n,1 E n,n + B +− n,2 E n,n−1 . (3.66) Consequently, the graded decomposition (3.6) of the blockΞ −− takes the form Along the same lines as in the analysis of Ξ ++ , one can easily verify that actually each term in the decomposition (3.67) vanishes, and soΞ −− = 0 n follows at once. As concerns the off-diagonal block Ξ +− (3.51), from (3.2), (3.11) and (3.24) we see (3.73) Now, looking back to the explicit form of the matrix entries, we arrive at Ξ −+ = 0 n . To sum up, we see that each block of Ξ (3.48) vanishes. Therefore, turning back to the defining equation (3.46), the Theorem follows.
At this point we are in a position to compare the content of Theorem 3.1 with the related earlier works on Ruijsenaars' n-particle open relativistic Toda chain. In this respect the most important contribution is due to Suris [Su91], who constructed certain n × n bidiagonal matrices obeying equations analogous to (3.44). Following Suris' terminology, we shall call (3.44) a 'Lax triad'. Admittedly, the bidiagonal matrix defined in equation (11) in [Su91] provided a firm guide in constructing our lower bidiagonal matrix L (3.28). However, finding the correct form of A (3.43) and B (3.40) required a much greater effort on our part. Indeed, while the analogous matrices A ± and B ± given below equations (13) and (14) in [Su91] are bidiagonal, our matrices A and B are pentadiagonal for κ > 0.
Due to Theorem 3.1, matrix A defined in (3.43) is a key player in our study of the deformed relativistic Toda system (1.7). Partly for this reason it is worthwhile to make it as explicit as possible. To this end, by conjugating A (3.32) with D (3.27) we obtain the block matrix form whereas the off-diagonal blocks are given by It is plain that Ω 2 = −1 N , whence Ω is invertible and

77)
(3.80) Note that for any N × N block matrix partitioned into four n × n blocks, we have (3.82) As concerns the relationship between the Z-grading (3.4) and the inner automorphism provided by the Ω-conjugation, from the above relation and (3.13) we infer Consequently, the Ω-conjugate of any upper triangular N × N matrix is lower triangular, and vice versa.
To continue, consider the block matrix with the n × n invertible diagonal blocks Remembering ∆ (2.49), one can check that g is an invertible matrix with inverse (3.86) Note that both g and g −1 are tridiagonal for κ > 0. By inspecting the action of the Hamiltonian vector field X H (2.5) on the matrix valued smooth function g, elementary matrix multiplication yields (3.87) Recalling ∆ (2.49), we see that the diagonal entries are actually trivial, since On account of (3.21) and (3.22) it is also immediate that whence we end up with the concise expression (3.90) Having equipped with matrices Ω (3.79) and g (3.84), we are now in a position to formulate an important technical observation.
has the block matrix structure whereas the off-diagonal blocks are given by Keeping in mind the matrix entries of A (3.43) and g (3.84), we proceed with examining the above four n × n blocks of Λ, one at a time. Besides the commutation relations (3.2), during the calculations we shall use frequently the relationships (3.11) and (3.12), too. Starting with Λ ++ (3.94), from (3.75) and (3.85) it is plain that (3.98) Continuing with the last three terms appearing on the right hand side of (3.94), from (3.77), (3.78) and (3.76) we obtain As concernsΛ −− (3.95), we can easily cook up a formula for the product g ++Â−− g ++ simply by replacing symbol A ++ in (3.98) at each place withÂ −− . Furthermore, from (3.75), (3.77), (3.78) and (3.85) we obtain g ++ RA −+ E n,n = A −+ 1,n ∆ 1 2 E n,n + A −+ 2,n E n−1,n , (3.101) Now, by plugging (3.75) and the above expressions into (3.95), the majority of the terms trivially cancel, and we are left witĥ (3.104) However, from the explicit form of the matrix entries we getΛ −− = 0 n . Turning to the off-diagonal block Λ +− (3.96), notice the following: Plugging the above formula into (3.96), we obtain (3.110) Remembering (3.20) and the matrix entries of A, we find Λ +− = 0 n . Finally, for the terms appearing in the formula of Λ −+ (3.97) we can write It readily follows which immediately leads to Λ −+ = 0 n . Summarizing, we see that each block of Λ (3.93) is zero, and so necessarily Λ = 0 N . Thus, remembering the definition (3.92), the proof is complete.
Having completed the preparations, at this point we define a Lax matrix for the deformed relativistic Toda system (1.7) by the matrix valued function L = LΩL −1 Ω −1 g. (3.117) Recall that the constituent matrix L (3.28) lower bidiagonal. Thus, on account of (3.83), the product ΩLΩ −1 is an upper bidiagonal matrix, and so its inverse is upper triangular. Consequently, H = LΩL −1 Ω −1 = L(ΩLΩ −1 ) −1 (3.118) is an upper Hessenberg matrix. Since for κ > 0 the matrix g (3.84) is tridiagonal, it is plain that the proposed Lax matrix L = Hg (3.119) has lower bandwidth 2. So, in the most interesting cases L is not Hessenberg, but its strictly lower triangular part is still quite spare. (3.120) In other words, the pair of matrices (L, A) provides a Lax pair for the Hamiltonian dynamics generated by H (1.7).
Proof. Remembering (3.117), Leibniz rule allows us to write (3.121) Therefore, by exploiting Theorem 3.1, we obtain and so by invoking Proposition 3.2, the Theorem follows.

Solution algorithm
Constructing solution algorithms from the Lax representation of the dynamics is an important aspect of the theory of integrable systems (for review see e.g. [Pe]). Since this is the usual precursor for a more geometric treatment in the framework of symplectic reductions, we find it highly motivated to present a 'projection method' for the van Diejen-Toda systems (1.7), too. We achieve this goal by adapting to our deformed systems the algebraic machinery available for Ruijsenaars' relativistic Toda chains (see Section 4 in [R90]). To make this approach work, in the first half of this section we establish some useful algebraic properties of the Lax pair (L, A) with members given in (3.117) and (3.74). During the calculations we shall utilize the Z-gradation (3.6) and the notations (3.7) and (3.8), too, but this time only for N × N matrices. 4.1. Algebraic relations. Recalling (3.119), we start our investigation with the upper Hessenberg matrix H (3.118), which is essentially built upon the lower bidiagonal matrix L (3.28). Now, giving a glance at D and W (see (3.27)), from (3.82) it is plain that (4.1) whence (4.2) Consequently, the definition of H (3.118) immediately leads to the formula However, since W (3.27) is manifestly nilpotent, obeying W N = 0 N , the inverse appearing in the above equation can be calculated as Consequently, for H we obtain the graded decomposition where whereas for 0 ≤ k ≤ N − 2 we have (4.7) Of course, we still need control over the diagonal matrix Straightforward calculations show that it has the block matrix structure (4.9) where (4.10) At this point, in principle, we could easily work out explicit formulas for the matrix entries of H, too. However, the good news is that later on we shall need control over only the homogeneous parts H −1 , H 0 and H 1 . Starting with the degree −1 part of H (4.5), from (3.27) and (4.6) we see that it can be partitioned into the block matrix form (4.11) where the non-trivial n × n blocks are given by (4.14) Incidentally, from the above formulas it is apparent that H is actually an unreduced upper Hessenberg matrix, meaning that all its entries in the first subdiagonal are nonzero. As concerns the diagonal part of H, from (4.7) and (4.8) we obtain Since we are dealing here with diagonal matrices, we find effortlessly that Finally, turning now to the degree 1 part of H, due to (4.7) we can write However, since we already have explicit formulas for the matrices T and H −1 , it comes without any difficulty that where the non-trivial n × n blocks are furnished by the slightly complicated expressions whilst the lower-right-hand corner can be recovered from (4.23) Related to upper Hessenberg matrix H (3.118), we find it convenient to introduce the conjugated matrixH = ΩHΩ −1 .
(4.24) By exploiting (3.82), it is evident thatH is an unreduced lower Hessenberg matrix, and a moment of reflection also reveals that (4.25) Thus, keeping in mind (3.83), from (4.5) it follows thatH has the graded decompositioñ (4.27) In particular, making use of (3.82), we can easily cook up explicit formulas for the matrices H −1 ,H 0 andH 1 from the above given expressions for H 1 , H 0 and H −1 , respectively.
As the last piece of preparation of this subsection, with the aid of the standard coordinate functions (1.2) we introduce the matrix valued smooth function Q = diag(q 1 , . . . , q n , −q n , . . . , −q 1 ) ∈ C ∞ (P, R N ×N ).
(4.28) Proof. We start by noticing that g (3.84) has the graded decomposition g = g −1 + g 0 + g 1 (4.31) with homogeneous parts (4.32) Therefore, remembering (3.119) and (4.5), for the lower triangular part of the Lax matrix we can write L ≤0 = L −2 + L −1 + L 0 , (4.33) where (4.34) Below we shall analyze each homogeneous part separately. Starting with the degree −2 part of the Lax matrix, from (4.11), (4.13) and (4.32) we find immediately that (4.35) Simply by comparing the above matrix with A (3.74), focusing in particular on the form of the off-diagonal block A −+ (3.78), we arrive at the conclusion Proceeding with the degree −1 part of L, from (4.34) we infer that (4.37) where the diagonal blocks are given by (4.38) whilst the only non-trivial off-diagonal block takes the form (4.39) By performing the matrix operations, from (3.85), (4.11) and (4.16) we obtain (4.42) Upon comparison with the corresponding blocks of A (3.74), from the explicit formulas appearing in (3.75), (3.76) and (3.78) it is clear that Now, merging this observation with (4.36), notice that the first relationship displayed in (4.29) arises at once. Turning to the diagonal part of L, from (4.34) we see that (4.44) where the diagonal blocks are given by (4.46) Now, by exploiting the equations (3.85), (4.11), (4.16) and (4.20), notice that the above diagonal matrices can be cast into the form Therefore, taking into account Lemma 2.6, it is immediate that (4.49) establishing the first relationship in (4.30) between the Hamiltonian function H (1.7) and the Lax matrix L (3.117).
Having completed the study of the lower triangular part of L, it is clear that the upper triangular part of L −1 can be analyzed by the same technique. Highlighting only the major steps of the calculations, let us observe that on account of (3.119) and (4.25) we can write (4.50) Remembering (3.86) and (4.32) we see that the inverse of g has the graded decomposition g −1 = −g −1 + g 0 − g 1 .
(4.51) Therefore, combining this observation with (4.5), it is plain that (4.52) where the diagonal part is given by whereas the strictly upper triangular part is built upon the homogeneous parts Inspecting each homogeneous part separately, notice that meanwhile from (4.54) it is also clear that for the degree 1 part we can write Now, recalling the matrices given in (4.12), (4.13) and (4.17), the second relation displayed in (4.29) can be confirmed easily. Spelling out (4.53), let us also observe that the diagonal part of L −1 takes the form (4.57) where the non-trivial blocks are given by The point is that, by subtracting the above diagonal matrices from the corresponding diagonal matrices given in (4.47) and (4.48), Lemma 2.7 immediately leads to the second relationship in (4.30) for X H [Q].

Projection method.
In order to make the presentation simpler, it proves convenient to introduce the shorthand notation Remembering Lemma 4.1, it is evident that for the diagonal part of Y we have (4.61) Moreover, due to Theorem 3.3 we can write Theorem 4.2. Take an arbitrary point ζ ∈ P and consider the maximal integral curve of the Hamiltonian vector field X H (2.5) satisfying the initial condition γ ζ (0) = ζ. (4.65) Then one can find smooth matrix valued functions subject to the conditions l ζ (0) = u ζ (0) = 1 N , (4.67) such that for all t ∈ R the matrix l ζ (t) is unit lower triangular, the matrix u ζ (t) is unit upper triangular, and most importantly e tY(ζ) = l ζ (t)e Q(γ ζ (t))−Q(ζ) u ζ (t). (4.68) Proof. Besides the algebraic developments we have made so far, the proof of this theorem hinges on two elementary facts from the theory of homogeneous linear systems of first order differential equations. First, the domain of every maximally defined solution of such system coincides with the common domain of the (continuous) coefficients. Second, if the coefficient matrix has a special (diagonal or triangular) form, then we can infer information about the structure of the fundamental matrix solution via the Peano-Baker series, which is basically the path-ordered matrix exponential. For a nice account on these facts see e.g. [Si,Corollary 3.3] and subsection 2.1.1 in [LR]. Continuing with the proof proper, take an arbitrary point ζ ∈ P and keep it fixed. Note that, on account of the completeness result formulated in Theorem 2.2, the domain of the maximally defined integral curve (4.64) does coincide with the set of real numbers. Utilizing this trajectory, below we shall introduce three time-dependent matrices as follows.
First, remembering the matrices L (3.117) and A (3.43), let be the (unique) maximal integral curve of the differential equatioṅ with initial condition D ζ (0) = 1 N .
(4.71) Since the coefficients in the above linear system are defined for all reals, the domain of the maximally defined solution D ζ is indeed R, as we anticipated in (4.69). Furthermore, since the coefficient matrix is diagonal, it is clear that D ζ (t) is an invertible diagonal matrix for all t ∈ R.
Second, keeping in mind D ζ (4.69) and Y (4.60), consider the maximally defined solution of the differential equationl satisfying the initial condition l ζ (0) = 1 N . Third, deploying Q (4.28), too, let be the maximal solution of the differential equatioṅ subject to the initial condition u ζ (0) = 1 N . Again, since the coefficients of the linear systems (4.73) and (4.75) are defined for all reals, so are the maximal solutions (4.72) and (4.74). Notice also that the coefficient matrix in (4.73) is strictly lower triangular, while the coefficient matrix in (4.75) is strictly upper triangular. Though we are dealing here with linear systems of variable coefficients, the representation of the fundamental matrix by the Peano-Baker series ensures that ∀t ∈ R the matrix l ζ (t) is unit lower triangular, whilst the matrix u ζ (t) is unit upper triangular. In particular, both l ζ (t) and u ζ (t) are invertible.
To proceed, for all t ∈ R define the invertible lower triangular matrix and also introduce Φ (4.77) It is clear that the dependence of Φ ζ (t) on t is smooth, and Leibniz rule yieldṡ (4.78) Now, due to (4.63), for the time derivative of Y along the curve γ ζ we can write whence (4.78) entails (4.80) As for the derivative of (4.76), from the differential equations (4.70) and (4.73) we find (4.81) Therefore, bringing into play the relations (4.29) displayed in Lemma 4.1, we can write Plugging this formula back into (4.80), from the definition (4.60) we conclude Now, looking back to (4.77), it is also clear that (4.84) The last two equations entail Combining this simple formula with the definition (4.77), we end up with This equation implies that the time evolution of Y (4.60) along any trajectory is isospectral. Of course, this is what we expect from the Lax equation (4.63). However, the essential point here is that we have full control over ρ ζ (4.76) via the differential equations (4.70) and (4.73).
Given an arbitrary N × N matrix X = [X k,l ] 1≤k,l≤N ∈ R N ×N , (4.94) for any j ∈ N N let π j (X) ∈ R denote its jth leading principal minor; that is, π j (X) = det([X k,l ] 1≤k,l≤j ). (4.95) Also, introduce the notations m 1 (X) = π 1 (X) = X 1,1 and m j (X) = π j (X) π j−1 (X) (2 ≤ j ≤ N). (4.96) Utilizing the above objects, we can draw important conclusions from the above Theorem. Indeed, the most important observation displayed in (4.68) can be interpreted by saying that the matrix e tY(ζ) has an LDU factorization (or Gauss decomposition) for all t ∈ R.
As is known from the theory of matrices (see e.g. [HJ,Corollary 3.5.6]), the factors appearing on the right hand side of (4.68) are unique and for the cth diagonal entry (c ∈ N n ) of the diagonal factor e Q(γ ζ (t))−Q(ζ) we can write that 0 < e qc(γ ζ (t))−qc(ζ) = m c Ä e tY(ζ) ä . (4.97) Thus, for the time evolution of the particle positions we obtain q c (γ ζ (t)) = q c (ζ) + ln Ä m c Ä e tY(ζ) ää , (4.98) whereas the time evolution of the rapidities can be recovered from the relationships (2.20), (2.21) and (2.22). To sum up, we see that finding the trajectories of the Hamiltonian dynamics generated by H (1.7) boils down to the computation of the leading principal minors of the exponential matrix flow t → e tY(ζ) . In this sense our solution algorithm is of purely algebraic nature. Note that the above discussion nicely harmonizes with Theorem 4.2 in [R90]. Furthermore, similarly to the translation invariant case, in our analysis we made critical use of the algebraic properties of the Lax pair, as formulated in Lemma 4.1. As a matter of fact, the real reason behind our choice of gauge for the Lax pair is to maintain these algebraic relationships.

Discussion
The main achievement of our work is the construction of a Lax representation for the dynamics generated by the van Diejen-Toda Hamiltonian (1.7). By exploiting the algebraic properties of the proposed Lax pair, we could provide a solution algorithm, too. By the very nature of the subject, the presentation of the material has an inevitable algebraic flavor. However, even the algebraic part of the story may have further surprises in store. The Lax representation of the dynamics guarantees that the spectral invariants of Lax matrix L (3.117) are first integrals of the dynamics, but we are still in debt to prove that they are in involution. In the light of the earlier developments on Ruijsenaars' relativistic Toda chains, the most satisfactory step would be to provide an appropriate r-matrix structure for L. Indeed, for the translation invariant models Suris succeeded in casting the tensorial Poisson bracket of the Lax matrix into a Sklyanin bracket form. Of course, this result was instrumental in developing a geometric picture, too, in the framework of the Poisson-Lie groups (see [Su91]).
Turning to questions requiring analytic considerations, it would be highly desirable to construct action-angle variables to the systems (1.7). Besides Ruijsenaars' pioneering work [R90], in this respect the paper [CKA] also warrants mention, in which Suris' bidiagonal matrices featuring the Lax triad are directly utilized to construct action-angle map for the relativistic open Toda chains. Just as in [R90], solving this problem for our deformed systems could shed light on their scattering properties, too. Moreover, taking the lead of [R90], it would be an equally important task to uncover the dual systems associated with (1.7) in the sense of Ruijsenaars. Since the clarification of the above problems is unavoidable to complete the study of the dynamical systems (1.7), we wish to come back to these issues in later publications.
However, the real challenge would be the construction of Lax matrices for the most general van Diejen-Toda systems. In this respect we must mention the closely related Ruijsenaars-Schneider models [RS,R88], too. Indeed, one of the main themes of [R90] is the transition from the hyperbolic Ruijsenaars-Schneider model to the relativistic Toda chain by taking the so-called 'strong coupling limit'. In fact, the multi-parametric van Diejen-Toda chains were derived in a completely analogous manner from the elliptic and the hyperbolic van Diejen systems [D94]. As concerns the construction of Lax matrices for the van Diejen-Toda models along the same lines, for about two decades after their inception the main obstacle had been the very limited knowledge about the Lax representation of the classical van Diejen systems. Nevertheless, the situation has greatly improved in the last couple of years. In our paper [Pu12] we provided a Lax representation for the rational van Diejen models with the maximal number of 3 parameters. Built upon this development, we worked out a complete theory for certain 2-parameter subfamily of hyperbolic van Diejen systems [PG,Pu18], too. However, in this research domain the real breakthrough is due to Chalykh. Indeed, among many other fascinating results, in the beautiful recent paper [Chal] a quantum Lax matrix is constructed for the elliptic van Diejen system with 9 coupling parameters.
At this point we must mention that our experience with the van Diejen models proved to be essential related to this paper, too. In fact, we arrived at our Lax matrix L (3.28) by a laborious brute-force approach based upon a close inspection of a conjectured Lax matrix for certain 3-parameter subfamily of hyperbolic van Diejen systems (see equation (6.5) in [PG]). Fortunately, by borrowing ideas mainly from [R90,Su90], our effort resulted in the Lax pair (L, A) having nice properties: L has lower bandwidth 2, whereas A is pentadiagonal for κ > 0. Therefore, it is an easy guess that such structured matrices will play role in the Lax representation of the most general van Diejen-Toda systems, too. However, remembering vividly our struggle with the 1-parameter subfamily of deformed relativistic Toda systems (1.7), we do not advocate the idea of a trial and error approach to construct Lax matrices with more coupling constants. Rather, a systematic approach is required, with more insight. As mentioned above, in this respect Chalykh's paper [Chal] may come to our salvation. Indeed, by taking appropriate 'strong coupling limits' of his Lax matrices, we expect that Lax representation will emerge for the most general van Diejen-Toda systems as well.