text:
# Long-lived states in synchronized traffic flow. Empirical prompt and dynamical trap model

## I Anomalous property of the synchronized mode

Although the motion of individual vehicles is controlled by the motivated driver behavior rather than by simple physical laws the car ensembles on highways exhibit phenomena widely met in physical systems. Namely, the existence of various states (self-sustained steady-state modes of traffic flow), their coexistence, the phase transitions, etc. (reviewed in KTGF99 ; KIsr ; KL3 ), actually forms a novel branch of physics. Following the pioneering works by Lighthill & Whitham LW , and Richards R  the state of traffic flow on a highway is specified by the vehicle density $`\rho `$ and the mean velocity $`v`$ or, equivalently, by the traffic flow rate $`q=\rho v`$. In other words, the vehicle density $`\rho `$ and the flow rate $`q`$ are regarded as the complete set of state variables, leading to a correspondence between a point in the $`\rho q`$-plane with a particular state of traffic flow. The main assumption adopted in the classical theories (reviews can be found in Nag96 ; Helbook ; Sh00 ; Hrev ; Hrev2 ) is the existence of a relationship between the mean velocity and vehicle density, $`v=\vartheta (\rho )`$. So, the curve $`q=\rho \vartheta (\rho )`$ in the $`\rho q`$-plane is called fundamental diagram which characterizes the possible states of traffic flow.

According to traffic flow data analyzed recently by Kerner & Rehborn K1 ; K2 ; K3 , Kerner K4 ; K5 , and also Neubert et al. Neub99  there are three distinctive states of multilane traffic flow, the free flow ($`F`$), the synchronized mode ($`S`$), and wide moving jams ($`J`$). This endows the multilane traffic flow with a variety of properties. In particular, it has been demonstrated that the self-formation of moving jams mainly proceeds via the sequence of two phase transitions $`FSJ`$ K4 . Both of them are of the first order, i.e. exhibiting breakdown, hysteresis, and nucleation effects K3 . The present paper focuses its attention on a unique property of the synchronized mode itself, its complexity K2 . In contrast to the free flow, the synchronized mode matches a two-dimensional domain on the $`\rho q`$-plane rather than a curve on it (Fig. 1). It means that a track made up of the empirical points $`\{\rho (t),q(t)\}`$ obtained at sequential time moments unpredictably fills a two-dimensional domain. In this regard the synchronized mode is also referred to as a widely scattered traffic state for which a fundamental diagram in the form of a one-dimensional curve did not exist.

The synchronized mode is characterized by strong correlations in the vehicle motion at different lanes ordering the car arrangement across the road. To describe this state a model has to take into account particular details of the multilane car interaction. Whence it follows that the vehicle density $`\rho `$ and the mean velocity $`v`$ averaged over all the lanes across the highway form an incomplete set of state variables. One of the possible ways to tackle this problem is to ascribe to each lane its own vehicle density and mean velocity. Several models for multilane traffic flow applying to such an approach have been proposed, in particular, gas-kinetic theory GK1 ; GK2 ; GK3 ; GK4 ; GK5 ; GK6  and compressible fluid model FL . For a more detailed survey, including also the cellular automata multilane models, we refer to Nag96 ; Helbook ; Sh00 ; Hrev ; Hrev2 . Here we only note that in all such models the cause of the free flow instability is related to a delay in the driver’s response to variations of the headway distance. This instability gives rise to the traffic flow pinch leading to the formation of a new congested phase of vehicle motion.

In spite of the variety of available models for multilane traffic flow up to now no generally accepted explanation for the widely scattered states of the synchronized mode has been proposed. There are various points of view on the mechanism responsible for this phenomenon. In particular, it is a mixture of different vehicle types (cars and trucks) 17 ; NN1 , a heterogeneity in the headway distance 18  as well as in the highway structure 18a , changes in the behavior of “frustrated” drivers 19 , anticipation effects 20 ; 21 ; 22 ; 23 , non-unique equilibrium solutions of the Prigogine-Herman kinetic equation 24 , the existence of certain plateaus in the dependence of the optimal velocity on the headway distance NN2 ; Peter1 , and the variety of oscillating metastable states with different wavelengths 25 ; 26 . The requirement of multilane models to describe synchronized traffic flow has already been stated in GK3 ; LLK98  using hydrodynamic as well as kinetic equations.

Commenting on these models we would like to say the following. First, most of them apply to a one-dimensional representation of traffic flow, discrete or continuum one. Therefore they do not take into account the strong multilane vehicle interaction which is an essential feature of the synchronized mode. Such a theory can explain the complex behavior of traffic flow on a single-lane road rather than on multilane highways or can describe the congested multilane traffic flow providing that the multilane synchronization is a minor effect. Second, according to the analysis of single-vehicle data in Neub99  the synchronized mode is mainly singled out by strong multilane correlations in the vehicle velocity rather than the formal equality of speed at which cars moving on different lanes. In particular, in the observed free flow regime the averaged velocities at various lanes differ from one another only slightly. Therefore the approach based on the description of multilane traffic flow in terms of individual streams along different lanes with their own mean velocities and vehicle densities seems to be doubtful.

Following the classification proposed by Kerner & Rehborn K2  there are three distinctive kinds of synchronized flow: (i) stationary and homogeneous states where both the average speed and flow rate are approximately constant during a fairly long time interval, (ii) states where only the average vehicle speed is stationary (‘homogeneous-in-speed-states’), and (iii) non-stationary and non-homogeneous states. It should be noted that on macroscopic scales congested traffic exhibits a large variety of different phenomena, causing a more detailed classification L2 ; L3 ; Ha1 ; Ha2  (see also review Hrev ). Recently, Kerner formulated several hypotheses about the synchronized mode KTGF99 ; KIsr ; KL3 ; K2 ; H1H2  to explain the properties of the widely scattered states. In particular, regarding the synchronized traffic of types (i) and (ii) he assumed it to contain a whole multitude (continuum) of possible congested traffic states K2 ; H1H2 , which are stable with respect to infinitesimal perturbations. This leads to continuous spatial-temporal transitions between these states KTGF99 ; KIsr ; KL3 ; H1H2 . However, simulation models and further empirical evidence for these hypotheses are still to be found Hrev .

According to the prevailing notion (see, e.g. Ref. Hrev, ) the synchronized mode is mainly identified by the absence of a direct relationship between the locally averaged vehicle velocity and density provided the case of the “stop-and-go wave” pattern has been excluded from consideration. Therefore, the main part of the theoretical models relates the synchronized mode to the instability of the homogeneous traffic flow caused by the driver delay in the response to variations in the headway distance. These models explain the existence of the widely scattered traffic states taking into account the complex behavior of the developed spatial structure. From our point of view this is a quite natural approach that describes only the synchronized mode of type (iii), called “oscillating congested traffic” Hrev . However, to justify Kerner’s hypotheses concerning types (i) and (ii) in the framework of such an approach a traffic flow model has to admit a large variety of stable nonhomogeneous solutions that either are fixed in the space or move as a whole with a constant speed. Such a model has been proposed in paper 25  based on the car-following approach where the existence of the variety of stable nonhomogeneous solutions differing in wavelength has been found numerically. It is an open question whether this variety is due to the periodic boundary conditions adopted in the model or is an inherent property. At least under random noise, the model demonstrates a certain selection of a single stable solution, so, it seems to describe actually type (iii) of the synchronized mode or a certain mixture of types (ii) and (iii). We would like to underline once more that the cause of such a congested traffic state is related to the delay in the driver response giving rise to the instability of homogeneous traffic flow and the multilane synchronization is a secondary effect only. Besides we note paper 26  that models the widely scattered state region by applying actually the idea used also by Kerner KTGF99 ; KIsr ; KL3 ; H1H2  in justifying his hypotheses. Namely it is assumed that in congested traffic drivers try to maintain a constant time to collision. The latter is possible within a wide interval of headway distances if the velocity difference between cars is sufficiently small. However this effect is a consequence of a certain degeneration of the adopted model. Introducing other terms they can again produce the oscillating behavior of the developed nonhomogeneous solution.

This paper presents a new approach to describe the widely scattered states for the synchronized mode of types (i) and (ii) where the multilane car correlations are fundamental in understanding the origin of congestion. The notion ‘synchronized’ traffic seems to embrace distinguished modes of congested multilane traffic flow. So we just deal with another type of the synchronized mode. We think that the role of the multilane correlations in the “oscillating congested traffic” is either of minor importance or is reduced to prevent the traffic flow going into the jammed “stop-and-go way” state. In the next section we will try to show that there is another type of synchronized mode further called ‘light’ synchronized traffic that meets the aforementioned Kerner’s hypotheses. Namely, in this mode the traffic flow dynamics contains fragments at which the vehicle velocity keeps its value approximately unchanged inside relatively long time intervals. The transitions between these fragments proceed through sharp jumps and the particular values of such a quasi-stationary vehicle velocity vary inside a wide interval. The subsequent sections will present the corresponding model developed from our approach formulated in papers W1 ; W2 . It should be noted that in the presented model these quasi-stationary states are not stationary at all, at least in the standard physical meaning, so we prefer to call them long-lived states of synchronized traffic.

## II Single-vehicle data and two types of the observed synchronized traffic

In this section we analyze single-vehicle data collected by loop detectors at the German highway A1 near Cologne between 6-th and 17-th June 1996. Figure 2 illustrates the analyzed section of the highway and the position of the detectors. A detailed description of this data set and the detection devices has been presented in Neub99 . Here, we only briefly recall their main features. The two sets of detectors D1 and D4 were placed nearby the busy intersection between the highway A1 and the highway A 57 (Köln-Nord), while the detectors D2 and D3 are located close to the on-ramps and exit-ramps of the junction Köln-Lövenich. In between there is another junction (AS Köln-Bocklemünd) but with a rather low usage. The most pronounced congested traffic was observed at the detector D1 where the number of lanes is reduced from three to two for cars passing the intersection towards Köln-Lövenich. Therefore, the data collected by the selected detector were mainly studied in Neub99 . The congested traffic observed near the detector D1 is characterized by a sharp fall in the mean velocity persisting within several hours. In connection with the following we prefer to call it heavy synchronized mode or “oscillating congested traffic”.

The detectors D2, D3, and D4 also recorded congested traffic that, however, seems to be of another nature that will be called light synchronized mode in the following. Exactly in the vicinity of the detectors D2, D3, and D4 traffic flow demonstrated the behavior meeting Kerner’s hypotheses in the sense discussed at the end of the previous section. The present paper does not pretend to a thorough analysis of these data, we only demonstrate the fact that the hypothetical behavior of traffic flow has been recorded directly and is likely not to be a misinterpretation of non-stationary inhomogeneous vehicle structures moving on the highway.

For each detector the data set is mainly composed of the sequence of numbers $`\{t_i,v_i,l_i,k_i\}`$ showing that $`i`$-th car of type $`k_i`$ (passenger car, truck, truck trailer etc.) passed the given detector at time $`t_i`$ moving on lane $`l_i`$ with speed $`v_i`$. It should be pointed out that the passing times $`\{t_i\}`$ were measured with second’s accuracy. Below, the time variations of $`q_l(t,\mathrm{\Theta })`$ and $`v_l(t,\mathrm{\Theta })`$ will be analyzed. Here, $`q_l(t,\mathrm{\Theta })`$ is the traffic flow rate and $`v_l(t,\mathrm{\Theta })`$ the vehicle velocity on lane $`l`$ at time $`t`$ (measured in seconds) averaged over the time interval $`2\mathrm{\Theta }`$ according to:

$`q_l(t,\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \underset{t_i}{}}w_\mathrm{\Theta }(tt_i),`$ (1)
$`v_l(t,\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{1}{q_l(t,\mathrm{\Theta })}}{\displaystyle \underset{t_i}{}}w_\mathrm{\Theta }(tt_i)v_i.`$ (2)
The weight coefficients $`w_\mathrm{\Theta }(t)`$ are given by:

$$w_\mathrm{\Theta }(t)=\{\begin{array}{cc}Z\mathrm{exp}\left(\left|t\right|/\mathrm{\Theta }\right)\hfill & \text{for }\left|t\right|2\mathrm{\Theta }\text{,}\hfill \\ 0\hfill & \text{for }\left|t\right|>2\mathrm{\Theta }\text{,}\hfill \end{array}$$
(3)
subject to the normalization condition

$$\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}w_\mathrm{\Theta }(t)=1,$$
(4)
For $`\mathrm{\Theta }1`$, the factor $`Z`$ is given by $`Z[2\mathrm{\Theta }(1e^2)]^1`$. According to this expression the characteristic time scale of averaging with the weight $`w_\mathrm{\Theta }(t)`$ is $`2\mathrm{\Theta }`$, that is why we will refer to the time interval of averaging the data set keeping in mind the value $`2\mathrm{\Theta }`$ rather than $`\mathrm{\Theta }`$.

As has been mentioned above, heavy congested traffic is a typical case in the vicinity of the detector D1 during the morning hours of workdays and is accompanied with a substantial fall in the mean vehicle velocity. The other detectors measured solely light congested traffic when the mean velocity, on the average, does not drop essentially but exhibits strong fluctuations only. To illustrate these types of traffic flow we have chosen, by way of example, the time series of velocity and flow collected by the detectors D1 and D2 on Th 6 June 1996 and Sa 15 June 1996, respectively. Figure 3 visualizes these data. The row 1 left and row 2 left windows show time variations in the vehicle velocities at different lanes averaged over 1 minute time interval according to Eq. (2). The visualization time intervals (from 9:30 to 12:30 a.m. for D2 and from 9:00 to 12:00 a.m. for D1) span all the time during which the congested flow was recorded on these days. In spite of the clearly visible difference in the time patterns at the detectors D2 and D1 the vehicle velocities at the left and right lanes are strongly correlated. So both the observed modes of congested traffic can be categorized as the synchronized traffic flow. It should be noted that the mean vehicle speed at the off-ramp lane near the detector D1 (Fig. 2) was at least twice as high as the mean speeds at the inner-lanes on 6 June 1996 and did not exhibit strong fluctuations during the given time interval. It seems that the traffic flow on this off-ramp was “free” in comparison with the main stream and did not affect it essentially.

The row 1 right and row 2 right windows demonstrate time variations in the vehicle velocity and the traffic flow rate within two-lane 4-minute averaging. The depicted time interval has been extended to make possible comparison of the congested state with the free flow. We see that the flow rate at the detector D2 did not practically show any visible dependence on the traffic state, whereas in the vicinity of the detector D1 it exhibits strong variations rather than a considerable fall. The latter feature prompts us that the congestion observed near the detector D1 is not the developed “stop-and-go wave” pattern. This statement, however, will be justified more carefully below.

Let us now demonstrate that the two time series recorded at D1 and D2 present actually two synchronized modes different in nature. To substantiate this statement we depicted in the row 3 left and right windows the same time series of the vehicle velocity but on smaller scales (from 10:00 to 11:00 a.m.) and averaged only over 30-second interval. First, we see that the multilane correlations hold also on these scales. Second, the difference between the two time series becomes pronounced. The time series collected by D2 (the row 3 left window) contains several quasi-plateaus in the time variations of the vehicle velocities (singled out in this window and also in the row 4 left one by dotted ellipses). In contrast, the velocity time series obtained by D1 seems to be no more than a collection of spikes varying in amplitude. To argue for the statement about the different nature of these congestions the row 4 left and right windows compare these 30-second averaging time series of the vehicle velocity and flow rate simultaneously for the left lanes. For the detector D2 we meet a visible distinction between the time pattern of vehicle velocity and that of the flow rate. For the detector D1 they are quite similar in appearance.

To make this difference more pronounced we have analyzed the variance for the time series of the vehicle velocity and flow rate collected by the detectors D1 and D2 at the left lane during the observation time from 10:00 to 11:00 a.m. Fig. 4 presents the corresponding relative deviations

$$\frac{\delta v}{v}|_{\text{D1}},\frac{\delta q}{q}|_{\text{D1}},\frac{\delta v}{v}|_{\text{D2}},\frac{\delta q}{q}|_{\text{D2}}$$

depending on the averaging interval $`2\mathrm{\Theta }`$. For the D1 time series we meet at first a weak dependence of the relative deviations on the averaging scale $`2\mathrm{\Theta }`$ for both the vehicle velocity and flow rate which, then, decrease fast keeping approximately their ratio fixed, i.e.

$$\frac{\delta v(\mathrm{\Theta })}{v}|_{\text{D1}}\frac{\delta q(\mathrm{\Theta })}{q}|_{\text{D1}}.$$

However, the question as to whether the latter estimate has a certain physical meaning or just is a simple coincidence requires an individual analysis. In contrast, the D2 time series demonstrates another behavior. The relative deviation of the vehicle velocity exhibits a weak decrease approximately from 22% to 13% as the averaging interval grows 16-fold from 30 seconds to 8 minutes. Conversely, the relative deviation of the flow rate is scaled practically as (straight solid line in Fig. 4)

$$\frac{\delta q(\mathrm{\Theta })}{q}|_{\text{D2}}12\%\times \left(\frac{0.5\text{min}}{2\mathrm{\Theta }}\right)^{{\scriptscriptstyle \frac{1}{2}}}.$$
(5)
Expression (5) enables us to suppose that the fluctuations of the flow rate observed at the detector D2 during one hour from 10 to 11 a.m. are due to a certain white noise. Moreover, the mean flow rate $`q|_{\text{D2, left lane}}`$ during the given time interval was about 1800 veh/h/l, so on average the detector D2 recorded one car passing it on the left lane per two seconds. Keeping in mind that the detectors measured the times of car passing within seconds we may treat each record of a car having passed the detector as a random event with the probability $`p0.5`$. Let us consider these events mutually independent and describe them with a sequence $`\{\gamma (i)\}`$ of random numbers such that take values equal to either 1 or 0 with the probability $`p`$ and $`1p`$, respective, i.e.

$$\gamma (i)=p\text{and}\gamma (i)\gamma (i^{})=p\delta _{ii^{}}+p^2(1\delta _{ii^{}}),$$
(6)
where $`\delta _{ii^{}}`$ is Kronecker’s symbol. In these terms Eq. (1) can be rewritten as

$$q_l(t,\mathrm{\Theta })=\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}w_\mathrm{\Theta }(ti)\gamma (i).$$

Whence we immediately get the following expression for the mean relative deviation of the random variable $`q_l(t,\theta )`$:

$$\frac{\left([q_l(t,\mathrm{\Theta })q_l(t,\mathrm{\Theta })]^2\right)^{{\scriptscriptstyle \frac{1}{2}}}}{q_l(t,\mathrm{\Theta })}=\sqrt{\frac{1p}{2p}\frac{1}{2\mathrm{\Theta }}}.$$
(7)
In obtaining Eq. (7) we have taken into account the normalization condition (4) and for the sake of simplicity extended the exponential dependence (3) of the weight coefficients $`w_\mathrm{\Theta }(t)`$ to all the values of its argument $`t`$. It should be recalled that the value $`\mathrm{\Theta }`$ appearing in (3) matches the time scale measured in seconds. For $`p=0.5`$ and the time scale $`\mathrm{\Theta }`$ measured in minutes Eq. (7) directly gives us expression (5) with the replacement of 12% by 13%.

Thereby, the time variations in the flow rate recorded by the detector D2 are likely to be due to the discreteness of traffic flow bearing no correlations at all. In no case such fluctuations should be ascribed to the state of traffic flow. Thus dealing with the congested traffic in the vicinity of the detector D2 we actually meet the traffic flow state where the vehicle velocity exhibits strong long-time fluctuations with certain time plateaus and the traffic rate, in contrast, has to be regarded as a constant value. These plateaus are separated by sharp jumps as illustrated in Fig. 3, the row 3 left and row 4 left windows with respect the plateaus singled out by the first two dotted ellipses. In the row 3 left windows it is also clearly visible that the plateau position with respect to the $`v`$-axis varies substantially. It should be noted that Neubert et al. Neub99  also have found strong long-time correlations in the vehicle velocity and their absence in the headway distribution but with respect to the D1 time series. They also have obtained the flow rate independence from the vehicle density for the synchronized mode but after averaging over all the D1 time series. Besides, there are cellular automata models predicting also a constant value of the flow rate inside a certain interval of the vehicle density which is due to local highway defects (see Ref. CA1, ; CA2, ; CA3, ; CA4, ; CA5,  and the review by Chowdhury et al. Sh00 ). However, in the case under consideration a constant value should be ascribed to the current flow rate. The latter enables us also to suppose that a possible spatial pattern of the vehicle distribution was fixed on the highway during the observation time. So the D2 time series of the vehicle velocity is likely to visualize the real dynamics of traffic flow caused by its internal local properties rather than by the pinch effects. In this connection it should noted that Lee et al. L2 ; L3  and Helbing et al. Ha1  preferred to place the congested traffic states matching spatial vehicle patterns fixed as a whole on highways into an individual group.

In order to demonstrate that the analyzed traffic states belong in fact to the synchronized mode we present the fundamental diagram (Fig. 5) reflecting the traffic flow states on the $`\rho q`$-plane. For a chosen averaging interval $`2\mathrm{\Theta }`$ the current values of the flow rate $`q(t,\mathrm{\Theta })`$ and the vehicle velocity $`v(t,\mathrm{\Theta })`$ were calculated according to Eqs. (1), (2), and then the effective vehicle density $`\rho `$ was calculated as $`\rho (t,\mathrm{\Theta })=q(t,\mathrm{\Theta })/v(t,\mathrm{\Theta })`$. The row 1 and row 2 windows exhibit the fundamental diagrams for the traffic flow on the left lane near the detectors D2 and D1, respectively, depending on the averaging interval $`2\mathrm{\Theta }`$. Red color has been used to make visible the states of traffic flow recorded within the time from 10 a.m. to 11 a.m. and thus to separate the free flow states from the those of the congested traffic. We see that, in fact, the congested traffic flow under consideration possesses the widely scattered state domain and so can be referred to as the synchronized mode. It is most clear for the time series averaged over 1-minute interval. On one hand, the widely scattered states cover a sufficiently large domain on the $`\rho q`$-plane. On the other hand, for the D2 time series the gap between this domain and the free flow region becomes quite visible. For the D1 time series a certain “tail” of the widely scattered state domain in the vicinity of the origin $`(\rho =0,q=0)`$ (see the row 2 left window) typical for the “stop-and-go way” pattern Neub99  disappears on these scales. So the “stop-and-go way” traffic flow was not, at least, the dominant phase state.

The windows in row 3 show the fundamental diagram for the traffic flow states recorded during the time from 10 to 11 a.m. These windows actually exhibit the structure of the widely scattered state domains visualized in the given time series. Roughly speaking, both of them comprise a large amount of straight lines describing quasi-free flow states characterized by the relationship $`q=\rho v`$ with a constant coefficient $`v`$ and the continuous transitions between them. The particular magnitudes of $`v`$, however, do not coincide with the free flow velocity being about 100 km/h (see Fig. 3, rows 1,2 right windows) but vary inside a wide interval from 35 km/h to 90 km/h for the D2 series and from 5 km/h to 45 km/h for the D1 series. Concerning the lower boundary of the latter interval, however, it should be pointed out that the result can depend on the way of calculating the vehicle density Neub99 . For the D2 time series this feature remains also after averaging over 1-minute interval, whereas for the D1 series it disappears and only the transition lines are visible.

The row 1 and row 2 windows again raise the question about the fundamental diagram which has been discussed in many papers (see, e.g., reviews Sh00 ; Hrev ), in particular, which averaging interval should be used. As can be seen in these windows the scattered state domain shrank considerably on scales about 8 minutes. Besides, keeping in mind the drawn conclusion about the nature of the flow rate fluctuations recorded by the detector D2 we can suppose that the widely scattered states of the light synchronized mode are due to the traffic flow discreteness. The vehicle velocity, however, exhibit strong fluctuations comprising many long-lived states and the continuous jump transitions between them.

## III Synchronized mode as a glass-like phase state

When vehicles move on a multilane highway without changing the lanes they interact practically with the nearest neighbors ahead only. Therefore, there should be no internal correlations in the vehicle flow at different lanes. In particular, the drivers that would prefer to move faster than the statistically mean driver will bunch up forming the platoons headed by a relatively slower vehicle. When the cars begin to change lanes for overtaking slow vehicles the car ensembles at different lanes will affect one another. The cause of this interaction is due to the fact that a car occupies during a lane change manoeuvre two lanes simultaneously, affecting the cars moving behind it at both lanes. Figure 6($`b`$) illustrates this interaction for cars 1 and 2 through car 4 changing the lanes. The drivers of both cars 1 and 2 have to regard car 4 as the nearest neighbor and, so, their motion will be correlated during the given manoeuvres. In the same way car 1 is affected by car 3 because the motion of car 4 directly depends on the behavior of car 3. The more frequently lane changing is performed, the more correlated traffic flow is on a multilane highway. Therefore, to characterize traffic flow on multilane highways it is reasonable to introduce an additional state variable, the order parameter $`h`$ W1 ; W2 . This variable is the mean density of such car triplets normalized to the maximum possible for the given highway and, so, can play the role of a measure of the multilane correlations in the vehicle flow.

On the other hand, the order parameter $`h`$ introduced in this way can be regarded as a measure of the vehicle arrangement regularity. Let us discuss this question in detail for the free flow and synchronized mode individually. In the free flow the feasibility of overtaking makes the vehicle arrangement more regular because of platoon dissipation. So, as the order parameter $`h`$ grows the free traffic becomes more regular. Nevertheless, in this case the density of the car multilane triplets remains relatively low, $`h1`$, and the vehicle ensembles should exhibit weak correlations. Figure 7 illustrates the $`FS`$ transition. As the car density grows in free flow, the “fast” drivers that at first overtake slow vehicles individually begin to gather into platoons headed by more “slow” cars among them but, nevertheless, moving faster than the statistically mean vehicle, see Fig. 7($`a`$). The platoons are formed by drivers preferring to move as fast as possible keeping short headway distance without lane changing. Such a state of the traffic flow should be sufficiently inhomogeneous and is the reason for the occurrence of high-flow states in free traffic. Therefore, even at a sufficiently high car density the free flow should be characterized by weak multilane correlations and not too big values of the order parameter $`h_f`$. The structure of these platoons is also inhomogeneous: they comprise cars whose drivers would prefer to move at different headways (for a fixed velocity) under comfortable conditions. So, when the density of vehicles attains sufficiently high values and their mean velocity decreases remarkably with respect to the velocity on the empty highway some of the “fast” drivers can decide that there is no reason to move so slowly at such short headway requiring strain. Then they can either overtake the car heading the current platoon by changing lanes individually or leave the platoon and take vacant places, compare Fig. 7($`a`$). The former has to increase the multilane correlations and, in part, to decrease the mean vehicle velocity because the other drivers should give place for this manoeuvres in sufficiently dense traffic flow. The latter also will decrease the mean vehicle velocity because these places are vacant from the standpoint of sufficiently “fast” drivers only but not from the point of view of the statistically mean ones preferring to keep longer headway in comparison with the platoon headway. That means, that the statistically mean drivers have to decelerate, decreasing the mean vehicle velocity. The two manoeuvre types make the traffic flow more homogeneous by dissipating the platoons and smoothing the headway distribution, see Fig. 7($`b`$). Besides, the single-vehicle data Neub99  show that the synchronized mode is characterized by long-distant correlations in the vehicle velocities. The headway fluctuations are correlated only on small scales. These findings justify the assumptions of the synchronized mode being a more homogeneous state than the free flow. We think W2  that the given scenario describes the synchronized mode formation which must be characterized by a large value of the order parameter, $`h_s>h_f`$, and a lower velocity in comparison with the free flow at the same vehicle density. It should be noted that this mechanism of the synchronized mode emergence is related to the hypothesis by Kerner KIsr ; KL3 ; NP  that the synchronized mode is caused by a “Z”-like form of the overtaking probability depending on the vehicle density.

Keeping in mind this scenario we have assumed the mean velocity $`v`$ of multilane traffic flow to be determined by both the vehicle density $`\rho `$ and the order parameter $`h`$, namely, $`v=\vartheta (h,\rho )`$, where the function $`\vartheta (h,\rho )`$ is regarded as a phenomenological relationship known beforehand. Its general properties have been discussed in Ref. W1, ; W2, , in particular, it has been assumed to be a decreasing function with respect to $`h`$ for the synchronized mode. Previously we have written the governing equation for the order parameter completely describing the traffic state evolution in terms of time variations in the mean velocity $`v`$ and the order parameter $`h`$, the vehicle density is treated as a fixed value. In this way we have obtained a simple and natural explanation of the observed sequence of phase transitions $`FSJ`$, with each of them being of the first order. In particular, this model predicts that the order parameter $`h`$ exhibits a sharp jump from the value $`h_F(\rho _c)1`$ corresponding to the free flow up to the order parameter $`h_S(\rho _c)1`$ of the synchronized mode when the vehicle density $`\rho `$ exceeds a certain critical value $`\rho _c`$. The further evolution of the synchronized mode as the vehicle density grows is given via the dependence $`v=\vartheta [h_S(\rho ),\rho ]`$.

However, to describe the widely scattered state region of the synchronized mode we have to regard the multilane car interaction into more detail. The matter is that for a car to be able to leave a platoon or to change a lane, the local vehicle arrangement at the neighboring lanes should be of a special form (Fig. 8). So the car rearrangement essentially depends also on the particular details of the neighboring car configuration exhibiting substantial fluctuations. These fluctuations, however, are of another nature than those in the free flow. The latter can be treated as a white noise and slightly blurs the dependence $`v=\vartheta (h,\rho )`$ for the free flow (Fig. 1). In the synchronized mode, by contrast, attaining the optimal conditions of driving, including also overtaking slow vehicles, can be frustrated for a certain time as illustrated in Fig. 8. For car 1 to be able to overtake car 2 the neighboring car 3 should provide a room for this manoeuvre. Otherwise the driver of car 1 has to wait and the local car arrangement will not vary substantially. In other words, changes in the particular realizations of the local car arrangement can be frozen for a certain time although the globally optimal car configuration is not attained at the current time moment. Due to this self-freezing effect the synchronized mode can comprise a great amount of locally metastable states and correspond to a certain two-dimensional region in the $`\rho q`$-plane rather than to a line $`q=\vartheta (\rho )\rho `$. This feature seems to be similar to that met in physical media with local order, for example, in glasses where phase transitions are characterized by wide range of controlling parameters (temperature, pressure, etc.) rather than their fixed values (see, e.g., Ziman ).

The synchronized mode consists of clusters of some ten cars moving along the road as a whole Neub99  which will be called for convenience the fundamental clusters of the synchronized mode. Thereby, the introduced order parameter $`h`$ is actually an averaged characteristics of such a cluster as whole and, so, it does not allow for local fluctuations in the car arrangement inside this cluster. To describe the latter another state variable, the order parameter $`a`$ has to be introduced. Its physical meaning is to allow for the fact that the optimal and comfortable way of driving which the drivers try to attain individually depends essentially on the current realization of the car arrangement inside the given fundamental cluster. In other words, we assume that the order parameter $`h_0(a,\rho )`$ matching this optimal way of driving is a function of the vehicle density $`\rho `$ as well as the order parameter $`a`$. Without lost of generality in such a phenomenological approach we may relate $`h_0(a=0,\rho )`$ to the global optimum of the driving manner attained after averaging over all the realizations of the car configuration, i.e. set $`h_0(a=0,\rho )=h_S(\rho )`$. Besides, it is reasonable to relate variations of the order parameter $`a`$ of amplitude about unity, $`a^2^{1/2}1`$, with all the possible realizations of the car configuration.

To complete the description of the long-lived state continuum of the synchronized mode we should specify the evolution of the order parameter $`a`$. Since the order parameter $`a`$ allows for microscopic details of the fundamental cluster structure its fluctuations will be treated as a random noise whose amplitude depend on the vehicle density only. In contrast, the rate $`da/dt`$ of time variations in the order parameter $`a`$ has to be affected substantially by the current value of the order parameter $`h`$. In fact, as the order parameter $`h`$ tends to the local optimum value $`h_0(a,\rho )`$ for the given $`a`$ the rate $`da/dt`$ should be depressed because all the drivers forming the fundamental cluster prefer to wait until a more comfortable car configuration arises therefore inhibiting the evolution of the fundamental cluster structure.

Now, the model can be formulated. However, before doing this it is worth noting that the present paper does not describes the formation and evolution of the synchronized mode inside the free flow. We construct a continuum model for the long-lived states. So we may confine our consideration solely to the local properties ignoring the interaction between fragments of the synchronized traffic flow separated in space. The latter is worthy of individual and detailed investigations.

## IV The dynamical trap model

Keeping in mind the aforesaid and following W1 ; W2 , we specify the local state of traffic flow on a multilane highway by the system of four phase variables, the mean velocity $`v`$, the density $`\rho `$ of vehicles, and two order parameters $`h`$ and $`a`$. The parameter $`h`$ describes the effect of multilane car interaction on the individual driver behavior, so, the mean velocity $`v=\vartheta (h,\rho )`$ is a direct function of the order parameter $`h`$ and the vehicle density $`\rho `$ given beforehand and meeting general properties described in W1 ; W2 . The other order parameter $`a`$ allows for fluctuations in the car arrangement inside a fundamental cluster of the synchronized mode and affects the relaxation process of the order parameter $`h`$, namely, we write

$$\tau \frac{dh}{dt}=\varphi (h,a,\rho )+ϵ\sqrt{\tau }\xi (t).$$
(8)
Here, $`\tau `$ is the average time drivers require to come to the decision to begin or to stop overtaking manoeuvres, the function $`\varphi (h,a,\rho )`$ describes the regular behavior of individual drivers depending on the traffic state as well as on the particular details of the local car arrangement. Finally, the term $`ϵ\sqrt{\tau }\xi (t)`$, a stochastic Langevin force, allows for a random component in the driver evaluation of the current situation on the road:

$$\xi (t)=0,\xi (t)\xi (t^{})=2\delta (tt^{}).$$
(9)
The factor $`ϵ`$ quantifies the size of the noise. As discussed in the previous section it is reasonable to normalize the order parameters $`h`$ and $`a`$ to their maximum values, which enables us to accept that in the free flow the averaged value of the order parameter $`h_F[\rho ]`$ is much less than unity, $`h_F[\rho ]1`$, whereas in the synchronized mode it is about one, $`h_S[\rho ]1`$, and to assume that the order parameter $`a`$ varies around zero, i.e. $`a^21`$.

In order to proceed, recall briefly the key features of the synchronized mode. When the vehicle density exceeds a certain critical value $`\rho _c`$ the free flow becomes unstable and the synchronized modes arises on the road through a first order transition. Our previous model describes this $`FS`$ transition as a sharp jump in the order parameter $`h`$ from $`h_F[\rho ]`$ to $`h_S[\rho ]`$ at a fixed value of the vehicle density $`\rho `$ which should be accompanied with stepwise jumps in the mean velocity $`\delta v=v_F[\rho ]v_S[\rho ]`$ and the traffic flow rate $`\delta q=`$ $`\rho \delta v`$. These jumps stem from the dependence of the steady state order parameter $`h_{\text{st}}[\rho ]`$ on the vehicle density $`\rho `$ taking the multivalued form for $`\rho >\rho _c`$. In the given case the points $`(\rho ,q_F[\rho ])`$ and $`(\rho ,q_S[\rho ])`$ on the $`\rho q`$-plane have to correspond to well determined states of traffic flow. However, the available empirical data (illustrated in Fig. 1) justify the latter assumption with respect to the free flow only, whereas the synchronized mode stretches over a whole region with, may be, blurred boundaries $`q_S^{(l)}[\rho ]`$ and $`q_S^{(u)}[\rho ]`$. Besides, this region turns out to be sufficiently thick, $`\left(q_S^{(u)}[\rho ]q_S^{(l)}[\rho ]\right)q_F[\rho _c]`$ and adjacent to the free flow branch, $`\left(q_F[\rho ]q_S^{(u)}[\rho ]\right)q_F[\rho _c]`$. We intend to explain the existence of the synchronized mode continuum by the effect of random fluctuations in the vehicle velocity $`v`$. However, it is not sufficient to introduce additional stochastic Langevin forces into governing equations, as done with respect to equation (8). Indeed, in this case such a stochastic force quantifies the random component of the individual driver behavior, so, it should demonstrate similar features in both the free flow and synchronized mode. Therefore, in particular, we consider the amplitude $`ϵ`$ of the introduced stochastic force as well as its correlation function (9) to be independent of traffic flow state. The free flow matches a line on the $`\rho q`$-plane, the free flow branch, which is slightly blurred due to random fluctuations. Besides, it exhibits no long-time correlation. Thereby, an additional Langevin force cannot explain the wide quasi-equilibrium state continuum of the synchronized mode. Moreover, these Langevin forces cannot be strong, so, we may set $`ϵ1`$ in equation (8).

Furthermore, we assume that in the synchronized mode where the multilane car interaction is essential the order parameter $`h`$ is substantially affected by the order parameter $`a`$ quantifying the random fluctuations in the local car arrangement. In this case the possibility of a self-freezing effect discussed above is responsible for the emergence of the long-lived state continuum. Therefore, the next step should specify the function $`\varphi (h,a,\rho )`$ in such a way that the resulting steady-state dependence $`h_{\text{st}}(\rho ,a)`$ on the order parameter $`a`$ is strong. Since we deal solely with the synchronized mode itself we can write down

$$\varphi (h,a,\rho )=hh_0(a,\rho )$$
(10)
without lost of generality. Here the function $`h_0(a,\rho )`$ can be approximated by the expression:

$$h_0(a,\rho )=h_S(\rho )\mathrm{\Delta }a^2$$
(11)
which assumes that a zero value of the order parameter $`a`$ corresponds to the most comfortable conditions for driving, to the highest order of the traffic flow, and to the slowest vehicle motion. We note that in the previous model W1 ; W2  the value $`h_S(\rho )`$ characterizes the stationary state of the synchronized mode at the given vehicle density $`\rho `$. The deviation of the car arrangement from the optimal configuration destroys the order of traffic flow, bringing the synchronized mode closer to the free flow in properties. Therefore, we may set $`\mathrm{\Delta }`$ approximately equal to unity, $`\mathrm{\Delta }1`$, and to regard the traffic state with $`hh_0(1,\rho )`$ as practically the free flow or, more rigorously, the state directly adjacent to the free flow.

Now let us specify the governing equation for the order parameter $`a`$. Local variations of the car arrangement seem to be determined solely by the vehicle density. Indeed, the order parameter $`h`$ is actually produced by averaging, first, over many particular realizations of the fundamental cluster structure similar in properties with respect to the possibility of the lane changing and, second, over all the cars making up the fundamental cluster. So it is sufficiently rough characteristics of the traffic flow state and the dynamics of the order parameter $`a`$ can be treated as an independent of $`h`$ random process of unit amplitude and induced by the white noise except for the rate $`\mathrm{\Omega }(h,a)`$ at which the parameter $`a`$ responds to the white noise. Namely, we write

$$\tau \frac{da}{dt}=\mathrm{\Omega }(h,a)a+\mathrm{\Omega }^{1/2}(h,a)\sqrt{\tau }\varsigma (t),$$
(12)
where again the noise $`\varsigma (t)`$ meets the conditions

$$\varsigma (t)=0,\varsigma (t)\xi (t^{})=0,\varsigma (t)\varsigma (t^{})=2\delta (tt^{})$$
(13)
and the factor $`\mathrm{\Omega }(h,a)`$ describes the self-freezing effect of fluctuations in the car arrangement. In particular, when the order parameter comes close to the local equilibrium, $`hh_0(a,\rho )`$, i.e. $`\varphi (h,a,\rho )0`$ the rate $`\mathrm{\Omega }(h,a)`$ should tend to zero. Therefore, by keeping only the leading term we get

$$\mathrm{\Omega }(h,a)\{\begin{array}{ccc}\varphi ^2(h,a,\rho )/\varphi _0^2& \text{if}& \left|\varphi (h,a,\rho )\right|\varphi _0\\ \mathrm{\Omega }_0& \text{if}& \left|\varphi (h,a,\rho )\right|\varphi _0\end{array}.\text{ }$$
(14)
Here, the constant coefficients $`\varphi _0\mathrm{\Omega }_01`$ because under the general conditions time variations of the behavior of individual drivers as well as of the car configuration of fundamental cluster should proceed essentially at the same rate. The form of the random force in equation (12) has been chosen so that the amplitude of fluctuations in variable $`a`$ remains unchanged provided the factor $`\mathrm{\Omega }(h,a)`$ is constant. In other words, the proximity of the fundamental cluster structure to a local equilibrium slows down the time variations in the order parameter $`a`$. In some sense the curve $`h=h_0(a,\rho )`$ on the $`ah`$-plane (for a fixed value of $`\rho `$) forms certain traps for the path $`\{a(t),h(t)\}`$, which has given the model its name. If the order parameter $`h`$ had been constant, governing equation (12) would describe the collapse of the random motion $`\{a(t)\}`$, with the Ito process corresponding to the maximum of the collapse intensity (see, e.g., Ref. Coll, ). Therefore we also ascribe to the stochastic equation (12) the Ito type, which completes the model formulation. Naturally, in the given model the vehicle density $`\rho `$ is a fixed constant. Time variations in the order parameter $`h(t)`$ can be mapped onto the evolution of the traffic flow rate $`q(t)=\rho \vartheta (h(t),\rho )`$. Therefore, in the following it is sufficient to regard the order parameters $`a`$ and $`h`$ only.

Here we present only the preliminary investigation of the stated model, namely, the results obtained by its numerical simulation, which demonstrate that this model in fact leads the desired long-lived state continuum. For the sake of continuance in numerical simulations we have converted to the dimensionless time $`tt/\tau `$, introduced the variable $`\eta =(h_S(\rho )h)/\mathrm{\Delta }`$, renormalized the constants $`ϵϵ/\mathrm{\Delta }`$ and $`\varphi _0\varphi _0\mathrm{\Delta }`$, set $`\mathrm{\Omega }_0=1`$. This leads to:

$`{\displaystyle \frac{d\eta }{dt}}`$ $`=`$ $`(\eta a^2)+ϵ\xi (t)`$ (15)
$`{\displaystyle \frac{da}{dt}}`$ $`=`$ $`\mathrm{\Omega }_0\varpi (\eta ,a)a+\mathrm{\Omega }_0^{1/2}\varpi ^{1/2}(\eta ,a)\varsigma (t)`$ (16)
where the noise $`\xi (t)`$ and $`\varsigma (t)`$ meets the conditions (9) and (13). Finally, the function $`\varpi (\eta ,a)`$ may be specified as:

$$\varpi (\eta ,a)=\{\begin{array}{ccc}(\eta a^2)^2/\varphi _0^2& \text{if}& \left|\eta a^2\right|\varphi _0\\ 1& \text{if}& \left|\eta a^2\right|>\varphi _0\end{array}.$$

Again, the stochastic equation (16) is assumed to be of the Ito type.

Let us now discuss the results obtained by simulating numerically the system of Eqs. (15), (16) for $`\mathrm{\Omega }_0=1`$, $`ϵ=0.1`$, and $`\varphi _0=0.5`$. First of all, Fig. 9 illustrates the evolution of the order parameters $`\eta `$ and $`a`$ when there is no self-freezing. In this case the order parameter $`a`$ exhibits the standard random pattern of Brownian movement inside a region of unit width. We see a collection of practically independent spikes of unit width. A similar pattern (see Fig. 9($`a`$)) is demonstrated by the dynamics of the order parameter $`\eta `$ provided the interaction of the parameters $`\eta `$ and $`a`$ have been ignored, i.e. by replacing $`(\eta a^2)`$ by $`\eta `$. Naturally, in this case the synchronized mode matches a line on the $`\rho q`$-plane for $`ϵ1`$. The dynamics of the order parameter $`\eta `$ affected by the variable $`a`$ but without the reciprocal influence is shown in Fig. 9($`b`$). Again a collection of spikes can be seen, whose amplitude as well as width has increased tenfold. In other words, if we take into account the effect of the fluctuations in the car arrangement (the order parameter $`a`$) on the individual driver behavior (the order parameter $`h`$) then it is possible to explain an essential blurring of the synchronized mode state.

The dynamics changes dramatically for the full problem, see Fig. 10. The time pattern takes a form corresponding to the long-lived state continuum. When the point $`\{a(t),\eta (t)\}`$ representing the current state of the synchronized mode wanders on the $`a\eta `$-plane and reaches the curve $`\eta =a^2`$ at any point it will be trapped for a certain time until it finally escapes from the trap due to the noise $`ϵ\xi (t)`$. After that the system again wanders in the $`a\eta `$-plane during a time interval about unity before being trapped for the next time. Since the characteristic duration of the trapping is much longer than unity the pattern looks like a certain collection of local metastable states of the synchronized mode. However, such prolonged stays of the system are not metastable in the rigorous meaning and we preferred to call them simply the long-lived states. Since each point of the curve $`\eta =a^2`$ is a trap, the long-lived states make up a certain continuum.

## V Remarks about governing equations for light synchronized traffic

In writing equations (15) and (16) we actually have considered the state of the traffic flow homogeneous. Thereby, the interactions between different fragments separated in space have been neglected. In order to analyze a more realistic case of the evolution of traffic flow with the long-lived states we have to deal with the fields $`\rho (x,t)`$, $`\eta (x,t)`$, and $`a(x,t)`$ describing the vehicle distribution along the highway, the order of the car arrangement as well as local fluctuations in the car arrangement at different highway parts. It should be noted that the given model has been proposed for the light synchronized mode where the driver’s delay to the headway variations seems to be of minor importance. This enables us to assume the local value of the vehicle velocity practically equal to the optimal velocity $`\vartheta (h,\rho )`$ for current values of the vehicle density $`\rho `$ and the order parameter $`h`$.

In this section we only touch a possible way of constructing the distributed model for traffic flow with long-lived states. For inhomogeneous traffic flow Eqs. (15), (16) actually describe the process in the frame locally “attached” to the traffic flow, i.e. moving at the same speed along the road. So, in order to return to the physical coordinate system we, first, should replace the time derivative by the material derivative, i.e.

$$\frac{d}{dt}_t+u(\eta ,\rho )_x,$$

where we again use the dimensionless form of time $`t`$ in the term $`u(\eta ,\rho )=\tau \vartheta [(h_S(\rho )\mathrm{\Delta }\eta ),\rho ]`$.

Second, as was discussed in Sec. III the order parameter $`h`$ characterizes the arrangement of the fundamental car cluster of synchronized traffic flow. So it cannot exhibit considerable variations on spatial scales about the characteristics size $`\mathrm{}`$ of this cluster. So the right-hand side of equation (15) has to be modified to allow for this effect. Following our previous paper W2  we introduce in equation (15) the term

$$\mathrm{}^2_x^2h+2^{1/2}\mathrm{}_xh,$$

smoothing the field $`h`$ on the scales $`\mathrm{}`$ and taking into account the asymmetry in the vehicles interaction. Conversely, the order parameter $`a`$ allows for local variations of the car arrangement inside the fundamental car cluster, so we may regard the field $`a(x,t)`$ as totally independent at different points along highway, at least, within such a mesoscopic description.

In this way the system of equations (15) and (16) can be extended as follows:

$`_t\eta +u(\eta ,\rho )_x\eta `$ $`=`$ $`\mathrm{}^2_x^2h+2^{1/2}\mathrm{}_xh`$ (17)
$`(\eta a^2)+ϵ\xi (t,x)`$
$`_ta+u(\eta ,\rho )_xa`$ $`=`$ $`\mathrm{\Omega }_0\varpi (\eta ,a)a`$ (18)
$`+\mathrm{\Omega }_0^{1/2}\varpi ^{1/2}(\eta ,a)\varsigma (t,x),`$
where the terms $`\xi (t,x)`$ and $`\varsigma (t,x)`$ are random sources uncorrelated in time and space. Naturally, these equations must be completed by the continuity equation:

$$_t\rho +_x\left[u(\eta ,\rho )\rho \right]=(D\tau )_x^2\rho ,$$
(19)
where $`D`$ is the effective diffusion coefficient.

The proposed distributed model for the light synchronized traffic flow is worthy of an individual investigation. Here we have only demonstrated the feasibility of its construction.

## VI Conclusion

We proposed a novel notion of the widely scattered states in the synchronized traffic stimulated by Kerner’s hypotheses K2 ; KTGF99 ; KIsr ; KL3 ; H1H2  about the existence of a multitude of metastable states in the homogeneous or “homogeneous-in-speed” synchronized traffic.

First, we have analyzed the single-vehicle data collected at the German highway A1 near Colonge from 6-th to 17-th of June 1996 in order to verify the presence of the hypothetical multitude of metastable states. As a preliminary result we have found two types of synchronized mode observed at different detectors. One of them can be classified as the “oscillating congested traffic” Hrev  and has been analyzed in detail in Neub99 . The behavior of the other synchronized mode differs essentially therefore we preferred to classify it as an individual type called light synchronized traffic. Exactly this type of traffic flow possesses the hypothetical multitude of quasi-stationary states. Namely, it has been found that the vehicle velocity dynamics consists of a collection of certain fragments inside which the vehicle velocity is approximately constant and the continuous transitions between them occur via sharp jumps. The particular values of the vehicle velocity in these fragments vary in a wide interval. In contrast, the traffic flow rate has to be regarded as a constant because, as we have shown, its fluctuations are mainly due to the traffic flow discreteness. In spite of the latter properties the fundamental diagram drawn for this traffic mode in the standard way contains the widely scattered state domain. Keeping in mind the obtained results we preferred to call this domain the long-lived state region, because it is made up of states that are not stationary at all, at least in the standard physical meaning.

Second, we have proposed a simple mathematical model explaining this complex behavior of light synchronized traffic in terms of dynamical traps. Advancing the idea proposed in W1 ; W2 , we specify the state of traffic flow by four phase variables, the mean velocity $`v`$, the vehicle density $`\rho `$, and two order parameters $`h`$ and $`a`$. The parameter $`h`$ describes the effect of multilane correlations in the vehicle motion on the behavior of individual drivers originating from lane change manoeuvres. On the other hand, it characterizes the order of traffic flow, so the transition from the free flow to the synchronized mode is represented as a sharp jump in the value of the order parameter $`h`$. We assume the existence of the relationship between the mean velocity $`v`$, the vehicle density $`\rho `$ and the order parameter $`h`$, i.e. $`v=\vartheta (h,\rho )`$, where $`\vartheta (h,\rho )`$ is a certain phenomenological function known beforehand. In these terms we can describe the local dynamics of the synchronized mode in terms of the evolution of the order parameter $`h`$, here the vehicle density $`\rho `$ is treated as a fixed constant. The latter is justified because the synchronized mode is characterized by long-distance correlations, so its basic features should stem from the main properties of a certain fundamental cluster, the minimal fragment of car ensemble that moves as a whole and can be regarded as an atomic element of the synchronized mode.

The dynamics of the order parameter $`h`$ is governed by the vehicle density and in addition, by local fluctuations in the car arrangement affecting the feasibility of changing lanes for overtaking. The latter aspect is the key point of the mechanism responsible for the emergence of widely scattered synchronized states. Due to different realizations of the car arrangement inside the fundamental cluster there exist a continuum of locally quasi-equilibrium states of the order parameter $`h`$. So, when such a local quasi-equilibrium state is attained, the individual drivers prefer to wait for a more comfortable car configuration to continue overtaking. In this case also the rate of changes in the car arrangement slows down and, as a result, the current state of traffic flow is ‘frozen’. Exactly this effect of self-freezing gives rise to a continuum of long-lived states in the light synchronized traffic.

We have written a particular mathematical description of self-freezing and justified these qualitative conclusions by solving the stated model numerically.

###### Acknowledgements.

The authors would like to thank the Ministry of Economy and Transport, and the Landschaftsverband Rheinland in the German state Northrhine-Westfalia for the support in obtaining this beautiful data-set. Also this work has been supported in part by Russian Foundation of Basic Research, Grant # 01-01-00389.


format: text

text:
# Tolman-Bondi collapse in scalar-tensor theories as a probe of gravitational memory

## I Introduction

Scalar-tensor (ST) theories of gravity provide a natural alternative to general relativity (GR). They describe gravity with not only a metric $`g_{ab}`$ but also a scalar field $`\varphi `$. Derivatives of $`\varphi `$ appear as source terms in the field equations and $`\varphi `$ itself satisfies a wave equation. The strength of the gravitational coupling is determined by the function $`\omega (\varphi )`$, where GR is recovered in the limits $`\omega \mathrm{}`$ and $`\omega ^3(d\omega /d\varphi )0`$.

ST theories can also be regarded as being equivalent to GR with a varying gravitational “constant” $`G`$. The most simple example of such a ST theory is Brans-Dicke theory BD , where $`\omega (\varphi )`$ is constant and $`G\varphi ^1`$. However, weak field experiments have shown that $`\omega >500`$ Will  and so the deviation from GR is small. For more general ST theories, where $`\omega `$ is not constant, it is possible that $`\omega `$ was much smaller at earlier times. So observations allow such theories to greatly deviate from GR in the early universe.

There has been a renewed interest in ST theories in recent years due to the effective low energy actions of string theory involving one or more scalar fields. These scalar fields enter the field equations in much the same way as the scalar field in ST theories Call . Also, the increasing popularity of inflation and quintessence suggests that scalar fields might need to be incorporated into cosmological models.

The purpose of this paper is to study the effect of an evolving scalar field on the formation and evolution of a primordial black hole. In an asymptotically flat spacetime it is well known that a black hole radiates away any inhomogeneities in the scalar field until it becomes a stationary solution with constant $`\varphi `$ Hawk . This is a consequence of the famous “no hair” theorem. However, in ST cosmological models the scalar field is evolving with time and this would modify how the black hole evolves during its lifetime.

Barrow Barr  was the first to examine this problem. He considered the two extreme possibilities: scenario (A), where the scalar field evolves everywhere homogeneously in the same way as the cosmological background; and scenario (B), where the black hole forces the scalar field to remain constant in some local region around it. The second scenario Barrow called gravitational memory because the black hole would locally preserve the value of the scalar field from when it formed. Barrow & Carr BC  studied the evolution of primordial black holes for these scenarios and found that either case results in a significant deviation from the usual GR analysis. They considered ST theories, where $`G(\varphi )\varphi ^1`$ as in the Brans-Dicke case. Since $`\varphi `$ increases with cosmological time, this implies that $`G(\varphi )`$ decreases, so black holes would take longer to radiate away their mass via Hawking evaporation than in the GR case. In both scenarios the black holes form when gravity is stronger, so the rate of evaporation is less, but in scenario (B) the strength never decreases and so the lifetime is even longer.

The above scenarios are the two extremes and the reality is probably somewhere in between. Two more general scenarios have been proposed GC . In scenario (C) the scalar field evolves faster at the event horizon (EH) than at the particle horizon (PH), so $`\dot{\varphi }_{EH}>\dot{\varphi }_{PH}`$. Eventually the black hole must reach a stage where $`\dot{\varphi }_{EH}=\dot{\varphi }_{PH}`$, but this does not necessarily mean the scalar field is homogeneous since there could be some lag between the asymptotic and local increase. In scenario (D), the scalar field is evolving locally but at a slower rate than asymptotically, so $`\dot{\varphi }_{EH}<\dot{\varphi }_{PH}`$. There is still some gravitational memory but not in the strict sense of scenario (B). In this scenario the gradient of the scalar field is increasing but one would expect there to be some limit, due to the influx of scalar gravitational waves.

Gravitational memory for black holes which are small compared to the cosmological scale (i.e. the particle horizon) has already been investigated by Jacobson Jaco . In this case, the scalar field evolution can be considered as an asymptotic perturbation to the Schwarzschild metric and the end-state of scenario (C) applies, with the lag being found to be small. This suggests that gravitational memory is virtually non-existent. However, this approximation may not apply for primordial black holes since these can have a size comparable to the particle horizon at formation CH . It is still not clear what would happen in this case. Therefore another way of investigating gravitational memory, without assuming that the black hole is small, is needed.

In GR there have been several attempts to study analytic solutions which represent black holes within a cosmological background. The earliest used the Einstein-Straus solution ES , which matches a Schwarzschild interior to a Friedmann exterior, and this approach has also been used to study gravitational memory Sak . However, in most circumstances it can be shown Whin  that such a matching is only possible if the scalar field is constant, which just gives the GR solution. Another method used the McVittie metric McV , but it has been shown Nola  that this has a scalar curvature singularity at the event horizon.

A more successful method uses the Tolman-Bondi metric TB  to represent the collapse of dust to a black hole in an asymptotically Friedmann background. However, this only works for dust and cannot in general be applied to ST theories due to the derivatives of the scalar field appearing as source terms in the field equations. In this paper we overcome this problem by assuming that the effect of the scalar field on the spacetime is small compared to that of the matter. This means that we can use the usual Einstein field equations to generate the spacetime and then use the wave equation for the scalar field to determine its evolution. This approximation was used by Harada et al. HCNN  to calculate the scalar gravitational radiation emitted by Oppenheimer-Snyder collapse in ST theory. Jacobson also uses this approximation when calculating the effect of an evolving scalar field in Schwarzschild spacetime. We note, however, that self-consistent numerical calculations of spherical gravitational collapse in ST theory in asymptotic flat spacetime, in which the effect of the scalar field on the spacetime is fully incorporated, have been considered by previous authors SNN ; SST .

If one makes this approximation to investigate gravitational collapse in a Tolman-Bondi spacetime, the solution is specified by two arbitrary functions: the energy and mass functions. To represent collapse in a Friedmann background, the energy function has to be negative within some radius $`r_0`$ and zero outside it. This results in the eventual gravitational collapse of all the matter within $`r_0`$, while the matter outside $`r_0`$ expands forever as in a flat Friedmann model. In choosing the mass function, or equivalently the density perturbation, we adopt the “compensated” model. This means that the overdense region in the centre is surrounded by an underdense region outside, so that the total mass at infinity is unaffected.

We solve the field equations numerically using the characteristic method. This method was first applied to an inhomogeneous and dynamical background spacetime by Iguchi et al. IHN . The characteristic method integrates over null hypersurfaces with the event horizon as a boundary. This means that one never needs to calculate anything inside the black hole, thereby avoiding any numerical problems associated with singularities. The output of the code shows the spatial and temporal variation of the scalar field. The figures produced show that the initial collapse results in a large gradient in the scalar field. However, as time increases, the scalar field becomes almost homogeneous. This suggests that, within the approximation used, gravitational memory is not possible. It remains to be seen whether the back reaction of the scalar field could alter this conclusion.

In Section II we describe ST theories in more detail and derive the field equations for the approximation in which the effect of the scalar field on the metric is neglected. In Section III we transform the Tolman-Bondi solution into null coordinates, giving the equations necessary to apply the characteristic method. In Section IV we specify the model giving rise to black hole formation. We present the numerical results in Section V and discuss their implications in Section VI.

## II Basic equations

For scalar-tensor theories of gravity with $`G(\varphi )\varphi ^1`$ the field equations are

$$G_{ab}=\frac{8\pi }{\varphi }T_{ab}+\frac{\omega (\varphi )}{\varphi ^2}\left(_a\varphi _a\varphi \frac{1}{2}g_{ab}\varphi ^c\varphi _c\right)+\frac{1}{\varphi }(_a_b\varphi g_{ab}^c_c\varphi ),$$
(1)
$$^c_c\varphi =\frac{8\pi T(d\omega /d\varphi )^c\varphi _c\varphi }{3+2\omega (\varphi )}$$
(2)
where $`T_{ab}`$ is the usual energy-momentum tensor, $`T`$ is its trace and we have set $`c=1`$. In Brans-Dicke theory Eq.(1) remains unchanged but $`d\omega /d\varphi =0`$ in Eq.(2). These equations are expressed in the Jordan frame but ST theories can also be expressed in a conformal frame known as the Einstein frame. The Einstein frame is related to the Jordan frame by the transformation

$$\overline{g}_{ab}=(G_0\varphi )g_{ab}\overline{T}_{ab}=(G_0\varphi )^1T_{ab},$$
(3)
where $`G_0`$ is the present value of the gravitational “constant” as measured in solar system experiments. It is called the Einstein frame because it can be expressed as GR with a scalar field. However, the Jordan frame will be used throughout this paper.

In GR the Tolman-Bondi solution is given by

$$ds^2=dt^2+A^2(t,r)dr^2+R^2(t,r)[d\theta ^2+\mathrm{sin}^2\theta d\psi ^2],$$
(4)
where $`r`$ is the comoving radial coordinate and $`R`$ is given by

$$tt_s(r)=\sqrt{\frac{R^3}{F}}G\left(\frac{fR}{F}\right).$$
(5)
Here $`G(y)`$ is a positive function given by

$$G(y)=\{\begin{array}{cc}\frac{\text{Arcsin}\sqrt{y}}{y^{3/2}}\frac{\sqrt{1y}}{y}\text{or}\frac{\pi \text{Arcsin}\sqrt{y}}{y^{3/2}}+\frac{\sqrt{1y}}{y}\hfill & (0<y1)\hfill \\ \frac{2}{3}\hfill & (y=0)\hfill \\ \frac{\text{Arcsinh}\sqrt{y}}{(y)^{3/2}}\frac{\sqrt{1y}}{y}\hfill & (y<0)\hfill \end{array},$$
(6)
$`t_s(r)`$ is a constant of integration, and $`F(r)=2G_0m(r)`$ with $`m(r)`$ being the mass within radius $`r`$. $`A`$ is given by

$$A^2(t,r)=\frac{R^{}(t,r)^2}{1+f(r)},$$
(7)
and $`R`$ satisfies

$$\dot{R}^2(t,r)=\frac{F(r)}{R}+f(r).$$
(8)
The density of the dust $`\rho `$ is given by

$$\rho =\frac{F^{}}{8\pi R^2R^{}}.$$
(9)
In the above, a dot denotes $`_t`$ and a prime denotes $`_r`$.

There are two arbitrary functions in this solution: the mass function $`m(r)`$ and the energy function $`f(r)`$. Investigating collapse to a black hole just requires the appropriate choice for these functions. In this paper the energy function is chosen such that

$$f(r)<0\mathrm{for}r<r_0$$
(10)
$$f(r)=0\mathrm{for}r>r_0$$
(11)
for some $`r_0`$. This means that when a perturbation is applied to the background dust, all the matter interior to $`r_0`$ will eventually collapse to form a black hole, while the exterior matter will expand forever as in a flat Friedman universe DL . The mass function $`m(r)`$ is determined by putting $`R=r`$ at $`t=t_s(r)`$ in Eq. (5).

The key to this approximation is that the back reaction of the scalar field is neglected. It is assumed that the effect of the scalar field on the spacetime is small compared to that of the matter. The initial configuration used is the general relativistic one with constant $`\varphi `$ and for simplicity Brans-Dicke theory is used. Then, to the lowest order, the evolution of $`\varphi `$ is determined by the wave equation:

$$\left[^c_c\right]_{TB}\varphi =\frac{8\pi }{3+2\omega }T_{TB}$$
(12)
where $`\left[^c_c\right]_{TB}`$ and $`T_{TB}`$ are determined for the Tolman-Bondi metric and the general relativistic solution. Using the Tolman-Bondi metric, the wave operator is given by

$$\left[^c_c\right]_{TB}\varphi =\ddot{\varphi }\frac{(AR^2)^{}}{AR^2}\dot{\varphi }+\frac{1}{A^2}\varphi ^{\prime \prime }+\frac{1}{AR^2}\left(\frac{R^2}{A}\right)^{}\varphi ^{}.$$
(13)
## III Characteristic method

The last equation needs to be rewritten in terms of a null coordinate suitable for the characteristic method. The retarded time coordinate $`u`$ is introduced such that $`u=\mathrm{constant}`$ is an outgoing null geodesic. In the original coordinates the outgoing null geodesic is given by

$$\frac{dt}{dr}=A,$$
(14)
so we can write

$$\frac{u^{}}{\dot{u}}=A.$$
(15)
The coordinate system is now transformed from $`(t,r)`$ to $`(u(t,r),\overline{r}(r))`$ using the relations

$$du=\frac{1}{\alpha }(dtAdr),d\overline{r}=dr,$$
(16)
where

$$\alpha \frac{1}{\dot{u}}.$$
(17)
The partial derivatives are then related by

$$_t=\frac{1}{\alpha }_u,_r=\frac{A}{\alpha }_u+_{\overline{r}}.$$
(18)
In this coordinate system the metric becomes

$$ds^2=\alpha ^2du^22\alpha A(u,\overline{r})dud\overline{r}+R^2(u,\overline{r})[d\theta ^2+\mathrm{sin}^2\theta d\psi ^2].$$
(19)
To use the characteristic method it is necessary to introduce the derivative along the ingoing radial null geodesic. In the original coordinates the ingoing null geodesic is given by

$$\frac{dt}{dr}=A,$$
(20)
which in the new coordinate system becomes

$$\frac{d\overline{r}}{du}=\frac{\alpha }{2A}.$$
(21)
Therefore the derivative along the ingoing null can be obtained:

$$\frac{d}{du}=_u+\frac{d\overline{r}}{du}_{\overline{r}}=\frac{\alpha }{2}\left(_t\frac{1}{A}_r\right).$$
(22)
The partial derivatives $`_t`$ and $`_r`$ can now be rewritten in terms of $`d/du`$ and $`_{\overline{r}}`$:

$`_t`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{d}{du}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{A}}_{\overline{r}},`$ (23)
$`_r`$ $`=`$ $`{\displaystyle \frac{A}{\alpha }}{\displaystyle \frac{d}{du}}+{\displaystyle \frac{1}{2}}_{\overline{r}}.`$ (24)
The wave operator then becomes

$$\left[^c_c\right]\varphi =\frac{2}{\alpha AR}\frac{d\phi }{du}\frac{A^{}}{A^3R}\phi +\frac{1}{AR}\left[(A\dot{R})^{}\left(\frac{R^{}}{A}\right)^{}\right]\varphi $$
(25)
where

$$\phi _{\overline{r}}(R\varphi ).$$
(26)
Here  and  refer to the operators given in Eq.(23) and Eq.(24) respectively. It is also necessary to obtain an equation for $`\alpha `$. This is achieved by using $`(\dot{u})^{}=(u^{})^{}`$, which gives

$$_{\overline{r}}\alpha =\dot{A}\alpha .$$
(27)
Applying the full Tolman-Bondi solution, the basic equations that must be solved in Brans-Dicke theory are to first order

$`{\displaystyle \frac{d\phi }{du}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{2}}{\displaystyle \frac{\sqrt{1+f}}{R^{}}}\left({\displaystyle \frac{f^{}}{2(1+f)}}{\displaystyle \frac{R^{\prime \prime }}{R^{}}}\right)\phi +{\displaystyle \frac{\alpha }{2}}{\displaystyle \frac{F}{R\sqrt{1+f}}}\left({\displaystyle \frac{F^{}}{2F}}{\displaystyle \frac{R^{}}{R}}\right)\varphi `$ (28)
$`+{\displaystyle \frac{\alpha }{2}}{\displaystyle \frac{1}{3+2\omega }}{\displaystyle \frac{F^{}}{R\sqrt{1+f}}},`$
$`_{\overline{r}}\alpha `$ $`=`$ $`\pm {\displaystyle \frac{\alpha }{2}}{\displaystyle \frac{1}{\sqrt{(1+f)\left(\frac{F}{R}+f\right)}}}\left({\displaystyle \frac{F^{}}{R}}{\displaystyle \frac{FR^{}}{R^2}}+f^{}\right),`$ (29)
where the upper and lower signs correspond to an expanding and collapsing phase respectively. For numerical purposes it is also convenient to take the parametrized form of the Tolman-Bondi solution. For $`f=0`$, $`R`$ is given by

$$R=\left(\frac{9F}{4}\right)^{1/3}(tt_s(r))^{2/3}.$$
(30)
For $`f>0`$, it is given by

$`R`$ $`=`$ $`{\displaystyle \frac{F}{2f}}(\mathrm{cosh}\eta 1),`$ (31)
$`tt_s(r)`$ $`=`$ $`{\displaystyle \frac{F}{2f^{3/2}}}(\mathrm{sinh}\eta \eta ).`$ (32)
For $`f<0`$, it is given by

$`R`$ $`=`$ $`{\displaystyle \frac{F}{2(f)}}(1\mathrm{cos}\eta ),`$ (33)
$`tt_s(r)`$ $`=`$ $`{\displaystyle \frac{F}{2(f)^{3/2}}}(\eta \mathrm{sin}\eta ).`$ (34)
Here the signs are chosen so that they correspond to the big bang universe.

## IV Models

We choose the background primordial black hole model so that the following conditions are satisfied. (1) The big bang occurs at the same time everywhere, i.e. $`t_s(r)=0`$ (constant). (2) The model is asymptotically flat Friedmann and compensated (i.e. the overdense region in the centre is surrounded by an underdense region outside in such a way that the total mass at infinity is unaffected). (3) The model is free of shell-focusing or shell-crossing naked singularities, at least within the calculated region. (4) The central region is bound, while the asymptotic region is marginally bound. (5) At the initial time $`t=t_0`$, the condition $`R^{}>0`$ is satisfied everywhere.

In order to satisfy the above conditions, we set $`t=t_0`$ and choose the energy function $`f(r)`$ to have the form

$$f(r)=\{\begin{array}{cc}\left(\frac{r}{r_\mathrm{c}}\right)^2& \text{for}r<r_\mathrm{w}\hfill \\ \left(\frac{r}{r_\mathrm{c}}\right)^2\mathrm{exp}\left(\left(\frac{rr_\mathrm{w}}{r_\mathrm{w}}\right)^4\right)& \text{for}rr_\mathrm{w}\hfill \end{array}$$
(35)
where $`r_\mathrm{c}`$ gives the curvature radius in the central closed Friedmann region and $`r_\mathrm{w}`$ gives the scale of the overdense region. Equation (35) means that the central region $`r<r_\mathrm{w}`$ is described by the exact closed Friedmann solution, which ensures that there is no shell-focusing naked singularity. More general situations in which shell-focusing is avoided have been discussed elsewhere Josh . The code continually monitors for shell-crossing to check that this never happens within the calculated region. We then determine $`F(r)`$ so that $`r`$ coincides with $`R`$ at the $`t=t_0`$ spacelike hypersurface.

The choice of the function $`f(r)`$ requires some justification since it is nowhere exactly zero in the calculated region, so strictly speaking the entire region would eventually collapse to a black hole if one waited long enough. In practice, this does not matter because our calculated region is so large that the value of $`f`$ is effectively identical to $`0`$ for $`r5r_\mathrm{w}`$ (i.e. at the outer boundary). Since the region is finite, we can therefore always make the matching to the Einstein-de Sitter universe outside the calculated region. It would be easy in principle to perform the calculations for a situation in which the matching occurs within the calculated region (i.e. with $`f=0`$ exactly at the edge of the region). It is clear that this would make no qualitative difference to our conclusions but, to confirm this, we also adopt the choice

$$f(r)=\{\begin{array}{cc}\left(\frac{r}{r_\mathrm{c}}\right)^2\left[1\left(\frac{r}{r_w}\right)^4\right]^4& \text{for}r<r_\mathrm{w}\hfill \\ 0& \text{for}rr_\mathrm{w}\hfill \end{array},$$
(36)
in which the compensation is satisfied explicitly within the calculated region.

Before integrating, we have to fix the initial data for $`\varphi `$. One can take this to be the homogeneous cosmological solution given by

$$\varphi _\mathrm{c}=\varphi _0\left(1+\frac{1}{3+2\omega }\frac{4}{3}\mathrm{ln}\frac{t}{t_0}\right).$$
(37)
We then set the initial null hypersurface as the null cone whose vertex is at $`(t,r)=(t_0,0)`$ and regard the cosmological solution as the initial data on this hypersurface. Although the value of the scalar field at the cosmological particle horizon must be given by this solution, the value in the perturbed region and the surrounding region may be different from this. To examine the sensitivity of the results to this alteration, we consider another form of the initial data which is different from the cosmological solution in the region $`t(110)\times t_0`$. We now choose

$$\varphi _\pm =\varphi _\mathrm{c}\left[1\pm \mathrm{exp}\left[\left(\frac{t}{5t_0}\right)^2\right]\right],$$
(38)
so that we have an ingoing wave in the perturbed region, and examine the evolution of the scalar field thereafter. The numerical code has been checked by the following non-trivial test calculation. In the flat Friedmann universe, the code must reproduce the cosmological evolution Eq. (37) from the initial data. There is agreement to within 0.05% accuracy.

## V Results

We denote the Hubble parameter in the Friedmann background (far from the perturbed region) as $`H_0`$ at $`t=t_0`$. Recall that $`r_\mathrm{c}`$ gives the amplitude of the density perturbation, while $`r_\mathrm{w}`$ gives the size of the perturbed region. For super-horizon scale perturbations, we cannot set the density perturbation to be very large else the overdense region closes up on itself and becomes disconnected from the rest of the universe CH . Actually, the requirement $`R^{}>0`$ imposes an even stronger condition since it implies $`r_\mathrm{w}r_\mathrm{c}`$. We have set the Brans-Dicke parameter to be $`\omega =5`$. If $`r_\mathrm{c}`$ is much increased, then the amplitude of the overdensity is much decreased and the resulting black hole becomes very small compared with the horizon scale at the formation time. If $`r_\mathrm{c}`$ is much decreased, the overdense region becomes a separate closed universe. Models and parameters are summarized in Table 1. The difference between models A, B, C and F is in the choice of the background perturbation. The difference between models A, D and E is in the choice of the initial condition for the scalar field. A change of $`\omega `$ only scales the variation of the scalar field from $`\varphi _0`$, as indicated by Eq. (37).

In Fig. 7, the energy function $`f(r)`$ and the initial density perturbation

$$\delta (t_0,r)\frac{\rho (t_0,r)\rho (t_0,\mathrm{})}{\rho (t_0,\mathrm{})}$$
(39)
are plotted. The trajectories of outgoing null geodesics are plotted in Fig. 7. It is seen that a nearly horizon-scale black hole is formed for model A, while the black hole is smaller than the horizon scale for models B, C and F. This is because the size of the black hole is always roughly $`20H_0^1`$, whereas the time at which it forms is about $`20H_0^1`$ in model A and closer to $`100H_0^1`$ in the other cases. The initial data for the scalar field are set on the initial null Cauchy surface (see Fig. 7). We prepare three sets of initial data, $`\varphi _\mathrm{c}`$, $`\varphi _+`$ and $`\varphi _{}`$, and these are plotted in Fig. 7. We have investigated all the models listed in Table 1 and the results are seen in Figs. 7-7. In Fig. 7, the profile of the scalar field is plotted for constant $`t`$. In Fig. 7, it is plotted for constant $`R`$. The reason why some of the curves come to an abrupt end below some value of $`R`$ is that the event horizon has formed near the center and we did not calculate the evolution of the scalar field inside the event horizon. The reason why some of them come to an abrupt end above some value of $`R`$ is due to the finiteness of the region in which the numerical calculation is done.

We note that collapse ensures that the scalar field is initially concentrated in the central regions and this means that it rises above the asymptotic cosmological value everywhere. However, this central concentration tends to fall due to the underdensity surrounding the black hole. This effect, coupled with the increase of the cosmological value, means that the scalar field necessarily falls below the cosmological value at sufficiently large values of $`R`$, at least for the models under consideration. Eventually it may do so at every plotted value of $`R`$. Strictly speaking, the issue of gravitational memory is concerned with the process whereby the scalar field is raised to the cosmological value once it has fallen below it rather than with the process whereby it initially falls to the cosmological value.

Generally, the configuration of the scalar field tends to spatially homogeneity with increasing time and the value of the scalar field around the black hole follows the cosmological evolution of the scalar field at each moment. To see this more clearly, we identify the value of the scalar field on the latest null ray within the calculation as the value on the black hole event horizon, $`\varphi _{\mathrm{BH}}`$, because the null ray is very close to the event horizon for each model (at least for $`H_0t300`$, see Fig. 7). Both $`\varphi _\mathrm{c}`$ and $`\varphi _{\mathrm{BH}}`$ are plotted in Fig. 7. The end of the curve $`\varphi _{\mathrm{BH}}`$ corresponds to the formation time of the event horizon. We can see that $`\varphi _{\mathrm{BH}}`$ follows the evolution of $`\varphi _\mathrm{c}`$. There is a small deviation of $`\varphi _{\mathrm{BH}}`$ from $`\varphi _\mathrm{c}`$ but this can be explained by the central overdensity and the surrounding underdensity. For $`H_0t100`$, the event horizon runs through the overdense region and hence the scalar field is amplified compared with $`\varphi _\mathrm{c}`$. Thereafter the event horizon runs through the underdense region and so the scalar field is reduced compared with $`\varphi _\mathrm{c}`$. The overall evolution of $`\varphi _{\mathrm{BH}}`$ is well described by $`\varphi _\mathrm{c}`$.

Although there are minor differences among the models, we can conclude that the configuration of the scalar field is nearly spatially homogeneous and well described by the cosmological solution $`\varphi _\mathrm{c}`$, at least for around ten initial Hubble times after the formation of the event horizon.

## VI Summary

We have calculated the evolution of the Brans-Dicke scalar field in the presence of a primordial black hole formed in a flat Friedmann background. We have found that the value of the scalar field at the event horizon almost always maintains the cosmological value. This suggests that primordial black holes “forget” the value of the gravitational constant at their formation epoch. In this sense, we confirm the result of Jacobsen Jaco , although he never carried out an explicit evolutionary calculation. While it is possible that some choices of the function $`f(r)`$ might lead to a different conclusion, we would claim that our conclusion must at least hold in compensated models. For we have tried many sets of parameters, corresponding to models which describe a wide variety of physical situations, and this always seems to be the case. However, it is not clear from our analysis what would happen for non-compensated models and we have also avoided models with shell-crossing singularities, so we do not claim that our result is completely general.

It should also be stressed that this result has only been demonstrated for a dust universe in which the scalar field does not appreciably affect the background curvature and it remains to be seen whether the same conclusion applies when this assumption is dropped. As can be seen from Figs. 7 and 7, both the radial gradient and the time derivative of $`\varphi `$ are large at early times and both of these act as source terms in the field equations. However, as long as $`\omega `$ is large, this is unlikely to stop the scalar field becoming homogeneous eventually, since it is clear that both the radial gradient and time derivative become small well after the initial collapse. Neglecting the back reaction should therefore be a reasonable assumption in this situation, so gravitational memory seems unlikely. However, this conclusion might not apply for $`\omega 1`$.

## Acknowledgments

We would like to thank T. Nakamura for helpful discussion. TH is supported by the Grant-in-Aid (No. 05540) from the Japanese Ministry of Education, Culture, Sports, Science and Technology. CG is supported by the UK Particle Physics & Astronomy Research Council.