Discrete integrable systems and random Lax matrices

We study properties of Hamiltonian integrable systems with random initial data by considering their Lax representation. Specifically, we investigate the spectral behaviour of the corresponding Lax matrices when the number $N$ of degrees of freedom of the system goes to infinity and the initial data is sampled according to a properly chosen Gibbs measure. We give an exact description of the limit density of states for the exponential Toda lattice and the Volterra lattice in terms of the Laguerre and antisymmetric Gaussian $\beta$-ensemble in the high temperature regime. For generalizations of the Volterra lattice to short range interactions, called INB additive and multiplicative lattices, the focusing Ablowitz--Ladik lattice and the focusing Schur flow, we derive numerically the density of states. For all these systems, we obtain explicitly the density of states in the ground states.


Introduction
In this manuscript, we study properties of Hamiltonian integrable systems with random initial data by analysing the spectral properties of their Lax matrices considered as random matrices.
One of the first investigations in this direction was made in [14,15], where the authors considered random Lax matrices of various Calogero-Moser systems, and computed their joint eigenvalues densities. More recently, this idea gained new attention thanks to the work of Spohn [56], where the author connected the spectrum of the random Lax matrix of the Toda lattice with the sparse matrix of the Gaussian β-ensemble [21] at high temperature [6]. Later, the Lax matrix of the defocusing Ablowitz-Ladik lattice [53,55] with random initial data was connected with the sparse matrix of the circular β-ensemble [43] in the high temperature limit [37] by two of the present authors [32], and, independently, by Spohn [57]. In the same work, Spohn also considered the defocusing Schur flow [31], and he connected it to the Jacobi β-ensemble [21] in the high temperature regime [7]. Further developments on this subject were also presented in [35,49]. We mention also the work [61], where the author studied the statistical properties of the energy level of a quantum integrable system analysing the eigenvalues of the Lax operator.
Classically, from the seminal paper of Liouville [45], an integrable system is understood as a Hamiltonian system possessing enough first integrals such that the motion of the system is regular and predictable [8]. In the modern geometrical setting [10,54], we say that given a Poisson manifold P of dimension n, with local coordinates a = (a 1 , . . . , a n ), and Poisson bracket {. , .} of rank 2r, the flow generated by a scalar smooth function H = H(a): is integrable if it possesses k = n − r functionally independent first integrals in involution with respect to the given Poisson bracket. Finding first integrals is often a complicated task, and during the past decades several algorithms to construct them have been developed. One of the most effective methods to produce first integrals of a given mechanical system is the so-called Lax pair representation 1 . The concept of Lax pair originates from the work of P. D. Lax on the theory of PDEs [44], see also [3,62]. In the finite dimensional setting, the construction runs as follows [9,11] is equivalent to the Hamiltonian flow (1). Then the matrices L and A are a Lax pair for the Hamiltonian system and the matrix L is called Lax matrix. The main consequence of the Lax equation (2) is that the eigenvalues of L are first integrals of the Hamiltonian flow (1). So, provided that we can prove these eigenvalues give enough functionally independent quantities in involution, we can infer the integrability of Hamilton's equations (1) through the Lax pair. The fact that in many cases an integrable system is equivalent to a matrix relation gives the connection with random matrix theory. Indeed, when the initial data a(0) is chosen randomly, the Lax matrix L = L(a(0)) becomes a random matrix. To define a random initial data, we consider invariant measures with respect to the Hamiltonian flow. In general, such objects have the form dµ = m (a) da 1 ∧ · · · ∧ da n , where the density m(a) is such that the measure dµ ∈ L 1 (P). Hamiltonian systems have a natural invariant measure, the so called Gibbs measure [42] obtained from the Hamiltonian itself so that dµ H becomes effectively a probability measure. If the measure is not normalizable in the whole phase space P, one needs to restrict the measure to a suitable submanifold M. Because of the integrability of the system, we can consider a more general probability measure, called generalized Gibbs measure dµ G = 1 Z G e − k j=1 β j H j (a) dµ, 1 Often L-A pair in the Russian literature.
where β 1 , . . . , β k are constants and H 1 , . . . , H k are the first integrals. As above, we assume that the normalization constant Z G is finite, In this manuscript, we explore the behaviour of the spectrum of the Lax matrices of various relevant integrable systems when the number of degrees of freedom N → ∞, and the initial data is sampled according to a properly chosen Gibbs measure. Our main result is that we can compute the density of states for the Lax matrix of the exponential Toda lattice and the Volterra lattice. This is done through a one-to-one correspondence with the Laguerre β-ensemble at high temperature and with the antisymmetric Gaussian β-ensemble at high temperature, respectively. These are two known classes of random matrix ensembles, see [48] and [22,28] respectively. We consider other relevant cases of integrable systems, namely the focusing Ablowitz-Ladik lattice [1,2], the focusing Schur flow, and a class of integrable generalization of the Volterra lattice to short range interactions, called the Itoh-Narita-Bogoyavleskii (INB) additive and multiplicative lattices [17]. In these cases the corresponding random Lax matrices are not symmetric nor self-adjoint and we derive numerically their density of states that has support in the complex plane. Interesting patterns of the density of states emerge as we vary the parameters of the system. Finally, for all the integrable systems analysed in this manuscript, we are able to compute the density of states in the low-temperature limit, namely in the ground state.
This manuscript is organized as follows. In section 2 we give an introduction to the basic tools of the integrable system theory needed in this manuscript. Next, we study the density of states of the random Lax matrices of the exponential Toda lattice, the Volterra lattice, the INB lattices, the Ablowitz-Ladik lattice, and the Schur flow in sections 3, 4, 5, 6, and 7 respectively. Finally, in section 8 we give some conclusions and an outlook for future developments on this topic.

Preliminaries
In this section, we recall some standard tools to study Hamiltonian integrable systems that we need throughout the manuscript . For further details, we refer to various textbooks and monographs [9][10][11]54].
A Poisson manifold is a pair (P, {. , .}) where P is a n-dimensional differentiable manifold and {. , .} is an antisymmetric bilinear operation on the space C ∞ (P) of smooth functions over P, such that for all functions f, g, h ∈ C ∞ (P), it satisfies: The operator {. , .} is called a Poisson bracket. When there is no risk of confusion, we simply denote a Poisson manifold by P, where the Poisson bracket is assumed to be fixed and given.
In local coordinates a = (a 1 , . . . , a n ) the Poisson bracket is specified by an antisymmetric (2, 0) tensor π ij (a), the Poisson tensor, acting on the coordinates as The Jacobi identity on the coordinates is equivalent to the relation where we are summing over repeated indices. In an open subset of P the Poisson tensor has a fixed even rank 2r ≤ n. By antisymmetry, it follows that the Poisson tensor can be non-degenerate, meaning that det π(a) = 0, only if the dimension n of the base space is even, namely n = 2N . Given a function H(a) ∈ C ∞ (P), it generates a set of so-called Hamilton's equations through the relationȧ The function H itself is called a Hamiltonian. The previous set of equations defines a continuum time flow from an initial condition a(0) ∈ P to its time evolution t > 0, namely Φ t : a(0) → a(t).
In this case the quantity K is called a first integral or a constant of motion. The notion of Liouville integrability is strictly related to the number of first integrals and the rank of the associated Poisson tensor.
in a dense subset of P.
As discussed in the Introduction, random initial data are obtained from an invariant measure dµ of the form (3). More precisely, this means that the measure of every subset S ⊂ P with respect to dµ is preserved under the time-evolution Φ t , Interpreting the evolution as a coordinate transformation, we have where Φ * t (dµ) is the pull-back of dµ through Φ t . This shows that the condition (6) is satisfied if Φ * t (dµ) = dµ. In coordinates, for a measure written in the form (3), namely dµ = m (a) da 1 ∧ · · · ∧ da n , the invariance of the measure with respect to the Hamiltonian flow is expressed by the condition where div is the usual euclidean divergence, see e.g. [36,Chapter 1]. The vector field f H is specified by the Hamiltonian H via the relation (f H ) i = {a i , H}. The condition (7) can be written in the form From formula (8) we immediately have two important consequences.
• If the Hamiltonian vector field is divergence free, like in the case of a canonical Poisson bracket, it follows that the Euclidean measure is an invariant measure.
• If K is a first integral and m is the density of an invariant measure, then from the Leibniz rule (4) it follows that is the density of another invariant measure for every scalar function f ∈ C ∞ (P).
In all the examples of this manuscript, the Hamiltonian vector fields are divergence free, so we will be allowed to consider the generalized Gibbs ensemble described in the Introduction where dµ = dµ 0 the Euclidean measure in (9).
3 Laguerre β-ensemble and the exponential Toda lattice In this section, we introduce an integrable model that we call exponential Toda lattice, since it resembles the well-know Toda lattice. We construct the Lax pair for this system, and we define its generalized Gibbs measure. Finally, we compute the mean density of states of the Lax matrix. The exponential Toda lattice is the Hamiltonian system on P = R 2N with canonical Poisson bracket described by the Hamiltonian We consider periodic boundary conditions and ∆ ≥ 0 is an arbitrary constant. The equations of motion are given in Hamiltonian form aṡ The system possesses two trivial constants of motion, the first one due to periodicity, the second one due to the translational invariance of the Hamiltonian (11). In order to obtain a Lax pair for this system we introduce, in the spirit of Flaschka and Manakov [26,27,46], the variables where we notice that N j=1 y j = e − ∆ 2 . In these variables, the Hamiltonian (11) and the constants of motion (13) transform into The Hamilton's equations (12) becomė where x N +1 = x 1 , y 0 = y N . One can explicitly construct a Lax pair for this system. Let us introduce the matrix E r,s , defined as (E r,s ) ij = δ i r δ j s . Set where, accounting for periodic boundary conditions, indices are taken modulo N , so that E i,j+N = E i+N,j = E i,j for all i, j ∈ Z. For example, the matrix L in (15) has the explicit form The system of equations (14) then admits the Lax representatioṅ Hence, the quantities H m = Tr L m−1 , m = 2, . . . , N + 1 are constants of motion as well as the eigenvalues of L. For the exponential Toda lattice, we define the generalized Gibbs ensemble as where β, η, θ > 0, the Hamiltonians H E , H 0 and H 1 are defined in (11) and (13) respectively, Z E N is the normalization constant, dr = dr 1 . . . dr N and analogously for dp. We notice that according to this measure, all the variables are independent, moreover all p j are identically distributed, and so are the r j . After introducing the variables (r, p) → (x, y), the previous measure turns into Let χ 2α be the chi-distribution, defined by its density here Γ is the classical gamma function [19, §5]. Then, the variables x j and y j in the Gibbs measure (17) are independent random variables with scaled chi-distribution, (16) becomes a random matrix when the entries are sampled according to (17). Such random matrix can be linked to the so-called Laguerre α-ensemble [48]. The connection is obtained noticing that the matrix L admits the following decomposition where B is the matrix transpose. On the other hand, the Laguerre α-ensemble is given by the set of matrices The variables x n , y n are distributed according to chi-distribution y n ∼ χ 2α n = 1, . . . , N − 1.
Thus, the following entry wise measure on the matrices B α,γ can be defined, We observe that the matrix B in (18) has the same form of B α,γ in (20), with the addition of the corner element y N E 1,N . Furthermore, the rescaling of the variables (x j , y j ) → 1 √ 2β (x j , y j ) in (17), amounts to the matrix rescaling B → 1 √ 2β B, and comparing with (21) we see that 1 √ 2β B is a rank one perturbation of the matrix B θ, θ η . We are interested in studying the density of states ν ET for the Lax matrix L when the entries are distributed according to the Gibbs measure dµ ET in (17). The density of states ν ET is obtained from the weak convergence of the empirical measure of the Lax matrix L, namely where λ j are the eigenvalues of L and δ x is the Dirac delta function centred at x. In order to study the density of states of the exponential Toda lattice, we need the following result proved in [48]. (21). Then, its mean density of states ν Lα,γ takes the form and here ψ(v, w; z) is the Tricomi's confluent hypergeometric function (see Appendix B).
Remark 3.2. The proof of theorem 3.1 has been obtained by comparing the Laguerre α-ensemble (19) with the Laguerre β-ensemble at high temperature [7].
The following corollary follows.
Proof. First, we notice that by virtue of general theory, see [12,Theorem A.43], we can restrict to the case y N = 0 in (17). As observed above, performing the change of variables (x j , y j ) → 1 √ 2β (x j , y j ), which amounts to rescale B → 1 √ 2β B, one has that the matrix entries of 1 √ 2β B are distributed as the matrix entries of B θ, θ η . Applying Theorem 3.1 we obtain the claim.

Parameter Limit
In this section, we examine the low-temperature limit of the Hamiltonian system (11). Namely, we want to compute the eigenvalues of the Lax matrix L in (15) in the limit β, θ, η → ∞, in such a way that η = ηβ, θ = θβ, where η and θ are in compact sets of R + .
Since all x j and y j are independent random variables, we just have to consider the weak limit of the rescaled chi-distributions, respectively We explicitly work out one of the cases above.
We consider a continuous and bounded function h : R + → R and evaluate the limit The last identity has been obtained by applying the Laplace method (see [50]) and observing that the minimizer of the term 2 η log(x) − x 2 in the exponent of the integral is x 0 = η.
As a consequence, we conclude that x j η and y j θ, j = 1, . . . , N as β → ∞, where with we denote the weak convergence. The previous limit implies that the measure ν ET in (24) weakly converges, in the low temperature limit, to the density of states of the matrix L ∞ Indeed, the fact that L is tridiagonal with iid entries along the diagonals, implies its k-th moment depends on a multiple of k number of variables only; specifically looking back at the Lax matrix L in (16) for some function f (·) of its entries. Then, passing to the density of states (22) and renaming the iid variables, the scaling factor N identically cancels out and moments converge, The eigenvalues being functions of moments, density of states converges as well.
In particular, this also shows that the two limit commute in taking the density of states at low temperature, since the limits can be passed directly to the variables x i , y i . The matrix L ∞ is a circulant matrix, so its eigenvalues can be computed explicitly [33] as From this explicit expression, it follows that the density of states of L ∞ is here 1 (c − ,c + ) is the indicator function of the set (c − , c + ). In particular, this measure belongs to the class of Arcsin distributions. Thus, we proved the following result.
Proposition 3.4. Consider the random Lax matrix L in (15) sampled from the Gibbs ensemble dµ ET (17) of the exponential Toda lattice (14). The density of states of the matrix L in the lowtemperature limit, i.e. when β, θ, η → ∞ in such a way that η = ηβ, θ = θβ, with η, θ in compact subsets of R + , is the Arcsin distribution given by (25).

Volterra lattice
The Volterra lattice, also known as the discrete KdV equation, describes the motion of N particles on the line with equationsȧ It was originally introduced by Volterra to study population evolution in a hierarchical system of competing species. It was first solved by Kac and van Moerbeke in [41] using a discrete version of inverse scattering due to Flaschka [26]. Equations (26) can be considered as a finite-dimensional approximation of the Korteweg-de Vries equation. The phase space is R N + and we consider periodic boundary conditions a j = a j+N for all j ∈ Z. The Volterra lattice is a reduction of the second flow of the Toda lattice [41]. Indeed, the latter is described by the dynamical systeṁ and equations (26) are recovered just by setting b j ≡ 0. The Hamiltonian structure of the equations follows from the one of the Toda lattice. On the phase space R N + we introduce the Poisson bracket and the Hamiltonian H 1 = N j=1 a j so that the equations of motion (26) can be written in the Hamiltonian formȧ An elementary constant of motion for the system is H 0 = N j=1 a j that is independent of H 1 . The Volterra lattice is a completely integrable system, and it admits several equivalent Lax representations, see e.g. [41,51]. The classical one readṡ where we recall that the matrix E r,s is defined as (E r,s ) ij = δ i r δ j s and E j+N,i = E j,i+N = E j,i . There exists also a symmetric formulation due to Moser [51], which assumes that all a j > 0. Furthermore, we point out that there exists also an antisymmetric formulation for this Lax pair, indeed a straightforward computation yields

Gibbs Ensemble
We introduce a Gibbs ensemble for the Volterra lattice (26) by observing that its vector field f j = a j (a j+1 − a j−1 ) is divergence free, due to the periodic boundary conditions. Therefore, an invariant measure can be obtained from (10). We use H 0 = N j=1 a j , and H 1 = N j=1 a j as constants of motion to construct the invariant measure where and Γ (η) is the Euler gamma function [19, §5]. We notice that according to this measure, all the variables are independent and identically distributed (i.i.d.).
Next we want to characterize the density of states of the antisymmetric Lax L 3 of the Volterra lattice given in Proposition 4.1. Among the three Lax matrices of the Volterra lattice, the matrix L 3 is particularly useful since it allows us to connect the Volterra lattice with a specific α-ensemble, namely the antisymmetric Gaussian α-ensemble. The antisymmetric Gaussian α-ensemble, see [28], is the family of random antisymmetric tridiagonal matrices where y i are i.i.d. random variables with density which is just a rescaled chi-distribution. Even though we use a different expression of the chidistribution with respect to section 3, we keep the same notation f 2α (y) for the density that will be used only in this section. This distribution induces a measure on the entries of the matrix L α , namely In [28] the authors studied this matrix ensemble in connection with the antisymmetric Gaussian βensemble introduced by Dumitriu and Forrester [22] in the high temperature regime, and computed explicitly its density of states ν Lα (x), defined as where λ j are the eigenvalues of L α . Since the matrix L α is antisymmetric with real entries, its eigenvalues are purely imaginary numbers.
Theorem 4.2 (cf. [28]). The density of states of the random matrix L α in (32), is explicitly given by here Γ(x) is the gamma function and W k,µ (z) is the Whittaker function (see Appendix B) .
Remark 4.3. The proof of Theorem 4.2 has been obtained in [28] by comparing the antisymmetric Gaussian α-ensemble (32) with the antisymmetric Gaussian β-ensemble at high temperature, which was considered in the same paper.
We notice that performing the change of coordinates a j = x 2 j , the Gibbs ensemble (30) reads: which, up to a rescaling x j → x j / √ β and for the extra term x N in the probability distribution, is exactly the distribution (33) of the matrix L α . Furthermore, the matrix L 3 is a 2 rank perturbation of the matrix L α . Therefore, by a corollary of [12,Theorem A.41] and Theorem 4.2, we obtain the following.  (30). Then, the density of states of the matrix L 3 is explicitly given by where θ α (x) is given in (34).

Parameter Limit
As for the case of exponential Toda (14), in this section we consider the low-temperature regime of the Volterra lattice, namely the limit η, β → ∞, in such a way that η = β η, with η in a compact set of R + , and we compute the density of states of the matrix L 3 (28) in this regime.
Applying the same techniques of Section 3.1, we conclude that the density of states of the matrix L 3 in the low-temperature limit coincides with the one of the matrix L ∞ , where Since the matrix L ∞ is circulant, we can readily compute its eigenvalues as From this explicit formula, it follows that the density of states of the matrix L ∞ reads Such measure coincides with the measure ν L α in the low-temperature limit, and it belongs to the class of Arcsin distributions. Thus, we just proved the following.

Generalization of the Volterra lattice: the INB k-lattices
The Volterra lattice (26) can be generalized in a variety of ways. The most natural ones are two families of lattices described in [17] (see also [16,40,52] ) which include short range interactions. The first family is called additive Itoh-Narita-Bogoyavleskii (INB) k-lattice and is defined by the equationsȧ The second family is called the multiplicative Itoh-Narita-Bogoyavleskii (INB) k-lattice and is defined by the equationṡ In both cases we consider the periodicity condition a j+N = a j holds. Setting k = 1, we recover from both lattices the Volterra one (26). Further generalizations of the INB lattice were recently considered in [25]. A crucial difference in the two models is that in the additive lattice (35) the interaction is on arbitrary number of points, but the non-linearity is still quadratic like the original Volterra lattice (26); on the other hand , the multiplicative lattice (36) admits non-linearity of arbitrary order. Moreover, both families admit the KdV equation as continuum limits, see [17].
As mentioned earlier, the additive INB k-lattice is an integrable system for all k ∈ N and i ∈ Z, since they all admit a Lax pair formulation (2). For the additive INB lattice (35), it reads we recall that we are always considering periodic boundary conditions, so for all j ∈ Z, a j+N = a j and E i,j+N = E i+N,j = E i,j . The constants of motion obtained through this Lax pair are in involution with respect to the Poisson bracket Then, the additive INB k-lattice (35) can be written aṡ where the Hamiltonian function H 1 = N j=1 a j is the same as in equation (27). In the same way, it is possible to prove that the function H 0 = N j=1 a j is a first integral for the additive INB k-lattice (35) as well.
Similarly, the multiplicative INB k-lattices can be endowed with a Lax Pair for all k ∈ N, therefore it is another example of integrable systems. Specifically, for the periodic case we presented in equation (36), the Lax pair reads We notice that both H 1 = N j=1 a j , and H 0 = N j=1 a j are constants of motion for these systems, for all k ∈ N.
Remark 5.1. For fixed k, there exists a transformation that maps the multiplicative INB k-lattice to the additive one. Namely, consider the system (36) and define the new set of variables b i := a i · · · · · a i+k−1 , i = 1, . . . , N, where the indices are taken modulo N . Then, it is immediate to see thaṫ which in turn, due to telescopic summations, implieṡ which is (35). The transformation (39) is invertible only when k and N are co-prime, for a more detailed discussion see [17].

Gibbs Ensemble
We want to introduce an invariant measure for the INB lattices (35) and (36). Since H 0 = N j=1 a j and H 1 = N j=1 a j are constants of motion for all the INB lattices, and both systems are divergence free, in view of (10) we can consider as invariant measure the same one that we used for the Volterra lattice, namely where the normalization constant Z INB N (β, η) has the value given in equation (31). For example the random Lax matrix L (+,k) takes the form where the χ 2 2η distribution has density a 2η−1 2 η Γ(η) e − a 2 . Unlike the Lax matrix of the Volterra lattice, the Lax matrices of these generalizations lack of a known random matrix model to compare with. For this reason, we present numerical investigations of the density of states for these random Lax matrices for several values of the parameters k, η and β, see Figures 2-3. We notice that, for both the additive lattice and the multiplicative one, the density of states seems to possess a discrete rotational symmetry. In this spirit, we prove the following  Proof. We prove the statement for the additive case, the proof in the multiplicative one is analogous.
We notice that, due to the structure of L (+,k) , if L (+,k) (1, i 1 ) · · · L (+,k) (i −1 , 1) is not zero, then either i s+1 = i s + 1 or i s+1 = i s − k modulo N . Now, consider paths in the Z 2 plane from the point (0, 0) to ( , 0), such that the only permitted steps are the up step (1, 1) and the down step (1, −k). Since these paths resemble the classical Dyck paths, we call them (1, k)-Dyck paths of length . Given a non-zero element of the product in (41), we can construct the corresponding path in the following way. We start at (0, 0), then if |i 1 − 1| = 1 we make an up step of height 1, otherwise we make a down step of height k, and so on. For each path, let n be the number of up steps and m the number of down steps, then m + n = , n − mk = 0 , since there is a total of step, and the path has to go back to height 0. Thus, we deduce that m(k + 1) = , and the claim is proven.  Another interesting feature of these measures is that their supports seem to be exponentially localized to one dimensional contours. Specifically, it appears that the supports are the two hypotrochoids γ +,k , γ ×,k , respectively This feature is highlighted in figures 2-3, where we plot the empirical density of states and the corresponding hypotrochoid. This characteristic is important since this type of curves are also related to the density of some cyclic digraph, see [5], and may serve as a link between these two topics.
All these observations lead us to formulate the following conjecture Moreover, the densities are exponentially localized in a neighbourhood of the two hypotrochoids γ +,k (t, η, β) and γ ×,k (t, η, β) in (42) respectively.

Parameter limit
As in the previous cases, although we are not able to give an explicit formula for the density of states of the INB lattices for general β, η, we can characterize this measure in the low-temperature limit. Specifically, we consider the limit as β, η → ∞ in such a way that η = ηβ, with η in a compact set of R + , and we compute the density of states of the matrices L (+,k) , L (×,k) , endowed with the probability inherited from dµ INB (40), in this limit.
The procedure is the same as in the case of Volterra (see Section 4.2). Indeed, following the same line, we can conclude that the densities of states ν ∞ INB,+,k and ν ∞ INB,×,k coincide with the densities of L (+) and L (×) respectively, where We notice that both matrices are circulant, thus we can compute their eigenvalues explicitly as here j = 1, . . . , N . Thus, in the large N limit, we deduce that the support of the measures ν ∞

INB+,k
and ν ∞ INB×,k are the hypotrochoids and the limiting eigenvalue densities are We summarize these results in the following Proposition.
Proposition 5.4. The densities of states of the Lax matrices L (+,k) and L (×,k) in (37) and (38) endowed with the Gibbs measure dµ INB in (40), in the low temperature limit, i.e. when η, β → ∞, in such a way that η = ηβ, with η in a compact set of R + , are given respectively by ν ∞ INB+,k and ν ∞ INB×,k in (43).

The focusing Ablowitz-Ladik lattice
The focusing Ablowitz-Ladik lattice is the following system of spatial discrete differential equations where a j ∈ C, j = 1, . . . , N , N ≥ 3, and we consider periodic boundary conditions a j+N = a j for all j ∈ Z. This equation was introduced by Ablowitz and Ladik [1,2], by searching integrable spatial discretization of the cubic non-linear Schrödinger Equation (NLS) for the complex function ψ(x, t), x ∈ R, t ∈ R + i∂ t ψ(x, t) + ∂ 2 x ψ(x, t) + 2|ψ(x, t)| 2 ψ(x, t) = 0 , In contrast with what happens in the defocusing case, the particles (a 1 , . . . , a N ) are free to explore the whole C N , which is the phase space of the system.
On the space C ∞ (C N ) we consider the Poisson bracket [24,29] {f, We notice that the phase shift a j (t) → e −2it a j (t) transforms the AL lattice (44) into the equatioṅ which we call the reduced AL equation. We remark that the quantity H 0 = 2 ln N j=1 ρ 2 j is the generator of the shift a j (t) → e −2it a j (t), while generates the flow (45). Therefore, we can rewrite the AL equation aṡ Moreover, it is straightforward to verify that {H 0 , H 1 } = 0. The Poisson bracket induces the symplectic form that is invariant under the evolution generated by the Hamiltonians H 0 and H 1 . Therefore, the volume form ω N = ω ∧ · · · ∧ ω, is also invariant. In view of these properties, we can define the Gibbs ensemble for the focusing Ablowitz-Ladik lattice on the phase space C N as where a = (a 1 , . . . , a N ), d 2 a = N j=1 (ida j ∧ da j ) and Z AL N (β) is the normalization constant of the system. We notice that according to this measure, all the variables are i.i.d.
Remark 6.1. The measure with density exp(−βH AL ) and β > 0 is not bounded nor normalizable on the whole phase space. For this reason, we have defined the Gibbs ensemble as in (47). Furthermore, we observe that the measure (47) has a finite number of moments, which implies that the corresponding density of states of the Lax matrix (see (49) below), if it exists, would have a finite number of moments.
The focusing AL lattice is a complete integrable system. Indeed it admits a Lax representation, first obtained by Ablowitz and Ladik from the discretization of the Zakharov-Shabat Lax pair for the focusing non-linear Schrödinger equation [63]. Gesztesy, Holden, Michor, and Teschl [30] found a different Lax pair for the infinite case of focusing AL lattice, and for its general hierarchy. To adapt their construction, we double the size of the lattice according to the periodic boundary conditions, thus we consider a chain of 2N particles a 1 , . . . , a 2N such that a j = a j+N for j = 1, . . . , N . Define and the 2N × 2N matrices

Now let us define the Lax matrix
that has the structure of a 5-band diagonal matrix  * * * * * * * * * * * * * * * * . . . . . . * * * * * * * * * * * * * * * * The N -periodic equation (45) is equivalent to the following Lax equation for the matrix E, where the two projections M + , M − are defined for a 2N × 2N matrix as We notice that the Lax matrix E has a similar structure to the one of the defocusing AL lattice obtained by Nenciu, and Simon [53,55]. The crucial difference is that while for the defocusing AL lattice the blocks Ξ j are unitary matrices, for the focusing lattice this is not the case since Ξ j Ξ † j = I 2 where I 2 = 1 0 0 1 and † stands for hermitian conjugate.
The measure dµ AL induces a probability distribution on the entries of the matrix E, thus it becomes a random matrix. As in the previous cases, one would like to connect the density of states for this random matrix to the density of states of some β-ensemble in the high temperature regime, but, as in the case of the INB lattices, we lack of a matrix representation of some β-ensemble with eigenvalues supported on the plane.
We make the following observations. The matrix Ξ j (48) is complex symmetric, and it can be factorized in the form we can rewrite the Lax matrix E (49) as where, defining c j := a j |a j | |a j |+ √ 1+|a j | 2 , the matrices Λ odd and Λ even are given by Since we are interested in the distribution of the eigenvalues of E, it follows from the factorization (51) that we can also consider the matrix Λ odd EΛ even E , where E = L M. The eigenvalues of Λ even , Λ odd come in pairs, such that if λ is an eigenvalue, then also −λ −1 is an eigenvalue. The matrix E is a periodic CMV matrix [18], thus its eigenvalues are on the unit circle. Thus, we are in a similar setting considered in [34,59,60]. Indeed in [59] the authors derived the eigenvalues distribution of U √ D where U is a Haar distributed unitary matrix and D is a fixed diagonal matrix with positive eigenvalues. They show that the density of states is rotational invariant and it is supported on a single ring whose radii r 1 < r 2 satisfy the constraint d min < r 1 < r 2 < d max , where d min and d max are the minimum and maximum eigenvalues of D. In [34], the authors considered a similar problem, namely the characterization of the density of states for a matrix of the form U T V , where U, V are independent unitary matrices Haar distributed, and T is a real diagonal matrix independent of U, V . They proved, under some mild conditions, that the density of states of the matrix U T V is radially symmetric and it is supported on a ring.
It is therefore reasonable to expect that the density of states of the random Lax matrix of the Ablowitz-Ladik lattice is rotational invariant, but with unbounded support, indeed the eigenvalues of Λ even , Λ odd could grow to infinity.
To confirm our expectations, we perform several numerical investigations of the eigenvalues of the random Lax matrix of the Ablowitz-Ladik lattice for various values of β (see Fig 5-6). In Fig  5 the eigenvalue density is shown. As expected, the density seems to be rotational invariant, and concentrated on a ring, exactly as in [34,59,60]. For this reason, we investigate the behaviour of the modulus of the eigenvalues, see Fig 6. They seem to be concentrated in a small region, but, in view of Remark 6.1, we expect that the tails should decay just polynomially fast.

Parameter Limit
Despite not being able to explicitly compute the density of states for general values of β, we can perform such an analysis in the low-temperature limit, namely when β → ∞.
We notice that, according to (47), all the a j are independent. Hence, in order to obtain the density of states in the low-temperature limit, we have to compute the weak limit of the density Proceeding as in the previous cases, it follows that the following limit holds for all bounded and continuous f : D → R: The previous limit implies that the density of states of the Ablowitz-Ladik lattice in the low temperature limit is equal to the one of E = L M, where L, M are 2N × 2N matrices and Ξ is defined as the unitary matrix To compute the density of states for the matrix E, we notice that both M, and L are permutation matrices. Specifically, we identify them with the following permutations   This implies that the eigenvalues of E are all double, and given by From this explicit expression of the eigenvalues, it is straightforward to prove that Thus, we have proved the following Proposition 6.2. Consider the Gibbs ensemble dµ AL (47) of the focusing Ablowitz-Ladik lattice in the low-temperature limit, i.e. β → ∞. Then, the density of states ν AL of the Lax matrix E (49) is given by

Schur flow
The focusing Schur flow, also known as discrete mKdV, is an integrable system deeply related to the Ablowitz-Ladik lattice. Its equations of motion arė where here we consider periodic boundary conditions, a j = a j+N for all j ∈ Z. Notice that if a j (0) ∈ R for all j = 1, . . . , N , then a j (t) ∈ R for all times, implying that R N is an invariant subspace for the dynamics.
Recalling the definition (46) for K (1) = N j=1 a j a j+1 and introducing the Poisson bracket we can rewrite the equations of motion (52) aṡ Notice that the quantity H 0 = −2 ln N j=1 ρ 2 j is a global first integral for the system. Moreover, one can deduce immediately from the equations of motion that R N is invariant for the dynamics. Thus, in view of the Hamiltonian representation and this invariance, we define the Gibbs measure for the Schur flow as where Z S N (β) is the normalization constant of the system, Remark 7.1. Similarly to the focusing AL case, the classical Gibbs ensemble is not well-defined on the whole phase space. Indeed, the measure with density e −βH S , β > 0 cannot be normalized on R N .
The Schur flow is a completely integrable system since it admits a Lax formulation. Namely, define the 2 × 2 matrix Ξ j and, similarly to the Ablowitz-Ladik case, the 2N × 2N matrices Then, the N -periodic equation (52) is equivalent to the following Lax equation for the matrix E: where the two projection +, − are defined in (50).
Carrying on with the approach of this article, we study the density of states ν S for the matrix E when the entries are distributed according to the measure (53). First, we notice that Remark 6.1 is valid also in the case of the focusing Schur flow. Moreover, since the variables a j are real, we can factorize the matrices Ξ j in the following way: where we note that U −1 0 = U 0 , so that the above is a similarity transformation. Thus, defining we can rewrite the Lax matrix of the Schur flow E as where Λ odd = diag c 1 , − 1 c 1 , c 3 , . . . , , c 2 , . . . , c 2N .
As in the case of the Ablowitz-Ladik lattice, since we are interested in just the distribution of the eigenvalues of E, we can consider the matrix Λ odd EΛ even E , where As in the AL case, the eigenvalues of Λ even , Λ odd come in pairs, such that if λ is an eigenvalue, then also −λ −1 is an eigenvalue. The main difference with the case of the focusing AL lattice is that in this case the matrix E is deterministic. Thus, one can be led to think that the eigenvalue distribution of the Schur flow would be similar to the one of the AL lattice, but it is not the case. Indeed, we perform several numerical investigations, reported in Fig 7, which shows that the behaviour of the eigenvalues is different in the two situations. We notice that a consistent part of the eigenvalues tend to stay close to the real axis, see Fig  7. This behaviour is also typical of the orthogonal Ginibre ensemble [23]. The main reason is that the eigenvalues of E are not evenly spaced on the unit circle, but they are constrained to the left semicircle, and are more dense nearby ±i (see Fig 8). Indeed we can give an accurate description of the spectrum of this matrix.
More precisely, we have the following.
counting the multiplicity.
Passing to the limit N → ∞, by standard methods, we can compute the limiting density of the argument ϕ of the eigenvalues as This behaviour is shown in Fig 8. Furthermore, as in the case of the AL lattice, for large β the eigenvalues tend to accumulate on the unit circle, see Fig 7. In a similar way as done for the focusing AL lattice we conclude that the density of states of the random Lax matrix E in (49) with probability distribution entries given by the Gibbs measure dµ S in (53), converge in the limit β → ∞ to the uniform measure on the unit circle.

Conclusions and outlook
In this manuscript, we underlined the deep relation between integrable systems and random matrix theory. Specifically, we showed that the exponential Toda lattice, and the Volterra one are related to the Laguerre β-ensemble at high temperature and the antisymmetric Gaussian β-ensemble at high temperature respectively. As we already mentioned in the introduction, these are not isolated and lucky relations, indeed the Toda lattice is related to the Gaussian β-ensemble at high temperature, and the Ablowitz-Ladik one to the Circular β ensemble. Thus, in view of the new relations that we obtained, we conclude that each classical β-ensemble in the high temperature regime is related to some integrable model, see Table 1.
Furthermore, we numerically investigate the additive and multiplicative INB lattices, which are a generalization of the Volterra one. Since the former lattice is related to the antisymmetric Gaussian β-ensemble, it may be possible that the INB lattice could be related to some generalization of the antisymmetric Gaussian β-ensemble themselves.
Finally, we considered the focusing Ablowitz-Ladik lattice and the focusing Schur flow. In view of our numerical investigations, we were able to formulate a precise conjecture regarding the eigenvalues distribution of the focusing Ablowitz-Ladik lattice. Considering the focusing Schur flow, we were not able to formulate a precise conjecture for its eigenvalues distribution. The expectation is that it would be related to the a possible Ginibre β-ensemble.
It would be fascinating to consider the generalized Gibbs ensemble for the exponential Toda lattice and the Volterra one, and prove that both their empirical measures satisfy a Large Deviation principle in the spirit of [35,49]. Moreover, it would be interesting to apply the theory of Generalized Hydrodynamic [13] to have some insight regarding the behaviour of the correlation functions for these systems. Recall that the matrix E is defined as E = L M where L and M are as in (54). It is a block circulant matrix, indeed we can write One can immediately check that λ = ±i are eigenvalues for E with eigenvectors v ±i = (∓i, 1, ∓i, 1, . . . , ∓i, 1) .
We now claim that, for fixed N , the remaining eigenvalues have multiplicity 2 and are the (N − 1) solutions to Ω (λ) N = I 2 , where we defined Such solutions are obtained by solving the equation For j = 0 the solutions to (57) are ±i which we already treated separately. Indeed Ω (±i) is not diagonalizable and Ω (±i) N = I 2 for every N greater than 0. For any other j ∈ {1, . . . , N 2 }, the solutions to (57) Since both the real and imaginary part are monotone functions of j, different j s will correspond to different solutions. Hence, if N is odd, we will have a total of N − 1 solutions coming from (58); if N is even one has N − 2 distinct solutions coming from the equation in (57)  Let us check the correctness of the claim. Write Ev 1 := (w 1 , . . . , w N ) , where w j are two-dimensional row vectors, then using the fact that Ω (λ) N = I 2 , one can compute for any k = 1, . . . , N , which shows that v 1 is an eigenvector with eigenvalue λ. The same proof clearly applies to the other eigenvector v 2 .

B.1 Tricomi Confluent Hypergeometric function
The Tricomi Confluent Hypergeometric function ψ(a, b; z), that we introduced in Theorem 3.1, is defined as the solution of the Kummer's equation such that ψ(a, b; z) ∼ z −a , z → ∞ , | arg(z)| ≤ 3 2 π − δ , here δ is a small positive constant. The origin is a regular singular point for this equation with indexes 0 and 1 − b, moreover, this equation has an irregular singularity at infinity of rank one. Usually, ψ(a, b; z) has a branch point at 0, whose principal branch is the same as z −a , thus it has a cut on the z plane along the interval (−∞, 0].