Algebraic Geometrization of the Kuramoto Model: Equilibria and Stability Analysis
Dhagash Mehta
Noah S. Daleo
Florian D\"orfler
Jonathan D.~Hauenstein
\noindent Finding equilibria of the finite size Kuramoto model amounts to solving a nonlinear system of equations, which is an important yet challenging problem. We translate this into an algebraic geometry problem and use numerical methods to find all of the equilibria for various choices of coupling constants $K$, natural frequencies, and on different graphs. We note that for even modest sizes ($N\sim 10-20$), the number of equilibria is already more than 100,000. We analyze the stability of each computed equilibrium as well as the configuration of angles. Our exploration of the equilibrium landscape leads to unexpected and possibly surprising results including non-monotonicity in the number of equilibria, a predictable pattern in the indices of equilibria, counter-examples to popular conjectures, multi-stable equilibrium landscapes, scenarios with only unstable equilibria, and multiple distinct extrema in the stable equilibrium~distribution as a function of the number of cycles~in~the~graph.
kuramoto_draft.tex
\documentclass[10pt,twocolumn,prl,superscriptaddress,english]{revtex4-1}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}

\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
%\documentclass[english]{article}

\bibliographystyle{plain}

% latexsym
\usepackage{babel}

%\usepackage{multicol} % open this to use the two-column format

%\usepackage{fullpage}

\widowpenalty=300
\clubpenalty=300

%opening

\makeatother

\usepackage{babel}


% commenting options
\usepackage{color}
\newcommand{\tb}[0]{\color{blue}}
\newcommand{\tr}[0]{\color{red}}
\newcommand\fdmargin[1]{\marginpar{\color{blue}\tiny\ttfamily From FD: #1}}
\newcommand\redmargin[1]{\marginpar{\color{red}\tiny\ttfamily #1}}
\newcommand{\todo}[1]{\par\noindent{\raggedright\color{red}\texttt{#1}\par\marginpar{$\star$}}}
\newcommand{\bC}{{\mathbb C}}
\newtheorem{conj}{Conjecture}

\begin{document}

\title{Algebraic Geometrization of the Kuramoto Model: Equilibria and Stability Analysis}

\author{Dhagash Mehta}
\email{dmehta@nd.edu}
\affiliation{Dept.~of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA.}

\author{Noah S. Daleo}
\email{nsdaleo@ncsu.edu}
\affiliation{Dept.~of Mathematics, North Carolina State University, Raleigh, NC 27695, USA.}

\author{Florian D\"orfler}
\email{{dorfler@ethz.ch}}
\affiliation{Automatic Control Laboratory, Swiss Federal Institute of Technology (ETH) Z\"urich, 8092 Z\"urich, Switzerland.}

\author{Jonathan D.~Hauenstein}
\email{hauenstein@nd.edu}
\affiliation{Dept.~of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA.}

\begin{abstract}
\noindent 
Finding equilibria of the finite size Kuramoto model amounts to 
solving a nonlinear system of equations, which is an important yet challenging problem. 
We translate this into an algebraic geometry problem and 
use numerical methods to find all of the equilibria for various 
choices of coupling constants $K$, natural frequencies, and on different graphs.
We note that for even modest sizes ($N\sim 10-20$), 
the number of equilibria is already more than 100,000.
We analyze the stability of each computed equilibrium as well as the configuration of angles. 
Our exploration of the equilibrium landscape leads to unexpected and possibly surprising results
including non-monotonicity in the number of equilibria, a predictable pattern in the 
indices of equilibria, counter-examples to popular conjectures, multi-stable equilibrium landscapes, scenarios with only unstable equilibria, and multiple distinct extrema in the stable equilibrium~distribution as a function of the number of cycles~in~the~graph.
\end{abstract}

\maketitle

\noindent\textit{Introduction.~-~}%\label{Sec:Introduction}
%\section*{Introduction}
The Kuramoto model is a fascinating model proposed in 1975
to study synchronization phenomena \cite{kuramoto1975self}
that has gained attention from various scientific communities, including
biology, chemistry, physics, and electrical engineering,
due to its applicability. 
This model has been used to study various phenomena including
neural networks, chemical oscillators, Josephson junctions and laser arrays, 
power grids, particle coordination, 
spin glass models, and rhythmic applause~\cite{acebron2005kuramoto,strogatz2000kuramoto,FD-FB:13b}.
%, and biological synchronization phenomena
%\fdmargin{replaced an old reference by \cite{FD-FB:13b}}
  
The Kuramoto model is defined as a system of autonomous ordinary differential equations as:
\begin{equation}\label{Eq:KuramotoStandard}
\frac{d\theta_{i}}{d t}=\omega_{i} - \frac{K}{N}\sum_{j=1}^{N}a_{i,j}\sin(\theta_{i}-\theta_{j}),\mbox{~for~}i=1,...,N,
\end{equation}
where $K$ is the coupling strength, $N$ is the number of oscillators, 
$\Omega=(\omega_{1},\dots,\omega_{N})$ is the vector of intrinsic natural frequencies, and 
%\fdmargin{You had $a_{i,j} \in \{-1,0,1\}$ before. But our analysis in this article is restricted to nonnegative weights. i changed it. But it would be an interesting next step to look at some negative weights}
$a_{i,j} \in \{0,1\}$ is the $(i,j)$th element of the 
adjacency matrix of the coupling graph.
The natural frequencies $\omega_i$ indicate how the system oscillates in 
the absence of any dissipation or exogenous forces.
%; these frequencies may be chosen from some distribution $g(\omega)$.

The equilibrium conditions are $\frac{d\theta_{i}}{d t}=0$ for all $i$.
This system of equations has an $O(2)$ freedom, i.e.,
for any \mbox{$\alpha\in(-\pi,\pi]$}, the equations are invariant 
under replacing all $\theta_i$ with $\theta_i+\alpha$.
This rotational symmetry leads to a continua of 
equilibria. To remove this $O(2)$ freedom resulting in 
finitely many equilibria, we fix one of the angles, say, 
$\theta_{N}=0$, and remove the equation 
$\frac{d \theta_{N}}{d t}=0$ from the system. The remaining 
system consists of $N-1$ nonlinear equations in $N-1$ angles.

Provided the coupling strength $K$ is strong enough, 
the oscillators will synchronize as $t\rightarrow\infty$. 
In this setup, a critical coupling $K_c(N)$ exists at which 
the number of stable equilibria switches from $0$ to a nonzero value. 
In the special case of $N\rightarrow\infty$ and long-range (all-to-all) coupling $a_{ij} = 1$, one may analytically 
compute $K_c(N)$.  
However, for the finite size Kuramoto model, such an analysis may turn out to be very difficult. 
In particular, finding all equilibria, analyzing stability, and 
finding $K_{c}(N)$ is known to be
prohibitively difficult for a finite but large oscillator population. % due to the nonlinearity of the equations.

We point out that the equilibria of system (\ref{Eq:KuramotoStandard}) can also be viewed as the stationary points of the potential energy landscape drawn by the 
mean-field XY model with an exogeneous perturbation term:
%\fdmargin{added a constant to the potential to make it locally quadratic in $(\theta_{i} - \theta_{j})$}
\begin{equation}
\label{eq: potential}
%{\tb
 V(\theta) = \frac{K}{2N} \sum_{i,j=1}^{N} a_{i,j}(1-\cos (\theta_{i} - \theta_{j})) - \sum_{i=1}^{N}\omega_{i} \theta_{i},
% }
\end{equation}
whose gradient %first derivatives with respect to $\theta$-angles 
reproduces the right-hand side of equation~\eqref{Eq:KuramotoStandard}. Hence, in the following,
we use the words equilibria and stationary points interchangeably.

All stationary points of the finite $N$ mean-field XY model (i.e., the Kuramoto model with homogeneous frequencies)
were identified in Ref.~\cite{Casetti:June2003:0022-4715:1091}. 
Building on Ref.~\cite{Casetti:June2003:0022-4715:1091},
all stationary points of the one-dimensional 
nearest-neighbour $XY$ model (i.e.,~the Kuramoto model with local coupling) 
for any given $N$ with either periodic or anti-periodic boundary conditions
have been found~\cite{Mehta:2009,Mehta:2010pe,vonSmekal:2007ns,vonSmekal:2008es}.

Using these solutions, a class of stationary points of the 2-dimensional 
nearest-neighbour $XY$ model \cite{Nerattini:2012pi} and the $XY$ model
with long-range interactions \cite{kastner2011stationary} were built and analysed (see also \cite{Hughes:2012hg,Mehta:2014jla}). 
In Ref.~\cite{Mehta:2009,Mehta:2010pe,Mehta:2009zv,Hughes:2012hg}, all of the stationary points for small lattices were found 
%\fdmargin{added two related references}
using algebraic geometry methods. Bounds on the number of equilibria  \cite{JB-CIB:82} as well as some counterintuitive 
examples to plausible conjectures \cite{AA-SS-VP:81} have been reported for the same model in the domain of power systems. 

In Ref.~\cite{Casetti:June2003:0022-4715:1091} for the finite $N$ mean-field $XY$ model, and
in Refs.~\cite{Mehta:2009,Nerattini:2012pi} for the nearest neighbour $XY$ models, it was 
shown that there were exponentially many isolated stationary solutions 
as $N$ increases.  Moreover, even after breaking the global $O(2)$ 
symmetry, there were continuous solutions at the maximum value of the 
energy which were independently observed and 
termed as {\em incoherent manifolds} in Ref.~\cite{strogatz1991stability}.

In Refs.~\cite{verwoerd2009computing,verwoerd2008global,dorfler2011critical}, 
necessary and sufficient conditions are given for fixed points to 
exist for the finite 
size Kuramoto model for complete and bipartite graphs, and explicit upper 
and lower bounds of $K_c(N)$ for these systems were also computed, followed 
by providing an algorithm to compute $K_c(N)$. For the complete graph, similar results were presented in \cite{DA-JAR:04,REM-SHS:05}, where it was additionally shown that there is exactly one single stable equilibrium for $K > K_c(N)$. 
In Ref. \cite{FD-MC-FB:11v-pnas}, an analogous result was shown for acyclic graphs, short cycles, and complete graphs as well as combinations thereof.
In the case of 
homogeneous natural frequencies networks with sufficiently high 
nodal degrees, the only stable fixed point is known to be the phase-synchronized solution \cite{taylor2012there}.
In Ref.~\cite{ochab2009synchronization}, 
all of the stable synchronized 
states were classified for the one-dimensional Kuramoto model on a ring 
graph with random natural frequencies in addition to computing a lower 
bound on $K_c(N)$ (see also Ref.~\cite{ermentrout1985behavior,strogatz1988phase, tilles2011multistable}). 

\noindent\textit{Algebraic Geometry Interpretation and Setup.~-~}%\label{Sec:Alg_geom}
%\section*{Algebraic Geometry Interpretation and Setup}
We initiate an approach to study the Kuramoto model by 
using an algebraic geometry interpretation of the 
equilibria and studying synchronization. 
Upon fixing $\theta_N = 0$, 
we use the identity 
$\sin(x-y)=\sin x \cos y - \sin y \cos x$
%\fdmargin{mention how you handle the case of the $n$th reference angle. do you set $\theta_{n} = 0$?}
%$$0=\omega_i+\frac{K}{N}\sum_{j=1}^{N} a_{i,j} \left( \sin \theta_i \cos \theta_j - \sin \theta_j \cos \theta_i \right)$$
and substitute $s_i = \sin\theta_i$ and $c_i = \cos\theta_i$
to transform the $N-1$ equations into polynomials
that are coupled by the Pythagorean identity
$s_i^2+c_i^2=1$.
This results in a system of $2(N-1)$ polynomials
in $2(N-1)$ variables:
\begin{equation}\label{Eq:PolySystem}
\begin{array}{l}
0 = \omega_i+\frac{K}{N}\sum_{j=1}^{N} a_{i,j} \left( s_i c_j - s_j c_i \right) \\
0 = s_i^2 + c_i^2 - 1
\end{array}
\end{equation}

Given a fixed integer $N$, our first goal is to find all real 
solutions of polynomial system (\ref{Eq:PolySystem}) for a given choice of 
$\Omega$, $K$, and adjacency matrix $A=[a_{ij}]$. 
We accomplish this using a technique from numerical algebraic geometry 
via parameter homotopies, which we briefly review below.
Once the system is solved, we 
determine which solutions (if any) are stable steady-state 
solutions by analyzing the eigenvalues of the Jacobian at each real solution.

\noindent\textit{Numerical Algebraic Geometry Methods.~-~}
%\section*{Numerical Algebraic Geometry Methods}
For a fixed $N$, we can interpret \eqref{Eq:PolySystem} 
as the system $F_N(s,c;K,\Omega,A) = 0$ with variables
$s$ and $c$ and parameters $K$, $\Omega$, and $A$.
This leads to using parameter homotopies \cite{CoeffParam}
for solving, which is a two-step process.
First, in the {\em ab initio} phase, one computes all solutions for a 
sufficiently random set of parameters.  
Then, in the parameter homotopy phase, one solves for given parameters 
by deforming from the parameters selected in the {\em ab initio} phase.
The {\em ab initio} phase is performed once, while the parameter homotopy phase
is performed for each set of parameters.
The following provides a short introduction with
more details provided in \cite[Ch.~7]{SW05} and \cite[Ch.~6]{BHSW13}.

This process is a generalization of other standard homotopy methods,
such as a \textit{total degree homotopy} or \textit{multi-homogeneous homotopy}
frequently used, e.g.,~\cite{SW05, BHSW13,Mehta:2011wj,Maniatis:2012ex,Kastner:2011zz,Mehta:2012wk,Mehta:2012qr,Greene:2013ida,Mehta:2013fza,MartinezPedrera:2012rs,He:2013yk}.
The method has also recently used in solving the power-flow systems \cite{konsta,chandra_mehta14}.
%The strategy of solving via numerical polynomial homotopy continuation is summarized as follows.
For example, consider solving a system of $m$ polynomial equations in $m$ variables, $f(x)=(f_1(x),\dots,f_m (x))=0$ where $x=(x_1,\dots, x_m)$.
For a total degree homotopy, the space of interest 
is the space of all polynomial systems $g(x)$ of $m$ polynomials in $m$ variables
where $d_i := \deg f_i = \deg g_i$.  Since a random element in this space
has the B\'ezout number of roots, namely $\prod_{i=1}^m d_i$,
the system $G = (G_1,\dots,G_m)=\textbf{0}$ where $G_i(x) = x_i^{d_i} - 1$
is sufficiently random in this space.  

In this case, the {\em ab initio} phase for solving $G = 0$ is trivial.
Then, in the parameter homotopy phase, one can consider the
homotopy $H(x,t) = (1-t)f(x) + \gamma t G(x)$ for a random $\gamma\in\bC$.
Now, the solutions of $H(x,t) = 0$ at $t=1$ are known 
and we want to compute the solutions of $H(x,0) = 0$.  
Using a numerical predictor-corrector method implemented 
in {\tt Bertini} \cite{Bertini}, we track each of the solution paths from $t=1$ to $t=0$. 
The number~$\gamma$ ensures that such path tracking will obtain \textit{all} of the isolated complex solutions of $f(x)=0$, from which the real and nonreal solutions can be identified. 

Now, returning to our case $F_N$, suppose, for simplicity, 
that we want to study the behavior of the solutions 
of $F_N = 0$ as a function of $K$, that is, we also fix $\Omega$ and $A$.  
In the {\em ab initio} phase, we pick a random $K = K_0\in\bC$ 
and compute the solutions $S_0$ of $F_N(s,c;K_0,\Omega,A) = 0$
using a total degree homotopy described above via {\tt Bertini}.

In the parameter homotopy phase, we consider solving for
various choices of parameters $K$, e.g., $K=2,\ldots,100$.  
For each choice of $K$, we simply track the solution paths of
$F_N(s,c;K\cdot (1-t) + t\cdot K_0,\Omega,A) = 0$ starting
at $t = 1$ with solutions $S_0$.  
Each of these computations is relatively inexpensive, thereby making 
it practical to compute the solutions at hundreds of different 
choices for $K$.

The resulting solutions are sorted based on whether they are real or not, and their stability
is investigated as follows. From the  sine $s_i$ and cosine $c_i$ values, we compute
the corresponding angle $\theta_i$.  The eigenvalues of the 
Jacobian of \eqref{Eq:KuramotoStandard} are used 
to analyze stability.  The {\em index} is the number of positive
eigenvalues so that real solutions with index zero correspond
to stable steady states.

%\fdmargin{you should stay general as long as possible}

%With this setup, if we fix $A$, then $K$ is the only unfixed parameter left in the system.
%Generalizing to other distributions of natural frequencies is straightforward. 

\noindent\textit{Results for the complete graph.~-~}%\label{Sec:Results}
%\section*{Results for the complete graph}
We start with the most prominent and well studied case corresponding
to the complete graph, namely $a_{i,j} = 1$.
We select $\Omega=(\omega_1\ldots\omega_N)$ to be $N$ equidistant numbers, 
namely $\omega_i=-1 + (2i -1)/N$.
%\fdmargin{indices correct. they don't seem to add up}
%{\tr namely $\omega_i=-1 + 2(i-1)/N$ for $i \in \{1,\dots,N\}$.}
%Though the techniques presented here apply to arbitrarily complex interaction topologies,
%since the complete graph, i.e., all $a_{i,j}=1$, is the most prominent and well studied case, 
%we begin with this case.
%\fdmargin{what is your choice for the ``reference polynomial'' $g(x)$? here}
%When computing solutions to polynomial system \eqref{Eq:PolySystem} for 
%given values of $N$ and $K$, we are particularly interested in the real 
%solutions.  In Figures \ref{figure:NumEquilibriaSmallK} and
%\ref{figure:NumEquilibriaLargeK}, we summarize the number 
%of real solutions for $N=3,\ldots,18$ at various values of $K$.
Figures~\ref{figure:NumEquilibriaSmallK} and~\ref{figure:NumEquilibriaLargeK}
summarize the number of real-valued solutions to the 
polynomial system~\eqref{Eq:PolySystem} for $N=3,\ldots,18$ at various values of $K$.
%
\begin{figure}[ht]
  \centering
  \includegraphics[width=0.5\textwidth]{figure_NumEquilibriaSmallK.eps}
  \caption{Number of equilibria for the case of equidistant natural frequencies on the complete graph at $1\leq K \leq 2$.}
\label{figure:NumEquilibriaSmallK}
\end{figure}
%
\begin{figure}[ht]
  \centering
  \includegraphics[width=0.5\textwidth]{figure_NumEquilibriaLargeK.eps}
  \caption{Number of equilibria for the case of equidistant natural frequencies on the complete graph at $2\leq K \leq 100$.}
\label{figure:NumEquilibriaLargeK}
\end{figure}
%
%\fdmargin{were the  $328$ and $326$ solutions distinct in this case?}
In general, we see that the number of distinct real solutions tends to grow as $K$ 
increases.  However, there are some surprising exceptions to this expected monotinicity property, e.g., the $N=9$ case has
$328$ real solutions at $K=60$ and $326$ real solutions at $K=70$.  Our 
results suggest that eventually the number of real solutions stabilizes 
to a constant for each $N$. %\fdmargin{as $K$ increases, the number of equilibria equals those in (3) with $\omega_{i}/K = 0$ (which is known to be a bounded constant)}
For $N=18$, we computed 140,356 real solutions when $K=50$; 
computing these without our polynomial formulation would be very difficult.

Once the real solutions are found, we turn our attention to determining 
which ones are stable.  Via Jacobian analysis we find that there 
is exactly a single stable steady state solution for $N=3,\ldots,18$ at each 
$K$ value tested, with the exception of small values of $K$ which 
result in no real solutions.  This result is consistent with theoretic findings for the complete graph \cite{DA-JAR:04,REM-SHS:05,dorfler2011critical,verwoerd2008global}.
{Next, for each $N$ we use the results of our stability analysis 
to determine numerical upper and lower bounds on $K_c(N)$ as follows. 
We acquire an upper bound on $K_c(N)$ by determining 
the smallest tested value of $K$ for which a real steady state is 
found.  Similarly, we find a lower bound on $K_c(N)$ by determining 
the largest tested value of $K$ for which no real steady states 
are found.}
These results are presented in Figure \ref{figure:EquidistantK_c}. 
We verify that these bounds for $K_c(N)$ are consistent 
with the known explicit bounds from \cite[Corollary 6.7]{FD-FB:13b}. 
%{\tb  We also observe a few inconsistent instances for $N \in \{2,3\}$ as a result of a finite numeric resolution of tested values for $K_c(N)$.}
%{\tr 
%In regard to the above sentence:
%I'm not sure if \textit{inconsistent} is the right word to use here. 
%For example at $N=3$, the explicit bounds are $1.1363 \leq K_c \leq 1.3333$ while the numeric bounds are $1.10 \leq K_c \leq 1.15$. 
%I don't view this as an \textit{inconsistency}, because the true value could be $K_c=1.14$ (for example), which agrees with both bounds.
%I would prefer to say something such as: 
%\begin{quote} 
For $N \in \{3,4\}$, our computed lower bound is less than the 
explicit lower bound due to the coarse resolution of tested numeric values for $K_c(N)$.
%\end{quote}
%}

%\fdmargin{In Figure 3 and 8, it is not clear why we report \underline{two} numerical values: a lower bound and an upper bound. Shouldn't there be only one numeric number? It would be helpful if we drop a sentence or two regarding how we obtain numeric lower and upper bounds in Fig. 3 and 8.}

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_EquidistantK_c.eps}
\caption{Bounds for $K_c(N)$ in the case of equidistant natural frequencies on the complete graph.  Known explicit bounds from \cite[Corollary 6.7]{FD-FB:13b} are shown for comparison.}
\label{figure:EquidistantK_c}
\end{figure}

In Ref.~\cite{AA-SS-VP:81}, it is conjectured that if there is a stable 
equilibrium, then there is a stable equilibrium satisfying 
$|\theta_{i}-\theta_{j}| < \pi/2$ for all neighbors $\{i,j\}$. 
We checked this conjecture and we quickly found 
counterexamples. 
For $N=3$ and $K=1.15$, there is one stable 
equilibria at 
$(\theta_1,\theta_2,\theta_3)\approx (1.6820,0.8410,0)$.
Since $\theta_1$ and~$\theta_3$ are coupled with
$|\theta_{1}-\theta_{3}| > \pi/2$, we reject the conjecture.

For larger $N$, the unique stable equilibria can have angles
even more spread out, such as
$(\theta_{1},\ldots,\theta_{15})\approx (2.0981,1.8867,\ldots,0.2114,0)$ 
when $N=15$ and $K=1.3$.
%\fdmargin{nice result. i could also analytically prove it when i saw it}
These interesting numerical observations prompt us to analytically reject the conjecture in \cite{AA-SS-VP:81}: in \cite{FD-FB:13b} it is shown that the worst-case equilibrium configuration for $K=K_{c}$ and $n>3$ corresponds to a tripolar distribution of all angles spanning a half-circle. From \cite{DA-JAR:04,REM-SHS:05}, we know that this is the only stable equilibrium. These two arguments suffice to reject the conjecture in \cite{AA-SS-VP:81}.

As mentioned earlier, our stability analysis is 
based on computing the index corresponding to each real solution.  
Figure~\ref{figure:IndexPlotN3To15} shows a histogram of 
these values for the case $K=100$.
\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_IndexPlot.eps}
\caption{Number of equilibria with given index for equidistant natural frequencies and a complete graph with $K=100$.}
\label{figure:IndexPlotN3To15}
\end{figure}
We observe that when $j$ is small enough relative to $N$, the number 
of real solutions with index $j$ is exactly $\binom{N}{j}$.  
In particular, 
for $3\leq N\leq 4$ this behavior occurs for $j\leq 1$ and,
for $5\leq N\leq 18$, this behavior occurs for 
$j\leq \lceil{\frac{2}{5}(N-1)}\rceil$.
Based on these results, we conjecture that this phenomenon occurs in general: for equidistant natural frequencies, 
if $0 \!\leq\! j \!\ll\! N$ and $K\gg 0$, then the number of equilibria with index $j$ is 
expected to be $\binom{N}{j}$.
This conjecture can be used when $N$ may be 
too large to compute all solutions.  
For example, we expect 75,287,520 real equilibria of 
index $5$ for $N = 100$.

\noindent\textit{Results for cyclic graphs.~-~}%\label{Sec:ResultsCyclic}
%\section*{Results for the cyclic graph}
Although our investigation initially focused on the complete graph being the most well studied case, we also performed a preliminary 
analysis of the coupling arrangement defined by an undirected cyclic graph. 
Cyclic graphs are known for having multi-stable equilibria \cite{FD-FB:13b}, entirely unstable equilibrium landscapes \cite{AA-SS-VP:81}, and thus a quite distinct behavior from acyclic or complete graphs \cite{FD-MC-FB:11v-pnas,wiley2006size,ermentrout1985behavior}.

For $N=10$, we used the same equidistant natural frequencies mentioned 
earlier with Figure \ref{figure:CyclicGraph} showing the number of equilibria and 
the number of stable equilibria at each integer $K=0,\ldots,100$.

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_CyclicGraph.eps}
\caption{Number of equilibria and number of stable equilibria for a cyclic coupling arrangement when $N=10$.}
\label{figure:CyclicGraph}
\end{figure}

For some values of $K$, the system possesses only unstable equilibria.
In particular, when $K$ is $13$, $14$, or $15$, there are $76$, $164$, and $260$ equilibria, respectively, all of which are unstable.
Since the critical points of the potential \eqref{eq: potential} correspond to power flow in a transmission network, a very interesting technological implication of these results is that there are power demands which can be met only with unstable equilibria. 
We find one stable equilibria at $K=16$ and conclude $15\leq K_c(10)\leq 16$. 
Figure~\ref{figure:CyclicGraph} also indicates multistability for some values of $K$, 
with at most {\em three} stable equilibria for each investigated sample. 
These results are in contrast to the complete graph case, in which exactly one stable equilibria exists whenever the system has real-valued solutions. 
%{\tb {\bf Comment:} the first result is very nice and surprising! The second multi-stability result is not so surprising. Multi-stability for $N>5$ is well known. You can make the finding more interesting if you (for fixed $N$) slowly increase $K$ and count the number of equilibria as well as the number of stable equilibria. I am not even sure of these curves should be monotone (??). I would like to see them.}

Next, we take a closer look at the geometric configuration of angles occurring at the stable equilibria shown in Figure \ref{figure:CyclicGraph}.
For relatively small $K$ such as $K=16,17,18$, the configuration of angles at stable equilibria shows no discernible structure.
For $K=35$ there are two stable equilibria, with one exhibiting phase sync, i.e., angles clustered around $0$, and one exhibiting a splay state, 
i.e., angles approximately uniformly distributed on $[0,2\pi)$.
For the splay state, the $\theta_i$ decrease on $[0,2\pi)$.

%\fdmargin{what do you mean by ``decrease'' and ``increase''? Do you mean ``sorted in decreasing/increasing order according to the natural frequencies'' ?}
For each $K=36,\ldots,100$ there are three stable equilibria, with each case consisting of one phase sync and two splay states.
%\redmargin{I tried to make this sentence clearer.}
In one of these splay states $\theta_1,\ldots,\theta_9$ are arranged in increasing order on $[0,2\pi)$, and in the other $\theta_9,\ldots,\theta_1$ are arranged in decreasing order on $[0,2\pi)$.
In Figure \ref{figure:CyclicConfiguration}, we depict the three stable equilibria at $K=100$. 
As $K$ increases, the steady state phase sync gradually becomes more tightly clustered, with the angle range decreasing from $\sim0.9432$ for $K=35$ to $\sim0.3256$ for~$K=100$. 
This is in accordance with the asymptotic result that exact phase sync is a critical point of the Kuramoto potential~\eqref{eq: potential} as $K\to \infty$ \cite{FD-FB:13b}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_CyclicConfiguration.eps}
\caption{Configuration of $\theta_1,\ldots,\theta_{10}$ for the three stable~equilibria occurring at $K=100$ for a cyclic graph with $N = 10$.}
\label{figure:CyclicConfiguration}
\end{figure}


\noindent\textit{Results for random graphs.~-~}%\label{Sec:Conclusion}
%\section*{Results for random graphs}
We now turn our attention to random coupling arrangements.
In these computations, we choose each $a_{i,j}$ as $0$ or $1$ according to a pre-determined probability $P$ while setting $a_{j,i}=a_{i,j}$.
In other words, $[a_{i,j}]$ is the adjacency matrix of a symmetric Erd\"{o}s-R\'enyi random graph with coupling probability~$P$, where we restrict ourselves to connected graphs.
For these numerical experiments, we use the same equidistant natural frequencies discussed\,earlier.

%\fdmargin{Well, I guess you didn�t do that, right? You probably constructed a bunch of random graphs and ordered them according to the cycles, right? How did you select exactly 100 for each case? I think this needs some clarifications.}
%\redmargin{I updated Figure 7 (previously we had 0-10 cycles) and added this clarifying sentence.}
First, we fix $N=8$ and $K=100$ and investigate the number of 
stable equilibria for sparse random graphs.
For each $c=0,1,\ldots,20$, we construct $100$ random graphs having exactly $c$ cycles. 
{We achieve this by generating adjacency matrices of symmetric Erd\H{o}s-R\'enyi random graphs until we have $100$ instances of connected graphs with the desired number of cycles.}
Depending on the graph, we find $1$, $2$, or $3$ stable equilibria.
The averages of these results are shown in Figure~\ref{figure:RandomGraphCyclesvsNumStable}
showing that the sparsest (i.e., acyclic) as well as sufficiently dense graphs have exactly one stable equilibrium confirming the analytic results for these two extremal  topologies \cite{FD-MC-FB:11v-pnas,FD-FB:13b,DA-JAR:04,REM-SHS:05,taylor2012there}. Since cycles are known to display multi-stability (see Figure~\ref{figure:CyclicGraph}), it is to be expected that there is more than one stable equilibrium for small cycle numbers.  Quite surprisingly, multiple distinct extrema can be observed in the stable equilibrium distribution in Figure~\ref{figure:RandomGraphCyclesvsNumStable}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_RandomGraphCyclesvsNumStable.eps}
\caption{Number of stable equilibria for random graphs according to number of cycles in the case of $N=8$ and $K=100$.}
\label{figure:RandomGraphCyclesvsNumStable}
\end{figure}

Next, we fix $N=10$ and investigate how $K_c(N)$ depends on the density of the graph. 
For each coupling probability $P=0.25,0.375,\ldots,0.875$, we generated $100$ random graphs. 
For each graph, we computed equilibria and determine stability at the following values of $K$: 
for $P=0.25$ we use $K=1,\ldots,20$; for $P=0.375$ we use $K=1,\ldots,15$; for $P=0.5$ we use $K=1,\ldots,10$; and for each $P=0.625,0.75,0.875$, we use $K=1,\ldots,5$.
In all cases, we find a value of $K$ such that at least one stable equilibrium occurs at $K$ while no stable equilibria occur at $K-1$, thereby allowing us to estimate bounds on $K_c(10)$. 
Figure~\ref{figure:RandomGraphCyclesvsK_c} shows these results sorted according 
to the number of cycles in the graph.

For undirected and connected graphs,
a theoretic lower and a conjectured theoretic upper bound are \cite{FD-FB:13b,FD-MC-FB:11v-pnas} 
$$N\cdot \max_i \left\{ \frac{|\omega_i|}{\textup{deg}_{i}} \right\} \leq K_c(N)\leq N\cdot \| B^T L^\dagger \omega \|_\infty,$$
where $\textup{deg}_{i} = \sum_j a_{ij}$ is the degree of node $i$, $B$ is the oriented incidence matrix, and $L$ is the network Laplacian matrix. 
These bounds are shown in Figure \ref{figure:RandomGraphCyclesvsK_c} for comparison. 
We verified that our $600$ individual results as well as the averages shown are consistent with these theoretic bounds and validate their accuracy. 
Figure~\ref{figure:RandomGraphCyclesvsK_c} indicates that the numerical upper bound is tighter than the theoretic bound on average, while the numerical lower bound tends to be weaker. 
Since the tightness of our computed bounds depends on the resolution of $K$ values tested, one could check more refined $K$ values to obtain tighter numerical bounds for particular graphs of interest.

%\begin{figure}[ht]
%\centering
%\includegraphics[width=0.5\textwidth]{figure_RandomGraphPvsK_c.eps}
%\caption{Average bounds on $K_c(N)$ for different coupling probabilities in the case of $N=10$.  Computations were performed for $100$ random graphs at each value of $P$.  For the $P=1$ result, we reuse data from the complete graph computations discussed earlier.}
%\label{figure:RandomGraphPvsK_c}
%\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{figure_RandomGraphCyclesvsK_c.eps}
\caption{Average bounds on $K_c(N)$ for random graphs according to number of cycles in the case of $N=10$.}
\label{figure:RandomGraphCyclesvsK_c}
\end{figure}

%{\tb Comment: these are interesting experiments even for small $N$. May I add a few suggestions to get more spectacular results:

%\begin{enumerate}

%	\item fix $K$ (pick a large number), fix $N$ (as large as your computation capabilities allow), and fix $\omega$. Consider the coupling probability $p = \textup{Prob}(a_{ij}=a_{ji}=1)$ as a variable. For $p=0$ the graph is disconnected and for $p = 1$, we have a complete graph. Now construct a lot of samples while varying $p$ from a small value, say $p=0.2.5$ (important: the graph needs to be connected) to $p=1$ and count the number of stable equilibria. I conjecture that the plot will be concave, with $1$ stable equilibrium for small $p$ (for which we will likely have a very sparse and acyclic graph), $1$ stable equilibrium for $p=1$ (case of a complete graph), and multiple stable equilibria in between. This would be really interesting to see.

%	\item same experiment as before for fixed $N$ and $\omega$. Vary $p$ from $p=0.25$ to $p=1$ and try to find the critical value $K_{c}$. For $p=0.25$ (likely acyclic graph) and $p=1$ (complete graph) you should get $K_{c} \approx \| B^{T} L^{\dagger} \omega\|_{\infty}$, where $L^{\dagger}$ is the Moore-Penrose pseudoinverse of the network Laplacian $L$ and $B$ is the incidence matrix. This is more or less a conjectured result in \cite{FD-MC-FB:11v-pnas}.

%\end{enumerate}
%}



In future work, we aim to extend our investigations to directed graphs, negative weights, 
cosine coupling, and statistical properties of random graph models, as well as verify conjectures for variations of the Kuramoto models
as studied in Ref.~\cite{wiley2006size}.

\noindent\textit{Acknowledgement.~-~}
%\section*{Acknowledgement} 
DM would like to thank Carlo Laing and Steven Strogatz for their feedback at the initial stages of this work. 
NSD, JDH, and DM were supported by DARPA Young Faculty Award
with NSD and JDH additionally supported by NSF DMS-1262428.
FD was supported by in part by ETH Z\"urich startup funds.

\bibliographystyle{unsrt}
\input{kuramoto_bib.bbl}
%\bibliography{kuramoto_references}


%\end{multicols} % open this to use the two-column format

\end{document}

kuramoto_bib.bbl
\begin{thebibliography}{10}

\bibitem{kuramoto1975self}
Y. Kuramoto.
%\newblock Self-entrainment of a population of coupled non-linear oscillators.
\newblock In {\em International symposium on mathematical problems in
  theoretical physics}, pages 420--422. Springer, 1975.

\bibitem{acebron2005kuramoto}
J.A. Acebr{\'o}n, L.L. Bonilla, C.J.~P{\'e}rez Vicente, F. Ritort, and R. Spigler.
%\newblock The {K}uramoto model: A simple paradigm for synchronization phenomena.
\newblock {\em Reviews of modern physics}, 77(1):137, 2005.

\bibitem{strogatz2000kuramoto}
S.H. Strogatz.
%\newblock From {K}uramoto to {C}rawford: exploring the onset of synchronization in populations of coupled oscillators.
\newblock {\em Physica D: Nonlinear Phenomena}, 143(1):1--20, 2000.

\bibitem{FD-FB:13b}
F.~D{\"o}rfler and F.~Bullo.
%\newblock Synchronization in complex oscillator networks: A survey.
\newblock {\em Automatica}, 50(6):1539--1564, June 2014.

\bibitem{Casetti:June2003:0022-4715:1091}
L.~Casetti, M.~Pettini, and E.G.D. Cohen.
%\newblock Phase transitions and topology changes in configuration space.
\newblock {\em Journal of Statistical Physics}, 111:1091--1123(33), June 2003.

\bibitem{Mehta:2009}
D.~Mehta.
\newblock {Lattice vs. Continuum: Landau Gauge Fixing and 't Hooft-Polyakov
  Monopoles}.
\newblock {\em Ph.D. Thesis, The Uni. of Adelaide, Australasian Digital Theses
  Program}, 2009.

\bibitem{Mehta:2010pe}
D. Mehta and M. Kastner.
%\newblock {Stationary point analysis of the one-dimensional lattice Landau gauge fixing functional, aka random phase XY Hamiltonian}.
\newblock {\em Annals Phys.}, 326:1425--1440, 2011.

\bibitem{vonSmekal:2007ns}
L. von Smekal, D. Mehta, A. Sternbeck, and A.G. Williams.
%\newblock {Modified Lattice Landau Gauge}.
\newblock {\em PoS}, LAT2007:382, 2007.

\bibitem{vonSmekal:2008es}
L. von Smekal, A. Jorkowski, D. Mehta, and A. Sternbeck.
%\newblock {Lattice Landau gauge via Stereographic Projection}.
\newblock {\em PoS}, CONFINEMENT8:048, 2008.

\bibitem{Nerattini:2012pi}
R. Nerattini, M. Kastner, D. Mehta, and L. Casetti.
%\newblock {Exploring the energy landscape of XY models}.
\newblock {\em Phys.Rev.}, E87(3):032140, 2013.

\bibitem{kastner2011stationary}
M. Kastner.
%\newblock Stationary-point approach to the phase transition of the classical xy chain with power-law interactions.
\newblock {\em Physical Review E}, 83(3):031114, 2011.

\bibitem{Hughes:2012hg}
C. Hughes, D. Mehta, and J.-I. Skullerud.
%\newblock {Enumerating Gribov copies on the lattice}.
\newblock {\em Annals Phys.}, 331:188--215, 2013.

\bibitem{Mehta:2014jla}
D. Mehta and M. Schr\"ock.
%\newblock Enumerating copies in the first gribov region on the lattice in up to four dimensions.
\newblock Phys.~Rev.~D \textbf{89} 094512 (2014).

\bibitem{Mehta:2009zv}
D. Mehta, A. Sternbeck, L. von Smekal, and A.G. Williams.
%\newblock {Lattice Landau Gauge and Algebraic Geometry}.
\newblock {\em PoS}, QCD-TNT09:025, 2009.

\bibitem{JB-CIB:82}
J.~Baillieul and C.I. Byrnes.
%\newblock Geometric critical point analysis of lossless power system models.
\newblock {\em IEEE Transactions on Circuits and Systems}, 29(11):724--737,
  1982.

\bibitem{AA-SS-VP:81}
A.~Araposthatis, S.~Sastry, and P.~Varaiya.
%\newblock Analysis of power-flow equation.
\newblock {\em International Journal of Electrical Power \& Energy Systems},
  3(3):115--126, 1981.

\bibitem{strogatz1991stability}
S.H. Strogatz and R.E. Mirollo.
%\newblock Stability of incoherence in a population of coupled oscillators.
\newblock {\em Journal of Statistical Physics}, 63(3-4):613--635, 1991.

\bibitem{verwoerd2009computing}
M. Verwoerd and O. Mason.
%\newblock On computing the critical coupling coefficient for the {K}uramoto model on a complete bipartite graph.
\newblock {\em SIAM Journal on Applied Dynamical Systems}, 8(1):417--453, 2009.

\bibitem{verwoerd2008global}
M. Verwoerd and O. Mason.
%\newblock Global phase-locking in finite populations of phase-coupled oscillators.
\newblock {\em SIAM Journal on Applied Dynamical Systems}, 7(1):134--160, 2008.

\bibitem{dorfler2011critical}
F. D{\"o}rfler and F. Bullo.
%\newblock On the critical coupling for kuramoto oscillators.
\newblock {\em SIAM Journal on Applied Dynamical Systems}, 10(3):1070--1099,
  2011.

\bibitem{DA-JAR:04}
D.~Aeyels and J.A. Rogge.
%\newblock Existence of partial entrainment and stability of phase locking behavior of coupled oscillators.
\newblock {\em Progress on Theoretical Physics}, 112(6):921--942, 2004.

\bibitem{REM-SHS:05}
R.E. Mirollo and S.H. Strogatz.
%\newblock The spectrum of the locked state for the {K}uramoto model of coupled oscillators.
\newblock {\em Physica D: Nonlinear Phenomena}, 205(1-4):249--266, 2005.

\bibitem{FD-MC-FB:11v-pnas}
F.~D{\"o}rfler, M.~Chertkov, and F.~Bullo.
%\newblock Synchronization in complex oscillator networks and smart grids.
\newblock {\em Proceedings of the National Academy of Sciences},
  110(6):2005--2010, 2013.

\bibitem{taylor2012there}
R. Taylor.
%\newblock There is no non-zero stable fixed point for dense networks in the homogeneous {K}uramoto model.
\newblock {\em Journal of Physics A: Mathematical and Theoretical},
  45(5):055102, 2012.

\bibitem{ochab2009synchronization}
J.~Ochab and P.F.~G{\'o}ra.
\newblock Synchronization of coupled oscillators in a local one-dimensional {K}uramoto model.
\newblock {\em arXiv preprint arXiv:0909.0043}, 2009.

\bibitem{tilles2011multistable}
P.~F.~C.~Tilles, F.~F.~Ferreira and H.~A.~Cerdeira.
%\newblock Multistable behavior above synchronization in a locally coupled Kuramoto model.
\newblock Phys.~Rev.~E \textbf{83} 066206 (2011).

\bibitem{CoeffParam}
A.P. Morgan and A.J. Sommese.
%\newblock Coefficient-parameter polynomial continuation.
\newblock {\em Appl. Math. Comput.}, 29(2):123--160, 1989.

\bibitem{SW05}
A.J. Sommese and C.W. Wampler.
\newblock {\em The Numerical Solution of Systems of Polynomials Arising in
  Engineering and Science}.
\newblock World Scientific Publishing, Hackensack, NJ, 2005.

\bibitem{BHSW13}
D.J. Bates, J.D. Hauenstein, A.J. Sommese, and C.W. Wampler.
\newblock {\em Numerically solving polynomial systems with Bertini}, volume~25.
\newblock SIAM, 2013.

\bibitem{konsta}
D.~{Mehta}, H.~{Nguyen}, and K.~{Turitsyn}, ``{Numerical Polynomial Homotopy
  Continuation Method to Locate All The Power Flow Solutions},'' arXiv:1408.2732 (2014).

\bibitem{chandra_mehta14}
S.~Chandra, D.~Mehta, and A.~Chakrabortty, ``Exploring the impact of wind
  penetration on power system equilibrium using a numerical continuation
  approach,'' arXiv:1409.7844 (2014).

\bibitem{Mehta:2011wj}
D. Mehta.
%\newblock {Numerical Polynomial Homotopy Continuation Method and String Vacua}.
\newblock {\em Adv.High Energy Phys.}, 2011:263937, 2011.

\bibitem{Maniatis:2012ex}
M. Maniatis and D. Mehta.
%\newblock {Minimizing Higgs Potentials via Numerical Polynomial Homotopy Continuation}.
\newblock {\em Eur.Phys.J.Plus}, 127:91, 2012.

\bibitem{Kastner:2011zz}
M. Kastner and D. Mehta.
%\newblock {Phase Transitions Detached from Stationary Points of the Energy Landscape}.
\newblock {\em Phys.Rev.Lett.}, 107:160602, 2011.

\bibitem{Mehta:2012wk}
D. Mehta, Y.-H. He, and J.D. Hauenstein.
%\newblock {Numerical Algebraic Geometry: A New Perspective on String and Gauge Theories}.
\newblock {\em JHEP}, 1207:018, 2012.

\bibitem{Mehta:2012qr}
D. Mehta, J.D. Hauenstein, and M. Kastner.
%\newblock Energy-landscape analysis of the two-dimensional nearest-neighbor $\phi^{4}$ model.
\newblock {\em Phys. Rev. E}, 85:061103, Jun 2012.

\bibitem{Greene:2013ida}
B. Greene, D. Kagan, A. Masoumi, D. Mehta, E.J. Weinberg, and Xiao Xiao.
%\newblock {Tumbling through a landscape: Evidence of instabilities in high-dimensional moduli spaces}.
\newblock {\em Phys.Rev.}, D88(2):026005, 2013.

\bibitem{Mehta:2013fza}
D. Mehta, D.A. Stariolo, and M. Kastner.
%\newblock {Energy landscape of the finite-size spherical three-spin glass model}.
\newblock {\em Phys.Rev.}, E87(5):052143, 2013.

\bibitem{MartinezPedrera:2012rs}
D. Martinez-Pedrera, D. Mehta, M. Rummel, and A. Westphal.
%\newblock {Finding all flux vacua in an explicit example}.
\newblock {\em JHEP}, 1306:110, 2013.

\bibitem{He:2013yk}
Y.-H. He, D. Mehta, M. Niemerg, M. Rummel, and A. Valeanu.
%\newblock {Exploring the Potential Energy Landscape Over a Large Parameter-Space}.
\newblock {\em JHEP}, 1307:050, 2013.

\bibitem{Bertini}
D.J. Bates, J.D. Hauenstein, A.J. Sommese, and C.W. Wampler.
\newblock Bertini: Software for numerical algebraic geometry.
\newblock Available at \url{http://bertini.nd.edu}.

\bibitem{wiley2006size}
%{The size of the sync basin},
D.~A.~Wiley, S.~H.~Strogatz and M.~Girvan.
\newblock Chaos: An Interdisciplinary J. of Non. Sci., \textbf{16} 015103 (2006).

\bibitem{strogatz1988phase}
%Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random intrinsic frequencies},
S.~H.~Strogatz and R.~E.~Mirollo. 
\newblock Physica D: Non. Phen., \textbf{31}, 143 (1988).

\bibitem{ermentrout1985behavior}
%The behavior of rings of coupled oscillators
G.~B.~Ermentrout.
\newblock J. Math.~Bio.,\textbf{23} 55 (1985).

\end{thebibliography}
