Cycle flows and multistabilty in oscillatory networks: an overview

The functions of many networked systems in physics, biology or engineering rely on a coordinated or synchronized dynamics of its constituents. In power grids for example, all generators must synchronize and run at the same frequency and their phases need to appoximately lock to guarantee a steady power flow. Here, we analyze the existence and multitude of such phase-locked states. Focusing on edge and cycle flows instead of the nodal phases we derive rigorous results on the existence and number of such states. Generally, multiple phase-locked states coexist in networks with strong edges, long elementary cycles and a homogeneous distribution of natural frequencies or power injections, respectively. We offer an algorithm to systematically compute multiple phase- locked states and demonstrate some surprising dynamical consequences of multistability.


I. FROM KURAMOTO OSCILLATORS TO POWER GRIDS
Coupled oscillator models are ubiquitous in science and technology, describing the collective dynamics of various systems in micro to macro scale. Research on coupled oscillators dates back to Christian Huygens, who noticed that two clocks synchronize when they are coupled [1]. One of the most important mathematical models was introduced by Kuramoto [2,3] and successfully applied to describe the collective dynamics of coupled Josephson junctions [4], neuronal networks [5], chemical oscillators [6], and a variety of other synchronization phenomena [7][8][9][10].
That model [3] describes the dynamics of N coupled limit cycle oscillators. The equations of motions for the phases θ j , j ∈ {1, . . . , N } are given by The coupling matrix is assumed to be symmetric, K j, = K ,j and the ω j are the natural frequencies of the oscillators. Throughout this article we consider systems where all K j, ≥ 0.
A similar model of second-order oscillators describes the collective phenomena of animal flocks [11,12] or human crowds [13] as well as the coarse-scale dynamics of power grids [14][15][16][17][18][19][20]. For power grids, for instance, the units j describe synchronous machines, generators or motors, whose state is completely described by their phase θ j and the phase velocityθ j relative to the reference frequency of the grid, typically rotating at 50 Hz or 60 Hz. The acceleration (deceleration) of the machines is proportional to the sum of the mechanical power P j generated (consumed) by the machine including damping and the electric power exchanged with the grid. The detailed equations of motion are given by where M j is an inertia term and D j a damping constant. The coupling constants K j, = U 2 B j, are determined by the voltage U of the grid, which is assumed to be constant, and the admittance B j, of the electrical transmission line joining node j and node . The flow of electric real power from node to node j is F j, = K j, sin(θ − θ j ) = K j, S j, .
It is useful to describe the interaction topology of the system as a weighted graph G(V, E), whose vertex set V is identical with the set of oscillators, and edge set E is given by the set of all inter-oscillator coupling pairs, i.e., all pairs with K ,j > 0. We use the term network [21] (rather than the term graph) for the entire system with given natural frequencies ω j or the powers P j .
Here we distinguish two types of synchronization in oscillator networks. Traditionally, the emergence of partial synchrony has received the most interest in the physics community [2,3,7,8]. In his seminal work, Kuramoto investigated a set of oscillators with global coupling, K j, = K/N and natural frequencies drawn at random from a unimodal symmetric distribution g(ω). If the coupling constant K exceeds a critical value K c , a fraction of the oscillators start to synchronize in the sense that they rotate at the same angular velocity, although their natural frequencies differ. In this state, called "frequency locking", the phases of parts of the oscillators are ordered, but they are not strictly phase-locked, such that the phase difference of two synchronized oscillators (θ j − θ ) is generally small, but not constant.
In this article, we analyze the properties of globally phase-locked states, where all oscillators synchronize and the phase differences (θ j − θ ) are constant for all pairs (j, ). These states are especially important for power grids, as they describe the regular synchronous operation of the grid [14][15][16][17][18]. If this state is lost due to local outages or accidents, the grid will fragment into asynchronous islands which can no longer exchange electric energy [22]. For instance, the European power grid fragmented into three asynchronous areas on November 4th 2006 after the shutdown of one transmission line in Northern Germany. As a result, south-western Europe suffered an under-supply on the order of 10 GW and approximately 10 million households were disconnected [23].
Without loss of generality, we take j ω j = 0 or j P j = 0 respectively, by invoking a transformation to a co-rotating frame of reference. The globally phase-locked states are then the fixed points of the system. Both for the Kuramoto model and the power grid model, these states are given by the solutions of the algebraic equations replacing P j by ω j for the Kuramoto model. In the following we analyze the influence of the network topology given by the coupling matrix K j, on the existence of a fixed point. All results below hold for both models, nevertheless our intuition heavily relies on the interpretation of F j, = K j, sin(θ − θ j ) as a flow which is inspired from the power grid model. The results generalize to arbitrary coupling functions f instead of the sine (see, e.g., [24,25]).
In the following we mostly restrict ourselves to the common sine coupling for the sake of clarity.

II. THE NATURE AND BIFURCATIONS OF FIXED POINTS
Both the Kuramoto system and the oscillator model of power grids share the same set of fixed points (4). It has been shown that the similarity between these two systems runs deeper, namely, the linear stability properties of those fixed points are identical [26]. In this section we briefly review some basic results on the stability of the fixed points.
We analyze the dynamical stability of a certain fixed point θ * = (θ * 1 , . . . , θ * N ) by defining the potential function The fixed points correspond to the local extrema of this potential, where ∂V ∂θ j = 0 for all j.
with the residual capacity has positive eigenvalues only. It is worth noting that M has one eigenvector v 1 = (1, 1, · · · , 1) with eigenvalue µ 1 = 0, because any fixed point θ * is arbitrary up to an additive constant c. As such a global phase shift does not affect the locking of the phases we can discard it in the following and concentrate on the stability transversely to the solution space {θ * + c(1, 1, · · · , 1)|c ∈ R}.
Lemma 1. Let the eigenvalues of M be ordered such that µ 1 = 0 and µ 2 ≤ · · · ≤ µ N . If for a given network topology and a given fixed point, then this fixed point is transversally asymptotically stable for both Kuramoto system and the power grid model system. If one of the µ k < 0, then the dynamical system is linearly unstable.
Using some results from bifurcation theory, it has been shown in [26] that a stable fixed point can only be lost by an inverse saddle-node bifurcation when one of the eigenvalues becomes zero, µ 2 = 0. At this point linear stability analysis is not sufficient to predict stability of the fixed point but it is expected that the fixed point is unstable [27].
More insights can be gained about the loss of a fixed point when the phase differences across all edges in the network are sufficiently small: holds for all edges (i, j) in the network. Then the network is said to be in "normal operation".
Proof. To this end, we first define a meta graph as follows.
Definition 1 (Meta graph). Given a graph G(V, E), and a set of flows F uv across each edge e(u, v), its meta graphG is an undirected graph with vertex set V and edge set E defined as follows. For all edges e(u, v) ∈ E, with weight K uv , ∃ an edge e(u, v) ∈ E with weight K red uv = K 2 uv − F 2 uv , as per (7).
Then the matrix M as defined in (6) is seen to be the Laplacian matrix of the meta graphG. The eigenvalues of a Laplacian of a connected undirected graph with positive edge weights are always non-negative [21] such that we obtain the result.
During normal operation an eigenvalue of M can become 0 only whenG disconnects into two (or more) components. Such a split-up happens only when K red j, = 0 for all the transmission lines connecting two certain parts (denoted by G 1 , G 2 ) of the network, meaning that these lines are completely saturated Another scenario for the loss of stability is that one or more transmission lines leave normal operation. Then the edge weights become effectively negative, such that a simple graphtheoretic interpretation of the bifurcation is no longer possible [26].

A. Flow conservation and the dynamics condition
It is instructive to divide the defining equation (4) of a fixed point into two parts. First, every fixed point has to satisfy a dynamic condition which is nothing but the conservation of the flow at every node of the network Here, K j, S j, is the sum of all flows from the neighboring nodes to the node j, while P j is a source or sink term, respectively. The second part of this condition reflects the fact that the transmission capacity of each link is bounded, such that the magnitude of the flow |F j, | cannot exceed the capacity K j, . The dynamic condition (11) holds for all flow networks including also DC networks (i.e. Kirchhoff's rules) and biological network models [28,29].
To obtain a better understanding of the possible solutions, we slightly rephrase the dynamic condition (11). In particular, we label all the L edges in the network with e ∈ {1, . . . , L}. As the flows are directed, we have to keep track of the ordering of the vertices connected by the edge e. That is, each e corresponds to a directed link (j, ) in the following.
The ordering is arbitrary but must be kept fixed. Then we write S e = S j, and F e = F j, for the flow over a link e = (j, ). Furthermore, we define the unweighted edge incidence matrix and the weighted edge incidence matrix K ∈ R N ×L with the components K je = K e I je .
The dynamic condition (11) then reads in terms of the sine factors. Here, F = (F 1 , . . . , F L ) T and S = (S 1 , . . . , S L ) T are vectors in R L . The matrix K has N rows, but its rank is only (N − 1). This is due to the fact that the sum of all rows is zero as j K j,e = 0, since each edge has exactly one head and one tail. Hence, the solutions of the linear set of equations (14a) span an affine subspace of R L whose dimension is (L − N + 1). This statement will later be rigorously proved in Lemma 2. In many important applications L is much larger than the number of nodes N , such that we have a high dimentional submanifold B of R L with every S ∈ B a solution of (14), and hence, a fixed point of (1) and (2). However, the set of solutions of the dynamical equations can also be empty if the capacities K j, are too small. In fact, the condition (14b) defines a bounded convex polytope in R L . The solution of the full dynamical conditions (14) are given by the intersection of this polytope and the (L − N + 1) dimensional affine subspace.
We can further characterize the solution of the dynamic conditions, by establishing that the homogeneous solutions of the system (14a) are just the cycle flows which do not affect flow conservation. As the number of fundamental cycles in a network is (L − N + 1), the dimension of the solution space is also given by (L − N + 1). Derivation of these results follows.
Definition 2 (Simple cycle). Given an undirected graph G(V, E), a closed path c = where no vertex apart from v 1 occurs twice is called a simple cycle [30, Definition 4 (Signed characteristic vector of a cycle). An arbitrary assignment of a direction to each edge of an undirected graph G, which results in a directed graph, is called an orientation G σ [31]. Given a graph G with L edges and N vertices, and one such orientation, there exists an injective mapping from the set C of all simple cycles of G to R L as follows: z c is called the signed characteristic vector of each cycle.
Now we show that any fixed point of the system can be uniquely specified by a cycle flow along each cycle belonging to a cycle basis of the underlying graph, alongwith an arbitrary solution of (13).
Definition 5 (Cycle flow). Given a simple cycle c = (v 1 , v 2 , · · · , v l , v 1 ) belonging to an i.e. it is a contant nonzero flow along the cycle. Let e L ) and F = (F e 1 , F e 2 , · · · , F e L ) be the flows for the fixed points θ (0) and θ, respectively. Then due to the result from graph theory that the flow space of an oriented graph G σ is spanned by the signed characteristic vectors (Definition 4) of its cycles [31, p 311]. Since by definition the cycles in B C forms a basis of the cycle space, the coefficients f c are guaranteed to be unique. This concludes the proof.

B. The winding number and the geometric condition
In addition to the dynamic condition, there is a geometric condition for the existence of a fixed point: a fixed point exists if the flows F j, = K j, S j, satisfy the dynamic condition (14) and if for all edges ( , j) : ∃ (θ 1 , . . . , θ N ) such that S j, = sin(θ − θ j ).
We now rephrase this condition in a more instructive way. To this end we assume that we have already obtained a solution of the dynamic condition (14). Then we can try to successively assign a phase θ j to every node j in the network. Starting at a node j 0 with an arbitrary phase θ j 0 , we assign the phases of all neighboring nodes j 1 such that sin(θ j 1 −θ j 0 ) = S j 0 ,j 1 . We then proceed in this way through the complete network to assign the phase of an arbitrary node j n , where (j 0 , j 1 , . . . , j n ) is an arbitrary path form j 0 to j n and we have used a solution of the equation for every edge (j, ).
In general, a given node j n can be reached from j 0 via a multitude of different paths.
To define a unique set of phases that satisfies the geometric condition (17), we must assure that Eq. (18) yields a unique phase regardless of which path is taken from j 0 to j n . This is equivalent to the condition, that the phase differences over every simple cycle (as defined in Definition 2) in the network must add up to an integer multiple of 2π.
where ∆ j, is a solution of equation (19). Furthermore, it is sufficient if (20) is satisfied by the cycles in the cycle basis of the network defined in Definition 3: it will then automatically be satisfied for all simple cycles of the network, since the simple cyles form a vector space.
However, there are two distinct solutions of equation (19) that satisfy ∆ ± j,l ∈ [−π, π). To consider both, we define a partition of the edge set Alternatively, we can define the two sets in terms of the nodal phases as We note that a fixed point where the plus sign is realized for all edges (E − = ∅) is guaranteed to be linearly stable according to corollary 1. We refer to it as normal operation.
To operationalize the geometric condition we now define the winding number (27), following the notation used by Ochab and Gora [32]. with The winding vector is defined as Using the winding number we can reformulate the conditions for the existence of a fixed point and establish a correspondence between the description of a fixed points in terms of nodal phases of edge flows.
Theorem 7. Consider a connected network with power injections P ∈ R N and coupling matrix K ∈ R N ×N . Then the following two statements are equivalent: 1. θ * is a fixed point, i.e., a real solution of equation (4).

F ∈ R L satisfies the dynamic condition (13) and
∈ Z L−N +1 for some partition Proof. We prove the theorem in two parts.
(1) =⇒ (2): If θ * is a fixed point, then the flows F satisfying the dynamical condition (13) as given by (3) are Let's partition the edge set into E + and E − by We note the identity that Combining this identity with the definition of ∆ e in (28) and our chosen set partition (31), results in Combining (33) and (34), we obtain ∆ jk = 2m jk π + (θ k − θ j ), m jk ∈ Z (choosing the + sign for 2m jk π without loss of generality).
Then for any simple cycle c = (v 1 , v 2 , · · · , v l , v 1 ) in the cycle basis B C , the winding number is thus completing the first part of the proof.
This concludes the proof.

C. Geometric frustration
The previous reasoning shows that we can face the following situation: Given an oscillator network characterized by the frequencies P j and the capacity matrix K j, , we can find a solution of the dynamical conditions, such that the flow is conserved at all nodes of the network. Nevertheless, no fixed point exists as these solutions are incompatible to the geometric conditions. In this case we say that phase locking is inhibited due to geometric frustration. We summarize this in a formal definition before giving some examples for the importance of this phenomenon. This definition seems unfamiliar at first glance, but is completely compatible to the common concept of geometric frustration in condensed matter theory [33][34][35]. In that context, a system with multiple state variables (x 1 , x 2 , · · · , x n ) is called geometrically frustrated [36] if there must exist certain pairwise correlations between those variables, and no steady state can exist because all these correlations cannot be satisfied simultaneously.
To further clarify the relation to condensed matter systems, we consider a spin lattice system with anti-ferromagnetic interactions [33] where the state variables are the orientation (up or down) of spins. To minimize the total energy of the system, adjacent spins must be aligned anti-parallel introducing perfect pair correlations. Whether this is possible depends on the geometry or topology of the lattice. It is impossible for triangular lattices, since two adjacent spins being antiparallel means the third one has to be parallel to one of those.
Such lattices are thus called frustrated and do not posses a unique minimum energy state [33,34]. In our case, the correlations are given by equation (21): The phases of two adjacent oscillators j and differ by a fixed value defined by ∆ ± j, as given by (21). A fixed point (θ 1 , . . . , θ N ) must satisfy all the correlations, see Eq. (17), otherwise the network is said to be frustrated. As before geometric frustration depends crucially on the topology of the network as we will show in detail in the following section.
In general, the problem of geometric frustration can be traced back to the fundamental cycles of a lattice or network. In condensed matter physics, frustration is classified by the use of plaquette variables, which reveal whether a cycle of the lattice contains incompatible correlations [34]. In oscillator networks an analog function can be defined for every simple cycle c starting from the left-hand side of equation (20) The geometric condition is satisfied for cycle c if Φ c = 0, whereas the cycle is frustrated for

IV. EXAMPLES AND APPLICATIONS
In this section we discuss the importance of geometric aspects for the fixed points of an oscillator network with different topologies. In particular, we analyze the number of fixed points and show that geometric frustration can inhibit phase locking, which may lead to counter-intuitive phenomena.
A. Trees do not suffer from frustration.
By definition, a tree does not contain any cycle such that the geometric condition (20) does not apply. Therefore, the calculation of an fixed point of the power grid oscillator

B. Multiple solutions in cycle
We now consider the simplest nontrivial topology of a cyclic network with only three nodes and three links with equal strength K. The dynamical condition for the existence of a fixed point then reads and |S j, | ≤ 1. In particular for P j = 0 any solution is a cycle flow (S 1,2 , S 2,3 , S 3,1 ) = S × (1, 1, 1).
Taking into account that there are two possible solutions for the phase difference along each edge as per (21)), and since in order to satisfy the geometric condition (20), the sum of phase differences along the cycle must equal 2mπ for some integer m ∈ Z, we see that all fixed points must satisfy Taking all combinations of either ∆ + or ∆ − and corresponding possible values of m, we see that there are three intersections corresponding to three fixed points. These fixed points are illustrated in Figure 1. This shows that stationary states are generally not unique, not even for the simplest cycle network. In the present case only one of the solutions is dynamically stable, but this is generally not true in larger cycles as we will show in the following.
We now extend the above example to a single cycle with an arbitrary number of nodes with the same power P j ≡ 0. All links have an equal strength K as above. For the sake of notational convenience we label the nodes as 1, 2, . . . , N along the cycle and identify the node 1 with N + 1 and 0 with N . In order to have a non-trivial closed cycle we need N ≥ 3.
The dynamic conditions for a fixed points are then given by We stress that the dynamic conditions have a continuum of solutions, i.e. all values F in the interval [−K, K] are allowed.
The phase difference along the edges (j + 1, j) is given by equation (21), leaving two possible solutions + and −. Choosing the minus sign for at least one edge ( + 1, ) yields In this case one can show that the Hesse matrix M is not positive semi-definite such that the fixed point must be unstable. Restricting ourselves to the dynamically stable states, we find that the phase differences are all equal and given by The geometric condition now yields which can be satisfied only for certain discrete values of F . The geometric condition thus induces a 'quantization' of the phase differences where · denotes the floor function. We note that solutions with (θ j+1 − θ j ) = ±π/2 have jacobian eigenvalues µ k = 0 for all k ∈ {1, . . . , N }. In this case linear stability analysis fails to determine dynamical stability properties (see Khazin und Shnol [27] for details). For two coupled oscillators it is rather easy to see that the fixed point is nonlinearly unstable. In total, we thus find 2 × (N − 1)/4 + 1 different stable stationary states.
This example is very simple but illustrates three important general results. First, there can be multiple stable fixed points in cyclic networks as previously noticed by [32,[37][38][39].
This fact is not fully recognized in power engineering, probably because most authors in this community concentrate of fully connected networks which arise after a Kron reduction [17,37]. Second, the oscillator model (2) allows for stable fixed points with a persistent current around a cycle. Interestingly, theses states are phase locked but not phase ordered in the sense that the phase order parameter [7] vanishes exactly for K > 0. Third, the geometric condition induces the discreteness of the phase differences although the dynamic condition allows for continuous values of cycle flows.

D. Braess' paradox
Here we introduce a special example which illustrates the paradoxical effects of geometric frustration most clearly. We consider the oscillator network depicted in Fig. 2  The dynamic condition for this network reads and |S j+1,j | ≤ 1, identifying node j = 5 with j = 1. For notational convenience, we define the vector The solutions of the linear system of equations (46) span a one-dimensional affine space parametrized by a real number , Evaluating the condition |S j+1,j | ≤ 1 yields a necessary condition for the existence of a fixed point K ≥ P.
For κ = 0 this condition is also sufficient for the existence of a stable fixed point. If the capacity of the upper link increases, κ > 0, geometric frustration inhibits phase locking.
A solution of the dynamical conditions always exists for K ≥ P , but this can become incompatible with the geometric condition. We illustrate this in the stability diagram in Figure 2 (b). A stable fixed point exists only in the parameter region above the white line.
As we see in Figure 2 (b), the minimum K required to maintain steady operation, the critical coupling K c , increases when κ is increased.
To further characterize the long-time behavior of the oscillator network we define ω ∞ as the average phase velocities of all the nodes in the limit of large time Therefore ω ∞ must be zero for steady operation to take place. As expected we find ω ∞ = 0 in the stable parameter region above the white line K > K c and ω ∞ > 0 in the unstable parameter region below the white line K < K c . Remarkably, ω ∞ is largest for small values of κ and, of course, K < K c (κ).
This leads to the paradoxical effect that an increase of local transmission capacity reduces the ability of the network to support a phase locked fixed point. This behavior can also be seen as an example of Braess' paradox [40,41] which has been first predicted for traffic networks [42].

V. MULTISTABILITY AND THE NUMBER OF FIXED POINTS
The previous examples show that there can be a large number of stable fixed points in a cyclic network. In the following we derive conditions for the existence and bounds for the number of stable fixed points depending on the network structure. We start with a deeper analysis of the dynamic condition for arbitrary networks, which is a necessary prerequisite for the existence of a stable fixed point. Then we turn to the geometric condition and derive bounds for the number for fixed points. The arguments depend heavily on the network structure such that we will start with trees and simple cycles before we turn to more complex topologies.

A. The dynamic condition
We first analyze whether the dynamic condition (13) admits a solution. The problem reduces to the Multi-source multi-sink maximum flow problem, which can be solved by a variety of different algorithms [43,44].
So let G = (V, E) be a connected graph with N nodes and L edges. Each edge is assigned a capacity given by K 1 , . . . , K L and each node has an in-or outflux given by P 1 , . . . , P N .
We define an extended graph G = (V , E ) by adding two vertices s and t to the vertex set, and adding directed links connecting s (t) to all nodes with positive (negative) power injection: The capacity of the newly added links is infinite. Then one finds the theorem: Alternatively, a sufficient condition for the existence of a solution can be found from dividing the graph into parts: Let (V 1 , V 2 ) an arbitrary partition of V and E(V 1 , V 2 ) the cutset induced by this partition. Then we definē We note that have assumed that j P j = 0, w.l.o.g, such that we always hatP 1 +P 2 = 0.
then there exists a solution of the dynamic condition (13).
Proof. The idea is to prove that: a solution of the dynamic condition (13a) and (13b).
Reversing arguments then yields the theorem. It remains to show that the statement "⇒" is true.
Let F be a solution of (13a). According to our assumption the set of overloaded edges is not empty. Now consider one overloaded edge e 0 = (u, v) ∈ E ov . We assume w.l.o.g that the flow is from u to v, i.e. that F u→v > K uv > 0. We define the weighted directed network G(V, E) with E = E\e 0 and coupling constants We determine the maximum flow pattern ∆F e , e ∈ E with the value ∆F max from u to v in the network G. According to the max-flow min-cut theorem there is a partition (V 1 , V 2 ) with u ∈ V 1 and v ∈ V 2 and the associated cutset E(V 1 , V 2 ) such that Now consider the flow pattern F defined by This is a new solution of the condition (13a). Basically we have rerouted the maximum possible flow from the edge e 0 = (u, v) to alternative paths from u to v. Furthermore we , which is a cut of the original graph G.
We now have to distinguish two cases: Case 1: The maximum flow value ∆F max < F e 0 − K e 0 . Then the edge e 0 is still overloaded, i.e. we have F e 0 > K e 0 . Summing up equations (13a) over the nodes in V 1 and V 2 , respectively, yieldsP However, we know that F e 0 > K e 0 and F e = K e for all other e ∈ E(V 1 , V 2 ) such that and the statement "⇒" follows.
does no longer contain e 0 , i.e. E ov ∈ E ov \ e 0 . However, this set cannot be empty as we have assumed that there is no solution of (13a) satisfying (13b). Then we can just restart the procedure, selecting an edge e 1 ∈ E ov and finding a max. flow between its adjacent vectors.
Finally we must arrive at the case 1 for which the statement "⇒" follows.

B. Tree network
In the previous section we have argued that multistability arises due to the possibility of cycle flows. In a tree there are no cycles and thus no multistability and we obtain the following result. Whether the fixed points exist or not can then be decided solely on the basis of the dynamical codition (11) resp. using theorem 10.
Proof. By definition, a tree has L = N − 1 edges such that the space of solutions of the linear system (14a) has dimension L − N + 1 = 0. That is, there is either zero or exactly one unique solution for the flows F j, . In the first case no fixed point exists. In the latter case there are 2 possible values for the phase difference for each of edge given by equation (21). Hence there are 2 L = 2 N −1 fixed points. Choosing the +-sign in equation (21) where K red , +1 < 0 and M 1 and M 2 are defined as in equation (6) Due to the structure of the matrix M 1 we have M 1 (1, . . . , 1) Thus, the Hesse matrix A is not positive semi-definite, i.e. it has at least one negative eigenvalue and the fixed point is unstable (cf. lemma 1).

C. Cycle flows and winding vector
In the following we want to operationalize the theorem (7), which characterizes fixed points in terms of the flows and winding numbers, to derive strict bounds for the number of fixed points in a network. Restricting ourselves to normal operation (E − = ∅) and using the decomposition (16), the definition of the winding numbers (27) reads using equation (16). The concept of winding numbers is particularly useful when they are unique. If we can find upper and lower bounds for the c , then we can simply count the number of solutions ∈ Z L−N +1 to obtain the number of fixed points. Uniqueness is rigorously shown for planar graphs in the following lemma.
A graph is called planar if it can be drawn in the plane without any edge crossings. Such a drawing is called a plane graph or a planar embedding of the graph and any cycle that surrounds a region without any edges is called a face of the graph [30]. For the sake of simplicity we adopt the convention that for plane graphs the cycle basis B C is built up from the faces in the following. Notably, many power grids and other supply networks are actually planar. Crossing of power lines are not forbidden a priori, but are rare.
Lemma 3. For a planar network, let θ and θ be two fixed points satisfying the "normal operation" criterion (9). If (θ) = (θ ) then both fixed points are the same, i.e. the phases differ only by an additive constant: In other words, no two different fixed points in planar networks can have identical winding defining two cycle flow vectors f and f . We write (f ) and (f ) in short-hand notation for the corresponding winding vectors. We show that (f ) = (f ) implies f = f and thus F = F . As we are assuming normal operation, we can reconstruct the phases via (18) and thus find θ = θ + c(1, 1, · · · , 1) T as we need to show.
So assume that (f ) = (f ) and f c = f c for at least one cycle c. We show that this leads to a contradiction such that the lemma follows. First consider the case that f c − f c is the same for all cycles, f c − f c = ∆f = 0 for all cycles c ∈ B C . Then choose a cycle k at the boundary. If ∆f > 0 we find k (f ) > k (f ) and if ∆f < 0 we find k (f ) < k (f ). This contradict the assumption and the lemma follows.
Otherwise, choose a cycle for which the quantity f c − f c is the largest. We can find a cycle k such that f k − f k > f n − f n for at least one cycle n adjacent to k.
or, equivalently, f k − f n > f k − f n for at least one cycle n adjacent to k.
We now exploit that any edge belongs to at most two cycles, according to Mac Lane's planarity criterion [45]. Choosing an edge e which is part of both the cycles k and n, we have z k e z k e = 1 and z k e z n e = −1. For all other cycles = k, n, we have z e = 0. Thus we find For every other edge e in the cycle k we find by the same procedure (using (72)) that Substituting these two inequalities in the definition (67) and using that the arcsin is monotonically increasing and point-symmetric about the origin such that arcsin(z k e x) = z k e arcsin(x), we find This contradicts our contrary assumption, which concludes the proof.

D. Simple cycles
For networks containing a single cycle, tight upper and lower bounds can be obtained for the number of fixed points satisfying cos(θ * i − θ * j ) > 0 for all edges (i, j). These states correspond to the normal operation of a power grid and are guaranteed to be stable by corollary 1. Other stable steady states can exist in particular at the border of the stable parameter region [26]. We label the nodes as 1, 2, . . . , N along the cycle, fixing the direction of counting in the counter-clockwise direction and identify the node 1 with N + 1 and 0 with N . Likewise we fix the orientation of the edges e ∈ {1, . . . , L} such that F e > 0 describes a counter-clockwise flow and F e < 0 a clockwise flow.
We will first calculate the exact number of fixed points counting the number of different allowed winding numbers. However, this result depends on one particular solution of the dynamic conditions (11), thereby limiting its applicability. We therefore also derive lower and upper bounds for the number of fixed points in terms of a few simple characteristics of the grid, in particular the maximum partial net power. These bounds do not depend on any particular solution of the dynamical condition.
Remark 11. For any ring network R N with N nodes, the cycle flow vector defined in (2) and the winding vector defined in (29) naturally reduce into single numbers. We refer to them as cycle flow f c and winding number c , following [32]. These two quantities will be crucial in establishing the results in the rest of this section.
Theorem 12. For a ring network R N , the number of normal operation fixed point (denoted by N ) is given by where · denotes the floor function and · denotes the ceiling function. F (0) ij is one particular solution to the dynamic condition (11) and Proof. Suppose we have one fixed point θ 0 with the flows F We emphasize that f c cannot be equal to f max c or f min c , because otherwise one edge would be fully loaded with cos(θ i − θ j ) = 0, contradicting our assumption.
Second, all fixed points have to satisfy the geometric condition (cf. theorem 7) Since we restrict ourselves to normal operation, the winding number for a single cycle reads Using the bound for the cycle flow strength (81), and the fact that the arcsin is a monotonically increasing function, we find that the winding number is also bounded by  any fragment F i,j , the partial net powerP ij is defined as and the maximal partial net power is defined as This concept is illustrated in Fig. 3 Furthermore we define the maximum and minimum transmission capacities Lemma 4. For any ring fragment F i,j , the partial net power is equal to the net outwards flow:P Proof. According to lemma 12 the number of fixed points N is given by We make use of the fact that the arcsin function is bounded, arcsin (x) ∈ [−π/2, +π/2], such This proves the first part (93) of the corollary. To proof the second part, we start from Now one can obtain upper and lower bounds for all terms in the sum using the trigonometric relation which holds for all −1 ≤ y ≤ x ≤ 1. This yields 1 2π where we define ∆f c = f max c − f min c . Furthermore, this quantity can be bounded as such that the fraction in equation (99) becomes In a similar way we find Substituting these bounds into equation (99) yields which combined with (97) completes the proof.
In particular, ring networks R N with N ≤ 4 do not have multiple stable fixed points. Ring network R N with N ≥ 7 nodes will have multiple stable fixed points (N ≥ 2) if These relations can be proven by simply evaluating the bounds in corollary 3.

E. Complex networks
Obtaining bounds for the number of fixed points is hard in general, as we cannot simply decompose a network into single cycles, unless no two cycles of a network share an edge.
The diffculty arises because cycle flows in two faces sharing one or more edges can cancel or enhance each other. That is why one cannot simply multiply the bounds for number of fixed points for each cycle to obtain a bound for the total number of fixed points for a network.
We demonstrate this using two examples.

Two cycle flows destroying each other
First, we show that even if all single cycles support (multiple) fixed points in case there are isolated, the full network may not have a single fixed point at all. This is illustrated in Figure 5 for a network consisting of just two cycles. The network motifs shown in panels So we have seen that getting a lower bound for number of fixed points of a general network is hard, as multiplying lower bounds for each cycle in a cycle basis does not yield a valid lower bound. Next we will show why obtaining a good upper bound is also hard.
Consider any of the two identical single loop networks in Figure 6. It consists of one generator and one consumer, generating and consuming 2.1K power respectively. Each edge has capacity K. None of the two single loop networks have any fixed point: simply because the network does not have enough capacity to transport the 2.1K amount of power from the generator to teh consumer. However, when those two are fused together, two cycle flows emerge and a stable fixed point with winding vector ω = (1, −1) comes into existence. This should not come as a surprise: fusing two cycles in this case ended up with an alternate pathway for the powerflow being created.

F. Planar networks
Although obtaining estimates for number of fixed points for general topologies is quite difficult, we now show that for planar topologies, it is possible to derive some analytical insights.

Upper bound
Theorem 14. Consider a finite planar network. Choose the faces of the graph as the cycle basis B C . Then the number of normal operation fixed points is bounded from above by where N c is the number of nodes in the cycle c.
Proof. In a planar network no two different fixed points can have the same winding vector (see lemma 3) such that we can just count the different allowed winding vectors. For each fundamental cycle c ∈ B C we have because −π/2 < ∆ e < +π/2 in normal operation. Counting the number of different possible values of the winding numbers 1 , . . . , L−N +1 respecting these upper and lower bounds yields the result.

Asymptotic behaviour
We have shown in subsection V E that it is not straightforward to obtain bounds for the number of fixed points N in complex networks, unlike simple cycles. However, in the limit of N 1, we can nevertheless derive the scaling behaviour for N . F min 12 ≤ min Then the possible cycle flows in each cycle are bounded inside a convex polygon D described by Then, for K 1, n 1 1, n 2 1, the number of fixed points converges to the area of the image set of D under the mapping .
N ≈ where the Jacobian J( ) for the change of variable can be computed from the expression for in (67), which yields Taking the limits Redefining In the last line we use the symmetry in f 1 , f 2 , both in the integrand, and the domain of integration. We can simplify even further, by using the following change of variables in the second integral: We note that the domain remains the same after this change of variable, and the determinant of the Jacobian det(J) = −1. This allows the simplification N ≈ (n 1 n 2 + (n 1 + n 2 )n 12 ) 1 to finally yield this scaling result lim n 1 ,n 2 →∞ N =γ (n 1 n 2 + (n 1 + n 2 )n 12 ) (122) To evaluate the accuracy of (122), we apply it to two special cases. First, we consider networks with n = n 1 = n 2 , n 12 = 1, i.e. two identical cycles sharing only one single edge.
Second, we consider networks with n = n 1 = n 2 , n 12 = n, i.e. two identical cycles sharing half of their edges. In this case (122) becomes N (n, n, n) ≈ 3γn 2 .
We see in Figure 8 that in both these two cases, the scaling relations are quite accurate even for not very large network sizes, such as n = 50. where each cycle has 2n edges and they share n edges between them.

General planar graphs
The scaling results for two-cycle networks can be extended to general planar graphs, to this end we define a few quantities.
Definition 15 (Loopy dual graph). Given a planar graph G(V, E) and an embedding, we choose a cycle basis B C consisting of the faces of the embedding. The loopy dual graph G dual (G) is a undirected multigraph whose vertex set is equal to B C . Its edge set E is as follows. For each edge e ∈ E, if it is shared between two cycles c 1 and c 2 , then an edge between c and c is added to E . If e is at the boundary and belongs to only one cycle c, then a self loop is added at node c. We illustrate this definition in Figure 9.
Now, consider a planar graph and an arbitrarily chosen fixed point with flows F e . Let's denote byL loopy the loopy laplacian of its meta graph, as defined in Definition 1. Then equation (118) generalizes to

VI. UNSTABLE FIXED POINTS
In principle, we can generalize the cycle flow approach to find fixed points which do not satify the normal operation condition, too. These fixed points are typically linearly unstable (cf. the discussion in [26]). However, most of the results on the number of fixed points cannot be generalized to this case. As an instructive example, consider again a homgeneous ring as in section IV C. We label the nodes as 1, 2, . . . , N along the cycle and assume that N is an integer multiple of 4. All nodes have a vanishing power injection P j ≡ 0 and all links have an equal strength K as before. Then it is easy to see that θ * = (0, δ, π, π + δ, 2π, 2π + δ, 3π, . . .) T is a fixed point of the dynamics for each value of δ ∈ [0, π). This class of fixed points represents a pure cycle flow, for all edges (j, j + 1). The winding number is = N/4 independent of the value of δ and the edges belong alternately to E + and E − : This simple example shows that two main assumptions made for the normal operation fixed points (where E − = ∅) do not longer hold: First, the set of fixed points is no longer discrete. Instead we find a continuum of solutions parametrized by the real number δ.
Second, different fixed points yield the same winding number. Thus we cannot obtain the number of fixed points by counting winding numbers in general.

VII. CALCULATING FIXED POINTS
The cycle flow approach yields a convenient method to calculate multiple fixed points for oscillator networks. Generally, it is hard to make sure that a numerical algorithm yields all solutions for a nonlinear algebraic equation. However, we have shown that the winding numbers are unique at least for normal operation fixed points in planar networks. Thus we can scan the allowed values of the winding numbers and try to find a corresponding solution.
This can be done by starting from an arbitrary solution of the dynamical condition and adding cycle flows until we obtain the desired winding numbers.
In particular, we can calculate all fixed points in normal operation for a planar network using the following algorithm: 1. Find a solution F (0) of the dyanmic condition.
where the winding numbers are given by equation (27). Dropping the assumption of a normal operation, we loose the the guarantee of uniqueness as discussed in section VI. Nevertheless the method can be readily adapted to find most of the unstable fixed points, at least if the number |E − | is small. This can be very useful, as a systematic calculation of such fixed points is generally not straightforward. The results can be applied, among other things, to assess the global stability of a stable fixed point by analyzing the stability boundary [46,47] or the stability in the presence of stochastic fluctuations [48]. In particular, we must add another The output of this algorithm is shown in Figure 10 for a small test network and |E − | ≤ 2.
For this small network we have only L − N + 1 = 3 fundamental cycles of which one is decoupled. Hence we can graphically check that we hvae obtained all fixed points.

VIII. DISCUSSION
Oscillator networks are ubiquitous in nature and technology. A lot of research in statistical physics starting from Kuramoto's seminal work [2] has been devoted to the onset of partial synchronization in large networks. However, in some applications global synchronization is required. In particular in electrical power grids all generators have to run with exactly the same frequency and have to be strictly phase-locked to enable stable power flows to the customers. A desynchronization generally has catastrophic consequences. An example is provided by the European power blackout in November 2006. Following a shutdown of one transmission line and unsuccessful attempts to restore stable operation, the European grid fragmented in three mutually asynchronous clusters [23]. In the end more than 10 million customers were cut from the power supply.
In this article we have analyzed the existence of stable fixed points in finite oscillator networks. The main methodological advancement is to split the calculation into two parts: First, we calculate the flows which satisfy the continuity equation at all nodes. Then we single out the the specific solution which leads to consistent phases of the oscillators. We thus move the focus of the calculation from the nodes (phases) to the edges (flows) and cycles. An immediate consequence is that several fixed points can co-exist which differ by cycle flows. Thus oscillator networks are in general multistable. Several aspects of multistability have been previously discussed in the literature. Multistability in isolated rings was discussed in [32]. The limits (93) were derived and the basins of attraction of the different fixed points was studied numerically. The case of a densely connected graph was analyzed by Taylor in [37]. He was able to show that there is at most one stable fixed point if the node degree is at least 0.9395 × (N − 1). Mehta et al. investigate multistability in complex networks numerically using a similar approach as the present paper [39]. They argue that the number of fixed points scales with the number of cycles as each cycle can accommodate cycle flows. While this is valid for many graphs, there are counterexamples ( Figure 5). Delabays et al [49] have recently reported their treatment of multistability using cycle flows. They have extended the upper bounds for fixed points in single rings by [32] to include also those stable fixed points with phase differences along edges > π/2. They have also derived upper bounds [50] for number of fixed points for planar graphs in case of uniform power injections at all nodes. Xi et al. [51] have numerically shown that spatical heterogeneity of power injections P j reduces the number of fixed points, which fits with our analytical result in Corollary 3. Intriguingly, they have also found that in heterogeneous ring topologies, nonlinear stability of fixed points decrease with the ring size N .
In this work, we have obtained a lower bound for the number of fixed points, and thereby provided a sufficient condition for existence of multistability. Furthermore, we have shown that the length of the cycles N c and the homogeneityP max are equally important for multistability and thereby arrived at tighter bounds for the number of fixed points than Ochab and Gora [32] and Delabays et al [49]. Morever, we have derived scaling laws at the limit of infinite transmission strengths that are much tighter than the upper bound results previously reported. We have shown the derived scaling behaviour to match numerically computed exact results for moderately sized networks.
Interestingly, Our results show that a previous highly recognized result presented by Jadbabaie et al in [52] is incorrect. The authors claim that for any network of Kuramoto oscillators with different natural frequencies, there exists a K u such that for Thus Banach's contraction theorem cannot be applied, which spoils the proof.

IX. CONCLUSION
In summary, taking cycle flows as a basis of flow patterns we analyzed existence and stability of phase locked states in networks of Kuramoto oscillators and second order phase oscillators modeling the phase dynamics of electric power grids. We demonstrated that such systems exhibit multistability. Intriguingly, multistability prevails even under conditions where unique stable operating points were believed to exist in both a power engineering text book and a major complex network reference on Kuramoto oscillators [52,53]. For classes of network topologies, we have established necessary and sufficient conditions for multistability, derived lower and upper bounds for the number of fixed points. We explained why generalizing those bounds for arbitrary topologies is hard. Nevertheless, we have derived asymptotic scaling laws at large loop limit that has been found to match closely numerically