# Partial synchronization phenomena in networks

of identical oscillators with non-linear coupling

###### Abstract

We study a Kuramoto-like model of coupled identical phase oscillators on a network, where attractive and repulsive couplings are balanced dynamically due to nonlinearity in interaction. Under a week force, an oscillator tends to follow the phase of its neighbors, but if an oscillator is compelled to follow its peers by a sufficient large number of cohesive neighbors, then it actually starts to act in the opposite manner, i.e. in anti-phase with the majority. Analytic results yield that if the repulsion parameter is small enough in comparison with the degree of the maximum hub, then the full synchronization state is locally stable. Numerical experiments are performed to explore the model beyond this threshold, where the overall cohesion is lost. We report in detail partially synchronous dynamical regimes, like stationary phase-locking, multistability, periodic and chaotic states. Via statistical analysis of different network organizations like tree, scale-free, and random ones, and different network sizes, we found a measure allowing one to predict relative abundance of partially synchronous stationary states in comparison to time-dependent ones.

###### pacs:

Lead paragraph: Large population of coupled oscillators synchronizes if the coupling is attractive and desynchronizes if it repulsive. However, the nature of coupling may depend on the strength of the compelling field: if the force on the oscillator becomes strong, it can switch its behavior from a “conformist” to a “contrarian”. We study such a population on a network. Here the oscillators connected to many others become contrarians first, so that synchrony breaks. We show that the regimes of appearing partial synchrony can be rather complex, with a large degree of multistability, and suggest a network measure which allows to predict relative abundance of static and dynamic regimes.

## I Introduction

In a seminal work Kuramoto (1975), aiming to understand synchronization phenomena, Kuramoto proposed a mathematical model of non-identical, nonlinear phase-oscillators, mutually coupled via a common mean field. Studying this system, he identified a synchronization transition to an oscillating global mode when the coupling strength is larger than a critical value, which is proportional to the range of the distribution of the natural frequencies. Over the time, subsequent outcomes based on Kuramoto propositions have shown that his approach can be used as a framework to several natural and technological systems where an ordered behavior (synchronization) emerges from the interactions of a large number of dynamical agents Acebrón et al. (2005); Strogatz (2000). Furthermore, works have shown that the Kuramoto model can be exploited as a building block to develop highly efficient strategies to process information Follmann et al. (2014); Vassilieva et al. (2011).

Recently, generalizations of the Kuramoto model toward interconnections between the elements more complex than the mean field one, have received considerable attention. Indeed, in many real-world problems, each dynamical agent interacts with a subset of the whole ensemble Dörfler et al. (2013); Leonard (2014); Sadilek and Thurner (2014), which can better described using a network topology. A variety of studies have analyzed the onset of the synchronization regime in this context. For a general class of linearly coupled identical oscillators, the Master Stability Function, originally proposed by Pecora and Carrol Pecora and Carroll (1998), allows one to determine an interval of coupling strength values that yields complete synchronization, as a function of the eigenvalues of Laplacian matrix of the coupling graph. For networks of oscillators with non-identical natural frequencies, Jadbabaie et al. Jadbabaie et al. (2004) were able to give similar bounds for the coupling strength of the Kuramoto model without the assumption of infinitely many phase-oscillators. Among related works, Ref. Franci et al. (2010) deals with a model whose natural frequency oscillators change with time, even when they are isolated. Ref. Papachristodoulou et al. (2010) explores the effects of delay in the communication between oscillators. Besides, Ref. Pecora et al. (2013) builds a bridge between graph symmetry and cluster synchronization.

Taking into consideration all of these previous results, one can roughly state that the Kuramoto transition to synchronization always happens if the coupling between oscillators is attractive; while, when it changes to be repulsive, this synchronization state is absent Tomov and Zaks (2014); Tsimring et al. (2005). However, the structure of the coupling can nontrivially depend on the level of synchrony itself. Such a dependence, called nonlinear coupling scheme, has been explored in recent theoretical Filatrella et al. (2007); Rosenblum and Pikovsky (2007); Pikovsky and Rosenblum (2009); Baibolatov et al. (2009) and experimental Temirbayev et al. (2012, 2013) studies dealing with setup of global coupling. The main effect here is the partial synchrony, which establishes at moderate coupling strengths, where the coupling is balanced between the attractive and repulsive one.

Here, we consider the effects of the non-linear coupling on a network: a set of identical oscillators, which communicate via a connected simple coupling graph. each element is forced by a (local) mean field, which encompasses the oscillators that are connected to it. The coupling function is tailored so that its influence is attractive, if the local acting field is small, or repulsive, otherwise. This coupling strategy implies that only nodes with a large enough number of connections may become repulsive. Thus, the hubs play a key role for the ensemble dynamics. A non-linear coupling parameter in the system tunes the critical quantity of connections and how cohesive this mean field must be in order to allow this transition. So, our approach can be considered as a dynamical generalization of the inhomogeneous populations of oscillators consisting of conformists and contrarians Hong and Strogatz (2011). However, the type of oscillator’s depends on the force acting on it.

Overall dynamics in the model can be qualitatively understood as follows: Let us assume initially that all the mean fields are small. Then, there are only attractive interactions (conformists) in the system. So, in a first moment, they start to mutually adjust their phases. Above a threshold, the most connected oscillators start to feel a repulsive effect that drives them away from the synchronous state. In other words, if an oscillator has a sufficiently large number of neighbors and if it suffers enough cohesive pressure from them, instead of attractiveness, it becomes a contrarian, wishing to be in anti-phase with the force. Then, due to the repulsiveness of some nodes, other mean fields may also become smaller. Finally, this tendency can shift nodes to attractiveness again. As a consequence, an intermediate configuration may emerge due to the balance these conflicting tendencies in the system.

Depending on the non-linear coupling parameter, we report a variety of qualitative dynamic behaviors. In general, for small values of the non-linear coupling parameter, we observed full synchronization and phase-locked states. When this parameter is increased, multistability, periodic and chaotic dynamics take place.

The paper is organized as follows. Initially, we discuss the basic details of the model in section II. In section III, the analytical result about the stability of full synchronization is presented. Numerical experiments in section IV illustrate different possible regimes that the present model can display. Finally, in section IV.3, we perform a numerical exploration to address the correlation of stationary phase locking states with partial synchronization with the network parameters, by exploring different network topologies and sizes.

## Ii Model of oscillator network with nonlinear coupling

Let us consider a system of identical phase-oscillators represented by coupled through a simple and connected undirected graph . The dynamics for the -th oscillator is given by the following ordinary differential equation

(1) |

where denotes the set of neighbors of in the coupling graph . Equations (1) are formulated in the reference frame rotating with the common frequency of the oscillators, so that the latter one does not appear in the equations. The time is normalized by the linear coupling strength. The main feature of our model (1) is the nonlinear coupling parameter, , which modifies the coupling at each node, with a standard setup of the Kuramoto model on a network is recovered. This parameter enters as a coefficient at the square of the norm of the local order parameter:

(2) |

which measures the magnitude of the force acting on oscillator with index . On the other hand, we represent the (global) order parameter by

(3) |

where is its norm and is its phase.

We stress that is not normalized (in the sense that there is no division by the number of the terms in the summation, like in ), as it measures the total action of the neighbors on the -th oscillator, which is called local mean field. Simple calculations show that , where denotes the degree of the -th vertex, that is, the number of incoming or outgoing connection, since the graph is undirected. Thus, a necessary condition for a node to suffer repulsive coupling, i.e. , is that .

The introduced order parameters and are maximal in the case of full synchronization , whiles they decrease when oscillators begin to move apart from each other.

If , where , then all oscillator will attract each other so that the full synchronization is established. Next, if becomes larger than , the corresponding oscillator begins to be repulsive related to its local mean field, and the full synchronization breaks. As a result, may decrease and switch again the node to be attractive. Depending on the coupling graph , on the initial conditions , and on the intensity of the nonlinear coupling parameter , numerical simulation reveals that model (1) can exhibit different qualitative behaviors.

If the largest degree in the coupling graph satisfies , with , then no node can be repulsive. We demonstrate in the next section via the Lyapunov analyses, that in fact this condition guarantees that the full synchronized state is stable.

## Iii Stability of full synchronization

The basic steps to obtain the results in this section follow the main ideas from Jadbabaie et al. (2004). We begin presenting some preliminary concepts, including elements of the graph theory needed, and a generalized norm of the order parameter to define our Lyapunov function.

Let be the directed incidence matrix of a graph . Thus, is a matrix with rows and columns, where is the number of directed edges of the matrix. The number of undirected edges, \ieignoring the direction, equals is . The columns of represent the edges of the graph: if the -th arrow (directed edge) of the graph goes from to , then the -th column of is zero, except at positions and , where and . Regarding the dynamics of the system, an arrow from node to node in the graph means that node influences node . Although the directed incidence matrix is generally defined for directed graphs, it must be emphasized that only undirected graphs are considered here. We abuse terminology and identify a graph with its adjacency matrix, which is an matrix where ; , if there is and edge between nodes ; and , otherwise. So, . Another common characterization of a graph is its Laplacian matrix, . One can check that . A simple illustration of these concepts is given at Fig. 1.

The usage of the directed incidence matrix allows us to rewrite model (1) in a vector form as follows:

(4) |

where , and stands for the matrix with the elements of a vector on the leading diagonal, and elsewhere.

The square of the global order parameter can be expressed as

However, to build our Lyapunov function, we define a generalized norm of order as

(5) |

Note that requires the sum of all with (for ), but its generalization takes the sum () only through the edges of the graph. In the case of full coupling graph, direct substitution yields that both global and generalized norm of the order parameter have the same expression .

For any connected symmetrical coupling graph, one can check that the maximum of
is the unit, and that if and only if this value is
achieved ^{1}^{1}1On the other hand, the minimum value of
does not necessarily correspond to ..

Let

(6) |

be a candidate Lyapunov function. It is clear that the minimum value of corresponds to the maximum value of , which is equivalent to the fully synchronized state.

In fact, algebraic manipulations reveal that

(7) |

and that the differential of is given by

(8) |

As a result, we synthesize in the next theorem the previous intuitively stated idea that if is small enough then full synchronization is a robust phenomenon related to small perturbations over initial conditions.

###### Theorem 1.

In Model (1), if is smaller than a critical value , then the synchronized stated () is Lyapunov stable.

###### Proof.

Consider the potential field defined in (6). So, using the vector form of the model (4) and the expression of the differential from (8), we have that equals to

If we set , then we have that is larger or equal than . Moreover, we can also define a lower bound for , since ; where is the algebraic connectivity of the graph. In the last inequality we used that and that both matrices and have the same non-trivial eigenvalues , where is strictly larger than zero because the coupling graph is connected Godsil and Royle (2001). Therefore,

As a result, implies that , then the fully synchronized state is stable. ∎

## Iv Dynamics of partially synchronous states

In this section numerical simulations are performed to illustrate the rich repertoire of behaviors that model (1) may exhibit, specially beyond the threshold , where Theorem 1 cannot be applied.

### iv.1 Quantification of dynamical regimes

The numerical integration scheme applied is a forth order Adams-Bashforth-Moulton Method (see Burden and Faires (2010)) with discretization time step . We calculate the partial synchronization metric s from Ref. Gómez-Gardeñes et al. (2007), which, for every two oscillators in the network, takes values , indicating how much the mean phase difference between and varies in the time interval , with . This metric is defined as

One can check that if for some constant , then the exponent in the previous integral is constant and . Nevertheless, if assumes every possible value over the unit circumference with not clear trend, then is close to zero. Now, we average contributions of all neighbor oscillators under a graph with nodes to write

where is the quantity of undirected edges in the graph.

To exclude transients and to detect the statistically stationary state, we adopted the following procedure. For all experiments the time interval is always considered as transient time. Then, the numerical integration is performed in the subsequent intervals , with , until the first such that , or is achieved. Only such a time interval is regarded as non-transient. For the subsequent analysis, we use values of the phases in the stationary time interval regime (whose beginning is shifted to without loss of generality) at points .

### iv.2 Examples of complex behaviors

As it was claimed before, in dependence on the network structure, very different types of the dynamics are possible. In order to give impression on it, we present simulations of model (1) with two different coupling graphs displayed as inserts in Fig. 2. Both networks have nodes and they differ only by the rewiring of a single edge. We performed simulations for random initial conditions chosen with uniform distribution over for each experiment. For all these initial conditions , the norm of the order parameter , according to (3), is computed from the time series. As explained in the previous section, in these calculations a transient time is eliminated and that a statistically stationary regime of units of time and points is considered. Then, also for each distinct initial condition, the maximum, average and minimum values of the associated norm of the order parameter are computed, respectively denoted by ; ; and . Of course, converges to a constant if and only if . Now, having different simulations for a fixed coupling graph, we evaluated the maximum, average and minimum value of the average value of the norm of the order parameters over this ensemble, respectively denoted by ; ; and . So, if the norm of the order parameter converges to the same value for all initial conditions simulated, then . For the cases where the norm of the order parameter does not converge over all initial conditions, it will be useful to examine the overall maximum and overall minimum values of the norm of the order parameter, respectively denoted by ; and . Thus, if there is no fixed phase synchronization for all the initial conditions simulated, but the norm of the order parameter presents only small deviations around a common value, then the gap between and is also small. Also notice that , since for all initial conditions. Finally, the maximum Lyapunov exponent for each initial condition is also computed, according to the algorithm in Alligood et al. (1997). The maximum Lyapunov exponent over all the chosen initial conditions is also analyzed.

We now describe different regimes observed in the networks, using also Fig. 3, where we depict time series of for some particular choices of , indicated as green letters in the upper panel from Fig. 2 (this is the case we choose for illustrating different regimes). Notice that in both cases, so Theorem 1 guarantees that for the full synchronization state, , is locally stable as illustrated in Fig. 3 (a) (with ).

Panel (a) in Fig. 3 illustrates full synchronization in the network for . For slightly bigger than , simulations suggest that a stationary regime of partial phase synchronization, where , is locally stable as shown in Fig. 3(b) (). Details of this state are clear from Fig. 4. There we show the that the synchronization between the individual oscillators is complete if measured by quantity , and all the oscillators have the same frequency. However, the oscillators are split into two groups with a constant phase shift between them; this division originates in the edge which connects the two largest hubs in the network (vertexes ).

For larger values of , the regimes are still static but with multistability. For instance, at (see Fig. 3 (c)) two stable configurations emerge with , with (black) or (blue), depending on the initial condition. Fig. 5, which is analogous to Fig. 4, shows the existence of three subgroups, whose members may vary according to the initial condition.

Other types of multistabilities appear for instance at and , as illustrate in Figs. 3 (d,e). For (panel d) some initial conditions do no converge to a fixed phase synchronization, but to a regime where the order parameter is periodic in time. For (panel c), the norm of the order parameter of all trajectories simulated becomes periodic. Fig. 6 provides an illustration of this regime.

Finally, for (Fig. 3 (f)), one observes a chaotic state with , the distribution of phases and frequencies is illustrated in Fig. 7.

If , we also obtained multistability, with the coexistence of solutions converging to phase-lock and irregular order parameter after the transient, similar to Fig. 3(d).

Now, we compare the results for two slightly different networks depicted in panels (a) and (b) in Fig. 2. The interval of values of with fixed phase synchronization for all initial conditions simulated is very similar for both networks, namely ; also multistability of static partial synchronous regimes have been observed in both cases.

When , contrary to case (a), we observed that the solution for all initial conditions converged to the same phase-lock regime, similar to Fig. 3 (b).

In conclusion of this section, in Fig. 8 we show simulation results for two other networks. Panel (a) shows a random network with nodes and undirected edges. Here predominantly static regimes are observed, only in small ranges of coupling constant chaos with a positive Lyapunov exponent appears. Static regimes, however, demonstrate a large degree of multistability. In panel (b) we show a scale-free network with nodes and undirected edges. Here static states are rare, typically irregular regimes with low values of the order parameter are observed.

### iv.3 Dependence of partial synchronization regimes on network structure

We have seen that partially synchronous states can be rather different even for similar networks. It is therefore difficult to make general predictions for a relation between the network properties and the dynamical behaviors. Here we attempt such a description, focusing on the property of abundance of static regimes in comparison to time-dependent ones. For this purpose we define the convergence index as the ratio of values of such that converges to a constant value, considering all the random initial conditions. So, both networks Fig. 2 have similar values of the index: in case (a) while in case (b). In contradistinction, network shown in Fig. 8(a) has very large value of the index , while that in Fig. 8(b) a rather low value .

In order to explore which features of the coupling graph are related with , we performed numerical experiments with three sets of graphs, with nodes. Each set consists in three common types of networks, each one with members, generated as: (i) random (Erdös-Rényi) graphs with edges; (ii) scale-free graphs, also with ; and (iii) tree graphs ( edges). The Barabási-Albert algorithm is applied for the last two types of networks (ii),(iii), with an initial clique of nodes and with other nodes been connected to existing ones. For the -edges scale-free graphs, we fixed and ; while for the tree graphs ( edges scale-free graphs), . We point out that all graphs created are connected and symmetrical. Additionally, three sets of initial conditions , with uniform distribution over and , have been explored. So, for each of the coupling graphs we computed its correspondent values by numerical integration of model (1) for .

In Table 1 we report the mean value and the standard deviation of for each topology and size of coupling graph. From these data we see that the mean value of increases if we go from tree to scale-free and to random graphs, respectively. However, these difference becomes less noticeable for larger values of . Both the mean value and the standard deviation of decrease with larger networks.

Network | |||
---|---|---|---|

Tree | 0.421 (0.260) | 0.016 (0.006) | 0.008 (0.004) |

Scale-free | 0.857 (0.029) | 0.050 (0.015) | 0.013 (0.003) |

Random | 0.872 (0.090) | 0.183 (0.063) | 0.077 (0.022) |

We have explored different networks metrics, searching for one mostly correlated with the convergence index . Let denote the Laplacian eigenvalues of the coupling graph Godsil et al. (2001). Recall that this graph is assumed to be simple and connected. We stress that these eigenvalues express fundamental characteristics of the graph. For instance, is related with graph diameter and with its largest degree size.

We found that the quantity , defined as the ratio between the maximum eigenvalue and the average of the non-trivial eigenvalues of the Laplacian matrix of the graph, is rather suitable for this purpose. Formally, it is defined as

In Fig. 9 a correlation plot between and measure for the correspondent graph is presented. From there, we observe a clear trend indicating that the greater the value of is, the smaller is the value of . Independently of the network type and size, static regimes of partial synchronization, full synchronization and phase-lock, are typical for values of , like in the experiments from Fig. 2. On the other hand, graphs with larger values of this measure yields more irregular dynamics, like time-dependent periodic and chaotic regimes, as the ones from Fig. 8.

## V Conclusion

In this work we introduced and studied a Kuramoto-like model of identical oscillators with non-linear coupling. Our main parameter was , which governs the coupling nonlinearity strength. It is clear that the most influence of nonlinearity in the coupling is on the hubs which experience strong forcing from many connected oscillators, while less connected nodes may still operate in a linear-coupling regime.

We proved that if this parameter is smaller than the inverse of the square of the maximum vertex degree in the network, then the full synchronized state is stable. Via numerical experiments, we showed that our model can display a variety of other qualitative behaviors of partial synchronization, like stationary phase locking, multistability, periodic order parameter variations, and chaotic regimes. We explored the relative abundance of stationary phase locking regimes under different network topologies. Statistical analysis performed suggests that tree graphs are much less likely to exhibit stationary phase locking in comparison with scale-free or random networks. In addition, this type of behavior becomes more rare if we increase network sizes, irrespective to the network topology. Finally, we also found a good correlation between the ration between the maximum eigenvalue and the average of the non-trivial eigenvalues of the Laplacian matrix of the graph, and the proportion of the repulsion parameter values which yield stationary phase locking. The greater this measure is, the smaller tend to be presence of stationary phase locking states in the system.

As a future research, we want to investigate analytical conditions and correlations involving other graph measures related to other forms of synchronization in the model.

###### Acknowledgements.

We would like to thank the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES (Process: BEX 10571/13-2) for financial support. AP V.V. thanks the IRTG 1740/TRP 2011/50151-0, funded by the DFG/FAPESP, CNPq, and the grant/agreement 02.B.49.21.0003 of August 27, 2013 between the Russian Ministry of Education and Science and Lobachevsky State University of Nizhni Novgorod.## References

- Kuramoto (1975) Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, Vol. 39, edited by H. Araki (Springer Berlin Heidelberg, 1975) pp. 420–422.
- Acebrón et al. (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
- Strogatz (2000) S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
- Follmann et al. (2014) R. Follmann, E. Macau, E. Rosa, and J. Piqueira, Neural Networks and Learning Systems, IEEE Transactions on PP, 1 (2014).
- Vassilieva et al. (2011) E. Vassilieva, G. Pinto, J. Acacio de Barros, and P. Suppes, Neural Networks, IEEE Transactions on 22, 84 (2011).
- Dörfler et al. (2013) F. Dörfler, M. Chertkov, and F. Bullo, Proceedings of the National Academy of Sciences 110, 2005 (2013).
- Leonard (2014) N. E. Leonard, Annual Reviews in Control 38, 171 (2014).
- Sadilek and Thurner (2014) M. Sadilek and S. Thurner, arXiv preprint arXiv:1409.5352 (2014).
- Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Physical Review Letters 80, 2109 (1998).
- Jadbabaie et al. (2004) A. Jadbabaie, N. Motee, and M. Barahona, in American Control Conference, 2004. Proceedings of the 2004, Vol. 5 (2004) pp. 4296–4301 vol.5.
- Franci et al. (2010) A. Franci, A. Chaillet, and W. Pasillas-Lépine, in Decision and Control (CDC), 2010 49th IEEE Conference on (IEEE, 2010) pp. 1587–1592.
- Papachristodoulou et al. (2010) A. Papachristodoulou, A. Jadbabaie, and U. Munz, Automatic Control, IEEE Transactions on 55, 1471 (2010).
- Pecora et al. (2013) L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, arXiv preprint arXiv:1309.6605 (2013).
- Tomov and Zaks (2014) P. Tomov and M. Zaks, in Nonlinear Dynamics of Electronic Systems (Springer, 2014) pp. 246–253.
- Tsimring et al. (2005) L. S. Tsimring, N. F. Rulkov, M. L. Larsen, and M. Gabbay, Phys. Rev. Lett. 95, 014101 (2005).
- Filatrella et al. (2007) G. Filatrella, N. F. Pedersen, and K. Wiesenfeld, Phys. Rev. E 75, 017201 (2007).
- Rosenblum and Pikovsky (2007) M. Rosenblum and A. Pikovsky, Phys. Rev. Lett. 98, 064101 (2007).
- Pikovsky and Rosenblum (2009) A. Pikovsky and M. Rosenblum, Physica D: Nonlinear Phenomena 238, 27 (2009).
- Baibolatov et al. (2009) Y. Baibolatov, M. Rosenblum, Z. Z. Zhanabaev, M. Kyzgarina, and A. Pikovsky, Phys. Rev. E 80, 046211 (2009).
- Temirbayev et al. (2012) A. A. Temirbayev, Z. Z. Zhanabaev, S. B. Tarasov, V. I. Ponomarenko, and M. Rosenblum, Phys. Rev. E 85, 015204 (2012).
- Temirbayev et al. (2013) A. A. Temirbayev, Y. D. Nalibayev, Z. Z. Zhanabaev, V. I. Ponomarenko, and M. Rosenblum, Phys. Rev. E 87, 062917 (2013).
- Hong and Strogatz (2011) H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106, 054102 (2011).
- (23) On the other hand, the minimum value of does not necessarily correspond to .
- Godsil and Royle (2001) C. Godsil and G. Royle, Algebraic Graph Theory, Graduate Texts in Mathematics., Vol. 207 (volume 207 of Graduate Texts in Mathematics. Springer, 2001).
- Burden and Faires (2010) R. L. Burden and J. D. Faires, Numerical Analysis (Brooks Cole, 2010).
- Gómez-Gardeñes et al. (2007) J. Gómez-Gardeñes, Y. Moreno, and A. Arenas, Phys. Rev. Lett. 98, 034101 (2007).
- Alligood et al. (1997) K. Alligood, T. Sauer, and J. Yorke, Chaos: An Introduction to Dynamical Systems, Chaos: An Introduction to Dynamical Systems (New York, NY, 1997).
- Godsil et al. (2001) C. D. Godsil, G. Royle, and C. Godsil, Algebraic graph theory, Vol. 207 (Springer New York, 2001).