# Equilibria in Kuramoto oscillator networks:
An algebraic approach

Tung T. Nguyen, Roberto C. Budzinski, Jacqueline Đoàn, Federico W. Pasini,
Ján Mináč, Lyle E. Muller  Brain and Mind Institute and Department of
Mathematics, The University of Western Ontario, London, ON, Canada, N6A 5B7
<EMAIL_ADDRESS>Brain and Mind Institute and Department of Mathematics, The
University of Western Ontario, London, ON, Canada, N6A 5B7<EMAIL_ADDRESS>Department of Mathematics, Western University, London, Ontario, Canada N6A 5B7
<EMAIL_ADDRESS>Brain and Mind Institute and Department of Mathematics, The
University of Western Ontario, London, ON, Canada, N6A 5B7<EMAIL_ADDRESS>Brain and Mind Institute and Department of Mathematics, The University of
Western Ontario, London, ON, Canada, N6A 5B7<EMAIL_ADDRESS>Brain and Mind
Institute and Department of Mathematics, The University of Western Ontario,
London, ON, Canada, N6A 5B7<EMAIL_ADDRESS>
###### Abstract.

Kuramoto networks constitute a paradigmatic model for the investigation of
collective behavior in networked systems. Despite many advances in recent
years, many open questions remain on the solutions for systems composed of
coupled Kuramoto oscillators on complex networks. In this article, we describe
an algebraic method to find equilibrium points for this kind of system without
using standard approximations in the limit of infinite system size or the
continuum limit. To do this, we use a recently introduced algebraic approach
to the Kuramoto dynamics, which results in an explicitly solvable complex-
valued equation that captures the dynamics of the original Kuramoto model.
Using this new approach, we obtain equilibria for both the nonlinear original
Kuramoto and complex-valued systems. We then completely classify all
equilibria in the case of complete graphs originally studied by Kuramoto.
Finally, we go on to study equilibria in networks of coupled oscillators with
phase lag, in generalized circulant networks, multi-layer networks, and also
random networks.

###### Contents

  1. 1 Introduction
  2. 2 Equilibrium points of Kuramoto models
  3. 3 Equilibrium points for circulant networks
    1. 3.1 Complete classification of equilibrium points for complete graphs
  4. 4 Equilibrium points for phase-lag coupled oscillators
  5. 5 Equilibrium points for generalized circulant networks
  6. 6 Equilibrium points for multilayer networks
  7. 7 Equilibrium points for random networks
  8. 8 Discussion and conclusions

##  1\. Introduction

A paradigmatic system for understanding the collective behavior of coupled
nodes is given by the Kuramoto model, which was introduced by Yoshiki Kuramoto
in 1975. This model can be described as a set of coupled phase-oscillators
that interact through a nonlinear function, usually considered to be a sine
function [1]. Studies of this mathematical model have revealed several,
previously unknown and non-trivial phenomena that arise from the dynamics and
pattern of connections in this system, such as phase synchronization, remote
synchronization, cluster synchronization, chimera states, and Bellerophon
states [2, 3, 4, 5, 6, 7, 8]. In this sense, this model has been used as a
theoretical approach for studying collective behavior and emergent phenomena
in different fields, spanning from social interactions to biology and physics
[6, 3, 4, 9, 10].

Many recent mathematical studies have focused on possible solutions of the
Kuramoto model – specifically, the equilibrium points of the system [11]. One
of the possible ways to mathematically analyze coupled Kuramoto oscillators is
by considering a mean-field approach in the thermodynamic limit, i.e.
$N\rightarrow\infty$, where $N$ is the number of oscillators in the system
[12, 13]. In this case, it is possible to define the critical coupling
strength where the system transitions to synchronization [14, 15]. Using this
approach, Medvedev has shown the existence of equilibrium points (twisted
states) for small-world networks, where non-local connections play an
important role [16]. Also, the mean-field approximation is an important
strategy in the analytical treatment of coupled oscillators, where the
dynamics of the entire system is reduced to a few variables [10, 17].
Furthermore, many other approaches have been used in order to explore the role
of topology in the synchronization of Kuramoto oscillators. Recent advances
have shown different synchronization phenomena in complex networks [5, 12, 6].

An important question in this context is the possibility of finding
equilibrium points given a specific network of Kuramoto oscillators [18, 19,
20, 21]. Many authors have explored the role of the connection architecture in
the existence and stability of equilibrium points, which affects the process
of transition to a global phase-synchronized state. Interestingly, Townsend et
al. have recently shown, somewhat counter-intuitively, some examples where the
stability of the synchronized state depends most strongly not on the density
of connections in the network, but on the specific pattern of connections,
highlighting the importance of a particular network’s structure on Kuramoto
dynamics [20]. Even though these results represent a great advance in the
understanding of the collective behavior, they reach an inherent difficulty:
the Kuramoto model is originally described by a nonlinear equation, which
imposes a limitation to the analytical treatment, where methods of
linearization and approximations are usually necessary.

In recent work [22], a novel analytical approach to Kuramoto oscillators has
been introduced, which results in an explicitly solvable complex-valued
equation that captures the dynamics of the original, nonlinear Kuramoto model.
This allows the analytical study of the transient behavior for any individual
realization, i.e. the exact state of the system can be obtained at each time
step and for any particular initial condition. While the original Kuramoto
model is defined in terms of real numbers, the analytical one is defined in
terms of complex numbers, where the argument corresponds to solution of the
original Kuramoto model.

In this work, we now use this analytical approach to study the equilibria in
coupled Kuramoto oscillator networks. Because the complex-valued analytical
form has an explicit solution, we can now approach this problem using very
straightforward techniques, where no approximations are needed. We show the
equilibria of the analytical model correspond to those of the original, real-
valued, and nonlinear KM, using algebraic techniques on the adjacency matrix
of an individual network model. The main ingredient of our analysis is the use
of Euler’s formula and some important properties of the matrix exponential.
The results from this work provide insight into both the structure of
equilibria in the original KM, in addition to providing new insight into the
complex-valued analytical form we have introduced in previous work.

This work is organized as follows. We first show that some eigenvectors of the
network adjacency matrix are equilibrium points for the novel analytical model
(Sec. 2). We next show that the equilibrium points for the analytical model
are also equilibrium points for the original (nonlinear) Kuramoto model.
Utilizing this result, we extend our analysis to the case where the network
can be understood as a circulant graph (Sec. 3), focusing also on the typical
case of all-to-all coupling (complete graph) (Sec. 3.1). With this result
obtained, we then extend our analysis of equilibria to four new cases: (1)
equilibria in networks of Kuramoto oscillators with phase lag, which leads to
an interplay between attractiveness and repulsiveness in the coupling term
(Sec. 4), (2) equilibria in generalized circulant networks (Sec. 5), (3)
equilibria in multi-layer networks (Sec. 6), and (4) a new method to design
equilibria in Erdős-Rényi random graphs by applying specific changes to the
adjacency matrix (Sec. 7). Taken together, these results extend the analysis
of equilibria in networks of Kuramoto oscillators to cases that could not be
considered previously, demonstrating the utility of our complex-valued
analytical approach.

##  2\. Equilibrium points of Kuramoto models

The original Kuramoto model can be defined as the dynamical system governed by
the equation

(2.1)
$\frac{d\theta_{i}}{dt}=\omega_{i}+\epsilon\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i}),$

where $\theta_{i}$ is the phase of the $i^{\mathrm{th}}$ oscillator,
$\omega_{i}$ is its natural frequency, $\epsilon$ is the coupling strength,
$N$ is the number of oscillators in the network and $\bm{A}=(a_{ij})$ defines
the adjacency matrix. We initially consider $a_{ij}=0$ if $i$ and $j$ are
unconnected and $1$ if connected, and later consider weighted adjacency
matrices with real-valued entries. In this paper, we deal with the homogeneous
Kuramoto model, where $\omega_{i}=\omega$ for all $i\in[1,N]$. In this case,
considering a suitable rotating frame of reference, without loss of generality
we can assume $\omega=0$.

Following [22], starting with the original (homogeneous) Kuramoto model, we
can extend this to a complex-valued model, obtaining a new dynamical system,
the analytical Kuramoto model, governed by the following equation:

(2.2)
$\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}\big{[}\sin(\theta_{j}-\theta_{i})-\mathrm{i}\cos(\theta_{j}-\theta_{i})\big{]}\,,$

Note this expression now implies $\theta_{i}\in\mathbb{C}$. As shown in [22],
the above equation is equivalent to

(2.3)
$\frac{d}{dt}\left(e^{\mathrm{i}\bm{\theta}}\right)=\epsilon\bm{A}e^{\mathrm{i}\bm{\theta}},$

where $\bm{A}$ is the adjacency matrix. By letting
$\bm{x}=e^{\mathrm{i}\bm{\theta}}$, we have

(2.4) ${\frac{d\bm{x}}{dt}}=\epsilon\bm{A}\bm{x}.$

Let
$\bm{\theta}=\bm{\theta}_{{\text{re}}}+\mathrm{i}\bm{\theta}_{{\text{im}}}$ be
the decomposition of $\bm{\theta}$ into the real and imaginary parts. Then, we
have

(2.5)
$\bm{x}=e^{\mathrm{i}\bm{\theta}_{{\text{re}}}-\bm{\theta}_{{\text{im}}}}=e^{-\bm{\theta_{{\text{im}}}}}e^{\mathrm{i}\bm{\theta}_{{\text{re}}}}.$

By Eq. (2.5), we can observe that $\bm{\theta_{{\text{re}}}}$ is the argument
of the solution $\bm{x}$. This naturally leads to the following definition.

######  Definition 1.

We say that $\bm{\theta}_{0}\in[-\pi,\pi]^{N}$ is an equilibrium point of the
analytical KM – Eq. (2.2) – if for all time $t\geq 0$,
$\bm{\theta}(t)=\bm{\theta}_{0}$. Equivalently

$\arg(\bm{x(t)})=\arg(e^{\epsilon
t\bm{A}}\bm{x}_{0})=\arg(\bm{x}_{0})=\bm{\theta}_{0},$

where $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$.

We also recall the definition of equilibrium points in the theory of
differential equations (see [23] for further discussions).

######  Definition 2.

Consider the following differential equation

(2.6) $\dot{\bm{\theta}}(t)=f(\bm{\theta}(t)),$

where $\bm{\theta}(t)=(\theta_{1}(t),\theta_{2}(t),\ldots,\theta_{N}(t)).$ We
say that $\bm{\theta}_{0}$ is an equilibrium point for Eq. (2.6) if
$f(\bm{\theta}_{0})=0.$ In other words, $\bm{\theta}(t)=\bm{\theta}_{0}$ is a
solution of the following differential equation with initial condition

$\begin{cases}\dot{\bm{\theta}}(t)=f(\bm{\theta}(t))\\\
\bm{\theta}(0)=\bm{\theta}_{0}.\end{cases}$

In the context of the original Kuramoto model Eq. (2.1) (note that we have set
$\omega=0$), $\bm{\theta}_{0}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is an
equilibrium point if and only if for all $1\leq i\leq N$

$\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i})=0.$

######  Remark 1.

In this article, we concentrate on the classical notion of equilibrium points
as discussed in Definition 2 and Definition 1. In our subsequent work, we will
discuss a more general notion of equilibria in Kuramoto networks, namely those
solutions in which $\theta_{j}(t)-\theta_{i}(t)$ is a constant for all time
$t$.

In the following discussion, we develop an effective method to find
equilibrium points of the KM – Eq. (2.1) – and the analytical KM – Eq. (2.2).

We start with the following lemma.

######  Lemma 1.

Let $\bm{B}$ be a matrix. Suppose $\lambda$ is an eigenvalue of a matrix
$\bm{B}$ and $\bm{v}$ is an associated eigenvector. Then

$e^{\bm{B}}v=e^{\lambda}\bm{v}.$

###### Proof.

By definition we have

$e^{\bm{B}}=\sum_{n=0}^{\infty}\frac{\bm{B}^{n}}{n!}.$

Applying both sides to $\bm{v}$ we have

$e^{\bm{B}}\bm{v}=\sum_{n=0}^{\infty}\frac{\bm{B}^{n}}{n!}\bm{v}=\sum_{n=0}^{\infty}\frac{\lambda^{n}}{n!}\bm{v}=e^{\lambda}\bm{v}.$

∎

A direct consequence of the above lemma is the following.

######  Proposition 1.

Suppose $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of
$\bm{A}$ associated with a real eigenvalue $\lambda$. Then $\bm{\theta}_{0}$
is an equilibrium point of the complex KM in the sense of Definition 1.

###### Proof.

By Lemma 1 for $\bm{B}=\epsilon t\bm{A}$, we have

$\bm{x}(t)=e^{\epsilon t\bm{A}}\bm{x}_{0}=e^{\lambda\epsilon t}\bm{x}_{0}.$

Because $\lambda\in\operatorname{\mathbb{R}}$ we have

$\arg(\bm{x}(t))=\arg(e^{\lambda\epsilon
t}\bm{x}_{0})=\arg(\bm{x}_{0})=\bm{\theta}_{0}.$

∎

Next, we show that the above equilibrium points are also equilibrium points
for the original KM Eq. (2.1).

######  Proposition 2.

Suppose $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of
$\bm{A}=(a_{ij})$ associated with a real eigenvalue $\lambda$. Then
$\bm{\theta}_{0}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is an equilibrium
point of the original KM.

$\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i}).$

###### Proof.

By Definition 2, we need to show that for all $1\leq i\leq N$

$\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i})=0.$

First of all, since $\bm{x}_{0}$ is an eigenvector, we know that for all
$1\leq i\leq N$

$\sum_{j=1}^{N}a_{ij}e^{\mathrm{i}\theta_{j}}=\lambda
e^{\mathrm{i}\theta_{i}}.$

Taking the conjugation of both sides and noting that
$a_{ij}\in\operatorname{\mathbb{R}}$ and
$\lambda\in\operatorname{\mathbb{R}}$, we have

$\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\theta_{j}}=\lambda
e^{-\mathrm{i}\theta_{i}}.$

We recall Euler’s formula

(2.7) $\sin(x)=\frac{e^{ix}-e^{-ix}}{2\mathrm{i}},$

thanks to which

$\sin(\theta_{j}-\theta_{i})=\frac{e^{\mathrm{i}(\theta_{j}-\theta_{i})}-e^{-\mathrm{i}(\theta_{j}-\theta_{i})}}{2\mathrm{i}}.$

Hence

$\displaystyle(2\mathrm{i})\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i})$
$\displaystyle=\sum_{j=1}^{N}a_{ij}\left[e^{\mathrm{i}(\theta_{j}-\theta_{i})}-e^{-\mathrm{i}(\theta_{j}-\theta_{i})}\right]$
$\displaystyle=e^{-\mathrm{i}\theta_{i}}\sum_{j=1}^{N}a_{ij}e^{i\theta_{j}}-e^{i\theta_{i}}\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\theta_{j}}$
$\displaystyle=e^{-i\theta_{i}}\lambda
e^{\mathrm{i}\theta_{i}}-e^{\mathrm{i}\theta_{i}}\lambda
e^{-\mathrm{i}\theta_{i}}$ $\displaystyle=\lambda-\lambda=0.$

This completes the proof. ∎

######  Remark 2.

We state many results assuming that equilibrium points have a zero first
coordinate, that is, the phase of the first oscillator is $0$. Since at
equilibrium points all oscillators have a constant phase, this assumption is
not a real loss of generality, but rather it is a normalization obtained as
part of the choice of the rotating frame of reference. However, we remind the
reader that adding a fixed phase to all coordinates of an equilibrium point
gives another equilibrium point.

##  3\. Equilibrium points for circulant networks

We apply Proposition 1 and Proposition 2 to find equilibrium points of the the
KM when the topological structure of the network is circulant. More precisely,
let $\bm{C}$ be a symmetric circulant matrix with first column
$[c_{0},c_{1},\ldots,c_{n-1}]$. Then the Circulant Diagonalization Theorem
ensures that, defining

(3.1) $\bm{\theta}^{(j)}=\left(0,\frac{2\pi j}{n},\ldots,\frac{2\pi
j(n-1)}{n}\right)^{T}\qquad\mbox{for }j=0,\dots,n-1,$

the vector $\bm{x}_{0}^{(j)}=e^{i\bm{\theta}^{(j)}}$ is an eigenvector of
$\bm{C}$ associated with the (real) eigenvalue

$\lambda_{j}=c_{0}+c_{n-1}\omega_{n}^{j}+c_{n-2}\omega_{n}^{2j}+\dots+c_{1}\omega_{n}^{(n-1)j}.$

See also [20].

By Proposition 1 and Proposition 2, we have

######  Proposition 3.

For $0\leq j\leq n-1$, let

$\bm{\theta}^{(j)}=\left(0,\frac{2\pi j}{n},\ldots,\frac{2\pi
j(n-1)}{n}\right)^{T}.$

Then

  1. (1)

$\bm{\theta}^{(j)}$ is an equilibrium point of the analytical KM – Eq. (2.4).

  2. (2)

$\bm{\theta}^{(j)}$ is an equilibrium point of the KM – Eq. (2.1).

Based on these results, Fig. 1 portrays a graphic representation of an
equilibrium point for a circulant network of Kuramoto oscillators, which is
also called as a “twisted state”. Here, a system with $N=50$ oscillators
connected with a ring network, where each node has $k=10$ connections in both
directions. This solution follows Eq. (3.1) with $j=1$ (a), which is
equivalent to take the $2^{\mathrm{nd}}$ eigenvector of the matrix $\bm{A}$.
This leads to an equilibrium point where the oscillators’ phases are equally
spaced between $0$ and $2\pi$ (represented in color-code). Moreover, we also
consider the case where the equilibrium point follows Eq. (3.1), but with
$j=3$ (b), which is equivalent to take the $4^{\mathrm{th}}$ eigenvector of
$\bm{A}$. In this case, the distribution of the phases changes, but it is
still an equilibrium point. These kind of solution is also known as “twisted
states”.

Figure 1. Graphic representation of an equilibrium point for a circulant
network following Eq. (3.1) with $j=1$ (a) and $N=50$ nodes. This is
equivalent to take the $2^{\mathrm{nd}}$ eigenvector of the matrix $\bm{A}$.
The phase of each oscillator is represented in color-code, which is equally
distributed between $0$ and $2\pi$. We also consider the case with $j=3$,
which is equivalent to take the $4^{\mathrm{th}}$ eigenvector of the matrix
$\bm{A}$. This kind of equilibrium point is also known as a “twisted state”.
Spatiotemporal dynamics for both the original Kuramoto model and the new
analytical one. Here, the initial conditions for panel (c) are given by the
equilibrium point represented in (a), while for panel (d) are given by (b). In
both cases, the difference between the oscillators’ phases remains as time
evolves.

Then, we use the equilibrium points represented in Figs. 1a and 1b as initial
conditions for the temporal evolution of the system. The spatiotemporal
dynamics for both the numerical (original) Kuramoto model and the new
analytical version is represented in (c) and (d), where the oscillators’
phases are represented in color-code, and the initial conditions are given,
respectively, by the solutions depicted in (a) and (b). We observe that the
initial configuration remains as time evolves, which corroborates the solution
as an equilibrium point for both cases. The oscillators’ natural frequency is
given by $\omega=20\pi$, and the coupling strength is given by $\epsilon=1.0$.

Furthermore, one can observe that a difference between the dynamics of these
two cases, where the “shape” of these waves differs as the equilibrium point
changes. In this case, the wave pattern in (c) depicts diagonal structures
with a lower slope (in comparison with the horizontal line) than (d) since the
solutions are given by the $2^{\mathrm{nd}}$ and $4^{\mathrm{th}}$
eigenvectors, respectively.

###  3.1. Complete classification of equilibrium points for complete graphs

A particular case of circulant networks is given by complete networks, in
which every oscillator is connected to any other. In other words, the
topological connection is given by the complete graph $K_{N}$. In this case,
based on the previous results, we can give an exhaustive characterization of
equilibrium points.

Let $\bm{A}_{N}$ be the adjacency matrix of $K_{N}$. We have the following
proposition.

######  Proposition 4.

Let $\bm{\theta}_{0}=(\theta_{1},\ldots,\theta_{N})$ be an initial condition
such that

$\sum_{k=1}^{N}e^{\mathrm{i}\theta_{k}}=0.$

Then $e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of $\bm{A}_{N}$
associated with the eigenvalue $\lambda=-1$. Consequently, $\bm{\theta}_{0}$
is an equilibrium point of the analytical KM – Eq. (2.2) – and the original KM
– Eq. (2.1).

###### Proof.

We have

$(\bm{A}_{N}+\bm{I})\begin{pmatrix}e^{\mathrm{i}\theta_{1}}\\\ \vdots\\\
e^{\mathrm{i}\theta_{N}}\end{pmatrix}=\begin{pmatrix}\sum_{k=1}^{N}e^{\mathrm{i}\theta_{k}}\\\
\vdots\\\
\sum_{k=1}^{N}e^{\mathrm{i}\theta_{k}}\end{pmatrix}=\begin{pmatrix}0\\\
\vdots\\\ 0\end{pmatrix}.$

In other words, $e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of
$\bm{A}_{N}$ associated with the eigenvalue $\lambda=-1$. By Proposition 1 and
Proposition 2, $\bm{\theta}_{0}$ is an equilibrium point of the analytical KM
– Eq. (2.2) – and the original KM – Eq. (2.1) – as claimed. ∎

We can go further to classify all equilibrium points of the classical Kuramoto
model – Eq. (2.1). Suppose $\bm{\theta}_{0}=(\theta_{1},\ldots,\theta_{N})$ is
an equilibrium point of the KM – Eq. (2.1)–not encompassed by Proposition 4,
that is, satisfying

$\sum_{i=1}^{N}e^{\mathrm{i}\theta_{i}}\neq 0.$

By taking the conjugate of both sides, this also implies that

$\sum_{i=1}^{N}e^{-\mathrm{i}\theta_{i}}\neq 0.$

Since $\bm{\theta}_{0}=(\theta_{1},\ldots,\theta_{N})$ is an equilibrium
point, for all $1\leq i\leq N$

(3.2) $\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})=0.$

Using Euler’s formula (2.7), Eq. (3.2) becomes

$\sum_{j=1}^{N}e^{\mathrm{i}(\theta_{j}-\theta_{i})}-\sum_{j=1}^{N}e^{-i(\theta_{j}-\theta_{i})}=0.$

We can rewrite this as

$e^{-\mathrm{i}\theta_{i}}\left(\sum_{j=1}^{N}e^{\mathrm{i}\theta_{j}}\right)-e^{\mathrm{i}\theta_{i}}\left(\sum_{j=1}^{N}e^{-\mathrm{i}\theta_{j}}\right)=0.$

We then see that

$e^{2\mathrm{i}\theta_{i}}=\dfrac{\sum_{j=1}^{N}e^{\mathrm{i}\theta_{j}}}{\sum_{j=1}^{N}e^{-\mathrm{i}\theta_{j}}}.$

As the right hand side does not depend on $i$, we conclude that

$e^{2\mathrm{i}\theta_{i}}=e^{2\mathrm{i}\theta_{j}},\forall 1\leq i,j\leq N.$

Assuming, up to a rotation of the frame of reference, that $\theta_{1}=0$ (see
Remark 2). Then the above equation implies that

$e^{2\mathrm{i}\theta_{j}}=1,\forall 1\leq j\leq N.$

In other words $\theta_{j}\in\\{0,\pi\\}$. It is straightforward to check that
if $\theta_{j}\in\\{0,\pi\\}$ then $\bm{\theta}_{0}$ is indeed an equilibrium
point of Eq. (2.1). Note that the property $\theta_{j}=0$ for all $j$ is a
special case of this condition. In summary, we just proved the following.

######  Proposition 5.

The point $\bm{\theta}_{0}=(\theta_{1},\ldots,\theta_{N})$ is an equilibrium
point of the KM – Eq. (2.1) – if and only if $\bm{\theta}_{0}$ satisfies one
of the following conditions.

  1. (1)

The $\theta_{i}$s differ by integer multiples of $\pi$.

  2. (2)

$\sum\limits_{i=1}^{N}e^{\mathrm{i}\theta_{i}}=0$.

##  4\. Equilibrium points for phase-lag coupled oscillators

In this section, we generalize the results in the previous section to find
equilibrium points for the Kuramoto networks where the oscillators are coupled
with phase-lag.

The original KM can be described, in this case, as:

(4.1)
$\frac{d\theta_{i}}{dt}=\omega_{i}+\epsilon\sum\limits_{j=1}^{N}a_{ij}\sin{(\theta_{j}-\theta_{i}-\phi)},$

where $N$ is the number of oscillators, $\epsilon$ is coupling strength,
$\bm{A}=(a_{ij})$ defines the adjacency matrix, and $\phi$ is the phase-lag.
Using this new parameter $\phi$ we can transition from a purely attractive
coupling ($\phi=0$) to a purely repulsive one ($\phi=\nicefrac{{\pi}}{{2}}$).
Here, following the same idea described in Sec. 2, we assume all intrinsic
frequencies $\omega_{i}$ to be equal and we consider the rotating frame, in
which $\omega_{i}=0\hskip 8.5359pt\left(\forall i\in[1,N]\right)$.

Furthermore, using the same technique of complexification, the new analytical
model can be described as:

(4.2)
$\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}\left(\sin(\theta_{j}-\theta_{i}-\phi)-\mathrm{i}\cos(\theta_{j}-\theta_{i}-\phi)\right).$

By Euler’s formula we have

$i\sin(\theta_{j}-\theta_{i}-\phi)+\cos(\theta_{j}-\theta_{i}-\phi)=e^{\mathrm{i}(\theta_{j}-\theta_{i}-\phi)}.$

The analytical KM becomes

$\mathrm{i}\frac{d\theta_{i}}{dt}=\epsilon
e^{-\mathrm{i}\theta_{i}}\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\phi}e^{\mathrm{i}\theta_{j}}.$

Equivalently

$\mathrm{i}e^{\mathrm{i}\theta_{i}}\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\phi}e^{\mathrm{i}\theta_{j}}.$

Let $x_{i}=e^{\mathrm{i}\theta_{i}}$. Then
$\frac{dx_{i}}{dt}=\mathrm{i}e^{\mathrm{i}\theta_{i}}\frac{d\theta_{i}}{dt}$.
Therefore the above equation becomes

$\frac{dx_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\phi}x_{j}.$

The general solution of this linear ODE is

(4.3) $\bm{x}=e^{t\bm{K}}\bm{x}(0),$

where $\bm{x}=(x_{1},\ldots,x_{N})$ and $\bm{K}=\epsilon
e^{-\mathrm{i}\phi}\bm{A}$.

Finally, we explain how to get a real solution out of the complex solution for
the analytical Kuramoto model described by Eq. (4.2). Let
$\bm{\theta}=\bm{\theta}_{{\text{re}}}+\mathrm{i}\bm{\theta}_{{\text{im}}}$ be
the decomposition of $\bm{\theta}$ into the real and imaginary parts. Then we
have

(4.4)
$\bm{x}=e^{\mathrm{i}\bm{\theta}_{{\text{re}}}-\bm{\theta}_{{\text{im}}}}=e^{-\bm{\theta}_{{\text{im}}}}e^{\mathrm{i}\bm{\theta}_{{\text{re}}}}.$

We see that $\bm{\theta}_{{\text{re}}}$ is thus the argument of the analytical
solution $\bm{x}$. In particular, we can take
$\bm{\theta}_{{\text{re}}}\in[-\pi,\pi]$.

We have the following definition which naturally generalizes Definition 1.

######  Definition 3.

We say that $\bm{\theta}_{0}\in[-\pi,\pi]^{N}$ is an equilibrium point of the
phase-lag analytical KM – Eq. (4.2) if for all times $t\geq 0$

$\arg(\bm{x(t)})=\arg(e^{\epsilon
e^{-\mathrm{i}\phi}\bm{A}t}\bm{x}_{0})=\arg(\bm{x}_{0})=\bm{\theta}_{0},$

where $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$.

We have the following proposition about equilibrium points of the analytical
KM.

######  Proposition 6.

Suppose $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of
$\bm{A}=(a_{ij})$ associated with the eigenvalue $\lambda$. Suppose further
that $\lambda e^{-\mathrm{i}\phi}\in\operatorname{\mathbb{R}}$. Then
$\bm{\theta}_{0}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is an equilibrium
point of the phase lag analytical KM – Eq. (4.2)

$\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}\left(\sin(\theta_{j}-\theta_{i}-\phi)-\mathrm{i}\cos(\theta_{j}-\theta_{i}-\phi)\right).$

###### Proof.

Because $\bm{x}_{0}$ is an eigenvector of $\bm{A}$ associated with $\lambda$,
we have $\bm{A}\bm{x}_{0}=\lambda\bm{x}_{0}$. Therefore, applying Lemma 1 with
$B=t\epsilon e^{-\mathrm{i}\phi}\bm{A}$, we have

$\bm{x}(t)=e^{t\epsilon e^{-\mathrm{i}\phi}\bm{A}}\bm{x}_{0}=e^{t\epsilon
e^{-\mathrm{i}\phi}\lambda}\bm{x}_{0}.$

Taking the argument of both sides we have

$\arg(\bm{x}(t))=\arg(e^{t\epsilon
e^{-\mathrm{i}\phi}{\color[rgb]{0,0,1}\lambda}}\bm{x}_{0})=\text{Im}(t\epsilon
e^{-i\phi}{\color[rgb]{0,0,1}\lambda})\arg(\bm{x}_{0})=\arg(\bm{x}_{0})=\bm{\theta}_{0}.$

Here we use the crucial assumption that $\lambda
e^{-i\phi}\in\operatorname{\mathbb{R}}$. We conclude that $\bm{\theta}_{0}$ is
an equilibrium point of the phase-lag analytical KM in the sense of Definition
3. ∎

Next, we show that the equilibrium points described in Proposition 6 are
equilibrium points of the original phase-lag KM (4.1) as well (note that we
have set $\omega=0$.) This is a direct generalization of Proposition 2.

######  Proposition 7.

Suppose $\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}$ is an eigenvector of
$\bm{A}=(a_{ij})$ associated with the eigenvalue $\lambda$. Suppose further
that $\lambda e^{-\mathrm{i}\phi}\in\operatorname{\mathbb{R}}$ (or
equivalently ${\color[rgb]{0,0,1}-}\phi+\arg(\lambda)$ is an integer multiple
of $\pi$). Then $\bm{\theta}_{0}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is
an equilibrium point of the phase lag KM – Eq. (4.1)

$\frac{d\theta_{i}}{dt}=\epsilon\sum\limits_{j=1}^{N}a_{ij}\sin{(\theta_{j}-\theta_{i}-\phi)}.$

###### Proof.

We need to show that for $1\leq i\leq N$

$\sum\limits_{j=1}^{N}a_{ij}\sin{(\theta_{j}-\theta_{i}-\phi)}=0.$

First of all, since $\bm{x}_{0}$ is an eigenvector, we know that for all
$1\leq i\leq N$

$\sum_{j=1}^{N}a_{ij}e^{\mathrm{i}\theta_{j}}=\lambda
e^{\mathrm{i}\theta_{i}}.$

Taking the conjugation of both sides and noting that
$a_{ij}\in\operatorname{\mathbb{R}}$, we have

$\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\theta_{j}}=\bar{\lambda}e^{-\mathrm{i}\theta_{i}}.$

By Euler’s formula (2.7)

$\sin(\theta_{j}-\theta_{i}-\phi)=\frac{e^{\mathrm{i}(\theta_{j}-\theta_{i}-\phi)}-e^{-\mathrm{i}(\theta_{j}-\theta_{i}-\phi)}}{2\mathrm{i}}.$

Hence

$\displaystyle(2\mathrm{i})\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i}-\phi)$
$\displaystyle=\sum_{j=1}^{N}a_{ij}\left[e^{\mathrm{i}(\theta_{j}-\theta_{i}-\phi)}-e^{-\mathrm{i}(\theta_{j}-\theta_{i}-\phi)}\right]$
$\displaystyle=e^{-\mathrm{i}\theta_{i}}e^{-\mathrm{i}\phi}\sum_{j=1}^{N}a_{ij}e^{i\theta_{j}}-e^{i\theta_{i}}e^{\mathrm{i}\phi}\sum_{j=1}^{N}a_{ij}e^{-\mathrm{i}\theta_{j}}$
$\displaystyle=e^{-i\theta_{i}}e^{-i\phi}\lambda
e^{\mathrm{i}\theta_{i}}-e^{\mathrm{i}\theta_{i}}e^{\mathrm{i}\phi}\bar{\lambda}e^{-\mathrm{i}\theta_{i}}$
$\displaystyle=\lambda e^{-\mathrm{i}\phi}-\bar{\lambda}e^{\mathrm{i}\phi}$
$\displaystyle=\lambda e^{-\mathrm{i}\phi}-\overline{\lambda
e^{-\mathrm{i}\phi}}=0.$

The last inequality comes from the assumption that $\lambda
e^{-i\phi}\in\operatorname{\mathbb{R}}$. This completes the proof. ∎

Note that by the same argument, we have the following slightly more general
statement.

######  Proposition 8.

Let
$\bm{x}_{0}=e^{\mathrm{i}\bm{\theta}_{0}}=(e^{\mathrm{i}\theta_{1}},\ldots,e^{\mathrm{i}\theta_{N}})$.
Suppose that for each $1\leq i\leq N$, there exists a pair
$(\lambda_{i},\phi_{i})$ such that
$\lambda_{i}e^{-\mathrm{i}\phi_{i}}\in\operatorname{\mathbb{R}}$ and

$\sum_{j=1}^{N}a_{ij}e^{\mathrm{i}\theta_{j}}=\lambda_{i}e^{\mathrm{i}\theta_{i}}.$

Then $\bm{\theta}_{0}=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is an
equilibrium point of the phase lag KM – Eq. (4.1)

$\frac{d\theta_{i}}{dt}=\epsilon\sum\limits_{j=1}^{N}a_{ij}\sin{(\theta_{j}-\theta_{i}-\phi_{i})},\forall
1\leq i\leq N$

######  Remark 3.

Proposition 7 is a special case of Proposition 8, i.e the case when
$\lambda_{i}=\lambda$ for all $1\leq i\leq N$.

We then use computational analyses to further illustrate the equilibrium
points for phase-lag Kuramoto oscillators. Figure 2 depicts a graphic
representation of this kind of solution for coupled Kuramoto oscillators.
Here, the network is given with the same configuration than in Fig. 1: $N=50$
oscillators coupled in a ring network with $k=10$, the coupling strength is
given by $\epsilon=1.0$, and the natural frequency is given by $\omega=20\pi$.
Furthermore, the equilibrium points represented here are given by the
$2^{\mathrm{nd}}$ eigenvector of the matrix $\bm{K}$, which has information
about the network topology, as well as, the phase-lag factor $\phi$. Here, we
consider $\phi=1.00$ (a), and $\phi=\pi/2$, where one can observe that the
role of this factor is given by a rotation in the solution in the circle.

Figure 2. Graphic representation of equilibrium points for phase-lag Kuramoto
oscillators. Here, a network with $N=50$ is considered, where the adjacency
matrix is given by ring network, where each node has $k=10$ connections in
both directions. Moreover, the solutions are given by the $2^{\mathrm{nd}}$
eigenvector of the matrix $\bm{A}$ with $\phi=1.00$ (a), and $\phi=\pi/2$. The
effect of the phase-lag is observed in the rotation of the solution.
Spatiotemporal dynamics for this system (phase-lag Kuramoto oscillators) is
represented in (c) and (d), where the initial conditions follow the
equilibrium points represented in (a), and (b), respectively. Both the
original Kuramoto model and analytical one depict the same behavior, where the
dynamics remains the same as time evolves, and the equilibrium points can be
understood as traveling waves.

In addition to that, Figs. 2c and 2d illustrate the spatiotemporal dynamics of
the networks considering the initial conditions given by the solutions
represented in panels (a) and (b). Here, the wave pattern is observed for the
whole analysis, in both the original KM and the analytical model, which
corroborates the findings about these equilibrium points. The only difference
between panels (c) and (d) is the rotation in the solution observed in panels
(a) and (b), therefore leading to a shifting in the oscillators’ phases.

######  Example 1.

Let us consider a network with $4$ nodes and their topological connection is
given by the following adjacency matrix

$\bm{A}=\begin{pmatrix}0&0&1&1\\\ 1&0&0&1\\\ 1&1&0&0\\\ 0&1&1&0\end{pmatrix}.$

This is a circulant matrix which is not symmetric. By the CDT, we can see that
$\bm{v}=(1,i,-1,-i)^{T}$ is an eigenvector of $\bm{A}$ associated with the
eigenvalue $\lambda=-(1+i)$. Note that we have
$\bm{v}=e^{\mathrm{i}\bm{\theta}_{0}}$ with
$\bm{\theta}_{0}=(0,\frac{\pi}{2},\pi,-\frac{\pi}{2})$. Additionally, we have

$\lambda=-(1+i)=\sqrt{2}e^{\mathrm{i}(\frac{-3\pi}{4})}.$

Let $\phi=\frac{\pi}{4}$ then $\lambda
e^{-i\phi}=-\sqrt{2}\in\operatorname{\mathbb{R}}$. By Propositions 7 and 6 we
know that $\bm{\theta}_{0}$ is an equilibrium point of the KM associated with
$\bm{A}$ with phase lag $\phi=\frac{\pi}{4}$. Concretely, the KM is described
by the following system of differential equations

$\begin{cases}\frac{d\theta_{1}}{dt}=\sin(\theta_{3}-\theta_{1}-\frac{\pi}{4})+\sin(\theta_{4}-\theta_{1}-\frac{\pi}{4})\\\
\frac{d\theta_{2}}{dt}=\sin(\theta_{1}-\theta_{2}-\frac{\pi}{4})+\sin(\theta_{4}-\theta_{2}-\frac{\pi}{4})\\\
\frac{d\theta_{3}}{dt}=\sin(\theta_{1}-\theta_{3}-\frac{\pi}{4})+\sin(\theta_{2}-\theta_{3}-\frac{\pi}{4})\\\
\frac{d\theta_{4}}{dt}=\sin(\theta_{2}-\theta_{4}-\frac{\pi}{4})+\sin(\theta_{3}-\theta_{4}-\frac{\pi}{4}).\\\
\end{cases}$

Figure 3 portrays a graphic representation of the equilibrium point for this
example (a). Moreover, panel (b) brings the spatiotemporal dynamics for this
network when the initial condition is given by the equilibrium point
represented in (a). In this case, the dynamics is given by a traveling wave
and corroborates the solution being an equilibrium point.

Figure 3. A network following the description of Example 1 is analyzed, where
panel (a) shows an equilibrium point for this system given by
$\theta_{0}=(0,\pi/2,\pi,-\pi/2)$. Panel (b) shows the dynamics of this
network when the equilibrium point represented in (a) is used as initial
condition.

##  5\. Equilibrium points for generalized circulant networks

In this section, we discuss equilibrium points for Kuramoto networks, where
the connection architecture is given by circulant matrices associated with a
group $G$. Here, we use the results obtained in the Secs. 3 and 4 to explore
the features of the equilibrium points. First, we recall their definitions
(see [24] for a more thorough discussion.)

######  Definition 4.

Let $G$ be a finite group. A matrix $\bm{C}$ is called $G$-circulant if it has
the following form

$\bm{C}=(c_{\tau^{-1}\sigma})_{\tau,\sigma\in G},$

where $c_{g}\in\mathbb{C}$ for $g\in G$.

######  Remark 4.

When $G=\operatorname{\mathbb{Z}}/n$, we recover the notion of “circulant
matrices” discussed in Section 3.

The following proposition is proved implicitly in [24, Section 1.2].

######  Proposition 9.

Let $\chi:G\to\mathbb{C}^{\times}$ be a $1$-dimensional representation of $G$.
Then $\bm{v_{\chi}}=(\chi(\sigma))_{\sigma\in G}^{T}$ is an eigenvector of any
$G$-circulant matrix $\bm{C}$. The corresponding eigenvalue is

$Y_{\chi}=\sum_{\sigma\in G}c_{\sigma}\chi(\sigma)$

###### Proof.

For all $\sigma\in G$, the $\tau$-component of the column vector
$\bm{C}\bm{v_{\chi}}$ is given by

$(\bm{C}\bm{v_{\chi}})_{\tau}=\sum_{\sigma\in
G}c_{\tau^{-1}\sigma}\chi(\sigma).$

Letting $g=\tau^{-1}\sigma$, we have $\sigma=\tau g$. Hence

$\displaystyle(\bm{C}\bm{v_{\chi}})_{\tau}$ $\displaystyle=\sum_{g\in
G}c_{g}\chi(\tau g)$ $\displaystyle=\sum_{g\in G}c_{g}\chi(g)\chi(\tau)$
$\displaystyle=\chi(\tau)\sum_{g\in G}c_{g}\chi(g)$
$\displaystyle=Y_{\chi}\chi(\tau)=Y_{\chi}(\bm{v}_{\chi})_{\tau}.$

Here we use the fact that $\chi(\tau g)=\chi(\tau)\chi(g)$ as $\chi$ is a
$1$-dimensional representation of $G$. Since this true for all $\tau\in G$, we
conclude that $\bm{v_{\chi}}$ is an eigenvector of $\bm{C}$ associated with
the eigenvalue $Y_{\chi}$. ∎

Let $n=|G|$. Then by Lagrange’s theorem, $\sigma^{n}=1$ for all $\sigma\in G$.
Therefore

$\chi(\sigma)^{n}=\chi(\sigma^{n})=\chi(1)=1.$

We conclude that $\chi(\sigma)$ must be an $n$-th root of unity. Therefore, we
can define

$\arg(\bm{v_{\chi}})=(\arg(\chi(\sigma))_{\sigma\in
G}\in\frac{2\pi}{n}{\operatorname{\mathbb{Z}}}^{n}.$

By Proposition 7, and Proposition 6, we have the following.

######  Proposition 10.

Suppose the topological connection of oscillators is given by a $G$-circulant
matrix $\bm{A}=(a_{ij})$. Let
$\chi,\bm{v_{\chi}},Y_{\chi},\arg(\bm{v_{\chi}})$ be as above and
$\phi=\arg(Y_{\chi})$. Then

  1. (1)

$\arg(\bm{v_{\chi}})$ is an equilibrium point of the original phase-lag KM

$\frac{d\theta_{i}}{dt}=\epsilon\sum\limits_{j=1}^{N}a_{ij}\sin{(\theta_{j}-\theta_{i}-\phi)}.$

  2. (2)

$\arg(\bm{v_{\chi}})$ is an equilibrium point of the associated phase-lag
analytical KM

$\frac{d\theta_{i}}{dt}=\epsilon\sum_{j=1}^{N}a_{ij}\left(\sin(\theta_{j}-\theta_{i}-\phi)-\mathrm{i}\cos(\theta_{j}-\theta_{i}-\phi)\right).$

##  6\. Equilibrium points for multilayer networks

In this section, we construct some examples of non-circulant network that have
interesting equilibrium points. We do so by applying a recent result on the
join of two circulant graphs [25]. First, we introduce the following
convention. In the following, the operator $\operatorname{\mathbin{\ast}}$
denotes vector concatenation:

$(x_{1},\dots,x_{m})^{T}\operatorname{\mathbin{\ast}}(y_{1},\dots,y_{n})^{T}=(x_{1},\dots,x_{m},y_{1},\dots,y_{n})^{T}$

We start with the following observation (see also [25, Section 2]).

######  Proposition 11.

Let $\bm{A}$ be an $(k_{1}+k_{2})\times(k_{1}+k_{2})$ matrix of the form

$\bm{A}=\begin{pmatrix}\bm{C}&[\alpha]_{k_{1},k_{2}}\\\
[\beta]_{k_{2},k_{1}}&\bm{D}\end{pmatrix}.$

Here $\bm{C}=circ(c_{0},\dots,c_{k_{1}})$ and
$\bm{D}=circ(d_{0},\dots,d_{k_{2}})$ are circulant matrices of respective
sizes $k_{1}\times k_{1}$ and $k_{2}\times k_{2}$, $[\alpha]_{k_{1},k_{2}}$ is
a $k_{1}\times k_{2}$ matrix whose entries are equal to
$\alpha\in\operatorname{\mathbb{R}}$, and $[\beta]_{k_{2},k_{1}}$ is a
$k_{2}\times k_{1}$ matrix whose entries are equal to
$\beta\in\operatorname{\mathbb{R}}$. For $1\leq j\leq k_{1}-1$ let

$\bm{w}_{j}=(1,\omega_{k_{1}}^{j},\omega_{k_{1}}^{2j},\ldots,\omega_{k_{1}}^{(k_{1}-1)j},0\ldots,0)^{T}=\bm{v}_{j,k_{1}}\operatorname{\mathbin{\ast}}\underbrace{(0,0,\ldots,0)^{T}}_{\text{$k_{2}$
zeros}}$

with

$\bm{v}_{j,k_{1}}=(1,\omega_{k_{1}}^{j},\omega_{k_{1}}^{2j},\ldots,\omega_{k_{1}}^{(k_{1}-1)j})^{T}.$

For $1\leq j\leq k_{2}-1$, let

$\bm{z}_{j}=(0\ldots,0,1,\omega_{k_{2}}^{j},\omega_{k_{2}}^{2j},\ldots,\omega_{k_{2}}^{(k_{2}-1)j})^{T}=\underbrace{(0,0,\ldots,0)^{T}}_{\text{$k_{1}$
zeros}}*\bm{v}_{j,k_{2}},$

with

$\bm{v}_{j,k_{2}}=(1,\omega_{k_{2}}^{j},\omega_{k_{2}}^{2j},\ldots,\omega_{k_{2}}^{(k_{2}-1)j})^{T}.$

We have the following

  1. (1)

$\bm{w}_{j}$ is an eigenvector of $\bm{A}$ associated with the eigenvalue

$\lambda_{j}^{\bm{C}}=c_{0}+c_{k_{1}-1}\omega_{k_{1}}^{j}+c_{k_{1}-2}\omega_{k_{1}}^{2j}+\dots+c_{1}\omega_{k_{1}}^{(k_{1}-1)j}$

  2. (2)

Similarly, $\bm{z}_{j}$ is an eigenvector associated with the eigenvalue

$\lambda_{j}^{\bm{D}}=c_{0}+d_{k_{2}-1}\omega_{k_{2}}^{j}+d_{k_{2}-2}\omega_{k_{2}}^{2j}+\dots+d_{1}\omega_{k_{2}}^{(k_{2}-1)j}$

###### Proof.

By the Circulant Diagonalization Theorem, $\bm{v}_{j,k_{1}}$ is an eigenvector
of $\bm{C}$ with respect to the eigenvalue $\lambda_{C}^{j}$. By definition,
we have

$\bm{A}\bm{w}_{j}=\bm{C}\bm{v}_{j,k_{1}}\operatorname{\mathbin{\ast}}\underbrace{(t_{j},t_{j},\ldots,t_{j})^{T}}_{\text{$n-k$
terms}}=\lambda_{j}^{\bm{C}}\bm{v}_{j,k_{1}}\operatorname{\mathbin{\ast}}\underbrace{(t_{j},t_{j},\ldots,t_{j})^{T}}_{\text{$n-k$
terms}}$

Here

$t_{j}=\beta\sum_{i=0}^{k-1}\omega_{k}^{ij}.$

By the assumption $1\leq j\leq k_{1}-1$, we can see that $t_{j}=0$. Therefore,
we conclude that

$\bm{A}\bm{w}_{j}=\lambda_{j}^{\bm{C}}\bm{w}_{j}.$

This proves the first statement. The second statement can be proved by the
same argument. ∎

Let us consider a special case when $\bm{C}=\bm{D}$ and $\bm{C}$ is a
symmetric circulant matrix (so in particular, $k_{1}=k_{2}=k$). In this case,
$\bm{w}_{j}$ and $\bm{z}_{j}$ are both eigenvectors with respect to the same
eigenvalue $\lambda_{j}^{\bm{C}}$. Note that by assumption $\bm{C}$ is
symmetric, so $\lambda_{\bm{C}}^{j}\in\operatorname{\mathbb{R}}$. Furthermore,
for any $\phi\in[0,2\pi)$, $\bm{w}_{j}+e^{\mathrm{i}\phi}\bm{z}_{j}$ is an
eigenvector with respect to $\lambda_{\bm{C}}^{j}$. We observe that

$\displaystyle\bm{w}_{j}+e^{\mathrm{i}\phi}\bm{z}_{j}$
$\displaystyle=\bm{v}_{j,k}*e^{\mathrm{i}\phi}\bm{v}_{j,k}$
$\displaystyle=(1,\omega_{k}^{j},\omega_{k}^{2j},\ldots,\omega_{k}^{(k-1)j},e^{\mathrm{i}\phi},e^{\mathrm{i}\phi}\omega_{k}^{j},e^{\mathrm{i}\phi}\omega_{k}^{2j},\ldots,e^{\mathrm{i}\phi}\omega_{k}^{(k-1)j})^{T}$
$\displaystyle=(e^{\mathrm{i}\theta_{1}},e^{\mathrm{i}\theta_{2}},\ldots,e^{\mathrm{i}\theta_{k}},e^{\mathrm{i}\theta_{k+1}},e^{\mathrm{i}\theta_{k+2}},\ldots,e^{\mathrm{i}\theta_{2k}})^{T},$

where

$\displaystyle\bm{\theta}_{0}^{(j)}$
$\displaystyle=(\theta_{1},\theta_{2},\ldots,\theta_{k},\theta_{k+1},\theta_{k+2},\ldots,\theta_{2k})^{T}$
$\displaystyle=\left(0,\frac{2\pi
j}{k},\ldots,\frac{2\pi(k-1)j}{k},\phi,\frac{2\pi
j}{k}+\phi,\ldots,\frac{2\pi(k-1)j}{k}+\phi\right)^{T}.$

By this argument, Proposition 2, and Proposition 1, we have the following.

######  Proposition 12.

Let $\bm{C}$ be a symmetric $k\times k$ circulant matrix, and $\alpha,\beta$
be two arbitrary real numbers. Consider the following matrix

$\bm{A}=\begin{pmatrix}\bm{C}&[\alpha]_{k,k}\\\
[\beta]_{k,k}&\bm{C}\end{pmatrix}.$

For each $1\leq j\leq k-1$, and for any $\phi\in[0,2\pi)$, let

$\bm{\theta}_{0}^{(j,\phi)}=\left(0,\frac{2\pi
j}{k},\ldots,\frac{2\pi(k-1)j}{k},\phi,\frac{2\pi
j}{k}+\phi,\ldots,\frac{2\pi(k-1)j}{k}+\phi\right)^{T}.$

  1. (1)

$\bm{\theta}_{0}^{(j,\phi)}$ is an equilibrium point of the analytical KM.

  2. (2)

$\bm{\theta}_{0}^{(j,\phi)}$ is also an equilibrium point of the original KM.

######  Remark 5.

Proposition 12 can be generalized to the case where we join $d$ identical
circulant networks. We refer interested reader to [25] for further details.

Based on the results depicted in this section, we perform computational
analyses for the example demonstrated in proposition 12. Then, Fig. 4(a) shows
the graphic representation of an equilibrium point for a network with $N=50$
nodes following $\bm{\theta}_{0}$ with $p=1$. Furthermore, the phases of the
first $N/2$ nodes is represented in bigger circle, while the phases of the
last $N/2$ nodes is given by the smaller ones. Here, the network is described
by a non-circulant graph represented by matrix $\bm{A}$, where $\bm{C}$
follows a ring network with $k=5$, $\alpha=0.25$ and $\beta=0.75$. A graphic
representation of this non-circulant matrix is given by Fig. 4b.

Figure 4. An equilibrium point for a non-circulant network with $N=50$ is
graphically represented in (a), where the phase of each node is depicted in
color-code. Here, the graph is given by the matrix represented in (b), which
is composed by circulant matrices. The spatiotemporal dynamics for this system
when the equilibrium point (a) is used as initial condition is represented in
panel (c) for both the original and analytical Kuramoto models. The matrix and
the equilibrium point follow the example in proposition 12.

Moreover, Fig. 4c represents the spatiotemporal patterns for both original and
analytical Kuramoto models. Here, the analyses use the equilibrium point
represented in Fig. 4a as initial condition, so we can observe the wave
pattern as time evolves.

##  7\. Equilibrium points for random networks

Our analytical approach allows us to “design” a twisted state (equilibrium
point) on an undirected random network, given by a Erdős-Réyni graph with
$N=100$, and $p=0.25$ (Fig. 5a). We first evaluate the eigenspectrum of this
random matrix numerically, and then decompose the matrix using its
eigenvectors and eigenvalues ($A=VDV^{\prime}$). Eigenvalues were arranged in
ascending order by their real part. We then modified the $2^{\mathrm{nd}}$ and
$3^{\mathrm{rd}}$ last eigenvectors by applying the sine and cosine functions,
respectively, to the phase given by Eq. (3.1). Finally, we set the eigenvalues
associated with these eigenvectors to be equal and scaled appropriately. Then,
using $VDV^{\prime}$, we create a modified matrix $\bm{A}^{\prime}$ (Fig. 5b),
which is used as the adjacency matrix for the simulations with the original
KM. With this modification $\cos(\bm{\theta}_{0})$ and $\sin(\bm{\theta}_{0})$
are eigenvectors of the new weighted adjacency matrix associated to the same
real eigenvalue $\lambda$. Consequently
$e^{i\bm{\theta}_{0}}=\cos(\bm{\theta}_{0})+\mathrm{i}\sin(\bm{\theta}_{0})$
is an eigenvector of the new matrix $\bm{A}^{\prime}$ associated with the
eigenvalue $\lambda$. By Proposition 2, we know that $\bm{\theta}_{0}$ is an
equilibrium point of the new system associated with $\bm{A}^{\prime}$.
Finally, it is important to note that these results were for these two
eigenvalues scaled by approximately one order of magnitude, which demonstrates
that this result does not depend on the modified eigenvalue-eigenvector pairs
dominating the resulting system dynamics.

Figure 5. We modify a random (Erdős-Renyi) matrix in order to produce an
equilibrium point given by a twisted state. Here, the original matrix (a) is
modified (b) due to changes in some of the eigenvectors. The equilibrium point
(c) is then used as initial condition to the simulation. The Kuramoto order
parameter ($R(t)$) as a function of time shows that for the random matrix, the
system reaches a phase synchronized state, while for the modified matrix, the
system stays in a wave (“twisted”) state. The spatiotemporal patterns
corroborate these features (e, and f).

The equilibrium point for these systems is then represented in Fig. 5c, which
is given by Eq. (3.1). Using this phase configuration as initial condition for
the simulation leads the systems to different states: in the case of the
random matrix, the system reaches a phase synchronized state ($R=1$); in the
case of the modified matrix, the system stays in a twisted state, which is a
phase-locked but not phase synchronized state ($R=0$) – see Fig. 5d. These
features can be observed in the spatiotemporal dynamics of these networks,
which are depicted in Figs. 5e and f, respectively.

##  8\. Discussion and conclusions

In this paper, we have analyzed equilibria in Kuramoto systems. To do so, we
have used a complex-valued version of the Kuramoto model [22], which allows us
to further investigate equilibria in the original nonlinear Kuramoto model. In
this context, we have shown that some of the eigenvectors of the adjacency
matrix are equilibrium points for the original Kuramoto model and also for the
complex-valued novel analytical model. These results thus indicate there is a
strong correspondence between the original Kuramoto model, which is given by
nonlinear differential equations, and the complex-valued one, which admits an
exact analytical solution for individual realizations of the system and on
finite graphs. Using this approach, we are able to extend the analysis of
equilibria in Kuramoto networks to new conditions not previously considered.

Based on the general result presented in Sec. 2, where we have shown that we
can find equilibrium points by using the eigenvectors of the adjacency matrix
for Kuramoto systems, we first considered equilibria for the well-known case
of circulant networks. We then moved on to study the case of global, all-to-
all connections (complete graph on $N$ nodes), which is the case first studied
by Kuramoto, where we have completely characterized the equilibrium points.

We then have analyzed equilibrium points for phase-lag oscillators, where the
imaginary part of the eigenvalues becomes important. This result allows us to
understand equilibria in Kuramoto system with an interplay between attractive
and repulsive coupling, due to the phase-lag parameter. We also have analyzed
the case of generalized circulant graphs, where the eigenvectors and
eigenvalues can be also imaginary, therefore affecting the equilibria of this
kind of system.

Furthermore, we have studied equilibria in the case of multilayer networks,
where the system can be analyzed as the join of circulant networks. Based on
[25], we have found equilibrium points for this kind of system, which opens
the possibility of application in several systems for the study of spreading
dynamics, neuroscience, synchronization, technical system, and others [26, 27,
28, 29].

Finally, we have used the ideas developed in this paper to “design” an
equilibrium point in a random network. In this case, we have shown a procedure
to change the adjacency matrix in order to create new equilibria in the
system. This shows the utility of our analytical approach and extends the
study of equilibria to a class of networks beyond the circulant graphs
considered previously. Based on this result, future work can extend these
ideas to the controlling Kuramoto systems, demonstrating the utility of
analyzing the original, nonlinear Kuramoto model through the lens of the
complex-valued approach [22].

Throughout this paper, we have shown a novel approach to investigate
equilibria in Kuramoto networks. Based on our complex-valued model for
Kuramoto oscillators, we can now study, analytically, network of Kuramoto
oscillators under new and more varied conditions. Our study reveals a strong
correspondence with the original, nonlinear, Kuramoto model, where we can now
track, analytically, individual realizations of oscillator networks and push
further investigations of the rich dynamics that this kind of system offers.

## Appendix - Computational analyses

The solution for the original Kuramoto model (KM) is given by the numerical
integration of Eq. (2.1), where we use Euler’s method with time step of
$10^{-4}$. On the other hand, the analytical solution is obtained through the
evaluating of Eq. (4.3) using a mathematical tool called Expokit [30], which
is designed to solve exponential matrix equations. Here, we used a windowed
approach to the propagation of the solution, where the final solution for each
window as the initial condition for the subsequent one. We use this approach
to obtain better computational stability for the solution.

## Acknowledgments

This work was supported by BrainsCAN at Western University through the Canada
First Research Excellence Fund (CFREF), the NSF through a NeuroNex award
(#2015276), and the Natural Sciences and Engineering Research Council of
Canada (NSERC) grant R0370A01. J.M. gratefully acknowledges the Western
University Faculty of Science Distinguished Professorship in 2020-2021.

## References

  * [1] Y. Kuramoto.  Chemical oscillations, waves, and turbulence.  Springer Science & Business Media, 2012.
  * [2] D. M. Abrams and S. H. Strogatz.  Chimera states for coupled oscillators.  Physical Review Letters, 93(17):174102, 2004.
  * [3] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou.  Synchronization in complex networks.  Physics Reports, 469(3):93–153, 2008.
  * [4] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, and C. S. Zhou.  The synchronization of chaotic systems.  Physics Reports, 366(1-2):1–101, 2002.
  * [5] F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti, and M. Perc.  Chimeras.  Physics Reports, 2020.
  * [6] F. A. Rodrigues, T. K. D. M. Peron, P. Ji, and J. Kurths.  The kuramoto model in complex networks.  Physics Reports, 610:1–98, 2016.
  * [7] S. H. Strogatz.  From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1-4):1–20, 2000.
  * [8] C. Xu, S. Boccaletti, S. Guan, and Z. Zheng.  Origin of bellerophon states in globally coupled phase oscillators.  Physical Review E, 98(5):050202, 2018.
  * [9] S. H. Strogatz.  Exploring complex networks.  Nature, 410(6825):268–276, 2001.
  * [10] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens.  Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review.  The Journal of Mathematical Neuroscience, 10:1–43, 2020.
  * [11] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler.  The kuramoto model: A simple paradigm for synchronization phenomena.  Reviews of Modern Physics, 77(1):137, 2005.
  * [12] A. Pikovsky and M. Rosenblum.  Dynamics of globally coupled oscillators: Progress and perspectives.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(9):097616, 2015.
  * [13] A. Pikovsky and M. Rosenblum.  Partially integrable dynamics of hierarchical populations of coupled oscillators.  Physical Review Letters, 101(26):264103, 2008.
  * [14] L. Basnarkov and V. Urumov.  Phase transitions in the kuramoto model.  Physical Review E, 76(5):057201, 2007.
  * [15] D. Pazó.  Thermodynamic limit of the first-order phase transition in the kuramoto model.  Physical Review E, 72(4):046211, 2005.
  * [16] G. S. Medvedev.  Small-world networks of kuramoto oscillators.  Physica D: Nonlinear Phenomena, 266:13–22, 2014.
  * [17] X. Hu, S. Boccaletti, W. Huang, X. Zhang, Z. Liu, S. Guan, and C-H Lai.  Exact solution for first-order synchronization transition in a generalized kuramoto model.  Scientific Reports, 4(1):1–6, 2014.
  * [18] J. Lu and S. Steinerberger.  Synchronization of kuramoto oscillators in dense networks.  Nonlinearity, 33(11):5905, 2020.
  * [19] R. Taylor.  There is no non-zero stable fixed point for dense networks in the homogeneous kuramoto model.  Journal of Physics A: Mathematical and Theoretical, 45(5):055102, 2012.
  * [20] A. Townsend, M. Stillman, and S. H. Strogatz.  Dense networks that do not synchronize and sparse ones that do.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8):083142, 2020.
  * [21] R. Yoneda, T. Tatsukawa, and J. Teramae.  The lower bound of the network connectivity guaranteeing in-phase synchronization.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(6):063124, 2021.
  * [22] L. Muller, J. Mináč, and T. T. Nguyen.  Algebraic approach to the kuramoto model.  Physical Reiew. E, 104:L022201, Aug 2021.
  * [23] L. Perko.  Differential Equations and Dynamical Systems.  Springer-Verlag, New York, 2001.
  * [24] S. Kanemitsu and M. Waldschmidt.  Matrices of finite abelian groups, finite fourier transform and codes.  Proc. 6th China-Japan Sem. Number Theory, World Sci. London-Singapore-New Jersey, pages 90–106, 2013.
  * [25] J. Đoàn, J. Mináč, L. Muller, T. T. Nguyen, and F. Pasini.  Join of circulant matrices.  submitted, 2021.
  * [26] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin.  The structure and dynamics of multilayer networks.  Physics Reports, 544(1):1–122, 2014.
  * [27] D. S. Bassett and O. Sporns.  Network neuroscience.  Nature neuroscience, 20(3):353–364, 2017.
  * [28] M. De Domenico, C. Granell, M. A. Porter, and A. Arenas.  The physics of spreading processes in multilayer networks.  Nature Physics, 12(10):901–906, 2016.
  * [29] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter.  Multilayer networks.  Journal of Complex Networks, 2(3):203–271, 2014.
  * [30] R. B. Sidje.  Expokit: A software package for computing matrix exponentials.  ACM Transactions on Mathematical Software (TOMS), 24(1):130–156, 1998.