A lattice gas model for generic one-dimensional Hamiltonian Systems

We present a three-lane exclusion process that exhibits the same universal fluctuation pattern as generic one-dimensional Hamiltonian dynamics with short-range interactions, viz., with two sound modes in the Kardar-Parisi-Zhang (KPZ) universality class (with dynamical exponent $z=3/2$ and symmetric Pr\"ahofer-Spohn scaling function) and a superdiffusive heat mode with dynamical exponent $z=5/3$ and symmetric L\'evy scaling function. The lattice gas model is amenable to efficient numerical simulation. Our main findings, obtained from dynamical Monte-Carlo simulation, are: (i) The frequently observed numerical asymmetry of the sound modes is a finite time effect. (ii) The mode-coupling calculation of the scale factor for the $5/3$-L\'evy-mode gives at least the right order of magnitude. (iii) There are significant diffusive corrections which are non-universal.


Introduction
ciently, as detailed below. It turns out that deviations from the scaling predictions are indeed strong for early times, but there is a clear indication of convergence for larger times, as discussed in the conclusions. The model allows for varying the strength of the mode coupling coefficients and can thus be employed to simulate scenarios that arise in other types of models which are less amenable to numerical analysis. In our model the conserved densities corresponding to the hydrodynamic densities of mass, momentum and energy are identified as linear combinations of the three particle densities in the three different lanes. Ergodicity guarantees the absence of further conservation laws.
2 The three-lane partial exclusion process 2

.1 Definition
Consider a three-lane asymmetric partial exclusion process with up to m particles per site, L sites per lane, and periodic boundary conditions. We denote by n λ k ∈ {0, 1, . . . , m} the occupation number on site k on lane λ with λ ∈ {−1, 0, 1}. Thus the time-dependent numbers n λ k (t) represent the time evolution of a single realization of the stochastic process.
The dynamics is Markovian. Particles hop randomly to nearest neighbour sites on the same lane with rates that depend not only on the occupation numbers of the departure site k and the target site k ± 1, but also on those of the neighbouring sites on the same and on the neighbouring lanes. The rates r λ k for a particle jump from site k to site k + 1 and λ k from site k to site k − 1 are The parameter range is a i ≥ 0, b i +d i ≥ −a i /m, b 1 +b 2 ≥ −a 0 /m to ensure positivity of all jump rates. The inter-lane coupling strength is given by the constants b 1,2 and we shall assume at least one of them to be non-zero to have interaction between the lanes. Complete or partial decoupling also takes place if one of the lanes is empty or completely filled. We exclude these trivial cases from our considerations. The case m = 1 corresponds to coupled exclusion processes. Due to the periodic boundary conditions the total number of particles N λ = k n λ k in each lane is conserved.

Steady state properties 2.2.1 Stationary distribution
For parameters z λ ≥ 0 the product measures with site marginals are a family of stationary distributions which one proves by observing that both the right-hopping process and the left-hopping process individually leave (7) invariant, which comes from a cancellation of the terms in the lattice sum over k, and similarly for the other hopping terms. Notice that the exclusion parameter m does not appear in these equations. The fugacity z λ parametrizes the density on lane λ obtained from (7). The total particle number N λ in each lane is conserved under the dynamics, but is a fluctuating quantity among realizations in the grand canonical ensemble defined by the invariant measure (7). The product form yields the diagonal density covariance matrix K with matrix elements with Expectations of time-independent functions f in the stationary distribution are denoted by f .

Currents
The so-called instantaneous currents for the bond (k, k+1) in lane λ are the functions of the time-dependent occupation numbers n λ k (t) which represent a single realization of the stochastic process. For an arbitrary initial distribution P 0 of the particles, these instantaneous currents yield the microscopic continuity equations for the timedependent expected local density ρ λ k (t) := n λ k (t) P 0 as Here the brackets denote averages over histories of the process and an arbitrary initial distribution P 0 of the particles. For the stationary currents j λ , which are translation invariant, one has j λ = j λ k . The product structure of the stationary distribution yields The subscript here corresponds to the superscript in the previous equations for the instantaneous current. Notice that the stationary currents depend only on the difference of the individual rates a 1,2 , b 1,2 , d 1,2 and are independent of a 0 . The exclusion parameter m only renormalizes the densities, the asymmetry parameter a and the overall time scale. From the stationary current density relation (13) -(15) one obtains the current Jacobian J with matrix elements We are specifically interested in the symmetric case ρ + = ρ − =: ρ. In this case the current Jacobian takes the form The Hessians H ν are defined by the matrix elements Defining one has for equal densities ρ ± = ρ By construction, the Hessians are symmetric.

Dynamical structure function
The (real-space) dynamical structure function is the matrixS k (t) with matrix ele-mentsS λµ Here the brackets without subscript denote an average over histories of the process and the stationary distribution. Because of the factorized stationary distribution one hasS λµ k (0) = κ λ δ λ,µ δ k,0 (25) and, due to particle number conservation, with the compressibility matrix defined in (9). From the continuity equation (12), the absence of stationary correlations, and the ensemble property (9) one also has for the infinite lattice the exact property where J is the current Jacobian (17). The matrix product JK is symmetric as has been proved in [15], thus linking the purely static compressibility K with the dynamics encoded in J.
The dynamical structure function describes the flow and broadening of density fluctuations in the steady state. To elucidate this property we diagonalize the current Jacobian J, i.e., we study RJR −1 = V =: diag (v − , v 0 , v + ) with eigenvalues v α . The three eigenvectors of J define the real-space eigenmodes which are fluctuation fields that travel with velocities v α . We normalize the transformation R by where K is the compressibility matrix defined in (9). With the transformation matrix R to normal modes, one obtains the dynamical structure function for the eigenmodes Thus, with the normalization (29) the transformed dynamical structure function S satisfies The first relation normalizes the "mass" of each mode and the second relation gives the propagation velocities of the modes. We stress that these are exact relations valid on the infinite lattice for all finite times t ≥ 0. The lattice Fourier transformŜ(p, t) of the dynamical structure function is defined byŜ with p of the form p = m/L with m integer. We can also consider the correlation between the Fourier-transformed density fieldŝ from which one constructsS Exploiting translation invariance one findsS(p, q, t) = LŜ(p, t)δ p+q,0 and thereforê Fourier eigenmodes are calculated similar to (28) aŝ and from (35) one concludeŝ Monte Carlo simulations are performed for a periodic system of very large length L = 10 6 with fixed particle numbers N λ . The densities are then given by ρ λ = N λ /L. Since the exclusion parameter m is immaterial from a theoretical perspective and we look for optimal numerical efficiency we choose full exclusion corresponding to m = 1. Notice that on switching to a system with fixed particle numbers, the stationary distribution modifies and each configuration with N λ occupied sites on lane λ becomes equally likely to be observed. After the drawing of an initial configuration from this uniform stationary distribution the system is evolved in time by using random sequential update.
Making use of translation invariance and ergodicity, we define the Monte-Carlo estimator for the structure functions where l + k has to be taken modulo L. Further, τ is the time between measurements and M is the total number of measurements. In order to compute the structure functions we generate P independent initial configurations yielding Monte Carlo estimatorsσ (p) L,k for each initial configuration. Averaging over the initial configurations then yields the numerical structure function which depends on the simulation parameters L, P, M, τ , on the space and time parameters k, t, and on the model parameters. On choosing P and M sufficiently large the Monte Carlo error for the numerical structure function becomes Gaussian distributed with zero mean and standard deviation scaling as ∼ 1/ √ P M . Our values for these simulation parameters are given in Sec. (4).

Canonical ensemble correction for finite-size effects
In the grand canonical ensemble the dynamical structure function vanishes rapidly at finite time t outside the "light cone" which is the region enclosed by the modes with the lowest and highest velocity. For numerical purposes, we define the light cone by the interval where v l (v h ) is the lowest (highest) mode velocity, z l (z h ) is the dynamical exponent of the corresponding mode and c l (c h ) is a constant chosen such that correlations outside are indeed vanishing within statistical measurement precision, after taking care of the following remark. For a finite system the constant numbers of particles N λ = Lρ λ introduce longrange correlations extending over the whole lattice even for t = 0. The canonical invariant measure is uniform, which for m = 1 yields the static structure function for fixed particle numbers N λ . Thus in a canonical ensemble at time t > 0 and k outside the light cone one has, within numerical measurement precision, Transforming to eigenmodes one gets Thus the structure function in the canonical ensemble has a constant offset of order L −1 compared to the grandcanonical structure function used above to describe correlations in the infinite system and which does not require this correction.
For the numerical precision of our Monte-Carlo data this offset is relevant and needs to be taken into account. Correcting for the finite size effects exposed in (42) we will use the quantities However, other finite size effects do appear in either ensemble. First of all, at times t > L/(v h − v l ) different modes will meet and add interactions over which we have no control. At even longer times, proportional to L z λ , the discreteness of the allowed set of values for the Fourier variable p becomes noticeable in the various mode coupling contributions. To avoid having to deal with this we will restrict our comparisons between simulations and theory to times shorter than this. The system size is chosen as large as numerically possible, but large enough to allow for neglecting within Monte-Carlo accuracy, the finite size effects discussed above.

Choice of rates
We recall that only the rate differences a, b, d defined in (16) occur in the currents and therefore in the mode-coupling matrices. For efficient simulation, avoiding irrelevant jump processes, we choose the jump rates of the three-lane model as follows: When choosing the rates in this way, the hopping rates simplify to and result in less storage accesses during the simulation. Additionally, for a ≥ 0 this choice guarantees to have all rates greater than or equal to 0 and we avoid having the processes r − k and l + k for a 1 = 0. Notice that we have chosen the rates such that l + = r − ; l − = r + ; ρ + = ρ − and the rates r 0 and l 0 depend on the occupations of the neighboring lanes in a symmetric way. This way we are guaranteed to find the velocities in the form v ± = ±v and v 0 = 0, just like for the sound modes and the heat mode in one dimensional Hamiltonian systems. Below we will use these terms to refer to the modes in our model system as well.

Predictions from mode coupling theory
In generic Hamiltonian dynamics as well as in the present model with our choices of parameters all three velocities v α are different [1] and this will be used throughout the discussion in this section. Under these conditions one expects the off-diagonal elements of the dynamical structure function (30) to decay fast and on large scales one is left with the diagonal terms S α (x, t) with continuous space coordinate x.
We follow our previous work [8] and use the mode coupling equations with memory term to predict the large scale-scale behaviour of S α (x, t). Here the G α are the mode coupling matrices obtained from the Hessians H α (21) by the transformation The diffusion coefficients D α turn out to be immaterial for the theoretical predictions that arise from (59). By construction, the mode coupling matrices are symmetric and related to the mode coupling matrices W α of [1] by W α = 2G α .

Dynamical scaling
In the scaling limit t → ∞ and x → ∞ with constant scaling variable u α = E α (x − v α t)t −1/zα with a non-universal scale parameter E α , non-universal mode velocity v α and universal dynamical exponent z α the mode coupling equations (59) with memory term (60) can be solved exactly for all modes α and any number of conservation laws [8]. Thus one can determine from them self-consistently for each mode the dynamical exponent z α and the corresponding scaling form of the dynamical structure function S α (x, t) = (E α t) −1/zα s α (u). We note that generally, in the case of distinct mode velocities (v α = v β for all α = β) and the absence of purely diffusive modes δ (for which G δ αα = 0 for all modes α), any mode α with non-vanishing self-coupling coefficient G α αα is expected to be in the KPZ universality class with dynamical exponent z α = 3/2. Rather than determining the corresponding scaling functions by the mode coupling equations (59) we use in our approach the exact Prähofer-Spohn scaling s α (u) = f P S (u) [5].
For Hamiltonian dynamics this argument applies to the two sound modes with velocities v ± = ±v so that their asymptotic behavior is expected to become and speed v = |v ± | = 0. The scaling function f PS is known exactly and satisfies There is no expression in closed form for f PS which, however, has been calculated with high precision and is tabulated in [6]. Using this data a more precise calculation of c PS can be found in Eq. (74) of [8].
In a setup with three modes, two symmetric KPZ-modes traveling with v + = −v − = v, identical self-coupling |G + ++ | = |G − −− |, and a mode 0 with vanishing self-coupling G 0 00 = 0 and mode-velocity, but non-vanishing symmetric coupling |G 0 ++ | = |G 0 −− | to the KPZ modes, Eq. (59) predicts dynamical exponent z 0 = 5/3 for this mode. The corresponding scaling function is then a symmetric 5/3-stable Lévy distribution given by and Here Γ(·) is the Gamma function. At this point we would like to emphasize that the mode coupling solution for the heat mode differs from [1] only in the prefactor a h,[1] = 1.6712. We recall that non-linear fluctuating hydrodynamics [2] predicts the presence of diffusive corrections to these asymptotic results. Defining the Fourier transform aŝ the Fourier representation of the asymptotic heat scaling function with diffusive correction included is given by Here D 0 is a phenomenological constant.

Eigenmodes and mode coupling matrices of the threelane model
The scenario described above is expected for generic Hamiltonian dynamics with short-range interactions. Specifically, one requires in the frame of vanishing centerof-mass velocity v 0 = 0, v ± = ±v (73) and the mode coupling symmetry [1] with G − −− = 0 and G 0 −− = 0. Hence, to verify these requirements, one needs to compute all diagonal mode coupling coefficients G α ββ for the three-lane exclusion process.
Indeed, the characteristic polynomial of the Jacobian (18) yields the eigenvalues We choose the diagonalizing matrix R such that Together with the normalization (29) this yields Notice the following symmetry which will play a role below.
With the sum of the Hessians (23) we find The mode coupling symmetry (74) is indeed satisfied without any fine-tuning of parameters, as can be seen as follows.
The Hessians (23) exhibit the symmetry One finds for the independent non-vanishing coefficients Some steps of the lengthy computation are presented in the appendix.
4 Monte Carlo results for the structure-function 4.1 Data fit using L 1 -distance The prediction for the heat mode and the exact Prähofer-Spohn scaling function involve the the non-universal scale parameters λ (given in (63)) and E 0 (given in (69)). In order to allow for a comparison with Monte-Carlo results we define analogously to [12] the fitted scaling parameters λ fit α as the minimum of the L 1 -distance where K is the set of calculated sampling points. For a model independent comparison of simulation results and theory we consider analogously to [12] the scale factor coefficients

Overview
With Monte Carlo simulation we aim at exploring the asymptotic behavior in numerically accessible times. Therefore, to guarantee a quickly vanishing mode overlap, the model parameters a, b, d and the particle densities ρ, ρ 0 should be chosen such that the speed v of the KPZ sound modes becomes as large as possible. Furthermore, the scale factors λ, E 0 for the sound and heat modes should be as large as possible so as to dominate any diffusive (z = 2) finite time contribution. However, according to (69) a strong KPZ selfcoupling G + ++ and large sound velocity v weakens the 5/3-Lévy-mode scale factor E 0 . Thus, to access the asymptotic regime of both modes we have to find a system with well balanced scale factors for both modes.
Below we present two examples showing the following properties: 1. Good KPZ-modes and a 5/3-Lévy-mode with diffusive finite time effects.
In both examples the convergence to the predicted asymptotic results becomes apparent after subtraction of the diffusive correction.
We first consider the fit of the scale parameters obtained from the L 1 distance, as shown in Fig. 1. We observe a monotone convergence of the numerical L 1 distance and of the scale parameters λ ±,0 and a s,h resp. to their theoretical values. At the largest time t = 2.33 · 10 5 the scale parameter a fit s of the KPZ sound mode is very close to the exact theoretical value (precision ≈ 1%). The scale parameter a fit h of the heat mode, however, is significantly off (≈ 34%) the theoretical prediction.
In order to explore the heat mode further we make a fit in Fourier space and include a diffusive correction. In this way we find where D 0,fit is the fitted diffusion coefficient of the heat mode. This result corresponds to a deviation of only ≈ 7% from the theoretical result. This indicates that diffusive corrections are significant. However, since t is not sufficiently large to be in the asymptotic regime one cannot tell from this result whether the 7% difference between numerics and theory is due to residual finite-time corrections or to an imprecision of the mode-coupling approximation. For a more detailed analysis we plot the value of the maximum of the sound modes together with the exact theoretical value (Fig. 2) as well as a scaling plot of the full dynamical structure function in position space (Fig. 3). One finds excellent convergence of the numerical data to the theoretical curve.   The absence of scaling suggests the presence of strong diffusive corrections. This is confirmed by inclusion of a diffusive correction in the Fourier transform, as shown in Fig. 5. The data collapse is excellent.
Notice, in comparison to the first example, the smaller value of the sound mode self coupling coefficients |G ± ±± | and the much larger value |G 0 ± | of the heat mode coupling coefficient which lead to correspondingly different values of λ and E 0 .
As Monte Carlo parameters we choose One finds as in example 1 a convergence of the scale parameters to the theoretical asymptotic values (Fig. 6). However, at the largest time t = 233000 there are still significant deviations from the asymptotic values, both for the KPZ sound mode (≈ 6%) and the Lévy heat mode (≈ 14% from below). The correction to a 0 for the heat mode is non-monotonic. By looking at a scaling plot for the full dynamical structure in position space one notices that the KPZ sound mode does not exhibit a good data collapse (Fig. 7). Also the amplitude at the maximum shows deviations from the exact theoretical result that are significantly larger than in the first example (Fig. 8).  On the other hand, the Lévy heat mode exhibits quite good data collapse (Fig. 9). According to the conclusions drawn from the first example this should be indicative of small diffusive corrections. This is confirmed by studying the scaling function in Fourier space with diffusive corrections. One finds as fit parameters

Conclusions
The main conclusions from this numerical study can be summarized as follows.
• We have found a three-lane lattice gas model that exhibits the same stationary fluctuations (two KPZ sound modes with velocities ±v = 0 and Lévy heat mode with v = 0) of the three conserved densities as generic Hamiltonian dynamics with short-range interaction and conservation of mass, energy, and momentum (and no additional conservation laws).
• The mode coupling matrices of this lattice gas model have the same structure of nonzero elements as mode coupling matrices of generic continuum Hamiltonian dynamics.
• No fine-tuning of model parameters is required to observe this behaviour.
• Choosing exclusion dynamics for the lattice gas model allows for long-time simulations and very large system size, using a canonical ensemble correction scheme introduced here (Sec. 2.3.2) and which is important to the numerical accuracy of the simulation results.
• Diffusive corrections generically persist for long times but their impact can be varied by changing system parameters.
The present work also reconfirms that mode coupling theory predicts correctly the predicted universal scaling form of heat mode, but the precision of non-universal scale factor that generally appears in Lévy modes is still an open problem. Also a rigorous mathematical proof of the KPZ/KPZ/5/3-Lévy scaling for this lattice model or suitably chosen variants remains a challenge. The remarkable rigorous results of [16] on the 3/2-Lévy scaling in the case of an anharmonic chain with two conservation laws [10,17] provide hope that these questions can be answered. =