Chaotic Ising-like dynamics in traffic signals

The green and red lights of a traffic signal can be viewed as the up and down states of an Ising spin. Moreover, traffic signals in a city interact with each other, if they are controlled in a decentralised way. In this paper, a simple model of such interacting signals on a finite-size two-dimensional lattice is shown to have Ising-like dynamics that undergoes a ferromagnetic phase transition. Probabilistic behaviour of the model is realised by chaotic billiard dynamics that arises from coupled non-chaotic elements. This purely deterministic model is expected to serve as a starting point for considering statistical mechanics of traffic signals.

them, assuming that each signal is controlled on the basis of locally observed information on traffic flows. Probabilistic behaviour resembling the Ising model arises from these deterministic interactions.

Results
Traffic signal model. Specifically, the model is formulated as follows. Traffic signals are located at every crossing (node) of a two-dimensional lattice (Fig. 1a) of size L 3 L with a periodic boundary condition. Roads exist between neighbouring nodes, and each road consists of two lanes, one in each direction. We regard vehicles as fluids that flow along the lanes; they flow into a traffic lane at a crossing and flow out at a neighbouring crossing. As shown in Fig. 1b, the signals at the ith node have one of the following two states: s i 5 11, in which vertical (north-south) traffic flow is allowed, and s i 5 21, in which horizontal (east-west) traffic flow is allowed. We assume that when s i 5 11, only the vehicles coming from the vertical lanes flow into the crossing at the rate of 1 (vehicle) per unit time; at the crossing, vehicles go straight at the rate of a, and they turn right or left at the rate of 1 2 a, where a g [0, 1] is a parameter of the model. Similarly, when s i 5 21, vehicles from the horizontal lanes are allowed to enter the crossing, where they behave in the same way.
Then, the number of vehicles q ij in the traffic lane from node j to node i evolves according to the following differential equation: where a stands for 2a 2 1, and s ij g {11, 21} denotes the direction of the lane from node j to node i (11 for vertical and 21 for horizontal). Traffic signals are controlled locally; essentially, they are changed to allow traffic from the lanes in which more vehicles are waiting at each crossing. This type of dynamic and decentralised control strategies, rather than conventional cyclic and centralised ones, has recently been considered as promising [1][2][3][4] . Let us consider the average numbers of vehicles coming from the vertical and horizontal lanes at the ith node. Their difference is given by x i~P j[N i ð Þ s ij q ij 2, where N(i) denotes the set of four neighbours of the ith node. When x i . 0, more vehicles are waiting in the vertical direction, thus the vertical signals should be green (s i 5 11), and when x i , 0, the horizontal signals should be green (s i 5 21). To avoid frequent switching, signals are controlled by the following on-off control rule: where h . 0 is the threshold parameter for the deadband, and we set h 5 1 without loss of generality in the following.
Equations (1) and (2) define the dynamics of the traffic signal model. However, instead of considering the numbers of vehicles q ij in all the lanes, it is sufficient to consider only the differences x i at all the nodes. The time evolution of x i is described by Note that because 21 # a # 1, x i never increases when s i 5 11, and never decreases when s i 5 21. Therefore, x i can be controlled within the interval [2h, 1h]. This also means that the values of q ij are bounded and ensures that the lanes become neither too much congested nor empty provided that there are sufficient vehicles at the initial time. Because the traffic flows are controlled by discrete signals, the model can be regarded as a switched flow system 2,19 . More generally, it is a hybrid dynamical system 20 , because it is a completely deterministic system that consists of both continuous variables {x i } and discrete variables {s i }. The state (x 1 , …, x N ) of N 5 L 3 L crossings moves straight in the N-dimensional hypercube [2h, 1h] N according to equation (3), and its direction changes only on the boundary of the hypercube according to equation (2). Therefore, its dynamics can be understood as a pseudo billiard 21 or a strange billiard 22,23 . Note that the dynamics is invertible. For a 5 1/2 (i.e., a 5 0), interactions between the x i 's vanish, and the dynamics reduces to a normal hypercubic billiard. The decoupled dynamics of the nodes can be considered as nonchaotic oscillators because they behave periodically with a period of 4h.
Ising-like dynamics. Typical dynamics of the model is shown in { P ij s i s j N, respectively, for several values of a. Apparently the system shows probabilistic behaviour, although each element is non-chaotic.
Neighbouring spins have ferromagnetic interactions for a . 0, since as equation (3) shows, x i moves more slowly when more neighbouring spins are aligned. Similarly, they have antiferromagnetic interactions for a , 0, and the model is symmetric with respect to a 5 0 (see Methods for mathematical properties of the model). The parameter a corresponds to the coupling strength between the signals, which controls the apparent temperature of the system; as a goes away from 0, the system is considered to be cooled down. As a approaches 1, the average cluster size grows (Fig. 2a) and the average energy decreases (Fig. 2c).
In the extreme case of a 5 1, the system dynamics has two absorbing states; as equation (3) shows, once all the spins are aligned, the state never changes, and the system freezes. Figure 3 shows the time required to reach one of the absorbing states starting from a random (high-temperature) initial state. In all the calculations for L # 32, the system dynamics reaches the absorbing states in a finite time t # 10 8 . This implies that the traffic signal model of finite size undergoes spontaneous symmetry breaking as it approaches a 5 1. However, the absorption time is likely to diverge in the thermodynamic limit, because in Fig. 3 the absorption time seems to increase exponentially as L increases.
Thermodynamic behaviour of the model as it approaches a 5 1 is shown in Fig. 4. For jaj , 1, no absorbing states exist, and similarly to the case of the Ising model, it can be shown that the long-term average of the magnetisation AEmae is strictly zero (see Methods). As  a approaches 1, the average absolute magnetisation AEjmjae increases from 0 to 1 (Fig. 4a), which suggests a phase transition from the paramagnetic phase to the ferromagnetic phase. In the corresponding range of a, the magnetic susceptibility x 5 N(AEm 2 ae 2 AEjmjae 2 ) and the specific heat c~N 2 h i{h i 2 À Á have peaks ( Fig. 4c and d). The phase transition point indicated by these peaks seems to converge to a 5 1 in the thermodynamic limit. This implies that the ordered behaviour exists only at a 5 1 in the thermodynamic limit. In this sense, this phase transition in the model can be observed only as a finite size effect.
More detailed analysis of larger systems is necessary for accurate estimation of the critical exponents. However, because the calculation of thermodynamic quantities of the model requires lengthy computation, these analyses do not seem tractable at present.
The thermodynamic behaviour is realised by nonlinear dynamics of coupled non-chaotic elements. The dynamics is not ergodic in a rigorous sense, because it can be shown that it has invariance (see Methods). However, this non-ergodicity may be subtle from a macroscopic point of view. The dynamics can be characterised by the piecewise linear Poincaré map defined on the boundary of the hypercube, which we call the spin-flip map (see Methods).
We numerically investigated the dynamics for the smallest case L 5 2. Figure 5 shows the largest Lyapunov exponents of the spin-flip map for initial values on a plane in the boundary. The largest Lyapunov exponents are non-negative. It can be observed that the instability of the orbits depends strongly on the invariant plane; namely, it depends on the choice of initial states. Moreover, even if we limit the dynamics to an invariant plane, there are typically multiple invariant sets (ergodic components). Therefore, even in the chaotic case, the state wanders in a limited region on the invariant plane. It should be investigated as a next step whether these properties for L 5 2 remain in larger systems.

Discussion
Chaotic dynamics is known to emerge even in a very simple switched flow system 19 . Therefore, it has been pointed out that city traffic as a switched flow system may exhibit chaotic behaviour 2 . In this regard, we have presented the first concrete traffic model described by a switched flow system that exhibits chaotic dynamics in this paper.
The probabilistic distribution of spin configurations realised by the chaotic behaviour is unknown and should be investigated in the future. However, it can be intuitively related to the Ising model. Let us consider the behaviour of the ith node in the lattice surrounded by four neighbouring nodes whose states are all fixed. Then according to  Eq. (3), x i moves at the absolute speed of (1{a P j s j 4) when s i 5 11, and (1za P j s j 4) when s i 5 21. Therefore, if we sample s i at a random time, the probability of observing s i 5 11 is given by 1za 2, which is determined by the states of the neighbouring nodes as in the Ising model. This intuitive discussion gives us a probabilistic Ising-like model that may be approximately followed by the determinstic model.
Because the traffic signal model is very simple, it can be extended in many directions and related to studies in various fields. External magnetic fields can easily be implemented by considering different speeds of flows in the horizontal and vertical directions. Extensions to a d-dimensional hyperlattice and other network structures can be considered by modifying equation (3). The traffic signal model can also be regarded as a coupled oscillator system in which, unlike the Kuramoto model 24 , the elements interact only through binarised phases. In the sense that interactions are discrete, the model may be related to neural networks; in fact, it is similar to a neural network model of simplified hysteresis neurons 25 .
Equation (3) can be generalised to the form _ , which in conjunction with equation (2) defines a new class of dynamical systems with a large number of degrees of freedom. This new class is expected to exhibit chaotic pseudo-billiard dynamics that arises from coupled non-chaotic elements. Beyond what is shown in this paper, it may have connections to various models such as spin models in statistical physics, coupled oscillators and CMLs in nonlinear physics, and cellular automata in computer science. Recently, we have shown that the Boltzmann machines can be simulated by chaotic dynamics in this class 26 .
It should be also noted that over the past half-century, ideas from statistical physics have been successfully applied to traffic flows 27 . In particular, cellular automaton models of city traffic such as the BML model 28 (11,11,11,11). Lines in the diagonal direction (x 1 2 x 2 5 const.) correspond to invariant planes. The largest Lyapunov exponents vary depending on the invariant plane. Even on the same invariant plane, non-chaotic (black) and chaotic (coloured) regions coexist. The dynamics depends strongly on the choice of initial values. two-dimensional lattices. However, while these studies focused on vehicles and characterised traffic jams as phase transitions, we focused on traffic signals, which themselves constitute a many-body system in cities. Although the model presented here is too simple for real city traffic, it is natural to assume that signals are interacting with each other if they are somehow controlled locally, as pointed out in the previous studies 1,2 . In this sense, our model is expected to serve as a starting point for considering statistical mechanics of traffic signals.

Methods
Symmetries in the model. To understand the dynamics of the model, it is important to understand symmetries inherent in the model. The time evolution of x i 's (Eq. (3)) can be written in a matrix form as follows: where x~x 1 , . . . , x N ð Þ > , s~s 1 , . . . , s N ð Þ > , and the matrix A is the adjacency matrix of the lattice. There is a trivial symmetry in the state space of the nodes as follows: There is also a symmetry as to time reversal as follows: Let L 1 and L 2 be the sets of the black and white sites in the checkerboard pattern on the lattice. Define which shows the symmetry between ferromagnetism and antiferromagnetism (1a and 2a).
Long-term average of the magnetisation. In the traffic signal model, it can be shown that the long-term average of the magnetisation is strictly zero, as in the Ising model on a finite-size lattice. By averaging the system dynamics from time 0 to T, we have In the limit T R ', the left-hand side of this equation goes to zero, because x is bounded. For jaj , 1, the matrix (2I 1 (a/4)A) is diagonal dominant and regular, so that AEsae is zero. Therefore, the long-term average of the magnetisation AEmae must be also zero.
Invariant hyperplanes. As one of the dynamical properties of the model, it can be shown that the system is not ergodic. Let us define C as follows: While there is no spin flips, the value of C does not change, because where N(i) denotes the set of four neighbouring nodes of node i. When a spin flips, the value of C changes only by multiples of two. When s i changes from 11 to 21, x i has to be 11, and similarly when s i changes from 21 to 11, x i has to be 21. Therefore, s i x i always changes from 11 to 21 when the ith spin flips. Therefore the system evolves within the invariant hyperplanes determined by the invariance of C modulo 2.
Spin-flip map. To analyse the dynamics of the model, it is important to analyse the Poincaré map defined on the boundary of the hypercube, which we call the spin-flip map. Let t 1 ƒt 2 ƒt 3 ƒ Á Á Á be the time sequence of spin flips. For convenience, define x[k] 5 x(t k ) and s[k] 5 s(t k 1 0). The kth time interval t k 5 t k11 2 t k of the spin sequence is given by where e i is the ith standard basis for N-dimensional Euclidean space and Let i[k] be the spin actually flipped at time t k . Then, we have the spin-flip map Therefore, the spin-flip map is piecewise linear.