ar
X
iv
:c
on
d-
m
at
/9
41
20
45
v1
9
D
ec
1
99
4
Discrete stochastic models for traffic flow
M. Schreckenberg1, A. Schadschneider1, K. Nagel2, and N. Ito3,∗
1 Institut fu¨r Theoretische Physik, Universita¨t zu Ko¨ln, D–50937 Ko¨ln, Germany,
schreck@thp.uni-koeln.de and as@thp.uni-koeln.de
2 Zentrum fu¨r Paralleles Rechnen ZPR, Universita¨t zu Ko¨ln, D–50923 Ko¨ln, Germany,
kai@zpr.uni-koeln.de
3 CISC JAERI, Tokai, Ibaraki 319-11, Japan, ito@catalyst.tokai.jaeri.go.jp
Abstract
We investigate a probabilistic cellular automaton model which has been in-
troduced recently. This model describes single-lane traffic flow on a ring and
generalizes the asymmetric exclusion process models. We study the equilib-
rium properties and calculate the so-called fundamental diagrams (flow vs.
density) for parallel dynamics. This is done numerically by computer simu-
lations of the model and by means of an improved mean-field approximation
which takes into account short-range correlations. For cars with maximum
velocity 1 the simplest non-trivial approximation gives the exact result. For
higher velocities the analytical results, obtained by iterated application of the
approximation scheme, are in excellent agreement with the numerical simula-
tions.
∗Present address: Dep. of Applied Physics, Faculty of Engineering, University of Tokyo, Tokyo
113, Japan
1

I. INTRODUCTION
Recently, there has been considerable interest in the investigation of traffic flow using
methods of statistical physics. Independently in [1] and [2] cellular automaton models [3]
for the description of traffic flow have been proposed. Due to their computational simplicity,
lattice gas automata [3] have already successfully been applied to other problems, e.g. the
simulation of fluids [4] (for further applications, see the book of Wolfram [3]).
A similar class of discrete models — which may be interpretated as traffic models or as
models for surface roughening — have also been used for the description of the so-called
asymmetric exclusion process (driven diffusion) [5–13]. Here several exact solutions have
been obtained for processes where the particles can move at most one lattice spacing per
update step. A natural generalization would allow particles to move over larger distances.
This is more realistic with regard to modelling traffic since one usually has a whole spectrum
of allowed car velocities. These generalized exclusion models are thus more appropriate for
comparison with ’experiments’ (i.e. measurements on freeway traffic [14]).
In addition there have been a number of publications dealing with traffic flow in the
framework of statistical mechanics, especially using cellular automata [1,2, 15–17]. In most
of these works two-dimensional traffic has been studied which corresponds to modelling the
complex situation of traffic in the street-network of a city. Again the cars were allowed to
move at most one lattice site per time step. It could be shown numerically [1] that relaxing
this constraint, i.e. the cars can have integer velocities up to an upper speed limit larger than
one, behaves qualitatively as well as quantitatively in good agreement with real traffic (with
an appropriate upper speed limit). It will be shown in this paper that, from an analytical
point of view, the situation changes drastically when turning from maximum velocity one
to higher speed limits. In this case real long-range correlations occur even in the stationary
state which is not true for models with maximum velocity one.
The exclusion models mentioned above may be classified according to the boundary
conditions and dynamics used. There are two relevant types of boundary conditions: pe-
2

riodic and open. Using periodic boundary conditions one considers a ring on which the
cars/particles can move (’Indianapolis situation’). Open boundaries correspond to the so-
called ’bottleneck situation’ where one imposes certain input and output flows at the chain
ends.
Basically one has to distinguish four types of dynamics. The dynamical variables may
be updated one after the other in a certain order (sequential update), one after another in
random order (random-sequential), in parallel for all sites of a given sublattice (sublattice
update) or in parallel for all sites (parallel update). In the case of asymmetric update rules
and sequential dynamics one has to distinguish at least two cases: update in or opposite
to the direction of the traffic motion. Note that for open boundary conditions sequential
and parallel dynamics are identical if the update direction is the same as the direction of
traffic flow. For periodic boundary conditions, however, they are not identical. Here parallel
dynamics correspond to sequential dynamics with a special site which closes the ring. This
site memorizes a car even if it moved away one time-step before. This creates an obstacle
in the ring giving rise to a lower flow through this site compared with sequential dynamics
in the direction of the motion.
In the present paper we use periodic boundary conditions and parallel update. In this
case not much is known exactly, since all of the previous works [5–13] used either random-
sequential or sublattice update. The advantage of parallel update with respect to sublattice
or sequential update is that all sites are equivalent which should be the case in a realistic
model with translational invariance. On the other hand parallel update enhances the possible
system sizes in numerical simulations, especially because parallel or vector computers can
be used easily.
The paper is organized as follows: In Sect. 2 we introduce the model. In Sect. 3 the results
of computer simulations are presented. We describe the different simulation techniques
which have been used and compare their performance. In Sect. 4 we present a mean-field
analysis of the model for arbitrary velocities vmax. In Sect. 5 we introduce the so-called
n−cluster approximation [18]. This improved mean-field method takes into account short-
3

range correlations. In the case vmax = 1 we show that the 2−cluster approximation already
gives the exact result. For vmax = 2 we compare the results of the n−cluster approximation
(n = 1, . . . 6) with the computer simulations of Sect. 3 and find an excellent agreement. In
App. A the exact solution for the case vmax = 1 is derived whereas in App. B the equivalence
of this case to an Ising model with repulsive interactions is shown.
II. THE MODEL
In the following we study single-lane traffic on a ring of length L with periodic boundary
conditions. The model which has been introduced in [1] is defined as follows:
Each of the L sites can either be empty or occupied by one vehicle with velocity v =
0, 1, . . . , vmax. At each discrete time-step t → t + 1 an arbitary arrangement of N cars is
updated according to the following rules:
1) Acceleration: If the velocity v of a vehicle is lower than vmax the speed is advanced
by one [v = v + 1].
2) Slowing down (due to other cars): If the distance d to the next car ahead is not
larger than v (d ≤ v) the speed is reduced to d− 1 [v = d− 1].
3) Randomization: With probability p, the velocity of a vehicle (if greater than zero)
is decreased by one [v = v − 1].
4) Car motion: Each vehicle is advanced v sites.
These rules are applied to all cars in parallel (parallel update). The rules ensure that the
total number N of cars is conserved under the dynamics (which is not true in the bottleneck-
situation). Note that even for parallel update the randomization yields non-deterministic
behaviour. For random-sequential update the probability p > 0 is not essential because it
only rescales the time axis [1].
4

In the simplest case vmax = 1 the cars are allowed to move only one step during an
update. For this situation several results are known [5–13]. Especially it can be shown
that for random-sequential update the mean field Ansatz yields the exact equilibrium state
[1,6] which is equivalent to the fact that for a fixed number of cars every arrangement of
cars occurs with the same probability. Therefore it is quite natural to take the mean-field
approach also as a starting point for the investigation of higher velocities vmax > 1 and
parallel update.
Our main interest will be the calculation of the so-called fundamental diagram (flow
q vs. density ρ = N/L). As described in [1] these results can be compared directly with
measurements of real traffic [14]. One expects a transition from laminar flow to start-stop
waves with increasing car density. For vmax = 1 it is easy to see that the fundamental
diagram is symmetric with respect to ρ = 1/2 due to particle-hole symmetry. This is not
true for realistic traffic where on finds a distinct asymmetry. i.e. the maximum of the flow
is shifted to lower values of ρ (∼ 0.2).
The rules given above take into account the basic features of real traffic. Firstly, they
allow for a spectrum of velocities which seems to be necessary in order to break the particle
hole symmetry of the model with only two velocities vmax = 1. The maximum velocity
still enters the system as a free parameter but it can be argued [1] that most realistic are
values around vmax = 5. Secondly, the acceleration of the cars is very smooth compared to
the deceleration which can occur in only one time step. For simplicity we assume the same
maximum velocity for all cars but this condition is not essential. Finally the randomization
step is necessary to avoid purely deterministic dynamics and to take into account natural
fluctuations of driving.
The order of the update rules given above is crucial. Changing it would also change
the performance of the model. The randomization step could e.g. be placed between accel-
eration and deceleration (1–3–2–4) but then its influence is lowered drastically since every
decelerating car does not feel any randomness. On the other hand one could also change
the starting step. This will not influence the properties of the model but can simplify the
5

calculations. If one begins with step 2 (as we will do for the n-cluster approximation) and
proceeds 3–4–1 then one has the advantage that no cars with velocity zero occur since all
cars were just accelerated by one unit. Therefore the number of possible states of a site is
reduced by one.
III. THE SIMULATIONS
An integral part of our research were computer simulations. Simulations allow to obtain
quantitative insight into a model in relatively short time. In this way, they complement
the analytical work: Simulations are first used to test a variety of models until one has
some overview and the most useful is found, and high quality simulation data are used
to confirm analytical results (see Figs. 6 and 7). Further on, simulations are used to go
beyond the analytically treated cases: either for variations/extensions of the model, or for
more complicated quantities and issues such as the life-time distribution of the simulated
traffic jams [19] or the possibility of self-organized criticality [20]. Last, but definitely no
least, this work are the first steps towards an ultrafast microscopic simulation model for
large traffic networks. We already have results for 2-lane-traffic [21–23] and for real world
network implementations including ramps, intersections, and junctions [21,24,25].
A. Simulation technicalities
Initially, we used a simple code on a workstation (and sometimes its replication on several
nodes of a parallel computer) for getting an overview over the model’s properties and for
estimating its relevance for real world traffic [1]. Later, we implemented vectorizing and/or
parallelized versions in order to obtain faster simulation speeds.
Besides the advantage of getting high quality data in relatively short time, we have
taken advantage of the simplicity of the model to implement it with different algorithms on
many different supercomputers. This gives us intuition on how to implement more complex
models [26], and to make predictions on how fast our simulations will be in comparison to
6

other microscopic traffic simulation models [27]. In the cases of 2-lane-traffic and of the
network implementation, our predictions were quite accurate.
For the practical coding, we considered three different approaches: site oriented, particle
oriented, and an intermediate scheme. Site oriented directly implements the CA: A street is
represented by a chain ~v of integers with values between −1 and vmax, where −1 means that
there is no particle at this site, whereas the other values denote a particle and its velocity. In
contrast, particle oriented means that two lists (xi)i=1,...,N and (vi)i=1,...,N contain position
xi and velocity vi of each particle i (i = 1, . . . , N). This is similar to a molecular dynamics
algorithm, except that particles are constrained to integer positions and velocities.
Obviously, the particle-oriented approach will always be faster than the site-oriented
one for sufficiently low vehicle densities. The particle-oriented approach is more flexible,
and since in single-lane simulations passing of vehicles is not possible, the particle lists are
always ordered, making efficient codes for all kinds of computers easy to write. In addition,
an extension to continuous position and velocity is straightforward [28].
On the other hand, for the site-oriented (cellular automaton) approach single-bit cod-
ing [29] is possible. This means that the model is formulated in logical variables, which
may be stored bitwise into computer words. Logical operations on computer words treat all
bits of the word simultaneously, giving a theoretical speedup of b, where b is the number
of bits per word (usually 32 or 64). However, the practical gain for traffic simulations on a
workstation is much lower because the bit-oriented approach cannot take advantage of the
fact that only a fraction of all sites is occupied by a particle. Nevertheless, we found that,
on a workstation, the single-bit algorithm is faster for densities above 0.05 (for vmax = 5).
In addition, the single-bit code runs very efficiently on a Thinking Machines CM-5 using
dataparallel CM-Fortran and on a NEC-SX/3 traditional vector computer. The simulation
data for the fundamental diagrams have been obtained this way.
Once passing of vehicles is allowed (multi-lane traffic), for the particle oriented approach
efficient memory allocation for parallel and/or vector processors becomes more difficult, and
single-bit coding for the site oriented approach becomes tiresome. These observations led
7

to a third, intermediate approach. As in the site oriented approach, each site is in one of
(vmax +2) states, but for the update only the relevant sites are considered. It turns out (see
below) that on parallel but not vectorizing computers this algorithm is about as fast as the
single-bit version.
In Table I, we give an overview of the computational speeds on selected computers (see
[27] for more details on most of these results). All values are valid for a vehicle density of
ρ = 0.1 and system size of L = 1 333 333 sites, corresponding to 10 000 single lane kilome-
ters. MUPS corresponds to Mega-Updates Per Second, i.e., the number of sites updated
per second divided by 106. These values are useful in order to compare with other imple-
mentations of similar cellular automata or particle hopping systems. The other number is
the extrapolated real time system size, i.e. the extrapolated system size (assuming linear
speed-up) where the computation would be just as fast as reality.
The most noteworthy features of the table are emphasized in boldface:
• On vectorizing computers such as the SX-3 and the CM-5, the single-bit algorithm has
a notable advantage over the other algorithms. On all other machines, the intermediate
algorithm is never more than a factor 2.5 slower than the single-bit algorithm.
• Already on a relatively modest machine such as an Intel Paragon with 64 nodes, our
real time limit is 280 000 single lane kilometers. Since, e.g., the whole freeway net-
work of Germany is about 60 000 single lane kilometers long (12 000 km × 2 directions
× 2.5 lanes), the real use of this computational speed will be (i) real time applica-
tions, where the traffic forecast has to be computed before the situation arrives; and
(ii) Monte Carlo simulations, where many runs are necessary.
• The by far fastest “realistic” traffic microsimulation world-wide to date is the PARAM-
ICS microsimulation project [30]. Their real time limit is ≈ 20 000 km on 16k CM-200,
a machine which is in terms of peak performance slightly faster than a 128-node-
Paragon. In other words, it seems that our not completely realistic car following logic
buys us about an order of magnitude in performance.
8

B. Simulation results
Figs. 1 and 2 show typical time evolutions of different versions of the model. Fig. 1 shows
parallel update with vmax = 1 on the left, and random sequential update with vmax = 5 on
the right, both at density c = 0.5. Random sequential update for vmax = 1 looks the same as
the left picture. Obviously, neither using parallel update for vmax = 1 nor taking a vmax > 1
for random sequential update does change the phenomenlogical behavior from the stochastic
asymmetric exclusion process. For vmax = 1, ρ = 0.5 corresponds to maximum flow, and
in consequence, wave structures do not move in space [31]. For vmax > 1, the point of
maximum throughput is shifted to lower densities, and in consequence, in the left plot the
waves are moving backwards.
The situation is different when one combines parallel update and vmax > 1 (Fig. 2). Here,
in the regime of maximum throughput (left plot) waves are only sparse, and they clearly
move backwards. And even at higher densities (right plot), waves are much more distinct
than in the random sequential case.
In short, one can divide the models between random sequential update with arbitrary
maximum velocity on one hand, and parallel update with maximum velocity vmax ≥ 2 on
the other hand. Parallel update with vmax = 1 phenomenologically is an intermediate case.
Interestingly, it will turn out that this structure is reflected in the analytical calculations
below: For random sequential update, the mean field solution is already exact. For parallel
update and vmax = 1, the situation is only slightly different because already the first step
beyond mean field is exact. The situation is completely different for higher velocities in
connection with parallel update, where the analytical approximations only converge slowly
towards the simulation result.
In addition, Fig. 2 allows an interesting comparison with fluid-dynamical models. Start-
ing from ordered initial conditions, one clearly sees how instabilities develop and produce
the start-stop-waves, very similar to results in [32]. Working out these connections is the
topic of current research [25].
9

Fig. 3 (a) shows current vs. density curves for maximum velocities vmax between 1 and
5, plus for a different fluctuation parameter p = 0.25 at vmax = 5. One clearly sees that
the maximum throughput increases with increasing vmax, whereas the density of maximum
throughput decreases. In reality, the density of maximum throughput lies between ρ =
0.15 and ρ = 0.2; it is given by the maximum speed of trucks which dominate the speed
distribution for traffic near capacity [21]. For that reason, having a higher speed limit for
passenger cars does not help throughput; in many cases, it actually makes things worse [33].
However, reducing the fluctuation parameter p increases throughput enormously. This
is mostly due to the better acceleration behavior in that case [34].
In general, one sees that by varying the parameters vmax and p, the fundamental diagram
can easily be adapted to real traffic situations, although some of the underlying vehicle
dynamics remain somewhat unrealistic: E.g. average acceleration from 0 to 100 km/h takes
place in 10 seconds. That fact that the fundamental diagram is nevertheless quite realistic
is due to the fact that the first time steps fo the acceleration matter most [35]. And here,
4 seconds for an acceleration from 0 to 40 km/h are far more realistic.
Another quantity of interest for traffic engineers are the velocity fluctuations
σ(vloc) :=
√
(vloc − vloc)2 ,
where vloc is the “local” velocity of vehicles crossing a fixed line. (The average of the lo-
cal velocity is different from the usual ensemble average. For example, cars with velocity
zero never enter the local average.) According to measurements and fluiddynamical argu-
ments [33], these fluctuations are a good indicator of traffic near capacity. And indeed do we
find for our model (Fig. 3 (b)) that these fluctuations very abruptly increase near capacity.
IV. MEAN-FIELD THEORY
The complete analytic solution of the traffic model is not possible except for the case of
maximum velocity vmax = 1 (see Sect. 5.1 and Appendix A). Therefore some approximation
10

is necessary when one tries to study this model analytically for vmax > 1. The following
mean-field type theory will be the first step of such an analytical approach.
The calculation starts from the stochastic description of the traffic model. Instead of
specifing the car position and their velocities, we analyze the probability distribution of each
site at each time step. We denote the probability that there is no car at site i (i = 1, 2, 3,
· · ·, L) at timestep t by d(i, t) and the probability that there is a car with velocity α (α = 0,
1, 2, · · ·, vmax) at site i and timestep t by cα(i, t). The c and d distributions together take
into account all possible states of the different sites. Therefore one has the normalization
condition for all sites and all timesteps
d(i, t) + c0(i, t) + c1(i, t) + c2(i, t) + c3(i, t) + · · · cvmax(i, t) = 1 . (4.1)
Denoting with c(i, t) the total probability for site i to be occupied at timestep t, i.e.
∑vmax
j=0 cj(i, t), one simply has d(i, t) + c(i, t) = 1.
According to the update rules in four stages (see Sect. 2) the time evolution of these
probability distributions can be described by the following four sets of equations:
• Acceleration Stage
c0(i, t1) = 0
cα(i, t1) = cα−1(i, t) , 0 < α < vmax (4.2)
cvmax(i, t1) = cvmax(i, t) + cvmax−1(i, t)
• Deceleration Stage
c0(i, t2) = c0(i, t1) + c(i + 1, t1)
vmax
∑
β=1
cβ(i, t1)
cα(i, t2) = c(i + α + 1, t1)
α
∏
j=1
d(i + j, t1)
vmax
∑
β=α+1
cβ(i, t1) (4.3)
+ cα(i, t1)
α
∏
j=1
d(i + j, t1), 0 < α < vmax
cvmax(i, t2) =
vmax
∏
j=1
d(i + j, t1)cvmax(i, t1)
11

• Randomization Stage
c0(i, t3) = c0(i, t2) + pc1(i, t2)
cα(i, t3) = qcα(i, t2) + pcα+1(i, t2), 0 < α < vmax (4.4)
cvmax(i, t3) = qcvmax(i, t2)
• Motion Stage
cα(i, t + 1) = cα(i− α, t3), 0 ≤ α ≤ vmax (4.5)
The time t is assumed to take on only integer values. t1, t2 and t3 denote the intermediate
time steps between t and t + 1 (sometimes defined as t + 1/4, t + 1/2 and t + 3/4, respec-
tively). These time-evolution equations conserve for periodic boundary conditions the total
number of cars at each stage . This formulation of the dynamics neglects spatial correlations
completely since one assumes that all expectation values factorize into local terms.
The variables c and d are real valued between 0 and 1. The stochastic description
originates from the stochastic nature of the third randomization step. The other three steps
are purely deterministic.
These time evolution equations are non-linear and further analysis of the full system has
not been successful up to now. However, in the stationary state, i.e. in the limit t → ∞, the
c and d distributions become homogeneous in space (for periodic boundary conditions) and
therefore, apart from the time dependence, also the site dependence can be omitted. Using
this and combining the four update steps one gets the following set of vmax + 1 equations:
c0 = (c + pd)c0 + (1 + pd)c
vmax
∑
β=1
cβ
cα = dα[qcα−1 + (qc + pd)cα + (q + pd)c
vmax
∑
β=α+1
cβ], (4.6)
0 < α < vmax − 1
cvmax−1 = dvmax−1[qcvmax−2 + (qc + pd)(cvmax−1 + cvmax)]
cvmax = qdvmax [cvmax−1 + cvmax ]
12

These equations essentially describe the motion of a single car with statistical (densitity
dependent) ’obstacles’. A remarkable feature of the equations is that they are linear when
one specifies the density c = 1 − d of cars. Therefore (4.6) can be recast in matrix form
as A~c = ~c. The matrix A can be read off from (4.6), ~c is the vector with elements cα,
α = 0, ..., vmax. For small vmax we can invert A to find the densities cα explicitly. For large
values of vmax it is more convenient to write down a recursion relation in order to obtain
the steady state solution.
From the first equation in (4.6) one can determine c0 directly without knowledge of the
other c′αs to give
c0 = c2
1 + pd
1− pd2 . (4.7)
Using this result and the second equation of (4.6) one can also write down the expression
for c1
c1 = qc2d
1 + d + pd2
(1− pd3)(1− pd2) . (4.8)
For α larger than one a recursion relation can be derived calculating cα − dcα−1 using again
the second equations of (4.6)
cα =
1 + (q − p)dα
1− pdα+2 dcα−1 −
qdα
1− pdα+2 cα−2 (4.9)
for α = 2, 3, · · ·, vmax−2. Therefore, starting with the expressions (4.7), (4.8) for c0 and c1,
one can estimate c2, c3, · · ·, cvmax−2 recursively. These results do not depend on the actual
value of vmax and thus are valid generally (provided vmax ≥ 2).
Finally, from the last two equations of (4.6) one gets
cvmax−1 =
1− qdvmax
1− dvmax−1(q + pd)qd
vmax−1cvmax−2 (4.10)
cvmax =
qdvmax
1− qdvmax cvmax−1. (4.11)
The vmax dependence only occurs in these two quantities. In Fig. 4 some of the densities are
shown for large vmax. The densities of the fast cars go to zero rapidly, since one expects an
13

exponentially fast decay from the recursion relations, (4.9)–(4.11). The contributions from
cars with high velocities therefore are neglegible. We will mainly be interested in the flow
f(c, p), defined by
f(c, p) =
vmax
∑
α=1
αcα. (4.12)
In the limit of vmax →∞ it is possible to carry on the analysis using the iteration equations
(4.7)–(4.9) and the generating function
g(z) =
∞
∑
α=0
zαcα. (4.13)
As usual one simply has g′(1) = f(c, p). The equations (4.7)–(4.9) now can be combined to
give one single equation for g(z)
g(z)[1− dz]− g(dz)d2(1− z)[p + qz] = c2 + pc2d(1− z). (4.14)
Since vmax is infinite one does not have to worry about the upper boundary equations (4.10)
and (4.11) for cvmax−1 and cvmax whose contribution is neglegible. One should notice that
the generating function g occurs with two different arguments (z and dz) which makes it
impossible to solve this (linear) equation for g(z) explicitly.
After differentiation an expression for g′(1) is obtained
g′(1) = d(1− pc)− g(d)d
2
c (4.15)
The problem therefore is reduced to find an expression for g(d). This can be found by
successive application of equation (4.14) with z = dn, n = 1, 2, 3, . . .. The final result is an
asymptotic expression for g′(1) = f(c, p)
f(c, p) = qcd
[
1 +
∞
∑
n=1
d2n
n−1
∏
l=0
(p + qdl)
]
(4.16)
In Fig. 5 the flow f is shown as a function of the concentration c of cars for various deceler-
ation probabilities p. The slope at the origin (c ∼ 0) of the fundamental diagram is infinite
whereas the slope for c ∼ 1 is −q. This mean-field result yields, compared with the simu-
lation data shown above, much too small values of the flow. This can easily be understood
14

since the reduction to a single car problem ignores all spatial correlations of the cars. Cars
with high velocities tend to be equidistant and can therefore maintain a high velocity with a
larger probability than in the mean-field system where it is much more difficult to accelerate
and stay at high velocities over a certain period.
Furthermore, even for vmax = 1, the mean-field solution does not coincide with the exact
result shown below. Note that for random sequential update the mean-field solution is exact
for vmax = 1.
V. IMPROVED MEAN-FIELD THEORIES
In order to improve the simple mean-field theory of the preceeding section we have to
take into account correlations between neighbouring sites [18].
We divide the lattice into segments or clusters of length n (n = 1, 2, . . .) such that two
neighbouring clusters have n − 1 sites in common. The probability to find in equilibrium
a cluster in state (σ1, . . . , σn) will be denoted by Pn(σ1, . . . , σn). Due to the translational
invariance of equilibrium state of the system with periodic boundary conditions one has not
to specify the actual location of the n-spin cluster. In order to simplify the calculations
we apply, as mentioned above, the four update-rules in the order 2–3–4–1 instead of 1–2–
3–4 [18]. This has the advantage that after one update cycle one ends up with step 1 and
therefore no car has velocity 0. It follows that every site j is in one of the vmax states
σj = 0, 1, . . . , vmax where now 0 denotes an empty site. This change in the ordering finally
has to be taken into account in the calculation of the flow.
The equilibrium probabilities for an n-site cluster are then determined by
Pn(σ1, . . . , σn) =
∑
{τj}
W (σ1, . . . , σn|τ−vmax+1, . . . , τn+vmax)×
× P2vmax+n(τ−vmax+1, . . . , τn+vmax) , (5.1)
where the probability P2vmax+n(τ−vmax+1, . . . , τn+vmax) for clusters of length 2vmax + n has to
be expressed through the n-cluster probabilities Pn(τ1, ..., τn). This enlargement of the clus-
15

ter length occurs since all cars which can drive into or out of the cluster (σ1, . . . , σn) within
the next timestep contribute to the transition rates W ({σj}|{τj}). Thus we have to take
into account not only the given cluster but also the vmax sites to its left (with state variables
τ−vmax+1, . . . , τ0) and the vmax sites to its right (with state variables τn+1, . . . , τn+vmax). The
decomposition of the (2vmax + n)-cluster can be carried out by introducing the conditional
probabilities
P (τi|τi+1, . . . , τi+n−1) =
Pn(τi, τi+1, . . . , τi+n−1)
∑
τ Pn(τ, τi+1, . . . , τi+n−1)
(5.2)
at the left side and
P (τi, . . . , τi+n−2|τi+n−1) =
Pn(τi, τi+1, . . . , τi+n−1)
∑
τ Pn(τi, . . . , τi+n−2, τ)
(5.3)
at the right side. With this definition we rewrite P2vmax+n (in the n-cluster approximation)
in the following form:
P2vmax+n(τ−vmax+1, . . . , τn+vmax) =
0
∏
i=−vmax+1
P (τi|τi+1, ..., τi+n−1)× (5.4)
× Pn(τ1, . . . , τn)
vmax
∏
i=1
P (τi+1, ..., τi+n−1|τi+n) .
The transition probability W (σ1, . . . , σn|τ−vmax+1, . . . , τn+vmax) depends on the probability p
and vanishes if the configuration (τ−vmax+1, . . . , τn+vmax) cannot evolve in one timestep into
(σ1, . . . , σn) according to the rules 1) − 4) of Sect. 2. If W is non-zero it is of the form
pn1qn2 with integers n1, n2 describing how many cars have to be decelerated (n1) through
the randomization step when total number of cars which can drive is n1 + n2. With the
approximation (5.4) it is possible to write down a closed system of equations for the n-
cluster probabilities Pn(σ1, . . . , σn). The number of the equations is given by (vmax + 1)n,
the total number of possible configurations of n site variable with vmax + 1 possible states1
(without change of the order of the update steps one would have (vmax + 2)n equations).
1In practice some of these equations turn out to be trivial so that the relevant number is less than
(vmax + 1)n.
16

Due to the exponential growth with respect to n one is, especially for larger vmax, restricted
to only small cluster lengths n (For the realistic value of vmax = 5 one has for the 2-cluster
approximation already 36 equations!)
The above approximation converges for n → ∞ to the exact solution for an infinite
system (i.e. in the thermodynamic limit L → ∞. This approximation scheme is well known
in the literature. It is analogous to the (n, n− 1)−cluster approximation of [36], the n− 1
step Markovian approximation of [37] or the local structure theory of [38]. It goes back to
the probability path method introduced by Kikuchi [39].
With the knowledge of the n-cluster probabilities Pn(σ1, . . . , σn) it is then easy to calcu-
late the fundamental diagram, i.e. the flow f as a function of the density c of cars. Since we
have changed the order of the update steps one has to take this into account by performing
the steps 2–3–4 at the end since the last step must be 4 (≡ car motion). After that one
simply can calculate the density cα of cars which will drive α sites in the next timestep
cα =
∑
σ2,...,σn
Pn(α, σ2, . . . , σn) (5.5)
and then one proceeds as in the mean-field approximation eqn. (4.12) and calculates f =
∑αcα.
A. vmax = 1
In the case vmax = 1 the site variables take on only the values σ = 0, 1 where σ = 0
means no car and σ = 1 a car with velocity one. In the 2-cluster approximation ((5.1) with
n = 2) one has, according to the above arguments, a system of 4 equations. This can be be
reduced to only one equation for P (1, 0) very easily with the help of the relations
P (1, 0) = P (0, 1)
P (0, 0) = 1− c− P (1, 0), (5.6)
P (1, 1) = c− P (1, 0).
17

The first equation is due to the ’particle-hole’ symmetry P (1, 0) = P (0, 1) (in a closed ring
one must have the same number of (0,1) and (1,0) pairs, therefore occuring with the same
probability). The other two equations describe the conservation of cars in the system. The
remaining equation for P (1, 0) reads
qP 2(1, 0)− P (1, 0) + c(1− c) = 0. (5.7)
In the thermodynamic limit we therefore find [18]
P (0, 1) = P (1, 0) =
1−
√
1− 4qc(1− c)
2q . (5.8)
Going to the 3- and higher-cluster approximations one finds that the solution remains the
same indicating that this is the exact result. In App. A we indeed prove that the solution
(5.8) is exact in the thermodynamic limit. For finite systems the prove is also valid but one
has to take into account an additional correction (normalization) due to the constraint of a
fixed number N of cars.
The correct result, valid for any system size, is is given by
P(N,L) = 1N
∑
{σ}
′
L
∏
j=1
P (σj , σj+1) (5.9)
where N denotes the normalization constant and the sum ∑′ runs over all configurations
with a fixed number N of cars (i.e. ∑Lj=1 σj = N).
In contrast to random-sequential dynamics parallel dynamics leads to an effective at-
traction between ’particles’ and ’holes’ (i.e. P (0)P (1) = c(1 − c) ≤ P (0, 1)) and thus to
a higher flow. In App. B the mapping of this model to an equivalent Ising-model with
antiferromagnetic next-nearest-neighbour interactions and an nonvanishing external field is
shown. Due to the antiferromagnetic interactions the system shows an effective attraction
of cars over two lattice sites, thus taking into account the well known effect of car bunching
or platooning [40–42].
In order to calculate the flow one first has to perform the steps 2–3–4 yielding the
probabilities P˜ (τ) to find a site in the state τ = −1, 0, 1 (where now an empty site is
denoted by τ = −1 ) after the fourth step of the updating procedure:
18

P˜ (−1) = 1− c, P˜ (0) = qP (1, 1), P˜ (1) = qP (0, 1) (5.10)
Therefore the fow is finally given by
f(p, c) = qP (1, 0) =
1−
√
1− 4qc(1− c)
2
. (5.11)
In Fig. 6 the exact result for the flow is shown (from the 2-cluster approximation) in com-
parison to numerical simulations and the mean-field approximation (p = 1/2). One can see
that the numerical and analytical data are in excellent agreement. In the deterministic case
p = 0 the flow is a linear function of c: f = (1− |1− 2c|)/2 [18].
The mean velocity per car v¯ is then
v¯ = 1ρ
vmax
∑
τ=0
τP˜ (τ) = P˜ (1)/ρ . (5.12)
In the deterministic case p = 1 this simplifies to
v¯ =



1 for 0 ≤ c ≤ 1/2
(1− c)/c for 1/2 ≤ c ≤ 1
(5.13)
which is the result of [2].
B. vmax = 2
The case vmax = 2 is qualitatively very different from the case vmax = 1. The flow
diagram is no longer symmetric around c = 1/2. The n−cluster approximation seems not
to become exact for any finite n, i.e. in this case long-ranged correlations are important.
We have calculated the fundamental diagram for n = 1, . . . , 5. As shown in Fig. 7 the
approximation converges fast to the results obtained by simulations. The difference between
n = 4 and n = 5 is less then 1%.
The observation that the approximation does not become exact for small n reflects the
fact that the physics for vmax ≥ 2 is distinctly different from the case vmax = 1. As explained
in Sec. 2 the regime looks qualitatively different both from vmax = 1 (arbitrary update) as
well as from random sequential update (arbitrary vmax). Moreover, in the case vmax ≥ 2
19

(parallel update), jams show characteristic branching behavior which is not observed for
vmax = 1 [20].
VI. CONCLUSIONS
We have introduced and investigated a statistical model capable to describe accurately
real traffic. Through the introduction of higher velocities it was possible to produce the
asymmetric fundamental diagrams and the characteristic start-stop waves observed in real
traffic. Through simulation, different regimes depending on the type of update and the max-
imum velocity have been identified. Realistic for traffic is a model with parallel stochastic
update and a suitable choice of the maximum velocity vmax = 5. Already for vmax = 1
the stationary state is more complicated than for random sequential update. An effective
’antiferromagnetic’ interaction between cars favors equal spacing and in consequence higher
throughput.
Taking into account two-point correlations the case vmax = 1 can be solved exactly.
This is no longer true for higher maximum velocties where correlations become long ranged.
Nevertheless, it could be shown that the n-cluster approximation converges fast to the
simulation data.
ACKNOWLEDGMENTS
This work was performed within the SFB 341 Ko¨ln-Aachen-Ju¨lich supported by the DFG.
KN is supported by the “Graduiertenkolleg Ko¨ln/St. Augustin”. The Center for Parallel
Computing ZPR Ko¨ln, the German National Research Center KFA Ju¨lich, the “Rechenzen-
trum” of the University of Cologne, and the “Institut fu¨r Angewandte Informatik” of the
University Wuppertal provided computer time.
20

APPENDIX A: EXACT SOLUTION FOR VMAX = 1
In this Appendix we show that the stationary state of the model with vmax = 1 is in fact
given by equation (5.6) with an appropriate normalization constant Z. The complete set of
evolution equations for parallel update reads:
Pt+1({σ}) =
∑
{τ }
W ({σ} | {τ}) · Pt({τ}), (A.1)
where Pt({σ}) denotes the probability for state {σ} = {σ1, . . . , σL} at time t. The transition
probability W ({σ} | {τ}) to move in one timestep from state {τ} to state {σ} factorizes
into local terms:
W ({σ} | {τ}) =
L
∏
i=1
W (σi, σi+1 | τi, τi+1). (A.2)
The only non-vanishing transition probabilities are given by:
W (0, 1 | 1, 0) = 1− p
W (1, 0 | 1, 0) = p
W (1, σ′|σ, 1) = 1 (A.3)
W (σ, 0|0, 0) = 1
W (0, σ|0, 1) = 1, σ, σ′ = 0, 1
whereas W (τ, τ ′ | σ, σ′) = 0 in all other cases. (Note that in eqn. (5) of ref. [1] due to a
misprint a factor (1− σi+1)/2 is missing on the right-hand side.) We now make the Ansatz
that the probability in the stationary state P ({σ}) factorizes into local two-site terms Pσiσi+1 :
P ({σ}) =
L
∏
i=1
Pσiσi+1 (A.4)
with periodic boundary condition σL+1 = σ1. We will define nσ,σ′(σ) as the number of pairs
of next neighbours (σ, σ′) in a particular state σ. Due to the particle-hole symmetry of the
system one has n01 = n10 and the following simple relations for the system size L and the
(conserved) number of particles N hold:
21

L = 2 · n01 + n11 + n00 (A.5)
N = n01 + n11. (A.6)
In the stationary state equation (A.1) becomes time-independent and the states {τ} in the
summation on the right-hand side can be classified according to the number of particles l
which have to be moved in order to obtain the new state σ:
(P01P10)n01P n1111 PL−2n01−n1100 = (A.7)
n01
∑
l=0
∑
∆
g(n01, l,∆)(1− p)lpn01−l−∆(P01P10)n01−∆P n11+∆11 PL−2n01−n11+∆00 .
The summation index ∆ is defined as ∆ = n11(τ )− n11(σ). Therefore the range of the ∆-
summation depends on the particular state σ. The function g(n01, l,∆) counts the number
of possible states τ leading to state σ given fixed values of n01, l and ∆ and reflects a kind
of degeneracy. Summing g over ∆ yields simply
∑
∆
g(n01, l,∆) =
(
n01
l
)
(A.8)
Dividing eqn. (A.7) by the left-hand side one gets
1 =
n01
∑
l=0
∑
∆
g(n01, l,∆)(1− p)lpn01−l
(
P00P11
pP01P10
)∆
. (A.9)
Setting the expression in the parenthesis on the right-hand side to 1, (A.9) reduces with the
help of (A.8) to an identity. Therefore the condition for P (σ, σ′) reads
P00P11 = pP01P10 (A.10)
Note that P (σ, σ′) is not normalized for a finite system and therefore a normalization con-
stant N has to be taken into account.
APPENDIX B: MAPPING TO AN EQUIVALENT ISING-MODEL
Introducing Ising-variables σi = ±1 instead of the lattice-gas variables τi = 0, 1 (σi =
2τi − 1) one can look at the steady state as the equilibrium distribution P (σ) ∼ e−βH(σ)
22

of an Ising model with Hamiltonian H = −J ∑i σiσi+1 + h
∑
i σi In order to determine the
coupling constant J and the external field h one interpretes the local probabilities Pσ,σ′ as
transfer matrices Pσ,σ′ = αe−Jσσ′−h(σ+σ′)/2, (β ≡ 1). From eqn. (A.10) it follows directly that
e4J = p or J = log(p)/4 < 0. Therefore the corresponding Ising model has antiferromagnetic
interactions. According to eqns. (A.5) and (A.6) one has in addition
1 = 2 · P01 + P11 + P00 (B.1)
ρ = P01 + P11. (B.2)
Dividing the two expressions on both sides of the two equations gives one eqn. without the
constant α to determine the external field as
eh =
√
1− 4(1− p)ρ(1− ρ)− (1− 2ρ)
2
√p(1− ρ) . (B.3)
Therefore the steady state corresponds to an Ising-model with antiferromagnetic (repulsive)
interactions. Due to the conservation of the total number of particles (cars) one has to
impose the constraint of a fixed magnetization to the Hamiltonian. The normalization is
then simply the partition function calculated under this constraint.
23

REFERENCES
[1] Nagel K and Schreckenberg M 1992 J. Phys. I (France) 2 2221
[2] Biham O, Middleton A A and Levine D 1992 Phys. Rev. A 46 R6124
[3] Wolfram S 1986 Theory and Applications of Cellular Automata (Singapore: World Sci-
entific)
[4] Frisch U, Hasslacher B and Pomeau Y 1986 Phys. Rev. Lett. 56 1505
[5] Janowsky S A and Lebowitz J L 1992 Phys. Rev. A 45 618
[6] Gwa L-H and Spohn H 1992 Phys. Rev. A 46 844
[7] Derrida B, Domany E and Mukamel D 1992 J. Stat. Phys. 69 667
[8] Derrida B and Evans M R 1993 J. Phys. I (France) 3 311
[9] Derrida B, Evans M R, Hakim V and Pasquier V 1993 J. Phys. A: Math Gen. 26, 1493
[10] Kandel D, Domany E and Nienhuis B 1990 J. Phys. A: Math Gen. 23 L755
[11] Schu¨tz G and Domany E 1993 J. Stat. Phys. 72 277
[12] Schu¨tz G 1993 Phys. Rev. E 47 4265
[13] Schu¨tz G 1993 J. Stat. Phys. 71 471
[14] Hall F L, Brian L A and Gunter M A 1986 Transp. Res. A 20 197
[15] Cuesta J A, Mart´inez F C, Molera J M and Sa´nchez A 1993 Phys. Rev. E 48 4175
[16] Nagatani T 1993 J. Phys. Soc. Jpn. 62 1085
[17] Nagatani T 1993 Physica A 198 108
[18] Schadschneider A and Schreckenberg M 1993 J. Phys. A: Math Gen. 26 L679
[19] Nagel K 1994 Int. Mod. Phys. C 5 567
24

[20] Nagel K and Paczuski M, Emergent Traffic Jams (in preparation)
[21] Rickert M 1994 Simulationen zweispurigen Autobahnverkehrs mit Zellularautomaten
Master Thesis University of Cologne
[22] Latour A 1993 Simulation von Zellularautomaten-Modellen fu¨r Mehrspurverkehr Master
Thesis University of Cologne
[23] Rickert M, Latour A , Nagel K, and Schreckenberg M (in preparation)
[24] Nagel K and Rasmussen S 1994, in: Brooks R, Maes P (eds.) Proceedings of Alife 4
MIT press Cambridge MA USA
[25] Nagel K 1994 Ph.D. thesis (in preparation)
[26] TRANSIMS—The TRansportation ANalysis and SIMulation System project, Los
Alamos National Laboratory
[27] Nagel K and Schleicher A 1994 Parallel Comput. 20 125
[28] Nagel K and Herrmann H J 1993 Physica A 199 254
[29] Stauffer D 1991 J. Phys. A: Math Gen. 24 91 909
[30] Wylie B J N, McArthur D and Brown M D PARAMICS parallelisation schemes EPPC-
PARAMICS-CT.10
[31] Krug J 1991 in: Riste T and Sherrington D (eds.) Spontaneous formation of space-time
structures and criticality Kluwer Academic Publishers Netherlands 37
[32] Kerner B S and Konha¨user P 1994 Phys. Rev. E 50 54
[33] Zackor H, Ku¨hne R and Balz W 1988 Untersuchungen des Verkehrsablaufs im Bere-
ich der Leistungsfa¨higkeit und bei instabilem Fluß Forschung Straßenbau und Straßen-
verkehrstechnik 524, Bundesminister fu¨r Verkehr Bonn–Bad Godesberg
[34] Nagel K 1994 Fast low fidelity microsimulation of vehicle traffic on supercomputers
25

Paper No. 94 09 01 Transportation Research Board Meeting Jan. 1994 Washington
D.C. USA
[35] Piper H P 1991 Internationales Verkehrswesen 43 489
[36] ben-Avraham D and Ko¨hler J 1992 Phys. Rev. B 45 8358
[37] Crisanti A, Paladin G and Vulpiani A 1993 Products of Random Matrices in Statistical
Physics (Berlin: Springer)
[38] Gutowitz H A, Victor D and Knight B W 1987 Physica D 28 18
[39] Kikuchi R 1966 Prog. theor. Phys., Suppl. 35 1
[40] May A D Traffic flow fundamentals Prentice Hall, Englewood Cliffs, NJ, 1990
[41] Ben-Naim E, Krapivsky P L, and Redner S Kinetics of clustering in traffic flows cond-
mat preprint 1994
[42] Nagatani T Bunching of cars in asymmetric exclusion models for freeway traffic (sub-
mitted PRE)
26

TABLES
s.bit (F77) particle (F77) intermed. (C) netw. (C)
Sparc10 4.0 MUPS 3.4 MUPS 1.9 MUPS 1.2 MUPS
30 000 km 14 000 km 8800 km
PVM 19.0 MUPS 8.9 MUPS
(5× Sp10) 140 000 km 65 000 km
SX-3/11(1) 533 MUPS 2.8 MUPS
1 CPN 4000 000 km 21 000 km
GCel-3 102 MUPS 211 MUPS 121 MUPS
1024 CPNs 750000 km 1 550000 km 900000 km
iPSC 83 MUPS 35 MUPS
32 CPNs 630000 km 260000 km
Paragon
64 CPNs 290 000 km
CM-5(1) 173 MUPS 30 MUPS
32 CPNs 1 300 300 km 220000 km
CM-5(1)
1024 CPNs [> 1.7 e 6 km]
(1) CPN(s) has/have vector units (SIMD instruction set)
(2) using data parallel Fortran (CMF)
(3) using message passing (CMMD)
27

TABLE I. Computational speed of different algorithms and different machines. The machines
are a SUN Sparc10 Workstation, 5 coupled Workstations Sparc10 under PVM, a NEC SX-3/11
traditional vectorcomputer, a Parsytec GCel-3 massively parallel computer with 1024 relatively
slow T805-CPU’s, an Intel iPSC/860 Hypercube with 32 nodes, and Intel Paragon XP-10/S with
64 nodes, and a Thinking Machines CM-5 with 32 nodes, each node containing one Sparc processor
and 4 vector units. — “s.-bit” refers to the single-bit coded site-oriented algorithm, “particle” and
“intermed.” to the particle-oriented and the intermediate one, and “netw.” refers to a network
implementation of the freeway network of Germany. — All values are valid for a vehicle density
of c = 0.1 and system size of L = 1333 333 sites, corresponding to 10 000 single lane kilometers.
MUPS corresponds to Mega-Updates Per Second, i.e., the number of sites updated per second
divided by 106. The other number is the extrapolated real time system size, i.e. the extrapolated
system size (assuming linear speed-up) where the computation would be just as fast as reality.
Values in brackets [] are estimates.
28

FIGURES
FIG. 1. Evolution of different automata rules from an ordered initial state of density c = 0.5.
Particles are moving to the right. Left: Same as (stochastic) asymmetric exclusion, except that the
update is parallel. The same plot for exact asymmetric exclusion looks similar. Note that waves
do not move in space. Right: Same as asymmetric exclusion, except that the maximum velocity
is 5. Waves are now moving backwards, indicating a density above maximum flow. — Neither
the parallel update nor the higher maximum velocity alone are sufficient to change the qualitative
dynamics of the asymmetric exclusion model.
FIG. 2. Evolution of our traffic model (maximum velocity vmax = 5, parallel update) for
density c = 0.3 (left) and c = 0.1 (right) from ordered initial conditions. Random sequential
update with vmax = 5 at density c = 0.3 (not shown) looks qualitatively similar to Fig. 1 left,
whereas at density c = 0.1 and random sequential update the waves are moving to the right (not
shown). The density of maximum throughput is much lower for the parallel update, and instead of
the waves switching from backward to forward motion at this point, they tend to vanish completely
(cf. [20]). Moreover, the right picture illustrates the instability mechanism similar to the continuous
model of [32].
FIG. 3. (a) Fundamental diagrams flow j vs. density c for maximum velocities vmax = 1, 2, . . .,
5, and for a different fluctuation parameter p = 0.25 instead of 0.5 at vmax = 5. (b) Fluctuations
of local speed as a function of density (vmax = 5 and p = 0.5).
FIG. 4. Partial densities c0(c), c1(c), . . ., c5(c) (i.e. for velocities 0 to 5) for mean field approx-
imation. p = 0.5
FIG. 5. Flow for vmax = ∞ in mean field approximation. From bottom to top, the randomiza-
tion parameter p is 0.1, 0.3, 0.5, 0.7, and 0.9.
FIG. 6. Convergence of the n-cluster approximations to the simulation result for the case
vmax = 1 and p = 0.5. Already the 2-cluster approximation is exact.
29

FIG. 7. Convergence of the n-cluster approximations to the simulation result for vmax = 2 and
p = 0.5. Already the 5-cluster approximation gives a good fit of the simulation data.
30