{\bm \epsilon}gin{document}

\title{Subcritical Bifurcations in the Kuramoto Model with Positive Coupling}
\author{Timothy Ferguson$^1$}

{\bm \epsilon}gin{abstract}
The Kuramoto model is a standard model for the dynamics of coupled oscillator networks. In particular, it is used to study long time behavior such as phase-locking where all oscillators rotate at a common frequency with fixed angle differences. It has been observed that phase-locking occurs if the natural frequencies of the oscillators are sufficiently close to each other, and doesn't if they aren't. Intuitively, the closer the natural frequencies are to each other the more stable the system is. In this paper, we study the Kuramoto model on ring networks with positive coupling between the oscillators and derive a criterion for a bifurcation in which two branches of phase-locked solutions collide as we increase the coupling strengths. Furthermore, we state a criterion for when one of these branches consists of stable phase-locked solutions. In this case stability is lost as the positive coupling is increased. We then apply our criteria to show that for any size of the ring network there always exists choices of the natural frequencies and coupling strengths so that this bifurcation occurs. (We require at least five oscillators for one branch to consist of stable phase-locked solutions.) Finally, we conjecture that generically our bifurcation is locally a subcritical bifurcation and globally an $S$-curve. (In the case where one branch consists of phase-locked solutions, an $S$-curve results in bistability i.e. the existence of two distinct stable phase-locked solutions.) Finally, we note that our methods are constructive.
\end{abstract}

\maketitle

\noindent \footnotesize\textbf{Keywords.} Kuramoto model, fixed points, bifurcations, stability
\\
\\

\noindent \footnotesize\textbf{AMS subject classifications.} 34C15, 34C23, 34D20

\footnotetext[1]{Department of Mathematics, Arizona State University, 901 S. Palm Walk, Tempe, AZ 85281, USA}

\section{Introduction}

The Kuramoto model is a model for the dynamics of oscillator networks, or more generally, networks of involving periodic components. In particular, it is used to study synchronization whereby the individual oscillators act in unison and rotate at a common frequency rather than rotate at their initial individual frequencies. The Kuramoto model was first defined by Y. Kuramoto in 1984 in \cite{MR762432} and has many practical applications to periodic phenomena. Examples from physics and engineering include generators in the power grid \cite{Filatrella2008} and laser arrays \cite{PhysRevA.52.4089} while examples from biology include pacemaker cells in the heart \cite{peskin75} and flashing colonies of fireflies \cite{Buck1988SynchronousRhythmicFlashing}. For a more complete discussion of synchronization as well as a history of the development of the Kuramoto model see \cite{Sync.book,MR1783382,Acebron.etal.05}. The Kuramoto model on a ring network with $n$ oscillators is the system of $n$ coupled non-linear differential equations
{\bm \epsilon}gin{align} \label{eq:model}
\frac{d\theta_i}{dt} = \omega_i + \sigma ({\bf g}amma_{i+1,i} \sin(\theta_{i+1} - \theta_i) + {\bf g}amma_{i-1,i} \sin(\theta_{i-1} - \theta_i)) \quad \text{for} \quad i \in \{1,\dots,n\}
\end{align}
where $\theta_i$ denotes the phase angle of the $i$th oscillator, $\omega_i$ the natural frequency of the $i$th oscillator, and ${\bf g}amma_{i-1,i} = {\bf g}amma_{i,i-1} > 0$ represents strength of the coupling between $\theta_{i-1}$ and $\theta_i$. (Throughout we will always interpret indices modulo $n$.) Note that the right hand side of \eqref{eq:model} naturally breaks up into two contributions. The first contribution is the natural frequency $\omega_i$ which determines the behavior of $\theta_i$ without any coupling. The second contribution is the term ${\bf g}amma_{i+1,i} \sin(\theta_{i+1} - \theta_i) + {\bf g}amma_{i-1,i} \sin(\theta_{i-1} - \theta_i)$ which determines the behavior of $\theta_i$ due to its interaction with $\theta_{i-1}$ and $\theta_{i+1}$. The parameter $\sigma > 0$ allows us to change the relative strengths of these two contributions. It is well-known that if $\sigma$ is sufficiently large, then the Kuramoto model \eqref{eq:model} has a unique stable fixed point to which all solutions converge. In this case all oscillators rotate at a common rate and we say that the system synchronizes. For results concerning the critical value of $\sigma$ for the onset of synchronization as well as estimates for specific network structures see \cite{Verwoerd.Mason.2008,Dorfler.Chertkov.Bullo.2013,Dorfler.Bullo.2011,Chopra.Spong.2009,Verwoerd.Mason.2009}. For smaller values of $\sigma$ clusters of the network may synchronize to different frequencies. See \cite{Brady.2018} for diagrams illustrating this phenomenon. In particular, these diagrams illustrate bifurcations that occur when two or more synchronized clusters become a single larger synchronized cluster.

It is natural to ask what kind of bifurcations can occur in the Kuramoto model as we change the parameter $\sigma$. In particular, we will focus on bifurcations involving phase-locked solutions i.e. solutions for which all oscillators rotate at a common frequency with fixed angle differences. In Figure \ref{fig:1a} we give an example plot of all phase-locked solutions as we vary $\sigma$. In order to visualize these plots we represent a phase-locked solution by its order parameter $r = \frac{1}{n} \bigg\rvert \sum_{i=1}^n e^{\theta_i \sqrt{-1}}\bigg\rvert$ which is a measure of the solutions level of synchronization. One might expect that increasing $\sigma$ will only create rather than destroy stability of phase-locked solutions. However, this is not always the case. For example, Terada, Ito, Aoyagi, and Yamaguchi \cite{Terada_2017} numerically discovered a subcritical saddle-node bifurcation whereby a stable phase-locked solution collides with a 1-saddle. They did this by choosing the natural frequencies from a bimodal distribution. In fact, their saddle-node bifurcation is part of a larger object called an $S$-curve which is discussed in \cite{golubitsky1985singularities} and illustrated in Figure \ref{fig:1b}. (Note that for values of $\sigma$ slightly less than one the $S$-curve in Figure \ref{fig:1b} results in bistability, namely, the existence of two distinct stable phase-locked solutions.)

{\bm \epsilon}gin{figure}[H]
\captionsetup{font=scriptsize}
\centering
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Example.pdf}
  \caption{} \label{fig:1a}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{CloseUp.pdf}
  \caption{} \label{fig:1b}
\end{subfigure}
\caption{(A) A plot of all phase-locked solutions of the Kuramoto model \eqref{eq:model} with $(\omega_1,\omega_2,\omega_3) = (1,0,-1)$ and $({\bf g}amma_1,{\bf g}amma_2,{\bf g}amma_3) = (1,1,1)$ for a range of values of $\sigma$. A phase-locked solution $(\theta_1,\theta_2,\theta_3)$ is represented by its order parameter $r$. (B) A plot of phase-locked solutions of the Kuramoto model \eqref{eq:model} with $(\omega_1,\omega_2,\omega_3,\omega_4,\omega_5) = (-0.69,0.69,0,0,0)$ and $({\bf g}amma_1,{\bf g}amma_2,{\bf g}amma_3,{\bf g}amma_4,{\bf g}amma_5) = (0.20,1,1,1,1)$ for a range of values of $\sigma$. The plot contains a local subcritical saddle-node bifurcation at $\sigma = 1$ which is a part of a global $S$-curve. For values of $\sigma$ slightly less than one there are two distinct stable phase-locked solutions.} \label{fig:1}
\end{figure}

In Section \ref{sec:notation} we give a criterion for the Kuramoto model \eqref{eq:model} to have a bifurcation in which two branches of phase-locked solutions collide as $\sigma$ increases. We conjecture that this bifurcation is generically a subcritical saddle-node bifurcation. In addition, we give a sufficient condition for one of these branches to consist of stable phase-locked solutions. We end the section by applying this criterion to prove that for every $n {\bf g}e 3$ there exists values of ${\bm \omega}$ and ${\bf g}ammab$ for which the Kuramoto model \eqref{eq:model} exhibits this bifurcation. We also show for $n {\bf g}e 5$ that this can be done so that one branch consists of stable phase-locked solutions. We end this section by presenting our methods as a constructive algorithm.

In Section \ref{sec:future} we discuss additional questions and future work. We defer proofs to the appendix.

\section{Notation and Main Results} \label{sec:notation}

In this section we introduce all relevant notations and definitions as well as state our main results. In general we defer the proofs of our results to the appendix.

Because we are only considering ring networks we let $e_i = (i,i+1)$ denote the $i$th edge in the network, $\eta_i = \theta_{i+1} - \theta_i$ the angle difference along the $i$th edge, and ${\bf g}amma_i = {\bf g}amma_{i+1,i}$ the coupling strength along the $i$th edge. For simplicity we further define the vector quantities ${\bm \omega} = (\omega_1, \dots, \omega_n)$, ${\bf g}ammab = ({\bf g}amma_1, \dots, {\bf g}amma_n)$, ${\bm \theta} = (\theta_1, \dots, \theta_n)$, and ${\bm \eta} = (\eta_1, \dots, \eta_n)$. With this notation we can rewrite \eqref{eq:model} as
{\bm \epsilon}gin{align} \label{eq:model2}
\frac{d\theta_i}{dt} = \omega_i + \sigma ({\bf g}amma_i \sin \eta_i - {\bf g}amma_{i-1} \sin \eta_{i-1}) \quad \text{for} \quad i \in \{1,\dots,n\}.
\end{align}
This leads to the following definitions.
{\bm \epsilon}gin{definition} \label{def:gJp}
Define the vector-valued function ${\bf g}b(\cdot,{\bf g}ammab) : \mathbb{R}^n \rightarrow \mathbb{R}^n$ by
{\bm \epsilon}gin{align*}
{\bf g}b({\bm \theta},{\bf g}ammab)_i := -{\bf g}amma_i \sin \eta_i + {\bf g}amma_{i-1} \sin \eta_{i-1}
\end{align*}
and define the matrix-valued Jacobian ${\bf J}b(\cdot,{\bf g}ammab) : \mathbb{R}^n \rightarrow \mathbb{R}^{n \times n}$ by
{\bm \epsilon}gin{align*}
{\bf J}b({\bm \theta},{\bf g}ammab)_{ij} := \frac{\partial}{\partial \theta_j} {\bf g}b({\bm \theta},{\bf g}ammab)_i =
{\bm \epsilon}gin{cases}
-{\bf g}amma_i \cos \eta_i & \mbox{if $j = i+1$,}
\\
{\bf g}amma_i \cos \eta_i + {\bf g}amma_{i-1} \cos \eta_{i-1}  & \mbox{if $j = i$,}
\\
-{\bf g}amma_{i-1} \cos \eta_{i-1} & \mbox{if $j = i-1$,}
\\
0 & \mbox{otherwise.}
\end{cases}
\end{align*}
Furthermore, let ${\bf P}$ be an $(n-1) \times n$ matrix with rows which form an orthonormal basis for the orthogonal complement ${\bf 1}^\perp$ where ${\bf 1} = (1,\dots,1)^\top$.
\end{definition}

With this notation we can rewrite \eqref{eq:model} as the vector equation
{\bm \epsilon}gin{align}
\frac{d{\bm \theta}}{dt} = {\bm \omega} - \sigma {\bf g}b({\bm \theta},{\bf g}ammab).
\end{align}
Note that we chose the sign of ${\bf g}b({\bm \theta},{\bf g}ammab)$ to simplify computations later on. Now a phase-locked solution has the form ${\bm \theta}(t) = {\overline{\omega}} {\bf 1} t + {\bm \theta}_0$ where ${\overline{\omega}} = \frac{1}{n} \sum_{i=1}^n \omega_i$, the mean of ${\bm \omega}$, is the common frequency of rotation and ${\bm \theta}_0$ represents the fixed angle differences. By rotating the reference frame we can assume that ${\overline{\omega}} = 0$ without loss of generality. Therefore a phase-locked solution is just a fixed point of the vector equation $\sigma {\bf g}b({\bm \theta},{\bf g}ammab) = {\bm \omega}$. We will call a phase-locked solution stable if it is stable under all mean-zero perturbations, namely, perturbation orthogonal to ${\bf 1}$. This is because ${\bf g}b(\cdot,{\bf g}ammab)$ is inherently invariant under translations by the vector ${\bf 1}$. From our choice of sign, we conclude that a phase-locked solution ${\bm \theta}$ is stable if the eigenvalues $\lambda_1,\dots,\lambda_n$ of the Jacobian matrix ${\bf J}b({\bm \theta},{\bf g}ammab)$ are all positive with the exception of $\lambda_1 = 0$ corresponding to the eigenvector ${\bf 1}$. Finally, we note that ${\bf J}b({\bm \theta},{\bf g}ammab)$ is a Laplacian matrix, namely, a symmetric matrix with mean zero rows. Laplacian matrices form an important class of matrices and have an important connection to networks. One of the most famous connections between Laplacian matrices and networks is the Matrix Tree Theorem \cite{MR3097651}. We now state Theorem \ref{thm:detred} which is the special case of the famous Matrix Tree Theorem for ring networks.

{\bm \epsilon}gin{theorem}[Matrix Tree Theorem] \label{thm:detred}
The reduced determinant satisfies the identities
{\bm \epsilon}gin{align*}
\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) := \det({\bf P} {\bf J}b({\bm \theta},{\bf g}ammab) {\bf P}^\top) = n \sum_{i=1}^n \prod_{j \ne i} {\bf g}amma_j \cos \eta_j = n \left( \prod_{j=1}^n {\bf g}amma_j \cos \eta_j \right) \sum_{i=1}^n \frac{1}{{\bf g}amma_i \cos \eta_i} = \prod_{i=2}^n \lambda_i
\end{align*}
where $\lambda_1, \dots, \lambda_n$ are the eigenvalues of ${\bf J}b({\bm \theta},{\bf g}ammab)$ and $\lambda_1 = 0$ corresponds to the eigenvector ${\bf 1}$.
\end{theorem}

We require one more definition before stating our first theorem.

{\bm \epsilon}gin{definition} \label{def:vn}
Define the vector fields ${\bf v}(\cdot,{\bf g}ammab), {\bf n}(\cdot,{\bf g}ammab) : \mathbb{R}^n \rightarrow \mathbb{R}^n$ by
{\bm \epsilon}gin{align}
{\bf v}({\bm \theta},{\bf g}ammab) := -{\bf P}^\top \adj({\bf P} {\bf J}b({\bm \theta},{\bf g}ammab) {\bf P}^\top) {\bf P} {\bf g}b({\bm \theta},{\bf g}ammab) \quad \text{and} \quad {\bf n}({\bm \theta},{\bf g}ammab) := \nabla_{\bm \theta} \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))
\end{align}
where $\adj$ is the adjugate matrix. Furthermore, define their inner product
{\bm \epsilon}gin{align}
\delta({\bm \theta},{\bf g}ammab) := {\bf n}({\bm \theta},{\bf g}ammab)^\top {\bf v}({\bm \theta},{\bf g}ammab).
\end{align}
\end{definition}

Now we can state our first main result Theorem \ref{thm:criterion1} which is a criterion for the Kuramoto model \eqref{eq:model} to have a bifurcation in which two branches of phase-locked solutions collide as we increase $\sigma$. We defer the proof of Theorem \ref{thm:criterion1} to the appendix.

{\bm \epsilon}gin{theorem} \label{thm:criterion1}
Fix ${\bm \theta}_0$ and ${\bf g}ammab$ satisfying
{\bm \epsilon}gin{align} \label{eq:hypothesis}
\det_\text{Red}({\bf J}b({\bm \theta}_0,{\bf g}ammab)) = 0 \quad \text{and} \quad \delta({\bm \theta}_0,{\bf g}ammab) < 0.
\end{align}
Furthermore, fix $\sigma_0$ and set ${\bm \omega} = \sigma_0 {\bf g}b({\bm \theta}_0,{\bf g}ammab)$. Then the Kuramoto model \eqref{eq:model} with parameters ${\bm \omega}$ and ${\bf g}ammab$ has a bifurcation at $({\bm \theta},\sigma) = ({\bm \theta}_0,\sigma_0)$ which involves two branches of phase-locked solutions colliding as $\sigma$ increases to $\sigma_0$. Locally one branch satisfies $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) > 0$ and the other $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) <0$. In addition, if we suppose that the component-wise vector $\cos {\bm \eta}_0$ has no zero components, then locally on each branch
{\bm \epsilon}gin{align}
\text{$\#$positive eigenvalues of ${\bf J}b({\bm \theta},{\bf g}ammab)$} = n_+(\cos {\bm \eta}_0) -
{\bm \epsilon}gin{cases}
1 & \mbox{if $(-1)^{n-n_+(\cos {\bm \eta}_0)} \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) > 0$,}
\\
0 & \mbox{if $(-1)^{n-n_+(\cos {\bm \eta}_0)} \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) < 0$,}
\end{cases}
\end{align}
where $n_+(\cos {\bm \eta}_0)$ is the number of positive components of $\cos {\bm \eta}_0$. In particular, if $n_+(\cos {\bm \eta}_0) = n-1$, then a branch of stable phase-locked solutions collide with a branch of 1-saddles.
\end{theorem}

We note that the bifurcation in Theorem \ref{thm:criterion1} is a subcritical saddle-node bifurcation if there are no other branches in the bifurcation. This occurs if for each sufficiently small $\epsilon > 0$ there exists a neighborhood $({\bm \theta}_0,\sigma_0) \in U_\epsilon \subseteq \mathbb{R}^n$ for which
{\bm \epsilon}gin{align} \label{eq:local}
\{ ({\bm \theta},\sigma) \in U_\epsilon : \sigma {\bf g}b({\bm \theta},{\bf g}ammab) = {\bm \omega} \} = \{ ({\bm \theta}(s), \sigma(s)) \in \mathbb{R}^n : |s| < \epsilon \}.
\end{align}
Now since ${\bf v}({\bm \theta}_0,{\bf g}ammab) \ne {\bf 0}$ we see that \eqref{eq:local} holds if and only if $({\bm \theta}_0,\sigma_0)$ is an isolated zero of the analytic variety $V = \{ ({\bm \theta},\sigma) \in \mathbb{R}^{n+1} : \sigma {\bf g}b({\bm \theta},{\bf g}ammab) = {\bm \omega} \text{ and } \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0 \}$. Of course we mean that $({\bm \theta}_0,\sigma_0)$ is an isolated zero up to translation of ${\bm \theta}_0$ by multiples of ${\bf 1}$. This is true in all of our examples and we conjecture that it is true in generically. Heuristically, once we remove translation of ${\bm \theta}$ by ${\bf 1}$, $V$ is actually an analytic variety in $n$ variables. Furthermore, $\sigma {\bf g}b({\bm \theta},{\bf g}ammab) = {\bm \omega}$ contributes $n-1$ linearly independent equations since ${\bf 1}^\top {\bf g}b({\bm \theta},{\bf g}ammab) = {\bf 1}^\top {\bm \omega} = {\bf 0}$ for all values of ${\bm \theta}$ and ${\bf g}ammab$. Therefore $V$ represents an analytic variety in $n$ variables and $n$ equations. For an analysis of the Kuramoto model using algebraic geometry see \cite{2015Chaos..25e3103M}.

Now Theorem \ref{thm:criterion1} is a consequence of Lemma \ref{lem:flow}.

{\bm \epsilon}gin{lemma} \label{lem:flow}
Fix ${\bm \theta}_0$, ${\bf g}ammab$, and $\sigma_0$ and set ${\bm \omega} = \sigma_0 {\bf g}b({\bm \theta}_0,{\bf g}ammab)$. Then there exists a unique analytic vector-valued function ${\bm \theta}(\cdot) = {\bm \theta}(\cdot, {\bm \theta}_0,{\bf g}ammab) : \mathbb{R} \rightarrow \mathbb{R}^n$ satisfying
{\bm \epsilon}gin{align} \label{eq:flow}
\frac{d {\bm \theta}}{ds} = {\bf v}({\bm \theta}(s),{\bf g}ammab)
\end{align}
with ${\bm \theta}(0) = {\bm \theta}_0$. Now define the analytic scalar-valued function $\sigma(\cdot) = \sigma(\cdot,{\bm \theta}_0,{\bf g}ammab, \sigma_0) : \mathbb{R} \rightarrow \mathbb{R}_{>0}$ by
{\bm \epsilon}gin{align} \label{eq:sigma}
\sigma(s) = \sigma_0 \exp \left(\int_0^s \det_\text{Red}({\bf J}b({\bm \theta}({\tilde{s}}),{\bf g}ammab)) d{\tilde{s}} \right).
\end{align}
Then we have the following properties:
{\bm \epsilon}gin{enumerate}
\item ${\bm \theta}(s)$ is a phase-locked solution of the Kuramoto model \eqref{eq:model} with fixed parameters ${\bm \omega}$ and ${\bf g}ammab$ and coupling $\sigma(s)$ i.e.
{\bm \epsilon}gin{align} \label{eq:fixed}
\sigma(s) {\bf g}b({\bm \theta}(s),{\bf g}ammab) = {\bm \omega}
\end{align}
for all $s \in \mathbb{R}$.
\item $\frac{d \sigma}{ds}$ has the same sign as $\det_\text{Red}({\bf J}b({\bm \theta}(s),{\bf g}ammab))$.
\end{enumerate}
\end{lemma}

Essentially, the conditions \eqref{eq:hypothesis} in Theorem \ref{thm:criterion1} cause $\sigma(s)$ in Lemma \ref{lem:flow} to have a local maximum at $\sigma(0) = \sigma_0$. We defer the details to the appendix, but give an example illustrating Theorem \ref{thm:criterion1} and Lemma \ref{lem:flow}. We start our example by setting
{\bm \epsilon}gin{align} \label{eq:3}
{\bm \theta}_0 = (0.86727,1.84593,3.88114) \quad \text{and} \quad {\bf g}ammab = (0.55284,1,1).
\end{align}
These values were obtained by using Mathematica to minimize $\delta({\bm \theta},{\bf g}ammab)$ subject to the conditions $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) \\ = 0$ and ${\bf g}ammab \in [0,1]^3$. In fact we have the value $\delta({\bm \theta}_0,{\bf g}ammab) = -0.37347$. In Figures \ref{fig:2a} and \ref{fig:2b} we plot the functions ${\bm \theta}(s)$ and $\sigma(s)$ from Lemma \ref{lem:flow}. In Figure \ref{fig:2c} we project the level set $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0$, the normal vector ${\bf n}({\bm \theta},{\bf g}ammab)$, and the orbit of ${\bm \theta}(s)$ onto the mean-zero plane ${\bf 1}^\perp$. Note that the orbit of ${\bm \theta}(s)$ just barely crosses the level set because it makes an obtuse angle (greater than a right angle) with the normal vector. This crossing results in the bifurcation and can also been seen in the local maximum of $\sigma(s)$ at $s = 0$.

{\bm \epsilon}gin{figure}[H]
\captionsetup{font=scriptsize}
\centering
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plott3.pdf}
  \caption{} \label{fig:2a}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plotsn3.pdf}
  \caption{} \label{fig:2b}
\end{subfigure}

{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plotdvt3.pdf}
  \caption{} \label{fig:2c}
\end{subfigure}
\caption{(A) A plot of ${\bm \theta}(s)$. (B) A plot of $\sigma(s)$. The dashed red lines indicate the number of positive eigenvalues of ${\bf J}b({\bm \theta}(s),{\bf g}ammab)$. The dashed green line is simply to help see that $\sigma(s)$ has a local maximum at $s = 0$.(C) We projected $\mathbb{R}^3$ onto the mean-zero plane ${\bf 1}^\perp$ via the change of variables ${\bm \phi} = {\bf P} {\bm \theta}$. The blue ``circles" are connected components of the level set $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0$, the green arrow is the normal vector ${\bf n}({\bm \theta},{\bf g}ammab)$, and the red curve is the orbit of ${\bm \theta}(s)$.} \label{fig:2}
\end{figure}

Now we give another example. This time however our bifurcation will have a branch consisting of stable phase-locked solutions. In fact, this is the bifurcation shown in Figure \ref{fig:1b}. Again we start by setting
{\bm \epsilon}gin{align} \label{eq:5}
{\bm \theta}_0 = (3.80063,1.47760,6.77075,5,78071,4.79067) \quad \text{and} \quad {\bf g}ammab = (0.20075,1,1,1,1).
\end{align}
Again, these values were obtained by using Mathematica to minimize $\delta({\bm \theta},{\bf g}ammab)$ but now subject to the conditions $a_1({\bm \theta},{\bf g}ammab) = 0$, $a_i({\bm \theta},{\bf g}ammab) {\bf g}e 0$ for $2 \le i \le 4$, and ${\bf g}ammab \in [0,1]^5$ where $\det(\lambda {\bf I}b - {\bf J}b({\bm \theta},{\bf g}ammab)) = \lambda^5 + \sum_{i=1}^4 a_i({\bm \theta},{\bf g}ammab) (-\lambda)^i$. In fact we have the value $\delta({\bm \theta}_0,{\bf g}ammab) = -0.70959$. Note that $a_0({\bm \theta},{\bf g}ammab) = \det({\bf J}b({\bm \theta},{\bf g}ammab)) = 0$ identically and that $a_1({\bm \theta},{\bf g}ammab) = 0$ is equivalent to $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0$. We also note that these conditions on the characteristic polynomial imply that ${\bf J}b({\bm \theta},{\bf g}ammab)$ has two zero eigenvalues (one trivial and one non-trivial) and that the rest are positive. In Figures \ref{fig:3a} and \ref{fig:3b} we again plot ${\bm \theta}(s)$ and $\sigma(s)$ from Lemma \ref{lem:flow}. Again we see that $\sigma(s)$ has a local maximum at $s = 0$ and that on one branch ${\bf J}b({\bm \theta},{\bf g}ammab)$ has four positive eigenvalues corresponding to stable phase-locked solutions. 

{\bm \epsilon}gin{figure}[H]
\captionsetup{font=scriptsize}
\centering
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plott5.pdf}
  \caption{} \label{fig:3a}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plotsn5.pdf}
  \caption{} \label{fig:3b}
\end{subfigure}
\caption{(A) A plot of ${\bm \theta}(s)$. (B) A plot of $\sigma(s)$. The dashed red lines indicate the number of positive eigenvalues of ${\bf J}b({\bm \theta}(s),{\bf g}ammab)$. The dashed green line is simply to help see that $\sigma(s)$ has a local maximum at $s = 0$.} \label{fig:3}
\end{figure}

For our next main result we derive a condition on ${\bm \theta}_0$ that guarantees the existence of a ${\bf g}ammab$ satisfying the conditions \eqref{eq:hypothesis} in Theorem \ref{thm:criterion1}. Before we can do this, however, we require some definitions and supporting lemmas. Again we defer the proofs to the appendix.

{\bm \epsilon}gin{definition} \label{def:h}
Define the vector-valued function ${\bf h} : \mathbb{R}^n \rightarrow \mathbb{R}^n$ by
{\bm \epsilon}gin{align} \label{eq:h}
{\bf h}({\bf x})_i = \prod_{j \ne i} x_j.
\end{align}
\end{definition}

{\bm \epsilon}gin{definition} \label{def:abc}
Define the three $n \times n$ matrices ${\bf A}({\bm \theta})$, ${\bf B}$, and ${\bf C}({\bm \theta})$ by
{\bm \epsilon}gin{align*}
{\bf A}({\bm \theta})_{ij} := n \frac{\partial}{\partial \theta_j} {\bf h}(\cos {\bm \eta})_i
\quad
{\bf B}_{ij} :=
{\bm \epsilon}gin{cases}
1 & \mbox{if $i = j+1$,}
\\
-1 & \mbox{if $i = j$,}
\\
0 & \mbox{otherwise,}
\end{cases}
\quad
{\bf C}({\bm \theta})_{ij} =
{\bm \epsilon}gin{cases}
n \left( \sum_{k \ne i} \frac{\sin \eta_k}{\cos \eta_k} \right) {\bf h}(\cos {\bm \eta})_j & \mbox{if $i = j$,}
\\
-n \frac{\sin \eta_i}{\cos \eta_i} {\bf h}(\cos {\bm \eta})_j & \mbox{if $i \ne j$.}
\end{cases}
\end{align*}
Note that ${\bf B}$ is the incidence matrix for a ring network with $n$ vertices and $n$ oriented edges $e_i = (i,i+1)$. Furthermore, define the two $n \times n$ matrices ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$ by
{\bm \epsilon}gin{align}
{\bf S}({\bm \theta}) = n^2 {\bf h}(\cos {\bm \eta})^\top {\bf h}(\cos {\bm \eta}) \quad \text{and} \quad {\bf T}({\bm \theta}) = {\bf A}({\bm \theta}) ({\bf B}^\top)^{-1} {\bf C}({\bm \theta})
\end{align}
where $({\bf B}^\top)^{-1}$ is the pseudo-inverse of ${\bf B}^\top$.
\end{definition}

Continuing our example for $n = 3$ we have that
{\bm \epsilon}gin{gather*}
{\bf A}({\bm \theta}) =
{\bm \epsilon}gin{pmatrix}
-\cos \eta _ 2 \sin \eta _ 3 & \cos \eta _ 3 \sin \eta _ 2 & \cos \eta _ 2 \sin \eta _ 3-\cos \eta _ 3 \sin \eta _ 2
\\
\cos \eta _ 3 \sin \eta _ 1-\cos \eta _ 1 \sin \eta _ 3 & -\cos \eta _ 3 \sin \eta _ 1 & \cos \eta _ 1 \sin \eta _ 3
\\
 \cos \eta _ 2 \sin \eta _ 1 & \cos \eta _ 1 \sin \eta _ 2-\cos \eta _ 2 \sin \eta _ 1 & -\cos \eta _ 1 \sin \eta _ 2
\end{pmatrix},
\\
{\bf B} =
{\bm \epsilon}gin{pmatrix}
-1 & 0 & 1
\\
1 & -1 & 0
\\
0 & 1 & -1
\end{pmatrix},
\\
{\bf C}({\bm \theta}) =
{\bm \epsilon}gin{pmatrix}
3 \cos \eta _ 3 \sin \eta _ 2+3 \cos \eta _ 2 \sin \eta _ 3 & -3 \cos \eta _ 3 \sin \eta _ 1 & -3 \cos \eta _ 2 \sin \eta _ 1
\\
-3 \cos \eta _ 3 \sin \eta _ 2 & 3 \cos \eta _ 3 \sin \eta _ 1+3 \cos \eta _ 1 \sin \eta _ 3 & -3 \cos \eta _ 1 \sin \eta _ 2
\\
-3 \cos \eta _ 2 \sin \eta _ 3 & -3 \cos \eta _ 1 \sin \eta _ 3 & 3 \cos \eta _ 2 \sin \eta _ 1+3 \cos \eta _ 1 \sin \eta _ 2
\end{pmatrix}.
\end{gather*}

The motivation for Definitions \ref{def:h} and \ref{def:abc} are the quadratic form identities in Lemma \ref{lem:quadratic}.

{\bm \epsilon}gin{lemma} \label{lem:quadratic}
We have the quadratic form identities
{\bm \epsilon}gin{align}
\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))^2 = {\bf h}({\bf g}ammab)^\top {\bf S}({\bm \theta}) {\bf h}({\bf g}ammab) \quad \text{and} \quad \delta({\bm \theta},{\bf g}ammab) = {\bf h}({\bf g}ammab)^\top {\bf T}({\bm \theta}) {\bf h}({\bf g}ammab).
\end{align}
Furthermore, the restriction ${\bf h} : \mathbb{R}_{>0}^n \rightarrow \mathbb{R}_{>0}^n$ is a bijection with inverse
{\bm \epsilon}gin{align} \label{eq:hinverse}
({\bf h}({\bf x}b)^{-1})_i = \frac{1}{x_i} \left( \prod_{j=1}^n x_j \right)^{\frac{1}{n-1}}.
\end{align}
\end{lemma}

Next, we combine our matrices ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$ in Definition \ref{def:r}.

{\bm \epsilon}gin{definition} \label{def:r}
Define the parameterized family of symmetric $n \times n$ matrices ${\bf R}_\tau({\bm \theta})$ by
{\bm \epsilon}gin{align}
{\bf R}_\tau({\bm \theta}) = {\bf S}({\bm \theta}) \tau + \frac{1}{2}({\bf T}({\bm \theta}) + {\bf T}({\bm \theta})^\top).
\end{align}
\end{definition}

The motivation for Definition \ref{def:r} is the identity
{\bm \epsilon}gin{align} \label{eq:quadratic}
{\bf h}({\bf g}ammab)^\top {\bf R}_\tau({\bm \theta}) {\bf h}({\bf g}ammab) = \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))^2 \tau + \delta({\bm \theta},{\bf g}ammab)
\end{align}
which easily follows from Lemma \ref{lem:quadratic}. We now make the observation that if ${\bm \theta}_0$ and ${\bf g}ammab$ satisfy \eqref{eq:hypothesis}, then ${\bf x}b^\top {\bf R}_\tau({\bm \theta}_0) {\bf x}b < 0$ for all $\tau$ where ${\bf x}b = {\bf h}({\bf g}ammab) \in \mathbb{R}_{>0}^n$ by Lemma \ref{lem:quadratic}. Conversely, we will show that if there exists vectors ${\bf x}b_\tau \in \mathbb{R}_{{\bf g}e 0}^n$ such that ${\bf x}b_\tau^\top {\bf R}_\tau({\bm \theta}_0) {\bf x}b_\tau / \| {\bf x}b_\tau \|_1^2 < \epsilon < 0$ for all sufficiently large $\tau$, then there exists a vector ${\bf x}b \in \mathbb{R}_{>0}^n$ such that ${\bf x}b^\top {\bf R}_\tau(\theta_0) {\bf x}b < 0$. By Lemma 4 we can choose ${\bf g}ammab$ such that ${\bf x}b = {\bf h}({\bf g}ammab)$ in which case \eqref{eq:hypothesis} is satisfied. In order to demonstrate the existence of our vectors ${\bf x}b_\tau \in \mathbb{R}_{{\bf g}e 0}^n$ we require Lemma \ref{lem:min} which is a slight modification of a result of Kreps in \cite{KREPS1984105}. Before stating Lemma \ref{lem:min} we first make Definition \ref{def:min}.

{\bm \epsilon}gin{definition} \label{def:min}
Let ${\bf R}$ be an $n \times n$ matrix and ${\bf I}c \subseteq \{1, \dots, n\}$. Define ${\bf R}_{\bf I}c$ to be the sub-matrix of ${\bf R}$ with rows and columns corresponding to the indices in ${\bf I}c$. In other words, if ${\bf I}c = \{ i_1 < \dots < i_m\}$, then ${\bf R}_{\bf I}c$ is the $m \times m$ matrix with entries $({\bf R}_{\bf I}c)_{k\ell} = {\bf R}_{i_k i_\ell}$. Furthermore, given $i \in {\bf I}c$, define ${\bf R}_{{\bf I}c,i}$ to be the matrix ${\bf R}_{\bf I}c$ with the $i$th row replaced by a row of all ones.
\end{definition}

{\bm \epsilon}gin{lemma} \label{lem:min}
For any symmetric $n \times n$ matrix ${\bf R}$ we have the identity
{\bm \epsilon}gin{align} \label{eq:min}
\min_{{\bf x}b \in \mathbb{R}_{{\bf g}e 0}^n} \frac{{\bf x}b^\top {\bf R} {\bf x}b}{\| {\bf x}b \|_1^2} = \min \left\{ \frac{\det({\bf R}_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i})} : {\bf I}c \subseteq \{1,\dots,n\} \text{ satisfies } \min_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i}) > 0 \right\}
\end{align}
where $\| {\bf x}b \|_1 := \sum_{i=1}^n |x_i|$.
\end{lemma}

We note that $\det({\bf R}_\tau({\bm \theta})_{\bf I}c)$ and $\det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i})$ are polynomials in $\tau$. In fact, they have degree at most one since ${\bf S}({\bm \theta})$ has rank one. Now in order to apply Lemma \ref{lem:min} we make Definition \ref{def:limit}.

{\bm \epsilon}gin{definition} \label{def:limit}
Define
{\bm \epsilon}gin{align*}
{\bf I}c({\bm \theta}) := \{ {\bf I}c \subseteq \{1,\dots,n\} : \min_{i \in {\bf I}c} \det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i}) > 0 \text{ for all sufficiently large $\tau$} \}
\end{align*}
and for each ${\bf I}c \in {\bf I}c({\bm \theta})$ define
{\bm \epsilon}gin{align*}
p_{\bf I}c({\bm \theta}) :=
{\bm \epsilon}gin{cases}
1 & \mbox{if $\lim_{\tau \rightarrow \infty} \frac{\det({\bf R}_\tau({\bm \theta})_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i})} > 0$,}
\\
0 & \mbox{if $\lim_{\tau \rightarrow \infty} \frac{\det({\bf R}_\tau({\bm \theta})_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i})} = 0$,}
\\
-1 & \mbox{if $\lim_{\tau \rightarrow \infty} \frac{\det({\bf R}_\tau({\bm \theta})_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i})} < 0$,}
\end{cases}
\end{align*}
and
{\bm \epsilon}gin{align*}
q_{\bf I}c({\bm \theta}) :=
{\bm \epsilon}gin{cases}
1 & \mbox{if $\det({\bf R}_\tau({\bm \theta})_{\bf I}c) > 0$ for all sufficiently large $\tau$,}
\\
0 & \mbox{if $\det({\bf R}_\tau({\bm \theta})_{\bf I}c) = 0$ for all sufficiently large $\tau$,}
\\
-1 & \mbox{if $\det({\bf R}_\tau({\bm \theta})_{\bf I}c) < 0$ for all sufficiently large $\tau$,}
\end{cases}
\end{align*}
and set
{\bm \epsilon}gin{align*}
{\bf D}elta_p({\bm \theta}) := \min \{ p_{\bf I}c({\bm \theta}) : {\bf I}c \in {\bf I}c({\bm \theta}) \} \quad \text{and} \quad {\bf D}elta_q({\bm \theta}) := \min \{ q_{\bf I}c({\bm \theta}) : {\bf I}c \in {\bf I}c({\bm \theta}) \}.
\end{align*}
\end{definition}

The motivation for the definitions of ${\bf D}elta_p({\bm \theta})$ and ${\bf D}elta_q({\bm \theta})$ is that they in some sense model the sign of \eqref{eq:min} for ${\bf R}_\tau({\bm \theta})$ for large values of $\tau$. We plot ${\bf D}elta_p({\bm \theta})$ and ${\bf D}elta_q({\bm \theta})$ in Figure \ref{fig:4} for $n = 3$. Now we can state our next main result Theorem \ref{thm:criterion2}.

{\bm \epsilon}gin{theorem} \label{thm:criterion2}
There exists a ${\bf g}ammab$ satisfying \eqref{eq:hypothesis} if ${\bf D}elta_p({\bm \theta}_0) = -1$. Conversely, \eqref{eq:hypothesis} is not satisfied for any ${\bf g}ammab$ if ${\bf D}elta_q({\bm \theta}_0) \in \{0,1\}$.
\end{theorem}

We make a brief comment on the strength of Theorem \ref{thm:criterion2}. Now
{\bm \epsilon}gin{align} \label{eq:exhaustive}
 \{ {\bm \theta} \in \mathbb{R}^n : {\bf D}elta_p({\bm \theta}) \ne {\bf D}elta_q({\bm \theta}) \} \subseteq \bigcup_{{\bf I}c \subseteq \{1, \dots,n\}} \{ {\bm \theta} \in \mathbb{R}^n : \ell_{\bf I}c({\bm \theta}) = 0\}
\end{align}
where $\ell_{\bf I}c({\bm \theta}) = \lim_{\tau \rightarrow \infty} \frac{1}{\tau} \det({\bf R}_\tau({\bm \theta})_{\bf I}c)$. This is because $\det({\bf R}_\tau({\bm \theta})_{\bf I}c)$ and $\det({\bf R}_\tau({\bm \theta})_{{\bf I}c,i})$ are polynomials of degree at most one in $\tau$ and therefore $p_{\bf I}c({\bm \theta}) = q_{\bf I}c({\bm \theta})$ if $\ell_{\bf I}c({\bm \theta}) \ne 0$. By properties of the determinant we know that $\ell_{\bf I}c({\bm \theta})$ is a sum of determinants of matrices with entries from ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$. Since the entries of ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$ are analytic we easily deduce that $\ell_{\bf I}c({\bm \theta})$ is also analytic. Therefore both sets in \eqref{eq:exhaustive} have $n$-dimensional Lebesgue measure zero if $\ell_{\bf I}c({\bm \theta})$ is not identically zero for all ${\bf I}c \subseteq \{1,\dots,n\}$. In particular, this would imply that the exceptional set of Theorem \ref{thm:criterion2} has $n$-dimensional Lebesgue measure zero. Unfortunately, we have yet to discover a proof that $\ell_{\bf I}c({\bm \theta})$ is not identically zero due to the complexity of the definition of $\ell_{\bf I}c({\bm \theta})$. However, we conjecture that this is the case and give some examples for $n = 3$. In particular,
{\bm \epsilon}gin{gather*}
\ell_{\{1\}}({\bm \theta}) =
\cos ^2\eta _2 \cos ^2\eta _3,
\\
\ell_{\{1,2\}}({\bm \theta}) =
-6 \sin \eta _1 \sin \eta _2 \cos \eta _1 \cos \eta _2 \cos ^4\eta _3-3 \sin ^2\eta _1 \cos ^2\eta _2 \cos ^4\eta _3 -3 \sin ^2\eta _2 \cos ^2\eta _1 \cos ^4\eta _3-
\\
3 \sin \eta _1 \sin \eta _3 \cos \eta _1 \cos ^2\eta _2 \cos ^3\eta _3-3 \sin \eta _2 \sin \eta _3 \cos ^2\eta _1 \cos \eta _2 \cos ^3\eta _3,
\\
\ell_{\{1,2,3\}}({\bm \theta}) =
9 \sin \eta _2 \sin ^3\eta _3 \cos ^3\eta _2 \cos \eta _3 \cos ^4\eta _1+9 \sin ^3\eta _2 \sin \eta _3 \cos \eta _2 \cos ^3\eta _3 \cos ^4\eta _1+
\\
9 \sin \eta _1 \sin ^3\eta _2 \cos \eta _2 \cos ^4\eta _3 \cos ^3\eta _1+9 \sin \eta _1 \sin ^3\eta _3 \cos ^4\eta _2 \cos \eta _3 \cos ^3\eta _1+
\\
9 \sin ^3\eta _1 \sin \eta _2 \cos ^3\eta _2 \cos ^4\eta _3 \cos \eta _1+9 \sin ^3\eta _1 \sin \eta _3 \cos ^4\eta _2 \cos ^3\eta _3 \cos \eta _1+
\\
18 \sin ^2\eta _2 \sin ^2\eta _3 \cos ^2\eta _2 \cos ^2\eta _3 \cos ^4\eta _1+18 \sin ^2\eta _1 \sin ^2\eta _2 \cos ^2\eta _2 \cos ^4\eta _3 \cos ^2\eta _1+
\\
18 \sin ^2\eta _1 \sin ^2\eta _3 \cos ^4\eta _2 \cos^2\eta _3 \cos ^2\eta _1+45 \sin \eta _1 \sin \eta _2 \sin ^2\eta _3 \cos ^3\eta _2 \cos ^2\eta _3 \cos ^3\eta _1+
\\
45 \sin \eta _1 \sin ^2\eta _2 \sin \eta _3 \cos ^2\eta _2 \cos ^3\eta _3 \cos ^3\eta _1+ 45 \sin ^2\eta _1 \sin \eta _2 \sin \eta _3 \cos ^3\eta _2 \cos ^3\eta _3 \cos ^2\eta _1,
\end{gather*}
and
{\bm \epsilon}gin{align*}
\ell_{\{1\}}({\bm \theta}_0) = \frac{1}{16} \ne 0, \quad \ell_{\{1,2\}}({\bm \theta}_0) = -\frac{27}{128} \ne 0, \quad \ell_{\{1,2,3\}}({\bm \theta}_0) = \frac{2187}{4096} \ne 0,
\end{align*}
where ${\bm \theta}_0 = (0,\pi/3,2\pi/3)$. Of course the remaining functions $\ell_{\{2\}}({\bm \theta})$, $\ell_{\{3\}}({\bm \theta})$, $\ell_{\{1,3\}}({\bm \theta})$, and $\ell_{\{2,3\}}({\bm \theta})$ cannot be identically zero by symmetry.

Finally, we show that the two sets
{\bm \epsilon}gin{align} \label{eq:set1}
\{ {\bm \theta} \in \mathbb{R}^n : {\bf D}elta_p({\bm \theta}) = -1 \}
\end{align}
and
{\bm \epsilon}gin{align} \label{eq:set2}
\{ {\bm \theta} \in \mathbb{R}^n : \text{${\bf D}elta_p({\bm \theta}) = -1$, $\cos {\bm \eta}$ has no zero components, and $n_+(\cos {\bm \eta}) = n-1$} \}
\end{align}
are non-empty for all $n {\bf g}e 3$ and $n {\bf g}e 5$ respectively. Therefore by applying Theorems \ref{thm:criterion1} and \ref{thm:criterion2} we deduce our final main result Theorem \ref{thm:existence}.

{\bm \epsilon}gin{theorem} \label{thm:existence}
For every $n {\bf g}e 3$, there exists an ${\bm \omega}$ and ${\bf g}ammab$ for which the Kuramoto model \eqref{eq:model} has a bifurcation which involves two branches of phase-locked solutions colliding as $\sigma$ increases. Furthermore, for every $n {\bf g}e 5$, there exists an ${\bm \omega}$ and ${\bf g}ammab$ for which the Kuramoto model \eqref{eq:model} has a bifurcation which involves a branch of stable phase-locked solutions colliding with a branch of 1-saddles as $\sigma$ increases.
\end{theorem}

We include the proof of Theorem \ref{thm:existence} here because the proof is an application of an important property of the matrices ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$ given in Lemma \ref{lem:st}.

{\bm \epsilon}gin{lemma} \label{lem:st}
Let $1 \le n \le {\tilde{n}}$ and for any ${\bm \theta} = (\theta_1, \dots, \theta_n) \in \mathbb{R}^n$ set ${\bm \theta}t = (\theta_1, \dots, \theta_n, \theta_1, \dots, \theta_1) \in \mathbb{R}^{{\tilde{n}}}$. Then
{\bm \epsilon}gin{align*}
{\bf S}t({\bm \theta}t) =
{\bm \epsilon}gin{pmatrix}
(\frac{{\tilde{n}}}{n})^2 {\bf S}({\bm \theta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
\quad \text{and} \quad
{\bf T}t({\bm \theta}t) =
{\bm \epsilon}gin{pmatrix}
(\frac{{\tilde{n}}}{n})^2 {\bf T}({\bm \theta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}.
\end{align*}
In other words, ${\bf S}({\bm \theta})$ and ${\bf T}({\bm \theta})$ are positive multiples of the upper left $n \times n$ submatrices of ${\bf S}t({\bm \theta}t)$ and ${\bf T}t({\bm \theta}t)$ respectively.
\end{lemma}

{\bm \epsilon}gin{proof}[Proof of Theorem \ref{thm:existence}]
We start with the first statement for $n {\bf g}e 3$ and set
{\bm \epsilon}gin{align} \label{eq:nge3}
{\bm \theta}t_0 = (0.86727,1.84593,3.88114,0.86727,\dots,0.86727) \in \mathbb{R}^n.
\end{align}
Note that ${\bm \theta}t_0$ reduces to ${\bm \theta}_0$ in \eqref{eq:3} when $n = 3$. By Lemma \ref{lem:st} we know that
{\bm \epsilon}gin{align*}
\det({\bf R}t_\tau({\bm \theta}t_0)_{\{1,2\}}) = \left(\frac{n}{3} \right)^4 \det({\bf R}_\tau({\bm \theta}_0)_{\{1,2\}}) = \left(\frac{n}{3} \right)^4 \left(-0.01174 n^4 \tau -0.01162 \right)
\end{align*}
and
{\bm \epsilon}gin{align*}
\det({\bf R}t_\tau({\bm \theta}t_0)_{\{1,2\},i}) = \left(\frac{n}{3} \right)^2 \det({\bf R}_\tau({\bm \theta}_0)_{\{1,2\},i}) = \left(\frac{n}{3} \right)^2
{\bm \epsilon}gin{cases}
0.55240 \tau -1.42429 & \mbox{if $i = 1$,}
\\
0.44330 \tau -1.53160 & \mbox{if $i = 2$.}
\end{cases}
\end{align*}
Therefore ${\bf D}elta_p({\bm \theta}t_0) = -1$ for each $n {\bf g}e 3$ and the result follows from Theorem \ref{thm:criterion2}.

In a similar way we prove the second statement for $n {\bf g}e 5$ by setting
{\bm \epsilon}gin{align} \label{eq:nge5}
{\bm \theta}t_0 = (3.80063,1.47760,6.77075,5,78071,4.79067,3.80063,\dots,3.80063) \in \mathbb{R}^n.
\end{align}
Again note that ${\bm \theta}t_0$ reduces to ${\bm \theta}_0$ in \eqref{eq:5} when $n = 5$. By Lemma \ref{lem:st} we know that
{\bm \epsilon}gin{align*}
\det({\bf R}t_\tau({\bm \theta}t_0)_{\{1,2\}}) = \left(\frac{n}{5} \right)^4 \det({\bf R}_\tau({\bm \theta}_0)_{\{1,2\}}) = \left(\frac{n}{5} \right)^4\left( -0.00024 \tau + 0.00191 \right)
\end{align*}
and
{\bm \epsilon}gin{align*}
\det({\bf R}t_\tau({\bm \theta}t_0)_{\{1,2\},i}) = \left(\frac{n}{5} \right)^2 \det({\bf R}_\tau({\bm \theta}_0)_{\{1,2\},i}) = \left(\frac{n}{5} \right)^2
{\bm \epsilon}gin{cases}
0.02296 \tau -0.18050 & \mbox{if $i = 1$,}
\\
0.01844 \tau -0.15351 & \mbox{if $i = 2$.}
\end{cases}
\end{align*}
Therefore ${\bf D}elta_p({\bm \theta}t_0) = -1$ for each $n {\bf g}e 5$ and $\cos {\bm \eta}t_0 = (-0.68327,0.54866,0.54866,0.54866,0.54866,1,\dots,1)$. Therefore the result follows from Theorems \ref{thm:criterion1} and \ref{thm:criterion2}.
\end{proof}

For $3 \le n \le 4$ we used Mathematica to minimize $\delta({\bm \theta},{\bf g}ammab)$ subject to the same conditions on the characteristic polynomial as those used to find ${\bm \theta}_0$ and ${\bf g}ammab$ in \eqref{eq:5}. However, we were unable to obtain a negative value for $\delta({\bm \theta},{\bf g}ammab)$. This gives numerical evidence that Theorem \ref{thm:existence} is best possible although we admit that we do not have a proof of this.

We end this section by briefly presenting Theorem \ref{thm:existence} as a constructive algorithm. For the first part of Theorem \ref{thm:existence} choose any $n {\bf g}e 3$ and complete the following steps:
{\bm \epsilon}gin{enumerate}
\item Set ${\bm \theta}_0$ as in \eqref{eq:nge3},
\item Find a vector ${\bf x}b \in \mathbb{R}_{>0}^n$ such that ${\bf h}(\cos {\bm \eta}_0) {\bf x}b = 0$ and ${\bf x}b^\top {\bf T}({\bm \theta}_0) {\bf x}b < 0$,
\item Set ${\bf g}ammab = {\bf h}({\bf x}b)^{-1}$, $\sigma_0 = 1$, and ${\bm \omega} = {\bf g}b({\bm \theta}_0,{\bf g}ammab)$.
\end{enumerate}
The existence of the vector in (2) is guaranteed by the proofs of Theorems \ref{thm:criterion2} and \ref{thm:existence}. Then our choice of ${\bf g}ammab$ results in the equations $\det_\text{Red}({\bf J}b({\bm \theta}_0,{\bf g}ammab)) = {\bf h}(\cos {\bm \eta}_0) {\bf x}b = 0$ and $\delta({\bm \theta}_0,{\bf g}ammab) = {\bf x}b^\top {\bf T}({\bm \theta}_0) {\bf x}b < 0$. Therefore the Kuramoto model with parameters ${\bm \omega}$ and ${\bf g}ammab$ has the desired bifurcation at $({\bm \theta},\sigma) = ({\bm \theta}_0,\sigma_0)$ by Theorem \ref{thm:criterion1}.

For the second part of Theorem \ref{thm:existence} choose any $n {\bf g}e 5$ and set ${\bm \theta}_0$ as in \eqref{eq:nge5}.

\section{Other Directions and Future Work} \label{sec:future}

As mentioned in the introduction our bifurcation in Figure \ref{fig:1a} is locally a subcritical saddle-node bifurcation and globally an $S$-curve. We conjecture that this is generically the case. (Note that we discussed the genericity of the subcritical saddle-node bifurcation heuristically after the statement of Theorem \ref{thm:criterion1}.

One may wonder if we can numerically generate higher order $S$-curves. This can occur if we have a hysteresis point which we then perturb. In our case a hysteresis point is a pair $({\bm \theta}_0,{\bm \gamma})$ for which the curve ${\bm \theta}(s)$ is tangential to the level set $\det_\text{Red}({\bf J}b({\bm \theta}_0,{\bm \gamma})) = 0$. In other words, the curve has a higher order of contact with the level set. If we can make this order of contact high enough we should be able to construct a higher order $S$-curve. To do this we set
{\bm \epsilon}gin{align*}
\det_\text{Red}({\bf J}b({\bm \theta}(s),{\bm \gamma}) = \sum_{k=0}^\infty a_k({\bm \theta}_0,{\bm \gamma}) s^k
\end{align*}
and deduce that
{\bm \epsilon}gin{align*}
a_0({\bm \theta},{\bm \gamma}) = \det_\text{Red}({\bf J}b({\bm \theta},{\bm \gamma})) \quad \text{and} \quad  a_{k+1}({\bm \theta},{\bm \gamma}) = \langle {\bf v}({\bm \theta},{\bm \gamma}), \nabla_{\bm \theta} a_k({\bm \theta},{\bm \gamma}) \rangle
\end{align*}
for $k {\bf g}e 0$. We say that the curve ${\bm \theta}(s)$ has order of contact $d$ with the level set $\det_\text{Red}({\bf J}b({\bm \theta},{\bm \gamma})) = 0$ at ${\bm \theta}_0$ if $a_d({\bm \theta}_0,{\bm \gamma}) \ne 0$ and $a_k({\bm \theta}_0,{\bm \gamma}) = 0$ for $0 \le k < d$.

We executed this strategy and set $a_0({\bm \theta}_0,{\bm \gamma}) = a_1({\bm \theta}_0,{\bm \gamma}) = a_2({\bm \theta}_0,{\bm \gamma}) = 0$ and numerically found the values
{\bm \epsilon}gin{align*}
{\bm \theta}_0 = (5.00250, 3.78647, 2.09311, 4.35320, 5.58644, 0.99502),
\\
{\bm \gamma} = (0.87815, 1.22758, 2.45939, 0.51472, 6.30989, 2.91474).
\end{align*}
We then perturbed the value of ${\bm \theta}_0$ by
{\bm \epsilon}gin{align*}
{\bm \epsilon}_1 = (0.04000, 0.02000, 0.01000, 0.01000, 0.01000, -0.02968),
\\
{\bm \epsilon}_2 = (0.02000, -0.02000, 0.04000, -0.03000, -0.01000, -0.16604),
\end{align*}
which leads to the values
{\bm \epsilon}gin{align*}
\det_\text{Red}({\bf J}b({\bm \theta}_0,{\bf g}ammab)) &= 0 \quad \text{and} \quad \delta({\bm \theta}_0,{\bf g}ammab) = 0,
\\
\det_\text{Red}({\bf J}b({\bm \theta}_0 + {\bm \epsilon}_1,{\bf g}ammab)) &= 0 \quad \text{and} \quad \delta({\bm \theta}_0 + {\bm \epsilon}_1,{\bf g}ammab) = -6.18641,
\\
\det_\text{Red}({\bf J}b({\bm \theta}_0 + {\bm \epsilon}_2,{\bf g}ammab)) &= 0 \quad \text{and} \quad \delta({\bm \theta}_0 + {\bm \epsilon}_2,{\bf g}ammab) = 12.7292,
\end{align*}
and the three figures in Figure \ref{fig:4}.

{\bm \epsilon}gin{figure}[H]
\captionsetup{font=scriptsize}
\centering
{\bm \epsilon}gin{subfigure}{5cm}
  \centering
  \includegraphics[width=4.5cm]{Plots4a.pdf}
  \caption{} \label{fig:4a}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{5cm}
  \centering
  \includegraphics[width=4.5cm]{Plots4b.pdf}
  \caption{} \label{fig:4b}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{5cm}
  \centering
  \includegraphics[width=4.5cm]{Plots4c.pdf}
  \caption{} \label{fig:4c}
\end{subfigure}
\caption{(A) A plot of $\sigma(s)$ for ${\bm \theta}t_0 = {\bm \theta}_0$. Note that $\sigma(s)$ is extra flat at $s = 0$. (B) A plot of $\sigma(s)$ for ${\bm \theta}t_0 = {\bm \theta}_0 + {\bm \epsilon}_1$. Note that $\sigma(s)$ develops a local maximum at $s = 0$. (C) $\sigma(s)$ for ${\bm \theta}t_0 = {\bm \theta}_0 + {\bm \epsilon}_2$. Note that $\sigma(s)$ develops a local minimum at $s = 0$. Again, the dashed green lines are there to help visualize the behavior of $\sigma(s)$ near $s = 0$. Note that all three plots have characteristics similar to the plots of certain quartic polynomials.} \label{fig:4}
\end{figure}

We note that the equations $a_k({\bm \theta}_0,{\bm \gamma}) = 0$ become complex at an alarming rate as $n$ or $k$ increases, and it is for this reason only that we do not push this method farther to search for a higher order $S$-curve.

Next, now that we know the sets \eqref{eq:set1} and \eqref{eq:set2} are non-empty, it is natural to ask if they have any discernible structure or if we can estimate their size. Although we do not provide results in this direction, we demonstrate an interesting relationship in Figure \ref{fig:5}. It would be interesting to determine if the sets \eqref{eq:set1} and \eqref{eq:set2} are closely related to the order parameter in general. Finally, another class of sets that would be interesting to study are
{\bm \epsilon}gin{align*}
\{ {\bf g}ammab \in \mathbb{R}_{>0}^n : \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0 \text{ and } \delta({\bm \theta},{\bf g}ammab) < 0 \}
\end{align*}
for fixed ${\bm \theta} \in \mathbb{R}^n$. In particular, estimating its measure would be an interesting problem.

{\bm \epsilon}gin{figure}[H]
\captionsetup{font=scriptsize}
\centering
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plot5a.pdf}
  \caption{} \label{fig:5a}
\end{subfigure}
{\bm \epsilon}gin{subfigure}{7cm}
  \centering
  \includegraphics[width=7cm]{Plot5b.pdf}
  \caption{} \label{fig:5b}
\end{subfigure}
\caption{(A) A density plot of ${\bf D}elta_p({\bm \theta}) = {\bf D}elta_q({\bm \theta})$. (Note that equality holds almost everywhere by our discussion following Theorem \ref{thm:criterion2}.) (B) A density plot of the order parameter $r$. In both plots we again projected $\mathbb{R}^3$ onto the mean-zero plane ${\bf 1}^\perp$ via the variables ${\bm \phi} = {\bf P} {\bm \theta}$. Note that the white regions in (A) correspond to the white and black regions in (B) and that the black region in (A) corresponds to the blue region in (B). Therefore the region $\{ {\bm \theta} \in \mathbb{R}^3 : {\bf D}elta_p({\bm \theta}) = {\bf D}elta_q({\bm \theta}) = -1\}$ is roughly the set of ${\bm \theta}$ for which $r$ isn't too large or too small.} \label{fig:5}
\end{figure}

\section{Conclusion}

We derived a criterion for the Kuramoto model on a ring network with positive coupling to exhibit a bifurcation in which two branches of phase-locked solutions collide as $\sigma$ increases. Furthermore, we derived a sufficient condition for one of these branches to consist of stable phase-locked solutions. We then applied our criteria to show that for any $n {\bf g}e 3$ we can choose ${\bm \omega}$ and ${\bm \gamma}$ so that this bifurcation occurs. (We require that $n {\bf g}e 5$ to guarantee the existence of ${\bm \omega}$ and ${\bm \gamma}$ for which one of the branches consists of stable phase-locked solutions.) We conjecture that generically our bifurcation is locally a subcritical bifurcation and globally an $S$-curve.

\section{Appendix}

We start by making a definition that significantly shortens the proofs of Lemma \ref{lem:quadratic} and Lemma \ref{lem:st}.

{\bm \epsilon}gin{definition} \label{def:d}
Let ${\bf u} \in \mathbb{R}^n$ and $f : \mathbb{R} \rightarrow \mathbb{R}$ be given. Define the diagonal matrix ${\bf D}_{\bf u} \in \mathbb{R}^{n \times n}$ and the vector $f({\bf u}) \in \mathbb{R}^n$ by
{\bm \epsilon}gin{align*}
({\bf D}_{\bf u})_{ii} = u_i \quad \text{and} \quad f({\bf u})_i = f(u_i).
\end{align*}
\end{definition}

Now we give proofs of all relevant theorems and lemmas.

{\bm \epsilon}gin{proof}[Proof of Theorem \ref{thm:criterion1}]
Note that ${\bf n}(\theta_0,{\bf g}ammab)$ is a normal vector for the level set $\{ {\bm \theta} \in \mathbb{R}^n : \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) = 0\}$ at ${\bm \theta}_0$ pointing in the direction of most rapid increase for $\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))$. From \eqref{eq:hypothesis} we conclude that
{\bm \epsilon}gin{align} \label{eq:detsigma}
\det_\text{Red}({\bf J}b({\bm \theta}(s),{\bf g}ammab))
{\bm \epsilon}gin{cases}
> 0 & \mbox{if $s < 0$ is small,}
\\
= 0 & \mbox{if $s = 0$,}
\\
< 0 & \mbox{if $s > 0$ is small,}
\end{cases}
\quad \text{and} \quad
\frac{d \sigma}{ds}
{\bm \epsilon}gin{cases}
> 0 & \mbox{if $s < 0$ is small,}
\\
= 0 & \mbox{if $s = 0$,}
\\
< 0 & \mbox{if $s > 0$ is small,}
\end{cases}
\end{align}
where ${\bm \theta}(s)$ and $\sigma(s)$ are given in Lemma \ref{lem:flow}. In other words, $\sigma(0) = \sigma_0$ is a local maximum.

Next, by Theorem 2.9 of \cite{doi:10.1137/15M1034258} we have that
{\bm \epsilon}gin{align*}
\text{$\#$positive eigenvalues of ${\bf J}b({\bm \theta},{\bf g}ammab)$} = \# \{ i : {\bf g}amma_i \cos \eta_i > 0 \} -
{\bm \epsilon}gin{cases}
1 & \mbox{ if $\sum_{i=1}^n \frac{1}{{\bf g}amma_i \cos \eta_i} > 0$,}
\\
0 & \mbox{ if $\sum_{i=1}^n \frac{1}{{\bf g}amma_i \cos \eta_i} < 0$.}
\end{cases}
\end{align*}
For all ${\bm \theta}$ in a sufficiently small neighborhood of ${\bm \theta}_0$ we have that $\# \{ i : {\bf g}amma_i \cos \eta_i > 0 \} = n_+(\cos {\bm \eta}_0)$ and that
{\bm \epsilon}gin{align*}
\sum_{i=1}^n \frac{1}{{\bf g}amma_i \cos \eta_i} = \frac{1}{n} \left( \prod_{i=1}^n \frac{1}{{\bf g}amma_i \cos \eta_i} \right) \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))
\end{align*}
has the same sign as $(-1)^{n-n_+(\cos {\bm \eta}_0)} \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab))$.
\end{proof}

{\bm \epsilon}gin{proof}[Proof of Lemma \ref{lem:flow}]
We will only prove \eqref{eq:fixed}. To do this set ${\bf r}(s) := {\bm \omega} - \sigma(s) {\bf g}b({\bm \theta}(s),{\bf g}ammab)$. Then
{\bm \epsilon}gin{align*}
{\bf P} \frac{d {\bf r}}{ds} &= -\sigma'(s) {\bf P} {\bf g}b({\bm \theta}(s),{\bf g}ammab) - \sigma(s) {\bf P} {\bf J}b({\bm \theta}(s),{\bf g}ammab) {\bf v}({\bm \theta}(s),{\bf g}ammab)
\\
&= -\sigma'(s) {\bf P} {\bf g}b({\bm \theta}(s),{\bf g}ammab) + \sigma(s) {\bf P} {\bf J}b({\bm \theta}(s),{\bf g}ammab) {\bf P}^\top \adj({\bf P} {\bf J}b({\bm \theta}(s),{\bf g}ammab) {\bf P}^\top) {\bf P} {\bf g}b({\bm \theta}(s),{\bf g}ammab)
\\
&= [-\sigma'(s) + \sigma(s) \det_\text{Red}({\bf J}b({\bm \theta}(s),{\bf g}ammab)] {\bf P} {\bf g}b({\bm \theta}(s),{\bf g}ammab) = {\bf 0}.
\end{align*}
But ${\bf r}$ is orthogonal to ${\bf 1}$ which spans the nullspace of ${\bf P}$. Therefore we conclude that ${\bf r}(s) = {\bf r}(0) = {\bf 0}$.
\end{proof}

{\bm \epsilon}gin{proof}[Proof of Lemma \ref{lem:quadratic}]

We first note that \eqref{eq:quadratic} follows from the identity
{\bm \epsilon}gin{align} \label{eq:btv}
{\bf B}^\top {\bf v}({\bm \theta},{\bf g}ammab) = \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) \left( \left( \frac{{\bf 1}^\top {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}}{{\bf 1}^\top {\bf D}_{\bf g}ammab^{-1} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}} \right) {\bf D}_{\bf g}ammab^{-1} - {\bf D}_{\sin {\bm \eta}} \right) {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} = {\bf C}({\bm \theta}) {\bf h}({\bf g}ammab)
\end{align}
since the restriction ${\bf B} : {\bf 1}^\perp \rightarrow {\bf 1}^\perp$ is invertible and ${\bf n}({\bm \theta},{\bf g}ammab)^\top = n {\bf h}({\bf g}ammab)^\top \nabla_{\bm \theta} {\bf h}(\cos {\bm \eta}) = {\bf h}({\bf g}ammab)^\top {\bf A}({\bm \theta})$. For reference later we note the identities
{\bm \epsilon}gin{align} \label{eq:gJ}
{\bf g}b({\bm \theta},{\bf g}ammab) = {\bf B} {\bf D}_{\bf g}ammab \sin {\bm \eta} \quad \text{and} \quad {\bf J}b({\bm \theta},{\bf g}ammab) = {\bf B} {\bf D}_{\bf g}ammab {\bf D}_{\cos {\bm \eta}} {\bf B}^\top.
\end{align}
We start from the equation $({\bf P} {\bf J}b({\bm \theta},{\bf g}ammab) {\bf P}^\top) {\bf P} {\bf v}({\bm \theta},{\bf g}ammab) = -\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) {\bf P} {\bf g}b({\bm \theta},{\bf g}ammab)$ which simplifies to ${\bf J}b({\bm \theta},{\bf g}ammab) {\bf v}({\bm \theta},{\bf g}ammab) \\ = -\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) {\bf g}b({\bm \theta},{\bf g}ammab)$ since ${\bf P} {\bf P}^\top = {\bf I}b$ and ${\bf P}^\top {\bf P} = {\bf I}b - \frac{1}{n} {\bf 1} {\bf 1}^\top$. By \eqref{eq:gJ} we have that
{\bm \epsilon}gin{align*}
{\bf B} {\bf D}_{\bf g}ammab {\bf D}_{\cos {\bm \eta}} {\bf B}^\top {\bf v}({\bm \theta},{\bf g}ammab) = -\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) {\bf B} {\bf D}_{\bf g}ammab {\bf D}_{\sin {\bm \eta}} {\bf 1}
\end{align*}
hence there exists a scalar $x$ such that
{\bm \epsilon}gin{align*}
{\bf B}^\top {\bf v}({\bm \theta},{\bf g}ammab) = -\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} + x {\bf D}_{\bf g}ammab^{-1} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}
\end{align*}
since the null space of ${\bf B}$ is spanned by ${\bf 1}$. In fact we find that
{\bm \epsilon}gin{align*}
x = \det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) \frac{{\bf 1}^\top {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}}{{\bf 1}^\top {\bf D}_{\bf g}ammab^{-1} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}}
\end{align*}
by noting that ${\bf 1}^\top {\bf B}^\top = {\bf 0}$. This verifies the first half of \eqref{eq:btv}. Next we note that
{\bm \epsilon}gin{align*}
\det_{\text{Red}}({\bf J}b({\bm \theta},{\bf g}ammab)) {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} = n {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} {\bf h}(\cos {\bm \eta})^\top {\bf h}({\bf g}ammab)
\end{align*}
and that
{\bm \epsilon}gin{align*}
\det_\text{Red}({\bf J}b({\bm \theta},{\bf g}ammab)) \left( \frac{{\bf 1}^\top {\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}}{{\bf 1}^\top {\bf D}_{\bf g}ammab^{-1} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1}} \right) {\bf D}_{\bf g}ammab^{-1} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} = n\left( \sum_{k=1}^n \frac{\sin \eta_k}{\cos \eta_k} \right) {\bf D}_{{\bf h}(\cos {\bm \eta})} {\bf h}({\bf g}ammab)
\end{align*}
by Theorem \ref{thm:detred} where
{\bm \epsilon}gin{align*}
n ({\bf D}_{\sin {\bm \eta}} {\bf D}_{\cos {\bm \eta}}^{-1} {\bf 1} {\bf h}(\cos {\bm \eta})^\top)_{ij} = n \frac{\sin \eta_i}{\cos \eta_i} {\bf h}(\cos {\bm \eta})_j
\end{align*}
and
{\bm \epsilon}gin{align*}
n \left( \sum_{k=1}^n \frac{\sin \eta_k}{\cos \eta_k} \right) ({\bf D}_{{\bf h}(\cos {\bm \eta})})_{ij} =
{\bm \epsilon}gin{cases}
n \left( \sum_{k=1}^n \frac{\sin \eta_k}{\cos \eta_k} \right) {\bf h}(\cos {\bm \eta})_j & \mbox{if $i = j$,}
\\
0 & \mbox{if $i \ne j$,}
\end{cases}
\end{align*}
which verifies the second half of \eqref{eq:btv}.
\end{proof}

{\bm \epsilon}gin{proof}[Proof of Lemma \ref{lem:min}]
Following the arguments in \cite{KREPS1984105} we easily find that
{\bm \epsilon}gin{align*}
\min_{{\bf x}b \in \mathbb{R}_{{\bf g}e 0}^n} \frac{{\bf x}b^\top {\bf R} {\bf x}b}{\| {\bf x}b \|_1^2} {\bf g}e t \text{ if and only if } \min_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i}) \le 0 \text{ or } \frac{\det({\bf R}_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i})} {\bf g}e t
\end{align*}
for every ${\bf I}c \subseteq \{1, \dots, n\}$. This is of course equivalent to
{\bm \epsilon}gin{align*}
\min_{{\bf x}b \in \mathbb{R}_{{\bf g}e 0}^n} \frac{{\bf x}b^\top {\bf R} {\bf x}b}{\| {\bf x} \|_1^2} < t \text{ if and only if } \min_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i}) > 0 \text{ and } \frac{\det({\bf R}_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_{{\bf I}c,i})} < t
\end{align*}
for some ${\bf I}c \subseteq \{1, \dots, n\}$. Taking the infimum over all such $t$ gives the result.
\end{proof}

{\bm \epsilon}gin{proof}[Proof of Theorem \ref{thm:criterion2}]
Suppose that ${\bf D}elta_p({\bm \theta}_0) = -1$. Then there exists an ${\bf I}c \in {\bf I}c({\bm \theta}_0)$ such that
{\bm \epsilon}gin{align*}
\lim_{\tau \rightarrow \infty} \frac{\det({\bf R}_\tau({\bm \theta}_0)_{\bf I}c)}{\sum_{i \in {\bf I}c} \det({\bf R}_\tau({\bm \theta}_0)_{{\bf I}c,i})} = \epsilon < 0.
\end{align*}
By Lemma \ref{lem:min} there exists vectors ${\bf x}b_\tau \in \mathbb{R}_{{\bf g}e 0}^n$ such that $\| {\bf x}b_\tau \|_1 = 1$ and
{\bm \epsilon}gin{align*}
{\bf x}b_\tau^\top {\bf R}_\tau({\bm \theta}_0) {\bf x}b_\tau = ({\bf h}(\cos {\bm \eta}_0)^\top {\bf x}b_\tau)^2 \tau + {\bf x}b_\tau^\top {\bf T}({\bm \theta}_0) {\bf x}b_\tau < \frac{\epsilon}{2}
\end{align*}
for all sufficiently large $\tau$. But the unit ball $\{ {\bf x}b \in \mathbb{R}^n : \| {\bf x}b \|_1 = 1\}$ is compact so there exists an infinite sequence $\tau_1 < \tau_2 < \dots$ with $\tau_n \rightarrow \infty$ such that $\lim_{n \rightarrow \infty} {\bf x}b_{\tau_n} = {\bf x}b$ exists. Now the second term ${\bf x}b_\tau^\top {\bf T}({\bm \theta}_0) {\bf x}b_\tau$ is bounded since $\| {\bf x}b_\tau \|_1 = 1$. Therefore ${\bf h}(\cos {\bm \eta}_0)^\top {\bf x}b = \lim_{n \rightarrow \infty} {\bf h}(\cos {\bm \eta}_0)^\top {\bf x}b_{\tau_n} = 0$. Also, the first term $({\bf h}(\cos {\bm \eta}_0)^\top {\bf x}b_\tau)^2$ is non-negative so that ${\bf x}b^\top {\bf T}({\bm \theta}_0) {\bf x}b = \lim_{n \rightarrow \infty} {\bf x}b_{\tau_n}^\top {\bf T}({\bm \theta}_0) {\bf x}b_{\tau_n} \le \epsilon/2$.

Clearly ${\bf x}b \in \mathbb{R}_{{\bf g}e 0}^n$ and $\| {\bf x}b \|_1 = 1$. Therefore ${\bf h}(\cos {\bm \eta}_0)^\top {\bf x}b = 0$ implies that ${\bf h}(\cos {\bm \eta}_0)$ has at least one positive component and at least one negative component. Therefore we can slightly modify ${\bf x}b$ to obtain a new vector ${\bf x}bt \in \mathbb{R}_{>0}^n$ such that ${\bf h}(\cos {\bm \eta}_0)^\top {\bf x}bt = 0$ and ${\bf x}bt^\top {\bf T}({\bm \theta}_0) {\bf x}bt < 0$. By Lemma \ref{lem:quadratic} we can choose ${\bf g}ammab \in \mathbb{R}_{>0}^n$ such that ${\bf x}bt = {\bf h}({\bf g}ammab)$. Then $\det_\text{Red}({\bf J}b({\bm \theta}_0,{\bf g}ammab)) = {\bf h}(\cos {\bm \eta}_0)^\top {\bf x}bt = 0$ and $\delta({\bm \theta}_0,{\bf g}ammab) = {\bf x}bt^\top {\bf T}({\bm \theta}_0) {\bf x}bt < 0$ so that \eqref{eq:hypothesis} is satisfied.

Now suppose that there exists a ${\bf g}ammab$ satisfying \eqref{eq:hypothesis}. Then the vector ${\bf x}b = {\bf h}({\bf g}ammab) \in \mathbb{R}_{>0}^n$ satisfies
{\bm \epsilon}gin{align*}
{\bf x}b^\top {\bf R}_\tau({\bm \theta}_0) {\bf x}b = \det_\text{Red}({\bf J}b({\bm \theta}_0,{\bf g}ammab))^2 \tau + \delta({\bm \theta}_0,{\bf g}ammab) < 0
\end{align*}
for all $\tau$. By Lemma \ref{lem:min} this implies that for each $\tau$ there exists an ${\bf I}c_\tau \in {\bf I}c({\bm \theta}_0)$ such that $\det({\bf R}_\tau({\bm \theta}_0)_{{\bf I}c_\tau}) < 0$. But ${\bf I}c({\bm \theta}_0)$ is finite so there exits an ${\bf I}c \in {\bf I}c({\bm \theta}_0)$ such that $\det({\bf R}_{\tau_n}({\bm \theta}_0)_{\bf I}c) < 0$ for a sequence $\tau_1 < \tau_2 < \dots $ with $\tau_n \rightarrow \infty$. However, $\det({\bf R}_{\tau_n}({\bm \theta}_0)_{\bf I}c)$ is a polynomial in $\tau$ so that in fact $\det({\bf R}_{\tau_n}({\bm \theta}_0)_{\bf I}c) < 0$ for all sufficiently large $\tau$. Therefore we obtain that ${\bf D}elta_q({\bm \theta}_0) = -1$.
\end{proof}

{\bm \epsilon}gin{proof}[Proof of Lemma \ref{lem:st}]
We first note that ${\bm \eta}t = (\eta_1,\dots,\eta_n,0,\dots,0) \in \mathbb{R}^{{\tilde{n}}}$ where ${\bm \eta} = (\eta_1,\dots,\eta_n) \in \mathbb{R}^n$. Therefore
{\bm \epsilon}gin{align*}
{\bf h}(\cos {\bm \eta}t)_i =
{\bm \epsilon}gin{cases}
{\bf h}(\cos {\bm \eta})_i & \mbox{if $1 \le i \le n$,}
\\
\prod_{j=1}^n \cos \eta_j & \mbox{if $n < i \le {\tilde{n}}$,}
\end{cases}
\end{align*}
demonstrating the first equality. Next we note that
{\bm \epsilon}gin{align*}
{\bf C}t({\bm \theta}t) =
{\bm \epsilon}gin{pmatrix}
\frac{{\tilde{n}}}{n} {\bf C}({\bm \theta}) & {\bf h}dots
\\
{\bf 0} & \ddots
\end{pmatrix}
\end{align*}
from Definition \ref{def:abc}. Also we have that
{\bm \epsilon}gin{align*}
{\bf H}t(\cos {\bm \eta}t) =
{\bm \epsilon}gin{pmatrix}
{\bf H}(\cos {\bm \eta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
\quad \text{and} \quad {\bf A}t({\bm \theta}t) = {\tilde{n}} {\bf H}t(\cos {\bm \eta}t) {\bf D}t_{\sin {\bm \eta}t} {\bf B}t^\top \quad \text{and} \quad {\bf A}({\bm \theta}) = n {\bf H}(\cos {\bm \eta}) {\bf D}_{\sin {\bm \eta}} {\bf B}^\top
\end{align*}
where ${\bf H}t$ and ${\bf H}$ are the Jacobians of ${\bf h}t$ and ${\bf h}$ respectively. Therefore
{\bm \epsilon}gin{align*}
{\bf T}t({\bm \theta}t) &= {\tilde{n}} {\bf H}t(\cos {\bm \eta}t) {\bf D}t_{\sin {\bm \eta}t} {\bf B}t^\top ({\bf B}t^\top)^{-1} {\bf C}t({\bm \theta}t) =
{\bm \epsilon}gin{pmatrix}
\frac{{\tilde{n}}}{n} {\bf H}(\cos {\bm \eta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
{\bf D}_{\sin {\bm \eta}} & {\bf 0}
\\
{\bf 0} & {\bf 0}
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
\frac{{\tilde{n}}}{n} {\bf C}({\bm \theta}) & {\bf h}dots
\\
{\bf 0} & \ddots
\end{pmatrix}
\\
&=
{\bm \epsilon}gin{pmatrix}
(\frac{{\tilde{n}}}{n})^2{\bf H}(\cos {\bm \eta}) {\bf D}_{\sin {\bm \eta}} {\bf C}({\bm \theta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
=
{\bm \epsilon}gin{pmatrix}
(\frac{{\tilde{n}}}{n})^2{\bf H}(\cos {\bm \eta}) {\bf D}_{\sin {\bm \eta}} {\bf B} ({\bf B}^\top)^{-1} {\bf C}({\bm \theta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
=
{\bm \epsilon}gin{pmatrix}
(\frac{{\tilde{n}}}{n})^2{\bf T}({\bm \theta}) & {\bf h}dots
\\
{\bf v}dots & \ddots
\end{pmatrix}
\end{align*}
completing the proof.
\end{proof}

\iffalse
{\bm \epsilon}gin{proof}[Proof of Lemmas \ref{lem:minpos} and \ref{lem:minall}]
Next, Lemma \ref{lem:minall} follows by computing the maximum over each octant which can be indexed by ${\bm \epsilon}$ using the inequalities $\epsilon_i x_i {\bf g}e 0$. Therefore we can apply Lemma $\ref{lem:minpos}$ to the matrix $\tilde{{\bf S}}$ where $\tilde{s}_{ij} = \epsilon_i \epsilon_j s_{ij}$. Then for any principal minor $\tilde{{\bf M}}$ of $\tilde{{\bf S}}$ we find that $\det(\tilde{{\bf M}}) = \det({\bf M})$ and $\tilde{{\bf M}}_{ij} = \epsilon_i \epsilon_j {\bf M}_{ij}$ where ${\bf M}$ is the principal minor of ${\bf S}$ corresponding to the principal minor $\tilde{{\bf M}}$ of $\tilde{{\bf S}}$.
\end{proof}
\fi

\iffalse
{\bm \epsilon}gin{proof}[Proof of Lemma \ref{lem:nonempty}]
Let $N {\bf g}e 3$ and set ${\bm \theta}_\epsilon = (\pi,0,\dots,0) + \epsilon {\bf u}$ for some vector ${\bf u}$ to be chosen later. Then ${\bm \eta}_\epsilon = (-\pi,0,\dots,0,\pi) +\epsilon {\bf v}$ where ${\bf v} = {\bf B}^\top {\bf u}$. Therefore
{\bm \epsilon}gin{align*}
\sin {\bm \eta}_\epsilon = (-v_1,v_2,\dots,v_{N-1},-v_N)^\top \epsilon + O(\epsilon^3) \quad \text{and} \quad \cos {\bm \eta}_\epsilon = (-1,1,\dots,1,-1)^\top + O(\epsilon^2).
\end{align*}
from which we conclude that
{\bm \epsilon}gin{align*}
{\bf h}(\cos {\bm \eta}_\epsilon) = (-1,1,\dots,1,-1) + O(\epsilon^2) \quad \text{and} \quad \frac{\sin \eta_i}{\cos \eta_i} = v_i \epsilon + O(\epsilon^3).
\end{align*}
Therefore
{\bm \epsilon}gin{align*}
{\bf h}(\cos {\bm \eta}_\epsilon) {\bf h}(\cos {\bm \eta}_\epsilon)^\top =
{\bm \epsilon}gin{pmatrix}
1 & -1 & {\bf h}dots
\\
-1 & 1 & {\bf h}dots &
\\
{\bf v}dots & {\bf v}dots & \ddots
\end{pmatrix}
+ O(\epsilon^2)
\end{align*}
and
{\bm \epsilon}gin{align*}
{\bf C}({\bm \eta}_\epsilon) = {\bf v} {\bf 1}^\top \epsilon + O(\epsilon^3)
\end{align*}
where we have used that ${\bf v}$ is mean zero to rewrite the diagonal terms. Therefore we have that
{\bm \epsilon}gin{align*}
({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}_\epsilon) = ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf B}^\top {\bf u} {\bf 1}^\top \epsilon + O(\epsilon^3) = {\bf u} {\bf 1}^\top \epsilon + O(\epsilon^3)
\end{align*}
where we have used that ${\bf v}$ is mean zero to rewrite the diagonal terms in the second identity. Also
{\bm \epsilon}gin{align*}
({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}_\epsilon) = ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf B}^\top {\bf u} {\bf 1}^\top \epsilon + O(\epsilon^3) = {\bf u} {\bf 1}^\top \epsilon + O(\epsilon^3)
\end{align*}
if ${\bf u}$ is chosen to be mean zero. It is also not hard to see that ${\bf A}({\bm \eta}_\epsilon) = {\bf A}_0 \epsilon + O(\epsilon^3)$ (see computation below) from which we conclude that
{\bm \epsilon}gin{align*}
{\bf A}({\bm \eta}_\epsilon) ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}_\epsilon) =
{\bm \epsilon}gin{pmatrix}
({\bf A}_0 {\bf u})_1 & ({\bf A}_0) {\bf u})_1 & {\bf h}dots
\\
({\bf A}_0 {\bf u})_2 & ({\bf A}_0) {\bf u})_2 & {\bf h}dots
\\
{\bf v}dots & {\bf v}dots & \ddots
\end{pmatrix}
\epsilon^2 + O(\epsilon^4).
\end{align*}
Now if
{\bm \epsilon}gin{align*}
{\bf S}_\tau({\bm \eta}_\epsilon) =
{\bm \epsilon}gin{pmatrix}
s_{11} & s_{12} & {\bf h}dots
\\
s_{21} & s_{22} & {\bf h}dots
\\
{\bf v}dots & {\bf v}dots & \ddots
\end{pmatrix}
\tau +
{\bm \epsilon}gin{pmatrix}
t_{11} & t_{12} & {\bf h}dots
\\
t_{21} & t_{22} & {\bf h}dots
\\
{\bf v}dots & {\bf v}dots & \ddots
\end{pmatrix},
\end{align*}
and ${\bf M}$ denotes the upper left $2 \times 2$ minor of ${\bf S}_\tau({\bm \eta}_\epsilon)$, we find that
{\bm \epsilon}gin{align*}
d = \det {\bm \epsilon}gin{pmatrix} s_{11} & t_{12} \\ s_{21} & t_{22} \end{pmatrix} + \det {\bm \epsilon}gin{pmatrix} t_{11} & s_{12} \\ t_{21} & s_{22} \end{pmatrix}
, \quad d_{11} = s_{22} , \quad d_{12} = -s_{21}, \quad d_{21} = -s_{12}, \quad d_{22} = s_{11}.
\end{align*}
Since
{\bm \epsilon}gin{align*}
{\bf S}_\tau({\bm \eta}_0) =
{\bm \epsilon}gin{pmatrix}
1 & -1 & {\bf h}dots
\\
-1 & 1 & {\bf h}dots
\\
{\bf v}dots & {\bf v}dots & \ddots
\end{pmatrix}
\tau
\end{align*}
we conclude that
{\bm \epsilon}gin{align*}
d_{11} + d_{12} = s_{22} - s_{21} = 1 - (-1) = 2 > 0 \quad \text{and} \quad d_{21} + d_{22} = -s_{12} + s_{11} = - (-1) + 1 = 2 > 0
\end{align*}
at $\epsilon = 0$. Of course these quantities are continuous in $\epsilon$ and therefore remain positive for all sufficiently small values of $\epsilon$. Therefore the minor ${\bf M}$ does contribute to the value of $D$ and it remains to show that $d < 0$. Now
{\bm \epsilon}gin{align*}
d = \det {\bm \epsilon}gin{pmatrix} 1 & ({\bf A}_0 {\bf u})_1 \\ -1 & ({\bf A}_0 {\bf u})_2 \end{pmatrix} \epsilon^4 + \det {\bm \epsilon}gin{pmatrix} ({\bf A}_0 {\bf u})_1 & -1 \\ ({\bf A}_0 {\bf u})_2 & 1 \end{pmatrix} \epsilon^4 + O(\epsilon^8) = 2(({\bf A}_0 {\bf u})_1 + ({\bf A}_0 {\bf u})_2) \epsilon^4 + O(\epsilon^8).
\end{align*}

Therefore it remains to compute the first two rows of ${\bf A}_0$. To do this
{\bm \epsilon}gin{align*}
{\bf h}(\cos {\bm \eta})_1 = \cos \eta_2 \dots \cos \eta_N
\end{align*}
so that
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_1} {\bf h}(\cos {\bm \eta})_1 = -\cos \eta_2 \dots \cos \eta_{N-1} \sin \eta_N = v_N \epsilon + O(\epsilon^3)
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_2} {\bf h}(\cos {\bm \eta})_1 = \sin \eta_2 \cos \eta_3 \dots \cos \eta_N = -v_2 \epsilon + O(\epsilon^3)
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} {\bf h}(\cos {\bm \eta})_1 = \left( \prod_{k \ne 1,j} \cos \eta_k \right) \sin \eta_j -\left( \prod_{k \ne 1,j-1} \cos \eta_k \right) \sin \eta_{j-1} = (-v_j + v_{j-1})\epsilon + O(\epsilon^3).
\end{align*}
for $j \ne 1,2$. Similarly
{\bm \epsilon}gin{align*}
{\bf h}(\cos {\bm \theta})_2 = \cos \eta_1 \cos \eta_3 \dots \cos \eta_N
\end{align*}
so that
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_2} {\bf h}(\cos {\bm \eta})_2 = -\sin \eta_1 \cos \eta_3 \dots \cos \eta_N = v_2 \epsilon + O(\epsilon^3)
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_3} {\bf h}(\cos {\bm \eta})_2 = \cos \eta_1 \sin \eta_3 \cos \eta_4 \dots \cos \eta_N = v_3 \epsilon + O(\epsilon^3)
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} {\bf h}(\cos {\bm \eta})_N = \left( \prod_{k \ne 2,j} \cos \eta_k \right) \sin \eta_j - \left( \prod_{k \ne 2,j-1} \cos \eta_k \right) \sin \eta_{j-1} = (v_j -v_{j-1})\epsilon + O(\epsilon^3)
\end{align*}
for $j \ne 2,3$. (Here we interpret $j-1$ at $j=1$ as $N$.) Therefore
{\bm \epsilon}gin{align*}
{\bf A}_0 =
{\bm \epsilon}gin{pmatrix}
v_N & -v_2 & -v_3 + v_2 & -v_4 + v_3 & {\bf h}dots & -v_N + v_{N-1}
\\
v_1 - v_N & v_2 & v_3 & v_4 - v_3 & {\bf h}dots & v_N -v_{N-1}
\\
{\bf v}dots & {\bf v}dots & {\bf v}dots & {\bf v}dots & \ddots & {\bf v}dots
\end{pmatrix}
\end{align*}
and
{\bm \epsilon}gin{align*}
({\bf A}_0 {\bf u})_1 + ({\bf A}_0 {\bf u})_2 = (v_1, 0, v_2, 0 , \dots ,0)^\top {\bf u} = (u_2 - u_1, 0 , u_3 - u_2, 0 , \dots ,0)^\top {\bf u} = (u_2-  u_1)u_1 +  (u_3 - u_2) u_3.
\end{align*}
If we choose the mean zero vector ${\bf u} = (u_1, -u_1 - u_3,0,\dots,0)^\top$, then we arrive at the quadratic form
{\bm \epsilon}gin{align*}
({\bf A}_0 {\bf u})_1 + ({\bf A}_0 {\bf u})_2 = 2(u_3^2 - u_1^2)
\end{align*}
which is of course negative if $|u_3| < |u_1|$.
\end{proof}
\fi

\iffalse
Let ${\bf u} \in {\bf 1}^\perp$ and set ${\bm \theta} = \epsilon {\bf u}$. Then $\sin {\bm \eta}_i = \epsilon(u_{i+1} - u_i) + O(\epsilon^3)$ and $\cos \eta_i = 1 + O(\epsilon^2)$. Therefore
{\bm \epsilon}gin{align*}
{\bf C}({\bm \theta})_{ij} =
{\bm \epsilon}gin{cases}
-n \epsilon \left( \sum_{k \ne i} u_{k+1} - u_k \right) + O(\epsilon^3) = n\epsilon(u_{i+1} - u_i) + O(\epsilon^3) & \mbox{if $i = j$,}
\\
n \epsilon(u_{i+1} - u_i) + O(\epsilon^3) & \mbox{if $i \ne j$,}
\end{cases}
\end{align*}
so that ${\bf C}({\bm \theta}) = n \epsilon {\bf B}^\top {\bf u} {\bf 1}^\top$. Now
{\bm \epsilon}gin{align*}
(-{\bf B} {\bf B}^\top {\bf u} {\bf 1}^\top)_{ij} = u_{j+1} - 2u_j + u_{j-1}
\end{align*}
and
{\bm \epsilon}gin{align*}
(-{\bf D}_{{\bf B}^\top {\bf u}} {\bf B}^\top)_{ij} =
{\bm \epsilon}gin{cases}
u_{j+1} - u_j & \mbox{if $i = j$,}
\\
-u_j + u_{j-1} & \mbox{if $i = j-1$,}
\\
0 & \mbox{otherwise,}
\end{cases}
\end{align*}

Similarly,
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} \prod_{k=1}^n \cos \eta_k = \left( \frac{\partial}{\partial \theta_j} \cos \eta_i \right) {\bf h}(\cos {\bm \eta})_i + \cos \eta_i \frac{\partial}{\partial \theta_j} {\bf h}(\cos {\bm \eta})_i
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} \cos \eta_i =
{\bm \epsilon}gin{cases}
\epsilon(u_{j+1} - u_j) + O(\epsilon^3) & \mbox{if $i = j$,}
\\
-\epsilon(u_j - u_{j-1}) + O(\epsilon^3) & \mbox{if $i = j-1$,}
\\
0 & \mbox{otherwise,}
\end{cases}
\end{align*}
and
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} \prod_{k=1}^n \cos \eta_k &= \left( \frac{\partial}{\partial \theta_j} \cos(\theta_{j+1} - \theta_j) \cos(\theta_j - \theta_{j-1}) \right) \prod_{k \ne j-1,j} \cos \eta_k
\\
&= \left( \sin(\theta_{j+1} - \theta_j) \cos(\theta_j - \theta_{j-1}) - \cos(\theta_{j+1} - \theta_j) \sin(\theta_j - \theta_{j-1}) \right) \prod_{k \ne j-1,j} \cos \eta_k
\\
&= \left( \frac{\sin(\theta_{j+1} - \theta_j)}{\cos(\theta_{j+1} - \theta_j)} - \frac{\sin(\theta_j - \theta_{j-1})}{\cos(\theta_j - \theta_{j-1})} \right) \prod_{k=1}^n \cos \eta_k
\\
&= \epsilon(u_{j+1} - 2u_j + u_{j-1}) + O(\epsilon^3).
\end{align*}
Therefore
{\bm \epsilon}gin{align*}
\frac{\partial}{\partial \theta_j} {\bf h}(\cos {\bm \eta})_i =
{\bm \epsilon}gin{cases}
\epsilon(u_{j+1} - u_j) + O(\epsilon^3) & \mbox{if $i = j$,}
\\
-\epsilon(u_j - u_{j-1}) + O(\epsilon^3) & \mbox{if $i = j-1$,}
\\
0 & \mbox{otherwise,}
\end{cases}
\end{align*}
\fi

\iffalse
By a theorem due to Hilbert we know that each $p_{\bm \theta}$ is a sum of at most three quadratic forms being a ternary quartic polynomial. Since $p_{\bm \theta}$ is homogenous with even degree we know that it is the sum of at most three squares of homogenous polynomials of degree two. Define $P$ to be the convex hull of $\{(1,1,0),(1,0,1),(0,1,1)\}$, then $\newt(p_{\bm \theta}) = 2P$. There $\newt(p_i) \subseteq P$. The only monomials that are homogenous of contained in the Newton polytope $P$ are the vertices of $P$ which represent ${\bf g}amma_{12} {\bf g}amma_{23}$, ${\bf g}amma_{23} {\bf g}amma_{31}$, and ${\bf g}amma_{31} {\bf g}amma_{12}$.

{\bm \epsilon}gin{proof}
If $G$ is a star network, then
{\bm \epsilon}gin{align*}
{\bf v}v({\bm \theta}) = ({\bf 1} \otimes {\bf 1}-(n-1) {\bf I}) \w({\bm \theta}) \quad \text{where} \quad w({\bm \theta})_i =
{\bm \epsilon}gin{cases}
0 & \mbox{if $i=1$,}
\\
\sin (\theta_i - \theta_1) \prod_{j \ne i} \cos (\theta_j - \theta_1) & \mbox{if $i > 1$.}
\end{cases}
\end{align*}
Notice that ${\bf M}$ is an invertible matrix since it has eigenvalues $1$ and $-(n-1)$ with multiplicity $1$ and $n-1$ respectively.
\end{proof}

{\bm \epsilon}gin{align*}
({\bf P}P^\top \adj({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top) {\bf P}P {\bf g}({\bm \theta}))_i &= \sum_{j=1}^n {\bf P}P_{ij}^\top (\adj({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top) {\bf P}P {\bf g}({\bm \theta}))_j
\\
&= \sum_{j,k=1}^n {\bf P}P_{ij}^\top \adj({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top)_{jk} ({\bf P}P {\bf g}({\bm \theta}))_k
\\
&= \sum_{j,k,\ell =1}^n {\bf P}P_{ij}^\top \adj({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top)_{jk} {\bf P}P_{k\ell} {\bf g}({\bm \theta})_\ell
\\
&= \sum_{j,k,\ell =1}^n (-1)^{j+k} {\bf P}P_{ij}^\top {\bf P}P_{k\ell} \det(({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top)_{kj}) {\bf g}({\bm \theta})_\ell
\end{align*}
and
{\bm \epsilon}gin{align*}
\det(({\bf P}P {\bf J}({\bm \theta}) {\bf P}P^\top)_{kj}) &= \det({\bf P}P_k {\bf J}({\bm \theta}) {\bf P}P_j^\top)
\\
&= \sum_{ \{ r,s \} \subseteq [n]} \det({\bf P}P_{k,\{r,s\}}) \det(({\bf J}({\bm \theta}) {\bf P}P_j^\top)_{\{r,s\}})
\\
&= \sum_{ \{ r,s \} \subseteq [n]} \det({\bf P}P_{k,\{r,s\}}) \sum_{\{t,u\} \subseteq [n]} \det({\bf J}({\bm \theta})_{\{r,s\},\{t,u\}}) \det({\bf P}P_{j,\{t,u\}})
\end{align*}

Define
{\bm \epsilon}gin{align*}
d_i =
{\bm \epsilon}gin{cases}
0 & \mbox{if $i=1$,}
\\
\prod_{j \notin \{i,1\}} \cos(\theta_j - \theta_1) & \mbox{if $i>1$.}
\end{cases}
\end{align*}

Conjecture
{\bm \epsilon}gin{align*}
\det({\bf J}({\bm \theta})_{\{r,s\},\{t,u\}}) = \pm \sum_{i \in \{r,s\} \cap \{t,u\}} d_i
\end{align*}
\fi

\iffalse
Let ${\bm \theta}'(s) = F({\bm \theta}(s))$ with ${\bm \theta}(0) = {\bm \theta}_0$ and define
{\bm \epsilon}gin{align*}
{\bm \theta}(s) = \alpha(s) {\bf u}_0 + {\bm \epsilon}ta(s) {\bf v}_0 \quad \text{where} \quad {\bf u}_0 = {\bm \epsilon}gin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} {\bf v}_0 \quad \text{and} \quad {\bf v}_0 = \nabla \det_{\text{Red}}({\bf J}({\bm \theta}_0)).
\end{align*}
Then
{\bm \epsilon}gin{align*}
\frac{d{\bm \epsilon}ta}{d\alpha} = \frac{{\bm \epsilon}ta'(s)}{\alpha'(s)} \quad \text{where} \quad
{\bm \epsilon}gin{pmatrix}
\alpha'(s)
\\
{\bm \epsilon}ta'(s)
\\
0
\end{pmatrix}
=
{\bm \epsilon}gin{pmatrix}
{\bf u}_0 & {\bf v}_0 & {\bf 1}
\end{pmatrix}^{-1}
F({\bm \theta}(s)).
\end{align*}
\fi

\iffalse
Define the order parameter
{\bm \epsilon}gin{align*}
r(t) = \biggr\rvert \frac{1}{n} \sum_{i=1}^n e^{\theta_i \sqrt{-1}} \biggr\rvert
\end{align*}
then
{\bm \epsilon}gin{align*}
r(s)^2 = \frac{1}{n^2} \sum_{i,j=1}^n e^{(\theta_i(s)-\theta_j(s)) \sqrt{-1}}.
\end{align*}
Therefore
{\bm \epsilon}gin{align*}
(r(s)^2)' &= 2 r(s) r'(s) = \frac{1}{n^2} \sum_{i,j=1}^n e^{(\theta_i(s)-\theta_j(s)) \sqrt{-1}} (\theta_i'(s) - \theta_j'(s)) \sqrt{-1}
\\
&= \frac{2}{n^2} \frac{Q(s) + \overline{Q(s)}}{2} = \frac{2}{n^2} \mathbb{R}e Q(s) = -\frac{2}{n^2} \sum_{i,j=1}^n \sin(\theta_i(s) - \theta_j(s)) \theta_i'(s) = -\frac{2}{n^2} \sum_{i=1}^n \left( \sum_{j=1}^n \sin(\theta_i(s) - \theta_j(s)) \right) \theta_i'(s)
\end{align*}
where
{\bm \epsilon}gin{align*}
Q(s) = \sum_{i,j=1}^n e^{(\theta_i(s)-\theta_j(s)) \sqrt{-1}} \theta_i'(s) \sqrt{-1} \quad \text{and} \quad \overline{Q(s)} = -\sum_{i,j=1}^n e^{(\theta_i(s)-\theta_j(s)) \sqrt{-1}} \theta_j'(s) \sqrt{-1}.
\end{align*}
Define
{\bm \epsilon}gin{align*}
v_i = -\frac{1}{n^2} \sum_{j=1}^n \sin(\theta_i(s) - \theta_j(s))
\end{align*}
and recall that
{\bm \epsilon}gin{align*}
\theta'(s) = \det(J) P^\top(P J(\theta) P^\top)^{-1} P g(\theta)
\end{align*}
so that
{\bm \epsilon}gin{align*}
r(s) r'(s) = \det(J(s)) v(\theta(s))^\top P^\top(P J(\theta(s)) P^\top)^{-1} P g(\theta(s))
\end{align*}
\fi

\iffalse
Consider
{\bm \epsilon}gin{align*}
\dot{\theta}_1 &= \omega_1 - \sigma \sin(\theta_1 -  \theta_2)
\\
\dot{\theta}_2 &= \omega_2 - \sigma \sin(\theta_2 -  \theta_1)
\end{align*}
If we define ${\bf v}arphi = \theta_2 - \theta_1$ and ${\bf D}elta = \omega_2 - \omega_1$ we arrive at
{\bm \epsilon}gin{align*}
\dot{{\bf v}arphi} = {\bf D}elta - 2\sigma \sin {\bf v}arphi.
\end{align*}
Therefore our fixed point equation becomes
{\bm \epsilon}gin{align*}
{\bf D}elta - 2\sigma \sin {\bf v}arphi = 0.
\end{align*}
Differentiating with respect to $s$ gives
{\bm \epsilon}gin{align*}
\sigma' \sin {\bf v}arphi + \sigma \cos {\bf v}arphi \dot{{\bf v}arphi}.
\end{align*}
from which we deduce that
{\bm \epsilon}gin{align*}
\dot{{\bf v}arphi} = -\frac{\sigma'}{\sigma} \frac{\sin {\bf v}arphi}{\cos {\bf v}arphi}
\end{align*}
We choose the flow
{\bm \epsilon}gin{align*}
\dot{{\bf v}arphi} = \sin {\bf v}arphi
\end{align*}

{\bm \epsilon}gin{example}
We give a general proof in the appendix but here consider the ``all-to-all" Kuramoto model with $N = 3$, namely,
{\bm \epsilon}gin{align*}
\theta_1'(t) &= \omega_1 - \sigma[ \sin(\theta_1(t) - \theta_2(t)) + \sin( \theta_1(t) - \theta_3(t)],
\\
\theta_2'(t) &= \omega_2 - \sigma[ \sin(\theta_2(t) - \theta_1(t)) + \sin( \theta_2(t) - \theta_3(t)],
\\
\theta_3'(t) &= \omega_3 - \sigma[ \sin(\theta_3(t) - \theta_1(t)) + \sin( \theta_3(t) - \theta_2(t)].
\end{align*}
This system can be reduced by a dimension via the change of variables $\psi_1 = \theta_2 - \theta_1$, $\psi_2 = \theta_3 - \theta_2$, ${\bf D}elta_1 = \omega_2 - \omega_1$, and ${\bf D}elta_2 = \omega_3 - \omega_2$ and becomes
{\bm \epsilon}gin{align*}
\dot{{\bf v}arphi}_1 &= {\bf D}elta_1 + \sigma[ \sin {\bf v}arphi_2 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_1],
\\
\dot{{\bf v}arphi}_2 &= {\bf D}elta_2 + \sigma[ \sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2].
\end{align*}
Therefore we arrive at the equation
{\bm \epsilon}gin{align*}
&{\bf D}elta_1 + \sigma(s)[ \sin {\bf v}arphi_2(s) - \sin({\bf v}arphi_1(s)+{\bf v}arphi_2(s)) -2\sin {\bf v}arphi_1(s)] = 0,
\\
&{\bf D}elta_2 + \sigma(s)[\sin {\bf v}arphi_1(s) - \sin({\bf v}arphi_1(s)+{\bf v}arphi_2(s)) -2\sin {\bf v}arphi_2(s)] = 0.
\end{align*}
Differentiating with respect to $s$ gives
{\bm \epsilon}gin{align*}
\sigma'
{\bm \epsilon}gin{pmatrix}
\sin {\bf v}arphi_2 - \sin({\bf v}arphi_1 + {\bf v}arphi_2) -2\sin {\bf v}arphi_1
\\
\sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2
\end{pmatrix}
+ \sigma
{\bm \epsilon}gin{pmatrix}
-\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_1 & \cos {\bf v}arphi_2 - \cos({\bf v}arphi_1 + {\bf v}arphi_2)
\\
\cos {\bf v}arphi_1 - \cos({\bf v}arphi_1 + {\bf v}arphi_2) & -\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_2
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
\dot{{\bf v}arphi}_1
\\
\dot{{\bf v}arphi}_2
\end{pmatrix}
=
{\bm \epsilon}gin{pmatrix}
0
\\
0
\end{pmatrix}.
\end{align*}
Now the determinant of the $2 \times 2$ matrix is
{\bm \epsilon}gin{align*}
\det = 3[ \cos {\bf v}arphi_1 \cos {\bf v}arphi_2 + \cos {\bf v}arphi_1 \cos({\bf v}arphi_1 + {\bf v}arphi_2) + \cos {\bf v}arphi_2 \cos({\bf v}arphi_1 + {\bf v}arphi_2)].
\end{align*}
Now suppose that this is not zero. Then
{\bm \epsilon}gin{align*}
{\bm \epsilon}gin{pmatrix}
\dot{{\bf v}arphi}_1
\\
\dot{{\bf v}arphi}_2
\end{pmatrix}
&=
\frac{\sigma'}{\sigma}
{\bm \epsilon}gin{pmatrix}
-\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_1 & \cos {\bf v}arphi_2 - \cos({\bf v}arphi_1 + {\bf v}arphi_2)
\\
\cos {\bf v}arphi_1 - \cos({\bf v}arphi_1 + {\bf v}arphi_2) & -\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_2
\end{pmatrix}^{-1}
{\bm \epsilon}gin{pmatrix}
\sin {\bf v}arphi_2 - \sin({\bf v}arphi_1 + {\bf v}arphi_2) -2\sin {\bf v}arphi_1
\\
\sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2
\end{pmatrix}
\\
&= \frac{\sigma'}{\sigma} \frac{1}{\det}
{\bm \epsilon}gin{pmatrix}
-\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_2 & -\cos {\bf v}arphi_2 + \cos({\bf v}arphi_1 + {\bf v}arphi_2)
\\
-\cos {\bf v}arphi_1 + \cos({\bf v}arphi_1 + {\bf v}arphi_2) & -\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_1
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
\sin {\bf v}arphi_2 - \sin({\bf v}arphi_1 + {\bf v}arphi_2) -2\sin {\bf v}arphi_1
\\
\sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2
\end{pmatrix}.
\end{align*}
This system of course has the same orbits as the system
{\bm \epsilon}gin{align*}
{\bm \epsilon}gin{pmatrix}
\dot{{\bf v}arphi}_1
\\
\dot{{\bf v}arphi}_2
\end{pmatrix}
=
\pm \frac{\sigma'}{\sigma}
{\bm \epsilon}gin{pmatrix}
-\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_2 & -\cos {\bf v}arphi_2 + \cos({\bf v}arphi_1 + {\bf v}arphi_2)
\\
-\cos {\bf v}arphi_1 + \cos({\bf v}arphi_1 + {\bf v}arphi_2) & -\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_1
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
\sin {\bf v}arphi_2 - \sin({\bf v}arphi_1 + {\bf v}arphi_2) -2\sin {\bf v}arphi_1
\\
\sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2
\end{pmatrix}.
\end{align*}
Therefore if I consider the system
{\bm \epsilon}gin{align*}
{\bm \epsilon}gin{pmatrix}
\dot{{\bf v}arphi}_1
\\
\dot{{\bf v}arphi}_2
\end{pmatrix}
=
{\bm \epsilon}gin{pmatrix}
-\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_2 & -\cos {\bf v}arphi_2 + \cos({\bf v}arphi_1 + {\bf v}arphi_2)
\\
-\cos {\bf v}arphi_1 + \cos({\bf v}arphi_1 + {\bf v}arphi_2) & -\cos ({\bf v}arphi_1 + {\bf v}arphi_2) - 2\cos {\bf v}arphi_1
\end{pmatrix}
{\bm \epsilon}gin{pmatrix}
\sin {\bf v}arphi_2 - \sin({\bf v}arphi_1 + {\bf v}arphi_2) -2\sin {\bf v}arphi_1
\\
\sin {\bf v}arphi_1 - \sin({\bf v}arphi_1+{\bf v}arphi_2) -2\sin {\bf v}arphi_2
\end{pmatrix}.
\end{align*}
This corresponds to
{\bm \epsilon}gin{align*}
\frac{\sigma'}{\sigma} =
{\bm \epsilon}gin{cases}
1 & \mbox{if $\det > 0$,}
\\
-1 & \mbox{if $\det < 0$.}
\end{cases}
\end{align*}
Define
{\bm \epsilon}gin{align*}
\sigma(s) = \frac{\| {\bf g}({\bm \phi}(s)) \|}{\| {\bf g}({\bm \theta}) \|}
\end{align*}
this is real analytic at any $s$ for which ${\bf g}({\bm \phi}(s)) \ne 0$. This is however always the case as $... {\bf g}e \| {\bm \omega} \| / \max ... > 0$.

In this case ... becomes
{\bm \epsilon}gin{align*}
\omega_1 - \sigma(s)[ \sin(\phi_1(s) - \phi_2(s)) + \sin( \phi_1(s) - \phi_3(s)] &= 0,
\\
\omega_2 - \sigma(s)[ \sin(\phi_2(s) - \phi_1(s)) + \sin( \phi_2(s) - \phi_3(s)] &= 0,
\\
\omega_3 - \sigma(s)[ \sin(\phi_3(s) - \phi_1(s)) + \sin( \phi_3(s) - \phi_2(s)] &= 0.
\end{align*}
If we differentiate this system of equations with respect to $s$ we get that
{\bm \epsilon}gin{align*}
\sigma'(s) ... + \sigma(s) .... . ... = ...
\end{align*}
\end{example}
\fi

\iffalse
\section{Introduction}

The Kuramoto model is a system of non-linear differential equations which is used to study many biological and physical phenomena. In particular, it is often used to model systems with periodic components such as coupled oscillator networks. For such systems a fundamental question is whether or not the system phase-locks, namely, does the system approach a state where all of the individual oscillators rotate at a common frequency with fixed angle differences. This phenomenon is also often referred to as synchronization and has many applications, for example, coupled laser arrays ... cells in the heart ... and even swarms of fireflies ...

If there are $N$ ocillators in the network, then the Kuramoto model is given by the system of differential equations
{\bm \epsilon}gin{align} \label{eq:kuramoto}
\frac{d\theta_i}{dt} = \omega_i - \sum_{\{i,j\} \in E} {\bf g}amma_{ij} \sin(\theta_i - \theta_j) \quad \text{for} \quad i \in \{1,\dots,N\}
\end{align}
where $\theta_i$ and $\omega_i$ is the angle and natural frequency of the $i$th oscillator. Furthermore, $E$ enumerates the links between coupled oscillators with $\{i,j\} \in E$ if the $i$th and $j$th oscillators are coupled together. We also let ${\bf g}amma_{ij} = {\bf g}amma_{ji} > 0$ denote the strength of this coupling. Note that the assumption that ${\bf g}amma_{ij} = {\bf g}amma_{ji}$ physically means that the $i$th oscillator affects the $j$th oscillator in the same way that the $j$th oscillator affects the $i$th oscillator. Also, the assumption that the couplings are positive means that the coupling encourages the two oscillators to behave in a similar way. For the remainder of this paper we will only consider ring networks which can be represented by $E = \{ \{1,2\} , \{2,3\}, \dots, \{N-1,N\},\{N,1\} \}$ with $N {\bf g}e 3$.

where the gradient $\nabla$ is with respect to ${\bm \theta}$. First, we note that $\nabla \det_{\text{Red}}({\bf J}({\bm \theta}_0))$ is the normal vector to the level set $\det_{\text{Red}}({\bf J}({\bm \theta}_0)) = 0$ pointing into the region with $\det_{\text{Red}}({\bf J}({\bm \theta}_0)) > 0$. Therefore we see that $\det_{\text{Red}}({\bf J}({\bm \theta}(t)))$ is positive, zero, and negative for small positive $t$, $t =0$, and small negative $t$ respectively. In particular, the $n_-$ must change at $t = 0$. Furthermore, by Lemma \ref{lem:flow} we know that $\| {\bf g}({\bm \theta}(t)) \|$ increases and decreases for small positive and negative $t$ respectively. Graphically, the plot of $\| {\bf g}({\bm \theta}(t)) \|$ near $t = 0$ is convex. Therefore if $n_-$ decreases as $t$ crosses zero from left to right, then ${\bf g}({\bm \theta}(t))$ is retracting linearly towards the origin for small negative $t$ approaching zero from the left. This then results in a decrease of stability at $t = 0$. Similarly, if $n_-$ increases as $t$ crosses zero from left to right, then ${\bf g}({\bm \theta}(t))$ is retracting linearly towards the origin for small positive $t$ approaching zero from the right. Again this results in a decrease of stability at $t = 0$. This setup is illustrated in Figures \ref{fig:generalexample} and \ref{fig:realdeal}   and motivates the following definition.

For our second criterion we will decouple the pair $({\bm \theta},{\bm \gamma})$ and derive a criterion in terms of just ${\bm \theta}$. To do this we recall that ${\bf h} : \mathbb{R}_{>0}^N \rightarrow \mathbb{R}_{>0}^N$ is a bijection by Lemma \ref{lem:h}. Therefore if we make the change of variables ${\bf x} = {\bf h}({\bm \gamma})$ we find that
{\bm \epsilon}gin{gather}
N {\bf h}(\cos {\bm \eta})^\top {\bf x} = \det_{\text{Red}}({\bf J}({\bm \theta})) = 0, \label{eq:criteriax1}
\\
{\bf x}^\top {\bf A}({\bm \eta}) ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}) {\bf x} = {\bf h}({\bm \gamma})^\top {\bf A}({\bm \eta}) ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}) {\bf h}({\bm \gamma}) > 0, \label{eq:criteriax2}
\end{gather}
where ${\bf x} \in \mathbb{R}_{>0}^N$. Therefore it remains to determine the set of ${\bm \theta}$ for which there exists a vector ${\bf x} \in \mathbb{R}_{>0}^N$ satisfying \eqref{eq:criteriax1} and \eqref{eq:criteriax2}. To do this define the matrix
{\bm \epsilon}gin{align*}
{\bf S}_\tau({\bm \theta}) := N^2 {\bf h}(\cos {\bm \eta}) {\bf h}(\cos {\bm \eta})^\top \tau - \sym ({\bf A}({\bm \eta}) ({\bf B} {\bf B}^\top)^{-1} {\bf B} {\bf C}({\bm \eta}))
\end{align*}
where $\sym$ denotes the symmetric part of the matrix, namely, the average of the matrix and its transpose. Recall that a matrix and its symmetric part represent the same bilinear form. With this definition we have that
{\bm \epsilon}gin{align} \label{eq:decomposition}
{\bf h}({\bm \gamma})^\top {\bf S}_\tau({\bm \eta}) {\bf h}({\bm \gamma}) = N^2 \det_{\text{Red}}({\bf J}({\bm \theta}))^2 \tau - \nabla \det_{\text{Red}}({\bf J}({\bm \theta})) {\bf P}^\top \adj({\bf P} {\bf J}({\bm \theta}) {\bf P}^\top) {\bf P} {\bf g}({\bm \theta}).
\end{align}
Since ${\bf h}(\cos {\bm \eta}) {\bf h}(\cos {\bm \eta})^\top$ is rank one, for any principal minor ${\bf M}$ of ${\bf S}_\tau({\bm \eta})$ we can write $\det({\bf M}) = d \tau + e$ and ${\bf M}_{ij} = d_{ij} \tau + e_{ij}$ where both $d$ and $e$ and $d_{ij}$ and $e_{ij}$ are independent of $\tau$ and only depend on ${\bm \eta}$ and the shape of the principal minor ${\bf M}$.

Note that the condition that ${\bf h}(\cos {\bm \eta})$ has non-zero components is a measure zero condition and furthermore isn't need in the case $D < 0$. Also, the zero set of $D$ is contained in the union of zero sets of $d$ as ${\bf M}$ ranges over all principal minors. Furthermore each $d$ is clearly analytic and therefore if each $d$ is not identically zero we have that the zero set of $D$ is a set of measure zero. In Lemma \ref{lem:nonempty} we show $d$ has positive values and hence is not identically zero for all $2 \times 2$ principal minors.

Note that Theorem \ref{thm:criteria2} gives a condition on ${\bm \theta}$ completely independent of ${\bm \gamma}$. For example, one can easily show that $D({\bf 0}) = -1$ and conclude that for any ring network with positive edge weights it is impossible for $n_-$ to decrease at ${\bm \theta} = {\bf 0}$ after the natural frequencies have been retracted linearly. This of course must also hold in some neighborhood of ${\bm \theta} = {\bf 0}$ since $D$ is continuous. On the other hand, in Lemma \ref{lem:nonempty} we show that there always exists ${\bm \theta}$ for which $D$ is positive and therefore for which there exists choices of ${\bm \gamma}$ such that $n_-$ decreases after the natural frequencies have been retracted linearly. To summarize, there exist some angles for which $n_-$ cannot decrease after the natural frequencies have been retracted linearly for any choice of ${\bm \gamma}$ and there exist some angles for which it does for an appropriate choice of ${\bm \gamma}$. In other words it appears that some angular configurations are inherently more stable than others independent of the couplings. Furthermore this inherent property is determined by the sign of $D$.
\fi

\iffalse

\section{Proof of Main Theorem and Numerical Experiments} \label{sec:main}

In this section we prove Theorem \ref{thm:main} which by Theorems \ref{thm:criteria1} and \ref{thm:criteria2} reduces to showing that $D$ takes positive values for every $N {\bf g}e 3$ which we do in Lemma \ref{lem:nonempty}. We then plot $D$ for $N = 3$ which in particular shows the structure of the set $D > 0$ and illustrates the strategy used in the proof of Lemma \ref{lem:nonempty}. 
We now state Lemma \ref{lem:nonempty}, which as mentioned above, completes the proof of Theorem \ref{thm:main}.

{\bm \epsilon}gin{lemma} \label{lem:nonempty}
For every $N {\bf g}e 3$ the quantity $D$ has positive values. In particular, let ${\bf u} \in \mathbb{R}^N$ be a vector such that ..., then $D({\bm \theta}_\epsilon) > 0$ for all sufficiently small $\epsilon \ne 0$ where ${\bm \theta}_\epsilon = (\pi,0,\dots,0) + \epsilon {\bf u}$.
\end{lemma}

The main idea is that the contribution to $D$ of the upper left $2 \times 2$ principal minor is zero at $(\pi,0,\dots,0)$. We then perturb in the direction of ${\bf u}$ to make the contribution of the principal minor positive. This amounts to showing that the sign of the lowest coefficient in the power series expansion with respect to $\epsilon$ is positive. We leave the details of the proof for the appendix. For $N = 3$, this setup is illustrated in Figure ... .

{\bm \epsilon}gin{figure}[H]
\centering
{\bm \epsilon}gin{tabular}{c}
\includegraphics[width=75mm]{PlotD.pdf}
\end{tabular}
\caption{A plot of $D$ for $N = 3$. The horizontal and vertical axes represent the variables $\theta_1$ and $\theta_2$ respectively while we set $\theta_3 = 0$.}
\label{fig:d} 
\end{figure}

{\bm \epsilon}gin{figure}[H]
\centering
{\bm \epsilon}gin{tabular}{c}
\includegraphics[width=37mm]{Plot3bg.pdf}
\end{tabular}
\caption{The plot shows the set of $\| {\bm \gamma} \|_2 = 1$ in the first octant (unit cube) for which \eqref{eq:criteria1} is positive for ${\bm \theta} = (0.39,-2.36,0.00)$ which satisfies $D > 0$ by Figure \ref{fig:d}.}

\label{fig:bg} 
\end{figure}

{\bm \epsilon}gin{figure}[H]
\centering
{\bm \epsilon}gin{tabular}{ccc}
\includegraphics[width=37mm]{Plot3bt.pdf} & \includegraphics[width=37mm]{Plot3vdetbt.pdf} & \includegraphics[width=37mm]{Plot3normnm.pdf}
\end{tabular}
\caption{The left plot shows the orbit ${\bm \theta}(t)$ which is the unique solution of \eqref{eq:flow} with ${\bm \theta}(0) = {\bm \theta}_0$. In addition we show the flow lines for \eqref{eq:flow} and the level set $\det_{\text{Red}}({\bf J}({\bm \theta})) = 0$. Since all relevant functions are invariant under the translation of ${\bm \theta}$ by a multiple of ${\bf 1}$ we can always set $\theta_3 = 0$ and plot $\theta_1$ and $\theta_2$ along the horizontal and vertical axes respectively. The right plot shows $\| {\bf g}({\bm \theta}(t)) \|$ and $n_-({\bf J}({\bm \theta}(t)))$.}
\label{fig:generalexample} 
\end{figure}

\iffalse
equation for lower bound of angle of incidence

{\bm \epsilon}gin{figure}[H]
\centering
{\bm \epsilon}gin{tabular}{c}
\includegraphics[width=37mm]{Chain.pdf}
\end{tabular}
\caption{A plot of the minimal angle of incidence of the flow ${\bm \theta}(t)$ in .. and the level set $\det_{\text{Red}}({\bf J}({\bm \theta})) = 0$.}
\label{fig:angle} 
\end{figure}
\fi
\fi

\iffalse

\section{Numerical Investigation of Loss of Stability} \label{sec:realdeal}

In the last section we showed that for every ring network with $N {\bf g}e 3$ that there exists $({\bm \theta},{\bm \gamma}) \in \mathbb{R}^N \times \mathbb{R}_{>0}^N$ for which $n_-({\bm \theta}(t))$ decreases. However this encompasses a lot of possible transitions such as a stable fixed point becoming a 1-saddle, a 1-saddle becoming a 2-saddle, etc. In this section we will consider the first case, namely a stable point becoming a 1-saddle. Forunately we can characterize these points quite easily using the characteristic polynomial of ${\bf P} {\bf J}({\bm \theta}) {\bf P}^\top$. It is easy to check that a polynomial has only nonpositive roots if and only if all of its coefficients are nonnegative. Therefore if we want ${\bf P} {\bf J}({\bm \theta}) {\bf P}^\top$ to have only negative eigenvalues except for a simple zero eigenvalue we need the coefficients of the characteristic polynomial of ${\bf P} {\bf J}({\bm \theta}) {\bf P}^\top$ to be all non-negative except for the constant and linear terms which should be zero and positive respectively. Bringing it all together we see that a point $({\bm \theta},{\bm \gamma}) \in \mathbb{R}^N \times \mathbb{R}_{>0}^N$ represents the transition of a stable fixed point to a 1-saddle if \eqref{eq:criteria1} is positive and the coefficients of the characteristic polynomial of ${\bf P} {\bf J}({\bm \theta}) {\bf P}^\top$ are all non-negative except for the constant and linear terms which are zero and positive respectively.

For $3 \le N \le 4$ we searched for such a pair $({\bm \theta},{\bm \gamma})$ but found nothing. We used Mathematica for this search, and our strategy was to start from a random initial point $({\bm \theta}_0,{\bm \gamma}_0)$ and then call FindMaximum in order to maximize \eqref{eq:criteria1} subject to the constraints on the coefficients of the characteristic polynomial of ${\bf P} {\bf J}({\bm \theta}) {\bf P}^\top$ and $0 < {\bf g}amma_i < 1$ starting from this initial point. We iterated this procedure 1000 times but never achieved a positive value for \eqref{eq:criteria1}. This strategy was chosen over the functions Maximize, FindInstance, and NSolve. As a result we conjecture that a stable fixed point cannot lose stability to become a 1-saddle for $3 \le N \le 4$.

However, things changed for $N = 5$. In this case we found such a pair $({\bm \theta},{\bm \gamma})$ and we conjecture that this in fact holds for all $N {\bf g}e 5$.

{\bm \epsilon}gin{figure}[H]
\centering
{\bm \epsilon}gin{tabular}{|c|c|c|}
{\bf h}line
${\bm \theta}_0$ & ${\bm \theta}(t)$ & $\| {\bf g}({\bm \theta}(t)) \|$ and $n_-({\bf J}({\bm \theta}(t)))$
\\
{\bf h}line
\includegraphics[width=37mm]{Plot5btring.pdf} & \includegraphics[width=37mm]{Plot5bt.pdf} & \includegraphics[width=37mm]{Plot5normnm.pdf}
\\
{\bf h}line
\end{tabular}
\caption{For $N = 5$ we set ${\bm \theta}_0 = (3.72, 1.39, 6.69, 5.70, 4.71)$ and ${\bm \gamma} = (0.20, 1.00, 1.00, 1.00, 1.00)$. On the left we plot ${\bm \theta}_0$ on a circle. In the middle we plot ${\bm \theta}(t)$ the unique solution of \eqref{eq:flow} with ${\bm \theta}(0) = {\bm \theta}_0$. On the right we plot both the norm $\| {\bf g}({\bm \theta}(t)) \|$ and $n_-({\bf J}({\bm \theta}(t)))$. It is difficult to see the slope of the plot of $\| {\bf g}({\bm \theta}(t)) \|$ but by Lemma \ref{lem:flow} we know that it increases, decreases and finally increases as $n_-$ equals 4, 3, and finally 4.}
\label{fig:realdeal} 
\end{figure}

We briefly give some modifications of lemmas that would be used if we wanted to remove the condition that ${\bm \gamma} \in \mathbb{R}_{>0}^N$ and replace it with just ${\bm \gamma} \in \mathbb{R}_{\ne 0}^N$.

{\bm \epsilon}gin{lemma} \label{lem:hgeneralize}
If $N$ is even, then the restriction ${\bf h} : \mathbb{R}_{\ne 0}^N \rightarrow \mathbb{R}_{\ne 0}^N$ of the map \eqref{eq:h} is a bijection with inverse given by
{\bm \epsilon}gin{align} \label{eq:inverseh}
{\bf h}^{-1}({\bf y})_i = \frac{1}{y_i} \left( \prod_{j=1}^N y_j \right)^{\frac{1}{N-1}}.
\end{align}
\end{lemma}

If $N$ is odd the image of ${\bf h}$ consists of vectors whose components have positive product. We can restrict our focus to these vectors by choosing the ${\bm \epsilon}$ in Lemma \ref{lem:minall} to also have components with positive product.

{\bm \epsilon}gin{lemma} \label{lem:minall}
For any symmetric matrix ${\bf S} \in \mathbb{R}^{n \times n}$ we have that
{\bm \epsilon}gin{align}
\min_{{\bf x} \in \mathbb{R}^n} \frac{{\bf x}^\top {\bf S} {\bf x}}{\| {\bf x} \|_1^2} = \min_{\{{\bf M}, {\bm \epsilon} : \sum_j \epsilon_i \epsilon_j {\bf M}_{ij} > 0\}} \frac{\det({\bf M})}{\sum_{i,j} \epsilon_i \epsilon_j {\bf M}_{ij}}
\end{align}
where $\| {\bf x} \|_1 := \sum_{i=1}^N |x_i|$, ${\bf M}$ denotes a principal minor of ${\bf A}$ with cofactors ${\bf M}_{ij}$, and ${\bm \epsilon}$ denotes a vector with components $\pm 1$. (Again we are using the convention that a singleton matrix has a single cofactor with value one.)
\end{lemma}
\fi

\iffalse

mention angle of incidence

give characteristic polynomials and value of ...

Discuss circle and possible mechanism ...

\fi

\iffalse
In conclusion, we have investigate the stability of fixed points of the Kuramoto model on ring network as the natural frequencies are brought closer together linearly. We have found that in spite of the general intuition that synchronization/stability is increased as the natural frequencies are brought closer together, that in fact this is not always the case. In particular, we showed that for any $N {\bf g}e 3$ there exists a choice of positive edge weights and a fixed point which decreases in stability (measured by the dimensions of the stable and unstable manifolds) as the natural frequencies are brought closer together. We did this by first deriving a criterion in terms of the angle variables and edge weight variables. We then simplified this criterion by decoupling these parameters and deriving a condition solely in terms of the angles. Finally, we numerically applied our criterion to search for a fixed point which although initially stable, loses its stability to become a 1-saddle as the natural frequencies are brought closer together. For $3 \le N \le 4$ we did not find any such point and therefore conjecture that a stable fixed point cannot become a 1-saddle in this case. However, for $N = 5$ we found such an example and in fact conjecture that there are such examples for all $N {\bf g}e 5$.

\fi

\iffalse
{\bm \epsilon}gin{proof}[Proof of Theorem ...]

Finally, let $\lambda_1(s) = 0$ denote the eigenvalue of ${\bf J}b({\bm \theta}(s))$ corresponding to the eigenvector ${\bf 1}$ and $\lambda_2(s) \le \dots \le \lambda_n(s)$ the remaining eigenvalues of ${\bf J}b({\bm \theta}(s))$. Then by assumption we have that $0 = \lambda_1(0) = \lambda_2(0) < \lambda_3(0) \le \dots \le \lambda_n(0)$. By continuity we know that $0 = \lambda_1(s) < \lambda_3(s) \le \dots \le \lambda_n(s)$ for $s$ in a neighborhood of zero. But $\det_\text{Red}({\bf J}b({\bm \theta}(s))) = \lambda_2(s) \dots \lambda_n(s)$ so that $\lambda_2(s) > 0$ for small $s > 0$ and $\lambda_2(s) < 0$ for small $s < 0$. Therefore the branch of the bifurcation for small $s > 0$ is stable.
\end{proof}
\fi

\end{document}