format: text

text:
# Heat conduction in a one-dimensional gas of elastically colliding particles of unequal masses

## Abstract

We study the nonequlibrium state of heat conduction in a one-dimensional system of hard point particles of unequal masses interacting through elastic collisions. A BBGKY-type formulation is presented and some exact results are obtained from it. Extensive numerical simulations for the two-mass problem indicate that even for arbitrarily small mass differences, a nontrivial steady state is obtained. This state exhibits local thermal equilibrium and has a temperature profile in accordance with the predictions of kinetic theory. The temperature jumps typically seen in such studies are shown to be finite-size effects. The thermal conductivity appears to have a very slow divergence with system size, different from that seen in most other systems.

Introduction: The problem of finding a one-dimensional system of interacting particles, evolving through Newtonian dynamics, in which Fourier’s law of heat conduction holds, has attracted considerable interest in recent years . One of the first studies in the field was the work of Rieder et al  who obtained the exact nonequilibrium steady state (NESS) of a chain of coupled harmonic oscillators. They found a constant temperature profile in the system and a heat current that was independent of system size. This result was not too surprising, considering that the harmonic chain is integrable and there is no scattering of phonons. Since then, a large number of studies have looked at the effects of introducing impurities, nonlinearity and external potentials \[$`38`$\]. The specific question asked in most of these studies is the following. Defining the thermal conductivity $`K`$ of a system through the linear response formula $`j=K(T)T`$, where $`j`$ is the the current and $`T`$ the local temperature, is there a one-dimensional model which gives a finite $`K`$ ? Contrary to initial expectations, it was found that adding impurities and/or nonlinearity into a system did not result in finite $`K`$. The first model in which a finite $`K`$ was found was the so-called ding-a-ling model  in which alternate particles on the line are bound to harmonic springs and the other particles move freely between pairs of the bound particles. Numerical studies of several other models have also given finite conductivity . A common feature in these models is the presence of external potentials which lead to momentum nonconservation. Recently it has been proved rigorously  that the conductivity as given by the Green-Kubo formula always diverges in one-dimensional momentum conserving systems with finite pressure. Two recent studies  have reported finite conductivity in momentum-conserving systems but both of these have vanishing pressure and so there is no contradiction.

While the proof of anamalous conductivity in momentum conserving systems has been a significant progress, there remain many issues that still need to be addressed. First of all, the proof uses the Kubo formula and this is not fully satisfactory since the validity of the Kubo formula, even in the limit of arbitrarily small temperature gradients, has never been established rigorously. It may be noted that any derivation of the Kubo formula for thermal conductivity essentially makes the following assumptions: (i) the NESS is in local thermal equilibrium (LTE) (ii) that Fourier’s law is valid and (iii) regression of equilibrium fluctuations occur in the same way as nonequilibrium processes . These assumptions are physically motivated but may not be true in all cases. Secondly we note that while the focus of most of the work on heat conduction has been addressed to obtaining Fourier’s law, the more general problem is one of understanding the nonequilibrium energy current carrying state of a many-body system. For example one important question is the existence or otherwise of LTE in the steady state. Thus one would like to know if it is possible to define a slowly varying temperature field which determines all other local properties in the system. This point has not attracted much attention, even though it is quite crucial even for stating Fourier’s law and using results of linear response theory. Naively one might expect that if thermal currents in the system vanish in the thermodynamic limit, then LTE should hold but this has been shown not to be true always . Another interesting question is the determination of the temperature profile itself. These questions are clearly of interest and can be asked independently of whether or not Fourier’s law holds. Finally we note that numerical studies of heat conduction are problematic for several reasons. One needs long time numerical solution of nonlinear differential equations which is time-consuming and not very accurate. This, in addition to long equilibration times typically occuring in such systems, restricts one to small system sizes. Also the treatment of boundary related problems, such as that of temperature jumps, is not straightforward. Hence it has often been difficult to arrive at correct conclusions from numerical studies and it is desirable to have more accurate studies.

In this paper we study heat conduction in a system of hard elastic particles of unequal masses moving on a one-dimensional line. The only interaction between the particles is through elastic collisions. This model was first considered by Casati  as a possible candidate for obtaining Fourier’s law but the numerical results were insufficient to draw any definite conclusions. More recently this model has been studied by Hatano  who obtained a diverging conductivity. The equal-mass case with dissipation was studied by Du et al  who obtained a rather surprising NESS which implied a breakdown of usual hydrodynamics.

Here we present extensive and accurate numerical work on this model and also an analytic formulation of the BBGKY-type. Our aim has been to give a more detailed characterization of the NESS than has been previosly done. The present model is particularly suitable for this purpose for two reasons: (i) simulations of this model do not require a numerical solution of nonlinear differential equations and it is possible to obtain very accurate results, (ii) analytically the BBGKY hierarchy has a comparatively simpler structure and some exact statements can be made. The most interesting result obtained is that the steady states for the case of equal masses and the case with arbitrarily small mass differences are completely different. In the thermodynamic limit the latter case exhibits LTE and the temperature profile approaches a form predicted by kinetic theory for a one-dimensional gas.This is surprising since kinetic theory predicts a finite conductivity with a $`T^{1/2}`$ dependence on temperature. On the other hand our model is momentum conserving and the proof for diverging Kubo conductivity holds. In our finite size studies we find a slow divergence of the conductivity ($`L^\alpha `$ with $`\alpha <0.2`$). Our work also clarifies some of the problems related to boundary effects. The jumps in temperature at the boundaries are shown to be finite-size effects which are studied systematically.

Definition of model: We consider $`N`$ point particles numbered $`i=1,\mathrm{}N`$ moving in a one-dimensional box extending from $`0`$ to $`L`$. The mass, position and velocity of the $`i`$th particle are denoted by $`m_i`$, $`x_i`$ and $`u_i`$. The only interaction between particles is through elastic collisions. After a collision between particle $`i`$ and $`(i+1)`$, the new velocities are obtained from momentum and energy conservation and given by the linear equations:

$`u_i^{}={\displaystyle \frac{m_im_{i+1}}{m_i+m_{i+1}}}u_i+{\displaystyle \frac{2m_{i+1}}{m_i+m_{i+1}}}u_{i+1}`$ (1)
$`u_{i+1}^{}={\displaystyle \frac{2m_i}{m_i+m_{i+1}}}u_i{\displaystyle \frac{m_im_{i+1}}{m_i+m_{i+1}}}u_{i+1}.`$ (2)
Between collisions the particles travel with constant velocity. The coupling to heat baths is implemented by using Maxwell boundary conditions. Thus whenever a particle of mass $`m`$ collides with a wall at temperature $`T`$, it is reflected back with a velocity chosen from the distribution $`P(u)=(mu/T)exp(mu^2/(2T))`$. The temperatures at the two ends are taken to be $`T_1`$ and $`T_2`$.

The case when all particles have equal masses can be solved easily and behaves similarly to the ordered harmonic springs case . The “temperature” profile is flat (with the value $`\sqrt{T_1T_2}`$), current is independent of system size and there is no LTE. The equal-mass problem is integrable and essentially reduces to a single-particle problem and so the results are not surprising. As soon as the masses are made different, the system becomes nonintegrable and is expected to have good ergodicity properties , and correspondingly the NESS should be very different.

BBGKY-type equations: In principle, a complete solution of the heat conduction problem could be obtained from the steady-state solution of the master equation for evolution of the $`N`$-point distribution equation $`\rho (\{x_l,u_l\},t)`$. In practice this is difficult and a simpler approach is to work with the so-called BBGKY hierarchy which deals with reduced distribution functions . Let us first make the following definitions: $`\rho _l(x,u,t)=\delta (xx_l)\delta (uu_l)`$, $`\rho _{l,l+1}(x_1,u_1;x_2,u_2)=\delta (x_1x_l)\delta (u_1u_l)\delta (x_2x_{l+1})\delta (u_2u_{l+1})`$, where $`A=\rho (\{x_l,u_l\},t)A(\{x_l,u_l\})_ldx_ldu_l`$. Further we define $`p_{l,l+1}(x,u_1,u_2)`$ as the number of collisions per unit time occuring at $`x`$ between the $`l`$th and $`(l+1)`$th particles with respective velocities $`u_1`$ and $`u_2`$. Clearly $`p_{l,l+1}(x,u_1,u_2)=\rho _{l,l+1}(x,u_1;x,u_2)(u_1u_2)\theta (u_1u_2)`$. In terms of these the BBGKY-type equations relating one-point functions to two-point functions are the following:

$`\rho _l(x,u,t)/t+u\rho _l(x,u,t)/x=`$ (3)
$`{\displaystyle p_{l1,l}(x,u_1,u_2)\delta (u_2^{}u)𝑑u_1𝑑u_2}{\displaystyle p_{l1,l}(x,u_1,u)𝑑u_1}+`$ (4)
$`{\displaystyle p_{l,l+1}(x,u_1,u_2)\delta (u_1^{}u)𝑑u_1𝑑u_2}{\displaystyle p_{l,l+1}(x,u,u_2)𝑑u_2}.`$ (5)
These equations hold for $`x0,L`$. At the boundaries, the distribution functions satisfy appropriate boundary conditions.

The physical observables that we will be interested in are the particle density $`n(x,t)`$, the energy density $`ϵ(x,t)`$ and the energy current density $`j(x,t)`$. These can be expressed in terms of the one-point functions $`\rho _l(x,u,t)`$. Thus we have:

$`n(x,t)`$ $`=`$ $`{\displaystyle \underset{l}{}}\delta (xx_l)={\displaystyle \underset{l}{}}{\displaystyle \rho _l(x,u,t)𝑑u}`$ (6)
$`ϵ(x,t)`$ $`=`$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{m_lu_l^2}{2}}\delta (xx_l)={\displaystyle \underset{l}{}}{\displaystyle \frac{m_l}{2}}{\displaystyle u^2\rho _l(x,u,t)𝑑u}`$ (7)
$`j(x,t)`$ $`=`$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{m_lu_l^3}{2}}\delta (xx_l)={\displaystyle \underset{l}{}}{\displaystyle \frac{m_l}{2}}{\displaystyle u^3\rho _l(x,u,t)𝑑u}`$ (8)
The temperature field $`T(x)`$ is defined as $`T(x,t)=2ϵ(x,t)/n(x,t)`$. Our first result is the following exact current conservation equation: $`ϵ(x,t)/t+j(x,t)/x=0`$. This result is easily obtained by multiplying Eq. (5) by $`m_lu^2/2`$, integrating over $`u`$, and summing over all $`l`$. There is a pairwise cancellation of all terms on the right-hand side. Similarly by multiplying Eq. (5) by $`m_lu/2`$, integrating and summing gives, in the steady state: $`ϵ(x)/x=0`$. Thus for any choice of masses $`\{m_l\}`$, the energy density in the steady state is constant in space. Physically the constancy of energy density follows from the constancy of pressure and the linear relation between the two quantites in an ideal gas. However this does not imply a constant temperature profile since the temperature also depends on the number density $`n(x)`$ which is not constant.

From the fact that the dynamics \[Eq. (2)\] is invariant under a constant scaling of the masses, it follows that the temperature profile does not change under $`m_i\nu m_i`$. Also from the boundary conditions it is easily shown that $`T(\nu T_1,\nu T_2,x)=\nu T(T_1,T_2,x)`$. From now on, we shall consider the dimer case where we consider only particles of two different masses $`m_1`$ and $`m_2`$ placed alternately on the line. Because of the above scaling properties the only independent variables are the ratios, $`m_2/m_1`$ and $`T_2/T_1`$, and $`N`$. We will henceforth consider the case $`m_1=1`$, $`m_2=m`$, $`T_1=2`$ and $`T_2=8`$.

Numerical results: In our numerical simulations we let the system evolve with the appropriate boundary conditions and compute time-averages of various quantities in the steady state. The time evolution does not require numerical solution of differential equations since the exact solution is essentially known. The system is evolved by computing successive collision times and updating velocities using the collision rules Eq. (2). The only errors are those due to round-off. We have verified that the simulations reproduce all the exact known results both for the equal and the unequal mass cases.

In our simulations we vary $`\delta =m1`$ and $`N`$. The mass values $`\delta =0.078,0.11,0.22,0.44`$ were studied for particle numbers $`N=41,81,161,321,641`$ and $`1281`$. The size of the boxes were changed so that in all cases the average density of particles was fixed at $`2`$. The number of particles is chosen to be odd so that at both ends the particles in contact with the bath have the same mass. The number of collisions over which the averaging is done was between $`10^910^{10}`$. In all cases we checked that increasing the time of averaging by a factor of $`10`$ did not significantly change the data.

In Fig. (1), we plot the steady state temperature profiles at different system sizes for $`\delta =0.22`$. The temperature has a smooth and continuously varying profile. There is a jump at the boundaries which decreases as the system size is increased and is expected to vanish in the thermodynamic limit when the current also becomes vanishingly small. We find that this is true for any non-zero $`\delta `$. For smaller $`\delta `$ one needs to go to larger system sizes to get the same temperature profile. Infact for small $`\delta `$ and large $`N`$ the temperature-profile depends on $`\delta `$ and $`N`$ only through the scaling combination $`\delta ^2N`$. This is illustrated in Fig. (2) where we plot data corresponding to five different values of $`\delta `$, each with a different $`N`$, chosen such that $`\delta ^2N`$ is the same. Note that for $`\delta =0`$, the temperature profile (which is flat and given by $`T=\sqrt{T_1T_2}`$) and the energy current are both independent of system size; thus $`\delta =0`$ is a singular point, while the steady state in the limit $`\delta 0,N\mathrm{}`$ with $`\delta ^2N`$ constant is quite different and nontrivial.

For fixed value of $`\delta `$, as we increase $`N`$ the temperature profile approaches a limiting form. Quite amazingly this limiting form seems to be exactly one that would be predicted by kinetic theory. We recall that kinetic theory for a one-dimensional gas predicts Fourier behaviour with $`KT^{1/2}`$ and hence a nonlinear temperature profile $`T_k(x)=[T_1^{3/2}(1x/L)+T_2^{3/2}x/L]^{2/3}`$. This has been plotted in Fig. (1). We find that the following scaling form gives a reasonable collapse of our data:

$`T(x,N,\delta )=T_k(x)+{\displaystyle \frac{1}{(\delta ^2N)^\gamma }}g(x).`$ (9)
The inset in Fig. (1) shows the collapse of data for $`\delta =0.22`$ obtained using the above scaling form with $`\gamma =0.67`$.

We now look at the dependence of $`K`$ on system size. In Fig. (3) we plot $`j`$ versus $`L`$ for $`\delta =0.22`$. It is clear from the data that the conductivity which is proportional to $`Lj`$ has a slow divergence given by $`KL^\alpha `$ with $`\alpha <0.2`$. We note that this is significantly different from the system size dependence of $`KL^{0.4}`$ found by Hatano in the same model  and also in other momentum-conserving systems like FPU and the diatomic Toda .

Finally we have checked for LTE: to do so we compute the steady state expectation $`u^{(4)}(x)=_lm_lu_l^4\delta (xx_l)`$. If there was LTE, this quantity would be determined by the local temperature $`T(x)`$. In Fig. (4) we plot $`u^{(4)}(x)`$, as determined directly by taking time averages and also the value predicted from the local temperature $`T(x)`$. We see that at large system sizes LTE is indeed achieved. We also find that for smaller values of $`\delta `$, one needs to go to larger system sizes to get LTE.

In summary we have studied heat conduction in the unequal mass problem which appears to be the simplest nontrivial deterministic system, in one dimension, for which a very detailed investigation of the NESS can be made. Our study shows that a meaningful hydrodynamic description of the steady state is possible even in a situation where (presumably) $`K\mathrm{}`$ in the thermodynamic limit. It is clear that further studies of this model can throw much light on the difficult problem of transition from the microscopic to macroscopic description in the context of nonequilibrium phenomena.

I thank D. Dhar and M. Rao for many valuable suggestions. I also thank J. Das, C. Dasgupta, R. Pandit, S. Ramaswamy and B. S. Shastry for discussions.

Also at the Poornaprajna Institute, Bangalore. dabhi@rri.res.in


text:
# Random walk theory of jamming in a cellular automaton model for traffic flow

## 1 Introduction

One of the most interesting phenomena of traffic flow is the appearance of traffic jams. Two classes of jams can be distinguished: Those induced by external influences, e.g., bottlenecks, lane reductions or intersections , and spontaneous jams, sometimes called phantom jams, caused by velocity fluctuations. The latter effect was first shown by Treiterer  who examined a series of aerial photographs of a highway. Further measurements on real traffic by Kerner and coworkers  have revealed the existence of phase separated jammed states and homogeneous metastable states with a very high throughput. Experimentally they observed several characteristic features of wide jams. Wide jams are compact stable structures which are clearly separated from moving vehicles. The outflow of a jam is reduced remarkably compared to the maximal flow in homogeneous traffic. The upstream velocity of the jam front is approximately constant about $`15km/h`$, independent of the road conditions. These quantities are nowadays regarded as important parameters of highway traffic which can be used for example to calibrate theoretical models. In addition to homogeneous states and wide traffic jams further phases of vehicular traffic on highways were found in recent years. Empirical observations  showed the existence of so-called synchronized traffic states and stop and go traffic. Recently Neubert et al.  analyzed single-vehicle data in detail. They showed that the cross-correlations between density and flow can be used as a criterion for the identification of the different traffic states. It should be mentioned that the nature of these traffic states is still under a vivid debate.

Besides the analysis of experimental data various theoretical approaches for the description traffic flow theory have been suggested in recent years . These approaches can be divided into fluid-dynamical descriptions and microscopic models where individual vehicles are represented by a particle. The first model with the ability to describe the formation and evolution of wide jams accurately was a macroscopic approach proposed by Kerner and Konhäuser . Besides these fluid-dynamical models several microscopic approaches  that are able to reproduce the dynamics of wide jams have been discussed. Nowadays one of the main interests in applications to real traffic is to perform real-time simulations of large networks with access to individual vehicles. Therefore cellular automata (CA) models have become quite popular for traffic flow simulations.

Due to their simplicity these microscopic models can be used very efficiently for computer simulations. A few years ago, Nagel and Schreckenberg (NaSch) proposed a probabilistic CA model for traffic flow . This model is the simplest known CA that can reproduce the basic phenomena encountered in real traffic, e.g., the occurance of phantom traffic jams. On the other hand, the NaSch model can not explain all experimental results. Therefore modifications have been suggested to obtain a description on a more detailed level. Here we concentrate our investigations on the NaSch model with velocity-dependent randomization (VDR model ). It exhibits metastable states and large phase separated jams whose upstream front velocities can be easily adjusted through its randomization parameters.

Recently Knospe et al.  have proposed a high fidelity cellular automaton model for highway traffic which is able to reproduce empirical single-vehicle data . In connection to the work presented in this paper it is interesting to note that the formation of wide jams in the model of  is implemented using VDR-like mechanisms. An additional interesting application of the VDR model, or models exhibiting metastability in general, is the flow optimization by stabilizing homogeneous states of traffic at bottlenecks . Due to the wide field of applications and the fact that the stochastic dynamics of traffic jams is similar in the different CA approaches there is a need for an exact understanding of the traffic jams that occur in the models or, at least, for good approximations.

## 2 NaSch model with velocity-dependent randomization VDR

The description of a physical system in terms of a cellular automata (CA) approach is in general an extreme simplification of the real world conditions. In the CA models for traffic the space, speed, acceleration and even time are treated as discrete variables. The motion of vehicles is realized through a simple set of rules. In the spirit of modeling complex phenomena in statistical physics Nagel and Schreckenberg (NaSch) have chosen a minimal set of rules for their model . The aim of the NaSch model is to describe basic phenomena of traffic flow and not to be very accurate on a microscopic level. Therefore additional rules have to be added or the standard set of rules has to be modified for a proper modeling of the fine-structure of traffic flow (see  and references therein).

In the NaSch model the road is divided into cells of length $`7.5m`$. Each cell can either be empty or occupied by at most one car. The speed $`v_n`$ of each vehicle $`n=1,2,\mathrm{},N`$ can take one of the $`v_{\mathrm{max}}+1`$ allowed integer values $`v_n=0,1,\mathrm{},v_{\mathrm{max}}`$. The state of the road at time $`t+1`$ can be obtained from that at time $`t`$ by applying the following rules for all cars at the same time (parallel dynamics):

* Step 1: Acceleration:
$`v_n\mathrm{min}(v_n+1,v_{\mathrm{max}})`$
* Step 2: Braking:
$`v_n\mathrm{min}(v_n,d_n1)`$
* Step 3: Randomization with probability $`p`$:
$`v_n\mathrm{max}(v_n1,0)`$
* Step 4: Driving:
$`x_nx_n+v_n`$

Here $`x_n`$ denotes the position of the $`n`$-th car and $`d_n=x_{n+1}x_n`$ the distance to the next car ahead. The density of cars is given by $`\rho =N/L`$, where $`L`$ is the length of the system, i.e. the number of cells. One time-step corresponds to approximately $`1s`$ in real time.

The NaSch model neither exhibits metastable states and hysteresis nor wide phase-separated jams. However, a simple generalization exists which is able to reproduce these effects. In the so-called VDR model  a velocity-dependent randomization (VDR) parameter $`p=p(v(t))`$ is introduced in contrast to the constant randomization in the original formulation. This parameter has to be determined before the acceleration step 1. One finds that the dynamics of the model strongly depends on the randomization. For simplicity we study the so-called slow-to-start case

$$p(v)=\{\begin{array}{ccc}p_0& \mathrm{for}& v=0,\\ p& \mathrm{for}& v>0,\end{array}$$
(1)
with two stochastic parameters $`p_0`$ and $`p`$ which already contains the expected features. In this case the randomization $`p_0`$ for standing cars is much larger than the randomization $`p`$ for moving cars, i.e., $`p_0p`$. An interesting fact is that an alternative choice of $`p(v)`$, e.g., $`p_0p`$, leads to a completely different behavior. Note that for $`p_0=p`$ the original NaSch model is recovered. Throughout the paper we will always assume that $`pp_0`$ if not stated otherwise.

Typical fundamental diagrams of the VDR model (see Fig. 1) show a density regime $`\rho _1<\rho <\rho _2`$ where the flow can take on two different states depending on the initial conditions. The homogeneous branch is metastable with an extremely long life-time and the jammed branch shows phase separation between jammed and free flowing cars. Neglecting interactions among free flowing vehicles the fundamental diagram of the model can be derived on the basis of heuristic arguments in good agreement to numerical results . The microscopic structure of the jammed states in the VDR model differs from those found in the NaSch model. While jammed states in the NaSch model contain clusters with an exponential size distribution , one finds phase separation in the VDR model. The reason for this different behavior is the reduction of the outflow of a jam compared to the maximal possible global flow. A large stable jam can only exist if the outflow from the jam is equal to the inflow. If the outflow is the maximal possible flow a stable jam can only exist at the corresponding density $`\rho _{\mathrm{max}}`$, but will easily dissolve due to fluctuations. If the outflow of a jam is reduced, the density in the free flow regime is smaller than the density $`\rho _{\mathrm{max}}`$ of maximal flow where interactions between vehicles lead to jams so that cars can propagate freely (for $`pp_0`$) through the low density part of the road. Therefore no spontaneous jam formation is observable. It is obvious that no phase separation can occur in the NaSch model due to the fact, that the outflow of a jam is maximal. Figure 1 shows the typical structure of the jammed state in the VDR model.

## 3 Jamming process

The fundamental diagram of the VDR model with periodic boundary conditions can be divided into three different regimes according to the jamming properties. For densities up to $`\rho _1`$ no jams with long lifetime appear and jams existing in the initial conditions dissolve very quickly. Obviously in the regime $`0\rho \rho _1`$ the outflow of a jam is greater than the inflow. This behavior has to be contrasted to the one found for densities above $`\rho _2`$. Here no homogeneous state without jams can exist. The most interesting regime lies between the two densities $`\rho _1`$ and $`\rho _2`$ where the system can be in two different states. One is a metastable homogeneous state with an extremely long life time where jams can appear due to internal fluctuations. Note that the homogeneous states can be destroyed by external perturbations, e.g., by stopping cars. The other state is a phase separated state with large jams which can be reached through the decay of the homogeneous state or directly owing to the initial condition. The microscopic structure of a spontaneous emerging jam can be seen in the space time plot in Fig. 1. The origin of the wide jam in the initially homogeneous state is a local velocity fluctuation that leads to a stopped car. A velocity fluctuation of a car can force the following car to brake if its headway becomes too small. If the local density is large enough this can lead to a chain reaction which finally forces a car to stop. Such local velocity fluctuations determine the typical density-dependent lifetime of the homogeneous states.

As one can see in Fig. 1, the density in the outflow regime of the jam is reduced compared to the average density. Therefore the jam length is growing approximately linear until outflow and inflow coincide due to the periodic boundary conditions. It is quite evident that interactions between vehicles in the outflow region of a jam are negligible so that no spontaneous jams can appear. A recent study by Pottmeier et al.  indicates that local defects can break the strong phase separation found in the VDR model in such a way that stop and go traffic can exist.

Furthermore Fig. 1 shows that strong fluctuations in the upstream jam front appear due to the fact that the outflow of a jam is determined by a stochastic parameter. As aforementioned the locally emerged jam will probably grow for a certain time because the mean inflow into the jam is larger than the mean outflow. At this point it should be stressed, however, that although the inflow is larger than the outflow these quantities are stochastic and therefore even in this growth regime a complete dissolution of the emerged jam is possible through fluctuations. The dynamics of this growth (dissolution) process is the main object of our investigations in this paper. Nonetheless, assuming that the locally emerged jam does not dissolve it will grow a certain time $`t_{\mathrm{stat}}`$ until the mean inflow and outflow are equal due to periodic boundary conditions. The average jam length then strongly fluctuates. This can also lead to a complete dissolution of the jam when these fluctuations are of the order of the jam length. The phase separation in the jammed state can be identified directly using the jam-gap distribution . It was shown that for $`p=0`$ the jam is compact in the sense that no holes between the jammed cars appear. For $`p>0`$ holes are formed due to velocity fluctuations of vehicles entering the jam. Nevertheless, it can be assumed that for $`p_0p`$ the jammed states are phase separated, i.e. the size of the jam is of the order of the system size.

A schematic representation of a single jam is depicted in Fig. 2 where jammed cars are represented by grey cells and empty cells are white. The parameter $`\alpha `$ denotes the outflow and is equal to $`p_0`$ in the studied model. Since for $`pp_0`$ car-car interactions can be neglected in the outflow region a car that leaves a jam can be treated as escaped for all times. Therefore the velocity of the upstream jam front is only determined by $`p_0`$, i.e., the waiting time of the first car. Neubert et al.  analyzed the jam velocity in the VDR model with an autocorrelation function and find an excellent agreement with the simple assumption made above. Note that this behavior is quite different from that of the NaSch model where even in the outflow region spontaneous jams occur for sufficiently large $`p`$ (stop and go traffic). The parameter $`\beta `$ is a measure for the inflow into the jam. This quantity is determined by the vehicle distribution on the road section upstream the jam (see Fig. 2).

In the next section we will introduce a stochastic theory of jamming based on random walk like arguments. Our approach gives the exact solution in the special case $`p=0`$ where the fluctuations of free flowing vehicles are completely suppressed. For the general case $`0<pp_0`$ the results are in good agreement with simulations.

## 4 Random walk theory of jamming

The main object of our investigation is a sequence of $`n`$ cars at rest forming a compact jam (see Fig. 2). In every timestep the first car can leave the megajam with probability $`\alpha =p_0`$ according to the acceleration and randomization steps in the update algorithm (see Sec. 2). Additionally a new car is able to enter the jam at its end with probability $`\beta `$ which is determined by the vehicle distribution behind the jam. In general $`\beta `$ will be generated by the outflow of a megajam (see below).

The outflow $`\alpha `$ realized through $`p_0`$ is an independent identically distributed (i.i.d.) random variable. For simplicity we start with the case that $`\beta `$ is also an i.i.d. random variable. This condition is fulfilled if fluctuations of free flowing vehicles are completely suppressed, i.e., for $`p=0`$, and when the vehicle distribution behind the considered jam is generated by the outflow of a megajam. In this case we will consider open boundary conditions with an (infinite) megajam at the left end and an empty system at the right end (Fig. 2). The density in the free flow regime is then determined by the waiting time of the first car in the megajam. In this special case the gaps between the vehicles are always a multiple of their maximum velocity $`v_{\mathrm{max}}`$ since the waiting times are discrete and no interactions between cars occur in the outflow region for $`p=0`$. For example a waiting time of $`t_w`$ timesteps leads to a gap of $`t_wv_{\mathrm{max}}`$ cells. Then a car with this gap will reach the end of the jam exactly $`t_w`$ timesteps after its predecessor. It is obvious that the waiting times of the first car in the megajam determine directly the inflow $`\beta `$ into the analyzed jam which therefore can be treated as an i.i.d. variable. Note that the randomization parameter $`\stackrel{~}{p}_0`$ of cars standing in the megajam can be considered as a control parameter for the vehicle distribution in the free flow regime, i.e., the inflow $`\beta =\beta (\stackrel{~}{p}_0)`$. Therefore it is in general different from the $`p_0`$ of the cars that have already left the megajam, i.e., the $`p_0`$ governing the resolution of the jam under consideration. In this way we can realize $`\alpha \beta `$ for the open system. The case of periodic boundary conditions can be treated as a special case of the problem where in the jammed state the inflow into the jam is determined by the outflow out of the jam and hence $`\alpha =\beta `$. A detailed description of the initial and boundary conditions considered here as well as a discussion of the influence of deviating gap distributions in the inflow region (e.g. if $`\beta `$ is no longer an i.i.d. variable) will be given below.

In the following we will map the jam dynamics onto a random walk problem. The number of standing cars $`n`$ in the jam determines the position $`n`$ of a random walker. The walker moves on a discrete lattice in discrete time. A car leaving (entering) the jam then corresponds to one step to the left (right). In the following we determine the probability $`\pi _{t,n}`$ that a jam of width $`n`$ resolves after $`t`$ timesteps. Here a jam is considered to have resolved when the last remaining car accelerates.

In random walk terminology this problem is equivalent to the calculation of the first passage time of a walker starting at position $`n_0=n(t=0)`$. $`\pi _{t,n}`$ is the probability that walker at position $`n`$ reaches the origin $`n=0`$ of the system in $`t`$ timesteps. Taking into account that $`\alpha `$ and $`\beta `$ are i.i.d. random variables one gets the following master equation for the duration of the process until the random walker reaches the origin:

$`\pi _{t+1,n>1}`$ $`=`$ $`\alpha (1\beta )\pi _{t,n1}+\beta (1\alpha )\pi _{t,n+1}+\left[(1\alpha )(1\beta )+\alpha \beta \right]\pi _{t,n},`$
$`\pi _{t+1,1}`$ $`=`$ $`\alpha \pi _{t,0}+\beta (1\alpha )\pi _{t,2}+(1\alpha )(1\beta )\pi _{t,1}.`$ (2)
It is obvious that $`\pi _{0,n_0>0}=0`$ and that by definition $`\pi _{1,1}=\alpha `$ since a jam of length $`n=1`$ resolves with probability $`\alpha `$ in one timestep independent of the inflow111This is due to the update procedure and will be discussed later.. Furthermore we assume an absorbing barrier at the origin $`\pi _{0,0}=1`$, i.e. the process stops when the random walker reaches the origin (the jam is dissolved). Note, that the variable $`n`$ covers the whole spektrum of possible positions of a random walker during his movement, while $`n_0`$ denotes the starting position. For a discussion of the mathematical aspects of first passage time problems we refer to .

As an example for the derivation of these equations consider a walker starting at position $`n_0>1`$. If the first trial results in a movement to the left the process continues as if the initial position had been $`n_01`$. A movement to the left means that a jam of width $`n_0`$ evolves into a jam of width $`n_01`$. This event occurs if the first car leaves the jam and no additional car enters it. The probability for this is $`\alpha (1\beta )`$. Similarly, if the first trial results in a movement to the right the process continues as if the initial position had been $`n_0+1`$. The probability for a movement to the right is given by $`\beta (1\alpha )`$, i.e., the first car remains in the jam and an additional car enters it. Furthermore the position of the walker, and thus the jam length, are unchanged if either one car leaves and one car enters the jam (probability $`\alpha \beta `$) or no car leaves and enters the jam (probability $`(1\alpha )(1\beta )`$). From the first equation of (2) we see that the random walk is symmetric (for $`n>1`$) if $`\alpha (1\beta )=(1\alpha )\beta `$. For $`\alpha (1\beta )>(1\alpha )\beta `$ the walker is biased to the left and the jam will dissolve quickly. For $`\alpha (1\beta )<(1\alpha )\beta `$ the bias is to the right and the jam will grow on average. But even in this case a complete resolution is possible through fluctuations.

In order to determine the probabilities we introduce the following generating functions:

$$\mathrm{\Pi }_{n_0}(z)=\underset{t=0}{\overset{\mathrm{}}{}}\pi _{t,n_0}z^t.$$
(3)
At this point it has to be taken into account that the case $`n=1`$ must be viewed separately. For a jam consisting of only one standing car the probability for resolving is equal to $`\alpha `$ (the standing car accelerates) even if a new car enters the minijam at its end. Therefore we will first look at the special case of a random walker starting at position $`n_0=1`$. A typical process in which the jam resolves after $`t>1`$ timesteps may be described as follows.
(I) The walker does not move for $`\mu `$ timesteps.
(II) Then the random walker moves one unit to the right.
(III) After the step to the right the walker first has to return to $`n=1`$ before the jam resolves, i.e. before the walker reaches $`n=0`$. The return to the initial position $`n=1`$ will take $`\nu 1`$ further timesteps ($`\nu =2,3,\mathrm{}`$).
(IV) Now there are $`t\nu \mu `$ further timesteps left for the walker to reach the origin at last.
These three events are mutually independent and so the probability of the simultaneous realization of the three events is given by the product of the single probabilities.

We first look at the case $`\mu =0`$. The cases $`\mu >0`$ will later be taken into account iteratively. We start with the solution of event (III). Since in this case the position of the walker is always larger than one only the first equation of the system (2) has to be considered. In contrast to the general solution for event (III) the hopping probabilities do not depend on the position anymore. A random walker has to return to the origin (in this case $`n=1`$) starting from a position shifted one step to the right (here $`n=2`$). The solution of this homogeneous first passage time problem is described in . Assuming an absorbing barrier and introducing a separate generating function for event (III),

$$\stackrel{~}{\mathrm{\Pi }}(z)=\underset{t=0}{\overset{\mathrm{}}{}}\stackrel{~}{\pi }_{t,1}z^t,$$
(4)
where $`\stackrel{~}{\pi }_{t,1}`$ is the probability that a walker starting at $`n=2`$ reaches $`n=1`$ for the first time after $`t`$ timesteps, the following solution for part (III) of the problem can be obtained:

$`\stackrel{~}{\mathrm{\Pi }}(z)`$ $`=`$ $`{\displaystyle \frac{1\left[\alpha \beta +(1\alpha )(1\beta )\right]z}{2\beta (1\alpha )z}}`$ (5)
$`\sqrt{\left[{\displaystyle \frac{1\left[\alpha \beta +(1\alpha )(1\beta )\right]z}{2\beta (1\alpha )z}}\right]^2{\displaystyle \frac{\alpha (1\beta )}{\beta (1\alpha )}}}.`$
The return probability, i.e. the probability that the walker reaches $`n=1`$ at an arbitrary time, is then given by

$$\stackrel{~}{\mathrm{\Pi }}(1)=\{\begin{array}{cc}\frac{\alpha (1\beta )}{\beta (1\alpha )}\hfill & \mathrm{if}\alpha (1\beta )\beta (1\alpha )\hfill \\ 1\hfill & \mathrm{if}\alpha (1\beta )\beta (1\alpha ).\hfill \end{array}$$
(6)
As already discussed above, for $`\alpha (1\beta )>\beta (1\alpha )`$ the walker is biased to the left and will always reach $`n=1`$. For $`\alpha (1\beta )<\beta (1\alpha )`$, on the other hand, the walker is biased to the right and will only return with probability $`\frac{\alpha (1\beta )}{\beta (1\alpha )}<1`$.

Now the complete event (I)–(IV) occurs for some $`\nu <t`$. Summing over all possible $`\nu `$ one gets:

$`\pi _{t,1}`$ $`=`$ $`\beta (1\alpha )[\stackrel{~}{\pi }_{1,1}\pi _{t2,1}+\stackrel{~}{\pi }_{2,1}\pi _{t3,1}+\mathrm{}..+\stackrel{~}{\pi }_{t2,1}\pi _{1,1}]`$ (7)
$`+(1\alpha )(1\beta )\pi _{t1,1}.`$
The last term takes into account event (I) for $`\mu 0`$ since $`(1\alpha )(1\beta )`$ is the probability that a walker at site $`n`$ will not move. In the first term, $`\beta (1\alpha )`$ is the probability of event (II). The quantity within the brackets is the probability of events (III) and (IV) for the allowed values of $`\nu `$. Note that it is the $`(t1)`$-st term of the convolution $`\{\pi _{t,1}\}\{\stackrel{~}{\pi }_{t,1}\}`$ (see ). After multiplying (7) with $`z^t`$ and summing over all times one finds an expression for the generating function (3) for $`n_0=1`$:

$$\mathrm{\Pi }_1(z)=\frac{\alpha z}{1(1\alpha )(1\beta )z\beta (1\alpha )\stackrel{~}{\mathrm{\Pi }}(z)z}.$$
(8)
The probability for the complete return is given by $`\mathrm{\Pi }_1(1)`$. Using equation (6) it is easy to obtain this quantity explicitly. The solution for the nontrivial case $`\alpha (1\beta )<\beta (1\alpha )`$ is given by

$$\mathrm{\Pi }_1(1)=\frac{\alpha }{\beta }.$$
(9)
We now turn to the general case where one has to deal with starting positions greater than one, i.e. initial conditions consisting of more than only one standing car, $`n_0>1`$. Similarly to the foregoing approach the process of reaching the origin can be seen as the realization of mutually independent events. For instance, the probability that a random walker starting at position $`n_0=2`$ reaches the origin is given by the event that the walker reaches first position $`n=1`$ and thereafter reaches position $`n=0`$.

The general resolving process of a jam with an initial width $`n_0`$ standing cars can be described as a chain of processes leading finally to the case $`n=1`$. Thus using the generating functions (3) and (4) this convolution of events can be expressed through

$$\mathrm{\Pi }_{n_0}(z)=\stackrel{~}{\mathrm{\Pi }}(z)^{n_01}\mathrm{\Pi }_1(z).$$
(10)
With equation (6) and (9) the following relation for the resolving probability of a jam of width $`n`$ can be obtained:

$$\mathrm{\Pi }_{n_0}(1)=\frac{\alpha }{\beta }\left[\frac{\alpha (1\beta )}{\beta (1\alpha )}\right]^{n_01}.$$
(11)
Besides the resolving probability there is also another interesting quantity that is directly accessible through the generating function (10), namely the average lifetime $`T_{n_0}`$ of a jam of initial length $`n_0`$:

$$T_{n_0}=\underset{t=0}{\overset{\mathrm{}}{}}t\pi _{t,n_0}=\mathrm{\Pi }_{n_0}^{}(1).$$
(12)
Using (10) one can obtain an explicit result for $`T_{n_0}`$. Figure 3 shows results for the resolving probability (11) and for the lifetime (12) of a jam for various initial widths $`n_0`$. The outflow parameter (fluctuation parameter of standing cars) is chosen such that $`\alpha =0.5`$ in both diagrams. In the left figure a strong dependence on the starting position of the walker (initial width of the jam) can be seen. Furthermore one can observe directly an outstanding difference between the case $`n_0=1`$ and $`n_0>1`$. While for $`n_0>1`$ the resolving probability fastly converges to zero for increasing $`\beta `$, this value is shifted to $`\alpha `$ for $`n_0=1`$. This shift can be explained through the update procedure of the model discussed above. A jam consisting of only one standing car resolves with probability $`\alpha `$ even if a new car enters the minijam at its end. The right part of Fig. 3 shows the mean disolving time for different $`n_0`$. It is obvious, that this quantity grows with increasing $`n_0`$ due to the fact, that the resolving process of a jam with $`n_0`$ standing can be described as a chain of resolving processes of smaller jams. Additionally it should be noted that a higher inflow $`\beta `$ leads to lower dissolving times, but it must be taken into account that the dissolution of a jam under a high inflow $`\beta `$ is a rather rare event.

## 5 Numerical results

In the previous section we have derived an analytical expression for the dissolution probability $`\mathrm{\Pi }_{n_0}(1)`$ of a jam with an initial width of $`n_0`$ standing cars. To compare the analytical predictions with simulation results we use a damage scenario by initializing a finite jam into an undisturbed system. The reaction of the system to such disturbances is characterized by the sensitivity $`S=1\mathrm{\Pi }_{n_0}(1)`$. The sensitivity is simply the probability that a cluster of $`n_0`$ standing cars causes a wide jam (for open boundary conditions) or leads to the jammed state of the system (for periodic boundary conditions). The following initial conditions are considered to produce an area of free flowing vehicles in the density regime $`\rho _1<\rho <\rho _2`$. Remember that the mean inflow into the induced minijam will always be greater than the mean outflow in this area.

(A) Here we set the randomization parameter $`p`$ of moving cars to zero, i.e. fluctuations in free flow are suppressed. We choose an open system with a sufficiently large megajam at the left boundary in this scenario. The gap distribution is than realized through the outflow of this megajam and fluctuations of free flowing vehicles are completely suppressed. The randomization parameter $`\stackrel{~}{p}_0`$ of standing cars in the megajam is used as control parameter for $`\beta `$. After a car has left the megajam it is therefore reset to $`p_0`$. Thus the distance between two consecutive cars is always a multiple of $`v_{\mathrm{max}}`$ whereby the inflow into the damage is an i.i.d. random variable in that case in correspondence to the theory.
(B) Also in the second case the free-flow car distribution is generated through a megajam, but fluctuations are permitted, i.e., $`p>0`$. Hence the gap distribution contains values deviating from $`nv_{\mathrm{max}}`$ due to velocity fluctuations or braking events.
(C) Finally, a homogeneous initialization of cars in a periodic system is considered. Simulations for this situation are performed until the system has relaxed into its free-flow steady state. As a result of the dynamics this steady state shows large deviations in the gap distribution in comparison to the megajam initializations.
Note that the inflow into the induced minijam is controlled in the cases (A)+(B) through the randomization parameter of standing cars at the left boundary (megajam) while in case (C) the inflow is controlled indirectly through the density of the homogeneous initialization. The inflow $`\beta `$ can be obtained easily in the simulations.

After the free-flowing vehicles are initialized according to the scenarios described above we induce a damage in the system by setting the velocity of a randomly chosen vehicle to zero. The acceleration of this stopped car is suppressed in the following update steps until the damage grows to $`n_0`$ cars. Now, the system is updated according to the update rules without any further outer influences. The damage considered can either grow to a large jam or dissolve.

In Fig. 4 we show the simulated sensitivity $`S=1\mathrm{\Pi }_{n_0}(1)`$ and the mean resolving time for the different initial conditions and compare these results to the analytical predictions of Sec. 4. It is clear that for scenario (A) the analytical results are exact. Even scenario (B) shows an excellent agreement with the analytical curve, but small deviations occur due to the fact that the inflow is no i.i.d. random variable anymore although the mean inflow is still identical to that of scenario (A). Furthermore it should be noted that for the considered dissolution of a small damage the time scale is small, so that local deviations in the gap distribution can play an important role. In the case of the homogeneous initialization (C) we found larger deviations from the analytical curve. The origin of this discrepancy is also the gap distribution which in contrast to the megajam initializations is not generated through a stochastic outflow parameter. Instead it is determined by vehicle interactions due to the model dynamics (simulation runs until relaxation). Hereby repulsive forces between the cars lower the probability of finding large gaps. Therefore the theory overestimates the dissolving probability in the case of a periodic system with homogeneous starting conditions (C). Note, that for large $`p`$ spontaneous jams can appear in the free flow region before a car is able to enter the induced minijam (starting condition (B)) or in the case of homogeneous initialization (C) jams can appear in the system before the steady state is reached. Therefore we do not consider initializations with higher $`p`$ in this work. Our aim was to analyze the dynamics of a single jam. Thus we chose the VDR model which exhibits phase separation but only for $`p_0p`$. Nonetheless we want to point out that the random walk approach for the dynamics of a single jam seems to be generic for various stochastic CA models for traffic flow.

## 6 Summary

We have analyzed the VDR model which exhibits phase separation and metastable states in order to obtain a deeper insight into the jamming behavior in cellular automata models for traffic flow. An analytical approach in terms of random walk theory has been suggested to determine characteristic quantities of wide jams, especially resolving probabilities and lifetimes. The analytical results reveal interesting peculiarities of the model. One finds, e.g., a shift in the convergence of the dissolving probability going from $`n_0=1`$ to greater values of $`n_0`$. The analytical predictions are compared to simulation results. For this purpose a damage scenario is considered by initializing a finite single jam into an undisturbed system. Different boundary conditions are assumed in our investigations. In one case we choose open boundary conditions whereby the free-flow car distribution is generated through the outflow of a sufficiently large megajam at the left boundary. Additionally, a periodic system with homogeneous initialization is considered. In this case the gap distribution is generated through the dynamics of the model (simulation are performed until the system has relaxed). For the megajam initialization we found a very good or even exact agreement between simulation and theory. The homogeneous initialization with periodic boundary condition shows larger deviations from the predicted curve but the overall agreement also rather good. We stress that the random walk approach renders the jamming dynamics of the model. Nagel and Paczuski  also analysed the lifetime of jams in another stochastic CA model for traffic flow, the cruise control limit of the NaSch model, and found good agreement with random walk theory. Therefore we assume that the jamming behavior is generic for a lot of the stochastic CA models for traffic flow.

Furthermore, the fact that one can interpret the emergence of wide jams in terms of probability theory points out a main difference between stochastic CA models for traffic flow and hydrodynamical approaches. To explain the major differences we want to focus on the hydrodynamical model introduced by Kerner and Konhäuser . Although the VDR model invokes most of the characteristic properties (i.e., linear growth of emerging jams, moving transition layer, …) of jams in fluid-dynamical models the formation of large jams due to local perturbations is completely different. In the hydrodynamical model a wide jam is formed from an initially homogeneous state through an external damage which exceeds a critical size. This critical size strongly depends on the density whereby damages below this size dissolve and damages above the critical size lead to wide jams. In contrast, in the VDR model an induced damage leads to a wide jam only with a certain probability. There is no quantity like a critical value for the formation of jams but the probability that a damage causes a jam of course depends on the size of the damage and the density of the initial homogeneous state. Additionally one finds spontaneous jam formation without external influences due to velocity fluctuations in the metastable branch of the model.

To conclude, the results presented here are of practical relevance for various applications of traffic flow using stochastic CA models. Due to its simplicity this class of models has become popular for large scale computer simulation (city or highway networks). Complex networks usually contain many bottlenecks such as crossings, lane reductions, traffic lights or traffic signs. Therefore induced jams often play an important role in realistic traffic scenarios and a proper understanding of the jamming process and dynamics is benefitable.