1 Introduction

It is by now well-established that one-dimensional short-ranged Hamiltonian systems with conserved particle density, energy, momentum generically exhibit universal time-dependent fluctuations of these quantities [1,2,3]. In the coordinate frame with zero center-of-mass velocity there are two oppositely moving sound modes. These are in the Kardar-Parisi-Zhang (KPZ) universality class [4] with dynamical exponent \(z=3/2\) and, in a comoving frame, given by symmetric Prähofer-Spohn scaling function [5, 6]. Moreover, mode coupling theory predicts a zero-velocity heat mode with dynamical exponent \(z=5/3\) with scaling form given by the 5/3-stable symmetric Lévy distribution [1, 2].

The fundamental assumptions underlying these predictions are that all slow variables of relevance for the long-time behavior of the time correlation functions are the long-wavelength Fourier components of the three conserved densities and that there are no further conservation laws. Thus one expects a broad universality of dynamical non-equilibrium phenomena in one dimension.

It should be noted, however, that these assumptions are not sufficient to guarantee the KPZ/Lévy/KPZ \((\frac{3}{2},\frac{5}{3},\frac{3}{2})\)-scaling scenario. In fact, it has been shown that the KPZ and 5/3-Lévy universality classes are members of an infinite family of dynamical universality classes whose dynamical exponents are – quite remarkably – given by the Kepler ratios of neighbouring Fibonacci numbers, beginning with \(z=2\) for diffusion and including also the limiting value, which is the golden mean \(\varphi =(1+\sqrt{5})/2\) [7, 8]Footnote 1. Specifically, for three conservation laws one can have among these Fibonacci universality classes normal diffusionFootnote 2 or marginal (logarithmic) superdiffusionFootnote 3 [10, 11], KPZ, modified KPZ (with unknown scaling function) [12] and Lévy (all with \(z=3/2\)), and Lévy modes with \(z=5/3\), \(z=8/5\) and \(z=\varphi \).

Which combinations of these eight different Fibonacci universality classes arise depends on the structure of the mode coupling matrices that can be obtained from the macroscopic current-density relation of the conserved densities [8, 12], and (only for marginal superdiffusion) on the presence of a cubic term in the underlying fluctuating hydrodynamic equation [10]. Moreover, notwithstanding the asymptotic predictions, one often observes numerically significant deviations from the predicted scaling behaviour which make the verification of the universality difficult in concrete applications [13, 14] or even put non-equilibrium universality into question [15].

Here we address these issues by introducing a lattice gas model for which mode coupling theory predicts the KPZ/Lévy/KPZ-scaling of generic Hamiltonian dynamics. One finds, without fine-tuning of parameters, the structure of the mode coupling matrices that is required for the \((\frac{3}{2},\frac{5}{3},\frac{3}{2})\)-scaling scenario, viz., non-vanishing self-coupling at quadratic order of the sound modes, vanishing self-coupling of the heat mode, but non-vanishing coupling of the heat mode to the sound modes.

This model is a three-lane exclusion process which can be simulated very efficiently, 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. The process is ergodic as can be seen by noting that the particle jumps are a sequence of permutationsFootnote 4 so that every configuration can be reached from any other after a finite number of jumps. Hence for any fixed set of particle numbers \(N_\lambda \) the invariant measure is unique which also implies 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^{\lambda } \in \{0,1,\dots ,m\}\) the occupation number on site k on lane \(\lambda \) with \(\lambda \in \{-1,0,1\}\). Thus the time-dependent numbers \(n_{k}^{\lambda }\left( t\right) \) 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\pm 1\), but also on those of the neighbouring sites on the same and on the neighbouring lanes. The rates \(r_k^{\lambda }\) for a particle jump from site k to site \(k+1\) and \(\ell _k^{\lambda }\) from site k to site \(k-1\) are

$$\begin{aligned} r^{-}_k= & {} \left\{ a_1 + \frac{b_1}{2} \left[ n^{0}_k + n^{0}_{k+1}\right] + \frac{d_1}{2} \left[ n^{-}_{k-1} + n^{-}_{k+2}\right] \right\} n^{-}_k \left( 1- \frac{n^{-}_{k+1}}{m}\right) \end{aligned}$$
(1)
$$\begin{aligned} \ell ^{-}_k= & {} \left\{ a_2 + \frac{b_2}{2}\left[ n^{0}_k + n^{0}_{k-1}\right] + \frac{d_2}{2} \left[ n^{-}_{k+1} + n^{-}_{k-2}\right] \right\} n^{-}_k \left( 1 - \frac{n^{-}_{k-1}}{m}\right) \end{aligned}$$
(2)
$$\begin{aligned} r^{0}_k= & {} \left\{ a_0 + \frac{b_1}{2} \left[ n^{-}_k + n^{-}_{k+1}\right] + \frac{b_2}{2} \left[ n^{+}_k + n^{+}_{k+1}\right] \right\} n^{0}_k \left( 1 - \frac{n^{0}_{k+1}}{m}\right) \end{aligned}$$
(3)
$$\begin{aligned} \ell ^{0}_k= & {} \left\{ a_0 + \frac{b_1}{2} \left[ n^{+}_k + n^{+}_{k-1}\right] + \frac{b_2}{2} \left[ n^{-}_k + n^{-}_{k-1}\right] \right\} n^{0}_k \left( 1 - \frac{n^{0}_{k-1}}{m}\right) \end{aligned}$$
(4)
$$\begin{aligned} r^{+}_k= & {} \left\{ a_2 + \frac{b_2}{2} \left[ n^{0}_k + n^{0}_{k+1}\right] + \frac{d_2}{2} \left[ n^{+}_{k-1} + n^{+}_{k+2}\right] \right\} n^{+}_k \left( 1- \frac{n^{+}_{k+1}}{m}\right) \end{aligned}$$
(5)
$$\begin{aligned} \ell ^{+}_k= & {} \left\{ a_1 + \frac{b_1}{2}\left[ n^{0}_k + n^{0}_{k-1}\right] + \frac{d_1}{2} \left[ n^{+}_{k+1} + n^{+}_{k-2}\right] \right\} n^{+}_k \left( 1 - \frac{n^{+}_{k-1}}{m}\right) . \end{aligned}$$
(6)

The parameter range is \(a_i \ge 0\), \(b_i+d_i \ge - a_i/m\), \(b_1+b_2 \ge - 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_\lambda = \sum _k n_{k}^{\lambda }\) in each lane is conserved.

2.2 Steady State Properties

2.2.1 Stationary Distribution

For parameters \(z_{\lambda } \ge 0\) the product measures with site marginals

$$\begin{aligned} \mathrm{Prob}\left[ \, {n_k^{\lambda } = n}\, \right] = \frac{z_{\lambda }^n}{(1+z_{\lambda })^m} {m \atopwithdelims ()n} \end{aligned}$$
(7)

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

$$\begin{aligned} \frac{1}{2}\left[ b_2 \left( n^{0}_k + n^{0}_{k+1}\right) + d_2 \left( n^{+}_{k-1} + n^{+}_{k+2}\right) \right] \left( n^{+}_k - n^{+}_{k+1}\right) \end{aligned}$$

and

$$\begin{aligned} \frac{1}{2}b_2 \left( n^{+}_k + n^{+}_{k+1}\right) \left( n^{0}_k - n^{0}_{k+1}\right) \end{aligned}$$

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 fugacityFootnote 5\(z_{\lambda }\) parametrizes the density

$$\begin{aligned} \rho _{\lambda } = \sum _{n=1}^{m} n\cdot \mathrm{Prob}\left[ \, {n_k^{\lambda } = n}\, \right] = m \frac{z_{\lambda }}{1+z_{\lambda }} \end{aligned}$$
(8)

on lane \(\lambda \) obtained from (7). The total particle number \(N_{\lambda }\) 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

$$\begin{aligned} K_{\lambda \mu } = z_\mu \frac{\partial \rho _\lambda }{\partial z_\mu } = \kappa _{\lambda } \delta _{\lambda \mu } \end{aligned}$$
(9)

with

$$\begin{aligned} \kappa _{\lambda }=\rho _{\lambda }/(1+z_{\lambda })=\rho _{\lambda }(1-\rho _{\lambda }/m). \end{aligned}$$
(10)

Expectations of time-independent functions f in the stationary distribution are denoted by \(\langle \, {f}\, \rangle \).

2.2.2 Currents

The so-called instantaneous currents for the bond \((k,k+1)\) in lane \(\lambda \) are the functions

$$\begin{aligned} j^{\lambda }_k(t) := r^{\lambda }_{k}(t) - \ell ^{\lambda }_{k+1}(t) \end{aligned}$$
(11)

of the time-dependent occupation numbers \(n^{\lambda }_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 time-dependent expected local density \(\rho ^{\lambda }_k(t):= \langle \, {n^{\lambda }_k(t)}\, \rangle _{P_0}\) as

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\rho ^{\lambda }_k(t) = \langle \, {j^{\lambda }_{k-1}(t)}\, \rangle _{P_0} - \langle \, {j^{\lambda }_k(t)}\, \rangle _{P_0}. \end{aligned}$$
(12)

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_\lambda \), which are translation invariant, one has \(j_{\lambda } = \langle \, {j^{\lambda }_k}\, \rangle \). The product structure of the stationary distribution yields

$$\begin{aligned} j_{-}= & {} - \left( a + b \rho _0 + d \rho _{-} \right) \rho _{-} \left( 1-\frac{1}{m} \rho _{-}\right) \end{aligned}$$
(13)
$$\begin{aligned} j_0= & {} b\left( \rho _{+}-\rho _{-}\right) \rho _0 \left( 1-\frac{1}{m} \rho _{0}\right) \end{aligned}$$
(14)
$$\begin{aligned} j_{+}= & {} \left( a + b \rho _0 + d \rho _{+} \right) \rho _{+} \left( 1-\frac{1}{m} \rho _{+}\right) \end{aligned}$$
(15)

with

$$\begin{aligned} a := a_2 - a_1, \quad b := b_2 - b_1, \quad d := d_2 - d_1. \end{aligned}$$
(16)

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

$$\begin{aligned} J_{\lambda \mu } = \frac{\partial j_\lambda }{\partial \rho _\mu }. \end{aligned}$$
(17)

We are specifically interested in the symmetric case \(\rho _+=\rho _-=:\rho \). In this case the current Jacobian takes the form

$$\begin{aligned} J = \left( \begin{array}{ccc} - \delta &{} - b\kappa &{} 0 \\ - b\kappa _0 &{} 0 &{} b\kappa _0 \\ 0 &{} b\kappa &{} \delta \end{array} \right) \end{aligned}$$
(18)

with

$$\begin{aligned}&\delta = a + b \rho _0 + 2\left( d - \frac{a+b\rho _0}{m}\right) \rho - \frac{3d}{m} \rho ^2, \end{aligned}$$
(19)
$$\begin{aligned}&\kappa = \rho \left( 1- \frac{\rho }{m}\right) , \quad \kappa _0 = \rho _0 \left( 1- \frac{\rho _0}{m}\right) . \end{aligned}$$
(20)

The Hessians \(H^\nu \) are defined by the matrix elements

$$\begin{aligned} H^\nu _{\lambda \mu } = \frac{\partial ^2 j^\nu }{\partial \rho _\lambda \partial \rho _\mu }. \end{aligned}$$
(21)

Defining

$$\begin{aligned} x = d-\frac{a+b\rho _0 + 3d \rho }{m}, \quad y = b\left( 1 - \frac{2\rho }{m}\right) , \quad z = b\left( 1 - \frac{2\rho _0}{m}\right) \end{aligned}$$
(22)

one has for equal densities \(\rho _\pm =\rho \)

$$\begin{aligned} H^{-} = - \left( \begin{array}{ccc} 2x &{} y &{} 0 \\ y &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{array} \right) , \,\, H^0 = z \left( \begin{array}{ccc} 0 &{} -1 &{} 0 \\ -1 &{} 0 &{} 1 \\ 0 &{} 1 &{} 0 \end{array} \right) , \,\, H^{+} = \left( \begin{array}{ccc} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} y \\ 0 &{} y &{} 2x \end{array} \right) . \end{aligned}$$
(23)

By construction, the Hessians are symmetric.

2.2.3 Dynamical Structure Function

The (real-space) dynamical structure function is the matrix \(\bar{S}_k(t)\) with matrix elements

$$\begin{aligned} \bar{S}^{\lambda \mu }_k(t) = \langle \, {(n^{\lambda }_k(t)-\rho _\lambda ) (n^{\mu }_0(0)-\rho _\mu )}\, \rangle . \end{aligned}$$
(24)

Here the brackets without subscript denote an average over histories of the process and the stationary distribution. Because of the factorized stationary distribution one has

$$\begin{aligned} \bar{S}^{\lambda \mu }_k(0) = \kappa _\lambda \delta _{\lambda ,\mu } \delta _{k,0} \end{aligned}$$
(25)

and, due to particle number conservation,

$$\begin{aligned} \sum _{k=0}^{L-1} \bar{S}_k\left( t\right) = K \quad \forall t \ge 0 \end{aligned}$$
(26)

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

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\sum _k k \bar{S}_k\left( t\right) = JK \quad \forall t \ge 0 \end{aligned}$$
(27)

where J is the current Jacobian (17). The matrix product JK is symmetric as has been proved in [16], 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 =: \text {diag}\left( v_-,v_0,v_+\right) \) with eigenvalues \(v_\alpha \). The three eigenvectors of J define the real-space eigenmodes

$$\begin{aligned} \phi _k^\alpha := \sum _\mu R_{\alpha \mu } (n_k^\mu - \rho _\mu ) \end{aligned}$$
(28)

which are fluctuation fields that travel with velocities \(v_\alpha \). We normalize the transformation R by

$$\begin{aligned} R K R^{T} = \mathbb {1}\end{aligned}$$
(29)

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

$$\begin{aligned} S^{\alpha \beta }_{k}\left( t\right) = \left( R\bar{S}_{k}\left( t\right) R^{T}\right) _{\alpha \beta } = \langle \, {\phi ^\alpha _k(t) \phi ^\beta _0}\, \rangle . \end{aligned}$$
(30)

Thus, with the normalization (29) the transformed dynamical structure function S satisfies

$$\begin{aligned} \sum _{k} S_k\left( t\right) = \mathbb {1}, \quad \frac{\mathrm {d}}{\mathrm {d}t}\sum _k k S_k\left( t\right) = V. \end{aligned}$$
(31)

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\ge 0\).

The lattice Fourier transform \(\hat{\bar{S}}(p,t)\) of the dynamical structure function is defined by

$$\begin{aligned} \hat{\bar{S}}^{\lambda \mu }(p,t) = \sum _{k=0}^{L-1} \mathrm {e}^{-2\pi i pk} \bar{S}^{\lambda \mu }_k(t), \end{aligned}$$
(32)

with p of the form \(p=l/L\) with l integer. We can also consider the correlation between the Fourier transformed density fields

$$\begin{aligned} \hat{n}^{\lambda }(p,t) = \sum _{k=0}^{L-1} \mathrm {e}^{-2\pi i pk} (n^{\lambda }_k(t)-\rho _\lambda ), \end{aligned}$$
(33)

from which one constructs

$$\begin{aligned} \tilde{\bar{S}}^{\lambda \mu }(p,q,t) = \langle \, {\hat{n}^{\lambda }(p,t)\hat{n}^{\mu }(q,0)}\, \rangle . \end{aligned}$$
(34)

Exploiting translation invariance one finds \(\tilde{\bar{S}}(p,q,t) = L \hat{\bar{S}}(p,t) \delta _{p+q,0}\) and therefore

$$\begin{aligned} \hat{\bar{S}}(p,t) = \frac{1}{L} \tilde{\bar{S}}(p,-p,t). \end{aligned}$$
(35)

Fourier eigenmodes are calculated similar to (28) as

$$\begin{aligned} \hat{\phi }^\alpha (p,t) := \sum _\mu R_{\alpha \mu } \hat{n}^\mu (p,t) \end{aligned}$$
(36)

and from (35) one concludes

$$\begin{aligned} \hat{S}^{\alpha \beta }\left( p, t\right) = \left( R\hat{\bar{S}}\left( p, t\right) R^{T}\right) _{\alpha \beta } = \frac{1}{L}\langle \, {\hat{\phi }^\alpha (p,t) \hat{\phi }^\beta (-p,0)}\, \rangle \end{aligned}$$
(37)

satisfying

$$\begin{aligned} \hat{S}(0,t) = \mathbb {1}. \end{aligned}$$
(38)

2.3 Monte Carlo Simulations for the Structure Function

2.3.1 Monte Carlo Algorithm

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\). Monte Carlo simulations are performed for a periodic system of very large length \(L=10^6\) with fixed particle numbers \(N_\lambda \). Hereby, large periodic systems allow to suppress finite size effects [17] while the translational invariance increases the quality of the estimator. Further, fixing the particle numbers allow a better control of finite size effects which are of order \(\mathcal {O}(1/L)\) instead of \(\mathcal {O}(1/\sqrt{L})\). The densities are then given by \(\rho _\lambda = N_\lambda /L\). Notice that on switching to a system with \(m=1\) and fixed particle numbers, the stationary distribution is modified and each configuration with \(N_\lambda \) occupied sites on lane \(\lambda \) becomes equally likely to be observed, as follows from Eqs. (1-6) in a similar way as for the grand canonical case. 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

$$\begin{aligned} \bar{\sigma }_{L,k}^{\lambda \mu }\left( M,\tau ,t\right)= & {} \frac{1}{M}\sum _{j=1}^{M}\frac{1}{L}\sum _{l=1}^{L}\left[ n_{l}^{\lambda }\left( j\tau \right) n_{l+k}^{\mu } \left( j\tau +t\right) -\rho _{\lambda }\rho _{\mu }\right] \end{aligned}$$
(39)

where \(l+k\) has to be taken modulo L. Further, \(\tau \) is the time between measurements and M is the total number of measurements. This allows a further variance reduction of the Monte Carlo estimator. The optimal value of \(\tau \) has to be chosen by minimizing the estimators variance for fixed computation time. For too short \(\tau \) subsequent measurements contain too little new information to compensate for the computational overhead and for too long \(\tau \) valuable information is lost. The optimal value for \(\tau \) is determined empirically, as it depends on the computation code and system parameters. In order to compute the structure functions we generate P independent initial configurations yielding Monte Carlo estimators \(\bar{\sigma }^{(p)}_{L,k}\) for each initial configuration. Averaging over the initial configurations then yields the numerical structure function

$$\begin{aligned} S^{\text {MC}}_{L,k}\left( P,M,\tau ,t\right) = \frac{1}{P}\sum _{p=1}^{P} R\bar{\sigma }_{L,k}^{\left( p\right) }\left( M,\tau ,t\right) R^{T} \end{aligned}$$
(40)

which depends on the simulation parameters \(L,P,M,\tau \), 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 \(\sim 1/\sqrt{PM}\). Our values for these simulation parameters are given in Sec. (4).

2.3.2 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

$$\begin{aligned} \mathbb {L}= [v_lt-c_lt^{1/z_l}, v_ht+c_ht^{1/z_h}], \end{aligned}$$
(41)

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_{\lambda }=L \rho _{\lambda }\) introduce long-range 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

$$\begin{aligned} \langle \, {n^\lambda _k n^\mu _0}\, \rangle _{N_+,N_0,N_-} - \rho _{\lambda }\rho _{\mu } = - \frac{1}{L-1} \rho _{\lambda } (1-\rho _{\lambda }) \delta _{\lambda \mu } \quad \text{ for } k \ne 0 \end{aligned}$$
(42)

for fixed particle numbers \(N_{\lambda }\). Thus in a canonical ensemble at time \(t>0\) and k outside the light cone one has, within numerical measurement precision,

$$\begin{aligned} \bar{S}_{\text {can},k}^{\lambda \mu }\left( t\right) = - \frac{\rho _{\lambda }\left( 1-\rho _{\lambda }\right) }{L-1} \delta _{\lambda \mu } . \end{aligned}$$
(43)

Transforming to eigenmodes one gets

$$\begin{aligned} S_{\text {can},k}^{\alpha \beta }\left( t\right) = -\frac{1}{L-1} \delta _{\alpha \beta }. \end{aligned}$$
(44)

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

$$\begin{aligned} \tilde{\sigma }_{L,k}^{\lambda \mu }\left( M,\tau ,t\right) = \bar{\sigma }_{L,k}^{\lambda \mu }\left( M,\tau ,t\right) + \frac{\rho _{\lambda }\left( 1-\rho _{\lambda }\right) }{L-1} \delta _{\lambda \mu }. \end{aligned}$$
(45)

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_\lambda }\), 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.

2.3.3 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:

$$\begin{aligned} a_{0}= & {} -\min \left( 0,b\right) \end{aligned}$$
(46)
$$\begin{aligned} a_{1}= & {} -\min \left( 0,b\right) -\min \left( 0,d\right) \end{aligned}$$
(47)
$$\begin{aligned} a_{2}= & {} a+a_{1} \end{aligned}$$
(48)
$$\begin{aligned} b_{1}= & {} 0 \end{aligned}$$
(49)
$$\begin{aligned} b_{2}= & {} b \end{aligned}$$
(50)
$$\begin{aligned} d_{1}= & {} 0 \end{aligned}$$
(51)
$$\begin{aligned} d_{2}= & {} d \end{aligned}$$
(52)

When choosing the rates in this way, the hopping rates simplify to

$$\begin{aligned} l_{k}^{-}= & {} \left\{ a_{2}+\frac{b}{2}\left( n_{k}^{0}+n_{k-1}^{0}\right) +\frac{d}{2}\left( n_{k+1}^{-}+n_{k-2}^{-}\right) \right\} n_{k}^{-}\left( 1-n_{k-1}^{-}\right) \end{aligned}$$
(53)
$$\begin{aligned} r_{k}^{-}= & {} a_{1}n_{k}^{-}\left( 1-n_{k+1}^{-}\right) \end{aligned}$$
(54)
$$\begin{aligned} l_{k}^{0}= & {} \left\{ a_{0}+\frac{b}{2}\left( n_{k}^{-}+n_{k-1}^{-}\right) \right\} n_{k}^{0}\left( 1-n_{k-1}^{0}\right) \end{aligned}$$
(55)
$$\begin{aligned} r_{k}^{0}= & {} \left\{ a_{0}+\frac{b}{2}\left( n_{k}^{+}+n_{k+1}^{+}\right) \right\} n_{k}^{0}\left( 1-n_{k+1}^{0}\right) \end{aligned}$$
(56)
$$\begin{aligned} l_{k}^{+}= & {} a_{1}n_{k}^{+}\left( 1-n_{k-1}^{+}\right) \end{aligned}$$
(57)
$$\begin{aligned} r_{k}^{+}= & {} \left\{ a_{2}+\frac{b}{2}\left( n_{k}^{0}+n_{k+1}^{0}\right) +\frac{d}{2}\left( n_{k-1}^{+}+n_{k+2}^{+}\right) \right\} n_{k}^{+}\left( 1-n_{k+1}^{+}\right) \end{aligned}$$
(58)

and result in less storage accesses during the simulation. Additionally, for \(a\ge 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^+\); \(\rho _+=\rho _-\) 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_{\pm }=\pm 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.

3 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_\alpha \) 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_\alpha (x,t)\) with continuous space coordinate x.

We follow our previous work [8] and use the mode coupling equations

$$\begin{aligned} \partial _t S_{\alpha }(x,t) = -v_\alpha \partial _x S_{\alpha }(x,t)+ D_{\alpha } \partial _x^2S_{\alpha }(x,t) + \int _0^t \mathrm {d}s \int _{-\infty }^\infty \mathrm {d}y\, S_{\alpha }(x-y,t-s) M_{\alpha \alpha }(y,s) \end{aligned}$$
(59)

with memory term

$$\begin{aligned} M_{\alpha \alpha }(y,s) = 2 \sum _{\beta ,\gamma } (G^{\alpha }_{\beta \gamma })^2 \partial _y^2\left[ S_\beta (y,s)S_\gamma (y,s) \right] \end{aligned}$$
(60)

to predict the large scale-scale behaviour of \(S_\alpha (x,t)\). Here the \(G^{\alpha }\) are the mode coupling matrices obtained from the Hessians \(H^\alpha \) (21) by the transformation

$$\begin{aligned} G^{\alpha } = \frac{1}{2}\sum _\beta R_{\alpha \beta } \left( R^{-1}\right) ^T H^\beta R^{-1}. \end{aligned}$$
(61)

The diffusion coefficients \(D_\alpha \) turn out to be immaterial for the asymptotic theoretical predictions that arise from (59). By construction, the mode coupling matrices are symmetric and related to the mode coupling matrices \(W^{\alpha }\) of [1] by \(W^{\alpha } = 2 G^{\alpha }\).

3.1 Dynamical Scaling

In the scaling limit \(t\rightarrow \infty \) and \(x\rightarrow \infty \) with constant scaling variable \(u_\alpha = E_\alpha (x-v_\alpha t)t^{-1/z_\alpha }\) with a non-universal scale parameter \(E_\alpha \), non-universal mode velocity \(v_\alpha \) and universal dynamical exponent \(z_\alpha \) the mode coupling equations (59) with memory term (60) can be solved exactly for all modes \(\alpha \) and any number of conservation laws [8]. Thus one can determine from them self-consistently for each mode the dynamical exponent \(z_\alpha \) and the corresponding scaling form of the dynamical structure function \(S_\alpha (x,t) = (E_\alpha t)^{-1/z_\alpha } s_\alpha (u)\).

We note that generally, in the case of distinct mode velocities (\(v_\alpha \not = v_\beta \) for all \(\alpha \not = \beta \)) and the absence of purely diffusive modes \(\delta \) (for which \(G^\delta _{\alpha \alpha } = 0\) for all modes \(\alpha \)), any mode \(\alpha \) with non-vanishing self-coupling coefficient \(G^\alpha _{\alpha \alpha }\) is expected to be in the KPZ universality class with dynamical exponent \(z_\alpha = 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_\alpha (u) = f_{PS} (u)\) [5].

For Hamiltonian dynamics this argument applies to the two sound modes with velocities \(v_\pm = \pm v\) so that their asymptotic behavior is expected to become

$$\begin{aligned} S_{\pm }\left( x,t\right) \simeq \left( \lambda t\right) ^{-2/3}\cdot f_{\text {PS}}\left( \left( x\mp vt\right) \cdot \left( \lambda t\right) ^{-2/3}\right) \end{aligned}$$
(62)

with

$$\begin{aligned} \lambda = 2^{3/2}\left| G_{++}^{+}\right| = 2^{3/2}\left| G_{--}^{-}\right| \end{aligned}$$
(63)

and speed \(v = |v_\pm | \ne 0\). The scaling function \(f_{\text {PS}}\) is known exactly and satisfies

$$\begin{aligned} f_{\text {PS}}\left( -x\right)= & {} f_{\text {PS}}\left( x\right) \end{aligned}$$
(64)
$$\begin{aligned} \max f_\text {PS}(x)= & {} f(0)=0.542461\ldots \end{aligned}$$
(65)
$$\begin{aligned} \intop _{-\infty }^{\infty }f_{\text {PS}}\left( x\right) \text {d}x= & {} 1 \end{aligned}$$
(66)
$$\begin{aligned} \intop _{-\infty }^{\infty }\left( f_{\text {PS}}\left( x\right) \right) ^{2} \text {d}x= & {} 0.389813\ldots =:c_{\text {PS}} . \end{aligned}$$
(67)

There is no expression in closed form for \(f_\mathrm {PS}\) which, however, has been calculated with high precision and is tabulated in [6]. Using this data a more precise calculation of \(c_\mathrm {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

$$\begin{aligned} S_{0}(x,t) = \frac{1}{2\pi }\intop _{-\infty }^{\infty }\exp \left( -E_{0}t\left| p\right| ^{5/3}\right) \mathrm {e}^{ipx}\text {d}p \end{aligned}$$
(68)

with [8]

$$\begin{aligned} E_{0} =\frac{a_h}{2} v^{-1/3} \left| G_{++}^{0}\right| ^{2} \left| G_{++}^{+}\right| ^{-2/3}. \end{aligned}$$
(69)

and

$$\begin{aligned} a_{h} = 4\Gamma \left( \frac{1}{3}\right) \sin \left( \frac{\pi }{3}\right) c_{\text {PS}} = 3.6175\ldots \end{aligned}$$
(70)

Here \(\Gamma (\cdot )\) 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\)Footnote 6

3.2 Eigenmodes and Mode Coupling Matrices of the Three-Lane Model

The scenario described above is expected for generic Hamiltonian dynamics with short-range interactions. Specifically, one requires in the frame of vanishing center-of-mass velocity

$$\begin{aligned} v_0 = 0, \quad v_\pm = \pm v \end{aligned}$$
(71)

and the mode coupling symmetry [1]

$$\begin{aligned} G^\gamma _{\alpha \alpha } = - G^{-\gamma }_{-\alpha -\alpha }. \end{aligned}$$
(72)

with \(G^-_{--} \ne 0\) and \(G^0_{--}\ne 0\). Hence, to verify these requirements, one needs to compute all diagonal mode coupling coefficients \(G^\alpha _{\beta \beta }\) for the three-lane exclusion process.

Indeed, the characteristic polynomial of the Jacobian (18) yields the eigenvalues

$$\begin{aligned} v_0 = 0, \quad v_{\pm } = \pm v \end{aligned}$$
(73)

with the strictly positive constant

$$\begin{aligned} v = \sqrt{\delta ^2 + 2 b^2 \kappa \kappa _0} . \end{aligned}$$
(74)

We choose the diagonalizing matrix R such that

$$\begin{aligned} R J R^{-1} = \text{ diag } (-v,0,v). \end{aligned}$$
(75)

Together with the normalization (29) this yields

$$\begin{aligned} R = \frac{1}{\xi } \left( \begin{array}{ccc} \frac{b\sqrt{\kappa _0}}{v-\delta } &{} \frac{1}{\sqrt{\kappa _0}} &{} -\frac{b\sqrt{\kappa _0}}{v+\delta } \\ \frac{1}{\sqrt{\kappa }} &{} -\frac{\delta }{b\sqrt{\kappa }\kappa _0} &{} \frac{1}{\sqrt{\kappa }} \\ -\frac{b\sqrt{\kappa _0}}{v+\delta } &{} \frac{1}{\sqrt{\kappa _0}} &{} \frac{b\sqrt{\kappa _0}}{v-\delta } \end{array} \right) , \, R^{-1} = \frac{1}{\xi } \left( \begin{array}{ccc} \frac{b\sqrt{\kappa _0}\kappa }{v-\delta } &{} \sqrt{\kappa } &{} -\frac{b\sqrt{\kappa _0}\kappa }{v+\delta } \\ \sqrt{\kappa _0} &{} -\frac{\delta }{b\sqrt{\kappa }} &{} \sqrt{\kappa _0} \\ -\frac{b\sqrt{\kappa _0}\kappa }{v+\delta } &{} \sqrt{\kappa } &{} \frac{b\sqrt{\kappa _0}\kappa }{v-\delta } \end{array} \right) \end{aligned}$$
(76)

where

$$\begin{aligned} \xi = \frac{v}{b\sqrt{\kappa \kappa _0}} = \sqrt{2+\frac{\delta ^2}{b^2\kappa \kappa _0}} > 0. \end{aligned}$$
(77)

Notice the following symmetry

$$\begin{aligned} R_{\alpha \beta } = R_{-\alpha -\beta }, \quad (R^{-1})_{\alpha \beta } = (R^{-1})_{-\alpha -\beta } \end{aligned}$$
(78)

which will play a role below.

With the sum

$$\begin{aligned} \tilde{H}^{\gamma }= & {} \sum _{\lambda } R_{\gamma \lambda } H^{\lambda } \end{aligned}$$
(79)
$$\begin{aligned}= & {} \left( \begin{array}{ccc} - 2xR_{\gamma -} &{} -(zR_{\gamma 0}+yR_{\gamma -}) &{} 0 \\ -(zR_{\gamma 0}+yR_{\gamma -}) &{} 0 &{} zR_{\gamma 0}+yR_{\gamma +} \\ 0 &{} zR_{\gamma 0}+yR_{\gamma +} &{} 2xR_{\gamma +} \end{array} \right) . \end{aligned}$$
(80)

of the Hessians (23) we find

$$\begin{aligned} G^\gamma _{\alpha \alpha }= & {} \frac{1}{2} \sum _{\mu \nu } (R^{-1})_{\mu \alpha } \tilde{H}^\gamma _{\mu \nu } (R^{-1})_{\nu \alpha } \end{aligned}$$
(81)
$$\begin{aligned}= & {} \frac{1}{2} \left[ \tilde{H}^\gamma _{--} \left( (R^{-1})_{-\alpha }\right) ^2 + \tilde{H}^\gamma _{++} \left( (R^{-1})_{+\alpha }\right) ^2 \right] \nonumber \\&+ \tilde{H}^\gamma _{0-} (R^{-1})_{0\alpha } (R^{-1})_{-\alpha } + \tilde{H}^\gamma _{+0} (R^{-1})_{+\alpha } (R^{-1})_{0\alpha } \end{aligned}$$
(82)
$$\begin{aligned}= & {} x \left[ R_{\gamma +} \left( (R^{-1})_{+\alpha }\right) ^2 - R_{\gamma -} \left( (R^{-1})_{-\alpha }\right) ^2 \right] \nonumber \\&+ (R^{-1})_{0\alpha } \left[ (z R_{\gamma 0} + yR_{\gamma +}) (R^{-1})_{+\alpha } - (z R_{\gamma 0} + y R_{\gamma -}) (R^{-1})_{-\alpha } \right] . \end{aligned}$$
(83)

The mode coupling symmetry (72) is indeed satisfied without any fine-tuning of parameters, as can be seen as follows.

The Hessians (23) exhibit the symmetry

$$\begin{aligned} H^\gamma _{\alpha \beta } = - H^{-\gamma }_{-\alpha -\beta }. \end{aligned}$$
(84)

which induces \(\tilde{H}^\gamma _{\alpha \beta } = - \tilde{H}^{-\gamma }_{-\alpha -\beta }\). Therefore, with \((R^{-1})^T_{\alpha \beta }=R^{-1}_{\beta \alpha }\) one gets

$$\begin{aligned} G^\gamma _{\alpha \alpha }= & {} \frac{1}{2} \sum _{\mu \nu } (R^{-1})_{\mu \alpha } \tilde{H}^\gamma _{\mu \nu } (R^{-1})_{\nu \alpha } \end{aligned}$$
(85)
$$\begin{aligned}= & {} - \frac{1}{2} \sum _{\mu \nu } (R^{-1})_{\mu \alpha } \tilde{H}^{-\gamma }_{-\mu ,-\nu } (R^{-1})_{\nu \alpha } \end{aligned}$$
(86)
$$\begin{aligned}= & {} - \frac{1}{2} \sum _{\mu \nu } (R^{-1})_{-\mu ,\alpha } \tilde{H}^{-\gamma }_{\mu \nu } (R^{-1})_{-\nu ,\alpha } \end{aligned}$$
(87)
$$\begin{aligned}= & {} - \frac{1}{2} \sum _{\mu \nu } (R^{-1})_{\mu ,-\alpha } \tilde{H}^{-\gamma }_{\mu \nu } (R^{-1})_{\nu ,-\alpha } \end{aligned}$$
(88)
$$\begin{aligned}= & {} - G^{-\gamma }_{-\alpha ,-\alpha } \end{aligned}$$
(89)

One finds for the independent non-vanishing coefficients

$$\begin{aligned} G^-_{--}= & {} - \frac{\kappa \kappa _0}{v^2} \sqrt{\kappa } \left[ x \left( \frac{\delta ^2}{\kappa \kappa _0} + \frac{b^2}{2}\right) + \frac{y b \delta }{\kappa } + z b^2 \right] \end{aligned}$$
(90)
$$\begin{aligned} G^-_{00}= & {} - \frac{\kappa \kappa _0}{v^2} b \sqrt{\kappa } \left( x b - \frac{y\delta }{\kappa }\right) \end{aligned}$$
(91)
$$\begin{aligned} G^-_{++}= & {} - \frac{\kappa \kappa _0}{v^2} b^2 \sqrt{\kappa } \left( \frac{x}{2} - z \right) \end{aligned}$$
(92)
$$\begin{aligned} G^0_{--}= & {} - \frac{\kappa \kappa _0}{v^2} b \sqrt{\kappa _0} \left( (x-z)\frac{\delta }{\kappa _0} + y b \right) . \end{aligned}$$
(93)

Some steps of the lengthy computation are presented in the appendix.

4 Monte Carlo Results for the Structure-Function

4.1 Data Fit

The prediction for the heat mode and the exact Prähofer-Spohn scaling function involve the non-universal scale parameters \(\lambda \) (given in (63)) and \(E_0\) (given in (69)). In order to allow for a comparison with Monte-Carlo results we fit the scale parameters using two different approaches handling finite time effects.

4.1.1 Time Dependent Scale Parameter Using a \(L_{1}\)-Distance Fit

Analogously to [13] we define a time dependent scale parameters \(\Lambda _\alpha (t)\) and study the convergence to their asymptotic theoretical value. \(\Lambda _\alpha \) is given as the minimum of the \(L_{1}\)-distance

$$\begin{aligned} \Lambda _{\alpha }(t) = \underset{\Lambda }{\text {argmin}} \left[ \sum _{k\in \mathcal {K}}\left| S_{\alpha }\left( k,t\right) -\left( \Lambda t\right) ^{-1/z_{\alpha }}f_{\alpha }^{\text {th}} \left( \left( \Lambda t\right) ^{-1/z_{\alpha }}\left( k-v_{\alpha }t\right) \right) \right| \right] \end{aligned}$$
(94)

where \(\mathcal {K}\) is the set of calculated sampling points, \(z_s=3/2\), \(f_s^{\text {th}}(x)=f_\text {PS}(x)\), \(z_h=5/3\) and \(f_h^{\text {th}}(x)=\frac{1}{2\pi }\int _{-\infty }^{\infty } \exp (-p^{5/3})e^{ipx} \mathrm {d}p\).

For a model independent comparison of simulation results and theory we consider analogously to [13] the scale factor coefficients

$$\begin{aligned} a_{s}&=\frac{\lambda }{\left| G_{++}^{+}\right| }=2^{3/2}=2.82842\ldots \end{aligned}$$
(95)
$$\begin{aligned} a_{s}^{\text {fit}}(t)&=\frac{\Lambda _{s}(t)}{\left| G_{++}^{+}\right| } \end{aligned}$$
(96)
$$\begin{aligned} a_{h}&=\frac{2E_{0}}{v^{-1/3}\cdot \left| G_{++}^{0}\right| ^{2}\cdot \left| G_{++}^{+}\right| ^{-2/3}} =4\Gamma \left( \frac{1}{3}\right) \sin \left( \frac{\pi }{3}\right) c_{\text {PS}} =3.61751\ldots \end{aligned}$$
(97)
$$\begin{aligned} a_{h}^{\text {fit}}(t)&=\frac{2\Lambda _{h}(t)}{v^{-1/3}\cdot \left| G_{++}^{0}\right| ^{2}\cdot \left| G_{++}^{+}\right| ^{-2/3}}. \end{aligned}$$
(98)

4.1.2 Scale Parameter Fit with Non-asymptotic Corrections

The strength of finite time effects depends on the chosen model parameter and may disguise the approach to the asymptotic regime within accessible times, leading to erroneous conclusions about scaling factors or universality [14]. The solution of the mode coupling equations (Eqs. (41) and (44) of [8]) quantifies how diffusion and coupling between modes (\(G^{\alpha }_{\beta \beta }\)) contribute to finite time effects. For the KPZ modes, we expect non-asymptotic corrections arising from diffusion and mode coupling to the other KPZ (\(G^{\pm }_{\mp \mp }\)) and heat mode (\(G^{\pm }_{00}\)). On defining the Fourier transform as

$$\begin{aligned} \hat{f}\left( p\right) = \intop _{-\infty }^{\infty }f\left( x\right) \mathrm {e}^{-ipx}\text {d}x \end{aligned}$$
(99)

the Fourier representation of the asymptotic KPZ scaling functions with diffusive and mode coupling corrections with phenomenological constants \(D_{\pm }\) can be approximated as

$$\begin{aligned} \hat{S}_{\pm }(p,t)\simeq e^{\pm i pvt}\hat{f}_{\text {PS}}\left( p(\lambda t)^{2/3}\right) \cdot e^{\left[ D_\pm p^2 + A_1(-iv_p^{\pm \mp })(G^\pm _{\mp \mp })^2 p^{5/3} + A_2(-iv_p^{\pm 0})(G^\pm _{00})^2p^{8/5} \right] t} \end{aligned}$$
(100)

with

$$\begin{aligned} v_p^{\alpha \beta }=(v_\alpha -v_\beta )\text {sgn}(p). \end{aligned}$$
(101)

The heat mode couples only to the KPZ modes and just an additional diffusive term appear in the corrected scaling function

$$\begin{aligned} \intop _{-\infty }^{\infty }S_{0}\left( x,t\right) \mathrm {e}^{-ipx}\text {d}x \simeq \exp \left( -E_{0}t\left| p^{5/3}\right| - D_{0}p^{2}t\right) . \end{aligned}$$
(102)

Here \(D_{0}\) is a further phenomenological diffusion constant. We determine the heat mode scaling factor \(E_0\) by performing a least square fit of the data to Eq. (102) with respect to \(E_0\) and \(D_0\) (see Figs. 5, 10 and related discussion). For the KPZ mode strong diffusive corrections are already known for single-component systems empirically [14]. A fit of the data to Eq. (100) with respect to \(A_{1,2}\) and \(D_+\) involves too many degrees of freedom and is unstable. Thus, we cannot analyze the finite time effects for the KPZ modes in more detail.

4.2 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 \(\rho ,\rho _0\) should be chosen such that the speed v of the KPZ sound modes becomes as large as possible. Furthermore, the scale factors \(\lambda \), \(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. 1.

    Asymptotic KPZ-modes and a 5/3-Lévy-mode with diffusive finite time effects.

  2. 2.

    Asymptotic 5/3-Lévy-mode and KPZ-modes with diffusive finite time effects.

In both examples the convergence to the predicted asymptotic results becomes apparent after subtraction of the diffusive correction.

4.3 Simulation Data

4.3.1 Example 1: Asymptotic KPZ-Modes and a 5/3-Lévy-Mode with Diffusive Finite Time Effects

We choose the model parameters

$$\begin{aligned} \rho= & {} 0.25 \end{aligned}$$
(103)
$$\begin{aligned} \rho _{0}= & {} 0.3 \end{aligned}$$
(104)
$$\begin{aligned} a= & {} 0.5031056 \end{aligned}$$
(105)
$$\begin{aligned} b= & {} 0.4968944 \end{aligned}$$
(106)
$$\begin{aligned} d= & {} 0 \end{aligned}$$
(107)

From this we obtain the sound velocities \(v_{\pm } = \pm 0.35465\) and the mode coupling matrices

$$\begin{aligned} G^{-}= & {} \left( \begin{matrix}0.2139 &{} 0.1269 &{} -0.0255\\ 0.1269 &{} 0.0509 &{} -0.0154\\ -0.0255 &{} -0.0154 &{} 0.0176 \end{matrix}\right) \end{aligned}$$
(108)
$$\begin{aligned} G^{0}= & {} \left( \begin{matrix}0.0854 &{} 0 &{} 0\\ 0 &{} 0 &{} 0\\ 0 &{} 0 &{} -0.0854 \end{matrix}\right) \end{aligned}$$
(109)
$$\begin{aligned} G^{+}= & {} \left( \begin{matrix}-0.0176 &{} 0.0154 &{} 0.0255\\ 0.0154 &{} -0.0509 &{} -0.1269\\ 0.0255 &{} -0.1269 &{} -0.2139 \end{matrix}\right) . \end{aligned}$$
(110)

The theoretical scale parameters are thus given by

$$\begin{aligned} \lambda = 0.605, \quad E_{0} = 5.209\cdot 10^{-2}. \end{aligned}$$
(111)

As Monte Carlo parameters defined in (39) and (40) we choose

$$\begin{aligned} P= & {} 1000 \end{aligned}$$
(112)
$$\begin{aligned} M= & {} 200 \end{aligned}$$
(113)
$$\begin{aligned} \tau= & {} 1000. \end{aligned}$$
(114)

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 \(a_{s,h}\) resp. to their theoretical values. At the largest time \(t=2.33 \cdot 10^5\) the scale parameter \(a_{s}^{\text {fit}}\) of the KPZ sound mode is very close to the exact theoretical value (precision \(\approx 1\%\)). The scale parameter \(a_{h}^{\text {fit}}\) of the heat mode, however, is significantly off (\(\approx 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

$$\begin{aligned} E_{0,\mathrm {fit}}= & {} \left( 1.07\pm 0.02\right) E_{0} \end{aligned}$$
(115)
$$\begin{aligned} D_{0,\mathrm {fit}}= & {} 0.95\pm 0.006 \end{aligned}$$
(116)

where \(D_{0,\mathrm {fit}}\) is the fitted diffusion coefficient of the heat mode. This result corresponds to a deviation of only \(\approx 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.

Fig. 1
figure 1

(Example 1) Comparison of the numerical scale factors \(a^{\text {fit}}_{s,h}\) with their theoretical values \(a_{s,h}\). Asterisks and dashed horizontal line: KPZ sound modes. Bullets and solid horizontal line: Heat mode. The measured scale parameter \(a^{\text {fit}}_s\) for the KPZ modes shows a nice convergence to its theoretical value \(a_s\). The dotted lines between data points are guides to the eye

Full size image

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.

Fig. 2
figure 2

(Example 1) Log-log plot for the value of the KPZ-mode maximum (data points with error bars) versus time compared to the asymptotic solution (62) (full curve). Statistical errors are in order of symbol size

Full size image
Fig. 3
figure 3

(Example 1) Data collapse in position space of the KPZ-mode compared for various times to the asymptotic Prähofer-Spohn scaling-function (62) (Black curve) both in lin-lin (left panel) and lin-log (right panel) representation. Statistical errors are of the order of symbol size

Full size image

For the heat mode the measured dynamical structure function does not collapse in the time range accessible to simulation, see Fig. 4, both for position space and momentum space.

Fig. 4
figure 4

(Example 1) Data collapse of the 5/3-Lévy-mode and scaling-exponent \(z=5/3\) compared to a fitted symmetric 5/3-Lévy-stable distribution described in Fourier representation by Eq. (102) and (69) with \(E_{0,\mathrm {fit}}=\left( 1.07\mp 0.02\right) \cdot E_{0}\) and \(E_{0}= 5.209\cdot 10^{-2}\). Upper panel: Momentum space. Lower panels: Position space in lin-lin (left) and lin-log (right) representation. Statistical errors are in order of symbol size

Full size image

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.

Fig. 5
figure 5

(Example 1) Fit for diffusive finite time effects in the 5/3-Lévy-mode in Fourier space. Data are compared to a fitted symmetric 5/3-stable distribution including diffusive finite time effects given by Eq. (102) with \(D_{0,\mathrm {fit}}=0.095\pm 0.006\), \(E_{0,\mathrm {fit}}=\left( 1.07\mp 0.02\right) \cdot E_{0}\) and \(E_{0}= 5.209\cdot 10^{-2}\)

Full size image

4.3.2 Example 2: Asymptotic 5/3-Lévy-Mode and KPZ-Modes with Finite Time Effects

Here we demonstrate that the presence of diffusive corrections is not a particular feature of the heat mode. They appear also in the sound modes and their impact has to do with the values of the scale parameters and hence with the mode coupling coefficients. To this end we take as model parameters

$$\begin{aligned} \rho= & {} 0.395 \end{aligned}$$
(117)
$$\begin{aligned} \rho _{0}= & {} 0.09 \end{aligned}$$
(118)
$$\begin{aligned} a= & {} 0.522066265 \end{aligned}$$
(119)
$$\begin{aligned} b= & {} 0.477933735 \end{aligned}$$
(120)
$$\begin{aligned} d= & {} 0. \end{aligned}$$
(121)

From this we obtain \(v_{\pm } = \pm 0.1517\) and

$$\begin{aligned} G^{-}= & {} \left( \begin{matrix}0.1487 &{} 0.1342 &{} -0.0318\\ 0.1342 &{} 0.0635 &{} -0.0358\\ -0.0318 &{} -0.0358 &{} 0.064 \end{matrix}\right) \end{aligned}$$
(122)
$$\begin{aligned} G^{0}= & {} \left( \begin{matrix}0.1556 &{} 0 &{} 0\\ 0 &{} 0 &{} 0\\ 0 &{} 0 &{} -0.1556 \end{matrix}\right) \end{aligned}$$
(123)
$$\begin{aligned} G^{+}= & {} \left( \begin{matrix}-0.064 &{} 0.0358 &{} 0.0318\\ 0.0358 &{} -0.0635 &{} -0.1342\\ 0.0318 &{} -0.1342 &{} -0.1487 \end{matrix}\right) \end{aligned}$$
(124)

and

$$\begin{aligned} \lambda = 0.4205, \quad E_{0} = 0.2927. \end{aligned}$$
(125)

Notice, in comparison to the first example, the smaller value of the sound mode self coupling coefficients \(|G^\pm _{\pm \pm }|\) and the much larger value \(|G^0_\pm |\) of the heat mode coupling coefficient which lead to correspondingly different values of \(\lambda \) and \(E_{0}\).

As Monte Carlo parameters we choose

$$\begin{aligned} P= & {} 1000 \end{aligned}$$
(126)
$$\begin{aligned} M= & {} 200 \end{aligned}$$
(127)
$$\begin{aligned} \tau= & {} 1000. \end{aligned}$$
(128)

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 (\(\approx 6\%\)) and the Lévy heat mode (\(\approx 14\%\) from below). The correction to \(a_0\) for the heat mode is non-monotonic.

Fig. 6
figure 6

(Example 2) Comparison of the numerical scale factors \(a^{\text {fit}}_{s,h}\) with their theoretical values \(a_{s,h}\). Asterisks and dashed horizontal line: KPZ sound modes. Bullets and solid horizontal line: Heat mode. The measured scale parameter \(a^{\text {fit}}_s\) for the KPZ modes shows a nice convergence to its theoretical value \(a_s\). The dotted lines between data points are guides to the eye

Full size image

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) but converges to its asymptotic prediction. Also the amplitude at the maximum (Fig. 8) shows deviations from the exact theoretical result that are significantly larger than in the first example (Fig. 2). As discussed in Sect. 4.1.2 the disagreement can be explained by diffusive contributions and a coupling to the other modes.

Fig. 7
figure 7

(Example 2) Data collapse in position space of the KPZ-mode compared for various times to the asymptotic Prähofer-Spohn scaling-function (62) (Black curve) both in lin-lin (left panel) and lin-log (right panel) representation. Statistical errors are of the order of symbol size

Full size image
Fig. 8
figure 8

(Example 2) Log-log plot for the KPZ-mode maximum versus time compared to the asymptotic solution Eq (62). Statistical errors are in order of symbol size

Full size image

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

$$\begin{aligned} E_{0,\mathrm {fit}}= & {} \left( 0.84\mp 0.01\right) E_{0} \end{aligned}$$
(129)
$$\begin{aligned} D_{0,\mathrm {fit}}= & {} 0.015\pm 0.005 \end{aligned}$$
(130)

corresponding to a small diffusion coefficient. The fit with the diffusive correction further improves the data collapse, see Fig. 10.

Fig. 9
figure 9

(Example 2) Data collapse in position space of the 5/3-Lévy-mode and scaling-exponent \(z=5/3\) compared to a fitted symmetric 5/3-Lévy-stable distribution described in Fourier representation by Eq. (102) and (69) with \(E_{0,\mathrm {fit}}=\left( 0.84\mp 0.01\right) \cdot E_{0}\) and \(E_{0}=0.2927\). Upper panel: Momentum space. Lower panels: Position space in lin-lin (left) and lin-log (right) representation. Statistical errors are in order of symbol size

Full size image
Fig. 10
figure 10

(Example 2) Fit for diffusive finite time effects in the 5/3-Lévy-mode in Fourier space. Data are compared to a fitted symmetric 5/3-stable distribution including diffusive finite time effects given by Eq. (102) with \(D_{0,\mathrm {fit}}=0.015\pm 0.005\), \(E_{0,\mathrm {fit}}=\left( 0.84\mp 0.01\right) \cdot E_{0}\) and \(E_{0}=0.2927\)

Full size image

5 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 \(\pm v \ne 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), which is important to the numerical accuracy of the simulation results.

  • Finite time effects: Diffusive corrections generically persist for long times. The sound modes are additionally affected by the coupling to different mode pairs, i.e. \(G^\pm _{\mp \mp }\) and \(G^\pm _{00}\). The impact of these effects 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 exact value of the 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 [18] on the 3/2-Lévy scaling in the case of an anharmonic chain with two conservation laws [12, 19] provide hope that a definite answer can be found.