# Critical behavior of a cellular automaton highway traffic model

## 1 Introduction

Cellular automaton (CA) models of traffic flow have attracted much interest since the publication of the Nagel-Schreckenberg (NS) model . The NS model is a probabilistic CA model of traffic flow on a one-lane highway. The road is represented by a lattice of L sites with periodic boundary conditions. Each site is either empty (in state $`e`$) or occupied by a car with velocity $`v=0,1,2,\mathrm{},v_{\mathrm{max}}`$ (in state $`v`$). If $`d_i`$ is the distance between car $`i`$ and car $`i+1`$ (cars are moving to the right), velocities are updated in parallel according to the following two subrules:

$`v_i(t+1/2)`$ $`=\mathrm{min}(v_i(t)+1,d_i(t)1,v_{\mathrm{max}})`$ (1)   
$`v_i(t+1)`$ $`=\mathrm{max}(v_i(t+1/2)1,0)\text{with probability}p,`$ (2)   
where $`v_i(t)`$ is the velocity of car $`i`$ at time $`t`$; then, cars are moving according to the subrule:

$$x_i(t+1)=x_i(t)+v_i(t+1),$$
(3)   
where $`x_i(t)`$ is the position of car $`i`$ at time $`t`$. The model contains three parameters: the maximum speed $`v_{\mathrm{max}}`$, which is the same for all cars, the braking probability $`p`$, and the car density $`\rho `$. For a clear presentation of the various approximate techniques to calculate relevant physical quantities of the NS model compared with results obtained from computer simulations see .

If $`p=0`$, the NS model is deterministic and the average velocity over the whole lattice is exactly given by

$$v=\mathrm{min}(v_{\mathrm{max}},\frac{1}{\rho }1).$$
(4)   
This expression shows that, below a critical car density $`\rho _c=1/(v_{\mathrm{max}}+1)`$, all cars move with a velocity equal to $`v_{\mathrm{max}}`$, while above $`\rho _c`$, the average velocity is less than $`v_{\mathrm{max}}`$. This transition from a free-moving regime to a congested regime is usually viewed as a second-order phase transition.

If $`p0`$, except in the case $`v_{\mathrm{max}}=1`$, no exact expression for the average velocity has been obtained, and there is no consensus concerning the existence of the phase transition (see, for instance  and  and references therein). One point that is particularly unclear is the definition of the order parameter.

If CA models of traffic flow exhibit second-order phase transitions, we have to understand the nature of the order parameter, show how it is related to symmetry-breaking, determine the symmetry-breaking field conjugate to the order parameter, define the analogue of the susceptibility, study the critical behavior, and, possibly, find scaling laws.

The purpose of this short note is to fulfill this program. In our discussion we shall consider a slightly simpler model of traffic flow. We shall assume that the acceleration, which is equal to 1 in the NS model, has the largest possible value (less or equal to $`v_{\mathrm{max}}`$) as in the Fukui-Ishibashi (FI) model , except that these authors apply random delays only to cars whose velocity is equal to $`v_{\mathrm{max}}`$. That is, in our model, we just replace (1) by

$$v_i(t+1/2)=\mathrm{min}(d_i(t)1,v_{\mathrm{max}}).$$
(5)   
## 2 Symmetry considerations

Deterministic CA rules modeling traffic flow on one-lane highways are number-conserving ($`\rho =\text{constant}`$). Limit sets of number-conserving CA have, in most cases, a very simple structure and, these limit sets are reached after a number of time steps proportional to the lattice size . The limit set of our model is identical to the limit set of the NS or FI models.

If $`\rho <\rho _c`$, any configuration in the limit set consists of “perfect tiles” of $`v_{\mathrm{max}}+1`$ cells as shown below

| $`e`$ | $`e`$ | $`\mathrm{}`$ | $`e`$ | $`v_{\mathrm{max}}`$ |
| --- | --- | --- | --- | --- |

in a sea of cells in state $`e`$. In particular, it can be shown  that the probability to find in the limit set a sequence of $`v_{\mathrm{max}}+1`$ empty sites is exactly given by

$$P(\underset{v_{\mathrm{max}}}{\underset{}{ee\mathrm{}e}})=1(v_{\mathrm{max}}+1)\rho .$$

This relation is a simple consequence of the following argument: If, in an empty road, you add $`L\rho `$ cars, in the limit set, for $`\rho <\rho _c`$, we have

$$LP(\underset{v_{\mathrm{max}}}{\underset{}{ee\mathrm{}e}})=L(v_{\mathrm{max}}+1)L\rho ,$$

where $`L`$ is the lattice size.

If $`\rho >\rho _c`$, a configuration belonging to the limit set only consists of a mixture of tiles containing $`v+1`$ cells of the type

| $`e`$ | $`e`$ | $`\mathrm{}`$ | $`e`$ | $`v`$ |
| --- | --- | --- | --- | --- |

where $`v=0,1,\mathrm{},v_{\mathrm{max}}`$. If $`\{\rho _vv=0,1,2,\mathrm{},v_{\mathrm{max}}\}`$ is the velocities distribution, we have

$`\rho `$ $`={\displaystyle \underset{v=0}{\overset{v_{\mathrm{max}}}{}}}\rho _v`$   
$`1`$ $`={\displaystyle \underset{v=0}{\overset{v_{\mathrm{max}}}{}}}(v+1)\rho _v`$   
$`v`$ $`={\displaystyle \frac{1}{\rho }}\left({\displaystyle \underset{v=0}{\overset{v_{\mathrm{max}}}{}}}v\rho _v\right).`$   
Note that (4) is a simple consequence of these relations.

If we introduce random braking, then, even at low density, some tiles become defective, which causes the average velocity to be less than $`v_{\mathrm{max}}`$. The random-braking parameter $`p`$ can, therefore, be viewed as a symmetry-breaking field, and the order parameter, conjugate to that field is

$$m=v_{\mathrm{max}}v.$$
(6)   
This point of view implies that the phase transition characterized by $`m`$ will be smeared out in the presence of random braking as in ferromagnetic systems placed in a magnetic field.

From (6) and (4), it follows that, for $`p=0`$,

$$m=\{\begin{array}{cc}0\hfill & \text{if }\rho \rho _c,\hfill \\ [0.5em]\frac{\rho \rho _c}{\rho \rho _c}\hfill & \text{otherwise}.\hfill \end{array}$$
(7)   
Hence, the critical exponent $`\beta `$ is equal to 1.

## 3 Approximate techniques

To determine the other critical exponents, we have used local structure approximations as described in  and numerical simulations.

To construct a local structure approximation, it is more convenient to represent configurations of cars as binary sequences, where zeros represent empty spaces and ones represent cars. Since for $`v_{max}=2`$ the speed of a car is determined by the states of, at most, two sites in front of it, the minimal block size to obtain nontrivial results is 3 (the site occupied by a car plus two sites in front of it). In what follows, we limit our attention to order-3 local structure approximation.

Using 3-block probabilities, we can write a set of equations describing the time evolution of these probabilities

$$P_{t+1}(b_2b_3b_4)=\underset{\genfrac{}{}{0pt}{}{a_i\{0,1\}}{i=0,1,\mathrm{},6}}{}w(b_2b_3b_4|a_0a_1a_2a_3a_4a_5a_6)P_t(a_0a_1a_2a_3a_4a_5a_6),$$
(8)   
where $`P_t(b_2b_3b_4)`$ is the probability of block $`b_2b_3b_4`$ at time $`t`$,   
and $`w(b_2b_3b_4|a_0a_1a_2a_3a_4a_5a_6)`$ is the conditional probability that the rule maps the seven-block $`a_0a_1a_2a_3a_4a_5a_6`$ into the three-block $`b_2b_3b_4`$. Letters $`a`$ denote the states of lattice sites at time $`t`$, while $`b`$ denote states at time $`t+1`$, so that, for example, $`a_3`$ is the state of site $`i=3`$ at time $`t`$, and $`b_3`$ is the state of the same site at time $`t+1`$. Conditional probabilities $`w`$ can be easily computed from the definition of the rule, although, since their number is quite large, we used a computer program for their determination.

Equation (8) is exact. The approximation consists of expressing the seven-block probabilities in terms of three-block probabilities using Bayesian extension. That is,

$`P_t(a_0a_1a_2a_3a_4a_5a_6)=`$   
$`\text{ }{\displaystyle \frac{P_t(a_0a_1a_2)P_t(a_1a_2a_3)P_t(a_2a_3a_4)P_t(a_3a_4a_5)P_t(a_4a_5a_6)}{P_t(a_1a_2)P_t(a_2a_3)P_t(a_3a_4)P_t(a_4a_5)}},`$   
where $`P_t(a_ia_{i+1})=P_t(a_ia_{i+1}0)+P_t(a_ia_{i+1}1)`$ for $`i=1,\mathrm{},4`$. Equations (8) and (3) define a dynamical system whose fixed point approximates three-block probabilities of the limit set of the CA rule. Due to the nonlinear nature of these equations it is not possible to find the fixed point analytically. We have determined it numerically.

From the knowledge of three-block probabilities $`P_t(a_2a_3a_4)`$, the car density is given by

$$\rho =P_t(1)=P_t(100)+P_t(101)+P_t(110)+P_t(111),$$
(10)   
and the average velocity by

$$v=\frac{2(1p)P_t(100)+pP_t(100)+(1p)P_t(101)}{\rho }.$$
(11)   
This last result gives the expression of the order parameter $`m=v_{max}v`$.

Note that $`\rho =P_t(1)`$ is a constant of motion for equation (8) whose fixed point depends, therefore, upon the initial car density $`P_0(1)=P_0(100)+P_0(101)+P_0(110)+P_0(111)`$. Consequently, using order-3 local structure approximation allows to determine the order parameter $`m`$ as a function of car density, and the critical exponents.

For most simulations, we used a lattice size equal to 1000 and our results are average of 1000 runs of 1000 iterations. For $`p=0.0005`$ we took a lattice of 10000 sites and averaged 500 runs of 10000 iterations.

## 4 Results

The susceptibility $`\chi _\rho `$ at constant $`\rho `$, defined by

$$\chi _\rho =\underset{p0}{lim}\frac{m}{p},$$
(12)   
cannot be calculated exactly. So we determined its value using local structure approximation as described in the preceding section and computer simulations for $`v_{\mathrm{max}}=2`$. Figures 1a and 1b show the results obtained from order-3 local structure approximation<sup>1</sup><sup>1</sup>1If $`v_{\mathrm{max}}>2`$, one should use $`(v_{\mathrm{max}}+1)^{\mathrm{th}}`$-order local structure approximation and simulations, respectively. In the limit $`p0`$, the susceptibility diverges as $`(\rho _c\rho )^\gamma `$ for $`\rho <\rho _c`$, and as $`(\rho \rho _c)^\gamma ^{}`$ for $`\rho >\rho _c`$. Using local structure approximation (Figures 2a and 3a), we found

$$\gamma =0.91\pm 0.03\text{and}\gamma ^{}=0.98\pm 0.03$$

while simulations (Figures 2b and 3b) yield

$$\gamma =0.86\pm 0.05\text{and}\gamma ^{}=0.94\pm 0.05.$$

Another exponent of interest is $`\delta `$. It characterizes the behavior of $`m`$ as a power of $`p`$ for $`\rho =\rho _c`$. Here again we have determined the value of

$$\underset{p0}{lim}\frac{m(\rho _c,0)m(\rho _c,p)}{p}$$

using order-3 local structure approximation and simulations. Our results are represented in Figure 4, and the values of $`\delta `$, obtained using order-3 local structure approximation and computer simulations, are, respectively, given by

$$1/\delta =0.51\pm 0.01\text{and}1/\delta =0.53\pm 0.02$$

In is interesting to note that the values

$$\beta =1,\gamma 1,\delta 2$$

obtained for the critical exponents are found in equilibrium statistical physics in the case of second-order phase transitions characterized by nonnegative order parameters above the upper critical dimensionality.

Close to the phase transition point, critical exponents obey scaling relations. If we assume that, in the vicinity of the critical point $`(\rho =\rho _c,p=0)`$, the order parameter $`m`$ is a generalized homogeneous function of $`\rho \rho _c`$ and $`p`$ of the form

$$m=|\rho \rho _c|^\beta f\left(\frac{p}{|\rho \rho _c|^{\beta \delta }}\right),$$
(13)   
where the function $`f`$ is such that $`f(0)0`$, then, differentiating $`f`$ with respect to $`p`$ and taking the limit $`p0`$, we readily obtain

$$\gamma =\gamma ^{}=(\delta 1)\beta ,$$
(14)   
which is verified by our numerical values obtained either using order-3 local structure approximation or computer simulations. Figures 5 and 6, clearly confirm the existence of a universal scaling function.

## 5 Conclusion

In this short note, we have derived the critical behavior of a CA traffic flow model using an order parameter breaking the symmetry of the jam-free phase. Random braking, which is thought to be an essential ingredient of any CA traffic flow model, appears to be the symmetry-breaking field conjugate to the order parameter. For $`v_{\mathrm{max}}=2`$, we have determined the critical exponents $`\beta `$, $`\gamma `$ and $`\delta `$ using order-3 local structure approximation and computer simulations. These critical exponents satisfy a scaling relation which can be derived assuming that the order parameter is a generalized homogeneous function of $`|\rho \rho _c|`$ and $`p`$ in the vicinity of the phase transition point.

## 6 Acknowledgements

One of the authors (H. F.) acknowledges financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada. He also wishes to thank The Fields Institute for Research in Mathematical Sciences and the Department of Mathematics and Statistics of University of Guelph for generous hospitality.