Non-linear network dynamics with consensus-dissensus bifurcation

We study a non-linear dynamical system on networks inspired by the pitchfork bifurcation normal form. The system has several interesting interpretations: as an interconnection of several pitchfork systems, a gradient dynamical system and the dominating behaviour of a general class of non-linear dynamical systems. The equilibrium behaviour of the system exhibits a global bifurcation with respect to the system parameter, with a transition from a single constant stationary state to a large range of possible stationary states. Our main result classifies the stability of (a subset of) these stationary states in terms of the effective resistances of the underlying graph; this classification clearly discerns the influence of the specific topology in which the local pitchfork systems are interconnected. We further describe exact solutions for graphs with external equitable partitions and characterize the basins of attraction on tree graphs. Our technical analysis is supplemented by a study of the system on a number of prototypical networks: tree graphs, complete graphs and barbell graphs. We describe a number of qualitative properties of the dynamics on these networks, with promising modeling consequences.


I. INTRODUCTION
Network dynamics are widely used as a natural way to model complex processes taking place in systems of interacting components. Within this framework, time-varying states are assigned to the nodes of a network and evolve according to interaction rules defined between neighbouring nodes. Sufficiently simple for theoretical investigations, the resulting dynamics may yet exhibit complex emergent behaviour of the global network state, making them suitable to model various real-world systems. Moreover, the interplay between the underlying network structure and the rich phenomenology of dynamics taking place on it make network dynamics a powerful tool to better understand and characterize the network itself. Some well-known examples of network dynamics include random walks [1,2], epidemic spreading [3], synchronization of oscillator systems [4,5], consensus dynamics and voter models [6,7] and power grids [8]. An overview of these applications and many other examples can be found in [9][10][11]. In this article, we propose a new non-linear dynamical system inspired by the pitchfork bifurcation normal form. Our choice of dynamical equations is supported by a number of different interpretations. We find that the system can be seen as (i) a set of interacting (1D) pitchfork systems, (ii) a gradient dynamical system for a potential composed of double-well potentials over the links of the network and finally (iii) as the dominating behaviour of a general class of non-linear dynamics with odd coupling functions. Qualitatively, the main property of the system is that it exhibits a bifurcation in the possible stationary states. In the first parameter regime, our system is essentially diffusive and evolves to a unique, uniform stationary state. In the second parameter regime the coupling function is a mixed attractive/repulsive force and the equilibrium is characterized by a large number of stationary states. We find an explicit description for (a subset of) these stationary states and analyse their stability using linear stability analysis. Our main technical result classifies the stability of these stationary states in terms of the effective resistance of certain links. The effective resistance is a central concept in graph theory with links to random walks [1,12], distance functions and embeddings [13][14][15], spectral sparsification [16] and many more. Its appearance as a determinant for (in)stability in our non-linear dynamical system is very surprising and at the same time a perfect example of the rich interplay between structure and function in network dynamics. Furthermore, analytical results are found for the basins of attraction (on tree graphs) of the stationary states, and an exact solution of the system is derived for certain types of graphs which include graphs with external equitable partitions. The latter result adds to a long list of interesting observations of dynamics on graphs with (external) equitable partitions and related symmetries [17][18][19][20][21][22]. Our technical analysis is supplemented by a detailed description of the system on complete and barbell graphs. On the complete graph, we find that a subset of the stable stationary states determine a balanced bipartition of the graph with each group corresponding to one of two existing state values, and neither group being too dominant (hence balanced). On the barbell graph, a similar balanced bipartition is observed within each of the complete components but with a non-zero difference between the average states of both components. We discuss how these observations might be interpreted in the framework of opinion dynamics.
Our choice to focus on a specific dynamical system is restrictive in various ways and our results only pertain to a small corner of the theory of non-linear dynamics on networks as a consequence. In a follow-up on the present work however, we found that our results generalize to a much broader class of non-linear systems [23], suggesting a potential wider relevance. Other works on this subject, notably the results of Golubitsky et al. [22,24] and Nijholt et al. [25,26] describe and characterize general classes of systems whose dynamics are constrained by a given underlying structure. Their results allow to determine which dynamical features (e.g synchronization conditions, bifurcations) are robust (generic) with respect to the network structure; in other words it details which features can be explained purely from the network structure irrespective of the specific choice of coupling functions. Our contributions are no attempt at such generality on the system level, but instead aim at developing a qualitative understanding of non-linear dynamics on graphs, starting from a basic 'toy' system and describing its interesting properties, with a focus on the influence of the network structure on these properties. A second relevant line of research is the recent work by Franci et al. [27] and Bizyaeva et al. [28] which study decision-making in (multi-option) opinion dynamics. They formulate opinion dynamics in a fully generalized setting, and show -independent of further system models -that this setting can exhibit a variety of rich nonlinear dynamical features such as consensus-dissensus bifurcations and opinion cascades. The model analysed in our article fits in the framework of Franci et al. as a particular two-option opinion dynamical system and consequently, certain features such as the global consensus-dissensus bifurcation and the observations in Section VI can be explained in this context. However, our contributions are complementary to those made in [27,28] as our particular model choice allows us to derive many other specific and interesting results, in particular related to the stability of stationary states and exact solutions in the presence of external equitable partitions. Furthermore, the main results of Franci et al. follow from a so-called equivariant analysis of the system which deduces properties of the system, starting form its symmetries. Our results (and those in [23]) follow from an algebraic and graph theoretic analysis instead and are valid for more general network structures as a result. A third body of related work is the well-developed field of coupled oscillator systems [4,5,21], where many similar questions are studied for non-linear (oscillator) systems on networks. In Section VII we briefly discuss the setup of coupled oscillator systems and highlight a particular result from [29] which closely relates to our stability result, Theorem 1.
The rest of this paper is organized as follows. Our dynamical system is introduced in Section II together with a number of interpretations of the system. Section III introduces the notion of stable and unstable stationary states, and describes the stability results for our system. Section IV describes some cases where the system equations can be solved exactly, and Section V deals with the characterization of basins of attraction. In Section VI finally, system (1) is studied on a number of prototypical networks with a focus on the qualitative behaviour of the solutions. A related result about synchronization in coupled oscillators is described in Section VII, and the article is concluded in Section VIII with a summary of the results and perspectives for future research.

II. THE NON-LINEAR SYSTEM
We will study a dynamical system defined by a set of non-linear differential equations that determine the evolution of a dynamical state x(t). This state is defined on a graph G where each of the N nodes has a corresponding state value x i (t) ∈ R which together make up the system state as x(t) = (x 1 (t), . . . , x N (t)). The dynamics of x(t) are determined at the node level by a non-linear coupling function between neighbouring nodes. For a node i with neighbours j ∼ i, the dynamics are described by where r is a scalar parameter, called the system parameter. Since the states are coupled via their differences, the average state value does not affect the dynamics and the state space of system (1) is thus equal to X = R N /1, i.e. with any two states x and y equivalent if x − y is constant for all nodes. In other words, the dynamics is translation invariant. When considering a specified initial condition x(0) = x 0 , we will also write the solution of system (1) as x(t, x 0 ). There are various ways to interpret the node states.
In the setting of (linear) consensus dynamics, as used frequently in the robotics and control community, the state variable x i (t) represents a real-valued parameter or measurement of an agent in a physical system and the goal is to coordinate these variables globally by following some local dynamics [6], similar to our system (1). In the setting of opinion dynamics [27,30] on the other hand, the node states z i (t) ∈ I in an interval (usually I = [0, 1]) reflect the commitment of an agent in the network to an option/belief A (z i = 0) or to an alternative B (z i = 1) instead [31]; the state dynamics then model the (social) processes by which agents update their opinions or beliefs. As shown in [28], there is a mapping of (forward invariant and bounded) dynamics in R N , i.e. system (1) on X, to opinion dynamics with state space Z. In this context, the system parameter r is sometimes interpreted as a measure of social attention or susceptibility to social influence. Our system can thus be seen as a non-linear generalization of consensus dynamics (see also [32]) or can be mapped onto a two-option opinion dynamics model; the further derivations in this article will be independent of these interpretations.
In what follows, we show how our system appears naturally in three different settings. Apart from suggesting different motivations for the study of our system, each perspective comes with a set of tools and results that will be used in our further analysis.
A. Three perspectives on the dynamical system

Pitchfork bifurcation normal form
The definition of system (1) is inspired by the socalled pitchfork bifurcation dynamical system. This 1dimensional system with state x(t) ∈ R is given by the non-linear differential equation where we will further also use the short-hand notation p(x) = rx − x 3 for the pitchfork function. System (2) is the prototypical form (i.e. normal form) for dynamical systems that exhibit a bifurcation from a single stationary state to three distinct stationary states [11]. This bifurcation occurs between a single stable stationary state x ⋆ = 0 when r < 0, and two stable states x ⋆ = ± √ r and one unstable state x ⋆ = 0 when r > 0. Figure 1 shows the solutions of the pitchfork system (see also Section IV) and illustrates the characteristic bifurcation diagram to which the system thanks its name.  1. (a) Exact solutions for the pitchfork system dx/dt = p(x) (as described in Section IV, Appendix C) for r = 1 and a range of initial conditions x0. These solutions illustrate the stable (x ⋆ = ± √ r) and unstable (x ⋆ = 0) stationary states for positive r. (b) Together with the stable stationary state x ⋆ = 0 for negative r, these solutions determine the characteristic, eponymous bifurcation diagram of system (1) The system studied in this article thus consists of a pitchfork bifurcation system for the state difference (x i − x j ) over each of the links, with interactions coming from the shared variables of links with common nodes. Unsurprisingly, the larger interconnected system exhibits more complex behaviour than each of the smaller systems added together. In particular, our main result Theorem 1 highlights that the stable stationary states of the interconnected system can differ greatly depending on the way in which the links are interconnected. Another way to see that system (1) is closely related to the pitchfork bifurcation normal form is by introducing the link variable y ℓ = (x i − x j ) for all links ℓ = (i, j) ∈ L with an orientation i j fixed by taking the difference (x i − x j ). The dynamics can then be rewritten as where two links m ∼ ℓ meet if they share a common node, and where the sign σ(ℓ, m) = ±1 of the interaction term depends on the relative orientation of the links; the matrix with entries σ(ℓ, m) for adjacent links and zero otherwise is also referred to as the edge adjacency matrix.

Dominating behaviour of odd coupling functions
System (1) is a specific example of a more general class of non-linear dynamical systems on a graph: An important property of this class of systems is that the average state x 1 N x i is always a conserved quantity [33] for the dynamics. If we furthermore assume f to be analytic, the dominating behaviour for systems of the form (4) around the consensus state can be studied by looking at the Taylor expansion of f around (x i −x j ) = 0 as A first-order approximation retrieves a simple, linear diffusion process. For the third-order approximation on the other hand, we see that by introducing the parameter r = −6( df dx d 3 f dx 3 )| x=0 and rescaling time as t ′ = −(6 d 3 f dx 3 )t we retrieve system (1). In other words, the analysis of system (1) is indicative for a general class of non-linear systems with odd coupling functions in the near-consensus regime [34]. In [32], systems of the form (4) are considered within the general problem of non-linear consensus and called relative non-linear flow. They are studied alongside absolute non-linear flow, of the form dx i /dt = (f (x i ) − f (x j )) and disagreement non-linear flow, of the form dx i /dt = f ( (x i − x j )). While some general results are found for the latter two, the discussion of relative flow systems in [32] is limited to the description of a number of small systems.

Gradient dynamical system
System (1) also has the strong property that it is a gradient dynamical system. This means that there exists a potential function V : X → R on the state space, such that the state dynamics are given by the negative gradient of this potential. For system (1), the potential takes the form from which the dynamics are retrieved as dx/dt = −∇V (x) with the gradient operator Interestingly, we see that the potential V (x) in (5) is composed of a separate potential term for each of the links. As illustrated in Figure 2 these terms are equal to a double-well potential, which are minimal at (x i − x j ) = ± √ r separated by a local maximum at x i = x j . As we will see later, the link differences at these local optima also appear as stationary solutions of the system.

FIG. 2.
A double-well potential is a symmetric potential function with two local minima (the 'wells') separated by a local maximum. In the case of our gradient system (1), the potential function V is composed of a double-well potential term for each of the link differences, as where the specific double-well function V dw following from (5) is illustrated above.
An important feature of gradient dynamical systems is that the potential is a decreasing function of time, i.e. the potential satisfiesV ≤ 0 with equality if and only if the system is at a stationary point. This means that system (1) is dissipative for the potential, in contrast with its conservation of the average value˙ x = 0. This feature restricts the possible evolution of a system for a given initial state as V (x(t, x 0 )) ≤ V (x 0 ) must always be satisfied.

III. STATIONARY STATES
Starting from the definition of the dynamics (1), we study a number of different aspects of the system. A first important characterization is the long-term behaviour of the dynamics: starting from some initial state at t = 0, in which states can we expect to observe the system after waiting sufficiently long? This question is answered by studying the stationary states x ⋆ of the system, which are equilibrium states where the system is at rest, i.e. characterized by dx ⋆ /dt = 0. Since these states are the only points in which the potential does not strictly decrease, the system is guaranteed to evolve to a stationary state eventually. From a practical perspective, the stronger notion of asymptotically stable stationary states is interesting. These are states for which the system is in a robust equilibrium, i.e. in the case of some perturbation ǫ, the state x ⋆ + ǫ will evolve back to x ⋆ . We start by determining (a subset of) the stationary states of system (1), followed by an analysis of their stability. A direct translation of the stationarity condition yields the following characterization of stationary states: Generally, finding a stationary state x ⋆ thus involves solving a (potentially large) system of cubic equations. However, the possible solutions for x ⋆ differ greatly depending on the value of r. When r < 0, only a single stationary state is possible: the consensus stationary state where each node state equals the same constant value x ⋆ i = c, equivalent to 0 ∈ X in the state space. In other words, for r < 0 system (1) is a (non-linear) form of diffusion or consensus dynamics. For r > 0 on the other hand, the equations for x ⋆ can have many different solutions. Consider the case where a pair of linked nodes i and j have a state difference equal to √ r. Then the coupling function between i and j will vanish, since p( √ r) = 0, and the same happens when this difference is equal to − √ r or 0. As a consequence, if the difference over all links equals one of these three values, all of the coupling functions will be inactive and the system will be in equilibrium. In other words, any state of the form is stationary. Since the global equilibrium in these states originates from a local equilibrium (balance) for each of the links, we refer to solutions of this type as detailed-balance stationary states [35]. Links with a zero difference will also be called consensus links and links with a ± √ r difference dissensus links. We note that the terms consensus and dissensus are often used to describe the global state of a system instead, while two neighbouring nodes having the same (different) state is then called (dis)agreement [6,27]. Our choice to use the consensus/dissensus terminology at the level of single links follows from our interpretation of system (1) as an interconnected collection of smaller systems, each of which can be in consensus, dissensus or another state (in non-detailed balance states). In certain cases these two uses of consensus and dissensus coincide, see e.g. Section VI B.
We recall from Section II A (and Figure 2) that a dissensus link corresponds to a minimum for the doublewell potential and a consensus link to a local maximum. This means that local stationary states are composed to form global stationary states. In principle, 3 L possible detailed-balance solutions exist, with each link independently taking one of the possible differences. When the graph contains cycles however, these differences must be consistent across each cycle which reduces the number of possible solutions, down to a minimum of just 2 N possible detailed-balance states for the (maximally cyclic) complete graph (see Section VI B). From the perspective of gradient dynamics, the potential of detailed-balance stationary states can be expressed compactly in terms of the number of dissensus links ℓ as In other words, the higher the fraction of dissensus links in a stationary state, the lower the corresponding potential. We use this result in Section V when describing the basins of attraction on tree graphs.
When solving equations (6) directly or simulating the system, other stationary states can be found. In the case of highly symmetrical graphs for instance, tools from equivariant dynamics can be used to find explicit descriptions of stationary states [36]. Generally, whenever a graph has cycles (i.e. it is not a tree as in Section VI A) solutions may exist which are not detailed-balance states. Such states are difficult to describe in general and might even be degenerate. On the 3-cycle graph for instance, all states x ⋆ on the circle In what follows, we focus exclusively on detailed-balance stationary states as they admit an explicit description. However, we have found in a follow-up investigation that the results on detailed-balance stationary states in the following (sub)sections fully generalize to all stationary states [23].

A. Stability conditions
As mentioned earlier, the stationary states of a dynamical system do not always correspond to a robust equilibrium. To characterize the stability of a state, we study how a perturbed state x ⋆ + ǫ evolves and in particular, whether it converges back to x ⋆ or not. To this end, we assume the perturbation ǫ to be sufficiently small such that the dynamics are determined by the linearized system around x ⋆ as where J(x ⋆ ) is the Jacobian of system (1) at x ⋆ . If this Jacobian is (positive) negative definite, it implies directly that the stationary state x ⋆ is linearly (un)stable, which in turn implies (no) asymptotic stability. As we only consider stability criteria following from the linearized system (10), we will further omit 'linear' and simply write stable and unstable. If the Jacobian is semi-definite instead, the linearization is not sufficient to determine the stability of x ⋆ and other techniques are required.
Restricting our analysis to the detailed-balance stationary states (7), we can further simplify the linearized system (9) and characterize certain stable stationary states.
Here, we present a first stability result for system (1): On any graph, the following states (if they exist) are stable stationary states of system (1) For r > 0 the full consensus state is unstable.
Proof: See Section III B. Proposition 1 only gives a rough picture of the stability of system (1), but it does illustrate clearly how the local dynamics are manifest in the global dynamics: the fact that dissensus is stable for each link (1-dimensional pitchfork system) locally while consensus is unstable, is observed globally as well. In the following section we refine this picture and show that the interconnected system also supports different types of stable states which are not simply inherited from the local dynamics. In particular, we find that for r > 0 in the range between full consensus (and thus instability) and full dissensus (and thus stability) there may be stable mixed states with both types of links present. As consensus links cannot exist stably for the local dynamics (see pitchfork dynamics, Section II A), the existence of these stable mixed states is necessarily a feature of the system as a whole.

B. Laplacian form of the linearized system
Before continuing our analysis, we introduce some more information about the graph (network) on which the system takes place. By G = (N , L) we will denote a graph with a set of N nodes N , and a set of L links L ⊆ N × N that connect pairs of distinct nodes, written as i ∼ j or (i, j) ∈ L. We assume the graph to be finite and connected, i.e. with at least one path connecting each pair of nodes. Any graph G has a corresponding (N × N ) Laplacian matrix Q, with entries defined by where the degree d i of a node i equals the number of neighbours of i in G. The Laplacian matrix is just one among several matrix representations, but it is known to have close relations to many important graph properties [37][38][39] and appears in the formulation of diffusion processes on a given graph [2].
In the case of our system, the Laplacian matrix appears when calculating the Jacobian J(x ⋆ ) of the system around some detailed-balance stationary state x ⋆ . From their definition in (9), we find that the entries of the Jacobian equal r} denote the links of G over which there is consensus, respectively dissensus in x ⋆ . Correspondingly, we define the Laplacians Q = and Q = of the subgraphs of G restricted to the consensus, respectively dissensus links [40]; these matrices satisfy Q = Q = + Q = since L = L = ∪ L = holds. This subgraph decomposition allows the Jacobian to be expressed as follows: Lemma 1 The Jacobian J(x ⋆ ) of system (1) at a detailed-balance solution x ⋆ with consensus and dissensus links L = and L = can be written as Proof: Identity (11) follows directly from the elementwise expression (10) and the definition of Laplacian matrices. Lemma 1 implies that the stability problem for detailedbalance stationary states comes down to characterizing the spectrum of a difference of Laplacian matrices and, in particular, the positivity/negativity of its spectrum. An important result about the Laplacian matrix of a connected graph is that it is positive semidefinite (i.e. non-negative eigenvalues) with a single zero eigenvalue corresponding to the constant eigenvector [37,38]. As the state space X is orthogonal to the constant vector (by conservation of average), the Laplacian is thus effectively positive definite. This observation leads to a direct proof of the stability result from Section III B. Proof of Proposition 1: If x ⋆ is the full consensus stationary state, we have that L = = L and thus J(x ⋆ ) = rQ, which is positive definite if r > 0 (x ⋆ is unstable) and negative definite when r < 0 (x ⋆ is stable). If x ⋆ is the full dissensus state on the other hand, we find L = = L such that J(x ⋆ ) = −2rQ which is negative definite for r > 0 (x ⋆ is stable). While Proposition 1 is a direct result of the relation (11) between the Jacobian and the Laplacian matrix of the graph on which the system takes place, the result does not depend on the specific structure of G but only on the properties of the Laplacian matrix in general. The specific structure will play an important role in the case of mixed stationary states.

C. Stability via effective resistances
Somewhat surprisingly, the stability of mixed stationary states can be characterized in terms of the effective resistance. The effective resistance was originally defined in the context of electrical circuit theory, but has found its way into graph theory through various applications such as random walks [12], distance functions [13], graph embeddings [14] and, more recently, graph sparsification [16]. The effective resistance ω ij between a pair of nodes i and j in a graph G can be defined as with Q † the pseudoinverse of the Laplacian of G. For more intuition into the effective resistance, we refer the readers to [8,41], where expression (12) is derived starting from the electrical circuit equations. One of the important properties of the effective resistance is that it determines a metric between the nodes of G [13], where a small effective resistance between a pair of nodes indicates that these nodes are essentially close and 'well connected', while a large effective resistance indicates the opposite. For instance, for a pair of linked nodes i ∼ j the extreme values for effective resistance correspond to ω ij = 2/N for the complete graph (i.e. very well connected) and ω ij = 1 for a tree graph (i.e. poorly connected).
We can now continue to characterize the stability of detailed-balance stationary states in the r > 0 regime. From Proposition 1 we know that in full dissensus the system is stable while full consensus is unstable. Here, we provide a partial answer to the stability question for mixed detailed-balance states with both consensus and dissensus links. In particular, we consider the case where a single consensus link is added to an otherwise full dissensus state; in this case, the stability depends on which link the consensus takes place: Theorem 1 (single consensus link stability) For system (1) with r > 0 on any graph G, the mixed stationary state x ⋆ with a single consensus link L = = {(i, j)} satisfies Proof: The proof is given in Appendix A and is based on Lemma 1 and a new approach to bound the eigenvalues of a difference of Laplacian matrices. Theorem 1 states that a single consensus link state can be stable, depending on the effective resistance of the consensus link i ∼ j. Importantly, the criteria in Theorem 1 are tight (except for a single point). If the effective resistance of the consensus link is high, i.e. if i and j are not well-connected, the state will not be stable. As mentioned before, an extreme example of the effective resistance is the case of tree graphs, where each pair of nodes has only a single link between them with no other possible paths such that ω ij = 1. Generally, a large effective resistance indicates 'bridge links', i.e. links between nodes which have few (or long) parallel paths between them (see example in Section VI C). Adding more parallel paths between i and j will gradually reduce the effective resistance until ω ij = 2/3 is crossed, at which point the corresponding mixed state turns stable. In other words, bridge-like links with few alternative paths in parallel cannot sustain consensus, while links with many alternative parallel paths can.
The answer to the initial question whether mixed stationary states can be stable is thus yes, with the important condition that the consensus occurs between wellconnected nodes. The proof of Theorem 1 is easily adapted to provide a condition for mixed states with several consensus links: Proposition 2 (mixed stationary state stability) For system (1) with r > 0 on any graph G, the mixed stationary state x ⋆ with consensus links L = satisfies Proof: See Appendix B. While Proposition 2 is applicable to all mixed stationary states, the stability criteria are not tight like the criteria of Theorem 1. Indeed, there are generally many detailed-balance states x ⋆ on a graph which satisfy neither of the criteria (13) and for which Proposition 2 thus does not apply. As discussed in Section VI C, one of the consequences of Theorem 1 and Proposition 2 seems to be that in networks with a community structure, the stable states will generally contain more dissensus links between different communities than within. This would result in a higher similarity of node states within each of the communities, compared to an expected bias between the communities, which is an attractive modeling feature e.g. in the context of social cleavage [42]. Crucially however, the results in Section VI B show that within each of the communities, a certain level of dissensus is still expected to occur -the so-called spontaneous symmetry breaking described in [27] -which can also be explained based on effective resistances in the graph, as shown in [23].
To summarize, we studied the stationary states of system (1) and identified the detailed-balance states (7) as a subset of all possible stationary states. The characterization of the Jacobian matrix around detailedbalance stationary states as a difference of Laplacian matrices (Lemma 1) enables a characterization of the stability in terms of the effective resistance. Most importantly, we find a tight stability condition for states with a single consensus link (Theorem 1) as well as more general, but less tight conditions for any mixed stationary state (Proposition 2). In follow-up work [23], we have found that all these results generalize to the setting of system (4) with any odd coupling function f , and for all stationary states (using a suitable reformulation).

IV. EXACT SOLUTIONS
On certain networks, the stationary states x ⋆ of system (1) can coincide with eigenvectors of the network Laplacian Q. As developed in detail in [43] for contagion dynamics, this allows for an exact solution of the state evolution. Applied to our system, we find the following result: Theorem 2 (Exact solution) If system (1) on a graph G has a stationary state x ⋆ ∈ X which is also a Laplacian eigenvector satisfying Qx ⋆ = µx ⋆ , then the exact solution for initial state x 0 = α 0 x ⋆ and r > 0 is given by In particular, the system will reach the stationary state lim t→∞ x(t, x 0 ) = x ⋆ .
Proof: see Appendix C.
In other words, Theorem 2 states that if the subspace Z ⊂ X spanned by an eigenvector z of the Laplacian matrix contains a stationary state of system (1), then any initial condition in Z allows for an exact solution [44]. Moreover, as x 0 ∈ Z implies that x(t, x 0 ) ∈ Z, the subspace Z is a positive invariant set for the dynamics. For r < 0, solution (14) still holds as long as |α 0 | < 1. The question remains for which graphs there exist stationary states of system (1) which are also Laplacian eigenvectors. In other words, we are looking for graphs for which there exists a state x ⋆ that satisfies for all i. We will further refer to the states that satisfy (15) as eigenstates of our system; regarding our system as a map φ t : x 0 → x(t, x 0 ), we find that φ t (x ⋆ ) = α t x ⋆ for these vectors, similar to the definition of eigenvectors for linear maps. Example: An elementary example of an eigenstate can be found for system (1) on a pair of connected nodes, i.e. G = K 2 . The corresponding Laplacian matrix 1 −1 −1 1 has a single non-constant eigenvector equal to z = (1, −1) T with corresponding eigenvalue µ = 2. Scaling this eigenvector as x ⋆ = √ r/2(1, −1) T yields a detailed-balance stationary state, indicating that x ⋆ is an eigenstate of system (1). Consequently, the system can be solved exactly for K 2 consistent with the fact that we can solve the pitchfork normal form exactly, as shown in Figure 1. In the following subsection we describe how simple examples like this two-node graph can be used as a starting point to construct new examples.

A. Graphs with external equitable partitions
In the study of network dynamics and Laplacian matrices, an important type of graph symmetry are equitable partitions [17,45]. A partition π of a graph divides the nodes of G into K disjoint groups N 1 , . . . , N K ⊆ N and is called an external equitable partition (EEP) if all nodes in a group have the same number of links d km to all other groups, in other words (16) for all i, j ∈ N k =m . If G has an external equitable partition π, its structure at the partition level can be summarized by the quotient graph G π . This quotient graph has node set {1, . . . , K} corresponding to the node groups of G and a set of weighted, directed links − → L that connect node group pairs (k, m) between which there exist links in G, and with link weights d km for the link going from k to m, and d mk for the link going from m to k. Some examples of equitable partitions and quotient graphs are given in Figure 3. For more details on equitable partitions and their relation to dynamical systems, we refer the reader to [17,45]. The concept of external equitable partitions will allow us to construct eigenstates of system (1) on graph G based on eigenstates on its quotient graph G π . Since G π is generally a directed and weighted graph, we generalize the definition of eigenstates to this setting as for all k. In Appendix D we show that if a vector y ⋆ satisfies (17) on G π then the corresponding vector x ⋆ with entries x ⋆ i = y ⋆ k for i ∈ N k will also satisfy (15) on G. As a result, we find that Proposition 3 For a graph G with external equitable partition π, any eigenstate y ⋆ of system (1) on G π has a corresponding eigenstate x ⋆ on G.
Proof: See Appendix D. Proposition 3 is a powerful tool for constructing examples of graphs with eigenstates. Indeed, starting from a (directed, weighted) graph G with an eigenstate y ⋆ we can construct many examples of graphs G ′ for which G is a quotient graph, i.e. G = G ′π with respect to an EEP π of G ′ , and for which there thus exists an eigenstate x ⋆ . In this construction, any node k in G can be replaced by a set N k of nodes in G ′ which can be interconnected in any desired way, and where the nodes from N k are then given d km links to nodes in N m , which requires that the identity |N k |d km = |N m |d mk holds for all pairs of partitions. This construction and the corresponding relation between eigenstates is illustrated in Figure 3.

FIG. 3. Illustration of external equitable partitions (EEPs)
, quotient graphs and the construction of eigenstates. In the first row, a partition π of the nodes of G (in two colours) is shown. Since each yellow node has one white neighbour and vice versa, this partition is an EEP and the corresponding quotient graph G π with directed, weighted links is given. In the third column, a number of other graphs G ′ with EEPs are give for which G π is again the quotient graph. In the second row, another instance of G, G π and G ′ is given, together with an example of how an eigenstate y ⋆ (satisfying (17)) on the quotient graph can be used to find eigenstates x ⋆ (satisfying (15)) on graphs G ′ with EEPs.
External equitable partitions arise from the notion (16) of equivalence between nodes of a network, which is based on (local) symmetries between the neighbourhoods of the nodes. This relation between symmetry and dynamics closely resembles the perspective on network dynamics developed by Golubitsky, Stewart et al. [22,46] and generalized in [47][48][49]. Their framework for general directed, labelled graphs and non-linear couplings focuses on local symmetries (which carry the structure of a so-called 'groupoid') of the graph, and studies how these symmetries are manifest in dynamics that respect the network structure. Our result that eigenstates y ⋆ on the quotient graph G π can be 'lifted' to eigenstates on the graph G can be directly understood in context of this framework as a relation between (EEP) symmetry and dynamics, as shown in [50]. These results are complementary to other works that study the effect of global symmetries (i.e. graph automorphisms which carry the structure of a 'group') on dynamical properties [21,36].

V. BASINS OF ATTRACTION
Another classical question in (non-linear) dynamics is to determine which initial conditions lead to which stationary states. More specifically, the problem consists of characterizing the basins of attraction W (x ⋆ ) of the stationary states x ⋆ , which are subsets of the state space defined as [51] W where we recall that x(t, x 0 ) is the system state at time t with initial condition x(0) = x 0 . The basins of attraction are positive invariant sets, since for each w ∈ W (x ⋆ ) we have x(t, w) ∈ W (x ⋆ ) for all t > 0, and determine a partition of the state space of system (1) as with any pair of distinct basins disjoint W (x ⋆ )∩W (x ′⋆ ) = ∅. Less formally, expression (18) captures the intuitive fact that for any initialization, system (1) will converge to some stationary state. In general, it can be difficult to determine the basins of attraction for a non-linear system, but in the case of system (1) we can use the additional properties of the dynamics (see e.g. Section II A) to find some partial characterization. For instance, using the non-increasing property of the potential V , we know that the basin of attraction of a stationary point x ⋆ can only contain states of a higher potential, i.e. that V (w) ≥ V (x ⋆ ) for each w ∈ W (x ⋆ ). When system (1) takes place on a tree graph T , we can say even more about the basins of attraction. In Propositions 5 and 6 in Section VI A, we will show that all stationary states are detailed-balance stationary states and that among these states only the full dissensus states are stable. By (8), this means that all stable states on T have the same minimal potential V min = −r 2 L/4. Moreover, from the non-increasing property of the potential we find that there is a critical potential V c −(L − 1)r 2 /4 which determines a transition between states (with potential V > V c ) which in principle could be in the basin of attraction of any stationary state, and states (with potential V < V c ) which can only be in the basin of attraction of a stable state. We find the following characterization of the basins of attraction in the sub-critical regime: Proposition 4 (attraction basins on trees) For system (1) on tree graphs, the state space region with potential lower than the critical potential X| V <Vc can be partitioned into basins of attraction of just the stable stationary states Furthermore, the basins of attraction in this region are given by Proof: See Appendix E.

VI. EXAMPLES AND MODELING OBSERVATIONS
In the previous sections, we focused on the technical analysis of system (1) and in particular on its stationary states. In the rest of the article, we study the system on a number of prototypical networks. We give a qualitative description of the system solutions and suggest how certain properties might be useful when considering our system as a complex systems model.

A. System on loopless networks
On a loopless network, or tree graph T , several of the earlier results are simplified or hold with less restrictions. Firstly, since the graph contains no loops, any assignment of {0, ± √ r} to the links of T is possible; this amounts to 3 N −1 possible detailed-balance stationary states (as any connected tree has N − 1 links). Moreover, condition (6) for stationarity implies (7) for the case of tree graphs, and thus: Proposition 5 For system (1) on tree graphs, all stationary states (satisfying (6)) are detailed-balance stationary states (satisfying (7)).
Proof: See Appendix F. Furthermore, the effective resistance between any pair of nodes of a tree graph is equal to the geodesic distance between these nodes [13] which means that for linked nodes i ∼ j we have ω ij = 1. Consequently, by Theorem 1 and Proposition 2 we find that the stability of system (1) is given by Proposition 6 For system (1) on tree graphs and r > 0, the full dissensus state is stable while all other stationary states are unstable.
Proof: From Proposition 5, the fact that ω ij = 1 for all links i ∼ j in T and Theorem 1, it follows that any stationary state x ⋆ with a consensus link, i.e. with L = non-empty, has max L= ω ij = 1 > 2/3 ⇒ x ⋆ is unstable. The stability of the full dissensus state follows from Proposition 1. One consequence of Proposition 6 is that the proportion of stable stationary states on a tree equals (2/3) N −1 which vanishes exponentially fast for larger trees. As discussed in the previous section in Proposition 4, we also have some information about the basins of attraction for tree graphs.

B. Balanced opinion formation in the complete graph
In the complete graph K N , every node is connected to all (N − 1) other nodes, making it the densest possible graph. Moreover, it means that the graph contains a high level of symmetry in the sense that no two nodes are distinguishable from their connections to other nodes, which greatly simplifies the description of detailed-balance stationary states. Since any three nodes in K N form a triangle, the only stationary states (x ⋆ i , x ⋆ j , x ⋆ k ) these nodes can achieve is some permutation of (0, 0, 0) or (0, 0, ± √ r).
As a consequence, the stationary states at the level of the full graph must be some permutation of the state (0, . . . , 0 with V entries equal to √ r and N − V entries equal to 0. In other words, all stationary states are parametrized by the number V ∈ [0, N ] of √ r-states; since there are N V ways to choose V such nodes, the complete graph has 2 N detailed-balance stationary states (by the binomial theorem). Moreover, if we use the degree of freedom provided by the average state to fix the state of an (arbitrary) reference node to x ⋆ i = 0, the (rescaled) state parameter v = V /N will be related to the average state value by For the stability of the stationary states, we find the following result: Proposition 7 For system (1) on the complete graph K N with r > 0, the detailed-balance stationary states x ⋆ satisfy v ∈ (1/3, 2/3) : This characterization of the (stable) stationary states on the complete graph is very interesting from a modeling point of view. First, in any detailed-balance stationary state the nodes of G are split into two groups with V and (N − V ) nodes, respectively. Within each of these groups the nodes are in consensus, while between the groups there is dissensus. Furthermore, Proposition 7 states that the size of the two groups needs to be balanced in stable states, i.e. the group sizes can differ by at most N/3 and neither of the groups can dominate the full graph. Figure 4 below illustrates the findings of Proposition 7 in the bifurcation diagram of K N . In contrast to loopless graphs, Proposition 7 shows that a non-vanishing proportion of 2/3 of all detailed-balance stationary states are stable on the complete graph.
Assuming the framework of opinion dynamics [7,27,28], where nodes play the role of individuals in a population with states x i (t) recording their preference in the range between a certain opinion A with x i = + √ r/2, or an opposing opinion B if x i = − √ r/2 (rivaling political party, competing product, etc.), we might interpret these results as follows: for r < 0 any difference between initial individual preferences will be disappear from the network until the population reaches a global consensus where all individuals agree. This qualitative behaviour is studied in various contexts like engineering [6] or social sciences [7,52], and can also be reproduced by the simpler diffusion dynamics dx i /dt = j∼i (x j − x i ). For r > 0 on the other hand, an atypical stationary distribution emerges in system (1) where instead of reaching global consensus, the population splits into two groups adhering to different opinions. Moreover, the stability condition v ∈ (1/3, 2/3) guarantees that neither of these groups can be too dominant in the population, i.e. that there is a balanced coexistence of opinions. This qualitative behaviour is observed in real social systems, where it is often called social cleavage or polarization [42]. As noted in [27] it is remarkable that a fully interconnected system with indistinguishable nodes (agents) can exhibit spontaneous symmetry breaking into a state with distinct groups of nodes. The analysis in [27] however explains how this behaviour is expected for a broad class of dynamical systems Finally, we remark that the complete graph should be seen as a prototype for more general 'dense graphs', and that our qualitative description should hold approximately for dense random graphs like, for instance, Erdős-Rényi random graphs with high link probability p, as a result of concentration of measure [53]. Importantly, the above description of the equilibrium behaviour of system (1) on K N does not take any non detailed-balance states into account, for which we might observe very different types of stable states.
An equivalent system consisting of pitchfork bifurcation normal forms on the complete graph was analysed by Aronson et al. in [54] as a model of coupled Josephson junctions. In contrast to our ad-hoc derivation, they make a principled equivariant analysis of the system dynamics and derive the two-group stationary states from the fact that this fully interconnected system has S n × 2 symmetry (permuting nodes × sign change). The same stability conditions as Proposition 7 are noted in [54] based on calculations of the system Jacobian for the pitchfork bifurcation normal form (similar to our proof), however there is no suggestion as to how this stability result might generalize to less symmetrical systems. In a sense, this broader view on the relation between (in)stability and structure is exactly what Lemma 1 and Theorem 1 in the present work and some (stronger) results in [23] build towards. This illustrates the complementarity between the approach and tools of [27,28] (symmetric graphs but general systems) and our approach (general graphs but specific system).

C. Biased opinion formation in the barbell graph
The barbell graph B 2N , illustrated in Figure 5, consists of two complete graphs joined by a single link i ∼ j. Similar to the complete graph, the high number of symmetries in each of the individual complete parts allows for a compact description of the stationary states. In fact, the detailed-balance stationary states on B 2N can be parametrized as the stationary state on two 'independent' complete graphs, i.e. with V A , V B ∈ [0, N ] denoting the number of nodes with a different value from x ⋆ i , x ⋆ j in the two complete graphs respectively. The average state value in the complete graphs is then related to their respective (scaled) parameters as Comparing the average state value between the two components yields When restricting to stable stationary states, we find that the bridge link in the barbell graph has effective resistance ω ij = 1 which implies that this bridge link must have dissensus in all stable states. Furthermore, by Proposition 7 we have that the stationary states in the complete graphs are stable for v A , v B ∈ (1/3, 2/3). Assuming that x ⋆ j > x ⋆ i we then find that the difference between the average stable state values in both components equals In Figure 5 this finding is illustrated in the bifurcation diagram of B 2N , which clearly shows the non-zero dif-ference that exists between both complete components for r > 0. The stable detailed-balance states on the barbell graph are again interesting from a modeling perspective.
In the setting of consensus dynamics, we again find that for r < 0 all individuals converge to a common opinion. For r > 0 on the other hand, a balanced coexistence of opinions will be established within each of the complete graphs separately but, importantly, with a non-zero bias existing on average between the components. In other words, within the dense subgraphs the opinions coexist without either opinion dominating the other, while a difference will exist between the subgraphs. This qualitative behaviour might seem interesting if we think of the barbell graph as a prototypical example of a graph consisting of dense groups of nodes which are sparsely interconnected in between, a structure commonly known as assortative communities. In this setting, one might expect such opinion biases to exist between communities rather than within due to a different level of coordination or communication, and a model similar to our system might help explain the underlying mechanisms.
We remark that the stationary states in this example could also be derived based on the (S N −1 ≀ 2 ) × 2 symmetry of the system (permuting non-bridge nodes within a complete component and/or interchanging components × sign change), where '≀' denotes the wreath product between groups [55]. From the theory of equivariant dynamics [27,36], we know that certain stationary states will be associated to subgroups of these system symmetries. This analysis would not provide any stability information however.

VII. RELATED RESULT: SYNCHRONIZATION OF PHASE OSCILLATORS
A famous example of non-linear dynamics on networks are systems of interacting phase oscillators. The underlying idea is that many natural (herds/shoals of animals, groups of neurons, etc.) and man-made (power grids, electrical oscillators, etc.) systems can be modeled effectively as a population of oscillators which establish some form of synchronization due to interactions [5,11]. The periodic behaviour of a single entity is abstracted as an oscillator whose state θ i (t) ∈ R/2π cycles periodically according to a natural frequency f i . These oscillators are then interconnected in a network G, with a coupling function h driving the phases of adjacent entities to a common value as with the coupling strength µ as a system parameter. The easiest example of a periodic, odd coupling function is the sine function. Similar to how our non-linear system (1) is the 3 rd order Taylor approximation for any odd coupling function f (on R), the sine function can be seen as the 1 st order Fourier expansion for any periodic odd coupling function h (on R/2π). System (21) with h(x) = sin(x) on the complete graph is also known as the Kuramoto model and is widely studied in the context of synchronization, see for instance [4,5,56]. One of the key features that motivates the study of interacting oscillator systems is that the oscillators exhibit synchronization for certain parameter ranges of µ and {f i } on certain graph structures. The onset of various types of synchronization (phase, frequency, chimera, etc.) has been studied extensively in these systems, and is used as a theoretical explanation for observed synchrony in many real-world systems. Here, we mention a specific result about coupled oscillators which is similar to our main result Theorem 1.
A particular notion of synchronization is "frequency synchronization with γ-cohesive phases", which is defined as the (rotating) state θ ⋆ where all oscillators rotate at the same instantaneous frequency dθ ⋆ i /dt = f ⋆ and where the phase differences between adjacent oscillators in the network satisfy (θ ⋆ i − θ ⋆ j ) mod 2π ≤ γ for all i ∼ j. We will call such a state γ-synchronized. In [29] the authors propose to study for which choices of natural frequencies f = (f 1 , . . . , f N ) T this type of synchrony can occur. Their interesting finding is that for many graphs (certain extremal graphs, and dense sets of random graphs) the following criterion is a sufficient condition for system (21) on a graph with Laplacian Q and sinusoidal coupling, to have a γsynchonized state. In particular, this implies the known result that that system (21) with a constant natural frequency f = α1 can have a γ-synchronized state, for any γ. Moreover, if the natural frequencies are equal for all but one pair of connected nodes i and j, which differ by |f i − f j | = c, then the synchrony criterion becomes c ≤ sin(γ)/ω ij , i.e. the difference c is upper-bounded by the inverse of the effective resistance. In other words, starting from the constant frequency distribution for which there is synchrony possible, and changing a single link to have a frequency difference of c, then synchrony is conserved depending on the effective resistance of the respective link. More specifically, a small (large) effective resistance will admit a large (small) phase difference. While the setting of [29] is very different, this result is reminiscent of Theorem 1, and a further investigation of this similarity might be worthwhile.

VIII. CONCLUSION
In this paper, we have introduced and studied a nonlinear dynamical system on networks inspired by the pitchfork bifurcation. Our analysis is motivated by different interpretations of the system as a collection of interdependent pitchfork systems, a gradient dynamical system for a potential composed of interacting doublewell potentials and finally as the dominating behaviour for more general non-linear systems with odd coupling functions. In a certain sense system (1) is the 'simplest' dynamical system of a broad class of non-linear systems (e.g. with general odd coupling functions f in (4), or general symmetric potentials V ). The choice to study equations (1) specifically is thus the outcome of a wish to implement a model with more complexity than simple linear models, while wielding Occam's razor. Our technical analysis mainly focused on the equilibrium behaviour of the system. The bifurcation from a single stationary state to a myriad of possible stationary states and in particular their stability, provides a clear picture of how the simple local dynamical rule in our system gives rise to interesting global phenomena. Specifically, as a main technical result (Theorem 1) we found stability conditions that depend on the full structure of the network, as captured by their dependence on the effective resistance. Our further analysis of the system includes the identification of exact solutions for certain graphs, which include graphs with external equitable partitions, and the description of basins of attraction for loopless networks. Finally, we looked at the system on a number of prototypical graphs and describe some interesting qualitative behaviour of the solutions. On the complete and barbell graph, our results suggest an interpretation of the system as an opinion dynamic model: in one parameter regime the system is driven to a global consensus state, while the stable states in the other regime are characterized by a balanced bipartition of opinions (states) in dense components and with an overall non-zero bias between sparsely connected components, that grows as the system goes deeper in the parameter regime. These results support system (1) as a rich model for complex systems allowing to identify unexpected bridges between network properties, like the effective resistance, and dynamical ones, which could trigger future advances in the more general study of non-linear systems on networks.
By Courants Minmax principle we can write for any (a, b) ∈ L = and with positive normalizing constant c = Q † (e a − e b ) −2 . Since (B1) holds for any consensus link, it is clear that max L= ω ab > 2/3 ⇒ λ max > 0 proving the first stability condition.
For the upper-bound, we proceed similar to Appendix A and let y be the normalised vector such that y T Jy = λ max and define P = y T Qy. We can then write where the maximum in the third line is achieved by The upper-bound (B2) shows that if L= ω ij < 2/3 ⇒ λ max < 0 proving the second stability condition.

Appendix C: Proof of exact solution
We prove Theorem 2 about exact solutions of system (1) for certain initial conditions. In line with the approach of [43], we will show that if at some time t, the state is of the form x(t) = α(t)x ⋆ with x ⋆ a stationary state of system (1) parallel to an eigenvector of the Laplacian matrix with eigenvalue µ (i.e. x ⋆ is an eigenstate), the dynamic equations (1) simplify to an equation for α(t). This proves that the subspace spanned by the vector (eigenstate) x ⋆ is a positive invariant subspace, as for any initial state x 0 = α 0 x ⋆ the solution is of the form x(t, x 0 ) = α(t)x ⋆ . Secondly, we show that the timedependent coefficient α(t) can be solved exactly as the solution of a 1-dimensional Bernoulli differential equation. Proof: The system equations for x(t) are given by If at some time t the state is of the form x(t) = α(t)x ⋆ , where x ⋆ is an eigenstate of the system satisfying conditions (15), these equations become Since x ⋆ is a stationary state of the dynamics, we have Next, as x ⋆ is an eigenvector of the Laplacian Q, i.e. with This shows that, for all t ′ > t the solution will be of the form x(t ′ , x(t)) = α(t ′ )x ⋆ and thus that {cx ⋆ for some c ∈ R} is a positive invariant set for the dynamics [58]. Furthermore, the equation for α(t) is a 1-dimensional Bernoulli differential equation which can be solved by introducing β = α −2 − 1 with dβ/dt = −2α −3 dα/dt. This yields a linear differential equation dβ dt = −2rµβ with solution β(t) = β 0 e −2rµt . Introducing the initial condition β 0 = α −2 0 −1 and changing the variable back to α(t) we find the solution which proves Theorem 2.
Appendix D: Proof of Proposition 3 Let y ⋆ be an eigenstate of the quotient graph G π , and x a (N × 1) vector with entries x i = y ⋆ k when i ∈ N k , defined on the nodes of G. For any k and i ∈ N k we can then write j∼i (x i − x j ) = where we used definition (16) of the link weight d km . This illustrates that x is an eigenvector of Q with eigenvalue µ. Similarly, for any k and i ∈ N k we can write which shows that x is also a stationary state of system (1). Together, these identities show that the vector x constructed from the eigenstate y ⋆ of G π is an eigenstate of G.
are determined by the equations Since leaf node l has degree d l = 1, the only equation where x ⋆ l appears is the balance equation for x ⋆ n . Moreover, introducing (x ⋆ n − x ⋆ l ) ∈ {0, ± √ r} in equation (F2) for x ⋆ n eliminates x ⋆ l from the equations for x ⋆ altogether. From this, it follows that the stationary state x ⋆ can be determined from the stationary states x ⋆′ of the tree graph T ′ = T \{l} as x ⋆ i = x ⋆′ i for all i = l (as they obey the same equations) and with (x ⋆ l − x ⋆ n ) ∈ {0, ± √ r} (by the solution of equation (F1)). By the induction hypothesis, x ⋆′ is the stationary state of system (1) on a K-node tree graph and is thus a detailed-balance state. From this follows that x ⋆ will be a detailed-balance stationary state as well.
Appendix G: Stability on the complete graph We prove Proposition 7 which states that stationary states x ⋆ with V ∈ (N/3, 2N/3) are stable by explicitly calculating the spectrum of the Jacobian J(x ⋆ ). Proof: The stationary state x ⋆ partitions the set of nodes N of K N into two disjoint sets V with |V| = V and V = N \V such that for all i ∈ V, j ∈ V we have (x ⋆ i − x ⋆ j ) = ± √ r. In other words, all consensus links go between nodes within a same set, while dissensus links go between nodes of a different set. If we order the nodes as V = 1, . . . , V and V = V + 1, . . . , N then the Laplacian matrices can be written as where I * and u * denote the identity matrix and all-one vector of dimensions indicated by * , and with P V = I V − u V u T V /V the projector on the space orthogonal to u V and similarly for N − V . The Jacobian of system (1) at x ⋆ can then be calculated as J = r(Q = − 2Q = ). We will show that J has four types of eigenvectors and, correspondingly, four types of eigenvalues. Type 1: Any vector of the form z = (z V , 0 N −V ) T with z V a V -dimensional vector that satisfies z T V u V = 0 and z V = 1, will give Jz = r (Q = z − 2Q = z) = r(V − 2(N − V ))z = 3r(V − 2/3N )z which shows that z is an eigenvector of J with eigenvalue 3r(V − 2/3N ). Since we can choose a basis of V − 1 vectors orthogonal to u V which all are of the form of z, the Jacobian will have V − 1 eigenvalues equal to 3r(V − 2/3N ). where superscripts denote the multiplicity of the eigenvalues. The first zero eigenvalue corresponds to the constant eigenvector u which is orthogonal to the state space X and thus does not influence the stability of state x ⋆ . The second zero eigenvalue is always negative for r > 0. Finally, if V ∈ (N/3, 2N/3) then all other eigenvalues are negative and if V / ∈ [N/3, 2N/3] they are positive, which proves the (in)stability of the stationary states x ⋆ with respect to V . When V ∈ {N/3, 2N/3}, the other eigenvalues become zero and the linearization method does not provide the necessary information to determine the stability of the corresponding states.