Kardar-Parisi-Zhang Universality of the Nagel-Schreckenberg Model

TL;DR

This paper proves that the Nagel-Schreckenberg traffic model belongs to the KPZ universality class for all maximum velocities, using hydrodynamics theory and large-scale simulations to confirm the theoretical predictions.

Contribution

It establishes the KPZ universality class for the NaSch model with arbitrary maximum velocities, extending previous results limited to the case v_max=1.

Findings

NaSch model belongs to KPZ class for v_max>1

Simulation results match KPZ asymptotic solutions

Early-time effects influence numerical exponent determination

Abstract

Dynamical universality classes are distinguished by their dynamical exponent $z$ and unique scaling functions encoding space-time asymmetry for, e.g. slow-relaxation modes or the distribution of time-integrated currents. So far the universality class of the Nagel-Schreckenberg (NaSch) model, which is a paradigmatic model for traffic flow on highways, was not known except for the special case $v_{\text{max}}=1$. Here the model corresponds to the TASEP (totally asymmetric simple exclusion process) that is known to belong to the superdiffusive Kardar-Parisi-Zhang (KPZ) class with $z=3/2$. In this paper, we show that the NaSch model also belongs to the KPZ class \cite{KPZ} for general maximum velocities $v_{\text{max}}>1$. Using nonlinear fluctuating hydrodynamics theory we calculate the nonuniversal coefficients, fixing the exact asymptotic solutions for the dynamical structure function and the distribution of time-integrated currents. Performing large-scale Monte-Carlo simulations we show that the simulation results match the exact asymptotic KPZ solutions without any fitting parameter left. Additionally, we find that nonuniversal early-time effects or the choice of initial conditions might have a strong impact on the numerical determination of the dynamical exponent and therefore lead to inconclusive results. We also show that the universality class is not changed by extending the model to a two-lane NaSch model with dynamical lane changing rules.

Equations

$$v_n\to \text{min}\left(v_n+1,v_{\text{max}}\right)$$

$$v_n\to \text{min}\left(v_n,d_n\right)$$

$$v_n\stackrel{p_s}{\to }\text{max}\left(v_n-1,0\right).$$

$$x_n\to x_n+v_n$$

$${\partial }_t\mathtt{u}(x,t)$$

$${\mathcal{B}}^2=2\kappa \mathcal{D},$$

$$\kappa =\int ⟨\mathtt{u}(0,t)\mathtt{u}(x,t)⟩\text{d}x$$

$${\partial }_{t}h=\nu {\partial }_x^{2}h+\frac{\lambda }{2}{\left({\partial }_{x}h\right)}^2+\sqrt{D}\xi$$

$$\nu =\mathcal{D},\text{leavevmode}\text{leavevmode}\lambda =-j^{\prime \prime }(\rho )\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\sqrt{D}=\mathcal{B}.$$

$$S(x,t)$$

$$\mathcal{P}(J,t)\simeq (\Gamma t{)}^{-\frac{1}{2z}}{F}_{\mathrm{BR}}\left(-J\cdot (\Gamma t{)}^{-\frac{1}{2z}}\right)$$

$${\mathtt{J}}_t$$

$$E$$

$$\Gamma$$

$$\kappa$$

$$=$$

$$S(x,t)=⟨n_{x,t}{n}_{0,0}⟩-{\rho }^2$$

$${\mathtt{j}}_{x,t}=\sum _{x^{\prime }=x+1}^{x+v_{\max }}\sum _{p=x^{\prime }-v_{x^{\prime },t}}^{x^{\prime }-1}{\delta }_{x,p}$$

$$J_{x,t}=\sum _{s=1}^t[{\mathtt{j}}_{x,s}-j(\rho )]-\sum _{x^{\prime }=0}^{\lfloor j^{\prime }(\rho )t\rfloor }[n_{x+x^{\prime },0}-\rho ].$$

$$\bar{P}\left(\mathbf{n}\right)=\prod _{x=-\infty }^{\infty }P\left(n_x,n_{x+1}\right),$$

$$P(0,0)$$

$$P(1,1)$$

$$P(0,1)$$

$$P(1,0)=\frac{1}{2(1-p_s)}\left[1-\sqrt{1-4(1-p_s)\rho (1-\rho )}\right].$$

$$j(\rho )$$

$$\kappa$$

$$T_{\mathrm{relax}}$$

$${\kappa }^{-1}(Et{)}^{2/3}S\left((Et{)}^{2/3}x+j^{\prime }(\rho )t,t\right)$$

$$d_n<\min \left(v_n+1,v_{\text{max}}\right)\text{leavevmode}\text{AND}\text{leavevmode}{d}_n<d_n^{\text{(a)}}$$

$$b_n>\min (v_n^{\text{(a)}}+1,v_{\text{max}})$$

$$S^{\lambda \mu }(x,t)$$

$$S(x,t)$$

$$\kappa$$

Full text

I Introduction

In statistical physics, nonequilibrium systems are divided into universality classes according to their dynamical behavior. Dynamical universality classes are distinguished by their dynamical exponent $z$ and unique scaling functions encoding space-time asymmetry for, e.g. slow-relaxation modes or the distribution of time-integrated currents. Systems in the same universality class show for large times of order $1\ll t\ll L^{z}$, where $L$ is the length of the system, identical statistical properties, while local interactions are coded in nonuniversal scaling factors. The two most prominent examples are the diffusive class with dynamical exponent $z=2$ and the superdiffusive Kardar-Parisi-Zhang (KPZ) class with $z=3/2$ KPZ ; Halp15 . A generic example for the latter is the totally asymmetric simple exclusion process (TASEP) asep ; Derrida ; Schuetz00 which describes the single-file motion of uni-directionally moving particles on a discrete one-dimensional lattice. Due to exclusion each lattice site can accommodate at most one particle.

The first indirect proof of KPZ-universality in the TASEP and its partially asymmetric generalization (ASEP) came from finite-size scaling analysis of the spectral gap of the Markov generator, using the Bethe ansatz Dhar ; Gwa92 ; Kim95 . These results yield the dynamical exponent $z=3/2$. Since then there has been remarkable progress that has led to a much more detailed understanding of the fluctuations in one dimensional systems using techniques from random matrix theory FP ; Spohn06 ; S07 ; Corwin2012 ; QS2015 . Several exact solutions for models in the KPZ universality class, such as the ASEP and the KPZ equation Johansson2000 ; BR ; Prae02 ; Prae04 ; FS ; TW2009a ; SS2010 ; ACQ2011 , have resulted in explicit expressions for universal distribution functions and correlations of physical quantities in appropriate scaling limits.

For applications to highway traffic the TASEP has been generalized to the Nagel-Schreckenberg (NaSch) model NaSch92 . Here the particles have an internal degree of freedom, called velocity, which determines their hopping range. The velocity changes dynamically and is limited by a maximum value $v_{\text{max}}$. In contrast to the standard TASEP the NaSch model is defined by a parallel updating scheme which leads to more realistic results. For the special case $v_{\text{max}}=1$ the NaSch model reduces to the TASEP with parallel dynamics.

For more than 20 years one has tried to determine the universality class of the NaSch model. For the case $v_{\text{max}}=1$ it was expected that it also belongs to the KPZ class, i.e. that the use of parallel dynamics does not change the universality class. This was confirmed with random matrix theory in Johansson2000 and was subsequently generalized to other parallel update schemes with $v_{\text{max}}=1$, using determinantal techniques derived from Bethe ansatz Rako05 . However, for general maximum velocities $v_{\text{max}}>1$ the universality class has remained under debate since the internal degree of freedom (i.e. the velocity) might lead to a different universality class and numerical studies were inconclusive CsanyiK95 ; SasvariK97 .

Here we will show that the NaSch model indeed belongs to the KPZ class for all parameter values. Using nonlinear fluctuation hydrodynamics theory we calculate the nonuniversal coefficients fixing the exact asymptotic solutions for the dynamical structure function and the distribution of time-integrated currents. Performing large-scale Monte-Carlo simulations we show that the simulation results match the exact asymptotic KPZ solutions without any fitting parameter left. Additionally, we find that nonuniversal early-time effects, or the choice of initial conditions might have an strong impact on the numerical dynamical exponent determination and therefore lead to inconclusive results. We also show that the universality class is not changed by extending the model to a two-lane NaSch model with dynamical lane changing rules. This implies that neither the use of random-sequential dynamics nor single-file behaviour are essential for the universality.

II Nagel-Schreckenberg Model

The model introduced by Nagel and Schreckenberg (NaSch) NaSch92 is by now regarded as a minimal cellular automation model for traffic flow on highways. It can be viewed as an extension of the TASEP with parallel dynamics to longer-range interactions. However, in contrast to other generalisations of the TASEP which allow for the movement of particles beyond the nearest-neighbour, the NaSch model has a kind of velocity memory controlled by an internal parameter $v$. $v_{n}(t)$ corresponds to the number of cells particle $n$ has moved forward in time step $t$. By the dynamical rules of the model, $v_{n}(t)$ can at most increase by one in the next timestep which mimics the limited acceleration properties of vehicles. This velocity makes the model at the same time more realistic, but also more difficult to analyse. Especially it is not clear whether it has an impact on the dynamical universality class of the model.

The dynamical rules for the NaSch model are given by four steps which are applied to all vehicles at the same time (parallel or synchronous update). The update rule for the $n$-th vehicle is:

Acceleration: If $v_{n}<v_{\text{max}}$, the speed of the $n$-th vehicle is increased by one, but $v_{n}$ remains unaltered if $v_{n}=v_{\text{max}}$, i.e

[TABLE] 2. 2.

Deceleration: If $v_{n}>d_{n}$, ($d_{n}$ is the headway of the $n$-th vehicle) the speed of the $n$-th vehicle is reduced to $d_{n}$, i.e.

[TABLE] 3. 3.

Randomization: If $v_{n}>0$, the speed of the $n$-th vehicle is decreased randomly by unity with probability $p_{s}$, i.e.

[TABLE]

With probability $1-p_{s}$ the velocity of the vehicle remains unchanged. The velocity $v_{n}$ does not change if $v_{n}=0$. 4. 4.

Vehicle movement: Each vehicle is moved forward according to its new velocity determined in 1.-3., i.e.

[TABLE]

These rules are minimal in the sense that the basic features of real highway traffic (e.g. spontaneous jam formation) are no longer reproduced if one rule is left out. Also the order of the rules is essential for realistic behaviour. Throughout this paper, numerical results will refer to an implementation of the NaSch model on a lattice with periodic boundary conditions.

In contrast to the TASEP, so far no closed solution for the stationary state of the NaSch model is known. The reviews SCNBook ; ChowdhurySS00 ; MaerivoetM05 give an overview over known results. Fig. 1 shows the fundamental diagram, i.e. the density-dependence of the stationary current, for different values of $v_{\text{max}}$. For $v_{\text{max}}>1$ the particle-hole symmetry of the TASEP is lost and as a result the function $j(\rho)$ is no longer symmetric around $\rho=1/2$.

III Nonlinear Fluctuating Hydrodynamics (NLFH)

Nonlinear fluctuating hydrodynamics Mori ; Swift ; Das is a powerful phenomenological tool to describe the large-scale behaviour behaviour of fluctuations of conserved quantities in many-body system both in and out of thermal equilibrium Spoh14 . Notably, it captures the large-scale properties of the dynamical structure function in the universality class of the KPZ equation Halp15 as well as an infinite discrete family of other universality classes Popk15b ; Popk16 . In particular, the theory allows for adopting exact results obtained for specific models Prae02 ; Prae04 ; Bern16 to models that are not exactly solvable but are within the same universality class.

In the following we use the non-linear fluctuating hydrodynamic equation for conservative driven diffusive lattice gases with one conserved density and thus establish the connection of the KPZ equation with the NaSch model. Significantly, we will present exact analytic predictions for the dynamical structure function and time integrated current distributions which will serve as a test of the KPZ-universality.

Particle conservation along with local stationarity and slow relaxation of the conserved modes implies that the long-time evolution of the NaSch model at large scales is described in terms of a conservation law $\partial_{t}\varrho(x,t)+\partial_{x}\mathtt{j}(x,t)=0$, where $\varrho(x,t)$ is the coarse-grained local density field, and $\mathtt{j}(x,t)$ is the associated current. Local stationarity ensures that the current $\mathtt{j}(x,t)$ depends on $x$ and $t$ only through the density $\varrho(x,t)$, i.e., one has $\mathtt{j}(x,t)=j(\varrho(x,t))$ with the stationary current-density relation $j(\varrho)$ Kipn99 . Thus $\partial_{x}\mathtt{j}=j^{\prime}(\varrho)\partial_{x}\varrho$ where $j^{\prime}(\varrho)=\mathbf{d}j(\varrho)/(\mathbf{d}\varrho)$. Evidently, $\varrho(x,t)=\rho$ with any constant $\rho$ in the physically permissible range is the stationary solution to this hydrodynamic equation. In this deterministic hydrodynamic description the effects of the noise disappear because of the spatial coarse-graining of the density and the Eulerian scaling of time in which the microscopic space and time scales are rescaled proportionally to a common scaling factor. In the case of lattice gas models these microscopic scales are the lattice constant and the time scale of particle jumps between lattice sites.

As next step one subtracts from the local density field $\varrho(x,t)$ its stationary background $\rho$ to obtain the fluctuation field $\mathtt{u}(x,t)=\varrho(x,t)-\rho$, and expands the current $j(\varrho(x,t))$ in $\mathtt{u}(x,t)$ around the constant $\rho$. To incorporate fluctuations and thus capture the effects of noise arising from the stochastic dynamics and to arrive at the fluctuating hydrodynamic description, a phenomenological diffusion term $\mathcal{D}\partial_{x}\mathtt{u}(x,t)$ and Gaussian white noise $\mathcal{B}\xi(x,t)$ are added to the current. To capture the universal behavior correctly, it suffices to expand the current-density relation up to second order Spoh14 . Possible logarithmic corrections, which arise from higher orders if $j^{\prime\prime}=0$ at some density, are neglected BT_vanB12 ; Delf07 . Thus we arrive at the nonlinear fluctuating hydrodynamics (NLFH) equation

[TABLE]

The noise magnitude $\mathcal{B}$ and the diffusion coefficient $\mathcal{D}$ are related by the fluctuation-dissipation theorem

[TABLE]

where

[TABLE]

is independent of time due to the global particle conservation. This quantity contains information about the system’s space correlations and is a nonequilibrium analogue of the thermodynamic compressibility.

Notice that performing a Galilean transformation $x\rightarrow x-v_{\mathrm{col}}t$ with $v_{\mathrm{col}}\equiv j^{\prime}(\rho)$ removes the drift term $j^{\prime}\left(\rho\right)\mathtt{u}(x,t)$ from the NLFH equation (1) which then by writing $\partial_{x}h(x,t)=-\mathtt{u}(x,t)$ turns into the originally proposed KPZ equation KPZ

[TABLE]

for the surface height $h(x,t)$ with parameters

[TABLE]

For the lattice model the substitution $\partial_{x}h(x,t)=-\mathtt{u}(x,t)$ is motivated by the exact mapping of the TASEP to a discrete surface growth process BARABASI_GROWTH ; HHZ ; KRUG_GROWTH ; KRUG_Book_GROWTH ; TASEP_GROWTH that is known as the single-step model. In Fig. 2 this mapping is generalized to a NaSch scenario resulting in a growing surface with diamonds of different size.

The universal large-scale properties of the KPZ equation are by now well-understood, see Halp15 for a recent review. The dynamical exponent that relates the scaling of space and time variables as $x\sim t^{1/z}$ takes the value $z=3/2$, as opposed to $z=1$ of the deterministic Eulerian scaling or $z=2$ for normal diffusion. Two prominent exact analytic results displaying the space-time symmetry with dynamical exponent $z=3/2$ are the asymptotic limit of the dynamical structure function

[TABLE]

with the Prähofer-Spohn scaling function $f_{\mathrm{PS}}$ Prae04 and the distribution

[TABLE]

of centered time-integrated currents

[TABLE]

with the Baik-Rains scaling function $F_{\mathrm{BR}}(x)=2F_{0}(2x)$ and $F_{0}(\cdot)$ defined in BR . The scaling parameters are

[TABLE]

The exactly known scaling functions $f_{\mathrm{PS}}$ and $F_{\mathrm{BR}}$ are given by solutions of of certain Painlevé II transcendent equations PJF03 , and cannot be expressed in closed form but are tabulated with high precision Prae04_DATA .

In order to check whether the NaSch model is in the KPZ universality class we first calculate the hydrodynamic quantities $\kappa$, $j$, $v_{\mathrm{col}}$ and $j^{\prime\prime}$ exactly for $v_{\text{max}}=1$ and from Monte-Carlo simulations for $v_{\text{max}}>1$. This allows us to fix the analytic predictions for the dynamical structure function and time-integrated current distribution for comparison with numerical data. We stress that the quantities $\kappa$, $v_{\text{col}}$ and $j^{\prime\prime}(\rho)$ are purely stationary quantities that do not require any knowledge about space-time symmetry. Thus, a comparison of simulation results for the dynamical structure function and current statistics with the analytical predictions (6) and (7) serves as a reliable check whether the NaSch model truly belongs the KPZ universality class.

As the discussion so far treats the dynamics as continuous in space and time we need to define for the NaSch model a discrete version of the hydrodynamic quantities $\rho$, $j(\rho)$, $\kappa$, the structure function and time-integrated currents. To this end a configuration $\mathcal{C}(t)=\{n_{x,t},v_{x,t}\}$ of the NaSch model at the end of an update cycle $t$ is expressed by a pair of occupation numbers $n_{x,t}\in\{0,1\}$ and its associated velocities $v_{x,t}\in\{0,\ldots,n_{x,t}\cdot v_{\text{max}}\}$ at site $x$. We limit our simulations to periodic systems of length $L$ with fixed particle density $\rho=\frac{1}{L}\sum_{x=1}^{L}n_{x,t}$. The current-density relation is calculated as $j(\rho)=\rho\left<v\right>$ where $\left<v\right>$ is the stationary average velocity of the cars. For $v_{\text{max}}>1$ the stationary state is unknown and the compressibility is calculated from space correlations as

[TABLE]

where the cutoff $K\ll L/2$ excludes exponentially decaying space correlation which can be neglected within statistical accuracy. With the hydrodynamic quantities at hand the dynamical structure function

[TABLE]

can be measured and compared with the scaling form (6). To define a discrete version of the centered time-integrated current we have to introduce a discrete version of the instantaneous current

[TABLE]

indicating if a particle passes between sites $x$ and $x+1$ during the update from $t-1$ to $t$. Finally, the discrete centered time-integrated current satisfying Eq. (7) is given as

[TABLE]

With these quantities we are in a position to probe in detail the dynamical universality class of the NaSch model.

IV NaSch Model with $v_{\text{max}}=1$

For $v_{\text{max}}=1$ the NaSch model corresponds to the TASEP with parallel dynamics. In this case, the stationary state is exactly known allowing to determine all nonuniversal scaling factors exactly. Therefore, this special case serves as a benchmark to show the convergence towards the asymptotic limit and allows to identify early time contributions. The latter will play a key role to point out that early time contribution might persist longer than expected. This leads to essential insights into finite-time and -size effects for simulations with unknown steady state ($v_{\max}>1$).

Interpreting $p=1-p_{s}$ as the probability that a vehicle will move, we have full equivalence to the TASEP with parallel update rule. Hereby the random-sequential update is included as a limiting case when taking in the limit of $p\rightarrow 0$ (time properly rescaled), while in the limit of $p\rightarrow 1$ the dynamics become deterministic. The exact stationary probability distribution $\bar{P}(\mathbf{n})$ to observe a configuration $\mathbf{n}$ factorizes into a two cluster form

[TABLE]

where $n_{k}\in{0,1}$ is the occupation number at site $k$. Using the Kolomogorov consistency relations

[TABLE]

with $P(1)=\rho$ and $P(0)=1-\rho$ one consequently has to solve the master equation for $P(1,0)$. Expressed in terms of $P(1,0)$ the stationary master equation reduces to a quadratic form and yields

[TABLE]

With the stationary distribution at hand, one calculates the current-density relation and its compressibility as

[TABLE]

fixing the nonuniversal scaling parameter. Knowing these hydrodynamic quantities exactly we are in the position to compare simulation results to the exact asymptotic predictions derived in Sec. III without any free parameter left. Fig. 3 shows a scaling plot with dynamical exponent $z=3/2$ of simulation data obtained for $\rho=1/2$ and various $p_{s}$. Additionally, in Meersoon_Schmidt_17 the parallel update TASEP has been shown to exhibit the Baik-Rains distribution (7) for current fluctuations. Remarkably, the data for the dynamical structure function and current distribution matches the predicted scaling form perfectly, although it is obtained for a model continuous in time and space.

V NaSch Model with $v_{\text{max}}>1$

As mentioned before the stationary state is unknown for $v_{\text{max}}>1$. Thus, the system has to be relaxed before one starts recording observables using the Metropolis sampler. Especially, it has been shown that the KPZ statistics are sensitive to initial conditions and might reveal different scaling functions Meersoon_Schmidt_17 ; Ferrari_Spohn . For a stochastic model the relaxation time $T_{\mathrm{relax}}$ can often be defined through the spectral gap $G_{\mathrm{S}}$ of the time evolution operator SCNBook which depends on the system size as $T_{\mathrm{relax}}\sim G_{\mathrm{S}}\sim L^{z}$, where $z$ is the dynamic exponent. For the ASEP with periodic and open boundaries the spectral gap was calculated exactly using Bethe ansatz methods Dhar ; Gwa92 ; Kim95 ; Gier_Essler .

Since the spectral gap is in general hard to calculate, one might use the dynamical structure function to define an equivalent relaxation time. The dynamical structure function carries the information about the slow relaxation mode and displays the evolution of a perturbation/fluctuation trough the system. The amplitude of the dynamical structure function will decay exponentially instead of $t^{-1/z}$ after the dynamical structure function has been spread over the whole system Prolhac_16 . The width $\sigma$ of the dynamical structure function scales with time as $\sigma\sim t^{1/z}$. In case of the KPZ universality class we define the width as $\sigma(t)\equiv(\sqrt{2\kappa}|j^{\prime\prime}|t)^{2/3}$ covering $\intop_{-0.5}^{0.5}f_{\text{PS}}(x)\text{d}x=50.057\ldots\%$ of the Prähofer Spohn KPZ scaling function, whereas one has $f_{\text{PS}}(0)=2f_{\text{PS}}(\pm 0.88046626\ldots)$. Thus, a lower boundary for a proper relaxation time is given when the structure function width $\sigma$ covers the whole system. Solving $L\lesssim(\sqrt{2\kappa}|j^{\prime\prime}|T_{\mathrm{relax}})^{2/3}$ we derive the relaxation time as

[TABLE]

It turns out that the system’s relaxation is the major computation bottleneck. A propper bound for a minimum required relaxation time allows a significant reduction of computation cost, and ensures relaxation artifacts to be absent. Because, the derived relaxation bound Eq. (23) applies for systems near the stationary state, we introduce a two level relaxation. Thus, we first initialize the system with equally spaced vehicles, velocity $v=v_{\text{max}}$, pre-relax the state according to Eq. (23) and store it in memory. In this way, the chosen initial condition prevents the system of being stuck in a jam which may have a long life-time Nagel_Jam_93 . Second, to generate a new independent state, the relaxed state is loaded and again independently propagated according to Eq. (23). However, there are various ways to initialize the system which may have advantages or disadvantages, depending on the observed quantities. We have tested our data for independence on the initial state by choosing different initial conditions and comparing the observables. Only for $T_{\mathrm{relax}}\gtrsim L^{3/2}/E$ these differences disappear.

Fig. 4 shows the compressibility and the second derivative of the current as function of the density for different values of $v_{\text{max}}$. The behaviour of the compressibility for $v_{\text{max}}>1$ differs clearly from that for $v_{\text{max}}=1$. In the latter case, $\kappa(\rho)$ increases monotonically with increasing density $\rho$ whereas for $v_{\text{max}}>1$ two local extrema exist in the interval $0<\rho<1$. The compressibility is strongly enhanced at higher densities, reflecting the formation of spontaneous traffic jams.

The data for the dynamical structure functions (Fig. 5) for $v_{\text{max}}>1$ and $\rho\approx\rho^{\star}$, where $\rho^{\star}:=\text{argmax}_{\rho\in(0,1)}E$ is the density for which the scaling parameter $E$ becomes maximal, collapse well and show a very good agreement with the asymptotic scaling function (6). The time collapse of the distribution for the time-integrated currents (15) in Fig. 6 shows a nice agreement with the asymptotic Baik-Rains distribution (7).

VI Early time dynamical structure function

NLFH has produced only asymptotic results so far. A full space-time solution of Eq. (1) would allow for a better comparison with simulation data and therefore a better identification of corrections which may arise from higher order corrections to Eq. (1). In this section we will take a closer look on simulation data for the early time dynamical structure function, showing a density-dependent asymmetry that vanishes with time. Non-asymptotic effects might have a strong impact on the identification of universal behavior and may lead to inconclusive results. Therefore, a qualitative understanding of the early time dynamical structure function asymmetry is crucial for the interpretation of simulation data.

In order to easily compare data for different models and parameters, the dynamical structure functions are rescaled to its scaling function as

[TABLE]

Due to the particle-hole symmetry of the TASEP the measured dynamical structure function is symmetric for $\rho=1/2$ and matches the symmetry prediction of the asymptotic solution Eq. (6). However, for densities $\rho\not=1/2$, the early time dynamical structure function shows an asymmetry which vanishes with increasing time. As shown in Fig. 7 the asymmetry is present both for $v_{\text{max}}=1$ (TASEP) and $v_{\text{max}}>1$. Thus the asymmetry is not a special feature of the NaSch model where particle velocities might be interpreted as an internal degree of freedom. Note that fitting the dynamical exponent from the maximum of the dynamical structure function in a non-asymptotic regime will lead to density-dependent and therefore inconclusive results as observed in CsanyiK95 ; SasvariK97 .

The skew of the dynamical structure function that we observe at early times disappears for $\rho\approx\rho^{\star}$. It is negative for $\rho<\rho^{\star}$ and positive for $\rho>\rho^{\star}$ (Fig. 7) and increases with $|\partial_{\rho}E|$. This indicates the role of cubic corrections to Eq. (1) for the full time solution of the dynamical structure function. On the other hand, the distribution of the time-integrated current does not show indications for higher order corrections (Fig. 8).

VII Two-Lane NaSch Model with Dynamical Lane Changes

In order to further understand the relevance of universal behavior for traffic-like models we will now relax the condition of single-file motion. We consider a one-dimensional system with two lanes and dynamical symmetric lane changing rules that allow overtaking on both lanes. For our purpose, we do not need rules that lead to a very realistic simulation of multilane traffic Rickert ; NagelWWS ; SCNBook , but represent only the basic aspects of lane changing.

Generically, a lane change decision is based on two criteria: The incentive criterion which tests for an improvement of the individual traffic situation, e.g. to move forward with their desired velocity, and the safety criterion where each vehicle considers a lane change based on the available backward gap in the desired lane Rickert ; NagelWWS ; SCNBook .

It is natural to split the multi-lane-model update into two substeps: In the first substep vehicles may change lanes and in the second substep vehicles move forward as in the single-lane NaSch model.

The investigated lane change protocol is designed as follows:

• •

Incentive criterion: If the headway $d_{n}$ in front of the $n$-th vehicle on lane $\lambda$ $(\lambda=1,2)$ is too small to travel with the desired speed in the ensuing NaSch update and the headway $d_{n}^{\text{(a)}}$ in the adjacent lane is larger, the vehicle considers a lane change. Otherwise, it stays in its actual lane, i.e.

[TABLE]

• •

Safety criterion: The $n$-th vehicle got a neighboring vehicle on the adjacent lane which might be next to or behind to it. This neighboring vehicle is moving with velocity $v_{n}^{\text{(a)}}$ and $b_{n}$ measures the backward gap. The backward gap is equal to zero if the vehicles are next to each other. To avoid conflicts due to lane changes in the following NaSch update, the backward gaps should be sufficiently large, so that neighboring cars won’t break due to the lane changes, i.e.

[TABLE]

• •

Randomization: If the criteria above are satisfied the vehicle performs a lane change with probability $p_{c}$

Lane changes are performed in parallel. Fig. 9 shows a typical lane change situation.

Note that all vehicles are identical and the system only conserves the overall vehicle density, therefore the system is expected to support a single KPZ-mode. Accounting for the symmetry of the model, the structure function and its hydrodynamic quantities can be defined as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In Fig. 10 a nice agreement between Monte-Carlo simulations and predicted asymptotic KPZ scaling behaviour (6) is shown. In order to reach the asymptotic regime within computation limits a vehicle density $\rho\approx\rho^{\star}$ was used to avoid an early time skew of the dynamical structure function.

The results of this section once again show the robustness of the KPZ universality class. To observe KPZ behavior, single-file motion is not a necessary condition.

VIII Discussion

We have provided strong numerical evidence that the NaSch model of traffic flow belongs to the KPZ universality class for all choices of the parameters $v_{\text{max}}$ and $p$. Previously, this was only known for the special limit $v_{\text{max}}=1$ and random-sequential dynamics where the model corresponds to the TASEP. Previous studies CsanyiK95 ; SasvariK97 were unable to determine the universality class conclusively because of the strong finite-size and finite-time effects (see Fig. 7).

The results presented here provide deeper insights both in the universality of driven diffusive systems and the dynamics of traffic flow. They indicate that neither the updating procedure (parallel update for NaSch vs. random-sequential update for TASEP) nor the internal degree of freedom (i.e. the velocity which is introduced for $v_{\text{max}}>1$) affect the universality. Furthermore, we have shown by considering a multi-lane version of the NaSch model that the universality class is also not changed by deviations from strict single-file motion, i.e. allowing for changes in the particle ordering.

To test for KPZ universality the dynamical structure function and the distribution of time-integrated currents were recorded for various times. Both observations show a nice agreement with the analytical predictions. The dynamical structure function shows a vanishing finite-time asymmetry which is likely to be universal and controlled by cubic corrections of the NLFH theory. The distribution of time-integrated currents do not show any indication for the relevance of cubic or higher corrections. Overall, we have found strong indications that the NLFH theory works properly for systems discrete in space and time (i.e. parallel update). Relevant quantities determining the asymptotic behavior are the current-density relation and the compressibility.

The slow relaxation modes are controlled by the universality class of the system. Monte Carlo simulations showed that, observables recorded in systems with insuffciently relaxed initial states show strong deviations to observables recorded in stationary systems. To overcome effects caused by insufficient relaxation, we have derived a relaxation criterion (23) for single species models that exhibit a nonlinear current-density relation. This criterion yields a precise estimate for the minimal time necessary to reach a state that can be considered stationary in simulations. It would be of interest to understand better the survival of universality in transient regimes of the NaSch model with time-dependent boundaries.

One expects NLFH to hold as well for multi-species models such as traffic models that incorporate cars as well as buses. In this case one expects fluctuations in the eigenmodes of NLFH to be described by explicitly known universal scaling functions, see Popk15b ; Popk16 for the general multi-species case. In CdGHS a more mathematical treatment of KPZ modes has been given for two-species models where NLFH is not postulated and universal distributions have been derived from first principles, confirming NLFH predictions.

Acknowledgements.

JdG would like to thank Tim Garoni for discussions and gratefully acknowledges financial support from the Australian Research Council. JS thanks the University of Melbourne, where parts of this work were done, for hospitality and support. He acknowledges ACEMS and the Bonn-Cologne-Graduate-School for covering travel expenses. This work was supported by Deutsche Forschungsgemeinschaft (DFG) under grant SCHA 636/8-2. AS, JS and GS acknowledge support by the German Excellence Initiative through the University of Cologne Forum ”Classical and Quantum Dynamics of Interacting Particle Systems”.

Appendix A Simulation method

In order to run efficient Monte Carlo simulations, it is recommendable to utilize translational invariance due to periodic boundary conditions and stationarity allowing for ergodic measurements by averaging over space and time, i.e.

[TABLE]

where $\tilde{f}$ and $f$ are Metropolis-Hastings Monte-Carlo estimators evaluating a single stationary Markov Chain $\{n_{x,t}\}_{L}$ of a system with $L$ sites. The evaluation points of interest are the positions $\vec{x}_{*}$ and times $\vec{t}_{*}$, their corresponding ones vectors $\vec{1}=(1,\ldots,1)^{t}$ shift these points in order to make use of the translational invariance and stationarity. The average of $\tilde{f}_{M,\tau}$ over $P\rightarrow\infty$ independent realisations $\{n_{x,t}\}_{L}$ guarantees the convergence to the desired quantity

[TABLE]

Note that, in case of stationarity and translational invariance one has $\mathbb{E}\left(\tilde{f}_{M,\tau}(\vec{x}_{*},\vec{t}_{*})\right)=\mathbb{E}\left(f(\vec{x}_{*},\vec{t}_{*})\right)$, whereas $\tilde{f}_{M,\tau}$ supports a significantly lower variance than $f$ and therefore consumes less computation time to reach the desired accuracy. Further, the time between two ergodic measures $\tau$ may serve as a variance reduction parameter allowing to minimize the uncertainty of the estimator $\tilde{f}_{M(\tau),\tau}$ under fixed computation cost.

E.g. the estimator $\tilde{f}_{M,\tau}$ for the single lane dynamical structure function $S(x,t)=\mathbb{E}(\tilde{f}_{M,\tau})$ (see Eq. (13)) is based on $f\left[\{n_{x,t}\}_{L}\right]((0,x)^{t},(0,t)^{t})=n_{0,0}n_{x,t}-\rho^{2}$.

Independent stationary Markov Chains $\{n_{x,t}\}$ are realised by using independent initial states $\{n_{x,0}\}$ drawn from stationary distribution, and propagated according to the systems update rules with independent sets of random numbers. In case of unknown stationary distribution ($v_{\mathrm{max}}>1$), we use the initial configuration where all vehicles are equally distributed and assigned to their maximum velocity. In order to reach the stationary limit, each configuration is independently propagated with at least $T_{\text{relax}}$ updates (see Eq. (23)).

All pseudo random numbers throughout this paper are generated by the Mersenne Twister generator, implemented in the C++ standard library random.