High-low pressure domain wall for the classical Toda lattice

We study the classical Toda lattice with domain wall initial conditions, for which left and right half lattice are in thermal equilibrium but with distinct parameters of pressure, mean velocity, and temperature. In the hydrodynamic regime the respective space-time profiles scale ballisticly. The particular case of interest is a jump from low to high pressure at uniform temperature and zero mean velocity. Thereby the scaling function for the average stretch (also free volume) is forced to change sign. By direct inspection, the hydrodynamic equations for the Toda lattice seem to be singular at zero stretch. In our contribution we report on numerical solutions and convincingly establish that nevertheless the self-similar solution exhibits smooth behavior.


Introduction
Over the past decade the out-of-equilibrium dynamics of many-body systems in one dimension has received a lot of attention, see for example the recent survey article [1]. Besides the laws governing the dynamics, also the initial conditions have to be specified. One popular choice is quantum quench [2,3], for which one prepares the ground (or thermal) state for a particular hamiltonian and then evolves in time according to some other translation invariant hamiltonian. A further much studied choice, the one discussed in our contribution, is domain wall: In the left/right half lattice one prepares thermal states of a given translation invariant hamiltonian. The product coupling of these states evolves then under the full hamiltonian. If in either half the thermodynamic parameters are identical, one deals with a small perturbation of a global equilibrium state. But, if the parameters differ, the time evolution is far from equilibrium. Besides being a physically natural choice, macroscopic properties can be predicted on the basis of a suitable hydrodynamic theory.
In our contribution we focus on the classical Toda lattice with initial domain wall. The dynamics is integrable and, in principle, generalized hydrodynamics (GHD) can be used to quantitatively describe the ballistic spreading [4,5]. The average current of GHD is of the form ν −1 (v eff −q 1 ), where the effective velocity, v eff , is the solution of a rate equation and q 1 is the average momentum [18]. ν denotes the average stretch (also free volume), which is defined through ν = q j+1 −q j with q j the position of the j-th Toda particle. In thermal equilibrium ν is independent of j. For low pressure the stretch is positive, ν > 0, since particles are far apart. The stretch decreases with increasing P and at some critical pressure, P c , the stretch vanishes, while ν < 0 for P > P c , compare with Fig. 1(d) below. In the latter regime the particles are typically anti-ordered, i.e. q j > q j+1 , whereas for ν > 0 they are typically ordered.
The low-high pressure domain wall refers to imposing a pressure P − < P c in the left half lattice and P + > P c in the right half, with uniform temperature and zero mean velocity throughout. As a consequence ν = ν − > 0 in the left half lattice and ν = ν + < 0 in the right one. On the basis of the hydrodynamic equations, the stretch is expected to reach a self-similar form as ν(x, t) = g(x/t) with some suitable scaling function g. Thus g(−∞) = ν − , g(∞) = ν + and the scaling function is forced to vanish at some intermediate spatial location, g(x c ) = 0, and hence ν(x c t, t) = 0. But at this location the average current seems to diverge because of the prefactor ν −1 . Thus a more detailed investigation of the self-similar solution is required, which is the topic of the present work. The case of a low-low pressure jump has been studied before [20], including low order quantum corrections. In this case g stays strictly positive and varies smoothly, as anticipated. Of course, a corresponding behavior is expected for a high-high pressure domain wall.
Already in the pioneering contributions [4,5] it was realized that the solution of the GHD domain wall problem has a very specific structure. In the position-(spectral parameter) plane one has to determine a contact line at which the spectral parameter jumps from its left to its right value. This contact line fully characterizes the solution, which brings us to a second motivation for our study. While the entire self-similar solution of generalized hydrodynamics is computed numerically, most commonly only the spatial dependence of density, momentum, and temperature are displayed. Rather than such low order moments, we regard the contact line as more informative and will numerically determine its shape.
To give an outline: We briefly review the domain wall solution for non-integrable chains. Generalized hydrodynamics of the Toda lattice is recalled. Our main result are numerical solutions of GHD for the low-high pressure domain wall. As a by-product we elucidate the structure of the density of states for the Lax matrix as P varies through P c .

Domain wall for non-integrable chains
We briefly recall the domain wall problem of conventional hydrodynamics, so to provide a contrast with integrable chains. We start by considering a wave field over Z, displacements denoted by {q j ∈ R, j ∈ Z}, which is governed bÿ Here U is a smooth function, bounded from below, with at least a one-sided linear growth at infinity. Note that the right side depends only on the increments of q. Physically one should think of nonlinear springs coupling neighboring displacements with U the corresponding potential and −U the force. The dynamics is hamiltonian and obtained from q j , p j position and momentum of the j-th particle. For the harmonic potential, i.e., U (x) = 1 2 x 2 , Eq. (2.1) becomes the standard discrete linear wave equation. For a non-linear force, rather naively one would expect the dynamics to become chaotic. But based on the Kolmogorov-Arnold-Moser (KAM) theorem a mixed phase space, partially chaotic and partially quasi-periodic, is the more likely scenario. On the other hand, the domain wall setting is concerned with the infinite lattice, and how much of the KAM structure persists is a poorly understood subject. We proceed here with a scheme which seems to work in practice.
It will be convenient to first introduce the increments with e j the translate of e 0 to lattice site j. e 0 is a smooth function of its arguments. The equations of motion (2.1) then turn intȯ From the dynamics one reads off three local conservation laws, namely stretch, momentum and energy, r j , p j , e j (2.6) and their currents − p j , −U (r j−1 ), −p j U (r j−1 ) .
While this derivation is straightforward, in principle there could be some other smooth function h 0 , its translate to site j denoted by h j , which is strictly local in the sense to depend only on finitely many r j 's and p j 's and moreover satisfies the discrete conservation law If such a current density J j exists at all, it has to be smooth and strictly local. As obvious from (2.8), the time change of j=m j=m h j is then only through the boundary of the interval [m, . . . , m ]. As an outstanding puzzle in the subject, a strict dichotomy seems to be valid. The nonlinear wave equation (2.1) has either three or a countably infinite number of local conservation laws, the former case being called non-integrable. Integrable chains are exceptional. The obvious candidate is the harmonic chain, U (x) = 1 2 x 2 , see for example the discussion in [10]. Since the equations of motion are linear, the superposition principle holds. Thus, beyond being integrable the harmonic chain is also non-interacting. The respective GHD is a system of uncoupled conservation laws. The Toda chain, to be discussed below, is also integrable. But now GHD is a nonlinear system of coupled conservation laws, reflecting that conserved fields interact.
To construct an initial domain wall state for a non-integrable chain, note that the dual thermodynamic parameters are (P, u, β), pressure, mean velocity, and inverse temperature. Since H of (2.4) is a sum of one-point functions, in thermal equilibrium the {r j , p j }'s are independent and identically distributed with single site probability density function Z the normalizing partition function and P > 0, β > 0. The thermal state does not change under the dynamics (2.5). For a domain wall state we prescribe left and right parameters, (P ± , u ± , β ± ). The (r, p)'s are still independent, but for j < 0 we choose the parameter set (P − , u − , β − ) and for j ≥ 0 the parameter set (P + , u + , β + ). One would expect that on a macroscopic scale the state can be well characterized by block averaged values of the conserved fields. Given these values, the local distribution is close to the corresponding equilibrium state. Assuming such local equilibrium, one arrives at the coupled set of hydrodynamic equations where the hydrodynamic currents are given by Here l(x, t), u(x, t), e(x, t) are the hydrodynamic fields of stretch, velocity, and total energy. P is the pressure depending on stretch and internal energy. From the microscopic model P is obtained by setting u = 0 in Eq. (2.9) and computing the averages r 0 P,0,β , e 0 P,0,β . By inverting these functions one arrives at P . The Euler equations have to be solved with domain wall initial conditions In our context of a coupled set of hyperbolic conservation laws, (2.10) -(2.12) are known mathematically as Riemann problem, which has been thoroughly studied, see [11] for a very readable account. Here, let us merely observe that the solution is self-similar in the sense that it depends on x, t only through the ratio x/t. The solution consists of spatial intervals either with constant profiles or smoothly varying profiles, the latter known as rarefaction waves. Possibly they are separated by jump discontinuities, called shocks. To our knowledge there are only a few studies which compare molecular dynamics of the chain with domain wall initial conditions to numerical solutions of the respective hydrodynamic equations. In [12,13] the program is carried out for a hamiltonian with point hard core potential and alternating masses. Actually in this investigation the goal was to elucidate fluctuations of time-integrated currents at a location inside a rarefaction wave. In fact, the most fascinating part of the study lies in observing how the microscopic motion manages to create and maintain a jump discontinuity. Along the same lines is the recent numerical study of blast formation [14].
The solution to the Riemann problem relies on properties of the Euler equation linearized at local equilibrium, in our case a 3 × 3 matrix with matrix elements depending on the thermodynamic parameters. As the location x of the self-similar solution is varied, the corresponding thermodynamic parameters change and, in order to satisfy the boundary conditions, may be forced to switch between two branches characterized by distinct eigenvalues. In a rough sense, this is the mathematical mechanism behind the formation of shocks. In contrast, for integrable systems the linearization operator has continuous spectrum and no mechanism for shock formation seems to be available.

GHD of the Toda lattice
The Toda lattice is an anharmonic chain, for which the interaction potential is specified as U (x) = e −x . Hence the hamiltonian is written as In terms of the Flaschka variables [15], the Lax matrix is the tridiagonal real symmetric matrix with matrix elements For the finite lattice [1, . . . , N ] with periodic boundary conditions the N eigenvalues of L are conserved. As functions on phase space they are non-local [16], their local version being tr[L m ], m = 1, 2, . . . [15]. In generalized hydrodynamics such conserved fields are usually called charges or conserved charges and it is convenient to follow this practice. The locally conserved charges of the Toda lattice have a strictly local density given by with j ∈ Z. In addition, the stretch is conserved with density The respective current densities are of the form where L ↓ denotes the lower triangular part of L, see [17,21,22]. By construction one arrives at the microscopic conservation laws Because of the extensive number of conserved charges, thermal equilibrium has to be extended to generalized Gibbs ensembles (GGE). Generalized hydrodynamics is then obtained by averaging (3.7) in a local GGE state. There is some freedom in the choice of coordinates. The most convenient one comes from the observation that the GGE expectation of Q [m] is the m-th moment of the density of states (DOS) of the Lax matrix. Therefore, the natural hydrodynamic fields are the DOS of the Lax matrix and in addition the stretch, ν, both depending on the macroscopic space-time point (x, t). More details can be found in [18][19][20] and the very recent lecture notes [22]. In standard notation, the DOS is written as νρ p (v), where ρ p is the "particle density in spectral parameter space". Then the generalized hydrodynamic equations are Here q 1 = ν R dwwρ p (w) is the average momentum. The definition of the effective velocity is more lengthy. One introduces the integral operator resulting from the two-particle scattering shift of the Toda lattice. Then, for given ρ p , the effective velocity is the solution of the linear integral equation where ρ p (v) is viewed as a multiplication operator. Writing GHD in this way, the domain wall problem looks inaccessible. Surprisingly, as a general property of GHD, one can transform to normal modes in such a way that the quasi-linear operator is diagonal. This transformation is accomplished by (3.11) and the resulting normal form of the hydrodynamic equations reads ν∂ t n + (v eff − q 1 )∂ x n = 0, (3.12) see [21] for a proof. Note that the two-system (3.8) has merged into a single equation. The symbol n(v) denotes the "number density in spectral parameter space". In [17,21,22], the same object is denoted by ρ µ . Below we need some further notions from GHD. The dressing transformation of a general function ψ is given by One also uses the notation [w m ] dr for the dressing of the m-th power of the linear function. Note that the dressing is relative to n. After a few algebraic transformations based on (3.10), (3.11) and (3.13), the thermodynamic quantities can be expressed in terms of n, which will be convenient for the domain wall problem discussed below. Specifically, one arrives at (3.14) A seemingly more explicit formula for the effective velocity is  with V (w) = 1 2 βw 2 . The "chemical potential" µ is a Lagrange multiplier which ensures normalization as R dw n(w) = P.
(3.18) µ depends on the thermodynamic parameters β and P , and turns out to be equal to the free energy F (β, P ), for which an explicit formula is available [17,22], F (β, P ) = log β/(2π) + P log β − log Γ(P ), (3.19) with Gamma function, Γ. The thermally averaged stretch is then the derivative of the free energy with respect to the pressure, ν = ∂ P F (β, P ) = log β − ψ(P ), (3.20) ψ denoting the digamma function. Fig. 1 visualizes these thermal functions. F (β, P ) is a concave function, and assumes its maximum at the critical pressure P c > 0. Thus the average stretch ν decreases with increasing P and crosses zero at P c . For low pressure particles are far apart and the density is small. As P is increased, the free volume between particles with neighboring index shrinks and vanishes at P c . Upon further increase the free volume turns negative. From the shape of F (β, P ) (and thus of µ) one concludes that for given µ there are two different values of P . In (3.17) only µ appears as parameter. Thus this equation has two solution branches, low and high pressure, which match at P c . As P is moved through P c , νρ p and n vary smoothly. On the other hand, ρ p has to diverge as P → P c , since ν = 0 at P c . Numerically we investigated only thermal equilibrium. But any other GGE is expected to show the same behavior.  Numerical method. -We first discuss the computation of the thermal n(w): TBA is solved numerically via Newton iteration (Mathematica's FindRoot), together with the quasi-energy parameterization n(w) = e −ε(w) to ensure n(w) > 0. In practice, this requires a good starting point, which we obtain via an associated Fokker-Planck equation [23]. As sketch of a derivation, first insert n(v) = P ρ s (v) into (3.17), which yields Differentiating this equation with respect to w and then multiplying by ρ s (w) results in ρ s can be interpreted as the stationary solution of the time-dependent nonlinear Fokker-Planck equation Requiring the normalization R dwρ s (w) = 1, Eq. (3.23) has a unique stationary solution. While TBA and Fokker-Planck are equivalent analytically, we observed that, numerically, optimal precision is obtained by first solving the nonlinear Fokker-Planck equation, and then using these data as initial value for the TBA Newton iteration.
Regarding GHD, the key task is to realize the dressing transformation (3.13) numerically. According to (3.13), this amounts to solving a linear system of equations after discretization. We have found a finite element discretization of n(v), ρ p (v), v eff (v), etc. via hat functions on a uniform grid to work well in practice. Specifically, we use the grid spacing h = 1 10 . The action of the integral operator (3.9) can be symbolically precomputed for these hat functions, such that the discretization of T is a symmetric matrix.
Starting from a given n(v), one subsequently calculates [1] dr (v) and [w] dr (v), then ρ p (v) via (3.14) and v eff (v) via (3.15), as well as ν and q 1 via (3.16). Due to the discretization on a uniform grid, the integration in (3.16) amounts to a simple summation of the hat function coefficients of the integrand.
Regarding the domain wall problem discussed in the following section, n φ (v) turns out to be a piecewise combination of thermal functions, see Eq. (4.4) below. For such a n φ (v), we use the algorithm just described to obtain [1]

Solution to the domain wall initial condition
We slightly rewrite Eq. (3.12) as The domain wall initial conditions are Instead of n ± , physically it might be more natural to prescribe the DOS of the Lax matrix and the average stretch. But mathematically the normal form (4.1) is more accessible.
Since the solution to (4.1), (4.2) scales ballisticly, we set n(x, t; v) = n(t −1 x; v) andṽ eff (x, t; v) =ṽ eff (t −1 x; v). Without loss of generality one adopts t = 1 and arrives at In the example of the harmonic chain, the parameter v corresponds to the wave number which varies only over a bounded interval. In this case the contact line function cannot be inverted. On the other hand, for hard rods [6], v varies over the entire real axis and the contact line is monotone increasing. Since the Toda lattice is closer to hard rods, we will assume that the contact line is given by the inverse ofφ, denoted then by φ.
With such an assumption, for every x there is a contact point φ(x). Then the solution ansatz reads The superscript φ will be used to generically indicate that in the TBA formalism n φ is substituted for n and similarly the subscript ± signals the substitution of n ± . For example, compare with (3.13), (3.16), Then the condition on the left side of (4.3) translates to x =ṽ eff (x, φ(x)).    (4.9) which means that G is the inverse of the contact line φ. While G(φ) itself has to be obtained numerically, the large |φ| asymptotics can still be argued analytically. We start from Setting we note that for large φ the second term is exponentially small and can be neglected. The second summand on the left of (4.10) then reads Similarly for the right side of (4.10) one finds a logarithmic increase with a prefactor c + , which could be determined by a further iteration. For φ → ∞ the term q φ 1 /ν φ converges to q 1,+ /ν + . Thus one arrives at the large φ asymptotics, For φ → −∞, the same asymptotics holds upon substituting ν − , c − for ν + , c + . If the boundary functions n ± have exponential decay, the error in (4.13) would be of the same order.
Exemplary numerical solution. -The considerations above hold for an arbitrary choice of boundary values n ± . For the simulation we impose thermal boundaries, characterized by zero average momentum, constant inverse temperature β = 1, and pressures P − = 1 2 , P + = 2, respectively, such that P − < P c < P + . The results are shown in Fig. 2. Starting from n φ (v) defined in (4.4), all the other hydrodynamic functions can be obtained via dressing transformations. While n φ (v) pointwise agrees with n ± (v) by construction, see Fig. 2a, this is no longer the case for the normalized particle density ν φ ρ φ p (v) shown in Fig. 2b. For comparison, the dashed curves in Fig. 2 visualize the functions computed for thermal n ± (v).
Only a single low-high pressure domain wall has been simulated in detail. We tested a few other values. As physically to be expected, the resulting plots look rather similar. Only the condition P − < P c < P + has to be respected.

Toda fluid and related work
So far we viewed the Toda lattice as a discrete nonlinear wave equation. A physically more immediate perspective would be to have particles moving on the real line, which is referred to as Toda fluid. Instead of stretch, one then considers the fluid density ρ f . For a homogeneous system |ν|ρ f = 1. Including space-time variations, the mapping between lattice and fluid is discussed in [18]. Rather unexpected features appear for the low-high pressure domain wall, however. We consider the initial state and set q 0 = 0. Then {q j+1 −q j , j ≤ −1} are i.i.d. random variables and so are {q j+1 − q j , j ≥ 0}. By assumption q j+1 − q j = ν − > 0 for j ≤ −1, while q j+1 − q j = ν + < 0 for j ≥ 0. A typical ordering is of the form · · · < q −1 < q 0 = 0 > q 1 > . . . . Thus initially the average particle density equals ν − + |ν + | on (−∞, 0] and decays exponentially on [0, ∞). Under the dynamics, the point at which the ordered domain touches the anti-ordered domain is moving in time, its label being denoted by κ(t). Close to q κ(t) particles pile up. The domain boundary acts as bottleneck for particles. The position of the bottleneck is x bn (t) = q κ(t) . Our numerical simulations suggest that κ(t) and x bn (t) change linearly in time, at least approximately. Also beyond the bottleneck position the particle density vanishes rapidly. As observed numerically, the scaling function g has a nonzero slope at x c . Thus, near the bottleneck the particles should be distributed according to equilibrium with a linearly varying pressure. From this property one infers that in physical space the particle density at the bottleneck diverges as an inverse square root. It would be of interest to better understand the precise particle statistics close to the bottleneck.
Within GHD, domain wall initial conditions have been studied numerically also for a discrete version of the sinh-Gordon model [26]. Most detailed studies are available for the XXZ model [27,28]. In this case the spectral parameter space has in addition the type of string states and the contact line carries extra indices. Because of momentum conservation, the contact line of the Toda lattice is supported on the full real line. For XXZ, and other discrete models, one usually observes a light cone, i.e., the contact line is supported on an interval and the boundary values are maintained up to an edge. The behavior near the edge often shows then an intricate oscillatory decay, which has been elucidated in considerable detail [29,30]. On the Euler scale the self-similar solution has a sharp step at the contact line. One expects that, because of dissipation, the step is actually broadened to an error-like function. For the Toda lattice, and in general for integrable systems, the structure of diffusive corrections are available [24]. But respective numerical simulations for the Toda lattice are still missing. On the other hand for the XXZ chain, both on the microscopic and GHD level, broadening has been convincingly observed [25].