# Subcritical Hopf bifurcations in a car-following model with reaction-time delay

Gábor Orosz , Gábor Stépán

## Abstract

A nonlinear car-following model of highway traffic is considered, which includes the reaction-time delay of drivers. Linear stability analysis shows that the uniform flow equilibrium of the system loses its stability via Hopf bifurcations and thus oscillations can appear. The stability and amplitudes of the oscillations are determined with the help of normal-form calculations of the Hopf bifurcation that also handles the essential translational symmetry of the system. We show that the subcritical case of the Hopf bifurcation occurs robustly, which indicates the possibility of bistability. We also show how these oscillations lead to spatial wave formation as can be observed in real-world traffic flows.

Keywords:

## 1. Introduction

The so-called uniform flow equilibrium of vehicles following each other on a road is a kind of steady state, where equidistant vehicles travel with the same constant velocity. Ideally, this state is stable. Indeed, it is the goal of traffic management that drivers choose their traffic parameters to keep this state stable and also to reach their goal, that is, to travel with a speed close to their desired speed. Still, traffic jams often appear as congestion waves travelling opposite to the flow of vehicles (Kerner 1999). The formation of these traffic jams (waves) is often associated with the linear instability of the uniform flow equilibrium, which should be a rare occurrence. However, it is also well known among traffic engineers that certain events, such as a truck pulling out of its lane, may trigger traffic jams even when the uniform flow is stable. We investigate a delayed car-following model and provide a thorough examination of the subcriticality of Hopf bifurcations related to the drivers' reaction-time delay. This explains how a stable uniform flow can coexist with a stable traffic wave.

The car-following model analysed in this paper was first introduced in Bando et al. (1995) without the reaction-time delay of drivers, and that was investigated by numerical simulation. Then, numerical continuation techniques were used in Gasser et al. (2004) and Berg & Wilson (2005) by applying the package Auto (Doedel et al. 1997). Recently, Hopf calculations have been carried out in Gasser et al. (2004) for arbitrary numbers of cars following each other on a ring.

The reaction-time delay of drivers was first introduced in Bando et al. (1998), and its importance was then emphasized by the study of Davis (2003). In these papers, numerical simulation was used to explore the nonlinear dynamics of the system. The first systematic global bifurcation analysis of the delayed model (Davis 2003) was presented in Orosz et al. (2004), where numerical continuation techniques, namely the package Dde-Biftool (Engelborghs et al. 2001), were used. The continuation results were extended to large numbers of cars in Orosz et al. (2005), where the dynamics of oscillations, belonging to different traffic patterns, were analysed as well. In this paper, we perform an analytical Hopf bifurcation calculation and determine the criticality of the bifurcation as a function of parameters for arbitrary numbers of vehicles in the presence of the drivers' reaction delay.

While the models without delay are described by ordinary differential equations (ODEs), presenting the dynamics in finite-dimensional phase spaces, the appearance of the delay leads to delay differential equations (DDEs) and to infinite-dimensional phase spaces. The finite-dimensional bifurcation theory that is available in basic textbooks (Guckenheimer & Holmes 1997; Kuznetsov 1998) have been extended to DDEs in Hale & Verduyn Lunel (1993), Diekmann et al. (1995), Kolmanovskii & Myshkis (1999) and Hale et al. (2002). The infinite-dimensional dynamics make the bifurcation analysis more abstract. In particular, the Hopf normal form calculations require complicated algebraic formalism and algorithms, as is shown in Hassard et al. (1981), Stépán (1986, 1989), Campbell & Bélair (1995), Orosz (2004) and Stone & Campbell (2004). Recently, these Hopf calculations have been extended for systems with translational symmetry in Orosz & Stépán (2004), which is an essential property of car-following models. This situation is similar to the S1 symmetry that occurs in laser systems with delay, see Verduyn Lunel & Krauskopf (2000) and Rottschäfer & Krauskopf (2004). In Orosz & Stépán (2004), the method was demonstrated in the over-simplified case of two cars on a ring. Here, these calculations are extended to arbitrarily many cars, providing general conclusions for the subcriticality of the bifurcations and its consequences for flow patterns. Our results are generalization of those in Gasser et al. (2004) for the case without reaction-time delay. In particular, we prove that this delay makes the subcriticality of Hopf bifurcations robust.

## 2. Modelling issues

The mathematical form of the car-following model in question was introduced and non-dimensionalized in Orosz et al. (2004). Here, we recall the basic features of this model. Periodic boundary conditions are considered, i.e. n vehicles are distributed along a circular road of overall length L; see figure 1. (This could be interpreted as traffic on a circular road around a large city, e.g. the M25 around London, even though such roads usually possess higher complexity.) As the number of cars is increased, the significance of the periodic boundary conditions usually tends to become smaller.

Figure 1

Vehicles flowing clockwise on a circular road.

We assume that drivers have identical characteristics. Considering that the ith vehicle follows the (i+1)th vehicle and the nth car follows the 1st one, the equations of motion can be expressed as(2.1)where the dot stands for time derivative. The position, the velocity, and the acceleration of the ith car are denoted by xi, and , respectively. The so-called optimal velocity function depends on the distance of the cars , which is usually called the headway. The argument of the headway contains the reaction-time delay of drivers which now is rescaled to 1. The parameter α>0 is known as the sensitivity and 1/α>0 is often called the relaxation time. Due to the rescaling of the time with respect to the delay, the delay parameter is hidden in the sensitivity α and the magnitude of the function V(h); see details in Orosz et al. (2004).

Equation (2.1) expresses that each driver approaches an optimal velocity, given by V(h), with a characteristic relaxation time of 1/α, but reacts to its headway via a reaction-time delay 1. The general features of the optimal velocity function V(h) can be summarized as follows:

1. V(h) is continuous, non-negative and monotone increasing, since drivers want to travel faster as their headway increases. Note that in the vicinity of the Hopf bifurcation points, V(h) is required to be three times differentiable for the application of the Hopf theorem (Guckenheimer & Holmes 1997; Kuznetsov 1998).

2. as h→∞, where is known as the desired speed, which corresponds, for example, to the speed limit of the given highway. Drivers approach this speed with the relaxation time of 1/α when the traffic is sparse.

3. V(h)≡0 for h∈[0,1], where the so-called jam headway is rescaled to 1, see Orosz et al. (2004). This means that drivers attempt to come to a full stop if their headways become less than the jam headway.

One might, for example, consider the optimal velocity function to take the form(2.2)as already used in Orosz et al. (2004, 2005). This function is shown together with its derivatives in figure 2. Functions with shapes similar to equation (2.2) were used in Bando et al. (1995, 1998) and Davis (2003).

Figure 2

(a) The optimal velocity function (2.2) and (b–d) its derivatives.

Note that the analytical calculations presented here are valid for any optimal velocity function V(h): it is not necessary to restrict ourself to a concrete function in contrast to the numerical simulations in Bando et al. (1998) and Davis (2003) and even to the numerical continuations in Orosz et al. (2004, 2005).

The dimensionless parameters α and v0 can be obtained from their dimensional counterparts as shown in Orosz et al. (2004). If we consider the reaction-time delay 1–2 s and the relaxation time 1–50 s, we obtain the sensitivity as their ratio: α∈[0.02,2]. Also, if we assume desired speeds in the range 10–40 m s−1 and jam headways 5–10 m, their ratio times the reaction time provides the dimensionless desired speed v0∈[1,20]. The linear stability diagram does not change qualitatively when v0 is varied in this range as shown in Orosz et al. (2005).

## 3. Translational symmetry and Hopf bifurcations

The stationary motion of the vehicles, the so-called uniform flow equilibrium, is described by(3.1)where(3.2)and(3.3)Note that one of the constants can be chosen arbitrarily due to the translational symmetry along the ring. Henceforward, we consider the average headway h*=L/n as a bifurcation parameter. Increasing h* increases the length L of the ring, which involves scaling all headways hi accordingly.

Let us define the perturbation of the uniform flow equilibrium by(3.4)Using Taylor series expansion of the optimal velocity function V(h) about up to third order of , we can eliminate the zero-order terms(3.5)where(3.6)At a critical/bifurcation point , the derivatives take the values , and . Now, and further on, prime denotes differentiation with respect to the headway.

Introducing the notation(3.7)equation (3.5) can be rewritten as(3.8)where . The matrices and the near-zero analytic function are defined as(3.9)Here, stands for the n-dimensional identity matrix, while the matrix and the functions are defined as(3.10)

System (3.8) possesses a translational symmetry, therefore, the matrices satisfy(3.11)that is, the Jacobian has a zero eigenvalue(3.12)for any value of parameter h*. Furthermore, the near-zero analytic function preserves this translational symmetry, that is,(3.13)for all c≠0 satisfying ; see details in Orosz & Stépán (2004).

The steady state y(t)≡0 of equation (3.8) corresponds to the uniform flow equilibrium (3.1) of the original system (2.1). Considering the linear part of (3.8) and using the trial solution with and , the characteristic equation becomes(3.14)According to equation (3.11), the relevant zero eigenvalue (3.12) is one of the infinitely many characteristic exponents that satisfy equation (3.14). This exponent exists for any value of the parameter b1, that is, for any value of the bifurcation parameter h*.

At a bifurcation point defined by , i.e. by , Hopf bifurcations may occur in the complementary part of the phase space spanned by the eigenspace of the zero exponent (3.12). Then, there exists a complex conjugate pair of pure imaginary characteristic exponents(3.15)which satisfies equation (3.14). To find the Hopf boundaries in the parameter space, we substitute into equation (3.14). Separation of the real and imaginary parts gives(3.16)The so-called wave number k can take the values (for even n) and (for odd n). The wave numbers are not considered since they belong to conjugated waves producing the same spatial pattern. Note that k=1 belongs to the stability boundary of the uniform flow equilibrium, while the cases k>1 result in further oscillation modes around the already unstable equilibrium. Furthermore, for each k the resulting frequency is bounded so that ; see Orosz et al. (2004, 2005).

Note also that the function shown in figure 2b is non-monotonous, and so a b1cr boundary typically leads either to two or to zero boundaries. For a fixed wave number k, these boundaries tend to finite values of h* as , as is shown in Orosz et al. (2005). Using trigonometric identities, equation (3.16) can be transformed to(3.17)which is a useful form used later in the Hopf calculation together with the resulting form(3.18)

With the help of the identity(3.19)we can calculate the following necessary condition for Hopf bifurcation as the parameter b1 is varied as(3.20)where(3.21)Since equation (3.20) is always positive, this Hopf condition is always satisfied. Now, using the chain rule and definition (3.6), condition (3.20) can be calculated further as the average headway h* is varied to give(3.22)This condition is fulfilled if and only if , which is usually satisfied except at some special points. For example, the function shown in figure 2c becomes zero at a single point over the interval . Notice that is zero for and for , but the critical headway hcr never takes these values for α>0. Conditions (3.20) and (3.22) ensure that the characteristic roots (3.15) cross the imaginary axis with a non-zero speed when the parameters b1 and h* are varied.

## 4. Operator differential equations and related eigenvectors

The DDE (3.8) can be rewritten in the form of an operator differential equation (OpDE) which reflects its infinite-dimensional dynamics. For the critical bifurcation parameter , we obtain(4.1)where the dot still refers to differentiation with respect to the time t and is defined by the shift , on the function space of continuous functions mapping . The linear and nonlinear operators , are defined as(4.2)(4.3)where the matrices and the nonlinear function are given by(4.4)

The translational symmetry conditions (3.11) and (3.13) are inherited, that is,(4.5)and(4.6)for all c≠0 satisfying (L+R)c=0.

In order to avoid singularities in Hopf calculations, we follow the methodology and algorithm of Orosz & Stépán (2004) and eliminate the eigendirection belonging to the relevant zero eigenvalue λ0=0 (equation (3.12)). The corresponding eigenvector satisfies(4.7)which, applying the definition (4.2) of operator , leads to a boundary value problem with the constant solution(4.8)One finds that(4.9)where each component of the vector is equal to 1. Here, is a scalar that can be chosen freely, in particular, we choose(4.10)

In order to project the system to s0 and to its complementary space, we also need the adjoint operator(4.11)where asterisk denotes either adjoint operator or transposed conjugate vector and matrix. The eigenvector of associated with the eigenvalue satisfies(4.12)This gives another boundary value problem, which has the constant solution(4.13)Here, we obtain(4.14)However, is not free, but is determined by the normality condition(4.15)Defining the inner product(4.16)condition (4.15) gives the scalar equation(4.17)from which we obtain(4.18)

Note that, as the vectors s0 and n0 are the right and left eigenvectors of the operator belonging to the eigenvalues and , similarly, the vectors S0 and N0 are the right and left eigenvectors of the matrix , belonging to the same zero eigenvalues. For the non-delayed model (Bando et al. 1995), the vectors S0 and N0 are the left and right eigenvectors of the Jacobian (L+R), that is, their structure is related the spatial structure of the system along the circular road.

We are now able to eliminate the singular dynamics generated by the translational symmetry so that we project the system to s0 and to its complementary space. With the help of the eigenvectors s0 and n0, the new state variables and are defined as(4.19)Using the derivation given in Orosz & Stépán (2004), we split the OpDE (4.1) into(4.20)Its second part is already fully decoupled, and can be redefined as(4.21)where the new nonlinear operator assumes the form(4.22)Considering given by equation (4.3), and using the expressions (3.9), (3.10) and (4.4), and the eigenvectors (4.9) and (4.14), we obtain(4.23)where the definition is applied. Note that the system reduction related to the translational symmetry changes the nonlinear operator of the system while the linear operator remains the same. This change has an essential role in the centre-manifold reduction presented below.

Let us consider a Hopf bifurcation at a critical point . First, we determine the real and imaginary parts of the eigenvector of the linear operator , which belongs to the critical eigenvalue (3.15), that is,(4.24)After the substitution of definition (4.2) of , the solution of the resulting boundary value problem can be written as(4.25)with constant vectors satisfying the homogeneous equation(4.26)

Using formula (3.17) and following the steps (A 1)–(A 3) of appendix A, one may solve (4.26) and obtain(4.27)where the scalar parameters u and υ can be chosen arbitrarily and the vectors C, are(4.28)with the wave number k used in equations (3.16) and (3.17). The cyclic permutation of the components in C and S results in further vectors S1 and S2 that still satisfy equation (4.26). This result corresponds to the symmetry of the system, that is, all the cars have the same dynamic characteristics. Choosing u=1 and υ=0 yields(4.29)

The real and imaginary parts of the eigenvector of the adjoint operator associated with are determined by(4.30)

The use of definition (4.11) of leads to another boundary value problem having the solution(4.31)where the constant vectors satisfy(4.32)

Applying equation (3.17) and following the steps (A 4)–(A 7) of appendix A, the solution of equation (4.32) is obtained as(4.33)The scalar parameters , are determined by the orthonormality conditions(4.34)which using the inner product definition (4.16), results in the two scalar equations(4.35)Substituting equations (3.17), (4.29) and (4.33) into (4.35) and using the second-order trigonometric identities (B 3)–(B 5) of appendix B, we obtain(4.36)with the solution(4.37)where is defined by equation (3.21). The substitution of equation (4.37) into (4.33) leads to(4.38)

As and are the right and left eigenvectors of the operator belonging to the eigenvalues and , the vectors and are similarly the right and left eigenvectors of the matrix belonging to the same eigenvalues. Note that for the non-delayed model (Bando et al. 1995) the vectors and are the left and right eigenvectors of the Jacobian , that is, their structure is again related the spatial structure of the system.

## 5. Centre-manifold reduction

Now, we investigate the essential dynamics on the two-dimensional centre manifold embedded in the infinite-dimensional phase space. Since the eigenvectors s1 and s2 span a plane tangent to the centre manifold at the origin, we use s1, s2 and n1, n2 to introduce the new state variables(5.1)where and . Using the derivation presented in Orosz & Stépán (2004), we can reduce OpDE (4.21) to the form(5.2)The scalar parameters are induced by the translational symmetry and their expressions are determined in Orosz & Stépán (2004) as(5.3)Since in our case , we obtain(5.4)

The power series form of (5.2) is also given in Orosz & Stépán (2004) as(5.5)where the subscripts of the constant coefficients and the vector ones refer to the corresponding jth and kth orders of z1 and z2, respectively. The terms with the coefficients and come from the linear combinations of s1(ϑ) and s2(ϑ). The translational symmetry only enters through the terms with coefficients , so that the terms with coefficients and refer to the structure of the modified nonlinear operator (equation (4.22)). Using the third- and fourth-order trigonometric identities (B 6)–(B 11) of appendix B, we can calculate these coefficients for wave numbers kn/2, kn/3 and kn/4 as given by equations (A 8) and (A 9) in appendix A.

Note that the cases , and result in different formulae for the above coefficients, but the final Poincaré–Lyapunov constant will have the same formula as in the case of general wave number k. The detailed calculation of these ‘resonant’ cases is not presented here.

Approximate the centre manifold locally as a truncated power series of w depending on the coordinates z1 and z2 as(5.6)There are no linear terms since the plane spanned by the eigenvectors s1 and s2 is tangent to the centre manifold at the origin. Third and higher order terms are dropped. The unknown coefficients can be determined by calculating the derivative of w and substituting that into the third equation of (5.5). The solution of the resulting linear boundary value problem, given in detail in Orosz & Stépán (2004), is(5.7)where the vectors satisfy(5.8)Here, we also used that(5.9)in accordance with equation (A 8).

One can find that the 2n-dimensional equation for H0 is decoupled from the 4n-dimensional equation for H1, H2 in equation (5.8). Since (L+R) is singular due to the translational symmetry (4.5), the non-homogeneous equation for H0 in (5.8) may seem not to be solvable. However, its right-hand side belongs to the image space of the coefficient matrix (L+R) due to the translational symmetry induced terms . We obtain the solution(5.10)with the undetermined parameter κ that has no role in the following calculations.

At the same time, the non-homogeneous equation for H1, H2 in equation (5.8) is not effected by the vectors . Using equation (3.17) and following the steps (A 10)–(A 14) of appendix A one can find the solution(5.11)where(5.12)

Now, using equation (5.11) in (5.6) and (5.7), we can calculate Rw(−1) which appears in the first two equations of (5.5). In this way, we obtain the flow restricted onto the two-dimensional centre manifold described by the ODEs(5.13)where the coefficients have already been determined by equation (A 8) in (5.5), while the coefficients(5.14)originate in the terms involving Rw(−1). To determine equation (5.14), the trigonometric identities (B 3)–(B 11) of appendix B have been used. The coefficients are not displayed, since they have no role in the following calculations.

We note that the coefficients of the second-order terms are not changed by the centre-manifold reduction. The so-called Poincaré–Lyapunov constant in the Poincaré normal form of equation (5.13) can be determined by the Bautin formula (Stépán 1989)(5.15)The bifurcation is supercritical for negative and subcritical for positive values of Δ. We found that Δ>0 is always true when which is the case for real traffic situations (many vehicles n with a few waves k). This can be proven as is detailed below.

The first part of the expression (5.15) of Δ in front of the parenthesis is always positive, since . Within the parenthesis in equation (5.15), the first term is positive since equation (3.16) implies , which yields critical headway values such that ; see figure 2b,d. The second term in the parenthesis in equation (5.15) contains the ratio of two complicated expressions, which by using equations (3.16) and (5.12), can be rearranged in the form(5.16)and(5.17)Since equation (5.17) is always positive, the sign of equation (5.16) is crucial for deciding the overall sign of Δ. According to equation (3.16) , that is, the realistic case implies the oscillation frequency .

Figure 3a shows the numerator for some particular values of α demonstrating that for . Note that if α→0 then may become negative (see figure 3a for ), but this is a physically unrealistic case where drivers intend to reach their desired speed extremely slowly.

Figure 3

Quantities defined by (5.16) and (5.17) as a function of the frequency ω, for representative values of parameter α. (a) Numerator and (b) ratio .

Moreover, the ratio of equations (5.16) and (5.17), , is not only positive for but also when ω→0 (i.e. when n→∞) as shown in figure 3b. This feature provides robustness for subcriticality. Note that subcriticality also occurs for optimal velocity functions different from equation (2.2), e.g. for those that are considered in Orosz et al. (2004).

Using definition (3.6), formulae (3.18) and (3.22) and expressions (5.15)–(5.17), the amplitude A of the unstable oscillations is obtained in the form(5.18)Thus, the first Fourier term of the oscillation restricted onto the centre manifold is(5.19)Since close to the critical bifurcation parameter , we have , equation (5.19) yields(5.20)where the vectors S1, S2 are given in equation (4.29).

The parameter enters Δ and A via the derivatives , and which are all proportional to v0. Consequently, Δ depends linearly on v0 causing no sign change and v0 disappears from A. Furthermore, v0 is embedded in the critical parameter which is determined from equation (3.16) by inverting . However, one may check that this is relevant for only, when small v0 may result in supercriticality as demonstrated in Orosz et al. (2004). In contrast, the realistic case leads to robust subcriticality as explained in §6 and also demonstrated in Orosz et al. (2005).

Note that zero reaction time delay results in as shown in Gasser et al. (2004). In that case, subcriticality appears only for extremely high values of the desired speed v0 when the term 6b3cr becomes greater than at the critical points (of the non-delayed model). Consequently, the presence of the drivers' reaction-time delay has an essential role in the robustness of the subcritical nature of the Hopf bifurcation. This subcriticality explains how traffic waves can be formed when the uniform flow equilibrium is stable, as is detailed in the subsequent section.

## 6. Physical interpretation of results

The unstable periodic motion given in equation (5.20) corresponds to a spatial wave formation in the traffic flow, which is actually unstable. Substituting (4.29) into (5.20) and using definition (3.7), one can determine the velocity perturbation as(6.1)The interpretation of this perturbation mode is a wave travelling opposite to the car flow with spatial wave number k (i.e. with spatial wavelength ). The related wave speed is(6.2)where the elimination of the frequency ω with the help of equation (3.18) leads to(6.3)Since the uniform flow equilibrium (3.1) travels with speed , the speed of the arising wave is(6.4)

By considering the optimal velocity function (2.2), we obtain cwave<0, that is, the resulting wave propagates in the opposite direction to the flow of vehicles. Note that the non-delayed model introduced in Bando et al. (1995) exhibits the same wave speed apart from some differences in the coefficient of the correction term . If one neglects this correction term, the wave speed becomes independent of n and k, which corresponds to the results obtained from continuum models; see e.g. Whitham (1999).

In order to check the reliability of the Poincaré–Lyapunov constant (5.15) and the amplitude estimation (5.18), we compare these analytical results with those obtained by numerical continuation techniques with the package Dde-Biftool (Engelborghs et al. 2001). In figure 4a, we demonstrate the subcriticality for n=9 cars and k=1 wave. The horizontal axis corresponds to the uniform flow equilibrium, that is, stable for small and large values of h* (shown by green solid line) but unstable for intermediate values of h* (shown by red dashed line) in accordance with formula (3.16) and figure 2b. The Hopf bifurcations, where the equilibrium loses its stability, are marked by blue stars. The branches of the arising unstable periodic motions given by equation (5.18) are shown as red dashed curves.

Figure 4

The amplitude A of velocity oscillations as a function of the average headway parameter h* (a) and the corresponding velocity profiles at h*=2.9 (b) and (c) for n=9 cars, k=1 wave and parameters α=1.0, v0=1.0. (a) Horizontal axis (A≡0) represents the uniform flow equilibrium and the analytical results are coloured: green solid and red dashed curves represent stable and unstable branches, respectively, and blue stars stand for Hopf bifurcations. Grey curves correspond to numerical continuation results: solid and dashed curves refer to stable and unstable states and grey crosses represent fold bifurcations. The points marked by and refer to the velocity profiles shown in panels (b) and (c), respectively. (b) Stable stop-and-go oscillations are shown for point P1 in green and for point in grey and (c) unstable oscillations are displayed for point P2 in red and for point in grey.

In §2, conditions (i)–(iii) require that the optimal velocity function is bounded so that ; see figure 2a. Maximum principles show that for t→∞ the velocity of any solution is contained by the interval [0, v0]. This suggests that ‘outside’ the unstable oscillating state there exists an attractive oscillating state which may include stopping, accelerating, travelling with the desired speed v0, and decelerating. The amplitude of this solution is defined as . In figure 4a the horizontal green line at represents this stable stop-and-go oscillation. The corresponding stop-and-go wave propagates against the traffic flow, since vehicles leave the traffic jam at the front and enter it at the rear; see figure 1. The above heuristic construction reveals the existence of bistability on each side of the unstable equilibrium between the Hopf bifurcation point and the point where the branches of unstable and stable oscillations intersect each other (outside the frame in figure 4a). For such parameters, depending on the initial condition, the system either tends to the uniform flow equilibrium or to the stop-and-go wave.

In figure 4a, we also displayed the results of numerical continuation carried out with the package Dde-Biftool (Engelborghs et al. 2001). Grey solid curves represent stable oscillations while grey dashed curves represent unstable ones. The fold bifurcation points, where the branches of stable and unstable oscillations meet, are marked by grey crosses. The comparison of the results shows that the analytical approximation of the unstable oscillations is quantitatively reliable in the vicinity of the Hopf bifurcation points. The heuristic amplitude v0/2 of the stop-and-go oscillations is slightly larger than the numerically computed ones. The analytically suggested bistable region is larger than the computed one (between the Hopf point and the fold point), since the third degree approximation is not able to predict fold bifurcations of periodic solutions. In order to find these fold bifurcation points, it is necessary to use numerical continuation techniques as presented in Orosz et al. (2005). Nevertheless, qualitatively the same structure is obtained by the two different techniques.

As was already mentioned in §3, the wave numbers k>1 are related to Hopf bifurcations in the parameter region, where the uniform flow equilibrium is already unstable. This also means that the corresponding oscillations for k>1 are unstable independently of the criticality of these Hopf bifurcations. Still, we found that these Hopf bifurcations are all robustly subcritical for any wave number k (except for large ). Consequently, the only stable oscillating state is the stop-and-go motion for k=1. On the other hand, several unstable solutions may coexist as is explained in Orosz et al. (2005). Note that analytical and numerical results agree better as the wave number k is increased because the oscillating solution becomes more harmonic.

To represent the features of vehicles' motions, the velocity oscillation profiles of the first vehicle are shown in figure 4b,c for the points and marked in figure 4a, for headway h*=2.9. Again, the coloured curves correspond to the analytical results, while the grey curves are obtained by numerical continuation. In figure 4b,c, the time window of each panel is chosen to be the period of the first Fourier approximation given by equation (3.18) (red curve in figure 4c) and the dashed vertical lines indicate oscillation periods computed numerically with Dde-Biftool.

Figure 4b shows the stop-and-go oscillations. The heuristic construction (green curve) is obtained by assuming that the stopping and flowing states are connected with states of constant acceleration/deceleration, which is qualitatively a good approximation of the numerical result (grey curve). Figure 4c compares the unstable periodic motions computed analytically (red curve) from the Hopf calculation with those from numerical continuation (grey curve).

For a perturbation ‘smaller’ than the unstable oscillation, the system approaches the uniform flow equilibrium. If a larger perturbation is applied then the system develops stop-and-go oscillations and a spatial stop-and-go travelling wave appears as demonstrated in figure 1. Since the period of the stable and unstable oscillations are close to each other, the stable stop-and-go wave travels approximately with the speed of the unstable travelling wave; see equation (6.4) for k=1.

## 7. Conclusion

A nonlinear car-following model has been investigated with special attention paid to the reaction-time delay of drivers. By considering the average headway as a bifurcation parameter, Hopf bifurcations were identified. In order to investigate the resulting periodic motions, the singularities related to the essential translational symmetry had to be eliminated. Then the Hopf bifurcations were found to be robustly subcritical leading to bistability between the uniform flow equilibrium and a stop-and-go wave. The appearing oscillations manifest themselves as spatial waves propagating backward along the circular road.

In the non-delayed model of Bando et al. (1995), subcriticality and bistability occur only for extremely high values of the desired speed v0, as it is demonstrated in Gasser et al. (2004). We proved that subcriticality and bistability are robust features of the system due to the drivers' reaction-time delay, even for moderate values of the desired speed. This delay, which is smaller than the macroscopic time-scales of traffic flow, plays an essential role in this complex system because it changes the qualitative nonlinear dynamics of traffic.

Due to the subcriticality, stop-and-go traffic jams can develop for large enough perturbations even when the desired uniform flow is linearly stable. These perturbations can be caused, for example, by a slower vehicle (such us a lorry) joining the inner lane flow for a short-time interval via changing lanes. It is essential to limit these unwanted events, for example, by introducing temporary regulations provided by overhead gantries. Still, if a backward travelling wave shows up without stoppings, it either dies out by itself or gets worse ending up as a persistent stop-and-go travelling wave. In order to dissolve this undesired situation, an appropriate control can be applied using temporary speed limits given by overhead gantries that can lead the traffic back ‘inside’ the unstable travelling wave and then to reach the desired uniform flow. For example, the MIDAS system (Lunt & Wilson 2003) installed on the M25 motorway around London is able to provide the necessary instructions for drivers.

## Acknowledgments

The authors acknowledge with thanks discussions with and comments of Eddie Wilson and Bernd Krauskopf on traffic dynamics and on numerical bifurcation analysis. This research was supported by the University of Bristol under a Postgraduate Research Scholarship and by the Hungarian National Science Foundation under grant no. OTKA T043368.

## Footnotes

• Alternative address: Research Group on Dynamics of Vehicles and Machines, Hungarian Academy of Science, PO Box 91, Budapest 1521, Hungary.