Self-organized criticality of traffic flow: There is nothing sweet about the sweet spot

This paper shows that the kinematic wave model exhibits self-organized criticality when initialized with random initial conditions around the critical density. A direct consequence is that conventional traffic management strategies seeking to maximize the flow may be detrimental as they make the system more unpredictable and more prone to collapse. Other implications for traffic flow in the capacity state are discussed, such as: (i) jam sizes obey a power-law distribution with exponents 1/2, implying that both its mean and variance diverge to infinity, and therefore traditional statistical methods fail for prediction and control, (ii) the tendency to be at the critical state is an intrinsic property of traffic flow driven by our desire to travel at the maximum possible speed, (iii) traffic flow in the critical region is chaotic in that it is highly sensitive to initial conditions, (iv) aggregate measures of performance are proportional to the area under a Brownian excursion, and therefore are given by different scalings of the Airy distribution, (v) traffic in the time-space diagram forms self-affine fractals where the basic unit is a triangle, in the shape of the fundamental diagram, containing 3 traffic states: voids, capacity and jams. This fractal nature of traffic flow calls for analysis methods currently not used in our field.


Introduction
There is mounting empirical evidence suggesting the hypotheses that urban networks exhibit self-organized criticality (SOC) , Zhang et al., 2019, Zeng et al., 2020, Bak et al., 1987, a celebrated paradigm in the 90's that has proved very useful in countless areas of science, from earthquakes to stock prices to forest fires to astrophysics (Dȃnilȃ et al., 2015, Aschwanden, 2014. Despite being around for 2 or 3 decades now, the theory and control consequences of SOC have not permeated the transportation literature. This is unfortunate because if it turns out that urban networks do exhibit SOC, a whole new paradigm emerges for understanding urban congestion, based on the lessons learned from complexity, catastrophe theory and non-equilibrium critical phenomena (Taleb, 2020). It would mean for example that most of our current control methods based on standard calculus and statistics techniques may no longer apply or may need a major overhaul. Notably, this paradigm is in stark contrast to the currently accepted view in our field that an optimal strategy is to meter the network so that it remains at the "sweet spot": near the critical density where flows are maximized (Geroliminis et al., 2012, Hajiahmadi et al., 2013, Haddad et al., 2013, Ramezani et al., 2015, Keyvan-Ekbatani et al., 2015, Guo and Ye, 2016, Ampountolas et al., 2017, Haddad, 2017, Ding et al., 2017, Fu et al., 2017, Kouvelas et al., 2017, Sirmatel and Geroliminis, 2017, Yang et al., 2017, Yildirimoglu et al., 2018, Zhong et al., 2018, Guo and Ban, 2020. But given the absence of a general theory of SOC to date (Aschwanden, 2014), it is not yet possible to formally prove or disprove this hypothesis. This paper is a step towards such a theory for traffic flow. The defining characteristic of a SOC system is that it is naturally driven to be in a "critical state" where it undergoes a phase transition (Halperin and Hohenberg, 1969) (e.g. from liquid to solid or from free-flow to congestion). At the critical point the characteristic length scale of the system, such as the size S of forest fires or the size of traffic jams as shown here, diverges to infinity, intimately related to dimensional analysis, commonly used in our field, where the celebrated Buckingham-Pi theorem (Buckingham, 1914) states that a physical system may be fully described in terms of dimensionless combinations of variables, and that typically the describing function is independent of the microscopic details of the problem.
As seen in the introduction, recently de Gier et al. (2019) showed that the NaSch CA model (Nagel and Schreckenberg, 1992) and multi-lane variants that include lane changes belong to the Kardar-Parisi-Zhang (KPZ) universality class of parameter 3/2. This parameter represents the exponent that relates the scaling of space and time variables as t ∼ x 3/2 . For example, the relaxation time of a road segment is proportional to (road length) 3/2 , and the correlation length is proportional to t 2/3 . This result is more general than it looks, because the NaSch model is closely related to many other traffic flow models. In particular, it is known to be a stochastic version of Newell's car following (Newell, 2002) in discrete space, which in turn corresponds to the the kinematic wave model of (Lighthill andWhitham, 1955, Richards, 1956) with triangular fundamental diagram. The fundamental diagram can be readily measured in the field (Newell, 2002); e.g., the critical density free-flow speed and the wave speed. The parameter v max in the NaSch model actually represents the ratio of free-flow speed to the wave speed. When v max = 1 the deterministic NaSch model is effectively equivalent to ECA-184, which also corresponds to the kinematic wave model with isosceles fundamental diagram.
Notice that the proof of universality in de Gier et al. (2019) was shown analytically for v max = 1 but only numerically for v max > 1. One can strengthen this latter result by recalling that Laval and Chilukuri (2016) showed that the slopes of the fundamental diagram are irrelevant. In fact, exploiting a skewing symmetry in the kinematic wave model they showed that the numerical values for these parameters are irrelevant when it comes to the calculation of normalized flows and/or delays. This means that, effectively, traffic flow on a single link requires no parameters (after proper normalization) and therefore is independent of the details such as the value of v max .
Given all of the above, the conjecture becomes that possibly all traffic flow models belong to the same universality class. After all, they all obey the conservation law and have some sort of nonlinear fundamental diagram. Second order effects, bounded accelerations, lane changes and other driving considerations become details under this macroscopic view of traffic models.

Criticality of traffic flow
In this section we show that at the critical density, traffic on a long single-lane road is fully characterized by a standard a Brownian motion in continuous time-space, or a random walk in discreet time-space. Based on our discussion in the previous section, in the remainder of the paper we use the simplest traffic flow model, i.e. the kinematic wave model with isosceles fundamental diagram with the following parameters: free-flow speed = wave speed = 1, jam density = 1, so that it is equivalent to ECA-184. These parameters imply that both the capacity and the critical density are 1/2, that distances are measured in units of one vehicle length and that time is measured in units of 1/2 of a critical headway.

Connection with Brownian motions
Here we exploit the fact that the kinematic wave model can be posed as a Hamilton-Jacoby PDE (Daganzo, 2005, Laval andLeclercq, 2013), where the solution to the initial value problem with data G(x), 0 ≤ x ≤ L, is given by the celebrated Hopfs-Lax formula (Hopf, 1970). For a generic time-space point P ≡ (t, x), this formula gives N(P), the cumulative vehicle number having crossed location x by time t, which in the case of the isosceles fundamental diagram (2) reads: subject to initial conditions N(0, x) = G(x), and where x U = x − t, x D = x + t are the spatial coordinates of points U and D in Fig. 1. Variational formulation (3) is solved by finding the minimum straight-like path from P to the initial boundary at (0, y), where the cost of each path, f P (y), is the sum of the terminal cost G(y) and the maximum number of vehicles that could possibly cross this path, (t − x + y)/2 in our case. With deterministic initial conditions at the critical density we have Figure 1: Correspondence between the area under the random walk and the space-time diagram: Sample spatiotemporal realization for ECA rule 184 under critical random initial conditions along with the cost function f P (x) which is a standard random walk. The colored bars represent the congestion clusters, the dark gray bars the voids and the light gray areas the capacity state.
so that the cost function becomes f P (y) = (L − (x − t))/2, which is independent of y. Therefore, all paths from P to the boundary are optimal (since they all have the same cost) and N det (P) ≡ (L − (x − t))/2 is the deterministic solution.
The important point here is that when G(y) is a stochastic process, at the critical density we still have that its expected value E {G(y)} = (L − y)/2, and that the expected cost is y-independent: This y-independence implies that the random process f P (y) is stationary and with 0 drift. Assuming independent increments of initial conditions leads directly to a Brownian motion process for f P (y), and will be the working assumption in this paper. Notice that the factor "1/2" in (3), (4) and (5) is the critical density and implies that f P (y) is normal with mean (L − (x − t))/2 and variance (L − (x − t))/4. The independent increments assumption is convenient not only for tractability, but because among all stationary stochastic processes, the Brownian motion is akin to a first-order stochastic approximation. Therefore, the insights revealed by this solution is likely to encode most of the important features of correlated initial conditions (as we have corroborated; see the discussion section). As shown next, the tractability afforded by this assumption means that one can leverage many important results concerning Brownian motions to derive key properties of traffic flow.

Marginal distribution of the N-surface
From the variational formulation (3) it becomes clear that the distribution of the vehicle number N(P) is given by the minimum of the Brownian motion over the finite interval x U < y < x D . Conditional on the Brownian motion being zero at the beginning of the interval, f P (U) = 0, this minimum is known to have the half-normal distribution with parameter √ π/2(x D − x U ) in our case. Since x D − x U = 2t it follows that the probability density of N(P) conditional on f P (U) = n U is a shifted half-normal distribution given by: The random walk can be seen as a slice of the deviation surface at t = 0, and the area under the random walk as the projection of the surface into the (∆N, x)-plane; note that state C is lost in this projection. with and one can see that the stochastic solution is below the deterministic one, i.e. E {N(P) | n U } < N det (P), as expected from Jensen's inequality, but with the discrepancy √ t/π increasing only as ∼ √ (6) is illustrated by the dashed lines in Fig. 2.
The unconditional marginal distribution for N(P) is obtained by taking expectation of (6) with respect to the distribution of n U , which is normal with mean ν = (L − (x − t))/2 and variance (L − (x − t))/4: Unfortunately, analytical expressions for the moments of this distribution are unavailable at this time. Fig. 2(Left) illustrates the shape of this distribution at a fixed location and for different times. It can be seen that the distribution is much flatter than the half-normal distribution (6) (dashed lines in the figure), as a consequence of large variance of n U . A realization of the N-surface is shown in Fig. 2(Right), which depicts a 3D representation this surface in oblique coordinates where the third axis corresponds to deviations around the deterministic mean of the N-surface: ∆N(P) = N(P) − N det (P). The random walk can be seen as a slice of this surface at t = 0, and the area under the random walk as the projection of the surface into the (∆N, x)-plane; note that state C is lost in this projection.

Flow distribution
A full description of the N-surface distribution, Pr {N(P 1 ) < n 1 , N(P 2 ) < n 2 , . . .}, has proven difficult to obtain, from which the distribution of flows, q(P), could be derived. However, we have verified that the mean of the conditional (6) and unconditional (8) distributions are numerically very similar, and thus the mean flow can be approximated taking the time-derivative in (7), to obtain: where we have used q(t) instead of q(P) to emphasize that the flow is location-independent. Notice that this function starts at 0 at t = 1/π, increases rapidly to 0.4 by t = 10 and flattens out to converge asymptotically to the capacity 1/2. (As a reference, recall that the critical headway is 2 time units.) That the flow is expected to be low at earlier times is a consequence of the system being farther from equilibrium at earlier times, which translates into an increased activity in jam interactions that produce little flow. As time goes by these interactions resolve to produce flows closer to capacity. This phenomenon can be seen in Fig. 1, where it is clear that jam interactions are much more prevalent near the base of the triangle then near the bottom vertex.
To gain additional insights about the distribution of flows, note from Fig. 1 that the flow at each time step can be either 1/2 in the capacity state or 0 in the void and congested states. As shown in the figure, the duration of both the void and congested states at a given location, τ, does not have a heavy-tail distribution and has a mean of approximately 2.7, independent of L. This number is surprisingly small compared to the total number of elementary jams for a particular realization. Fig. 3 shows the (oblique) cumulative counts and flow time series at 3 different locations for a typical rollout. The roughness of the cumulative curve is apparent as a result of the "spiky" variations in the flow around the capacity value 1/2. As a first approximation to describe the stochastic processes for flows, a fractional Gaussian noise process was fitted for 2000 realizations of flows timeseries. Also, a fractional Brownian motion process was estimated for the cumulative counts. Each panel shows the estimated Hurst exponent and the histograms summarize the results across rollouts. The average Hurst exponent for the flows H = 0.81 is in line with the empirical results reported for Australian motorways during morning rush hour and nights (Chand et al., 2017).
These high Hurst exponents indicate a long memory in the time series as a result of long-range spatial correlations. In this case these correlations are due to the spatial interactions brought about by the traffic model, which imposes long periods of constant capacity flow interrupted by bursts of null flow. Notice that the fractal dimension of the flow time series can be approximated by D = 2 − H = 2 − 0.81 ≈ 1.2.

Elementary jams
Let us define an elementary jam as the collection of occupied cells connecting a point P located on the back of the queue and its corresponding downstream point D on the boundary; see Fig. 1. This is akin to a congested characteristic line connecting P to D but in discrete units. Points P in the back of the queue are characterized by a transition from free-flow to congestion at that point, which means that the cost of minimum paths P-U and P-D both solve (3) and therefore have identical cost. This means that the duration D − U in Fig. 1 is given by the first return time of our random walk f P (y), which are known to have a power law distribution with exponent 1/2 (see e.g. Schroeder, 2009, Kostinski andAmir, 2016, p.152). From Fig. 1 it is clear that the duration, H, of an elementary jam is H = (D − U)/2 and therefore retains the same power-law distribution with exponent 1/2: The consequences of this are far-reaching: both the mean and variance of elementary jam durations (and distance traveled) diverge to infinity, which renders standard statistical tools inadequate. Most notably, the law of large numbers does not apply: taking a sample of jam durations to estimate mean and variance is futile because these estimates will not coverage regardless of the sample size.

Cluster size distribution
We consider here a road segment divided into L cells, where the initial "density" in each cell i = 1, 2, . . . L can be either 0 or 1 with probabilities 1-p and p, respectively, independently of other cells (by the independent increments assumption); i.e., a Bernoulli process. Notice that p can also be interpreted as the mean density in the segment.
As mentioned to the introduction, (Laval, 2021) recently showed numerically that ECA 184 exhibits power-law behavior at the critical density, with cluster sizes, S , following a power-law distribution with an exponent close to one half: Pr This was shown by evolving ECA-184 for T = 5, 000 time periods starting with an initial row of L = 10, 000 cells, each with a probability p = 1/2 of being occupied. The time-space area A = LT/2 for cluster measurement is of triangular shape as in Fig. 1, where L = 2T to ensure that all points in A are able to reach the boundary. The tail distribution of the site percolation cluster size (Stauffer and Aharony, 2018) were plotted on a log-log plot so that the slope of the linear segment (if any) gives an idea of the exponent of the power law; see inset in Fig. 1 for p = 1/2. Piecewise linear regression, with 2 pieces, was used to estimate the power-law exponent in order to account for finitesize effects; see the 2 estimated lines in Fig. 4 and note that the power-law exponent is taken as the slope of the uppermost linear piece, the orange line in the figure.
To study the behavior near the critical density, here we expand the experiment in (Laval, 2021) to include a range of mean initial densities p. Fig. 4 shows the relationship between the estimated power-law exponent α and p in our experiment for L = 1, 000. The downwards trend is fairly linear around the critical density, as in second-order phase transitions, and can be described approximately as The sizable but decreasing interval width of the transition, ≈ ±0.1 at the critical density, is a strong indication that the exponent that best describes a particular realization is highly sensitive to the initial conditions. This makes prediction challenging, especially considering that L = 1, 000 is a fairly large road segment of around 10 kilometers, and that the interval width increases for smaller segments, e.g. ≈ ±0.2 for L = 100 as seen in Fig. 4. It is important to note that only near the critical density these power-law exponent are meaningful. Fig. 4(Bottom) shows 3 detailed outputs for selected values of p, where one can see that for p << 1/2 there is a clear power law relationship but with a very limited range, because in the free-flow jams are very small. For for p >> 1/2 the cluster distribution describes a piecewise-linear relationship, as a consequence of several clusters "percolating" to the opposite side of the triangle in the time-space plane.

Intuitive explanation of Eqn. (11)
In terms of elementary jams, a jam can be characterized by a collection of, say n, contiguous elementary jams. The size S of a jam can be seen as the sum of n correlated power laws with exponents 1/2: It is well known that for independent H i 's and for large values of S the power law approximation for S is accurate (Zaliapin et al., 2005), and maintains the same exponent, 1/2 in this case, regardless of the value of n. In light of (11) one concludes, perhaps surprisingly, that the correlations among the H i 's are not a problem that would alter the scaling exponent of cluster sizes. This can be understood intuitively thanks to (i) the correspondence between the random walk and the jam clusters unveiled in this paper, and (ii) the property of power laws that the conditional distribution Pr {H i > h|H i > a} ∼ h −1/2 is also a power law with the same exponent. To see the implication of (i), take the duration H i in Fig. 1 that corresponds to points P and D in the figure. One can see that the duration of the elementary jam above it, which is part of the same jam, is smaller than H i and therefore correlated, but the one below it, which is part of another jam, is far greater and uncorrelated given that this particular H i happens to begin/end at a local minimum of the random walk. We conclude that within a jam the H i 's are correlated, across jams they are not. This means that jam cluster sizes S are uncorrelated.
To see the implication of (ii) above, notice that the number of contiguous elementary jams n can be seen as a random stopping time for the sun in (13), and that, starting with the smallest element in that sum, each consecutive element will be greater than the previous one. Therefore, given property (ii) above, the sought result follows.

Aggregate measures of performance
Let Ψ ≡ A k(t, x) dA and Φ ≡ A q(t, x) dA be the total time traveled and total distance traveled, respectively, over a (discrete) time-space region A. By letting the symbol A also denote its area, Edie's (Edie, 1965) average density, flow and speed in A, k A , q A and v A , respectively, are given by: These general definitions can be greatly simplified in our case because the spatiotemporal solution to our problem only involves 3 traffic states: capacity state C with (flow, density) = (1/2, 1/2), the empty or void state O with (0,0) and the jam state J with (0, 1). The total area that each state occupies in A is indicated with a subscript, and they satisfy: Airy distribution Figure 5: PDF of the Airy distribution, g(a), 0 ≤ a ≤ 1, in blue. Although it is analytical, there is no closed-form expression for this distribution, only a complicated infinite sum (Majumdar and Comtet, 2005), which is not implemented by default in most libraries/packages (the first 24 terms in this sum were considered for making this figure). The shaded area corresponds to the area under the Inverse Gaussian Distribution with mean parameter 0.625416 and shape parameter 9.82985, which can be seen to provide an accurate (and fast) approximation.
The total travel time Ψ = A C 1 2 · dA + A J 1 · dA = A C /2 + A J , and the total distance traveled Φ = A C /2. Using (15) we have: and therefore A J < A/2 is bounded. Also, where p J = A J /A correspondence to the probability that a time-space point in A is congested. It can be seen that the total time-space area in congestion, A J , turns out to be the key measure of performance from which all others can be derived.

Moments and probability densities
The main result here is that for an area A = L 2 /4 of triangular shape as in Fig. 1, the area A J of congested clusters corresponds to half the area under a Brownian excursion, e(x): A Brownian excursion is a conditioned one dimensional Brownian motion over the time interval (0, L) such that its path starts and ends at the origin and it is constrained to stay positive in between. It is important to note that the area under the excursion is also given by the area between a Brownian bridge and its minimum value in (0, L). This second interpretation is precisely our context: (i) exactly at the critical density the number of empty and unoccupied cells are equal so that in our random walk f P (0) = f P (L), which corresponds to a Brownian bridge, and (ii) from Fig. 1 it can be seen that (half) the total area under the random walk and its minimum value is A J , since: where m is the total number of elementary jams in the system. It is important to highlight that this equation remains valid so long as the time-space region of analysis is as in Fig. 1, i.e. an isosceles triangle in the shape of the fundamental diagram, with area A = L 2 /4. This triangle also corresponds to the domain of dependence of the continuum traffic model, and contains all time-space points that can be affected by perturbations at the boundary. Unlike (13), the correlations among the H i 's in (19) play a significant role here because the sum across all elementary jams in the system is limited by the size of the system itself, and therefore its distribution will not be a power law, as shown next.
For the unit length road segment the area under a Brownian excursion, Z = 1 0 e(x)dx is known to have the Airy distribution (Janson, 2007, Majumdar and Comtet, 2005, Agranov et al., 2020; let g(a), 0 ≤ a ≤ 1 be its PDF. Fig.  5 plots the PDF of the Airy distribution, g(a), 0 ≤ a ≤ 1. Notice that it is not a heavy-tail distribution; (Janson, 2007) shows that its mean and standard deviation are given by E {Z} = √ π/2/2 ≈ 0.627, SD {Z} = √ 5/12 − π/8 ≈ 0.155. Although the Airy distribution is not exactly symmetric, one can approximate the interval where most of the realizations of Z would fall: Z ∈ 0.627 ± 0.155 = (0.472, 0.781).
Given the self-similarity properties of Brownian motions, we have that A J /L 3/2 has the same distribution as Z/2 (Majumdar and Comtet, 2005), and therefore is an invariant quantity in our problem. Most of the realizations of A J /L 3/2 would fall in the interval (0.236, 0.391). Similarly, the invariant probability of congestion is p J L 1/2 and has the same distribution as 2Z. Importantly, it also follows that: To correctly interpret these insights, which are not contradictory, recall that the analysis area A = L 2 /4 changes with the size of the segment. Therefore, doubling the size of the road increases the analysis area by a factor of 4, but the congested area A J only by a factor of 2 3/2 ≈ 2.8, implying economies of scale. It is not surprising then, that the probability of congestion actually decreases by a factor of 2 −1/2 ≈ 2.8/4 ≈ 0.71.
To obtain the PDF of A J one can use the above self-similarity arguments to show that The distribution of all other traffic flow variables will be given by different scalings of the Airy distribution. It follows from (17) that the PDFs for these variables are given by which only depend upon L 1/2 . It can be seen from (22a) that the Airy distribution can be interpreted as the likelihood that a particular point in A = L 2 /4 is congested on a segment of length L = 4.

Discussion
This paper has presented overwhelming evidence that traffic flow in individual road segments exhibits all the characteristics of critical systems subject to SOC. We have shown that this critical behavior can be traced back to the deep connection with Brownian motions formalized here, a theory well understood by now that you have important implications for traffic flow going forward; we have only scratched the surface here.
This critical behavior added to the tendency of the system to be in the critical state is what makes the case for SOC so compelling. This tendency is a consequence of our driving behavior as we seek to maximize travel speeds-subject to safety-and thereby maximizing the flow out of a jam. The corresponding SOC "sandpile" model is one where vehicles are added one at a time to an initially empty ring road and observing whether or not a slowdown is produced. If not, another vehicle is added until this happens. The first slow down gives an indication of the beginning of the critical region, and congestion clusters of small size are expected. We keep adding vehicles to the system until jams of all sizes are observed, which correspond to "avalanches" in SOC theory.
Starting with random initial conditions at the critical density, we have shown that on long segments congestion clusters propagate for a time (or distance) duration that follows a power law with exponent 1/2. Exponents less than unity are problematic because no moments of the random variable exist, such as the average jam propagation distance (or time) and its variance, which makes prediction virtually impossible as traditional analysis tools are bound to fail. In fact, jams (and voids) of any duration may be expected for no particular reason at the critical density, in the same way earthquakes of any magnitude are expected in some countries. An effective control method is therefore to meter the demand to prevent the density from reaching the critical region, and not to maximize throughput as in current practice. Traffic flow near the critical density is therefore a fractal object, a property that has not been exploited in the current literature. In the (t, x)-plane the basic fractal "generator" unit is a triangle similar in shape to the fundamental diagram, containing 3 traffic states: voids, capacity and jams; see Fig. 1. Within each such triangle all the properties and formulas put forward in this paper apply so long as L is interpreted as the base of the triangle. Similarly, due to the symmetry property of the kinematic wave model, all the results in this paper for the jam state J transfer identically to avoid state O. For example, the distribution of the duration of voids at a given location is identical to the distribution of the duration of jams at the same location; the total area of voids is also half of the area under a Brownian excursion, etc.
This paper has taken this symmetry one step further, by exploiting the connection with the theory of random walks, best illustrated using the traffic flow N-surface as in Fig. 2 (right). Notably, one can see that the area under the random walk can be interpreted as the projection of this surface into the (∆N, x)-plane, and that state C is lost in this projection. It becomes clear then, that delays at the critical density can be calculated simply by integration of the initial data, without the need to solve the traffic flow model explicitly.
In the case of flow timeseries, we found compelling evidence for a Hurst exponent H > 0.5 (so flow fractal dimension D < 3/2) indicating long-range correlations for both cumulative counts and flow timeseries, settling the discrepancies in the existing literature between conflicting results (Chand et al., 2017) and (Krause et al., 2017) in favor of the former. It remains to be established, however, that a fractional process approach is justified for traffic flow. From Fig. (3) it is clear that the flow time series is far from a random noise symmetric around its mean value, as required by the fractional process theory. A better description of cumulative accounts and flow time-series can be obtained through fractal geometry, in particular, the relationship E {A J } ∼ L 3/2 unveiled here. This relationship implies a fractal dimension of 3/2 for the cluster areas in the time-space plane, and therefore a fractal dimension of 3/2+1=5/2 for the N-surface, and a fractal dimension of 5/2-1=3/2 for any slice of this surface, such as the cumulative curves or flow time-series.
As mentioned in the introduction, early work (Nagel andPaczuski, 1995, Paczuski and conjectured that SOC in single-lane models emerges only in "cruise-control" type traffic flow models, where random perturbations can only take place in congested traffic states. Here we have found that SOC appears to be an intrinsic property of traffic flow, giving a strong indication that the results presented here should apply in general, as a consequence of our conjecture of universality of the kinematic wave model/ECA-184. These references also pointed out the connection with random walks, based on numerical results that hinted to a power nap exponent close to one half for the duration of a jam. Here we have put this results on firmer grounds thanks to the connection with variational theory presented in section 2.1. A discrepancy found with this earlier work is in the distribution of jem sizes (11), where they found α = 1/3 as opposed to one half as shown here. Their result is a direct consequence of their conjecture that the number of vehicles under congestion scales with time as t 1/2 , for which we found no support here.
In this paper we have concentrated on initial conditions of the "white noise" type, which is the simplest paradigm in stochastic analysis as it assumes stationary and independent increments. The next step is to analyze more general initial conditions that allow for spatial correlations, such as fractional Gaussian noise processes with Hurst exponent H, where H = 1/2 corresponds to the Brownian white noise. It is important to note that this parameter also determines the occupation probability p. This is shown in Fig. 6 (top), which shows the relationship between the occupation probability p and Hurst exponent H obtained numerically. It can be seen that for H < 1/2 the occupation probability follows a narrow normal distribution around a mean of 1/2, and that for H > 1/2 this distribution widens as H → 1. This is as expected since long-range dependencies imply that two neighboring cells will tend to have the same value (either occupied or not) as H → 1. The bottom part of the figure shows the power-law exponents and α obtained with these fractional initial conditions and the probation probability p as a function of H. The main conclusion here is that the critical region seems to be wider than in the random walk case, i.e. that apart from the mean initial density p, the spatial correlations in the initial density are also an important factor when it comes to critical behavior. This is currently under investigation by the author. So, is it time for paradigm change? Most likely yes, at least in the critical region. We have seen that prediction of "local" variables such as individual jam propagation and its effects on flow time-series is virtually impossible in the critical region. The main difficulty here is that traffic flow is "chaotic" in the sense that jam propagation patterns are highly sensitive to initial conditions, i.e. to the individual vehicle arrangement at the beginning of the analysis period. Importantly, trying to predict when a particular road segment will spill back is a lost cause with standard loopdetector technology, but might be achievable with vision technology tracking each individual vehicle. But even then it is possible that random fluctuations will prevent accurate predictions. In the meantime, a robust control strategy at the network level appears to be perimeter control to prevent the system for reaching the critical region, as opposed to current practice to keep vehicle accumulation at the "sweet spot" to maximize flows. It turned out that there is nothing sweet about the critical region! That is where catastrophes take place seemingly out of nowhere: e.g., a vehicle not willing to reroute due to a spill back in front, or another one unable to change lanes in time and forced to stop, or a delivery truck suddenly double parking, all might trigger severe congestion in the whole network.
We are therefore in a position now to accept Nagel and Paczuki's 1995 conjecture that traffic management strategies seeking to maximize flow are bound to cause more harm than good. The main reason why flow-maximizing strategies have not proven problematic in the literature might be explained because of the traffic models used, typically some variation of the cell transmission model. While it is true that this model is a numerical solution method for the kinematic wave model, it is unable to capture the nuances of the sensitivity to initial conditions exposed here, which are wiped-out by the mesoscopic averaging over typically long cells and by the considerable numerical error of this model. This argument makes the case for the use of microscopic models (where vehicles do not disappear when unable to take an exit!), but at the expense of increased computation times. Another research direction would be mesoscopic models that capture critical behavior, which appeared to be nonexistent today.