Kinetic theory of quantum and classical Toda lattices
Vir B. Bulchandani Xiangyu Cao Joel E. Moore



^1 Department of Physics, University of California, Berkeley, Berkeley CA 94720, USA

^2 Materials Science Division, Lawrence Berkeley National Laboratory, Berkeley CA 94720, USA

 We consider the kinetic theory of the quantum and classical Toda lattice models. A kinetic equation of Bethe-Boltzmann type is derived for the distribution function of conserved quasiparticles. Near the classical limit, we show that the kinetic theory depends smoothly on Planck's constant, and explicitly characterise the leading quantum corrections to classical behaviour. The classical kinetic theory is compared with direct numerical simulations and shows excellent agreement. Finally, we connect the Bethe-Boltzmann approach with the classical inverse scattering method, by identifying conserved quasiparticles with the spectrum of the Lax matrix.





§ INTRODUCTION
  

The derivation of hydrodynamics from microscopic many-particle dynamics is an enduringly difficult mathematical problem. Despite the technical hurdles involved, rigorous progress has been made for some simple classical many-body systems <cit.>.  However, beyond the relatively tractable case of weakly interacting Bose gases <cit.>, the hydrodynamics of quantum many-body systems remains poorly understood.

Recently, a promising new avenue of inquiry has opened up with the development of “generalized hydrodynamics” of quantum integrable systems in one dimension <cit.>, which has been accompanied by rigorous results on the local equilibrium dynamics of one-dimensional conformal field theories <cit.>. The Bethe-Boltzmann kinetic equation for quasiparticle densities, which captures some essential features of generalized hydrodynamics with an infinite set of conserved quantities and Euler equations, has yet to be derived rigorously, but there is a widespread belief in the field, supported by a wealth of numerical and analytical tests <cit.> (and a recent experiment <cit.>), that this equation describes the ballistic part of local equilibrium dynamics in such models exactly.

The kinetic equation can be understood as a generalization of earlier kinetic theories of classical soliton gases <cit.>, with classical phase-shifts replaced by the quantum mechanical phase-shifts between Bethe quasiparticles.  The concept that the dynamics of Bethe quasiparticles could capture far-from-equilibrium transport in quantum integrable spin chains was discussed early on <cit.>.  Given that the kinetic approach amounts to describing quantum dynamics by an essentially classical equation, a natural question that arises is the following: how does a kinetic theory of quantum solitons differ from that of classical solitons? The integrable Toda lattice, whose thermodynamic Bethe Ansatz (TBA) is understood both classically and quantum mechanically <cit.>, provides an ideal setting in which to analyse this important question.

In the present work we derive Bethe-Boltzmann kinetic equations for the quantum and classical Toda lattices (Section <ref>). Hydrodynamic predictions obtained for the classical Toda lattice are compared against microscopic numerical simulations of “two-reservoir” initial conditions  and show excellent quantitative agreement (Section <ref>). Our kinetic equation for the classical Toda lattice can be understood as the semiclassical limit of the kinetic theory of the quantum system, allowing for a direct, quantitative analysis of the effect of introducing a small and non-zero ħ in the Bethe-Boltzmann equation.  Hence the Toda problem illustrates, in a single model with variable quantum parameter ħ, how the same type of kinetic equation describes the hydrodynamic limit of quantum and classical integrable systems.

Finally, we make a connection between equilibrium quasiparticle distribution functions of the classical Toda lattice and classical inverse scattering theory, interpreting the local distribution functions of quasiparticle density and current directly in terms of the spectrum of the Lax matrix. This “Bethe-Lax correspondence” turns predictions of the kinetic theory to mathematical statements about tridiagonal random matrix ensembles. They are convincingly supported by our numerical results, which indicate that the conjectured equations of state for the quasiparticle density and current distribution functions in Section <ref> hold exactly for equilibrium states. The Bethe-Lax correspondence offers a new qualitative insight into the origin of this powerful hydrodynamic formalism. In the classical Toda case, an intriguing outcome would be a rigorous derivation of our equation of state for the quasiparticle currents from Dumitriu-Edelman β-random matrix theory <cit.>, as was very recently achieved for the quasiparticle density of states <cit.>.
We note that the hydrodynamics of classical Toda lattice is treated in an extensive recent work <cit.> using a different approach.



§ KINETIC EQUATIONS: QUANTUM AND CLASSICAL

  



 §.§ Quantum TBA and kinetic theory

The quantum Toda Hamiltonian can be defined as

    H = ∑_n -ħ^2/2m∂_n^2 + V(x_n+1-x_n)  ,  
    V(r) = mω^2/γ^2 e^-γ r  .

Here γ quantifies the anharmonicity of the potential and ω is the harmonic frequency as γ→ 0. Using Sutherland's technique of “asymptotic Bethe ansatz”, one can deduce the equations of thermodynamic Bethe ansatz <cit.> for this model. At temperature T and chemical potential μ, the Yang-Yang equation for excitation energies is given by

    ϵ_k = ħ^2k^2/2m - μ + 1/β∫_-∞^∞ dk'   K_k,k'ln(1+e^-βϵ_k')  ,

while the total density of states ρ^t_k satisfies

    ρ_k^t = 1/2π -∫_-∞^∞ dk'   K_k,k'θ_k'ρ^t_k'  ,

where θ_k = (1+e^βϵ_k)^-1 denotes the  Fermi factor. The differential phase-shift is given explicitly by

    K_k,k' = 1/2πγ[2 lnC - ψ(1-i(k-k')/γ) - ψ(1+i(k-k')/γ)]  .

Here, the dimensionless parameter C=mω/ħγ^2 and ψ(x) denotes the digamma function. Abstracting to general local equilibrium states, with arbitrary Fermi factors θ_k : ℝ→ [0,1], we find that the derivatives of quasiparticle energies and momenta satisfy

    ϵ'_k +∫_-∞^∞ dk'   K_k,k'θ_k'ϵ'_k'   = ħ^2k/m 
    p^dr'_k +∫_-∞^∞ dk'   K_k,k'θ_k'p^dr'_k'   = ħ  .


From these formulae, one can write down a “Bethe-Boltzmann” equation, modelling the ballistic part of locally equilibrated time-evolution in the quantum Toda lattice. The resulting formalism is identical to that for the Lieb-Liniger gas <cit.>, and the initial value problem for smoothly varying temperature and chemical potential profiles β(x) and μ(x) can be conveniently summarized in terms of local Fermi factors θ_k(x,t), which evolve under the system of equations

    ∂_t θ_k + v_k[θ]∂_x θ_k = 0,    v_k[θ] = (1̂ + K̂θ̂)^-1[ħ k'/m]_k/(1̂ + K̂θ̂)^-1[1]_k  ,

with initial data determined locally in terms of {β(x),μ(x)} via thermodynamic Bethe Ansatz. Here, the notation Ô indicates the integral operator corresponding to a given integral kernel O_k,k'. In particular, (1̂+K̂θ̂)^-1 stands for the inverse of the integral operator with kernel  δ(k-k')+K_k,k'θ_k', that appears in Eq.  <ref>. Thus (1̂ + K̂θ̂)^-1[f]_k ∫ dk'   U_k,k'f_k', with U_k,k' determined by the integral equation

    U_k,k' + ∫ dk” K_k,k”θ_k”U_k”,k' = δ(k-k').




 §.§ The semiclassical limit of TBA


We now turn to the semiclassical limit ħ→ 0 of the quantum TBA and hydrodynamics outlined above, following closely the analysis of Ref. <cit.>. We proceed in two stages: we first express the quantum TBA equations in terms of classical momentum p = ħ k, before taking the semiclassical limit as ħ→ 0. In order to guarantee a finite free energy in the semiclassical limit <cit.>, we define classical quasiparticle energies and classical chemical potentials, as

    ϵ^cl_p = ϵ_p + 1/βlnβħω , μ^cl = μ- 1/βlnβħω  .

The Yang-Yang equation for ϵ^cl_p reads

    ϵ^cl_p = p^2/2m - μ^cl + 1/βħ∫_-∞^∞ dp'   K_p/ħ,p'/ħln(1+βħω e^-βϵ^cl_p')  ,

the Fermi factor is given by

    θ_p = βħω e^-β^cl_p/1+βħω e^-βϵ^cl_p  ,

and the dressing equations read

    ϵ^cl'_p +∫_-∞^∞ dp'   K_p/ħ,p'/ħ(θ_p'/ħ)ϵ^cl'_p'   = p/m   , 
    p^dr'_p +∫_-∞^∞ dp' K_p/ħ,p'/ħ(θ_p'/ħ)p^dr'_p'   = 1   .

The rescaled density of states, defined by ρ_p^t dp = ρ^t_k dk, satisfies 2πħρ^t_p = p^dr'_p. Note that equations Eq. <ref> through <ref> are merely a rewriting of quantum TBA, and are exact for ħ>0. Later, they will allow us to quantify the quantum corrections to classical hydrodynamics. Let us first consider passing to the semiclassical limit ħ→ 0 of quantum TBA. 

In this limit, the phase-shift is replaced by its leading asymptotic,

    K^cl_p,p'   = lim_ħ→ 0 K_p/ħ,p'/ħ = -1/πγ(ln|γ/mω| + ln|p-p'|)  ,

the Yang-Yang equation becomes

    ϵ^cl_p    = p^2/2m - μ + ∫_-∞^∞ dp'   K^cl_p,p'ω e^-βϵ^cl_p'  ,

and the effective Fermi factor appearing in the dressing equations is given by

    θ_p = lim_ħ→ 0θ_p/ħ = βω e^-βϵ^cl_p  .

In particular, the density of occupied states reads

    ρ_p = lim_ħ→ 0ρ^t_p θ_p = θ_p(1+K̂_clθ̂)^-1(1/2π)_p   ,

and from this expression, the thermal expectation values of conserved charges of the classical Toda lattice, ⟨ Q_n ⟩_ρ = ∫ dp  ρ_p p^n, may be evaluated. The relation Eq. <ref> is verified numerically in Section <ref>, by direct comparison with spectrum of the Lax matrix.



 §.§ Classical kinetic theory


For any non-zero ħ, the quasiparticle distribution function of the full quantum TBA includes both “phonon” (particle-like) and “soliton” (hole-like) degrees of freedom. However, in the semiclassical limit at non-zero temperature, the soliton degrees of freedom become infinitely heavy and do not contribute to the thermodynamics of the model <cit.>. Thus for the purposes of considering local equilibrium dynamics in the Toda lattice, it suffices to restrict attention to phononic degrees of freedom; this will be justified more carefully below. The phonon dispersion relation is given by <cit.>

    Δ E^ph_p    = ħω e^-βϵ^cl_p,   Δ P^ph_p = 2πħ∫_p^∞ dp'  ρ_p',

and implies an interaction-dressed phonon group velocity

    v_p = ΔE^ph'_p/ΔP^ph'_p = ϵ^cl_p'/2πħρ_p^t,

for a state with occupied density of states ρ_p. From this expression, we can write down a Bethe-Boltzmann equation for the local phonon density of states, given by

    ∂_t ρ_p(x,t) + ∂_x (ρ_p(x,t) v_p[ρ](x,t)) = 0.

As usual, this may be expressed in advection form as the system

    ∂_t θ_p + v_p[θ]∂_x θ_p    = 0,   v_p[θ]    = (1̂+K̂^clθ̂)^-1[p'/m]_p/(1̂+K̂^clθ̂)^-1[1]_p,

with initial data determined in terms of initial temperature and (classical) chemical potential via the classical TBA equations Eqs. <ref> and <ref>. Solutions to this system can be compared directly with classical numerical simulations of Toda lattice dynamics. Such a comparison is provided in Section <ref>.

While the above discussion was physically motivated, notice that the classical hydrodynamic equation Eq. <ref> could have been obtained directly from the quantum hydrodynamic equations Eq. <ref> by taking the semiclassical limit described above (the possibility of obtaining classical hydrodynamics as a semiclassical limit was noted for the sinh-Gordon model in Ref. <cit.>). Since the quantum distribution function ρ_k does not distinguish between phonon and soliton degrees of freedom, this indicates that solitons can indeed be neglected at the level of classical hydrodynamics. This conclusion is supported by the numerical results of Section <ref>.



 §.§ Quantum corrections to classical kinetic theory


We now consider the quantum corrections to the kinetic theory of the classical Toda lattice obtained above. The quantum TBA equations written in terms of classical momenta, Eqs. <ref> - <ref>, allow for a direct comparison between quantum and classical results, yielding predictions at a given temperature and (classical) chemical potential that tend smoothly to the semiclassical results as ħ→ 0. This is shown in Fig. <ref> for the two-reservoir non-equilibrium steady state studied in Section <ref>. As the dimensionless scale m ω / ħγ^2 falls below 1, deviations from the classical limit become apparent. The main effect of increasing the relative magnitude of ħ appears to be a reduction in particle density in each reservoir, with a concomitant decrease in energy density.

Let us now quantify these corrections analytically. It turns out that the main subtlety arises from the quantum correction to the phase shift, which we focus on at first. Then, as an application, we calculate explicitly the first-order correction to the dressed energy. Our results can be straightforwardly extended to higher-order and other quantities.

The quantum correction to the phase shift is naturally defined as the difference between quantum and classical phase shifts. By Eqs. <ref> and <ref>, we have:

    δ K_p,p' = K_p/ħ,p'/ħ - K_p,p'^cl
    
        = 1/πγ(ln|p-p'|/γħ - Re{ψ(1+i  |p-p'|/γħ)}) .
   
We wish to analyze the semiclassical expansion of this quantity about ħ = 0. The final result will be given in Eq. <ref>. However, we find it instructive to first consider a wrong approach. Indeed, one might be tempted to perform an asymptotic expansion in the limit Δ k := Δ p / ħ→∞, which corresponds to the limit ħ→ 0 with fixed Δ p ≠ 0. Using the asymptotic expansion <cit.>

    Re{ψ(1+ik)}∼ln k + ∑_n=1^∞(-1)^n-1/2n k^2nB_2n,   k →∞,

of the digamma function (B_2n denote the Bernoulli numbers), one would conclude that leading quantum correction is of order ħ^2:

    δ K_p,p' = -ħ^2 γ/12π1/|p-p'|^2 , for   fixed   
    |p-p'|>0  .

However, the correction kernel on the right hand side is pathological because of the divergence at p-p' → 0. In fact, the expansion <ref> gives diverging kernels at all orders ħ^2n, and must be unphysical. This indicates that the quantum correction comes essentially from static (“soft”) scattering p-p' → 0 (more precisely, |p-p'|∼ħγ).   
 
This is indeed the case. In fact, we claim that δ K_p,p' is concentrated at p-p'=0 at any order of ħ. More precisely, 

    δ K_p,p' = ∑_n=0^∞κ_2n/(2n)!δ^(2n)(p-p') γ^2nħ^2n+1 , 
    where κ_2n =  1/π∫_-∞^∞ d u   u^2n[ ln |u|  - Re{ψ(1+i u)}]  ,

and δ^(s) is the s-th derivative of δ. Curiously, the semiclassical expansion of the phase shift has only odd positive powers of ħ, although the phase shift is an even function of (p-p')/ħ. In particular, the leading order is of order ħ (its coefficient can be explicitly evaluated):

    δ K_p,p' = ħκ_0 δ(p-p') + 𝒪(ħ^2)    , κ_0 = -1/2 .


To show Eq. <ref>, we apply the kernel δ K_p,p' to an arbitrary smooth function f_p of p, write the integral in variables k = p/ħ , k' = p' / ħ at the quantum scale, and then Taylor expand f_p' around p:

    ∫_-∞^∞  dp' δ K_p,p' f_p'    = 
     ħ∫_-∞^∞  dk' δ K_k,k' f_k'ħ
        = ħ∫_-∞^∞  dk' δ K_k,k'∑_n=0^∞1/n! f^(n)_kħħ^n (k'-k)^n

where f^(n) = d^n f_p/ dp^n. Such an expansion is physically justified because f_kħ is a very slowly varying function of k. At each order of ħ, the integral over k' is decoupled from f_p by a change of variable (u = (k'-k)γ)

    ∫_-∞^∞  dk' (k'-k)^nδ K_k,k' 
        = ∫_-∞^∞ d u/π[ ln |u|  - Reψ(1+i u) ] u^nγ^n  = κ_nγ^n ,

where we used Eq. <ref>. Combining Eqs. <ref> and <ref>, we obtain

    ∫_-∞^∞  dp' δ K_p,p' f_p' =
       ħ∑_n=0^∞1/n! f^(n)_p (ħγ)^n  κ_n ,

which is equivalent to Eq. <ref>, since Eq. <ref> vanishes for odd n by parity.



Eq. <ref> allows the semiclassical expansion of the kinetic theory to any order. To illustrate this, we compute the leading quantum correction to the dressed energy. To this end, we expand it as follows:

    ϵ^cl = ϵ^cl^(0) + δϵ^cl + 𝒪(ħ^2)  .







It follows from the Yang-Yang equation in the form Eq. <ref> that the perturbation δϵ^cl_p satisfies the integral equation

    δϵ^cl_p + ∫_-∞^∞ dp' K^cl_p,p'θ^(0)_pδϵ^cl_p'
    
    =     -1/2(βħω)∫_-∞^∞ dp' K^cl_p,p'ω e^-2βϵ^cl^(0)_p' + ∫_-∞^∞  dp' δ K_p,p'ω 
    e^- βϵ^cl^(0)_p'

up to O(ħ) corrections. Applying Eq. <ref>, we obtain the integral equation that determines the leading quantum correction to ϵ^cl_p:

    δϵ^cl_p + ∫_-∞^∞ dp' K^cl_p,p'δϵ^cl_p'
    
        =    -ħβω/2∫_-∞^∞ dp' K^cl_p,p'ω e^-2βϵ^cl^(0)_p' - 1/2ħω e^-βϵ^cl^(0)_p  + O(ħ^2)  .

This expression is exact and can be verified numerically, see Fig. <ref> (inset) for an example. Observe that in Eq. <ref>, the only kernel that relates unequal momenta p - p' > 0 is the classical one K_p,p'^cl. This is a direct consequence of Eq. <ref> and holds at any order and for all other quantities (e.g., the velocity functional). We also note that the first term of the semiclassical expansion Eq. <ref> is controlled by the dimensionless parameter λ = βħω, so is valid at fixed nonzero temperature. 


at order 𝒪(ħ). There are generically two types of quantum correction to the classical kinetic equation Eq. <ref>. The first arise at the level of initial conditions; these are 𝒪(ħ) corrections arising from the leading quantum corrections to TBA, and are controlled by 

    λ_1 = βħω ,
 
where β denotes a typical inverse temperature scale. Such corrections include the coupling to gradients of temperature in the two-reservoir initial condition. The second type of corrections arise due to quantum corrections to the classical phase-shift, and are more subtle to analyze. and alter the dependence of the velocity on the Fermi factor. The size of such corrections is controlled by the velocity scale 

    λ_2 = (ħγ/ m)(βħω)  .
 
Therefore, the dominant quantum corrections to the classical kinetic theory of the Toda lattice are the “static” 𝒪(ħ) corrections to the thermodynamics of local equilibrium states; the “dynamical” corrections occur at one order higher.






As an illustrative example, we analyze the structure of 𝒪(ħ) corrections to the non-equilibrium steady state arising from two-reservoir initial conditions. Recall that for initial Fermi factors of the form

    θ_p(x,0) = θ_L + (θ_R - θ_L)Θ(x),

where Θ(x) denotes the Heaviside step function, the non-equilbrium steady state satisfies <cit.>

    θ_p (x,t) = θ_L + (θ_R - θ_L)Θ(x-v_p[θ(x,t)]t).

After a brief calculation, which makes use of the linear degeneracy of the velocity functional <cit.>, we find that the linear response perturbation to this expression in ħ is given by

    δθ_p(ζ) = δθ^eq_p(ζ) - (θ^(0)_R - θ^(0)_L) (∫_-∞^∞ dp'  δθ^eq_p'(ζ) δ v_p/δθ_p'(ζ))δ(p-p^*(ζ))/|v'_p*(ζ)|,

where ζ = x/t, p^*(ζ) is defined by[This assumes monotonicity of the dressed group velocity in p, which holds for bare group velocities in the Toda lattice.] v_p^*(ζ)[θ^(0)(ζ)] = ζ and

    δθ^eq_p(ζ) = δθ^L_p +(δθ^R_p-δθ^L_p)Θ(ζ - v_p[θ^(0)(ζ)]).

The first term in Eq. <ref> is the leading quantum correction to TBA in each individual reservoir, while the second term arises from differentiating the Heaviside function (with a regulator ϵ=0^+ implicit). These terms can be viewed respectively as “bulk” and “gradient” corrections to the two-reservoir steady state. Both terms are linear in δθ^eq, and by the above results, scale with ħ through the dimensionless coupling λ_1 = βħω.





§ NUMERICAL STUDIES

In this section, we numerically test the classical kinetic theory in two ways.  In Section <ref>, we test the kinetic theory directly in a non-equilibrium two-reservoir setting. In Section <ref>, we connect the quasiparticles in kinetic theory to the spectrum of the Lax matrix. For convenience, we set γ = m = ω = 1 throughout this section.



 §.§ Out-of-equilibrium dynamics

We now compare predictions from the classical kinetic theory with microscopic numerical simulations of the classical Toda lattice. While performing direct microscopic simulations on the quantum Toda lattice is not feasible, previous comparison of results from integrability to microscopic DMRG simulations on the XXZ chain <cit.> suggest that the quantum kinetic equation is comparably accurate in the hydrodynamic limit.

For the sake of simplicity, we focus on time-evolution from “two-reservoir” initial conditions, for which the solution to the hydrodynamic equations <ref> is well-known by now <cit.>. To obtain hydrodynamic predictions, the reservoir occupation numbers θ^L/R_p are computed from the classical TBA equations Eq. <ref>, from which the occupation numbers θ_p(x,t) for the hydrodynamic steady-state are obtained by solving the system

    θ_p (x,t) = θ^L_p + (θ^R_p - θ^L_p)Θ(x-v_p[θ(x,t)]t)

by numerical iteration, as usual.


We then obtain, via Eq. <ref>, profiles of particle density and energy density

    q_0(x,t) = ∫ dp  ρ_p(x,t)  ,  q_2(x,t) =  1/2∫ dp   p^2 ρ_p(x,t)  .

These predictions can be compared to direct numerical simulations. The equations of motion are integrated using the Flaschka variables <cit.>, see Eq. <ref> below. The thermal equilibrium of each reservoir can be sampled directly, since the momenta and positions have the following factorized distribution:

    P({r_n, p_n }) = 1/Z∏_n e^-βp_n^2/2exp( -β e^-r_n - β P  r_n )  ,  r_n = x_n+1 - x_n  ,

where Z is the partition function and P>0 is the pressure. β and P determine the classical thermal state, of which the chemical potential and average stretch are
    
    μ^cl =  -1/β( lnΓ(β P) - ln√(2π) + (1/2 - β P) lnβ)  , 
             a =  ∂μ^cl/∂ P  ,

respectively. For a fixed β, μ^cl(P) is not a one-to-one function (see inset of Fig. <ref>), so there is an ambiguity in defining a thermal state by β and μ^cl, as in TBA. This is lifted by further requiring 
    a ≥ 0  ,

which selects the low-pressure branch. We observe that the latter solution is systematically selected by the solution (by numerical iteration) of the classical TBA equation <ref>. Intuitively, this is because a state with a < 0 is inaccessible from the asymptotic Bethe Ansatz that we assumed initially (see Section <ref> above).

To sample a two reservoir initial condition with 2N particles, we generate p_1, …, p_N, r_1, …, r_N using Eq. <ref> with β_L, μ_L, and p_N+1, …, p_2N, r_N+1, …, r_2N-1 with β_R, μ_R. The overall translation freedom is fixed by requiring (x_N + x_N+1)/2 = 0.
A result of the comparison is shown in Fig. <ref>. We find an excellent agreement already at a moderate space-time scale, x, t ∼ 10^2. At shorter times, the effect of higher-derivative corrections becomes visible, especially for the particle density.





   
  



 §.§ Lax matrix and quasiparticles

It is well-known that the equation of motion of the classical Toda lattice can be expressed in terms of a Lax pair <cit.>, as

    L̇ = [B, L]   .

For an open chain, B and L are the following tridiagonal matrices, which we shall regard as operators on the Hilbert space spanned by the orthonormal basis | n ⟩, n = 1, …, N:

    L = ∑_n=1^N  b_n | n ⟩⟨ n |  + 
        ∑_n=1^N-1 a_n (| n ⟩⟨ n + 1 | + 
        | n + 1⟩⟨ n  |)  , 
    
        B =  1/2∑_n=1^N-1 a_n 
        (| n + 1⟩⟨ n  | - | n ⟩⟨ n + 1 | )  .

Here,

    b_n = p_n  ,  a_n = e^-r_n/2

are the Flaschka variables (up to a factor 2). Eq. <ref> implies that the spectrum of the Lax matrix L is conserved, as are the quasimomenta  in the kinetic theory. Can they be identified? Here we show that this is the case; moreover, the quasiparticles appear as eigenstates of the Lax matrix. 

First, we compare the quasiparticle density ρ(k) with the density of states P(λ) of the Lax matrix for a thermal state at inverse temperature β and pressure P. The former is determined by the TBA equation <ref>. P(λ) is the averaged density of states for the random Lax matrix ensemble whose matrix elements a_n, b_n are obtained from the thermal distribution <ref> (via Eq. <ref>). In terms of the Green's function, we have

    P(λ) =  1/i π.1/NTr[G(λ - iϵ)] |_ϵ→ 0_+ ,  G(z) := (z - L)^-1 .

where […] denotes the ensemble average. Note that L is Hermitian and has a real spectrum. We computed P(λ) by exact diagonalization for a few thermal states, and found that it agrees perfectly with ρ_p (this depends on the normalization of Eq. <ref>), up to a normalization constant given by the average stretch:

    P(λ) = a  ρ_p , where  λ = p  .

One example is shown in Fig. <ref> (top panel). 

We next turn to the quasiparticle currents, given by

    j_p = v_p[ρ]ρ_p = 1/2πθ_p ϵ^cl_p'
  
in the kinetic theory. To obtain their (conjectural) expression in terms of the Lax matrix, consider an eigenstate |λ⟩ of L with eigenvalue λ. It evolves under

    d/dt|λ⟩ = B |λ⟩ .

We want to interpret this eigenstate as a quasiparticle. Its position is given by an expectation value

    x(λ) = ∑_n⟨λ| n ⟩ x_n ⟨ n |λ⟩ ,  x_n = r_n-1 + r_n-2 + … + r_1 + x_1  ,

which is not well-defined in terms of the Flaschka variables (because of the overall displacement). However, the time derivative of x(λ) can be calculated without ambiguity using Eq. <ref>:

    ẋ(λ) =  ⟨λ| V |λ⟩ , 
    
        V := ∑_n  [ b_n | n ⟩⟨ n | 
         +  1/2 a_n ln a_n ( | n ⟩⟨ n + 1 | + | n + 1 ⟩⟨ n |)   ]  .

Summing over eigenstates, we have

    j_p = 1/i π.1/NTr[ V  G(p - iϵ)] |_ϵ→ 0_+ .

We numerically evaluate the RHS for thermal states by exact diagonalization and compare it to the kinetic theory prediction, obtaining again an excellent agreement, see Fig. <ref> (bottom panel). 


The identification of quasiparticles with Lax eigenstates (which may be called a Bethe-Lax correspondence) is conceptually appealing: it explains why the former can be regarded as point-like objects at the level of kinetic theory. While this is straightforward in simple systems such as hard rods (a quasiparticle is literally a hard rod, modulo permutations), it is not a priori clear for the classical Toda lattice. An intuitive justification at non-zero temperature is as follows: since L is a tridiagonal matrix with random elements of order unity, it can be viewed as a one-dimensional, disordered tight-binding Hamiltonian. It is then well-known that its eigenstates are all Anderson localized <cit.>. Therefore, the quasiparticles have a position uncertainty of order unity, and can be treated as points at the hydrodynamic scale. This is illustrated in Fig. <ref>, inset, where we show a space-time snapshot of Toda lattice evolution, from the viewpoint of scattering quasiparticles.




§ CONCLUSION

We have studied the kinetic theory of the quantum and classical Toda lattices, from the perspective of an interacting gas of particles on a line. While the quantum kinetic theory was obtained from TBA in a manner that is straightforward at this point, and the classical kinetic theory can be obtained directly <cit.>, our analysis of the semiclassical (ħ→ 0) limit is physically instructive and allows us to quantify the size of quantum corrections through a systematic expansion in ħ; this shows in particular that the classical limit is non-singular. This may not be surprising given that the quantum Toda lattice can be viewed as an infinite-spin limit of SU(2) integrable spin chains <cit.>, and so is already “classical” in some sense. 
Here, we found a hydrodynamical consequence of this fact, in terms of the semiclassical expansion of the phase shift (Eq. <ref>): at any order in ħ, the quantum corrections only affect scatterings in the classically static limit. We suspect that this is a general phenomenon in the semiclassical hydrodynamics of integrable systems. For instance, this holds for the Lieb-Liniger model, for which K^cl_p,p' = 0 is moreover trivial. It will be interesting to analyze the implications of this for the dynamics of the trapped δ-Bose gas, as was realized in the famous “Quantum Newton's Cradle” experiment <cit.>.

Our numerical studies validate the classical kinetic theory by a direct comparison with microscopic non-equilibrium dynamics, as well as justifying the underlying equations of state from the viewpoint of a Bethe-Lax correspondence. This relates predictions from kinetic theory with mathematical conjectures about certain tridiagonal matrix ensembles, which are closely related to the β-random matrices <cit.>. This relation was recently exploited <cit.> in order to construct the generalized Gibbs ensemble of the classical Toda lattice from a random-matrix point of view. A similar approach might be useful in developing a kinetic theory for the Calogero-Sutherland-Moser systems, which admit a similar Lax-pair formulation. It should be clear that mathematically, the Bethe-Lax correspondence is nothing more than the inverse scattering method. However, in the context of the hydrodynamics of integrable systems, it has not previously been realized in a concrete way. Doing so for quantum integrable systems (including the quantum Toda lattice <cit.>) might significantly elucidate the connection between kinetic theory and microscopic quantum dynamics, and represents an exciting challenge for future study.



Acknowledgments. We thank Herbert Spohn for illuminating discussions and for sharing his results prior to their publication in <cit.>. We acknowledge support from the DARPA DRINQS program (VBB), the ERC synergy Grant UQUAM, DOE grant DE-SC0019380 (XC), and the Center for Novel Pathways to Quantum Coherence in Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, and a Simons Investigatorship (JEM). Numerical computations are carried out on the cluster of Laboratoire de Physique Théorique et Modèles Statistiques (Orsay, France).






