Dynamics of hard rods with initial domain wall state

TL;DR

This paper analyzes the non-equilibrium dynamics of classical hard rods with an initial domain wall state, providing solutions to hydrodynamic equations, exploring contact discontinuities, and discussing viscous effects and entropy production.

Contribution

It offers the first explicit solution to the Euler equations for hard rods with domain wall initial conditions and compares classical and quantum transport approaches.

Findings

Solution composed of contact discontinuities for each velocity

Viscous corrections lead to broadening of discontinuities

Entropy production analyzed in viscous regime

Abstract

Inspired by recent results on the non-equilibrium dynamics of many-body quantum systems, we study the classical hard rod problem in one dimension with initial domain wall condition. Hard rods are an integrable system, in the sense that for each velocity the density of particles is locally conserved. It was proven by Boldrighini, Dobrushin and Sukhov (1983) that on the hydrodynamic space-time scale, the fluid of hard rods satisfies Euler-type equations which comprise all conservation laws. We provide the general solution to these equations on the line, with an initial condition where the left and right halves are, asymptotically, in different states. The solution is interpreted as being composed of a continuum of contact discontinuities, one for each velocity. This is a classical counterpart of the transport problem solved recently in quantum integrable systems. We provide two independent derivations, one paralleling that in the quantum case, the other using explicitly the microscopic dynamics. Finally, we discuss the Navier-Stokes (viscous) corrections, and study its effect on the broadening of the contact discontinuity and on entropy production.

Equations

$$Z^{-1}\chi\big{(}\{|q_{i}-q_{j}|\geq a,1\leq i<j\leq N\}\big{)}\prod_{j=1}^{N}h(p_{j})\,\tfrac{1}{N!}dq_{j}dp_{j}$$

$$\mathfrak{n}(x;v)=\sum _{j\in \mathbb{Z}}\delta (x-q_j)\delta (v-p_j).$$

$$\mathfrak{n}(x,t;v)=\sum _{j\in \mathbb{Z}}\delta (x-q_j(t))\delta (v-p_j),$$

$${\partial }_t\mathfrak{n}(x,t;v)+{\partial }_x\mathfrak{j}(x,t;v)=0$$

$$\mathfrak{j}(x,t;v)=\sum _{j\in \mathbb{Z}}{\dot{q}}_j(t)\delta (x-q_j(t))\delta (v-p_j),$$

$$\dot{q}_{j}(t)=p_{j}+a\sum_{j\neq k\in{\mathbb{Z}}}|p_{k}-p_{j}|\big{(}\chi(p_{k}<p_{j})\delta(q_{k}(t)-q_{j}(t)-a)-\chi(p_{k}>p_{j})\delta(q_{k}(t)-q_{j}(t)+a)\big{)}$$

$${\int }_t^{t+\tau }ds{\int }_{\Delta }dv\mathfrak{j}(x,t;v),$$

$$⟨\mathfrak{n}(x,t;v){⟩}_f=f(v),\frac{⟨\mathfrak{j}(x,t;v){⟩}_f}{⟨\mathfrak{n}(x,t;v){⟩}_f}=v+a(1-a\rho {)}^{-1}{\int }_{\mathbb{R}}dw(v-w)f(w),$$

$$\rho (x,t)={\int }_{\mathbb{R}}dvf(x,t;v),u(x,t)=\frac{1}{\rho (x,t)}{\int }_{\mathbb{R}}dvvf(x,t;v).$$

$$\partial_{t}f(v)+\partial_{x}\big{(}v_{[f]}^{\mathrm{eff}}(v)f(v)\big{)}=0.$$

$$v_{[f]}^{\mathrm{eff}}(v)=v+a(1-a\rho {)}^{-1}{\int }_{\mathbb{R}}dw(v-w)f(w)=v+\frac{a\rho (v-u)}{1-a\rho },$$

$$v_{[f]}^{\mathrm{eff}}(v)=\frac{v-a\rho u}{1-a\rho }.$$

$${\partial }_t{u}_{\alpha }+{\partial }_x{j}_{\alpha }(\vec{u})=0.$$

$${\partial }_t{u}_{\alpha }+\sum _{\beta =1}^n{A}_{\alpha \beta }(\vec{u}){\partial }_x{u}_{\beta }=0,A_{\alpha \beta }={\partial }_{u_{\beta }}{j}_{\alpha }.$$

$$n(v)=(1-a\rho {)}^{-1}f(v).$$

$$\partial_{t}(1-a\rho)+\partial_{x}\big{(}v_{[f]}^{\mathrm{eff}}(v)(1-a\rho)\big{)}=-a\big{(}\partial_{t}\rho+\partial_{x}(u\rho)\big{)}=0$$

$${\partial }_{t}n(v)+v_{[f]}^{\mathrm{eff}}(v){\partial }_{x}n(v)$$

$$\displaystyle=(1-a\rho)^{-1}\big{(}\partial_{t}f(v)+v_{[f]}^{\mathrm{eff}}(v)\partial_{x}f(v)\big{)}-(1-a\rho)^{-2}f(v)\big{(}\partial_{t}(1-a\rho)+v_{[f]}^{\mathrm{eff}}(v)\partial_{x}(1-a\rho)\big{)}$$

$$=-(1-a\rho {)}^{-1}f(v){\partial }_x{v}_{[f]}^{\mathrm{eff}}(v)+(1-a\rho {)}^{-1}f(v){\partial }_x{v}_f^{\mathrm{eff}}(v)$$

$${\partial }_{t}n(v)+v_{[f]}^{\mathrm{eff}}(v){\partial }_{x}n(v)=0.$$

$$f_0(x,v)=f_{-}(v)\mathrm{for}x<0,f_0(x,v)=f_{+}(v)\mathrm{for}x>0.$$

$$u_{\alpha }(x,t)={\mathsf{u}}_{\alpha }(x/t).$$

$$\Big{(}\sum_{\beta=1}^{n}\big{(}\xi\delta_{\alpha\beta}-A_{\alpha\beta}(\!\vec{\,\mathsf{u}})\big{)}\Big{)}\partial_{\xi}\mathsf{u}_{\beta}=0.$$

$$f(x,t;v)=\mathsf{f}(\xi ;v),n(x,t;v)=\mathsf{n}(\xi ;v),v_{[f]}^{\mathrm{eff}}(v)=v^{\mathrm{eff}}(\xi ;v).$$

$$(\xi -v^{\mathrm{eff}}(\xi ;v){\partial }_{\xi }\mathsf{n}(\xi ;v)=0.$$

$$\mathsf{n}(\xi ;v)=n_{-}(v)\chi (v>v_{\ast }(\xi ))+n_{+}(v)\chi (v<v_{\ast }(\xi )),$$

$$v^{\mathrm{eff}}(\xi ;v_{\ast })=\xi .$$

$$\kappa (v)={\int }_v^{\infty }dw{n}_{-}(w)+{\int }_{-\infty }^{v}dw{n}_{+}(w),\nu (v)={\int }_v^{\infty }dww{n}_{-}(w)+{\int }_{-\infty }^{v}dww{n}_{+}(w)$$

$$\xi =v_{\ast }-a{\int }_{\mathbb{R}}dv(v-v_{\ast })\mathsf{n}(\xi ;v)$$

$$\displaystyle=v_{*}-a\Big{(}\int_{v_{*}}^{\infty}dv\,vn_{-}(v)+\int^{v_{*}}_{-\infty}dv\,vn_{+}(v)\Big{)}+av_{*}\Big{(}\int_{v_{*}}^{\infty}dv\,n_{-}(v)+\int^{v_{*}}_{-\infty}dv\,n_{+}(v)\Big{)}$$

$$=v_{\ast }-a\nu (v_{\ast })+a{v}_{\ast }\kappa (v_{\ast })=g(v_{\ast }).$$

$$g^{\prime }=1-a{\nu }^{\prime }+a\kappa +av{\kappa }^{\prime }=1+a\kappa .$$

Full text

Dynamics of hard rods with initial domain wall state

Benjamin Doyon∗ and Herbert Spohn*†*

∗* Department of Mathematics, King’s College London, Strand, London WC2R 2LS, U.K.

† Physik Department and Zentrum Mathematik, Technische Universität München, Boltzmannstrasse 3, 85748 Garching, Germany*

Inspired by recent results on the non-equilibrium dynamics of many-body quantum systems, we study the classical hard rod problem in one dimension with initial domain wall condition. Hard rods are an integrable system, in the sense that for each velocity the density of particles is locally conserved. It was proven by Boldrighini, Dobrushin and Sukhov (1983) that on the hydrodynamic space-time scale, the fluid of hard rods satisfies Euler-type equations which comprise all conservation laws. We provide the general solution to these equations on the line, with an initial condition where the left and right halves are, asymptotically, in different states. The solution is interpreted as being composed of a continuum of contact discontinuities, one for each velocity. This is a classical counterpart of the transport problem solved recently in quantum integrable systems. We provide two independent derivations, one paralleling that in the quantum case, the other using explicitly the microscopic dynamics. Finally, we discuss the Navier-Stokes (viscous) corrections, and study its effect on the broadening of the contact discontinuity and on entropy production.

1 Introduction

The study of non-equilibrium, deterministic dynamics of many body systems has received a great amount of attention in recent years. This is particularly true in the quantum context, where investigations into the unitary quantum evolution of extended systems has led to a wealth of new fundamental theoretical principles. For instance, in so-called homogeneous quantum quenches, a coupling parameter of the model is suddenly changed, and the question of thermalization and, for integrable systems, generalized thermalization has been much studied [1, 2, 3, 4, 5].

Inhomogeneous states and dynamics are, however, more difficult to approach. Of particular interest is the partitioning protocol for generating non-equilibrium steady states from an initial inhomogeneity, already studied in the 1970ies [6]. In the partitioning protocol, two halves of an infinite system are independently thermalized (a “domain wall” initial condition), then connected and let to evolve unitarily for a large time. If the dynamics admits ballistic processes, then at large times, local observables are in a stationary state with nonzero currents, for instance of energy or particles – this is a non-equilibrium steady state. In the quantum context, Ruelle’s work [7] provided a general mathematical framework, and rigorous results for non-equilibrium currents and some other observables in quantum chains admitting free-fermionic descriptions were obtained [8, 9, 10, 11, 12], see also [13]. The protocol received renewed attention after an exact general solution at one-dimensional quantum critical points was given [14, 15], with numerical verification [16], experimental evidence [17], and a mathematically rigorous proof in the context of algebraic conformal field theory [18]. See [19, 20, 21, 22, 23, 24, 26, 27, 28, 25, 29, 30, 31, 32, 33, 34, 35, 36, 37] and the reviews [38, 39].

The problem of finding exact solutions to the domain-wall initial-value problem in the context of interacting quantum integrable models stayed unsolved for some time, although various conjectures were given [40, 41, 42]. A breakthrough came [43, 44] when the powerful ideas of emerging hydrodynamics, whereby a many-body quantum system receives an effective description as a fluid [45, 46, 47, 48], were generalized to the infinite number of conservation laws present in integrable models. In this context, the hydrodynamic theory, called generalized hydrodynamics, has become a fundamental tool in order to study the partitioning protocol and more generally situations with inhomogeneous state or in inhomogeneous force fields, see [49, 50, 51, 52] and related ideas [53, 54].

Despite the large amount of work in the context of quantum integrable models, relatively few results exist in classical integrable models. In this note we study the domain wall problem for an integrable classical dynamics: that of hard rods located on the one-dimensional line. In the usual parlance the model is interacting, thus behaving very differently from an ideal gas and the harmonic chain, for which the exact solution in the case of a domain wall initial state is written down with ease. We mention that a detailed study for Fermi-Pasta-Ulam-type chains with domain wall initial state has been carried out recently [55]. But these are non-integrable systems with a dynamical behavior very different from the integrable case discussed here.

Hard rods are segments, say of length $a$, that move freely in one dimension except for elastic collisions whereby they exchange their velocities. Clearly, velocities are conserved implying an integrable dynamics. Under certain natural conditions, as proven in [56], at the Eulerian scale, a fluid of hard rods satisfies dynamical equations which account for all local conservation laws. Using these hydrodynamic equations, we provide an exact solution to the domain wall problem, for any initial distribution on the two halves. As is usual in the domain wall problem, this solution is self-similar, and at large times, it is valid for any initial distribution with the same asymptotics on the far left and right. Paralleling the case of quantum integrable models [43, 44], this solution is composed of infinitely many contact singularities, one for each velocity. We provide a derivation inspired by the recent discoveries in quantum integrable models as well as a complementary derivation that takes into account the microscopic dynamics.

One important problem in the study of non-equilibrium systems using hydrodynamic principles is to establish the viscous corrections to the Euler equations and their effects on the dynamics. This problem is still largely unsolved in the quantum context, see however the recent study [51]. For hard rod dynamics, the Navier-Stokes equation has been derived in [57], offering us the opportunity of studying its effects on the transport problem. We discuss such viscous corrections, and how they affect the broadening of the contact singularities and lead to the production of entropy. The non-vanishing entropy production signals that, albeit being integrable, hard rods are interacting.

The full relation between the generalized hydrodynamic equations found in quantum integrable systems and the classical hard rod dynamics will be developed in separate works [58, 59]. Here we simply mention that the hard rod problem can be seen, at the Euler scale, as a generalized hydrodynamic problem, but with a particular choice of the differential scattering phase appearing in the thermodynamic Bethe ansatz [60, 61].

The paper is organized as follows. In Section 2, we review the hard rods dynamics and Euler scale hydrodynamics. In Section 3, we derive the exact solution to the hydrodynamics with domain wall initial state. In Section 4, we discuss diffusive corrections and entropy production. The appendix provides comparison of the exact solution proposed with simulations of the hard rod gas.

2 Hard rod dynamics

2.1 Microscopic dynamics

Consider a one-dimensional fluid of hard rods. The hard rods have length $a>0$, positions $q_{j}\in{\mathbb{R}}$, and momenta $p_{j}\in{\mathbb{R}}$. They move freely except for elastic collisions conserving momentum and energy, whereby two hard rods exchange their momenta upon contact. For convenience we set their mass to $1$, and later will mostly refer to $p_{j}$ as velocity. In a collision, we label the particles so as to maintain their velocities. Thus particle $j$ moves along a straight line, $\dot{q}_{j}=p_{j}$, interspersed with jumps back and forth by $a$ due to collisions. Clearly, such a system with $N$ rods has $N$ conservation laws labeled by velocities.

Let us first discuss the generalized Gibbs ensembles (GGEs). A GGE is specified by fixing the volume, say the finite interval $[-L,L]$, the density $\rho=N/(2L)$ with $0<a\rho<1$, and the single-particle velocity probability density $h(v)\geq 0$ with $\int_{\mathbb{R}}dv\,h(v)=1$; that is, $h(p_{j})dp_{j}$ is the probability for particle $j$ to have a velocity between $p_{j}$ and $p_{j}+dp_{j}$. Then the canonical ensemble is given by the probability measure

[TABLE]

on phase space $([-L,L]\times\mathbb{R})^{N}$ (where $Z$ is the normalization factor). Here $\chi$ denotes the indicator function of the set specified. Under this equilibrium measure, velocities are independent of the positions and the velocities themselves are independent. The positions are however correlated because of the hard core repulsion. The infinitely extended system is obtained through the limit $N,L\to\infty$ at fixed $\rho$. The state is fully characterized by $f(v)=\rho h(v)$, and the corresponding equilibrium average will be denoted by $\langle\cdot\rangle_{f}$. In the following we will consider directly the infinite hard rod system. For the infinite system it is known that in the class of sufficiently regular measures, the only time-stationary measures are given by the above generalized equilibrium [56].

Let us now discuss the dynamics. It is convenient to define densities which will have well-behaved hydrodynamic limits. The particle number density, at time $t=0$, is

[TABLE]

That is, $\mathfrak{n}(x;v)dxdv$ is the number of particles in the square $[x,x+dx]\times[v,v+dv]$ of one-particle phase-space. This is a conserved field thanks to conservation of velocities under collisions, and our somewhat unusual notation emphasizes that $v$ is the label of the conservation law. The time-evolved fields are

[TABLE]

where we recall that the $p_{j}$’s are independent of $t$ according to our choice of labelling. The particle number density satisfies the local conservation law

[TABLE]

with current $\mathfrak{j}(x,t;v)$ explicitly given by

[TABLE]

where

[TABLE]

represents the linear motion with jumps at collisions.

The observables $\mathfrak{n}(x,t;v)$ and $\mathfrak{j}(x,t;v)$ are generalized functions on phase space. Integrating against a smooth localized function $g$ on phase space, as $\int dxdv\,g(x,v)\mathfrak{n}(x,t;v)$ or $\int dxdv\,g(x,v)\mathfrak{j}(x,t;v)$, yields functions on phase space. Currents can be defined also at a given location, $x$, by integrating in time and velocities as

[TABLE]

which counts the signed number of particles with velocities in the interval $\Delta$ crossing the given point $x$ during the time interval $[t,t+\tau]$. The averages of these observables in GGEs are given by

[TABLE]

see [62] Sec. I.3.3 for more details.

2.2 Eulerian hydrodynamics and normal modes

We now focus on the Eulerian scale. This is the large-distance, long-time scale, which emerges when conserved fields are slowly varying in space, and locally the system is approximately in a generalized equilibrium state. An analysis of the viscous correction (see sub-section 4.1) indicates that the Eulerian scale is controlled by the rod length $a$ as follows: it emerges, in some large region, when the space between rods throughout this region, in mean with respect to the local generalized equilibrium, times the mean variation length of the rod density, is much greater than $a^{2}$. We note that the local equilibrium may be understood in various ways: one may average over space, time, or space-time “fluid cells”, large enough to give negligible fluctuations, but small enough so that local states are slowly varying; and one may additionally average over realizations of the dynamics with different initial conditions determined by a given distribution, thus allowing the size of fluid cells to be reduced arbitrarily.

We therefore upgrade $f(v)$ to the space-time dependent function $f(x,t;v)$, which fully characterizes the local state at $x,t$. The evolution equation for $f$ is formally obtained by averaging (4) in a local equilibrium state. To write down the Euler type evolution equation let us first introduce the local density and velocity

[TABLE]

The $x,t$ argument will be omitted if obvious from the context. Using the local averages as provided in (8), the dynamics of $f$ is governed by

[TABLE]

Here

[TABLE]

which can also be written

[TABLE]

A tracer particle (in our context also called pulse) with “bare” velocity $v$ in fact moves with a corrected, or effective, velocity $v_{[f]}^{\mathrm{eff}}(v)$ arising from the collisions with other pulses of velocities $w$. The velocity $v_{[f]}^{\mathrm{eff}}(v)$ depends explicitly on $v$ and is a functional of $f$. For hard rods the dependence is only through the first two moments of $f$ and given by a linear rational function. For integrable quantum models, as for instance the Lieb-Liniger delta-Bose gas, the dependence on $f$ is more complicated [43, 44]. As a parenthesis we recall that (10) holds not only on average. Dobrushin et al. [56] prove a law of large numbers in the sense that typical initial configurations will follow the hydrodynamic evolution with a probability very close to one.

Let us briefly deviate from our topic by recalling the Eulerian hydrodynamic problem in the case of a finite number, say $k$, conservation laws, $u_{\alpha}$, $\alpha=1,...,k$, governed by

[TABLE]

Since currents are functions of densities through the equations of state, this can equivalently be written in quasilinear form as

[TABLE]

For small deviations from a uniform background, one has to diagonalize $A$. The corresponding normal modes propagate convectively with velocities given by the spectrum of $A$. In many other problems the quasilinear version is the easiest to deal with.

Above, the hydrodynamic variables are the $f(v)$’s. Accordingly, we may bring Eq (10) into the form (14), but the resulting linearized operator $A$ is of a complicated form. Surprisingly, there is a simple transformation diagonalizing $A$. We first introduce the “free density”

[TABLE]

As to be stressed there is no relation to the microscopic particle density $\mathfrak{n}$ introduced above. Then

[TABLE]

by conservation of mass. Hence

[TABLE]

and therefore $n(v)$ has vanishing convective derivative,

[TABLE]

The fields $n(v)$ are the normal modes of the hard rod fluid equations, and the convective equation identifies their propagation velocities as $v_{[f]}^{\mathrm{eff}}(v)$ with $f$ related to $n$ by (15).

The free density is closely related to the occupation function of the thermodynamic Bethe ansatz [60, 61] in the quantum problem; see [58, 59] for the details of this relation.

3 Domain wall initial state, Euler scale

We are interested in solving Eq. (10) with the domain wall initial conditions

[TABLE]

As before, let us briefly recall the Eulerian hydrodynamic problem in the case of $k$ conservation laws governed by (13). The domain wall initial state is invariant under spacial dilations, and the hydrodynamic equation (13) is invariant under simultaneous dilations of space and time. Therefore, if no spontaneous breaking of this symmetry occurs, the solution is of scaling form,

[TABLE]

Setting $\xi=x/t$, one obtains

[TABLE]

Thus the domain wall solution is related to the eigenvalues of $A$. Since $A$ has discrete spectrum, the actual construction of $\vec{\!\mathsf{\,u}}$, or of the normal modes, requires considerable efforts and we refer to [63] for a most readable introduction. As a result of such analysis, generically $\vec{\,\mathsf{u}}$ will have flat pieces interrupted by shocks, contact discontinuities, and rarefaction waves. Relevant for us are the contact discontinuities, which are defined by the property that the velocity of the jump equals to one of the eigenvalues of $A$. In our case, the $A$ matrix is infinite-dimensional and has only continuous spectrum. As will be discussed, there are only contact discontinuities.

3.1 Exact solution

We first provide a derivation that closely parallels that given in the quantum context [43, 44]. By scale invariance of the initial condition, one expects the solution to Eqs. (10), (19) to be self-similar, of the form

[TABLE]

Using (2.2), $\mathsf{n}$ satisfies

[TABLE]

Setting $n_{\pm}(v)=(1-a\rho_{\pm})^{-1}f_{\pm}(v)$, this equation together with the boundary condition (19) is solved by

[TABLE]

where $v_{*}$ is defined implicitly through

[TABLE]

Let us first make sure that (25) has a unique solution. Introducing

[TABLE]

and using $v^{\mathrm{eff}}(v)=v+a\rho(1-a\rho)^{-1}(v-u)$, one obtains

[TABLE]

Then

[TABLE]

Hence $g^{\prime}>0$, $g$ is invertible, and $v_{*}=g^{-1}(\xi)$. For large $\xi$, one has $v_{*}=\xi$. Thus the solution (24) indeed satisfies the boundary condition (19).

The final step is to switch back to the original field. Observing that $(1-a\rho)^{-1}=1+a\int_{\mathbb{R}}dv\,\mathsf{n}(v)$ one concludes

[TABLE]

Hence

[TABLE]

Reading the identity (25), the solution (24) is naturally interpreted as a family of contact discontinuities. This is the same type of solution as was found in [43, 44] in the case of generalized hydrodynamics for quantum systems. See Appendix A for comparisons between this exact solution and numerical simulations of the hard rod gas. Note that in the quantum context, unicity of the solution to (25) was only verified numerically. Here, since the effective velocity has a simpler form, unicity could be proved. The explicit solution (30) also is particular to the hard rod problem, but a generalization to other differential scattering phases can be obtained as integral equations [64].

3.2 Microscopic derivation

The domain wall solution (30) can be guessed by considering the microscopic hard rod dynamics and appealing to a suitable law of large numbers. So let us see how this works. At $t=0$ we shrink $a$ to zero. Then $f_{0}(x;v)$ goes over to $\tilde{f}_{0}(y;v)=(1-a\rho_{-})^{-1}f_{-}(v)\chi(y<0)+(1-a\rho_{+})^{-1}f_{+}(v)\chi(y>0)$, which evolves under the free dynamics to

[TABLE]

We have to compute how a small volume element at $y$ with velocity $v$ is shifted to $x$ upon expanding the hard rods back to size $a$. This shift has two contributions. Firstly, under the free dynamics the initial density jump at [math] is shifted to $vt$. Upon expanding back to the original hard rod length, for $vt<y$ the particles in $[vt,y]$ yield the shift $(1-a\rho_{+})^{-1}\rho_{+}(y-vt)a$, while for $vt>y$ the particles in $[y,vt]$ yield the shift $-(1-a\rho_{-})^{-1}\rho_{-}(vt-y)a$. Hence in sum

[TABLE]

A second contribution comes from other particles, say with velocity $w$, crossing the positions $y-v(t-s)$ at times $0\leq s\leq t$, which causes a shift $-a$ upon crossing from left to right and shift $a$ from right to left. Thus we have to compute $a$ times the total number of signed crossings as

[TABLE]

Using that $\tilde{f}_{0}$ is a step function one divides into the sets $\{y-wt<0\}$ and $\{y-wt>0\}$, respectively $\{y-vt<0\}$ and $\{y-vt>0\}$. Working out the integrals,

[TABLE]

The second summand precisely cancels the contribution from Eq. (32). We arrived at the expected ballistic scaling and may set $t=1$. Then the point $x$ is shifted to $x=y+a(-\nu(y)+y\kappa(y))=g(y)$ with volume element $dx=g^{\prime}(y)dy$. Inserting in (31) we arrive at the domain wall solution (30).

A re-statement of the above derivation is as follows. Consider the coordinate $y(x,t)=\int_{-x_{0}}^{x}dx^{\prime}\,(1-a\rho(x^{\prime},t))+x_{0}$. This is the free length from $x_{0}$ to $x$, after the rods have been shrunk to points. In terms of this coordinate, the solution to the domain-wall problem becomes quite transparent. We choose $x_{0}$ far enough on the left so that it lies within the left bath up to large enough times. Then, integrating (29) over $x$ at $t$ fixed, we have

[TABLE]

Since $v_{*}$, as defined in (25), is equal to $g^{-1}(\xi)$, and since $v_{*}\to\xi$ as $\xi\to-\infty$, then we may use $g^{-1}(x_{0}/t)=x_{0}/t$ for $|x_{0}|$ large enough. Therefore,

[TABLE]

That is, the function $g^{-1}(\xi)$, which determines the velocity at which a jump occurs in the solution (24) for $n(\xi;v)$, is nothing else but the ray in the $y$ coordinate. Such geometric ideas are at the basis of the viewpoint provided in [64].

4 Diffusive corrections

4.1 Navier-Stokes correction

On the Eulerian time scale, for fixed $v$, the domain wall solution has a jump discontinuity at $g(v)$. Physically one expects a broadening of the step, which analytically should be accessible through a suitable Navier-Stokes correction of the Euler equations (10). To figure out such a correction we follow the standard recipe for computing transport coefficients in equilibrium, as applied to hard rods in [65, 66].

For the velocity distribution we recall that $f(v)=\rho h(v)$ with $\int dv\,h(v)=1$, and we set $\int dv\,h(v)v=0$, for convenience. On large scales the static fields are $\delta$-correlated in $x$ and correlated in velocity space as

[TABLE]

with

[TABLE]

The second summand results from the integral over the static density-density correlation. Linearizing Eq. (10) as $\rho h(v)+f(v)$ yields the linearized Euler operator with integral kernel

[TABLE]

Hence

[TABLE]

Clearly $AC=CA^{\mathrm{T}}$, in agreement with linear fluctuating hydrodynamics [62].

In our context, the most interesting feature is the total current correlation function, which turns out to be given by

[TABLE]

The constant term is the Drude weight. As valid in great generality it equals $ACA^{\mathrm{T}}$, compare with the discussion in [67]. To obtain the diffusion matrix, $D$, we use the Green-Kubo formula according to which $DC$ equals the integral over the total current correlation with the Drude weight subtracted. Hence

[TABLE]

which by using (38) implies

[TABLE]

Remark. In Eq. I.7.67 of [62] it is claimed that

[TABLE]

which is incorrect. This oversight leads to a form of the nonlinear Navier-Stokes equations slightly different from the correct version given below in (47). Thereby a longstanding discrepancy as regards to [57] is resolved.

The operator $-D$ is the forward generator of a jump process on $\mathbb{R}$ with jump rate $|v-v^{\prime}|h(v^{\prime})$ for the transition from $v$ to $v^{\prime}$. The invariant measure is $h(v)dv$ and the rates satisfy detailed balance with respect to $h$. In particular, the quadratic form of $D^{\mathrm{T}}$ in $L^{2}(\mathbb{R},h(v)dv)=\mathcal{H}$ reads

[TABLE]

Hence $D^{\mathrm{T}}=(D^{\mathrm{T}})^{*}$ in $\mathcal{H}$ and the eigenvalues of $D^{\mathrm{T}}$, thus also of $D$, are nonnegative. The space-time correlator is well approximated by

[TABLE]

An explicit error bound can be found in [66].

On this basis one guesses that the Navier-Stokes correction is given by

[TABLE]

where the dissipative term is

[TABLE]

A mathematical proof is established by Boldrighini and Suhov [57]. Note that linearizing (47) as $\rho h(v)+f(v)$ one obtains indeed (4.1). As a further simple consistency check, according to (47) one has

[TABLE]

which must hold, since it is a property for each microscopic trajectory.

4.2 Entropy production

For usual hydrodynamics with a few conservation laws, there is a balance equation for the entropy which consists of a local entropy flow and a non-negative entropy production proportional to the transport coefficients. In fact, on the Euler scale, even starting with smooth initial data, shocks may develop spontaneously at some finite time and entropy is produced at shocks and possibly other non-smooth parts of the solution. For hard rods, most likely such a mechanism does not work. One evidence comes from [56], where the uniqueness of solutions for (10) is established without an entropy condition, which would be required in the case of a finite number of conservation laws. Another argument in support of the absence of the spontaneous formation of stable entropy-producing structures comes from the presence of a continuous set of conservation laws [43]. These offer an escape route permitting, on the Euler scale, entropy to be conserved at all times (except possibly for discrete sets of space-time points where infinitesimal amounts of entropy may be produced). For instance, in the domain-wall problem, where the eigenvalue equation (23) must be solved, the set of conservation laws provides a continuous infinity of eigenvectors, allowing for continuous fluid states111Since a fluid state is a function of $v$, continuity must be taken in an appropriate topology, which should be that induced by an appropriate completion of the set of local observables. to join any two reservoirs (an equivalent expression is that a continuum of contact singularities form – no entropy is produced at such singularities). This is in sharp contrast to the finite case, where usually no continuous solution exists and shocks must form. Thus, in the hard rod hydrodynamics, we expect no entropy-producing structures to spontaneously develop, and entropy production to result only from dissipation.

The thermodynamic entropy has two contributions. Firstly, there is the positional part including the hard core repulsion with entropy $s_{\mathrm{hr}}(\rho)=\rho\log(\rho^{-1}(1-a\rho))$ [65]. Secondly, since velocities are independent, the momentum part yields the entropy

[TABLE]

Hence the local entropy equals

[TABLE]

Note that this can be rewritten as

[TABLE]

in terms of the free density (15), which provides a natural interpretation: if the rods were shrunk to points, the resulting density of free particles yields the entropy.

We remark that this entropy function is different from the so-called Yang-Yang entropy function occurring in generalized equilibrium states of fermionic quantum systems [60, 61], even if both cases are described by similar equations of states. In the quantum case, there would be an additional term counting the hole configurations. This, of course, points to the fact that the entropy function is an additional characterization of the fluid, not a property of the Euler scale. Different fluids with the same Euler equations may have different entropy functions, with positive entropy production under different viscosity terms.

In Eulerian hydrodynamics, because of (18), the free density $n(x,t;v)$ has zero convective derivative, hence so does any function of it. Further, by (16) the function $(1-a\rho(x,t))$ is a conserved density. As a consequence, $(1-a\rho(x,t))n(x,t;v)\log n(x,t;v)$ is also a conserved density, with the same current. Integrating over $v$, the Eulerian contribution to $\partial_{t}s$ can be immediately evaluated. Taking into account dissipation, one therefore obtains

[TABLE]

where the dissipative term has been simplified using (49) and $\int_{\mathbb{R}}dv\,\mathcal{N}_{[f]}(v)=0$.

We rewrite this as

[TABLE]

The first summand is the additional flow term

[TABLE]

The second summand is the entropy production

[TABLE]

Thus the entropy balance reads

[TABLE]

with the entropy production defined in Eq. (56) and the entropy flow

[TABLE]

The velocity part of the entropy is transported by the effective velocity $v_{[f]}^{\mathrm{eff}}$ and the spatial part of the entropy by the mean velocity $u$. In addition, there is a contribution proportional to the spatial gradient, which is comparable to the entropy flow for compressible Navier-Stokes containing a term proportional to the temperature gradient.

4.3 Domain wall including dissipation

It seems rather unlikely to still find an explicit domain wall solution once the dissipative term is included. But on a more qualitative level the dynamical behavior can be easily guessed. The discontinuity at $\xi=g(v)$ will be smeared out diffusively to become an error function. The diffusion is correlated and nonlinear. But if one just assumes constant diffusion a simple picture emerges. Let us consider the density $\rho(x,t)=\rho_{\mathrm{s}}(x/t)$ with $\rho_{\mathrm{s}}$ the $v$-integral of Eq. (30). Dissipation generates then a convolution of this profile with a Gaussian of width $\sqrt{t}$. If $h_{\pm}$ are smooth, then $\rho(x,t)$ is smooth because of the $v$-integration. For long times the convolution will hardly change the profile. In fact the most visible change should be detected at short times when $\rho$ has still large gradients. The same observation applies to the other physical fields, such as velocity and energy.

5 Conclusions

Inspired by recent developments in quantum many-body integrable models, we provided the exact solution to the general domain-wall initial value problem for the hard rod gas on the Euler scale. The solution was obtained in two independent ways, one following the solution strategy of [43, 44], the other based on the explicit microscopic dynamics of the hard rod gas. We then studied viscosity and entropy, establishing the exact entropy current equation with a source controlled by viscosity, and showed positivity of entropy production.

It would be interesting to further understand the effects of viscosity on the domain-wall solution. Hard rod hydrodynamics is of course applicable beyond the exact self-similar limit of the solution to the domain initial value problem, and it would be instructive to study further its range of applicability. Finally, it would be very interesting to study the large deviation theory of the fluctuating currents, perhaps via an extension of the macroscopic fluctuation theory [68].

Appendix A Molecular dynamics simulations

In this appendix, we compare the exact solution (30) with numerical simulations of the hard rod dynamics. The domain wall initial state is taken to be at uniform temperature with a jump in the density. Another popular choice is the converse, uniform density and jump in the temperature, for which we expect similar results. Our simulation is performed at infinite volume. Initially, besides the discontinuity at the origin, the density jumps to zero at some point to the right and also to the left. Thus there are two depletion zones moving inwards. The exact solution looses validity once it merges with the depletion zones. For the exact solution, we use the formulae $\rho=(1-1/g^{\prime}(v_{*}))/a$ and $\rho u=\nu(v_{*})/g^{\prime}(v_{*})=(v_{*}-g(v_{*})/g^{\prime}(v_{*}))/a$ with $v_{*}=g^{-1}(x/t)$.

1. Rod length $a=0.001$, initial condition:

[TABLE]

In this initial condition, the value of $a\rho\in[0,1]$ is approximately $0.18$ on the left and $0.09$ on the right, which corresponds to $\langle N\rangle=4050$. The fluid cells are taken of length $0.1$, and an average over 200 samples is done. Excellent agreement is found for both density and current, see Figs 1 and 2. Note that even after seemingly small times, the system is well into the Euler scale, because the rod length and the rod separations are very small as compared to the scale at which the graph is displayed.

2. Rod length $a=0.001$, initial condition:

[TABLE]

In this initial condition, the value of $a\rho$ is approximately $0.36$ on the left and $0.12$ on the right, which corresponds to $\langle N\rangle=2400$. Thereby one probes more deeply into the hard rod dynamics, further away from the free-particle limit. The fluid cells are again taken of length $0.1$, and an average over 50 samples is done. See Fig 3. Since the segments of non-zero density are shorter, the depletion zone reaches the central region at an earlier time when compared to 1.

3. In order to see deviations from the exact solution at Euler scale, we investigate the case of a rod length $a=0.5$, with initial condition:

[TABLE]

In this initial condition, the value of $a\rho$ is approximately $0.55$ on the left and $0.25$ on the right. On average there are 25 rods in the fluid with an inter-particle distance of the order of unity. The fluid cells are again taken of length $0.1$, but there are so few rods that it is essential to average over a large number of samples, 120000 in our simulation. See Fig 4. We see observe a slight departure from the exact solution. It is surprising that, even with so few rods, the Euler-scale solution is still fairly accurate.