The Biham-Middleton-Levine traffic model for a single junction

In the Biham-Middleton-Levine traffic model cars are placed with some density p on a two dimensional torus, and move according to a (simple) set of predefined rules. Computer simulations show this system exhibits many interesting phenomena: for low densities the system self organizes such that cars flow freely while for densities higher than some critical density the system gets stuck in an endless traffic jam. However, apart from the simulation results very few properties of the system were proven rigorously to date. We introduce a simplified version of this model in which cars are placed in a single row and column (a junction) and show that similar phenomena of self-organization of the system and phase transition still occur.

The BML traffic model
The Biham-Middleton-Levine (BML) traffic models was first introduced in [3] published 1992. The model involves two types of cars: "red" and "blue". Initially the cars are placed in random with a given density p on the N × N torus. The system dynamics are as follows: at each turn, first all the red cars try to move simultaneously a single step to the right in the torus. Afterwards all blue cars try to move a single step upwards. A car succeeds in moving as long as the relevant space above/beside it (according to whether it is blue/red) is vacant.
The basic properties of this model are described in [3] and some recent more subtle observations due to larger simulations are described in [4]. The main and most interesting property of the system originally observed in simulations is a phase transition: for some critical density p c one observes, that while filling the torus randomly with cars in density p < p c the system self organizes such that after some time all cars flow freely with no car ever running into another car (see figure 1), by slightly changing the density to some p > p c not only does system not reach a free flow, but it will eventually get stuck in a configuration in which no car can ever move (see figure 2).
Very little of the above behaviour is rigorously proven for the BML model. The main rigorous result is that of Angel, Holroyd and Martin [2], showing that for some Figure 1: Self organization for p = 0.3: the initial configuration in the left organizes after 520 steps to the "diagonal" configuration on the right which flows freely Figure 2: Traffic jam for p = 0.35: the initial configuration on the left gets to the stuck configuration on the right after 400 steps, in which no car can ever move fixed density p < 1, very close to 1, the probability of the system getting stuck tends to 1 as N → ∞. In [1] one can find a study of the very low density regime (when p = O( 1 N )).
First, we introduce a slight variant of the original BML model by allowing a car (say, red) to move not only if there is a vacant place right next to it but also if there is a red car next to it that moves. Thus, for sequence of red cars placed in a row with a single vacant place to its right -all cars will move together (as oppose to only the rightmost car in the sequence for the original BML model). Not only does this new variant exhibits the same phenomena of self-organization and phase transition, they even seem to appear more quickly (i.e. it takes less time for the system to reach a stable state). Actually the demonstrated simulations (figures 1, 2) were performed using the variant model. Note that the results of [2] appear to apply equally well to the variant model.
In the following, we will analyze a simplified version of the BML model: BML on a single junction. Meaning, we place red cars in some density p on a single row of the torus, and blue cars are placed in density p on a single column. We will show that for all p the system reaches optimal speed 1 , depending only on p. For p < 0.5 we will show the system reaches speed 1, while for p > 0.5 the speed cannot be 1, but the system will reach the same speed, regardless of the initial configuration. Moreover, at p = 0.5 the system's behaviour undergoes a phase transition: we will prove that while for p < 0.5 the stable configuration will have linearly many sequences of cars, for p > 0.5 we will have only O(1) different sequences after some time. We will also examine what happens at a small window around p = 0.5.
Note that in the variant BML model (and unlike the original BML model) car sequences are never split. Therefore, the simplified version of the variant BML model can be viewed as some kind of 1-dimensional coalescent process.
Much of the proofs below rely on the fact that we model BML on a symmetric torus. Indeed, the time-normalization of section 2.2 would take an entirely different form if the height and width were not equal. We suspect that the model would exhibit similar properties if the height and width had a large common denominator, e.g. if they had a fixed proportion. More importantly, we believe that for height and width with low common denominator (e.g. relatively prime), we would see a clearly different behaviour. As a simple example, note that in the relatively prime case, a speed of precisely 1 cannot be attained, no matter how low is p, in contrast with corollary 3.14. This dependence on the arithmetic properties of the dimensions is also apparent in [4].

Basic model
We start with the exact definition of our simplified model. On a cross shape, containing a single horizontal segment and a single vertical segment, both of identical length N , red cars are placed in exactly pN randomly (and uniformly) chosen locations along the row, and blue cars are similarly placed in pN locations along the column. p will be called the density of the configuration. For simplicity and consistency with the general BML model, we refer to the cars placed on the horizontal segment as red and those on the vertical segment as blue. For simplicity we may assume that the junction, i.e. the single location in the intersection of the segments, is left unoccupied. The segments are cyclic -e.g. a red car located at the rightmost position (N − 1) that moves one step to the right re-emerges at the leftmost position (0).
At each turn, all the red cars move one step to the right, except that the red car that is just left of the junction will not move if a blue car is in the junction (i.e. blocking it), in which case also all the red cars immediately following it will stay still. Afterwards, the blue cars move similarly, with the red and blue roles switched. As in the original BML, we look at the asymptotic speed of the system, i.e. the Figure 3: On the left: a junction configuration and the analogous configuration beneath it. The junction is marked with an asterisk. On the right: same configuration after 3 turns in both views (asymptotic) average number of steps a car moves in one turn. It is easily seen that this speed is the same for all blue cars and all red cars, since the number of steps any two cars of the same color have taken cannot differ by more than N . It is somewhat surprising, perhaps, that these two speed must also be the same. For instance, there is no configuration for which a blue car completes 2 cycles for every 1 a red car completes.

Time-normalized model
Though less natural, it will be sometimes useful to consider the equivalent description, in which two rows of cars -red and blue -are placed one beneath the other, with a "special place" -the junction -where at most one car can be present at any time. In every step first the red line shifts one to the right (except cars immediately to the left of the junction, if it contains a blue car) and then the blue line does the same. Furthermore, instead of having the cars move to the right, we can have the junction move to the left, and when a blue car is in the junction, the (possibly empty) sequence of red cars immediately to the left of the junction moves to the left, and vice verse. Figure 3 illustrates the correspondence between these models.
From the discussion above we get the following equivalent system, which we will call the time-normalized junction: 1} N representing the red and blue cars respectively, i.e. R i = 1 iff there's a red car at the i-th place. We require that R i = B i = p · N , and at place S (= N − 1) there is at most one car in both rows 2 .
2. In each turn: 2 The last requirement is that the junction itself can contain only one car at the beginning • If place S contains a blue car, and place S − 1 contains a red car (if B S = R S−1 = 1), push this car one step to the left. By pushing the car we mean also moving all red cars that immediately followed it one step to the left, i.e. set R S−1 = 0, R S−i = 1 for i = min j≥1 [R S−j = 0]. • set S = S − 1 Note that indeed the dynamics above guarantee that after a turn place (S − 1) -the new junction -contains at most one car. Generally, as long as cars flow freely, the time-normalized system configuration does not change (except for the location of the junction). Cars in the time-normalized system configuration actually move only when some cars are waiting for the junction to clear (in the non-time-normalized system).

Analysis of an (N, p) junction
Our analysis of the junction will argue that for any p, regardless of the initial configuration, the system will reach some optimal speed, depending only on p. First let us state what is this optimal speed: Theorem 3.1 For a junction with density p, the maximal speed for any configuration is min(1, 1 2p ).
Proof Obviously, the speed cannot exceed 1. Observe the system when it reaches its stable state, and denote the speed by s. Thus, at time t, a car has advanced ts(1+ o(ts)) steps (on average) which means that it has passed the junction ts/N (1 + o(1)) times. As only one car can pass the junction at any time, the total number of cars passing the junction at time t is bounded by t. Therefore, we have 2pN × ts/N ≤ t which implies s ≤ 1/2p.
We will now show that the system necessarily reaches these speeds from any starting configuration.

The case p < 0.5
We begin by proving that for p < 0.5 the junction will eventually reach speed very close to 1. A system in a stable state is said to be free-flowing if no car is ever waiting to enter the junction. Proof Obviously, this is just a reformulation of free-flowing.
We will now turn to show that for p < 0.5, following the system dynamics of the time-normalized junction will necessarily bring us to a free flowing state, or at worst an "almost" free flowing state, meaning a state for which the system speed will be arbitrarily close to 1 for large enough N .
For this let us consider some configuration and look at the set of "violations" to lemma 3.2, i.e. places that either contain both a blue and a red car, or that contain a blue car and a red car is in the place one to the left. As the following lemma will show, the size of the set of violations is very closely related to the system speed, and posses a very important property: it is non-increasing in time.
More formally, For a configuration R, B we define two disjoint sets for the two types of violations: . For a junction with some initial configuration R, B, let R t , B t be the system configuration at time t, and let V t = V B t ∪ V R t be the set of violations for this configuration, and X t be the corresponding indicator vector. Similarly, let S t denote the junction's position at time t.
For any t, the system speed is at least Proof Property (1) follows from the system dynamics. To see this, examine the three possible cases for what happens at time t: 1. If in turn t place S t does not contain a violation then the configurations do not change during the next turn, i.e.
2. If S t ∈ V B t , then the configuration B t+1 changes in two places: . This may affect X t+1 (S t − i), and X t+1 (S t − i + 1). However, for place S t − i + 1, by changing B t S t −i from 0 to 1 no new violation can be created (since by definition of i, B t+1 For other indices X t+1 (i) = X t (i) since R, B do not change, so between times t and t + 1, we have that X t (S t ) changes from 1 to 0, and at worst only one other place -X t (S t − i + 1) changes from 0 to 1, so 3. Similarly, if place S t ∈ V R t then the configuration R t+1 changes in two places: . However for place S t −i changing R S t −i does not affect whether this place is a violation or not, so at worst X t (S t −i+1) changed from 0 to 1.
By the same argument we get |V t+1 | ≤ |V t |.
For property (2) we note that in the time-normalized system, following a specific car in the system, its "current speed" is N N +k = (1 + k N ) −1 where k is the number of times the car was pushed to the left during the last N system turns. We note that if a car at place j is pushed to the left at some time t, by some violation at place S t , this violation can reappear only to the left of j, so it can push the car again only after S t passes j. Hence any violation can push a car to the left only once in a car's cycle (i.e. N moves). Since by (1) at any time from t onwards the number of violations in the system is at most |V t |, then each car is pushed left only |V t | times in N turns, so its speed from time t onwards is at least ( With lemma 3.3 at hand we are now ready to prove system self organization for p < 0.5. We will show that for p < 0.5, after 2N system turns |V t | = O(1), and hence deduce by part (2) of lemma 3.3 the system reaches speed As the junction advances to the left it pushes some car sequences, thus affecting the configuration to its left. The next lemma will show that when p < 0.5, for some T < N , the number of cars affected to the left of the junction is only a constant, independent of N . Lemma 3.4 Consider a junction with density p < 0.5. There exists some constant C = C(p) = p 1−2p , independent of N , for which: From any configuration R, B with junction at place S there exist some 0 < T < N such that after T turns: Proof First let us consider T = 1. For T = 1 to not satisfy the lemma conclusions, there need to be a car sequence (either red or blue) of length exceeding C, which is pushed left by the junction as it moves. Progressing through the process, if, for some T < N − C, the sequence currently pushed by the junction is shorter then C, then this is the T we seek. Therefore, for the conclusions not to hold, the length of the car sequence pushed by the junction must exceed C for all 0 < T < N . If this is the case, then leaving the junction we see alternating red and blue sequences, all of lengths exceeding C and one vacant place after any blue sequence and before any red one.
However, if this is the case for all 0 < T ≤ T then the average number of cars per location in {S − T , . . . , S} at time T must be at least 2C 2C+1 (at least 2C cars between vacant places). Therefore, the total number of cars in {S − T , . . . , S} at time T is more than 2C 2C+1 T = 2pT (Recall that C = 2p 1−2p ). Since there are only 2pN cars in the system, this cannot hold for all T up to N . Thus, there must be some time for which the conclusions of the lemma are satisfied.
We are now ready to easily prove the main result for this section.
Theorem 3.5 A junction of size N with density p < 0.5 reaches speed of 1 − C(p) N from any initial configuration.
Proof Let R, B be some initial configuration, with S = N − 1 and let V the corresponding set of violations and X the matching indicators vector. By lemma 3.4 there exist T 0 > 0 for which , and also, as long as Proceeding in this manner until T = T i ≥ N we will get that after T turns, X T (i) = 0 for all but at most C(p) places, hence by lemma 3.3 the system speed from this time onward is at least 1 − C(p) N .
We remark one cannot prove an actual speed 1 (rather than 1 − o(1)) for the case p < 1 2 since this is definitely not the case. Figure 4 (a) demonstrates a general construction of a junction with speed 1 − 1 N for all N . Finally we remark that a sharper constant C (p) can be provided, that will also meet a provable lower bound for the speed (by constructions similar to the above). This C (p) is exactly half of C(p). We will see proof of this later in this paper as we obtain theorem 3.5 as a corollary of a different theorem, using a different approach.
3.2 Number of segments for p < 0.5 The proof in the previous section were combinatorial in nature and showed that for a density p < 0.5 the system can slow-down only by some constant number independent of N regardless of the configuration. As we turn to examine the properties of the stable configuration of course we can not say much if we allow any initial configuration. There are configurations in which the number of different segments in However, we will show that for a random initial configuration, the system will have linearly many different segments of cars with high probability. Theorem 3.6 A junction of size N with density p < 0.5, started from a random initial configuration, will have Θ(N ) different segments at all times with high probability (w.r.t. N ).
Proof As we already seen in the proof of 3.4, as the system completes a full round, since every place (but a constant number of places) contains at most a single car there must be (1 − 2p)N places in which no car is present. Each two such places that are not adjacent must correspond to a segment in the cars configuration.
It is evident by the system dynamics, that the number of places for which R i = B i = R i−1 = B i−1 = 0 is non-increasing. More precisely, only places for which R i = B i = R i−1 = B i−1 = 0 in the initial configuration can satisfy this in the future. In a random initial configuration with density p, the initial number of these places is expected to be (1 − p) 4 N , and by standard CLT, we get that with high probability this number is at most ((1 − p) 4 + ε)N . Thus, the number of different segments in the system configuration at any time is at least ((1 − 2p) − (1 − p) 4 − ε)N . However for p very close to 1 2 this bound may be negative, so this does not suffice.
To solve this we note that similarly also for any fixed K, the number of consecutive K empty places in a configuration is non-increasing by the system dynamics, and w.h.p. is at most ((1 − p) 2K + ε)N for an initial random state. But this guarantees at least ( 3.3 The case p > 0.5 The proofs of speed optimality and segment structure for p > 0.5 will rely mainly of a the combinatorial properties of a stable configuration. A stable configuration for the system is a configuration that re-appears after running the system for some M turns. Since for a fixed N the number of possible configurations of the system is finite, and the state-transitions (traffic rules) are time independent, the system will necessarily reach a stable configuration at some time regardless of the starting configuration.
We will use mainly two simple facts that must hold after a system reached a stable configurations: (a) |V t | cannot change -i.e. no violation can disappear. This is clear from lemma 3.3; (b) Two disjoint segments of car cannot merge to one (i.e. one pushed until it meets the other), since clearly the number of segments in the system is also non-increasing in time.
These two facts alone already provide plenty of information about the stable configuration for p > 0.5. We begin with the following twin lemmas on the stable state.
Proof Since it is easy to lose the idea in all the notations, a visual sketch of the proof is provided in figure 5.
By the assumptions: To get (1), we note that with the assumption B 0 = 0, the blue sequence will be pushed to the left in the next s R turns, so by restriction (b), B i = 0 for i ∈ [N − s B − s R − 1, N − s B − 1] since otherwise 2 disjoint blue segments will merge while the sequence is pushed. Thus following the system, after s R turns we will get:  If we follow the system for s B more steps we note that any red car in [N − max(s R , s B ) − 1, N − s R − 1] will be pushed left until eventually hitting the red car already proven to be present at N − max(s R , s B ) − 2], thus The following lemma is completely analogous when reversing the roles of R, B, and can be proven the same way.  Putting lemmas 3.7, 3.8 together we get the following characterization for stable configurations: Proof Figure 6 outlines the proof, without the risk of getting lost in the indices. For the full proof, first get from lemma 3.7 that: From (1)

From (2),(3) we get
contain a sequence of length min(s R , s B ) of red cars with no blue cars in parallel to it. This sequence is possibly a part of a larger sequence of length . Now running the system for max(s R , s B ) turns, we will have the junction at place S = N − max(s R , s B ) − 1, B S = 1, followed by sequences of s R reds and s B = min(s B , s R )(≤ s R ) blues. Applying lemma 3.8 for the system (rotated by max(s R , s B ) i.e. for N = N − max(s R , s B ) − 1): In particular from (3) we get that no blue cars are in parallel to the entire red Furthermore by (2) we get that a sequence of blue cars which begins from place N − max(s R , s B ) − s R − 2 with no red cars in parallel to it by (1). Note that N − max(s R , s B ) − s R − 2 is exactly to the left of N − max(s R , s B ) − s R − 1 where the red sequence ended. Now clearly choosing K 1 = s R ; K 2 = min(s R , s B ) we get claims 3,4,5 in the lemma.
Putting it all together we can now get a very good description of a stable state. 1. Each place contains at most one type of car, red or blue.
2. Place N −M −1 is empty. Each empty place, is followed by a sequence of places containing red cars immediately left to it, which is followed by a sequence of places containing blue cars immediately left to it.
3. Any sequence of red or blue cars is of length at least min(s R , s B ) Proof This is merely applying lemma 3.9 repeatedly. Applying lemma 3.9 we know that there exist K 1 , K 2 ≥ min(s R , s B ) such that: Place N − M − 1 is empty, followed by K 1 consecutive places with only red cars and K 2 consecutive places with only blue cars left to it, thus the assertion holds for the segment [T, , and this is enough to advance the system for K 1 + M turns. The s B blue segment is pushed left s R places, further pushing the K 1 red sequence min(s R , s B ) places to the left such that its last min(s R , s B ) cars now overlap with the K 2 blue sequence.
So after K 1 +M turns the system evolves to a state where S = N −K 1 −M , B S = 0, and left to S there are K 2 > K 2 consecutive blue cars and exactly min(s R , s B ) consecutive red cars. Noting that this time M = max(K 2 , min(s R , s B )) = K 2 , once again we can deduce from lemma 3.9 that: there are no additional cars in [N −M −K 1 −K 2 , N −M −K 1 −1] (thus we are assured that the entire blue segment of length K 2 does not have red cars parallel to it), Place N − M − K 1 − K 2 − 1 is empty, followed by some K 3 consecutive places with only red cars and K 4 consecutive places with only blue cars left to it, for K 3 , K 4 ≥ min(K 1 , K 2 ) = min(s R , s B ) thus the assertion holds for the segment [T , Repeatedly applying lemma 3.9 as long as T > 0, we repeatedly get that the assertion holds for some [T, N − M − 1] for T strictly decreasing, so the assertion holds in [0, N − M − 1] We have worked hard for theorem 3.10, but it will soon turn to be worthwhile. Let us obtain some useful corollaries. On the other hand, Since any sequence of a red car or a blue car is of length at least m, an empty place in [0, N − m] can occur only once in 2m + 1 places, and other places contain one car. Thus the number of cars is lower bounded by The last property follows the fact that any sequence of cars is of length at least m, and by (3) total number of cars is at most N + m, thus number of different sequences is at most N +m m = N m + 1. Theorem 3.12 All cars in the system have the same asymptotic speed.
Proof As we have seen, when it reaches a stable state, the system consists of alternating red and blue sequences of cars. Obviously, the order of the sequences cannot change. Therefore, the difference between the number of steps two different cars have taken cannot be more then the length of the longest sequence, which is less then N . Thus, the asymptotic (w.r.t. t) speed is the same for all cars.
With these corollaries we can now completely characterize the stable state of a junction with p > 0.5, just by adding the final simple observation, that since the number of cars in the model is greater than N , there are violations at all times, including after reaching a stable state. Now let us look at some time when the junction reaches a violation when the system is in stable state. At this point the conditions of theorem 3.10 are satisfied, thus: Theorem 3.13 A junction of size N and density p > 0.5 reaches speed of 1 2p −O( 1 N ) (i.e. arbitrarily close to the optimal speed of 1 2p , for large enough N ), and contains at most a bounded number (depending only on p) of car sequences.
Proof We look at the system after it reached a stable state. Since 2pN > N at some time after that conditions of theorem 3.10 are satisfied for some s R , s B ≥ 1. Let m = min(s R , s B ) at this time. Using claims (3),(4) in corollary 3.11 we get: and reusing m ≥ (2p − 1)N on the left hand size we get: (1)). So m = (2p − 1)N + K, for K ≤ C. Now by claim (2) in 3.11, system speed is at least But by theorem 3.1 system speed is at most 1 2p , thus system speed is exactly 1 proving the first part of the theorem.
By (5) we get total number of segments in the system is at most N m + 1, applying m ≥ (2p − 1)N we get the number of segments is bounded by: Thus the second part proven.

p < 0.5 revisited
The characterization in 3.10 can be also proven useful to handle p < 0.5. Actually the main result for p < 0.5, theorem 3.5, can be shown using similar technique, and even sharpened. 3 Corollary 3.14 For p < 0.5 the junction reaches speed of at least 1 − C(p) N , for C(p) = p 1−2p . In particular, for p < 1 3 the junction reaches speed 1, for any initial configuration.
Proof Let R, B be any initial configuration. Looking at the configuration after it reached the stable state, if the system reached speed 1 we have nothing to prove. Assume the speed is less than 1. Since in this case violations still occur, at some time the stable configuration will satisfy theorem 3.10. As before, letting m = min(s R , s B ) at this time, by claim (4) in corollary 3.11 we have: In particular 2p − 1 + 1 2m+1 > 0, rearranging we get m < p 1−2p , and since m ∈ Z, m ≤ p 1−2p = C(p). m must be positive, thus for p < 1 3 , having C(p) = p 1−2p = 0 we get a contradiction, thus the assumption (speed < 1) cannot hold. For p ≥ 1 3 by claim (2) in corollary 3.11 the system speed is at least ( 3.5 The critical p = 0.5 Gathering the results so far we get an almost-complete description of the behaviour of a junction. For junction of size N and density p: • If p < 0.5 the junction will reach speed 1 − o(1) (asymptotically optimal), and contain linearly many different segments in the stable state.
• If p > 0.5 the junction will reach speed 1 2p − o(1) (asymptotically optimal), and contain constant many segments in the stable state.
From the description above one sees that the junction system goes through a sharp phase transition at p = 0.5, as the number of segments of cars as the system stabilizes drops from being linear to merely constant. The last curiosity, is what happens at p = 0.5. Once again by using the powerful theorem 3.10 we can deduce: Proof For p = 0.5 we have exactly N cars in the system. As we reach stable state, violations must still occur (since a system with exactly N cars must contain at least one violation), thus at some time theorem 3.10 is satisfied. For m = min(s R , s B ) at this time, by claim (4) in corollary 3.11 we have: From here by claim (2) the system speed is at least (1 + If S is the number of segments, then by theorem 3. Simulation results show that these bounds are not only tight, but typical, meaning that a junction with a random initial configuration with density p indeed has O( √ N ) segments in the stable state, with the largest segment of size near n 1/2 . This suggests that the system undergoes a second order phase transition.

Simulation results
Following are computer simulation results for the junction for critical and near critical p, demonstrating the phase transition. The columns below consist of N , the average asymptotic speed, the average number of car segments in the stable state, the average longest segment in the stable state, and the average number of segments divided by N and √ N .

p = 0.48
For large N , system reaches speed 1 and the average number of segments is linear (approx. 0.037N ).

Summary
The fascinating phenomena observed in the BML traffic model are still far from being completely understood. In this paper we showed a very simplified version of this model, which, despite its relative simplicity, displayed very similar phenomena of phase transition at some critical density and of self-organization, which in our case both can be proven and well understood.
We used two approaches in this paper: The first one was to use some sort of "independence" in the way the system evolves, as in the way we handle p < 0.5. We showed that the system self-organizes "locally" and with a bounded affect on the rest of the configuration, thus it will eventually organize globally. The second approach is the notion of the stable configuration, i.e. we characterize the combinatorial structure of any state that the system can "preserve" and use it to show it is optimal (in a way saying, that as long as we are not optimal the system must continue to evolve).
Can these results be extended to handle more complicated settings than the junction? Possibly yes. For example, considering a k-lanes junction (i.e. k consecutive red lines meeting k consecutive blue rows), one can look at a time-space-normalized version of the system, as shown in figure 9, with the junction now being a k × k parallelogram traveling along the red and blue lines, and "propagating violations" (which now have a more complicated structure depending on the k × k configuration within the junction). Stable states of this configurations seem to have the same characteristics of a single junction, with a red (blue) car equivalent to some red (blue) car in any of the k red (blue) lines. Thus a zero-effort corollary for a k-lanes junction is that it reaches speed 1 for p < 1 2k , but for k nearing O(N ) this bound is clearly non-significant. It is not surprising though, since combinatorics alone cannot bring us too far, at least for the complete BML model -even as little as 2N cars in an N 2 size torus can be put in a stuck configuration -i.e. reach speed 0.


Why is a new Journal of Informetrics needed?

In our study we analysed 3.889 records which were indexed in the Library and Information Science Abstracts (LISA) database in the research field of informetrics. We can show the core journals of the field via a Bradford (power law) distribution and corroborate on the basis of the restricted LISA data set that it was the appropriate time to found a new specialized journal dedicated to informetrics. According to Bradford's Law of scattering (pure quantitative calculation), Egghe's Journal of Informetrics (JOI) first issue to appear in 2007, comes most probable at the right time.

Introduction
In the new Journal of Informetrics (JOI 2 ) recent call for papers Leo Egghe mentioned that Elsevier (Oxford, UK) accepted his proposal for the foundation of this new scientific journal. It will be the first journal worldwide which carries 'informetrics' in the journal title. He explained in an editorial why this is necessary (Egghe, 2005). One main reason is the growth of the field informetrics. A broad definition of the field informetrics can be found in Egghe's introductory paper in the first special issue on informetrics in Information Processing and Management Vol. 41 (6): "we will use the term "informetrics" as the broad term comprising all-metrics studies related to information science, including bibliometrics (bibliographies, libraries, ...), scientometrics (science policy, citation analysis, research evaluation, ...), webometrics (metrics of the web, the Internet or other social networks such as citation or collaboration networks), ..." (Egghe, 2005(Egghe, p. 1311).
Egghe states several studies to show that the number of papers and authors in informetrics are constantly growing. In terms of multidisciplinarity, the field informetrics is also expanding, this can be shown by new research topics like webometrics and other web based informetric methods (see a review in Björneborn & Ingwersen, 2004). In the following sections we will show briefly that the distribution of papers published in the broad field of informetrics can be used to consider if a new forum of scientific communication is necessary. For establishing a new scholarly journal publishers generally require, on average, not only hundred papers per year from a hundred scientists (Umstätter, 2003), they also need to have a few hundred interested readers.

Informetric papers in the LISA database
To generate a data set of papers in the field of informetrics we used the Library and Information Science Abstracts (LISA) database 3 . The query bibliomet* OR cybermet* OR infomet* OR informet* OR scientomet* OR webomet* generated a total of 3.889 records that are indexed in LISA (1976-2004. The records are indexed under the following terms (  First we produced a bradfordized list to show the core journals which published the most papers in the field (see Figure 1). By using Bradford's Law of scattering (see e.g. Garfield, 1980;White, 1981;Nicolaisen & Hjørland, 2007 to appear) it is possible to extrapolate the number of publications in the world if we know the core journals. In principle, if a new journal arises the number of scientific publications, has to be about 2 to 3 times higher than 100, because Bradford's Law of scattering means, that all other journals in the world are also attractive for such publications, but with a decrease that follows the power law. It is simple to see that the core journal of the field Scientometrics 5 is ~2 to 4 times higher than expected by the Bradford distribution (see Figure 1 and Table  2). The core journals in our Bradford analysis are the following (compare to the journals Egghe mentioned in his editorial (2005)).  It seems to be clear that some journals are highly represented due to a special indexing of LISA. Beyond that we extracted the most productive authors in the informetrics field (see Table 3).

Growth and competition in informetrics
Starting with "statistical bibliography", in 1922 we can understand E. Windham Hulm as a pioneer. Interesting works from Bradford (1934), Lotka, Otlet, Pritchard, Zipf and others followed. If we look to the expansion of this field in LISA the growth rate seems to be some 10 years (see Figure 2). A rough estimation of papers in the field of informetrics in LISA shows the doubling rate for a 10 years period to be twice the normally seen value of 20 years. Under the assumption that LISA (dotted line in Figure 2), like similar databases, has only ~40% of the whole amount, we can calculate a number of ~500 published papers in 2002 (straight line in Figure 2).

Growth of publications in "informetrics"
Calculating the core journals only for the year 2002 (without the comparable too high value of Scientometrics), then we can extrapolate it to the total number of some 100.000 running journal titles in the world (see Figure 3). In this case we get 430 papers in 2002. By a growth rate of 7% per year (t 2 = 10 years), we have approximately 600 to 700 papers in 2006. But this value is also based on LISA, and additionally we know that some of the core journals are not indexed adequately. Another possibility is to calculate Bradford's Law of scattering using as of all the years included in LISA (see Table 4 as an idealized example) to reconstruct the publication growth for all journals with a doubling time of 9.6 years.  Then we can estimate roughly the number of 1.000 publications in 2005 in this field (see Table 5). The factor 2.5 means, that empirical assessments have shown repeatedly, that the search in one database has only a recall ratio of roughly 40%. Searchers know the motto, since the seventies of the last century, "one database is not enough".  A comparison of editorial boards is warranted as theses scientists often use their journals as a platform for their own publications (see Nourmohammadi & Umstätter, 2004) and collaborations 6 . Regardless, new journals are often founded too early. The consequence can be observed very easily by the death of roughly 50% of new journals (Umstätter & Rehm, 1984). The reason for this is also very simple; scientists on the search front always want be the first.

Conclusion
There has to be improved standards of quality and increasing specialization has to be found in the growing area of informetrics. Better measurements of information and knowledge would be one of the great desiderata. If we remember that in 1963 a first very rough estimation of the amount of information in the Library of Congress was calculated with 10 13 bit 7 .
The new Journal of Informetrics with its focus on "good mathematical (probabilistic) models and explanations of informetric regularities (in the broad sense) and/or papers in which interesting and important data-gathering" comes at the right time, according to our study. The basis of a well known publisher is empirically also very important.


extfieldsofstudy:
Computer Science
Mathematics

text:
Sample size cognizant detection of signals in white noise

The detection and estimation of signals in noisy, limited data is a problem of interest to many scientific and engineering communities. We present a computationally simple, sample eigenvalue based procedure for estimating the number of high-dimensional signals in white noise when there are relatively few samples. We highlight a fundamental asymptotic limit of sample eigenvalue based detection of weak high-dimensional signals from a limited sample size and discuss its implication for the detection of two closely spaced signals. This motivates our heuristic definition of the 'effective number of identifiable signals.' Numerical simulations are used to demonstrate the consistency of the algorithm with respect to the effective number of signals and the superior performance of the algorithm with respect to Wax and Kailath's"asymptotically consistent"MDL based estimator.

INTRODUCTION
The observation vector, in many signal processing applications, can be modelled as a superposition of a finite number of signals embedded in additive noise. Detecting the number of signals present becomes a key issue and is often the starting point for the signal parameter estimation problem. When the signals and the noise are assumed to be samples of a stationary, ergodic Gaussian vector process, the sample covariance matrix formed from m observations has the Wishart distribution.
The proposed algorithm uses an information theoretic criterion, motivated by the approach taken by Wax and Kailath (henceforth WK) in [1], for determining the number of signals in white noise by performing inference * Thanks NSF DMS-0411962 and ONR N00014-07-1-0269. on the eigenvalues of the resulting sample covariance matrix. The form of the estimator is motivated by the distributional properties of moments of the eigenvalues of large dimensional Wishart matrices [2].
The proposed estimator was derived by explicitly accounting for the blurring and fluctuations of the eigenvalues due to sample size constraints. Consequently, there is a greater theoretical justification for employing the proposed estimator in sample starved settings unlike the WK estimators which were derived assuming that the sample size greatly exceeds the number of sensors. This is reflected in the improved performance relative to the "asymptotically consistent" WK MDL based estimator.
Another important contribution of this paper is the description of a fundamental limit of eigen-inference, i.e., inference using the sample eigenvalues alone. The concept of effective number of identifiable signals, introduced herein, explains why, asymptotically, if the signal level is below a threshold that depends on the noise variance, sample size and the dimensionality of the system, then reliable detection is not possible. This paper is organized as follows. The problem is formulated in Section 2. An estimator for the number of signals present that exploits results from random matrix theory is derived in Section 3. The fundamental limits of sample eigenvalue based detection and the concept of effective number of signals are discussed in Section 4. Simulation results are presented in Section 5 while some concluding remarks and directions for future research are presented in Section 6.

PROBLEM FORMULATION
We observe m samples ("snapshots") of possibly signal bearing n-dimensional snapshot vectors x 1 , . . . , x m where for each i, x i ∼ N n (0, R) and x i are mutually independent. The snapshot vectors are modelled as where z i ∼ N n (0, σ 2 I), denotes an n-dimensional (real or complex) Gaussian noise vector where σ 2 is generically unknown, s i ∼ N k (0, R s ) denotes a k-dimensional (real or complex) Gaussian signal vector with covariance R s , and A is a n × k unknown non-random matrix.
Since the signal and noise vectors are independent of each other, the covariance matrix of x i can hence be decomposed as where with ′ denoting the conjugate transpose. Assuming that the matrix A is of full column rank, i.e., the columns of A are linearly independent, and that the covariance matrix of the signals R s is nonsingular, it follows that the rank of Ψ is k. Equivalently, the n − k smallest eigenvalues of Ψ are equal to zero. If we denote the eigenvalues of R by λ 1 ≥ λ 2 ≥ . . . ≥ λ n then it follows that the smallest n − k eigenvalues of R are all equal to σ 2 so that Thus, if the true covariance matrix R were known apriori, the dimension of the signal vector k can be determined from the multiplicity of the smallest eigenvalue of R. The problem in practice is that the covariance matrix R is unknown so that such a straight-forward algorithm cannot be used. The signal detection and estimation problem is hence posed in terms of an inference problem on m samples of n-dimensional multivariate real or complex Gaussian snapshot vectors.
A classical approach to this problem, developed by Bartlett [3] and Lawley [4], uses a sequence of hypothesis tests. Though this approach is sophisticated, the main problem is the subjective judgement needed by the practitioner in selecting the threshold levels for the different tests. This was overcome by Wax and Kailath in [1] wherein they propose an estimator for the number of signals (assuming m > n) based on the eigenvalues is the matrix of observations (samples). The Akaike Information Criteria (AIC) form of the estimator is given bŷ while the Minimum Descriptive Length (MDL) criterion is given bŷ is the geometric mean of the n − k smallest sample eigenvalues and a(k) = 1 n−k n j=k+1 l j is their arithmetic mean. It is known [1] that the AIC form inconsistently estimates the number of signals, while the MDL form estimates the number of signals consistently. The simplicity of the estimator, and the large sample consistency are among the primary reasons why the Kailath-Wax MDL estimator continues to be employed in practice [5]. In the two decades since the publication of the WK paper, researchers have come up with many innovative solutions ( [6,7,8] to list a few) for making the estimators more robust by exploiting some type of prior knowledge.
The most important deficiency of the WK and related estimators that remains unresolved occurs when the sample size is smaller than the number of sensors, i.e., when m < n. In this situation, the SCM is singular and the estimators become degenerate. Practitioners often overcome this in an ad-hoc fashion by, for example, restricting k in (7) to integer values in the range 0 ≤ k < min(n, m). Since large sample, i.e., m ≫ n, asymptotics were used to derive the estimators in [1], there is no rigorous theoretical justification for such a reformulation even if the simulation results suggest that the WK estimators are working "well enough." Other sample eigenvalue based solutions found in the literature that exploit the sample eigenvalue order statistics [6], or employ a Bayesian framework by imposing priors on the number of signals [9] are computationally more intensive and do not address the sample size starved setting in their analysis or their simulations. Particle filter based techniques [8], while useful, require the practitioner to the model the eigenvectors of the underlying population covariance matrix as well; this makes them especially sensitive to model mismatch errors that are endemic to high-dimensional settings. This motivates our development of an sample eigenvalue based estimator with a computational complexity comparable to that of the WK estimators.

ESTIMATING THE NUMBER OF SIGNALS
Given an observation y = [y(1), . . . , y(N )] and a family of models, or equivalently a parameterized family of probability densities f (y|θ) indexed by the parameter vector θ, we select the model which gives the minimum Akaike Information Criterion (AIC) [10] defined by where θ is the maximum likelihood estimate of θ, and k is the number of free parameters in θ. We derive an AIC based estimator for the number of signals by exploiting the following distributional properties of the moments of eigenvalues of the (signal-free) SCM. (1) with k = 0, λ = 1 then as m, n → ∞ and c m = n/m → c ∈ (0, ∞), then 

Theorem 1 (Dumitriu-Edelman [2]) Assume R is formed from m snapshots modelled as
where D denotes convergence in distribution, β = 1 (or 2) when x i is real (or complex) valued, respectively, and Q = 2 β c 2c (c + 1) 2c (c + 1) 2c (2c 2 + 5c + 2) . PROOF. This follows from applying the delta method [11] to the results in Theorem 1. When k > 0 signals are present and assuming k ≪ n, then the distributional properties of the n − k "noise" eigenvalues are closely approximated by the distributional properties of the eigenvalues given by Theorem 1 of the signal-free SCM, i.e., k = 0. Hence, by evaluating the statistic in Proposition 2 over a sliding window, and using the normal approximation for the statistic from Proposition 2 with c ≈ n/m and k + 1 free parameters in the AIC formulation in (8) where Here β = 1 if x i ∈ R n , and β = 2 if x i ∈ C n . When the measurement vectors represent quaternion valued narrowband signals, then we set β = 4. Quaternion valued vectors arise when the data collected from vector sensors is represented using quaternions as in [12].

FUNDAMENTAL LIMIT OF DETECTION
The following result exposes when the "signal" eigenvalues are asymptotically distinguishable from the "noise" eigenvalues.
PROOF. This result appears in [13] for very general settings. A matrix theoretic proof for when c < 1 for the real case may be found in [14] and an interacting particle system interpretation appears in [15].
Motivated by Proposition 3, we define the effective number of signals as k eff (R) = # eigs. of R > σ 2 1 + n m .

Identifiability of closely spaced signals
Suppose there are two uncorrelated (hence, independent) signals so that R s = diag(σ 2 S1 , σ 2 S2 ). In (1), let A = [v 1 v 2 ]. In a sensor array processing application, we think of v 1 ≡ v(θ 1 ) and v 2 ≡ v 2 (θ 2 ) as encoding the array manifold vectors for a source and an interferer with powers σ 2 S1 and σ 2 S2 , located at θ 1 and θ 2 , respectively. The covariance matrix given by has the n − 2 smallest eigenvalues λ 3 = . . . = λ n = σ 2 and the two largest eigenvalues respectively. Applying the result in Proposition 3 allows us to express the effective number of signals as In the special situation when v 1 = v 2 = v and σ 2 S1 = σ 2 S2 = σ 2 S , we can (in an asymptotic sense) reliably detect the presence of both signals from the sample eigenvalues alone whenever Equation (16) captures the tradeoff between the identifiability of two closely spaced signals, the dimensionality of the system, the number of available snapshots and the cosine of the angle between the vectors v 1 and v 2 . It may prove to be a useful heuristic for experimental design.

SIMULATIONS
Assume the covariance matrix R has n−2 "noise" eigenvalues with σ 2 = 1, and two "signal" eigenvalues with λ 1 = 10 and λ 2 = 3. When m = 4n samples are available, Figure 1(a) shows that the proposed estimator consistently detects two signals while the WK MDL estimator does not. However, when m = n/4, Figure 1(a) suggests that neither estimator is able to detect both the signals present. A closer examination of the empirical data presents a different picture. For the covariance matrix considered, when m = n/4, then from (12), k ef f = 1. Figure 1(b) shows that for large n and m = n/4, the new estimator consistently estimates one signal, as expected. The WK MDL estimator detects no signals. We conjecture that the new estimator consistently estimates k ef f in the n, m → ∞, n/m → c sense.

CONCLUSIONS
An estimator for the number of signals in white noise was presented that exhibits robustness to high-dimensionality, and sample size constraints. The concept of effective number of signals described provides insight into the (asymptotic) regime in which reliable detection with sample eigenvalue based methods, including the proposed method, is possible. This helps identify scenarios where algorithms that exploit any structure in the eigenvectors of the signals, such as the MUSIC and the Capon-MVDR [5] algorithms in sensor array processing, might be better able to tease out lower level signals from the background noise. It is worth noting that the proposed approach remains relevant in situations where the eigenvector structure has been identified. This is because eigen-inference methodologies are inherently robust to eigenvector modelling errors that are endemic to high-dimensional settings. Thus the practitioner may use the proposed estimator to complement and "robustify" the inference provided by algorithms that exploit the eigenvector structure.


doi:10.1155/2007/57054 Research Article Distributed Time Synchronization in Wireless Sensor Networks with Coupled Discrete-Time Oscillators

In wireless sensor networks, distributed timing synchronization based on pulse-coupled oscillators at the physical layer is currently being investigated as an interesting alternative to packet synchronization. In this paper, the convergence properties of such a system are studied through algebraic graph theory, by modeling the nodes as discrete-time clocks. A general scenario where clocks may have different free-oscillation frequencies is considered, and both time-invariant and time-variant network topologies (or fading channels) are discussed. Furthermore, it is shown that the system of oscillators can be studied as a set of coupled discrete-time PLLs. Based on this observation, a generalized system design is discussed, and it is proved that known results in the context of conventional PLLs for carrier acquisition have a counterpart in distributed systems. Finally, practical details of the implementation of the distributed synchronization algorithm over a bandlimited noisy channel are covered.

INTRODUCTION
Distributed timing synchronization refers to a decentralized procedure that ensures the achievement and maintenance of a common time-scale (frequency and phase) for all the nodes of the network [1]. This condition enables a wide range of applications and functionalities of a sensor networks, including complex sensing tasks (distributed detection/estimation, data fusion), power saving (all nodes sleep and wake-up at coordinate times), and medium access control for communication (e.g., time division multiple access and cooperative communications).
Conventional design of distributed algorithms for timing synchronization prescribes the exchange of local time information through packets carrying a time-stamp to be appropriately elaborated by the transmitting and receiving nodes [1]. Packet-based synchronization has been widely studied, especially in the context of wireline networks. However, the specific features and requirements of wireless sensor networks call for alternative methods that improve both the computational complexity (and therefore energy efficiency) and scalability. Toward this goal, physical layer-based synchronization protocols are currently being investigated that exploit the broadcast nature of radio propagation. The idea is to build distributed algorithms based on the exchange of pulses at the physical layer, thus avoiding the need to perform complex processing at the packet level.
Physical layer-based synchronization was studied in [2] using a mathematical framework developed in [3] in order to model the spontaneous establishment of synchronous periodic activities in biological systems, such as the flashing of fireflies. In [2,3], nodes are modeled as integrate-and-fire oscillators coupled through the transmission of pulses. Convergence is proved under the assumption of an all-to-all interconnection among the nodes. The model was later extended in [4], by explicitly including constraints on the transmission range of each node. In particular, the authors derived a bound on the velocity of convergence by using algebraic graph theory [5]. An implementation of distributed synchronization on a real sensor network testbed was reported in [6]. A related work is [7], where a generalization of the model in [3] is proposed and the regime of an asymptotically dense network is investigated. As a final remark, it should be noted that the framework of physical layer-based timing synchronization has been recently interpreted as a means to achieve distributed estimation/detection [8,9] or data fusion [10].

EURASIP Journal on Wireless Communications and Networking
In this paper, we reconsider physical layer-based synchronization by modeling the sensors as coupled discrete-time oscillators. Basically, each node modifies its current clock based on a weighted average of the residual differences of timing phases as measured with respect to other nodes. The synchronization algorithms proposed in [11] in the context of interbase station communication and [12] for intervehicle transmission can be seen as instances of this general model. The analytical framework is at the same time a generalization and an application of the literature on discrete-time consensus problems for networks of agents (see, e.g., [13]). In particular, differently from [13], here we address the case of clocks with generally different free-oscillation frequencies, and account for the specific features of a wireless network, namely channel reciprocity and randomness (fading). Analysis of convergence of the synchronization process is carried out by algebraic graph theory as in [4], allowing to relate global convergence properties to the local connectivity of the network. The results are first derived for a time-invariant scenario, and then extended to the case where the network topology (or fading) varies with time, building on the results presented in [14].
A central contribution of this paper is the observation that the distributed synchronization system at hand can be modeled as a set of coupled discrete-time phase locked loops (PLLs). The system can thus be seen as a discrete-time version of the network synchronization scheme of [15], that is based on continuously-coupled analog PLLs. This fact allows us to generalize the system design by introducing the concept of loop order. Moreover, we prove that known results about the convergence of conventional PLLs for carrier acquisition have a counterpart in distributed systems. In particular, it is shown that, under appropriate conditions on the interconnections between sensors, (i) a system of first-order distributed PLLs is able to recover perfectly a phase mismatch among the clocks; (ii) in case of a frequency error, first-order loops are able to recover the frequency gap, but at the expense of an asymptotic phase mismatch; (iii) this asymptotic phase mismatch can be reduced by considering second-order loops.
Finally, the analysis is complemented by addressing the issue of a practical implementation of the distributed synchronization algorithm over a bandlimited Gaussian channels.

SYSTEM MODEL AND MAIN ASSUMPTIONS
Let the wireless network be composed of K sensors, where each node, say the kth, has a discrete-time clock with period T k . If the nodes are left isolated, the timing clock of the kth sensor evolves as t k (n) = nT k +τ k (0), where 0 ≤ τ k (0) < T k is an initial arbitrary phase and n = 1, 2, . . . runs over the periods of the timing signal. Two synchronization conditions are of interest. We say the K clocks are frequency synchronized if for each k and for sufficiently large n, where 1/T is the common frequency. A more strict condition requires full frequency and phase synchronization 1 : We remark that the network is said to fractionate into, say, two clusters of synchronization if there exist a permutation function on the nodes' labels, π(i) : [1, . . . , n] → [1, . . . , n] such that for n large enough where the number of nodes in the two clusters is r and K − r, respectively. The definition above generalizes naturally to more than two clusters. Towards the goal of achieving synchronization, the clocks of different sensors can be coupled by letting any node radiate a timing signal as the one sketched in Figure 1. A pulse 2 is transmitted at times t k (n) by the kth node and received through independent flat fading channels by the other sensors. It is assumed that all the nodes transmit with the same power, and that the power P ki received on the wireless link between the ith and the kth user reads where C is an appropriate constant that depends on the transmitted power (assumed here to be the same for all nodes), d ki (n) = d ik (n) is the distance between node i and node k at the nth period, G ki (n) is a random variable accounting for the fading process, and γ is the path loss exponent (γ = 2 ÷ 4). Notice that the fading channel is reciprocal (all transmissions use the same carrier frequency), which implies that G ik (n) = G ki (n) and P ik (n) = P ki (n) for i / = k [16]. As detailed in the following, each node (at any period n) processes the received signal in order to estimate the time difference between its own clock t k (n) and the corresponding "firing" instant of other nodes, that is, t i (n)−t k (n), i / = k, and, based on this measure, it updates its own clock.

The synchronization algorithm
In this section, we consider the synchronization procedure under the ideal assumptions that any node, say the kth, is able to measure exactly the time differences t i (n) − t k (n) and the powers P ki (n) of other nodes (i / = k) based on the received signal. This model is elaborated upon in the first part of the Figure 1: Clock t k (n) of the kth node. τ k (n) is the timing phase in the nth period of the clock.
paper. A practical implementation of the system that alleviates the said assumptions (and in particular, does not require estimation of time of arrivals) is then discussed in Section 7. At the nth period, the kth node updates its clock t k (n) according to a weighted sum of timing differences Δt k (n + 1) 3 : where ε is the step-size (0 < ε < 1) and the coefficients α ki (n) are selected so that α ki (n) ≥ 0 and K i=1, i / =k α ki (n) = 1. The updating rule (5) generalizes the algorithms of [11,12] (and the consensus algorithms, see, e.g., [13]) to a frequencyasynchronous scenario. In this paper, we focus on the following choice for the coefficients α ki (n): The selection of the weighting coefficients (6) is inspired by the algorithms proposed in [11,12]. The rationale of this design is that time differences measured over more unreliable (i.e., low-power) channels should be weighted less when updating the clock, thus rendering the algorithm robust against measurement errors over the fading channels (see also Section 7). Notice that by using (5b) we are implicitly neglecting the propagation delays among nodes, that are assumed to be smaller than the timing resolution. A method to handle propagation delays is described in [11]. As a final remark, we notice that the dynamic system (5) updates the clock t k (n+1) as a convex combination of the times {t i (n)} K i=1 [14].
By defining the vector containing the clocks of all nodes as t(n) = [t 1 (n) · · · t K (n)] T and the vector of clock periods T = [T 1 · · · T K ] T , we can express (5) as the difference vector equation where A(n) is a K × K matrix such that we have [A(n)] ii = Notice that even though we assume channel reciprocity, matrix A(n) is not symmetric. Moreover, by construction, matrix A(n) is nonnegative and stochastic since the sum of the elements on each row sums to one, or equivalently

TIME-INVARIANT FREQUENCY-SYNCHRONOUS NETWORK
In this section, we study the convergence properties of the distributed synchronization algorithm (5) under the following assumptions: (i) frequency-synchronous network, that is, all the clocks share the same period T = T 1 = · · · = T K ; (ii) the network is time-invariant, that is, P ki (n) = P ki for any n and k / = i. From assumption (i), the clock of the kth node can be expressed as where τ k (n) is the timing phase 0 ≤ τ k (n) < T of the kth node in the nth period (see Figure 1). Moreover, by substituting (9) into (5a) and using assumption (ii), it easily follows that the synchronization algorithm (5) can be written in terms of the phases τ k (n) as τ k (n + 1) = τ k (n) + ε · Δτ k (n + 1), (10a) with coefficients α ki : Finally, by defining the vector containing the timings of all nodes as τ(n) = [τ 1 (n) · · · τ K (n)] T , the vector model (7) becomes where A is a K × K matrix such that we have [A] ii = 1 − ε on the main diagonal and [A] i j = ε · α i j for i / = j. Model (12) resembles the one considered in the literature on multiagent coordination (see, e.g., [13]). The goal of this section is to determine the conditions under which the system (12) converges to a unique cluster or to multiple clusters of synchronization for a fixed realization of the fading variables G ki in (4), that is, matrix A is assumed to be deterministic. We will define the conditions of convergence in terms of the properties of the graph associated to the wireless network under study, or equivalently in terms of the system matrix A.

The associated graph and useful definitions
The synchronization algorithm defines a weighted directed graph G = (V, E , A) of order K on the sensor network, where V ={1, . . . , K} is the set of nodes and E ⊆ V × V is the set of edges weighted by the off-diagonal elements of the K × K adjacency matrix [A] i j = α i j . The edge connecting the ith and the jth nodes, i / = j, belongs to E if and only if α i j > 0. Notice that the graph is directed (α i j / = α ji for i / = j), even though fading links are reciprocal (P i j = P ji for i / = j). Moreover, notice that the system matrix reads where L is the graph Laplacian of the network that is defined as [13]: The main result of this section (Theorem 1) relates the convergence properties of the distributed synchronization procedure in (10) with the connectivity of the graph G associated to the sensor network. We need the following definitions.

Definition 1.
A graph G is said to be strongly connected if there exists a path (i.e., a collection of edges in E ) that links every pair of nodes.
It can be proved that strong connectivity of graph G is equivalent to the irreducibility of matrix A [18].

Definition 2.
A K × K matrix A is said to be reducible if there exists a K × K permutation matrix P and an integer r > 0 such that where The degree of irreducibility of a matrix A, or equivalently of strong connectivity of the associated graph G, can be measured by the following quantity where the minimum is taken over all nonempty proper sub- It can be shown that σ = 0 if and only if the matrix A is reducible, or the associated graph G is not strongly connected [19]. 4

Convergence properties
The main result of this section can be now stated as follows.

Theorem 1. (i)
The distributed synchronization (10) converges to a unique cluster of synchronized nodes, τ 1 (n) = · · · = τ K (n) = τ * for n → ∞, if and only if the associated weighted directed graph G is strongly connected, or equivalently if system matrix A is irreducible. (ii) In this case, the system (12) converges to (for n → ∞) An immediate consequence of Theorem 1 is that the timing vectors converge to the average of their initial values τ(0) if and only if the system matrix A is doubly stochastic (i.e., if A T is stochastic as well). In fact, in this case A T 1 = 1 and vector v in (16) reads v = 1/K·1. This condition occurs in balanced networks [13], where i / = j α i j = 1 = i / = j α ji . In sensor networks, this result is of interest in applications where the steady state value of synchronization is used in order to infer the status of the process monitored by the sensor [8,9,20].
Proof. The proof of part (i) of Theorem 1 is available in the literature for applications where the graph G associated to the dynamic system (12) is undirected [5]. In the case of a directed graph, strong connectivity can generally be proved to be only a sufficient condition for synchronization. However, in a wireless fading case with reciprocal channels, the result can be proved as shown in the following. The second part (ii) of Theorem 1 follows from a result derived, among the others, in [13].
As explained above, in order to prove Theorem 1, we only need to show that strong connectivity is also a necessary condition for synchronization. As a by-product, the proposed proof brings insight into the formation of multiple clusters of synchronization (3). Let us assume that A is reducible (or equivalently the associated graph G is not strongly connected). Then, by definition, there exists a permutation matrix P and an integer r > 0 such that (14) holds. But if α i j = 0 in A, then for reciprocity P i j = P ji = 0 and then α ji = 0 (i / = j). This property is sometimes referred to as bidirectionality of the graph (i.e., α i j = 0 if and only if α ji = 0 but α i j and α ji need not to be equal [14]). Therefore, the (14) has all zero entries. Since the permuted matrix P T AP is nonnegative and stochastic, so are submatrices B and D. By applying the permutation function π(k) = P k [1 · · · K] T , where P k is the kth row of matrix P, to the nodes' labels, we can write the system (12) as where τ(n) = Pτ(n). Therefore, the set of r nodes {π(1), . . . , π(r)} evolves independently from the remaining nodes {π(r + 1), . . . , π(K)}. Now, if either B or D are reducible, the reasoning above can be iterated bringing to the formation of multiple independent set of nodes evolving separately. At the end of this procedure, the system matrix can be O. Simeone and U. Spagnolini Figure 2: The rectangular topology considered in the example in Section 3.3.
written as a block matrix with irreducible stochastic blocks on the diagonal. Without loss of generality, let us then assume that B and D are irreducible. From the first part of the proof (see also Appendix A), it follows the two cluster of r and (K − r) nodes synchronize among themselves according to (3). Moreover, the steady state values of the timing vectors depend on the left eigenvectors of B and D according to (16): where is the r × 1 vector collecting the first r entries of τ(n) and The convergence of the dynamic system at hand could be also studied in terms of the subdominant eigenvalue of matrix A, similarly to approach commonly adopted in the context of the analysis of Markov chains [21]. In particular, the following results can be proved relating convergence to the multiplicity of eigenvalue 1.

Theorem 2.
The distributed synchronization (10) converges to a unique cluster of synchronized nodes as in (2) if and only if the subdominant eigenvalue λ 2 / = 1.
Proof. By recalling Theorem 1, it is enough to prove that: (i) if λ 2 = 1, then the graph is not strongly connected; (ii) if the graph is not strongly connected then, λ 2 = 1. Part (i) can be proved similarly to [13]; however, in Appendix A we give an alternative proof based on the measure σ in (15) of irreducibilty of A. Part (ii) does not hold in general for problems with directed graphs but it is easily shown under the reciprocity assumption similarly to Theorem 1.

Numerical results
Here, we present a numerical example to corroborate the analysis discussed above. A network of K = 4 nodes is considered where the nodes are divided into two groups,  (11), the performance is not affected by the value of C in (4) and it only depends on relative distances. Figure 3 shows the timing vector

τ(n) versus n.
After a transient where the nodes tend to synchronize in pairs within the two groups, the system reaches the steady state to the average value τ * /T = 0.475, as stated in Theorem 1, since the system matrix is easily shown to be doubly stochastic for this specific example.
In order to quantify the rate of convergence, from Theorem 2, we notice that the convergence of the synchronization protocol (10) depends on the subdominant eigenvalue λ 2 . In particular, as it is well known from the theory of linear difference equations, the rate of convergence is ruled by a term proportional to |λ 2 | n (see, e.g., [22]). If we define a threshold λ o , we could say that the protocol reaches the steady state condition at the time instant n o for which |λ 2 | no = λ o : n o = log λ o / log |λ 2 |. Therefore, we can take v = − log λ 2 (19) as a measure of the rate of convergence of the algorithm. Figure 4 shows the rate of convergence v versus the normalized distance D/d for ε = 0.3, 0.7. As expected the rate v decreases with increasing D/d and decreasing ε. Along with v, Figure 4 shows the measure of irreducibility (or strong connectivity) σ (15) as dashed lines. It is interesting to note that the rate of convergence v and the measure of irreducibility σ have the same behavior as a function of D/d and ε. This confirms that convergence is strictly related to the connectivity properties of the associated graph, as proved in Theorem 1.

Effect of fading: an example
In this section, the effect of fading on the rate of convergence v is investigated via simulation for linear, ring, and star topologies (see Figure 5). Rayleigh fading is assumed, that is, the fading amplitude G ki in (4) accounts for Rayleigh fading with unit average power. Fading is assumed to be constant for any n during the evolution of the algorithm.  three networks versus the number of nodes K (ε = 0.3). Notice that for K = 2 the three networks coincide and recall that only relative distances are of concern for the behavior of the system (10). As it is expected, the star topology has the largest rate of convergences whereas the linear network yields the slowest convergence.

TIME-VARYING FREQUENCY-SYNCHRONOUS NETWORK
Here, we reconsider the performance of the synchronization procedure (5) by removing the assumption of timeinvariance underlying the analysis of the previous section. However, we still assume a frequency-synchronous network. Overall, the system (5) can be written in vector form in terms of phases as (recall (12)) with the definition of the system matrix A(n) in Section 2.
In this case, the sensor network can be described by a sequence of directed graphs G(n) = (V, E (n), A(n)) of order K, defined similarly to Section 3. In particular, A(n) is the adjacency matrix [A(n)] i j = α i j (n) and the edge connecting the ith and the jth nodes, i / = j, belongs to E (n) if and only if α i j (n) > 0. At each time instant n, the dynamic system describing the synchronization process evolves as where A(n) = I−εL(n) with L(n) being the graph Laplacian at time n (see Section 3). Study of convergence of a family of algorithms encompassing (10) has been recently attempted in a few works (see [13,14] and references therein). In particular, adapting a result first presented in [14] to our case, we are able to relate the convergence of dynamic system (20) to the connectivity properties of the associated sequence of graphs G(n). We need the following definition.

Theorem 3. The distributed synchronization
Proof. Theorem 3 can be proved by specializing the proof of [14, Theorem 3] to our scenario. The basic idea is to exploit the convexity, or the contractive property, of transformation (10). An interesting remark is that reciprocity of fading plays here a key role as it did in the proof of Theorem 1 for the case of fixed topology. Reciprocity of fading translates into bidirectionality of the associated graphs (see Section 3). As proved in [14], in presence of unidirectional communication among nodes (i.e., nonreciprocal fading in our scenario), convergence of synchronization is not necessarily guaranteed if each sensor communicates to every other sensor (either directly or via intermediate nodes) in an interval [n 0 , ∞). On the contrary, in order to guarantee convergence in a unidirectional graph, a limit should be imposed on the time it takes for the graph to become strongly connected, that is, the interval in Theorem 3 should be modified as [n 0 , n 0 + T] where T ≥ 0 finite.

FREQUENCY-ASYNCHRONOUS NETWORK
In the previous sections, it was assumed that all the nodes have the same clock period T (frequency synchronous network). However, in practice, different nodes might have different frequencies {1/T k } K k=1 , and the question arises of whether or not the physical layer-based scheme (5) is still able to achieve synchronization on a strictly connected graph. For a time-invariant scenario (i.e., P ki (n) = P ki for any n and k / = i), it will be shown below that, in presence of a frequency mismatch, the scheme (10) is able to synchronize the clock periods of the nodes (recall (1)), but not their timing phases, so that the full synchronization condition (2) is not achieved.
In this regard, it should be noted that, while perfect synchronization (2) is necessary for many applications, in other scenarios having nodes with synchronized frequency is the only requirement (i.e., to ensure equal sensor duty cycles).
For a frequency-asynchronous time-invariant network, the considered synchronization scheme (5) reads or, in vector form (see (7)): Let us now denote a possible common value for the clock period of all nodes as T (to be determined) as in (1). It follows that the clock of the kth sensor can then be written for sufficiently large n as or equivalently, in vector form, as t(n) = nT · 1 + τ(n) with τ(n) = [τ 1 (n) · · · τ K (n)] T . We are interested in determining if such common frequency 1/T exists and, if so, whether eventually the phases τ(n) converge to the same value for n → ∞. The main conclusion is summarized in the theorem below, whose proof is inspired by the analysis of the convergence of coupled analog oscillators in [23].

Theorem 4.
With reference to (23), under the assumption that the graph G is strictly connected, the system (21) synchronizes the clocks of the K nodes to the common period where v is the normalized left eigenvector of matrix A corresponding to eigenvalue 1: However, the timing phases τ(n) remain generally mismatched and given for n → ∞ by with (·) † denoting the pseudoinverse and the definitions The theorem above states that, in presence of a frequency mismatch, the algorithm (21) is able to synchronize the frequencies of different nodes to the common clock period T in (24). However, the system does not lead to phasesynchronous clocks, and the phase error is determined by the frequency (period) mismatch ΔT according to (25). Notice that, if the network is such that the system matrix A is doubly stochastic (as in the example of Section 3.3), the eigenvector v reads 1/K and the common period T is in this case the average T = 1/K K k=1 T k . Moreover, with doubly stochastic matrix A, condition (25) simplifies as η = 1/K K k=1 τ k (0) since 1 T L † = 0 (see proof below for further details). Finally, we remark that if the frequency mismatch is ΔT = 0 (or equivalently T k = T), Theorem 4 follows from Theorem 1.
Proof. Under the assumption of a connected graph (or irreducible matrix A), according to Theorem 2 or [13], the Laplacian L is easily shown to have rank K − 1, where the one-dimensional null subspaces are defined by the relationships Using the latter equality, recalling (13) and the definition of common clock period T and phases τ(n) in (23), the vector difference equation (22) can be written as An equilibrium state τ(n + 1) = τ(n) = τ * for the difference equation (29) satisfies τ(n + 1) − τ(n) = 0, which yields the condition From (28), it follows that (i) in order for (30) to be feasible (i.e., for an equilibrium point to exist), the common clock period T must satisfy v T ΔT = 0 or equivalently (24); (ii) an equilibrium phase vector τ * must read τ * = (L † /ε)ΔT+η1 where η is an arbitrary constant. It remains to show that the system actually converges for n → ∞ to the equilibrium point τ * determined above, and to evaluate the constant η. Toward the goal of studying convergence, let us define τ (n) = τ(n) − (L † /ε)ΔT. With this change of variables, the difference equation (29) boils down to

EURASIP Journal on Wireless Communications and Networking
Timing error detection Loop filter where we used the relationship LL † ΔT = ΔT, which easily follows from the definition of pseudoinverse and (24). As a consequence of (31), as per Theorem 1, we have τ (n)→ v T τ (0). This expression is equivalent to (25), thus proving the theorem. Notice that from (31) the rate of convergence is the same as in the case of no frequency mismatch.
As a final remark, we notice that the study of timevarying frequency-asynchronous networks is a challenging task and is left for future work.

DISTRIBUTED TIME SYNCHRONIZATION AS COUPLED DISCRETE-TIME PLLs
The purpose of this section is to discuss the distributed synchronization algorithm investigated throughout the paper by casting it into the framework of discrete-time phase locked loops (PLLs) [24]. The discussion is not only be beneficial for a better understanding of the system, but it also allows us to generalize the system design. In order to appreciate the similarity with a discrete-time PLL, Figure 7 depicts the synchronization procedure (5) carried out at each sensor as a linear dynamic feedback system. Similarly to a discrete-time PLL, the adder at the input evaluates the timing error Δt k (n), that is then multiplied by ε and then fed to a voltage controlled clock (VCC) that updates the clock according to (5a). The constant ε plays the role of the loop filter in a discrete-time PLL. The procedure (5) can then be interpreted as a firstorder discrete-time PLL since the loop filter is a trivial pure gain [25]. From the discussion above, it is clear that the second or third order discrete PLLs 5 can be obtained by introducing a loop filter ε(z), with one or two poles respectively, instead of the constant ε in the synchronization system of Figure 7. For instance, in the case of a second-order loop, we can introduce a pole μ in the loop by setting ε(z) = ε/(1 − μz −1 ) with 5 As discussed in [25], loops are never built with order larger than three. 0 < μ < 1. The corresponding update rule (5a) modifies as t k (n + 1) = t k (n) + ε · Δt k (n + 1) The updating rule (32) essentially corrects the local period T k by the estimate of the common clock period t k (n) − t k (n − 1). The convergence analysis of the second-order loop (32) can be carried out similarly to the previous section where a firstorder loop was considered. In particular, the following results hold.

Theorem 5.
If the network of distributed PLL is strictly connected and the system (32) converges, then it synchronizes the clocks of the K nodes to the common period (24). However, under the same conditions, the timing phases τ(n) remain generally mismatched and given for n → ∞ by with (·) † denoting the pseudoinverse, and with definitions and (27).
Proof. The proof is along the lines of the proof of Theorem 4 (see Appendix B for details).
Comparing the statement of the previous theorem with the results derived for a first-order loop (Theorem 4), it can be seen that introducing a pole μ in the loop causes a reduction in the steady state phase error by a factor 1−μ. However, this reduction comes at the expense of decreased margins of stability. In fact, convergence cannot be guaranteed for all values of 0 < ε < 1 and 0 < μ < 1. We refer to Appendix C for further analysis on this point. Here, we illustrate this issue by means of an example. Consider a network with two nodes. In this case, we have α 12 = α 21 = 1 and the graph is connected. Figure 8 shows the four eigenvalues of the system matrix associated with (32) (see Appendices B and C for further details): for different values of the pole μ and ε = 0.9. Notice that the system matrix (35) is 4×4 since (32) is a system of two second order difference equations [22]. Moreover, one eigenvalue of A is 1 irrespective of the value of μ. The absolute value of the remaining eigenvalues tends to one for μ → 1, showing that increasing the value of the pole leads to lack of stability of the equilibrium point (33). Moreover, the value of μ at which a couple of eigenvalues acquire a nonzero imaginary part can be calculated exactly as a function of the spectrum of matrix A, as shown in Appendix C (see (C.5)). In order to corroborate the conclusions above, Figure 9 shows the standard deviation ξ(n) of the timing vector t(n) O. Simeone and U. Spagnolini To conclude, it is interesting to revisit the results of Theorems 1, 4, and 5 in the light of the analogy with conventional PLLs drawn above. It has been shown that in a strictly connected network: (i) a phase error is perfectly recovered for n → ∞ by the distributed synchronization algorithm (5) (Theorem 1); (ii) a frequency error is perfectly recovered at the expense of a phase mismatch for n → ∞ (Theorem 4); (iii) the residual phase mismatch caused by a frequency error can be reduced by introducing a pole in the control loop (Theorem 5). All these results can be read as the counterpart of known facts in the analysis of linearized PLLs, which assert that a first-order loop is indeed able to (i) recover phase errors and (ii) to achieve a constant phase error (referred to as static phase error) in case of a frequency mismatch [25]. Moreover, in this second case, it is interesting to notice how the phase errors (25) and (33) depend on the frequency mismatch ΔT exactly as the static phase error of a PLL [25]. Further results on large-scale randomly deployed networks can be found in [26].

IMPLEMENTATION OF DISTRIBUTED COUPLED DISCRETE-TIME OSCILLATORS
In the previous sections, it was assumed that each node is able to measure time differences and powers of other nodes so as to calculate the phase update Δt k (n). Here, we remove this assumption by presenting a practical scheme to implement the phase detector over a bandlimited noisy channel. Since the algorithm is based only on instantaneous power measurements by different nodes, it applies to both a time-invariant and time-variant scenario. The scheme is inspired by the proposal in [12]. A carrier frequency is dedicated to the synchronization channel, where each node, say the kth, transmits a bandlimited waveform g(t) (say a square-root raised cosine pulse) centered at times t k (n) with symbol period 1/F s (i.e., the time between peak and first zero). The symbol period 1/F s defines the timing resolution of the system. Each node works in an half-duplex mode and measures the received signal on a interval of duration T k around the current timing instant t k (n). Due to the half duplex constraint and the finite switching time between transmitting and receiving mode, each sensor is not able to measure the received signal in an interval of (unilateral) size θ around the firing instant t k (n). It follows that the observation win- Figure 10(a) illustrates a block diagram of the operations performed at the receiver side by each sensor. The receiver performs baseband filtering matched to the transmitted waveform and than samples the received signal at some multiple L of the symbol frequency F s , that is, LF s with L ≥ 1. Based on the N = LF s T samples received in the nth observation window, the kth node computes the update (21a) in case a first-order loop is employed, or (32) if a second-order loop is considered. Not knowing the exact timings and powers of other nodes, t i (n) and P ki with i / = k, the kth sensor cannot directly calculate the updating term Δt k (n) from (5b) and (6). Instead, it estimates these quantities from the received samples, as explained below. After matched filtering and sampling, the discrete-time baseband signal received by the kth node in the nth time period reads (sampling index m ranges within −N/2 < m ≤ N/2 with m = 0 corresponding to the firing time t k (n) of the kth node): where the average energy per symbol reads E ki = C/(d γ ki · F s ) (recall (4)); β ki denotes the Rayleigh fading coefficient, that is a zero-mean and unit-power complex (circularly symmetric) Gaussian random variable with |β ki | 2 = G ki ; and w(n, m) is the additive Gaussian noise with zero mean and power N 0 . Notice that the sample in the interval −ΔLF s ≤ m ≤ ΔLF s is not measured (i.e., zero) due to the half-duplex constraint and the switching time between receive and transmit mode of node k. A sketch of a possible realization of the received signal (36) is provided in Figure 10(b) using an arbitrary waveform g(t).
A simple estimate of ΔT k (n) can then be obtained as where J is the subset of time instants m ∈ (−N/2, −θLF s ) (θLF s , N/2] for which the received signal |y k (n, m)| 2 is above a given threshold as in [12]. The threshold is a system parameter that can be optimized. Being based solely on instantaneous power measurements (i.e., on samples |y k (n, m)| 2 ), the practical scheme (37) has the advantage that it does not require any a priori knowledge at the nodes about network topology or average received powers. For the example of Section 3.3 with no fading (β ki = 1 for every i and k), Figure 11 shows the standard deviation ξ(n) of the timing vector t(n) versus n, where ξ 2 (n)

and expectation E[·]
is taken with respect to noise. We are considering equal clock periods T k = T = 1 for k = 1, . . . , K. Moreover, it is assumed that all nodes transmit the same power and the signal-tonoise ratio for transmission to the closest node (e.g., from 2 to 1) is set to SNR = E 12 /N 0 = 15 dB. Other parameters are as follows: ε = 0.9; the threshold is set for simplicity to zero (see discussion below); distances satisfy D/d = 2; the normalized timing resolution is 1/F s = 0.01; the waveform g(t) is a raised cosine with roll-off δ = 0.2; the switching time is set to θ = 1/F s . 6 We first consider the first-order loop (21a) (or equivalently (32) with pole μ = 0) with different oversampling factors L = 1, 2, 5, 10, 15. It can be seen that the finite resolution of the system produces a performance floor for increasing n, that can be lowered by increasing the oversampling factor L. In any case, an upper bound on the accuracy of synchronization is set by the finite switching time θ = 0.01. This bound is reached for n and L sufficiently large. 7 Adding a pole in the loop as in (32) can increase the convergence speed as shown in Figure 11 for μ = 0.2, 0.4, 0.6. Notice that convergence speed could also be improved by setting an appropriately chosen threshold in (37) (not shown). Finally, Figure 11 shows that an upper bound on the performance of the practical implementation discussed here is set by the performance of the dynamic system (10), where the performance gap is due to noisy observations, finite resolution and finite switching time.

CONCLUSIONS
An increasing number of applications of sensor networks requires the availability of a common time reference to all the nodes. Aiming at scalability and complexity reduction, physical-layer based synchronization qualifies as a valid alternative to conventional packet-based synchronization. In this work, the convergence properties of such an approach have been investigated by modeling the network as a set of discrete-time coupled oscillators and relying on the analytical framework of algebraic graph theory. It has been shown that the system can be equivalently studied as a set of distributed discrete-time PLLs. This observation allowed to generalize the design of the synchronization process, and prove that well-established results on conventional PLLs extend naturally to distributed systems. Finally, practical details of the implementation of the synchronization algorithm have been discussed on a bandlimited Gaussian channel.

Lemma 1. Every eigenvalue of A satisfies
Moreover, if G is connected, only one eigenvalue λ = 1.
From (C.3) and (C.4), it can be seen that the eigenvalue λ of A (with the positive sign in (C.3)) corresponding to λ = 1 equals one. Moreover, by substituting λ = 1 in (C.3), this eigenvalue is unique. In order to prove stability of the system at hand, we then need to show that all the other 2K −1 eigenvalues λ(λ, μ) different from λ = 1 are such that | λ(λ, μ)| < 1. In general, it is not easy to find general conditions on ε and μ that guarantee the condition | λ(λ, μ)| < 1 for any λ in (C.4) (except λ = 1). However, a numerical study can be easily carried out from (C.3) and (C.4). For illustration, consider an undirected graph G. In this case, the eigenvalues of matrix A are real and from (C.4) satisfy −1 < λ < 1 (recalling that 0 < ε < 1). The absolute value | λ(λ, μ)| (C.3) (with positive sign, similar result holds for the negative sign) is plotted versus μ and λ in Figure 12. It can be seen that unless μ → 1, the stability condition is satisfied. Moreover, by studying the sign of the polynomial in (C.3), it is also possible to evaluate the value of the pole μ, say μ 0 (λ), such that for μ ≤ μ 0 (λ) the eigenvalues λ corresponding to λ are real and for μ > μ 0 (λ) are complex conjugate. This value is easily found to be μ 0 (λ) = 2 − λ − 2 1 − λ, ( C . 5 ) which is shown in Figure 12.
As a final remark, we notice that, for μ = 1, the absolute value | λ(λ, μ)| equals one irrespective of λ, leading to an unstable system.