Nonlinear Network Dynamics with Consensus–Dissensus Bifurcation

Abstract

We study a nonlinear 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 nonlinear 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 modelling consequences.

1 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 makes network dynamics a powerful tool to better understand and characterize the network itself. Some well-known examples of network dynamics include random walks (Lovász 1993; Van Mieghem 2014), epidemic spreading (Pastor-Satorras et al. 2015), synchronization of oscillator systems (Strogatz 2000; Arenas et al. 2008), consensus dynamics and voter models (Olfati-Saber and Murray 2003; Castellano et al. 2009) and power grids (Dörfler et al. 2018). An overview of these applications and many other examples can be found in Barrat et al. (2008), Porter and Gleeson (2016) and Strogatz (2018).

In this article, we propose a new nonlinear 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 nonlinear 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 (Lovász 1993; Doyle and Snell 1984), distance functions and embeddings (Klein and Randić 1993; Fiedler 2011; Devriendt and Van Mieghem 2019), spectral sparsification (Spielman and Teng 2011) and many more. Its appearance as a determinant for (in)stability in our nonlinear 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 (Schaub et al. 2016; Pecora et al. 2014; Bonaccorsi et al. 2015; Devriendt and Van Mieghem 2017; Ashwin and Swift 1992; Golubitsky and Stewart 2006).

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 nonzero 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 nonlinear 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 nonlinear systems (Homs-Dones et al. 2020), suggesting a potential wider relevance. Other works on this subject, notably the results of Golubitsky and Stewart (2006), Gandhi et al. (2020) and Nijholt (2018), Nijholt et al. (2019), 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 nonlinear 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. (2020) and Bizyaeva et al. (2020) 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 Sect. 6 can be explained in this context. However, our contributions are complementary to those made in Franci et al. (2020) and Bizyaeva et al. (2020) 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 from its symmetries. Our results (and those in Homs-Dones et al. 2020) 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 (Strogatz 2000; Ashwin and Swift 1992; Arenas et al. 2008), where many similar questions are studied for nonlinear (oscillator) systems on networks. In Sect. 7, we briefly discuss the setup of coupled oscillator systems and highlight a particular result from Dörfler et al. (2013) which closely relates to our stability result, Theorem 1.

The rest of this paper is organized as follows. Our dynamical system is introduced in Sect. 2 together with a number of interpretations of the system. Section 3 introduces the notion of stable and unstable stationary states, and describes the stability results for our system. Section 4 describes some cases where the system equations can be solved exactly, and Sect. 5 deals with the characterization of basins of attraction. In Sect. 6 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 Sect. 7, and the article is concluded in Sect. 8 with a summary of the results and perspectives for future research.

2 The Nonlinear System

We will study a dynamical system defined by a set of nonlinear differential equations that determine the evolution of a dynamical state $\mathbf {x}(t)$. This state is defined on a graph G where each of the N nodes has a corresponding state value $x_i(t)\in \mathbb {R}$ which together make up the system state as $\mathbf {x}(t)=(x_1(t),\dots ,x_N(t))$. The dynamics of $\mathbf {x}(t)$ are determined at the node level by a nonlinear coupling function between neighbouring nodes. For a node i with neighbours $j\sim i$, the dynamics are described by

$$\begin{aligned} \frac{\mathrm{d}x_i}{\mathrm{d}t} = \sum _{j\sim i} r(x_i-x_j) - (x_i-x_j)^3 \end{aligned}$$

(1)

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=\mathbb {R}^N/\mathbf {1}$, i.e. with any two states $\mathbf {x}$ and $\mathbf {y}$ equivalent if $\mathbf {x}-\mathbf {y}$ is constant for all nodes. In other words, the dynamics is translation invariant. When considering a specified initial condition $\mathbf {x}(0)=\mathbf {x}_0$, we will also write the solution of system (1) as $\mathbf {x}(t,\mathbf {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 (Olfati-Saber and Murray 2003), similar to our system (1). In the setting of opinion dynamics (Degroot 1974; Franci et al. 2020) on the other hand, the node states $z_i(t)\in \mathcal {I}$ in an interval (usually $\mathcal {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¹; the state dynamics then model the (social) processes by which agents update their opinions or beliefs. As shown in Bizyaeva et al. (2020), there is a mapping of (forward invariant and bounded) dynamics in $\mathbb {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 nonlinear generalization of consensus dynamics (see also Srivastava et al. (2011)) 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.

2.1 Three Perspectives on the Dynamical System

2.1.1 Pitchfork Bifurcation Normal Form

The definition of system (1) is inspired by the so-called pitchfork bifurcation dynamical system. This 1-dimensional system with state $x(t)\in \mathbb {R}$ is given by the nonlinear differential equation

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t} = rx - x^3 \text { with parameter }r, \end{aligned}$$

(2)

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 (Strogatz 2018). This bifurcation occurs between a single stable stationary state $x^\star =0$ when $r<0$, and two stable states $x^\star =\pm \sqrt{r}$ and one unstable state $x^\star =0$ when $r>0$. Figure 1 shows the solutions of the pitchfork system (see also Sect. 4) and illustrates the characteristic bifurcation diagram to which the system thanks its name.

Fig. 1

a Exact solutions for the pitchfork system $\mathrm{d}x/\mathrm{d}t=p(x)$ (as described in Sect. 4, “Appendix C”) for $r=1$ and a range of initial conditions $x_0$. These solutions illustrate the stable $(x^\star =\pm \sqrt{r})$ and unstable $(x^\star =0)$ stationary states for positive r. b Together with the stable stationary state $x^\star =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_{{\ell }}=(x_i-x_j)$ for all links $\ell =(i,j)\in \mathcal {L}$ with an orientation $i\rightsquigarrow j$ fixed by taking the difference $(x_i-x_j)$. The dynamics can then be rewritten as

$$\begin{aligned} \frac{\mathrm{d} y_{\ell }}{\mathrm{d}t} = 2p(y_{\ell }) + \sum _{{m}\sim {\ell }}\sigma ({\ell },{m})p(y_{m}) \end{aligned}$$

(3)

where two links ${m}\sim {\ell }$ meet if they share a common node, and where the sign $\sigma ({\ell },m)=\pm 1$ of the interaction term depends on the relative orientation of the links; the matrix with entries $\sigma (\ell ,m)$ for adjacent links and zero otherwise is also referred to as the edge adjacency matrix.

2.1.2 Dominating Behaviour of Odd Coupling Functions

System (1) is a specific example of a more general class of nonlinear dynamical systems on a graph:

$$\begin{aligned} \frac{\mathrm{d}x_i}{\mathrm{d}t} = \sum _{j\sim i}f(x_i-x_j) \text { with odd function }f. \end{aligned}$$

(4)

An important property of this class of systems is that the average state $\langle \mathbf {x}\rangle \triangleq \tfrac{1}{N}\sum x_i$ is always a conserved quantity² 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

$$\begin{aligned} \frac{\mathrm{d}x_i}{\mathrm{d}t} = \sum _{j\sim i}\frac{\mathrm{d}f}{\mathrm{d}x}\Big \vert _{0}(x_i-x_j) + \frac{1}{6}\frac{\mathrm{d}^3f}{\mathrm{d}x^3}\Big \vert _{0}(x_i-x_j)^3 + O((x_i-x_j)^5). \end{aligned}$$

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(\tfrac{\mathrm{d}f}{\mathrm{d}x}\big /\tfrac{\mathrm{d}^3f}{\mathrm{d}x^3})\vert _{\mathbf {x}=0}$ and rescaling time as $t' = -(6\big /\tfrac{\mathrm{d}^3f}{\mathrm{d}x^3})t$ we retrieve system (1). In other words, the analysis of system (1) is indicative for a general class of nonlinear systems with odd coupling functions in the near-consensus regime.³

In Srivastava et al. (2011), systems of the form (4) are considered within the general problem of nonlinear consensus and called relative nonlinear flow. They are studied alongside absolute nonlinear flow, of the form $\mathrm{d}x_i/\mathrm{d}t = \sum \left( f(x_i)-f(x_j)\right)$ and disagreement nonlinear flow, of the form $\mathrm{d}x_i/\mathrm{d}t = f(\sum (x_i-x_j))$. While some general results are found for the latter two, the discussion of relative flow systems in Srivastava et al. (2011) is limited to the description of a number of small systems.

2.1.3 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\rightarrow \mathbb {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

$$\begin{aligned} V(\mathbf {x}) = \frac{1}{4}\sum _{i\sim j}(x_i-x_j)^2((x_i-x_j)^2-2r), \end{aligned}$$

(5)

from which the dynamics are retrieved as $\mathrm{d}\mathbf {x}/\mathrm{d}t = -\nabla V(\mathbf {x})$ with the gradient operator $\nabla =\left( \frac{\partial }{\partial x_1},\frac{\partial }{\partial x_2},\dots ,\frac{\partial }{\partial x_N}\right) ^T$. Interestingly, we see that the potential $V(\mathbf {x})$ in (5) is composed of a separate potential term for each of the links. As illustrated in Fig. 2, these terms are equal to a double-well potential, which are minimal at $(x_i-x_j)=\pm \sqrt{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 $V(\mathbf {x})=\sum _{j\sim i}V_{\text {d}w}(x_i-x_j)$ where the specific double-well function $V_{\text {d}w}$ 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 satisfies $\dot{V}\le 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 $\dot{\langle \mathbf {x}\rangle }=0$. This feature restricts the possible evolution of a system for a given initial state as $V(\mathbf {x}(t,x_0))\le V(\mathbf {x}_0)$ must always be satisfied.

3 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 $\mathbf {x}^\star$ of the system, which are equilibrium states where the system is at rest, i.e. characterized by $\mathrm{d}\mathbf {x}^\star /\mathrm{d}t=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 $\varvec{\epsilon }$, the state $\mathbf {x}^\star +\varvec{\epsilon }$ will evolve back to $\mathbf {x}^\star$. 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:

$$\begin{aligned} \frac{\mathrm{d}\mathbf {x}^\star }{\mathrm{d}t}=0 \Leftrightarrow \sum _{j\sim i} r(x^\star _i - x^\star _j) - (x^\star _i - x^\star _j)^3=0,~\forall i. \end{aligned}$$

(6)

Generally, finding a stationary state $\mathbf {x}^\star$ thus involves solving a (potentially large) system of cubic equations. However, the possible solutions for $\mathbf {x}^\star$ 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}^\star _i = c$, equivalent to $\mathbf {0}\in X$ in the state space. In other words, for $r<0$ system (1) is a (nonlinear) form of diffusion or consensus dynamics. For $r>0$ on the other hand, the equations for $\mathbf {x}^\star$ can have many different solutions. Consider the case where a pair of linked nodes i and j have a state difference equal to $\sqrt{r}$. Then the coupling function between i and j will vanish, since $p(\sqrt{r})=0$, and the same happens when this difference is equal to $-\sqrt{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

$$\begin{aligned} \mathbf {x}^\star \text { with } (x_i^\star - x_j^\star )\in \left\{ 0,\pm \sqrt{r}\right\} \quad \text {for all } i\sim j \end{aligned}$$

(7)

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.⁴ Links with a zero difference will also be called consensus links and links with a $\pm \sqrt{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 (Olfati-Saber and Murray 2003; Franci et al. 2020). 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. Sect. 6.2.

We recall from Sect. 2.1 (and Fig. 2) that a dissensus link corresponds to a minimum for the double-well 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 Sect. 6.2).

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 $\ell$ as

$$\begin{aligned} V(\mathbf {x}^\star ) = -\tfrac{1}{4}r^2\ell . \end{aligned}$$

(8)

In other words, the higher the fraction of dissensus links in a stationary state, the lower the corresponding potential. We use this result in Sect. 5 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 (Golubitsky and Stewart 2015). Generally, whenever a graph has cycles (i.e. it is not a tree as in Sect. 6.1) 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 $\mathbf {x}^\star$ on the circle

$$\begin{aligned} \left\{ \mathbf {x} \in X : (x_1-\langle \mathbf {x}\rangle )^2+(x_2-\langle \mathbf {x}\rangle )^2+(x_3-\langle \mathbf {x}\rangle )^2 = {2r/3}\right\} \end{aligned}$$

in X are stationary. 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 (Homs-Dones et al. 2020).

3.1 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 $\mathbf {x}^\star +\varvec{\epsilon }$ evolves and in particular, whether it converges back to $\mathbf {x}^\star$ or not. To this end, we assume the perturbation $\varvec{\epsilon }$ to be sufficiently small such that the dynamics are determined by the linearized system around $\mathbf {x}^\star$ as

$$\begin{aligned} \frac{\mathrm{d}(\mathbf {x}^\star +\varvec{\epsilon })}{\mathrm{d}t} \approx J(\mathbf {x}^\star )\varvec{\epsilon }\quad \text {with }(J(\mathbf {x}^\star ))_{ij} = \frac{\mathrm{d}^2x_i}{\mathrm{d}t\mathrm{d}x_{j}}\Big \vert _{\mathbf {x}=\mathbf {x}^\star } \end{aligned}$$

(9)

where $J(\mathbf {x}^\star )$ is the Jacobian of system (1) at $\mathbf {x}^\star$. If this Jacobian is (positive) negative definite, it implies directly that the stationary state $\mathbf {x}^\star$ 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 $\mathbf {x}^\star$ 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):

Proposition 1

(Full consensus/dissensus stability) On any graph, the following states (if they exist) are stable stationary states of system (1)

$$\begin{aligned} {\left\{ \begin{array}{ll} (r<0):\text { the full consensus state } (x^\star _i-x^\star _{j\sim i})=0\\ (r>0):\text { the full dissensus state }(x^\star _i-x^\star _{j\sim i})=\pm \sqrt{r} \end{array}\right. } \end{aligned}$$

For $r>0$ the full consensus state is unstable.

Proof

See Sect. 3.2. $\square$

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, Sect. 2.1), the existence of these stable mixed states is necessarily a feature of the system as a whole.

3.2 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=(\mathcal {N},\mathcal {L})$ we will denote a graph with a set of N nodes $\mathcal {N}$, and a set of L links $\mathcal {L}\subseteq \mathcal {N}\times \mathcal {N}$ that connect pairs of distinct nodes, written as $i\sim j$ or $(i,j)\in \mathcal {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\times N)$ Laplacian matrix Q, with entries defined by

$$\begin{aligned} (Q)_{ij} \triangleq {\left\{ \begin{array}{ll} d_i &{}\text { if } i=j\\ -1 &{}\text { if } i\sim j\\ 0 &{}\text { otherwise}, \end{array}\right. } \end{aligned}$$

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 (Mohar et al. 1991; Merris 1994; Chung 1997) and appears in the formulation of diffusion processes on a given graph (Van Mieghem 2014).

In the case of our system, the Laplacian matrix appears when calculating the Jacobian $J(\mathbf {x}^\star )$ of the system around some detailed-balance stationary state $\mathbf {x}^\star$. From their definition in (9), we find that the entries of the Jacobian equal

$$\begin{aligned} (J(\mathbf {x}^\star ))_{ij} = {\left\{ \begin{array}{ll} d_ir - 3\sum _{j\sim i}(x^\star _i-x^\star _j)^2 &{}\text { if }i=j\\ -r + 3(x_i^\star -x_j^\star )^2 &{}\text { if }i \sim j\\ 0 &{}\text { otherwise}. \end{array}\right. } \end{aligned}$$

(10)

We let $\mathcal {L}_{=}\triangleq \lbrace i\sim j: x^\star _i-x^\star _j=0\rbrace$ and $\mathcal {L}_{\ne }\triangleq \lbrace i\sim j : x^\star _i-x^\star _j = \pm \sqrt{r} \rbrace$ denote the links of G over which there is consensus, respectively, dissensus in $\mathbf {x}^\star$. Correspondingly, we define the Laplacians $Q_{=}$ and $Q_{\ne }$ of the subgraphs of G restricted to the consensus, respectively, dissensus links⁵; these matrices satisfy $Q=Q_{=}+Q_{\ne }$ since $\mathcal {L}=\mathcal {L}_{=}\cup \mathcal {L}_{\ne }$ holds. This subgraph decomposition allows the Jacobian to be expressed as follows:

Lemma 1

The Jacobian $J(\mathbf {x}^\star )$ of system (1) at a detailed-balance solution $\mathbf {x}^\star$ with consensus and dissensus links $\mathcal {L}_{=}$ and $\mathcal {L}_{\ne }$ can be written as

$$\begin{aligned} J(\mathbf {x}^\star )&= r(Q-3Q_{\ne }) \nonumber \\&= r(3Q_=-2Q)\nonumber \\&= r(Q_{=}-2Q_{\ne }) \end{aligned}$$

(11)

Proof

Identity (11) follows directly from the elementwise expression (10) and the definition of Laplacian matrices. $\square$

Lemma 1 implies that the stability problem for detailed-balance 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 (Mohar et al. 1991; Merris 1994). 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 Sect. 3.2.

Proof of Proposition 1

If $\mathbf {x}^\star$ is the full consensus stationary state, we have that $\mathcal {L}_==\mathcal {L}$ and thus $J(\mathbf {x}^\star ) = rQ$, which is positive definite if $r>0$ ($\mathbf {x}^\star$ is unstable) and negative definite when $r<0$ ($\mathbf {x}^\star$ is stable). If $\mathbf {x}^\star$ is the full dissensus state on the other hand, we find $\mathcal {L}_{\ne }=\mathcal {L}$ such that $J(\mathbf {x}^\star )=-2rQ$ which is negative definite for $r>0$ ($\mathbf {x}^\star$ is stable). $\square$

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.

3.3 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 (Doyle and Snell 1984), distance functions (Klein and Randić 1993), graph embeddings (Fiedler 2011) and, more recently, graph sparsification (Spielman and Teng 2011). The effective resistance $\omega _{ij}$ between a pair of nodes i and j in a graph G can be defined as

$$\begin{aligned} \omega _{ij} = (e_i-e_j)^TQ^\dagger (e_i-e_j), \end{aligned}$$

(12)

with $Q^\dagger$ the pseudoinverse of the Laplacian of G. For more intuition into the effective resistance, we refer the readers to Ghosh et al. (2008) and Dörfler et al. (2018), 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 (Klein and Randić 1993), 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\sim j$ the extreme values for effective resistance correspond to $\omega _{ij}=2/N$ for the complete graph (i.e. very well connected) and $\omega _{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 $\mathbf {x}^\star$ with a single consensus link $\mathcal {L}_{=}=\lbrace (i,j)\rbrace$ satisfies

$$\begin{aligned} {\left\{ \begin{array}{ll} \omega _{ij}<2/3:\mathbf {x}^\star \text { is stable}\\ \omega _{ij}>2/3: \mathbf {x}^\star \text { is unstable} \end{array}\right. } \end{aligned}$$

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. $\square$

Theorem 1 states that a single consensus link state can be stable, depending on the effective resistance of the consensus link $i\sim 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 $\omega _{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 Sect. 6.3). Adding more parallel paths between i and j will gradually reduce the effective resistance until $\omega _{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 well-connected 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 $\mathbf {x}^\star$ with consensus links $\mathcal {L}_{=}$ satisfies

$$\begin{aligned} {\left\{ \begin{array}{ll} \sum _{\mathcal {L}_{=}}\omega _{ij}<2/3: \mathbf {x}^\star is stable \\ \max _{\mathcal {L}_{=}}\omega _{ij}>2/3: \mathbf {x}^\star is unstable \end{array}\right. } \end{aligned}$$

(13)

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 $\mathbf {x}^\star$ on a graph which satisfy neither of the criteria (13) and for which Proposition 2 thus does not apply. As discussed in Sect. 6.3, 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 modelling feature e.g. in the context of social cleavage (Friedkin 2015). Crucially, however, the results in Sect. 6.2 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 Franci et al. (2020)—which can also be explained based on effective resistances in the graph, as shown in Homs-Dones et al. (2020).

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 detailed-balance 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 (Homs-Dones et al. 2020), 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).

4 Exact Solutions

On certain networks, the stationary states $\mathbf {x}^\star$ of system (1) can coincide with eigenvectors of the network Laplacian Q. As developed in detail in Prasse and Van Mieghem (2020) 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 $\mathbf {x}^\star \in X$ which is also a Laplacian eigenvector satisfying $Q\mathbf {x}^\star =\mu \mathbf {x}^\star$, then the exact solution for initial state $\mathbf {x}_0 = \alpha _0\mathbf {x}^\star$ and $r>0$ is given by

$$\begin{aligned} \mathbf {x}(t,\mathbf {x}_0) = \mathbf {x}_0\left( \alpha _0^2 - \left( \alpha _0^2-1\right) e^{-2\mu r t}\right) ^{-1/2} \end{aligned}$$

(14)

In particular, the system will reach the stationary state $\lim _{t\rightarrow \infty } \mathbf {x}(t,\mathbf {x}_0)=\mathbf {x}^\star$.

Proof

See “Appendix C”. $\square$

In other words, Theorem 2 states that if the subspace $Z\subset X$ spanned by an eigenvector $\mathbf {z}$ of the Laplacian matrix contains a stationary state of system (1), then any initial condition in Z allows for an exact solution.⁶ Moreover, as $\mathbf {x}_0\in Z$ implies that $\mathbf {x}(t,\mathbf {x}_0)\in Z$, the subspace Z is a positive invariant set for the dynamics. For $r<0$, solution (14) still holds as long as $\vert \alpha _0\vert <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 $\mathbf {x}^\star$ that satisfies

$$\begin{aligned} \mu {x}^\star _i = \sum _{j\sim i}(x^\star _i-x^\star _j) \text { and } \sum _{j\sim i}r(x^\star _i-x^\star _j) = \sum _{j\sim i}(x^\star _i-x^\star _j)^3 \end{aligned}$$

(15)

for all i. We will further refer to the states that satisfy (15) as eigenstates of our system; regarding our system as a map $\phi _t : \mathbf {x}_0 \mapsto \mathbf {x}(t,\mathbf {x}_0)$, we find that $\phi _t(\mathbf {x}^\star ) = \alpha _t \mathbf {x}^\star$ 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 $\left( {\begin{matrix} 1&{}-1\\ -1&{}1\\ \end{matrix}}\right)$ has a single non-constant eigenvector equal to $\mathbf {z}=(1,-1)^T$ with corresponding eigenvalue $\mu =2$. Scaling this eigenvector as $\mathbf {x}^\star = \sqrt{r}/2(1,-1)^T$ yields a detailed-balance stationary state, indicating that $\mathbf {x}^\star$ 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 Fig. 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.

4.1 Graphs with External Equitable Partitions

In the study of network dynamics and Laplacian matrices, an important type of graph symmetry are equitable partitions (Schaub et al. 2016; O’Clery et al. 2013). A partition $\pi$ of a graph divides the nodes of G into K disjoint groups $\mathcal {N}_1,\dots ,\mathcal {N}_K\subseteq \mathcal {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

$$\begin{aligned} \vert \lbrace v\sim i:v\in \mathcal {N}_m\rbrace \vert =\vert \lbrace v\sim j:v\in \mathcal {N}_m\rbrace \vert \triangleq d_{km} \end{aligned}$$

(16)

for all $i,j\in \mathcal {N}_{k\ne m}$. If G has an external equitable partition $\pi$, its structure at the partition level can be summarized by the quotient graph $G^{\pi }$. This quotient graph has node set $\lbrace 1,\dots ,K\rbrace$ corresponding to the node groups of G and a set of weighted, directed links $\overrightarrow{\mathcal {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 Fig. 3. For more details on equitable partitions and their relation to dynamical systems, we refer the reader to Schaub et al. (2016) and O’Clery et al. (2013).

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^{\pi }$. Since $G^{\pi }$ is generally a directed and weighted graph, we generalize the definition of eigenstates to this setting as

$$\begin{aligned} \sum _{m\sim k}d_{km}(y^\star _k-y^\star _m)&= \mu y^\star _k\quad \text { and}\nonumber \\ \sum _{m\sim k} d_{km}r(y^\star _k-y^\star _m)&= \sum _{m\sim k}d_{km}(y^\star _k-y^\star _m)^3. \end{aligned}$$

(17)

for all k. In “Appendix D”, we show that if a vector $\mathbf {y}^\star$ satisfies (17) on $G^{\pi }$ then the corresponding vector $\mathbf {x}^\star$ with entries $x_i^\star =y_k^\star$ for $i\in \mathcal {N}_k$ will also satisfy (15) on G. As a result, we find that

Proposition 3

For a graph G with external equitable partition $\pi$, any eigenstate $\mathbf {y}^\star$ of system (1) on $G^{\pi }$ has a corresponding eigenstate $\mathbf {x}^\star$ on G.

Proof

See “Appendix D”. $\square$

Proposition 3 is a powerful tool for constructing examples of graphs with eigenstates. Indeed, starting from a (directed, weighted) graph G with an eigenstate $\mathbf {y}^\star$ we can construct many examples of graphs $G'$ for which G is a quotient graph, i.e. $G=G'^{\pi }$ with respect to an EEP $\pi$ of $G'$, and for which there thus exists an eigenstate $\mathbf {x}^\star$. In this construction, any node k in G can be replaced by a set $\mathcal {N}_k$ of nodes in $G'$ which can be interconnected in any desired way, and where the nodes from $\mathcal {N}_k$ are then given $d_{km}$ links to nodes in $\mathcal {N}_m$, which requires that the identity $\vert \mathcal {N}_k\vert d_{km} = \vert \mathcal {N}_m\vert d_{mk}$ holds for all pairs of partitions. This construction and the corresponding relation between eigenstates are illustrated in Fig. 3.

Fig. 3

Illustration of external equitable partitions (EEPs), quotient graphs and the construction of eigenstates. In the first row, a partition $\pi$ 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^{\pi }$ with directed, weighted links is given. In the third column, a number of other graphs $G'$ with EEPs are given for which $G^{\pi }$ is again the quotient graph. In the second row, another instance of $G,G^{\pi }$ and $G'$ is given, together with an example of how an eigenstate $\mathbf {y}^\star$ (satisfying (17)) on the quotient graph can be used to find eigenstates $\mathbf {x}^\star$ (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 and Stewart (2006), Stewart et al. (2003) and generalized in Golubitsky et al. (2005), DeVille and Lerman (2015) and Nijholt et al. (2016). Their framework for general directed, labelled graphs and nonlinear 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^\star$ on the quotient graph $G^{\pi }$ 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 Aguiar and Dias (2018). 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 (Ashwin and Swift 1992; Golubitsky and Stewart 2015).

5 Basins of Attraction

Another classical question in (nonlinear) dynamics is to determine which initial conditions lead to which stationary states. More specifically, the problem consists of characterizing the basins of attraction $W(\mathbf {x}^\star )$ of the stationary states $\mathbf {x}^\star$, which are subsets of the state space defined as⁷

$$\begin{aligned} W(\mathbf {x}^\star ) \triangleq \left\{ \mathbf {x}_0\in X: \lim _{t\rightarrow \infty }\mathbf {x}(t,\mathbf {x}_0)=\mathbf {x}^\star \right\} , \end{aligned}$$

where we recall that $\mathbf {x}(t,\mathbf {x}_0)$ is the system state at time t with initial condition $\mathbf {x}(0)=\mathbf {x}_0$. The basins of attraction are positive invariant sets, since for each $\mathbf {w}\in W(\mathbf {x}^\star )$ we have $\mathbf {x}(t,\mathbf {w})\in W(\mathbf {x}^\star )$ for all $t>0$, and determine a partition of the state space of system (1) as

$$\begin{aligned} X = \bigcup _{\text {stationary }\mathbf {x}^\star }W(\mathbf {x}^\star ) \end{aligned}$$

(18)

with any pair of distinct basins disjoint $W(\mathbf {x}^\star )\cap W(\mathbf {x}'^\star )=\emptyset$. 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 nonlinear system, but in the case of system (1) we can use the additional properties of the dynamics (see e.g. Sect. 2.1) 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 $\mathbf {x}^\star$ can only contain states of a higher potential, i.e. that $V(\mathbf {w})\ge V(\mathbf {x}^\star )$ for each $\mathbf {w}\in W(\mathbf {x}^\star )$.

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 Sect. 6.1, 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^2L/4$. Moreover, from the non-increasing property of the potential we find that there is a critical potential $V_c \triangleq -(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\vert _{V<V_c}$ can be partitioned into basins of attraction of just the stable stationary states

$$\begin{aligned} \left( X = \bigcup _{\text {stable }\mathbf {x}^\star } W(\mathbf {x}^\star )\right) \Bigg \vert _{V<V_c} \end{aligned}$$

(19)

Furthermore, the basins of attraction in this region are given by

$$\begin{aligned} \left( W(\mathbf {x}^\star ) = \Big \lbrace \mathbf {x}\in X : (x_i-x_j)(x^\star _i-x^\star _j)>0 ~\forall i\sim j\right\} \Big )\Big \vert _{V<V_c}. \end{aligned}$$

(20)

Proof

See “Appendix E”. $\square$

6 Examples and Modelling 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.

6.1 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 $\lbrace 0,\pm \sqrt{r}\rbrace$ 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 G”. $\square$

Furthermore, the effective resistance between any pair of nodes of a tree graph is equal to the geodesic distance between these nodes (Klein and Randić 1993) which means that for linked nodes $i\sim j$ we have $\omega _{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 $\omega _{ij}=1$ for all links $i\sim j$ in T and Theorem 1, it follows that any stationary state $\mathbf {x}^\star$ with a consensus link, i.e. with $\mathcal {L}_{=}$ non-empty, has $\max _{\mathcal {L}_{=}}{\omega _{ij}} = 1 >2/3\Rightarrow \mathbf {x}^\star$ is unstable. The stability of the full dissensus state follows from Proposition 1. $\square$

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.

6.2 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^\star ,x_j^\star ,x_k^\star )$ these nodes can achieve is some permutation of

$$\begin{aligned} (0,0,0) \text { or }(0,0,\pm \sqrt{r}). \end{aligned}$$

As a consequence, the stationary states at the level of the full graph must be some permutation of the state

$$\begin{aligned} (\underbrace{0,\dots ,0}_{\#(N-V)},\underbrace{\sqrt{r},\dots ,\sqrt{r}}_{\# V}) \end{aligned}$$

with V entries equal to $\sqrt{r}$ and $N-V$ entries equal to 0. In other words, all stationary states are parametrized by the number $V\in [0,N]$ of $\sqrt{r}$-states; since there are ${{N}\atopwithdelims (){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^\star _i=0$, the (rescaled) state parameter $v=V/N$ will be related to the average state value by $\langle \mathbf {x}^\star \rangle = \pm v \sqrt{r}$. 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 $\mathbf {x}^\star$ satisfy

$$\begin{aligned} {\left\{ \begin{array}{ll} v\in (1/3,2/3) : \mathbf {x}^\star \text { is stable}\\ v\notin [1/3,2/3] : \mathbf {x}^\star \text { is unstable} \end{array}\right. } \end{aligned}$$

Proof

See “Appendix G”. $\square$

This characterization of the (stable) stationary states on the complete graph is very interesting from a modelling 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 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.

Fig. 4

Bifurcation diagram of system (1) on the complete graph $K_N$ with $N=75$ nodes. For $r<0$ the consensus state is the only stable stationary state. For $r>0$, the stable detailed-balance stationary states are parametrized by $v\in (1/3,2/3)$ with corresponding average state value $\langle \mathbf {x}^\star \rangle = \pm v\sqrt{r}$ when fixing an arbitrary node to $x^\star _i=0$. The diagram resembles the bifurcation diagram (Fig. 1) of the pitchfork bifurcation normal form

Assuming the framework of opinion dynamics (Castellano et al. 2009; Franci et al. 2020; Bizyaeva et al. 2020), 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=+\sqrt{r}/2$, or an opposing opinion B if $x_i=-\sqrt{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 (Olfati-Saber and Murray 2003) or social sciences (Castellano et al. 2009; Proskurnikov and Tempo 2017), and can also be reproduced by the simpler diffusion dynamics $\mathrm{d}x_i/\mathrm{d}t=\sum _{j\sim 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\in (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 (Friedkin 2015). As noted in Franci et al. (2020) 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 Franci et al. (2020), 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 (Bandeira 2018). 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. (1991) 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\times \mathbb {Z}_2$ symmetry (permuting nodes $\times$ sign change). The same stability conditions as Proposition 7 are noted in Aronson et al. (1991) 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 Homs-Dones et al. (2020) build towards. This illustrates the complementarity between the approach and tools of Franci et al. (2020) and Bizyaeva et al. (2020) (symmetric graphs but general systems) and our approach (general graphs but specific system).

6.3 Biased Opinion Formation in the Barbell Graph

The barbell graph $B_{2N}$, illustrated in Fig. 5, consists of two complete graphs joined by a single link $i\sim 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\in [0,N]$ denoting the number of nodes with a different value from $x^\star _i,x^\star _j$ in the two complete graphs, respectively. The average state value in the complete graphs is then related to their respective (scaled) parameters as $\langle \mathbf {x}^{\star }\rangle _A = x^\star _i \pm v_A\sqrt{r}$ and similarly for $\langle \mathbf {x}^\star \rangle _B=x^\star _j\pm v_B\sqrt{r}$. Comparing the average state value between the two components yields

$$\begin{aligned} \langle \mathbf {x}^{\star }\rangle _B-\langle \mathbf {x}^{\star }\rangle _A = (x^\star _i-x^\star _j) \pm (v_B\pm v_A)\sqrt{r}. \end{aligned}$$

When restricting to stable stationary states, we find that the bridge link in the barbell graph has effective resistance $\omega _{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\in (1/3,2/3)$. Assuming that $x^\star _j>x^\star _i$ we then find that the difference between the average stable state values in both components equals

$$\begin{aligned} \langle \mathbf {x}^\star \rangle _B - \langle \mathbf {x}^\star \rangle _A = \left[ 1\pm (v_B\pm v_A)\right] \sqrt{r} \text { with }v_A,v_B\in (1/3,2/3). \end{aligned}$$

In Fig. 5, this finding is illustrated in the bifurcation diagram of $B_{2N}$, which clearly shows the nonzero difference that exists between both complete components for $r>0$.

Fig. 5

Difference in average stationary state value of the two complete subgraphs for system (1) on the barbell graph $B_{2N}$ with $2N=48$ nodes. For $r<0$ the consensus state is the only stable stationary state. For $r>0$, the bridge link $i\sim j$ needs to be a dissensus link for stability (by Theorem 1), which allows to fix $x^\star _i=-\sqrt{r}/2$ and $x^\star _j=\sqrt{r}/2$. The stable stationary states are then parametrized by $v_A,v_B\in (1/3,2/3)$ for each of the complete components independently, and the forced bias over the bridge link results in a nonzero average difference of $\sqrt{r}$ between the state averages

The stable detailed-balance states on the barbell graph are again interesting from a modelling 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 nonzero 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}\wr \mathbb {Z}_2)\times \mathbb {Z}_2$ symmetry of the system (permuting non-bridge nodes within a complete component and/or interchanging components $\times$ sign change), where ‘$\wr$’ denotes the wreath product between groups (Wells 1976). From the theory of equivariant dynamics (Golubitsky and Stewart 2015; Franci et al. 2020), we know that certain stationary states will be associated with subgroups of these system symmetries. This analysis would not provide any stability information, however.

7 Related Result: Synchronization of Phase Oscillators

A famous example of nonlinear 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 modelled effectively as a population of oscillators which establish some form of synchronization due to interactions (Strogatz 2018; Arenas et al. 2008). The periodic behaviour of a single entity is abstracted as an oscillator whose state $\theta _i(t)\in \mathbb {R}/2\pi$ 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

$$\begin{aligned} \frac{\mathrm{d}\theta _i}{\mathrm{d}t} = f_i + \mu \sum _{j\sim i}h(\theta _i-\theta _j) \end{aligned}$$

(21)

with the coupling strength $\mu$ as a system parameter. The easiest example of a periodic, odd coupling function is the sine function. Similar to how our nonlinear system (1) is the $3^{\text {rd}}$-order Taylor approximation for any odd coupling function f (on $\mathbb {R}$), the sine function can be seen as the $1^{\text {st}}$ order Fourier expansion for any periodic odd coupling function h (on $\mathbb {R}/2\pi$). 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 Strogatz (2000), Arenas et al. (2008) and Dörfler and Bullo (2014).

One of the key features that motivates the study of interacting oscillator systems is that the oscillators exhibit synchronization for certain parameter ranges of $\mu$ and $\lbrace f_i\rbrace$ 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 $\gamma$-cohesive phases”, which is defined as the (rotating) state $\varvec{\theta }^\star$ where all oscillators rotate at the same instantaneous frequency $\mathrm{d}\theta _i^\star /\mathrm{d}t=f^\star$ and where the phase differences between adjacent oscillators in the network satisfy $(\theta ^\star _i-\theta ^\star _j){\text {mod}} 2\pi \le \gamma$ for all $i\sim j$. We will call such a state $\gamma$-synchronized. In Dörfler et al. (2013) the authors propose to study for which choices of natural frequencies $\mathbf {f}=(f_1,\dots ,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

$$\begin{aligned} \max _{i\sim j}\left| (e_i-e_j)Q^\dagger \mathbf {f}\right| \le \sin (\gamma ) \end{aligned}$$

is a sufficient condition for system (21) on a graph with Laplacian Q and sinusoidal coupling, to have a $\gamma$-synchronized state. In particular, this implies the known result that that system (21) with a constant natural frequency $\mathbf {f}=\alpha \mathbf {1}$ can have a $\gamma$-synchronized state, for any $\gamma$. Moreover, if the natural frequencies are equal for all but one pair of connected nodes i and j, which differ by $\vert f_i-f_j\vert =c$, then the synchrony criterion becomes $c\le \sin (\gamma )/\omega _{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 Dörfler et al. (2013) is very different, this result is reminiscent of Theorem 1, and a further investigation of this similarity might be worthwhile.

8 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 double-well potentials and finally as the dominating behaviour for more general nonlinear systems with odd coupling functions. In a certain sense system, (1) is the ‘simplest’ dynamical system of a broad class of nonlinear 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 nonzero 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 nonlinear systems on networks.

Appendix

1.1 Appendix A: Proof of Single Consensus Link Stability

We prove Theorem 1 about the stability of mixed detailed-balance stationary states $\mathbf {x}^\star$ with a single consensus link $\mathcal {L}_{=}=\lbrace i,j\rbrace$.

Proof

By Lemma 1, we know that the system Jacobian at $\mathbf {x}^\star$ equals

$$\begin{aligned} J(\mathbf {x}^\star ) = 3r\left( (\mathbf {e}_i-\mathbf {e}_j)(\mathbf {e}_i-\mathbf {e}_j)^T - 2/3Q\right) \end{aligned}$$

with unit vectors $(\mathbf {e}_i)_k=1$ iff $k=i$ and where the first term equals $Q_{=}$ which is the (rank-one) Laplacian of a single-link graph. We introduce the abbreviations $J(\mathbf {x}^\star )=J$ and $\lambda _{\max }=\lambda _{\max }(J)$. In what follows, we show that the sign of $\lambda _{\max }$ equals the sign of $(\omega _{ij}-2/3)$ which by the linearization of system (1) then implies Theorem 1.

By the Courant Minmax principle for real symmetric matrices, we know that

$$\begin{aligned} \lambda _{\max } \ge \mathbf {x}^TJ\mathbf {x} \text { for all } \mathbf {x}\in \mathbb {R}^N\text { with }\Vert \mathbf {x}\Vert = 1. \end{aligned}$$

(22)

In particular, we thus know that for $\tilde{\mathbf {x}}=Q^\dagger (\mathbf {e}_i-\mathbf {e}_j)$ and $\mathbf {x}=\tilde{\mathbf {x}}/\Vert \tilde{\mathbf {x}}\Vert$ we have

$$\begin{aligned} \lambda _{\max } \ge 3r\mathbf {x}^T\left( (\mathbf {e}_i-\mathbf {e}_j)(\mathbf {e}_i-\mathbf {e}_j)^T - 2/3Q\right) \mathbf {x} \end{aligned}$$

(23)

Taking definition (12) of the effective resistance into account, this can be rewritten as

$$\begin{aligned} \lambda _{\max } \ge 3r\omega _{ij}/\Vert \tilde{\mathbf {x}}\Vert ^2 (\omega _{ij}-2/3). \end{aligned}$$

(24)

Since the factor $3r\omega _{ij}/\Vert \tilde{\mathbf {x}}\Vert ^2$ is positive for $r>0$ we have that $(\omega _{ij}-2/3)\ge 0 \Rightarrow \lambda _{\max } \ge 0$.

To upper-bound the largest eigenvalue, we start again from the Courant Minmax principle (22) and use the fact that equality is attained for some vector $\mathbf {y}$, i.e. such that

$$\begin{aligned} \exists \mathbf {y}\text { with }\Vert \mathbf {y}\Vert = 1\text {, such that } \lambda _{\max }=\mathbf {y}^TJ\mathbf {y} \end{aligned}$$

Writing out the full expression for the Jacobian, we find

$$\begin{aligned} \lambda _{\max } = 3r \left( \left[ (\mathbf {e}_i-\mathbf {e}_j)^T\mathbf {y}\right] ^2 - 2/3\mathbf {y}^TQ\mathbf {y}\right) \end{aligned}$$

If we denote the quadratic product by $P=\mathbf {y}^TQ\mathbf {y}$, then we know that $\mathbf {y}$ satisfies $\Vert \mathbf {y}\Vert =1, \mathbf {y}^TJ\mathbf {y} = \lambda _{\max }$ and $\mathbf {y}^TQ\mathbf {y} = P$. If we relax the first two conditions on $\mathbf {y}$, we get an upper-bound of the form

$$\begin{aligned} \lambda _{\max } \le 3r\max _{\mathbf {z}: \mathbf {z}^TQ\mathbf {z} = P}\left\{ \left[ (\mathbf {e}_i-\mathbf {e}_j)^T\mathbf {z}\right] ^2-2/3\mathbf {z}^TQ\mathbf {z} \right\} . \end{aligned}$$

(25)

This is a valid upper-bound since $\mathbf {y}$ is in the domain of $\mathbf {z}$, but the domain for $\mathbf {z}$ is larger than $\mathbf {y}$ alone and can thus potentially achieve a larger objective value. Next, we introduce the decomposition of the (positive semidefinite) Laplacian as $Q=S^TS$ and its pseudoinverse as $Q^{\dagger } = S^{\dagger T}S^{\dagger }$, where the $(N-1)\times N$ matrices $S,S^\dagger$ have entries $(S)_{ij}=(z_j)_i\sqrt{\mu _j}$ and $(S^\dagger )_{ij} = (z_j)_i\sqrt{1/\mu _j}$ and thus satisfy $S^T S^{\dagger }=I-\mathbf {u}\mathbf {u}^T/N$ and $S^{\dagger }S^T=I_{N-1}$,⁸ with all-one vector $\mathbf {u}$. With the change of variable $\tilde{\mathbf {z}} = S\mathbf {z} \Leftrightarrow \mathbf {z}=S^{\dagger T}\tilde{\mathbf {z}}$ we can then rewrite (25) as

$$\begin{aligned} \lambda _{\max } \le 3r\max _{\tilde{\mathbf {z}}:\Vert \tilde{\mathbf {z}}\Vert ^2=P}\left\{ \left[ \tilde{\mathbf {z}}^T S^{\dagger }(\mathbf {e}_i-\mathbf {e}_j)\right] ^2 -2/3\Vert \tilde{\mathbf {z}}\Vert ^2 \right\} , \end{aligned}$$

which is solved by taking $\tilde{\mathbf {z}}$ parallel to $S^{\dagger }(\mathbf {e}_i-\mathbf {e}_j)$, i.e. as $\tilde{\mathbf {z}} = S^{\dagger }(\mathbf {e}_i-\mathbf {e}_j)\sqrt{P/\omega _{ij}}$. Introducing the definition of the effective resistance in the form $\omega _{ij}=\left[ S^\dagger (\mathbf {e}_i-\mathbf {e}_j)\right] ^2$ we then find

$$\begin{aligned} \lambda _{\max } \le 3rP(\omega _{ij}-2/3). \end{aligned}$$

(26)

Since the factor 3rP is positive for $r>0$, we have that $(\omega _{ij}-2/3)\le 0\Rightarrow \lambda _{\max }\le 0$. Combining inequalities (24) and (26), we indeed find that ${\text {sign}}(\lambda _{\max })={\text {sign}}(\omega _{ij}-2/3)$ which implies Theorem 1. $\square$

1.2 Appendix B: Proof of Proposition 2

As in “Appendix A”, we will bound the eigenvalues of the Jacobian matrix $J(\mathbf {x}^\star )$ starting from Lemma 1. For consensus links $\mathcal {L}_{=}$ we have

$$\begin{aligned} J(\mathbf {x}^\star ) = 3r\sum _{(i,j)\in \mathcal {L}_{=}}(\mathbf {e}_i-\mathbf {e}_j)(\mathbf {e}_i-\mathbf {e}_j)^T - 2rQ. \end{aligned}$$

By Courants Minmax principle, we can write

$$\begin{aligned} \lambda _{\max }&\ge c(\mathbf {e}_a-\mathbf {e}_b)^TQ^\dagger JQ^\dagger (\mathbf {e}_a-\mathbf {e}_b)\nonumber \\&= 3rc\left( \omega _{ab}^2 - 2/3\omega _{ab} \right) \nonumber \\&\quad +\, \sum _{(i,j)\in \mathcal {L}\backslash \lbrace (a,b)\rbrace }c\left( (\mathbf {e}_a-\mathbf {e}_b)^TQ^\dagger (\mathbf {e}_i-\mathbf {e}_j)\right) ^2\nonumber \\&\ge 3rc\omega _{ab}(\omega _{ab}-2/3) \end{aligned}$$

(27)

for any $(a,b)\in \mathcal {L}_{=}$ and with positive normalizing constant $c=\Vert Q^\dagger (\mathbf {e}_a-\mathbf {e}_b)\Vert ^{-2}$. Since (27) holds for any consensus link, it is clear that $\max _{\mathcal {L}_{=}}\omega _{ab}>2/3\Rightarrow \lambda _{\max }>0$ proving the first stability condition.

For the upper-bound, we proceed similar to “Appendix A” and let $\mathbf {y}$ be the normalized vector such that $\mathbf {y}^TJ\mathbf {y}=\lambda _{\max }$ and define $P=\mathbf {y}^TQ\mathbf {y}$. We can then write

$$\begin{aligned} \lambda _{\max }&= \mathbf {y}^TJ\mathbf {y} = 3r\mathbf {y}^TQ_{=}\mathbf {y} - 2r\mathbf {y}^TQ\mathbf {y}\nonumber \\&\le 3r \max _{\mathbf {z}:\mathbf {z}^TQ\mathbf {z}=P}\left\{ \sum _{(i,j)\in \mathcal {L}_{=}}\left[ (\mathbf {e}_i-\mathbf {e}_j)^T\mathbf {z}\right] ^2 - 2/3P\right\} \nonumber \\&\le 3r\sum _{(i,j)\in \mathcal {L}_{=}}\left\{ \max _{\mathbf {z}:\mathbf {z}^TQ\mathbf {z}=P}\left[ (\mathbf {e}_i-\mathbf {e}_j)^T\mathbf {z}\right] ^2 - 2/3P\right\} \nonumber \\&= 3rP\left( \sum _{(i,j)\in \mathcal {L}_{=}}\omega _{ij}-2/3\right) \end{aligned}$$

(28)

where the maximum in the third line is achieved by $\mathbf {z}=S^{\dagger }(\mathbf {e}_i-\mathbf {e}_j)(\sqrt{P/\omega _{ij}})$ for each consensus link $i\sim j$. The upper-bound (28) shows that if $\sum _{\mathcal {L}_{=}}\omega _{ij}<2/3\Rightarrow \lambda _{\max }<0$ proving the second stability condition.$\square$

1.3 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 Prasse and Van Mieghem (2020), we will show that if at some time t, the state is of the form $\mathbf {x}(t)=\alpha (t) \mathbf {x}^\star$ with $\mathbf {x}^\star$ a stationary state of system (1) parallel to an eigenvector of the Laplacian matrix with eigenvalue $\mu$ (i.e. $\mathbf {x}^\star$ is an eigenstate), the dynamic equations (1) simplify to an equation for $\alpha (t)$. This proves that the subspace spanned by the vector (eigenstate) $\mathbf {x}^\star$ is a positive invariant subspace, as for any initial state $\mathbf {x}_0=\alpha _0\mathbf {x}^\star$ the solution is of the form $\mathbf {x}(t,\mathbf {x}_0)=\alpha (t)\mathbf {x}^\star$. Secondly, we show that the time-dependent coefficient $\alpha (t)$ can be solved exactly as the solution of a 1-dimensional Bernoulli differential equation.

Proof

The system equations for $\mathbf {x}(t)$ are given by

$$\begin{aligned} \frac{\mathrm{d}x_i}{\mathrm{d}t} = \sum _{j\sim i}r(x_i-x_j) - (x_i-x_j)^3 \text { for all }i. \end{aligned}$$

If at some time t the state is of the form $\mathbf {x}(t)=\alpha (t)\mathbf {x}^\star$, where $\mathbf {x}^\star$ is an eigenstate of the system satisfying conditions (15), these equations become

$$\begin{aligned} \frac{\mathrm{d}[\alpha (t)x^\star _i]}{\mathrm{d}t} = \sum _{j\sim i}r\alpha (t)(x^\star _i-x^\star _j) - \alpha ^3(t)(x^\star _i-x^\star _j)^3. \end{aligned}$$

(29)

Since $\mathbf {x}^\star$ is a stationary state of the dynamics, we have that $\sum _{j\sim i}(x_i^\star -x_j^\star )^3=r\sum _{j\sim i}(x^\star _i-x^\star _j)$, and (29) simplifies to

$$\begin{aligned} \frac{\mathrm{d}\alpha }{\mathrm{d}t}x^\star _i = r\left( \alpha (t)-\alpha ^3(t)\right) \sum _{j\sim i}(x^\star _i-x^\star _j). \end{aligned}$$

Next, as $\mathbf {x}^\star$ is an eigenvector of the Laplacian Q, i.e. with $\mu x^\star _i = \sum _{j\sim i}(x^\star _i-x^\star _j)$ we can rewrite this as

$$\begin{aligned} \frac{\mathrm{d}\alpha }{\mathrm{d}t}x^\star _i = r\mu \left( \alpha (t)-\alpha ^3(t)\right) x^\star _i. \end{aligned}$$

This shows that, for all $t'>t$ the solution will be of the form $\mathbf {x}(t',\mathbf {x}(t)) = \alpha (t')\mathbf {x}^\star$ and thus that $\left\{ c\mathbf {x}^\star \text { for some } c\in \mathbb {R}\right\}$ is a positive invariant set for the dynamics.⁹

Furthermore, the equation for $\alpha (t)$ is a 1-dimensional Bernoulli differential equation

$$\begin{aligned} \frac{\mathrm{d}\alpha }{\mathrm{d}t} = r \mu (\alpha (t)-\alpha ^3(t)), \end{aligned}$$

which can be solved by introducing $\beta =\alpha ^{-2}-1$ with $\mathrm{d}\beta /\mathrm{d}t = -2\alpha ^{-3}\mathrm{d}\alpha /\mathrm{d}t$. This yields a linear differential equation $\frac{\mathrm{d}\beta }{\mathrm{d}t} = -2r\mu \beta$ with solution $\beta (t)=\beta _0e^{-2r\mu t}$. Introducing the initial condition $\beta _0=\alpha _0^{-2}-1$ and changing the variable back to $\alpha (t)$ we find the solution

$$\begin{aligned} \alpha (t)&= \left[ 1+\left( \frac{1}{\alpha _0^2}-1\right) e^{-2\mu rt}\right] ^{-1/2}\nonumber \\&=\alpha _0\left( \alpha _0^2 - \left( \alpha _0^2-1\right) e^{-2\mu rt}\right) ^{-1/2}. \end{aligned}$$

(30)

which proves Theorem 2. $\square$

Note that this solution method only works when $(\beta (t)+1)^{-1/2}$ is well defined. When $r>0$ or $r<0$ and $\vert \alpha _0\vert <1$ this is satisfied; when $r<0$ and $\vert \alpha _0\vert >1$ on the other hand, we find that $t\ge \tfrac{1}{2\mu \vert {r}\vert }\ln \left( \frac{\alpha _0^2}{\alpha _0^2-1}\right) \Rightarrow \beta (t)\le -1$ in which case solution (30) will not apply.

1.4 Appendix D: Proof of Proposition 3

Let $\mathbf {y}^\star$ be an eigenstate of the quotient graph $G^\pi$, and $\mathbf {x}$ a $(N\times 1)$ vector with entries $x_i=y^\star _k$ when $i\in \mathcal {N}_k$, defined on the nodes of G. For any k and $i\in \mathcal {N}_k$ we can then write

$$\begin{aligned} \sum _{j\sim i}(x_i-x_j)&= \sum _{m=1}^K\sum _{\begin{array}{c} j\sim i\\ j\in \mathcal {N}_m \end{array}}(x_i-x_j)\\&= \sum _{m\ne k}\sum _{\begin{array}{c} j\sim i\\ j\in \mathcal {N}_m \end{array}}(x_i-x_j)\\&= \sum _{m\ne k}(y^\star _k-y^\star _m)\sum _{\begin{array}{c} j\sim i\\ j\in \mathcal {N}_m \end{array}}1\\&= \sum _{m\sim k}d_{km}(y^\star _k-y^\star _m) = \mu y^\star _k = \mu x_i \end{aligned}$$

where we used definition (16) of the link weight $d_{km}$. This illustrates that $\mathbf {x}$ is an eigenvector of Q with eigenvalue $\mu$. Similarly, for any k and $i\in \mathcal {N}_k$ we can write

$$\begin{aligned} \sum _{j\sim i}(x_i-x_j)^3&= \sum _{m\ne k}\sum _{\begin{array}{c} j\sim i\\ j\in \mathcal {N}_m \end{array}}(x_i-x_j)^3\\&= \sum _{m\ne k}(y_k^\star -y_m^\star )^3\sum _{\begin{array}{c} j\sim i\\ j\in \mathcal {N}_m \end{array}}1\\&= \sum _{m\sim k}d_{km}(y^\star _k-y^\star _m)^3 \\&= \sum _{m\sim k}rd_{km}(y^\star _k-y^\star _m)=\sum _{j\sim i}r(x_i-x_j) \end{aligned}$$

which shows that $\mathbf {x}$ is also a stationary state of system (1). Together, these identities show that the vector $\mathbf {x}$ constructed from the eigenstate $\mathbf {y}^\star$ of $G^\pi$ is an eigenstate of G.

$\square$

1.5 Appendix E: Proof of Attraction Basins on Trees

We prove Proposition 4, which characterizes the basins of attraction for system (1) in the region of state space with potential below the critical potential $V_c=-(L-1)r^2/4$.

Proof

The non-increasing property of the potential $\dot{V}\le 0$ implies that any state $\mathbf {x}$ with $V(\mathbf {x})<V_c$ can only evolve to stationary states $\mathbf {x}^\star$ with $V(\mathbf {x}^\star )\le V(\mathbf {x})<V_c$. Since this only leaves the stable stationary states (i.e. with $\ell (\mathbf {x}^\star )=L$ and $V(\mathbf {x}^\star )=-Lr^2/4$), we have that any $\mathbf {x}$ with $V(\mathbf {x})<V_c$ can only converge to a stable stationary state. Furthermore, by partition (18) of the full state space, every state is in the basin of attraction of exactly one stationary state, which thus means that $\mathbf {x}$ is in the basin of attraction of exactly one stable stationary state. Since this holds for every state with sub-critical potential, we get decomposition (19) of the sub-critical region of state space into basins of attraction of the stable states.

Next, we derive the specific form of the basins of attraction in the sub-critical region. We start by showing that for states $\mathbf {w}$ with $V(\mathbf {w})<V_c$ we have ${\text {sign}}(x_i(t,\mathbf {w})-x_j(t,\mathbf {w}) = {\text {sign}}(w_i-w_j)$ for all links $i\sim j$ and $t>0$. In other words, the state difference over a link cannot ‘flip’ its sign. Assume for contradiction that there is a link $i\sim j$ for which the sign is flipped. Since the state differences satisfy the differential equation (3), $(x_i-x_j)(t)$ must be continuous and thus we have that

$$\begin{aligned}&\text {if }(x_i-x_j)(t_0)>0\text { and }(x_i-x_j)(t_1)<0,\\&\text {then } \exists \tau \in [t_0,t_1]\text { s.t. }(x_i-x_j)(\tau ) = 0 \end{aligned}$$

The potential of the state $\mathbf {x}=\mathbf {x}(\tau )$ at this time can then be lower-bounded by

$$\begin{aligned} V(\mathbf {x}) =&\frac{1}{4}\sum _{\begin{array}{c} a\sim b \\ \ne i\sim j \end{array}}(x_a-x_b)^2((x_a-x_b)^2-2r) \\&+ \frac{1}{4}\underbrace{(x_i-x_j)^2((x_i-x_j)^2-2r)}_{=0} \\ \ge&\sum _{\begin{array}{c} a\sim b \\ \ne i\sim j \end{array}}\frac{-r^2}{4} = -\frac{(L-1)r^2}{4} = V_c \end{aligned}$$

where we used the fact that $z^2(z^2-2r)\ge -r^2$. This bound contradicts the fact that $V(\mathbf {x}(t,\mathbf {w}))\le V(\mathbf {w})<V_c$ which means our assumption must be false. Thus, since the stable stationary states have $(x^\star _i-x^\star _j)=\pm \sqrt{r}$ and we know that ${\text {sign}}(x_i(t,\mathbf {w})-x_j(t,\mathbf {w})) = {\text {sign}}(w_i-w_j)$ for all $t>0$ and for $t\rightarrow \infty$ in particular, we know that the link-difference signs of a state in the sub-critical region of X determines its corresponding stationary state. This proves (20) in Proposition 4. $\square$

1.6 Appendix F: Stationary States on Tree Graphs

To prove Proposition 5, we will make use of the fact that each (non-singleton) tree graph T has at least one leaf node l (a node with degree $d_l=1$), and that removing a leaf node from T results in another tree $T'=T\backslash \lbrace l\rbrace$.

Proof of Proposition 5

For $N=2$, the only tree T consists of a pair of connected nodes $i\sim j$. System (6) which determines the stationary states reduces to a single equation $r(x^\star _i-x^\star _j)=(x_i^\star -x^\star _j)^3$ which solves to a detailed-balance stationary state (7). We will prove the general case of Proposition 5 by induction on the number of nodes N, with base case $N=2$. Assuming the induction hypothesis holds for $N=K$, i.e. for all trees on K nodes the stationary states are detailed-balance stationary states, we will now show that it holds for $N=K+1$ as well.

Let T be a tree graph on $N=K+1$ nodes and l one of its leaf nodes connected to one other node n. The stationary states $\mathbf {x}^\star$ of system (1) on T are found from $\mathrm{d}\mathbf {x}^\star /\mathrm{d}t=0$ which, for the leaf node, yields the equation

$$\begin{aligned} \frac{\mathrm{d}x_l^\star }{\mathrm{d}t} = 0 \Leftrightarrow r(x^\star _l-x^\star _n) - (x^\star _l-x^\star _n)^3=0, \end{aligned}$$

(31)

which is only satisfied if $(x^\star _l-x^\star _n)\in \lbrace 0,\pm \sqrt{r}\rbrace$. In other words, for any (not necessarily detailed-balance) stationary state $\mathbf {x}^\star$ of system (1) on T, the link difference $(x^\star _l-x^\star _n)$ must be either a consensus link or a dissensus link. The stationary state values on the other nodes $i\ne l$ are determined by the equations

Since leaf node l has degree $d_l=1$, the only equation where $x^\star _l$ appears is the balance equation for $x^\star _n$. Moreover, introducing $(x_n^\star -x^\star _l)\in \lbrace 0,\pm \sqrt{r}\rbrace$ in equation () for $x^\star _n$ eliminates $x^\star _l$ from the equations for $\mathbf {x}^\star$ altogether. From this, it follows that the stationary state $\mathbf {x}^\star$ can be determined from the stationary states $\mathbf {x}^{\star '}$ of the tree graph $T'=T\backslash \lbrace l\rbrace$ as $x^\star _i=x^{\star '}_i$ for all $i\ne l$ (as they obey the same equations) and with $(x^\star _l-x^\star _n)\in \lbrace 0,\pm \sqrt{r}\rbrace$ (by the solution of equation (31)). By the induction hypothesis, $\mathbf {x}^{\star '}$ is the stationary state of system (1) on a K-node tree graph and is thus a detailed-balance state. From this follows that $\mathbf {x}^\star$ will be a detailed-balance stationary state as well. $\square$

1.7 Appendix G: Stability on the Complete Graph

We prove Proposition 7 which states that stationary states $\mathbf {x}^\star$ with $V\in (N/3,2N/3)$ are stable by explicitly calculating the spectrum of the Jacobian $J(\mathbf {x}^\star )$.

Proof

The stationary state $\mathbf {x}^\star$ partitions the set of nodes $\mathcal {N}$ of $K_N$ into two disjoint sets $\mathcal {V}$ with $\vert \mathcal {V}\vert =V$ and $\overline{\mathcal {V}}=\mathcal {N}\backslash \mathcal {V}$ such that for all $i\in \mathcal {V},j\in \overline{\mathcal {V}}$ we have $(x^\star _i-x_j^\star ) = \pm \sqrt{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 $\mathcal {V}={1,\dots ,V}$ and $\overline{\mathcal {V}}={V+1,\dots ,N}$ then the Laplacian matrices can be written as

$$\begin{aligned}&Q_{=} = \begin{pmatrix} VP_V &{} 0 \\ 0 &{} (N-V)P_{N-V} \end{pmatrix}\\&Q_{\ne } = \begin{pmatrix} (N-V)I_V &{} -\mathbf {u}_V\mathbf {u}_{N-V}^T\\ -\mathbf {u}_{N-V}\mathbf {u}_V^T &{} VI_{N-V} \end{pmatrix} \end{aligned}$$

where $I_{*}$ and $\mathbf {u}_{*}$ denote the identity matrix and all-one vector of dimensions indicated by $*$, and with $P_V = I_V-\mathbf {u}_V\mathbf {u}_V^T/V$ the projector on the space orthogonal to $\mathbf {u}_V$ and similarly for $N-V$. The Jacobian of system (1) at $\mathbf {x}^\star$ can then be calculated as $J = r(Q_{=}-2Q_{\ne })$. We will show that J has four types of eigenvectors and, correspondingly, four types of eigenvalues.

Type 1: Any vector of the form $\mathbf {z} = (\mathbf {z}_V, \mathbf {0}_{N-V})^T$ with $\mathbf {z}_V$ a V-dimensional vector that satisfies $\mathbf {z}_V^T\mathbf {u}_V=0$ and $\Vert \mathbf {z}_V\Vert = 1$, will give

$$\begin{aligned} J\mathbf {z}&= r\left( Q_{=}\mathbf {z} - 2Q_{\ne }\mathbf {z}\right) \\&= r(V - 2(N-V))\mathbf {z}\\&= 3r(V-2/3N)\mathbf {z} \end{aligned}$$

which shows that $\mathbf {z}$ is an eigenvector of J with eigenvalue $3r(V-2/3N)$. Since we can choose a basis of $V-1$ vectors orthogonal to $\mathbf {u}_V$ which all are of the form of $\mathbf {z}$, the Jacobian will have $V-1$ eigenvalues equal to $3r(V-2/3N)$.

Type 2: The same approach based on vectors of the form $\mathbf {z}=(\mathbf {0}_{V},\mathbf {z}_{N-V})^T$ with $\mathbf {z}_{N-V}^T\mathbf {u}_{N-V}=0$ gives $(N-V-1)$ eigenvalues of J equal to $-3r(V-N/3)$.

Type 3: The third type of vector is given by $\mathbf {z}=(-\mathbf {u}_V/V,\mathbf {u}_{N-V}/(N-V))^T$ and has corresponding eigenvalue 0, as $J\mathbf {z}=-2rN$.

Type 4: Finally, the fourth type of vector is simply the constant vector $\mathbf {u}=(\mathbf {u}_V,\mathbf {u}_{N-V})$ with eigenvalue 0, as $J\mathbf {u}=0$.

By construction, we furthermore have that these N vectors of the form $(\mathbf {z}_V,\mathbf {0}_{N-V}), (\mathbf {0}_V,\mathbf {z}_{N-V})$, $(-\mathbf {u}_{V}/V,\mathbf {u}_{N-V}/(N-V))$,$(\mathbf {u}_{V},\mathbf {u}_{N-V})$ are orthogonal and and thus determine eigenvectors of J. The eigenvalues of the Jacobian $J(\mathbf {x}^\star )$ are thus equal to

$$\begin{aligned}&\lbrace 0\rbrace ^{(1)}, \lbrace -2rN\rbrace ^{(1)},\left\{ -3r(V-N/3)\right\} ^{(N-V-1)}, \\&\lbrace 3r(V-2N/3)\rbrace ^{(V-1)} \end{aligned}$$

where superscripts denote the multiplicity of the eigenvalues. The first zero eigenvalue corresponds to the constant eigenvector $\mathbf {u}$ which is orthogonal to the state space X and thus does not influence the stability of state $\mathbf {x}^\star$. The second zero eigenvalue is always negative for $r>0$. Finally, if $V\in (N/3,2N/3)$ then all other eigenvalues are negative, and if $V\notin [N/3,2N/3]$ they are positive, which proves the (in)stability of the stationary states $\mathbf {x}^\star$ with respect to V. When $V\in \lbrace N/3,2N/3\rbrace$, the other eigenvalues become zero and the linearization method does not provide the necessary information to determine the stability of the corresponding states. $\square$