# Broadcasting solutions on multilayer networks of phase oscillators

Tung T. Nguyen<EMAIL_ADDRESS>Roberto C. Budzinski<EMAIL_ADDRESS>Federico W. Pasini  Department of Mathematics, Western University, London, ON,
Canada  Brain and Mind Institute, Western University, London, ON, Canada
Western Academy for Advanced Research, Western University, London, ON, Canada
Robin Delabays  Center for Control, Dynamical Systems, and Computation,
University of California, Santa Barbara, CA, USA  Ján Mináč  Department of
Mathematics, Western University, London, ON, Canada  Western Academy for
Advanced Research, Western University, London, ON, Canada  Lyle E. Muller
Department of Mathematics, Western University, London, ON, Canada  Brain and
Mind Institute, Western University, London, ON, Canada  Western Academy for
Advanced Research, Western University, London, ON, Canada

###### Abstract

Multilayer networks have been used to model and investigate the dynamical
behavior of a variety of systems. On one hand, this approach offers a more
realistic and rich modeling option than single-layer networks. On the other
hand, it leads to an intrinsic difficulty in analyzing the system. Here, we
introduce an approach to investigate the dynamics of Kuramoto oscillators on
multilayer networks by considering the join of matrices, which represent the
intra-layer and inter-layer connections. This approach provides a “reduced”
representation of the original, multilayered system, where both systems have
an equivalent dynamics. Then, we can find solutions for the reduced system and
broadcast them to the multilayer network. Moreover, using the same idea we can
investigate the stability of these states, where we can obtain information on
the Jacobian of the multilayer system by analyzing the reduced one. This
approach is general for arbitrary connection scheme between nodes within the
same layer. Finally, our work opens the possibility of studying the dynamics
of multilayer networks using a simpler representation, then leading to a
better understanding of dynamical behavior of these systems.

##  I Introduction

Networked systems have been used to model and understand collective behaviors
of several systems [1, 2, 3]. Examples can be found spanning from neuroscience
[4] to ecology [5] to physical systems, like power-grids [6, 7, 8]. This
approach can be understood as single, undirected small networks, or network of
networks, or even networks with higher-order interactions [9, 10]. In this
context, one important approach to improve the modeling of real systems is
given by multilayer networks, where two levels of connections are considering:
connections within each layer (intra-layer), and connections between layers
(inter-layer) [11, 12]. This approach allows the study of important systems
such as broadcasting processes [13], animal behavior [14], neural dynamics
[15, 4], and synchronization phenomena [16, 17]. On the other hand, the
addition of an extra level of connections makes the analysis process more
complicated, especially for analytical and mechanistic insights. In this
paper, we introduce an approach that offers a simplified way to analyze the
dynamics of multilayer networks, which allows us to obtain equilibrium points
and analyze the stabiliy for a large class of such systems.

Based on the framework developed in [18, 19, 20], we can represent the
multilayer system as a join of matrices containing information about the
intra-layer and inter-layer connnections. We then introduce a reduced
representation of the system, which has a clear connection with the join
matrix that represents the multilayer network. Thus, this reduced system
allows us to study the dynamics of the mutilayer system in a simpler way.
Then, we can find solutions for the reduced system and broadcast them to the
multilayer networks. This construction offers a direct way to find equilibrium
points in multilayer networks with an arbitrary number of layers and arbitrary
connection scheme for nodes within the same layer (intra-layer connections).

Here, the individual dynamics of each node is modelled by the Kuramoto model
[21, 22, 23], so each node is represented by a phase oscillator. The Kuramoto
model is a central model for the studying of synchronization and the
collective behavior of several systems. It has been used in problems spanning
from social interactions [24] to biology [25, 26] and physics [27, 28].
Particularly, a variety of synchronization phenomena and spatiotemporal
patterns has been observed in Kuramoto systems: phase synchronization [23,
29], remote synchronization [30], twisted states [31, 32], chimera states [33,
34], and Bellerophon states [35].

Furthermore, by considering the properties of joins of matrices and the idea
of representing the multilayer network by the reduced system, we can
investigate into the stability of equilibrium points for Kuramoto oscillators
on multilayer networks. When the solutions for the multilayer systems are
obtained through the broadcasting from the reduced system, we have found a
relation between the Jacobian of both systems. In this sense, we can evaluate
the spectrum of a simpler matrix - representing the reduced systems - and
obtain information on the stability of the equilibrium points for the
multilayer network.

In this paper, we first introduce Kuramoto oscillators on multilayer networks
and expose our approach with the reduced system and also show how we broadcast
solutions to the multilayer system - Sec. II; we extend the approach to the
linear stability analysis of multilayer networks - in Sec. III; in addition to
this, we show numerical simulations and examples of Kuramoto multilayer
networks and their reduced version, which corroborate the main approach
introduced in this paper - Sec. IV; finally, we discuss the results presented
here and their implications, and draw our conclusions - Sec. V.

##  II Kuramoto oscillators on multilayer networks

The Kuramoto model [21, 22, 23] can be defined as the dynamical system
governed by the equation

$\frac{d\theta_{i}}{dt}=\omega_{i}+\sum_{j=1}^{\mathcal{N}}A_{ij}\sin(\theta_{j}-\theta_{i}),$
(II.1)

where $\theta_{i}\in[-\pi,\,\pi]$ is the phase of the $i^{\mathrm{th}}$
oscillator, $\omega_{i}$ is its natural frequency, $\mathcal{N}$ is the number
of oscillators in the system, and $A_{ij}$ are the elements of the weighted
adjacency matrix that represents the connections of this system. That is,
$A_{ij}=0$ if nodes $i$ and $j$ are not connected, and $A_{ij}>0$ if nodes $i$
and $j$ are connected. Here, we consider the case where all oscillators have
the same natural frequency $\omega_{i}=\omega$ for all
$i\in[1,\,\mathcal{N}]$.

###  II.1 Multilayer Kuramoto models

For multilayer systems, we can consider two levels of interactions: intra-
layer (within each layer) and inter-layer (between layers) connections. In
this case, the Kuramoto model for the $i^{\mathrm{th}}$ oscillator in the
$l^{\mathrm{th}}$ layer can be written as:

$\frac{d(\bm{\theta}_{l})_{i}}{dt}=\omega+\underbrace{\sum_{j=1}^{N}(\bm{A}_{ll})_{ij}\sin((\bm{\theta}_{l})_{j}-(\bm{\theta}_{l})_{i})}_{\mathrm{intra-
layer}}+\underbrace{\sum_{k=1,k\neq
l}^{M}\sum_{j=1}^{N}(\bm{A}_{lk})_{ij}\sin((\bm{\theta}_{k})_{j}-(\bm{\theta}_{l})_{i})}_{\mathrm{inter-
layer}},$ (II.2)

where $i,j\in[1,\,N]$ and $l,k\in[1,\,M]$, $M$ being the number of layers and
$N$ the number of nodes in each one. Here, the coupling strength between two
oscillators can assume different values based on the layers these oscillators
are in. Moreover, it is important to emphasize that the intra-layer coupling
considers the coupling between oscillators within a given layer; the inter-
layer coupling considers the coupling between oscillators in different layers.

In this paper, on the other hand, we consider a different representation of a
multilayer system. Here, we use properties of joins of matrices to describe
Kuramoto oscillators on multilayer networks as:

$\frac{d\theta_{i}}{dt}=\omega+\sum_{j=1}^{NM}A_{ij}\sin(\theta_{j}-\theta_{i}),$
(II.3)

where now $i\in[1,\,NM]$, such that:

$\bm{\theta}=(\underbrace{\theta_{1},\theta_{2},\cdots,\theta_{N}}_{1^{\mathrm{st}}\mathrm{layer}},\underbrace{\theta_{N+1},\theta_{N+2},\cdots,\theta_{2N}}_{2^{\mathrm{nd}}\mathrm{layer}},\cdots,\underbrace{\theta_{N(M-1)+1},\theta_{N(M-1)+2},\cdots,\theta_{NM}}_{M^{\mathrm{th}}\mathrm{layer}}).$
(II.4)

The adjacency matrix describing the multilayer system is the $NM\times NM$
matrix:

$\bm{A}=\left(\begin{array}[]{c|c|c|c}\bm{A}_{11}&\bm{A}_{12}&\cdots&\bm{A}_{1M}\\\
\hline\cr\bm{A}_{21}&\bm{A}_{22}&\cdots&\bm{A}_{2M}\\\
\hline\cr\vdots&\vdots&\ddots&\vdots\\\
\hline\cr\bm{A}_{M1}&\bm{A}_{M2}&\cdots&\bm{A}_{MM}\end{array}\right),$ (II.5)

where $\bm{A}_{lk}$ is a $N\times N$ matrix, which describes the connections
between the oscillators in the $l^{\mathrm{th}}$ and $k^{\mathrm{th}}$ layers.
In case $l=k$, the block represents the intra-layer connections in the
$l^{\mathrm{th}}$ layer. In case $l\neq k$, the block represents the inter-
layer connections between oscillators in layers $l$ and $k$. Note that in our
treatment these matrices absorb the information about the coupling strength.
For simplicity, we consider here the case in which the coupling strength
between oscillators in two fixed layers is constant. However, this coupling
need not be the same for all layers.

For details on the joins of matrices and on the mathematical background of our
approach, see [18, 19, 20].

###  II.2 Reduced Kuramoto models and broadcased solutions

In this paper we focus on networks satisfying a regularity condition: namely
that the off-diagonal blocks $\bm{A}_{lk}$ for $l\neq k$ have the same row sum
(that is, they are row-regular, or semimagic, matrices, [19]). In terms of the
Kuramoto model, this assumption means that, for any pair of layers, each
oscillator in the first layer is connected with the same number of oscillators
in the other layer. To investigate the dynamics of such a system, we introduce
an auxiliary reduced $M\times M$ network, in which each layer of the original
system is condensed into one global oscillator. More concretely, the reduced
system is described by the adjacency matrix

$\bar{\bm{A}}=\begin{pmatrix}r_{\bm{A}_{11}}&r_{\bm{A}_{12}}&\cdots&r_{\bm{A}_{1M}}\\\
r_{\bm{A}_{21}}&r_{\bm{A}_{22}}&\cdots&r_{\bm{A}_{2M}}\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{\bm{A}_{M1}}&r_{\bm{A}_{M2}}&\cdots&r_{\bm{A}_{MM}}\end{pmatrix},$ (II.6)

where $r_{\bm{A}_{lk}}$ is the row sum of $\bm{A}_{lk}$. At this point, it is
important to emphasize that these matrices have information about the coupling
strength for the intra-layer and inter-layer case, thus the line sum captures
this information. We also remark that the row sums $r_{\bm{A}_{lk}}$ capture
the information about the inter-layer and intra-layer coupling strengths; this
information is therefore present in the reduced system.

We can now define the reduced Kuramoto model:

$\frac{d\bar{\theta}_{i}}{dt}=\omega+\sum_{j=1}^{M}\bar{A}_{ij}\sin(\bar{\theta}_{j}-\bar{\theta}_{i}),$
(II.7)

where now $i\in[1,\,M]$. We emphasize that the main diagonal of $\bar{\bm{A}}$
does not affect the dynamics of the oscillators, since
$\sin(\bar{\theta}_{i}-\bar{\theta}_{i})=0$. Therefore, in the following, we
often make the harmless identification

$\bar{\bm{A}}=\begin{pmatrix}r_{\bm{A}_{11}}&r_{\bm{A}_{12}}&\cdots&r_{\bm{A}_{1M}}\\\
r_{\bm{A}_{21}}&r_{\bm{A}_{22}}&\cdots&r_{\bm{A}_{2M}}\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{\bm{A}_{M1}}&r_{\bm{A}_{M2}}&\cdots&r_{\bm{A}_{MM}}\end{pmatrix}\simeq\begin{pmatrix}0&r_{\bm{A}_{12}}&\cdots&r_{\bm{A}_{1M}}\\\
r_{\bm{A}_{21}}&0&\cdots&r_{\bm{A}_{2M}}\\\ \vdots&\vdots&\ddots&\vdots\\\
r_{\bm{A}_{M1}}&r_{\bm{A}_{M2}}&\cdots&0\end{pmatrix}.$ (II.8)

The Kuramoto network on the multilayer system - Eq. (II.3) - represented by
$\bm{A}$ \- Eq. (II.5) - and the Kuramoto network on the reduced system - Eq.
(II.7) - represented by $\bar{\bm{A}}$ \- Eq. (II.6) - have the same dynamics.
Thus, we can use the reduced system to better understand the dynamics on the
more complicated, multilayer network.

More precisely, solutions for the reduced system can be broadcast to the
original multilayer network: given a solution of the reduced Kuramoto model:

$\bar{\bm{\theta}}^{\ast}=(\bar{\theta}_{1}^{\ast},\bar{\theta}_{2}^{\ast},\cdots,\bar{\theta}_{M}^{\ast}),$
(II.9)

then

$\bm{\theta}^{\ast}=(\underbrace{\bar{\theta}_{1}^{\ast},\bar{\theta}_{1}^{\ast},\cdots,\bar{\theta}_{1}^{\ast}}_{1^{\mathrm{st}}\mathrm{layer}},\underbrace{\bar{\theta}_{2}^{\ast},\bar{\theta}_{2}^{\ast},\cdots,\bar{\theta}_{2}^{\ast}}_{2^{\mathrm{nd}}\mathrm{layer}},\cdots,\underbrace{\bar{\theta}_{M}^{\ast},\bar{\theta}_{M}^{\ast},\cdots,\bar{\theta}_{M}^{\ast}}_{M^{\mathrm{th}}\mathrm{layer}}),$
(II.10)

is the corresponding broadcasting solution of the multilayer system described
by $\bm{A}$.

Moreover, we note the dependence of the dynamics on the coupling strength. The
inter-layer coupling controlled by $\epsilon$ is also presented in the reduced
system, see Eq. (II.6). So, if we consider the same initial state, both
systems evolve to the same final state. However, the stronger the coupling,
the faster is the transition, which we summarize in the following proposition
whose proof follows from direct calculations (valid for nonzero coupling
only):

######  Proposition II.1.

Suppose
$\bar{\bm{\theta}}^{\ast}=(\bar{\theta}_{1}^{\ast},\bar{\theta}_{2}^{\ast},\cdots,\bar{\theta}_{M}^{\ast})$
is a solution of the Kuramoto model

$\dot{\theta}_{i}=\epsilon\sum_{j}A_{ij}\sin(\theta_{j}-\theta_{i}),$

with intial condition $\bm{\theta}(0)=\bm{\theta}_{0}.$ Let
$\bm{\psi}(t)=\bar{\bm{\theta}}^{\ast}\big{(}\frac{\epsilon^{\prime}}{\epsilon}t\big{)}$.
Then $\bm{\psi}(t)$ is a solution of the Kuramoto model

$\dot{\psi}_{i}=\epsilon^{\prime}\sum_{j}A_{ij}\sin(\psi_{j}-\psi_{i}),$

with the same initial condition $\bm{\psi}(0)=\bm{\theta}_{0}.$

######  Remark II.2.

We note that the idea presented in this section about the correspondence
between the multilayer network and the reduced system is valid even when the
layers have a different number of nodes. In this case, $\bm{A}$ would be a
$\mathcal{M}\times\mathcal{M}$ matrix, where
$\mathcal{M}=N_{1}+N_{2}+\cdots+N_{M}$. We presented here the version where
the layers have the same number of nodes $N$ for a simplicity purpose.

###  II.3 Kuramoto order parameter

In order to characterize the dynamical behavior of the system, we use the
Kuramoto order parameter [23, 22]. The order parameter for the multilayer
system is defined as

$R(t)=\frac{1}{NM}\left|\sum\limits_{j=1}^{NM}\exp{(\mathrm{i}\theta_{j}(t))}\right|,$
(II.11)

where $\theta(t)$ is given by Eq. (II.3). Here, $R(t)=1$ means that all
oscillators in all layers have the same phase at a given time $t$, which is
defined as phase synchronization. If the oscillators are not phase
synchronized, the order parameter is less than $1$. For the asynchronous
behavior, $R$ assumes residual values. When the system is on a phase-locking
solution, also called a twisted state, $R(t)=0$.

We also measure the level of synchronization of the Kuramoto network on the
reduced representation using the Kuramoto order parameter. In this case, the
order parameter is written as

$\bar{R}(t)=\frac{1}{M}\left|\sum\limits_{j=1}^{M}\exp{(\mathrm{i}\bar{\theta}_{j}(t))}\right|,$
(II.12)

where $\bar{R}(t)=1$ means the reduced system is phase synchronized at time
$t$, an asynchronous behavior leads to a residual value of $\bar{R}$, and if
the reduced system has a phase-locking solution (“twisted state”), then
$\bar{R}=0$.

By the definition of the broadcast solution and direct calculations, we have
the following proposition.

######  Proposition II.3.

Suppose that $\bar{\bm{\theta}}^{\ast}$ is a solution of the Kuramoto model on
$\bar{\bm{A}}$ and $\bm{\theta}^{\ast}$ is the corresponding broadcasting
solution of the Kuramoto model on $\bm{A}.$ Then the two systems have the same
Kuramoto order parameter for all times: $R(t)=\bar{R}(t),\,t\geq 0.$

##  III Stability of equilibrium points

We now propose to analyze the existence and stability of some equilibria for
multilayer networks through our knowledge of the existence and stability
properties of equilibria in the reduced system. Of course, there is loss in
complexity when considering the reduced system and we cannot guarantee that
its analysis provides a full picture of the whole multilayer system. However,
we will show that some crucial information can be obtained from the reduced
system by simply broadcasting a known solution thereof to the multilayer
network.

In this section, we require that connections between each pair of connected
layers are dense, which we approximate by a mean-field (or all-to-all)
coupling. In other words, we assume that if layers $l$ and $k$ are connected,
then each node in layer $l$ is connected to each node in layer $k$. Whereas
these assumptions formally restrict the applicability of our method, they
still cover a rather wide class of networked systems. It is furthermore quite
standard to consider identical coupling strength within a networked system and
to approximate unknown or uncertain couplings by mean-field interactions [36,
37, 38, 39, 40]. Formally, the above assumptions translate as:

  1. 1.

$\bm{A}_{ll}$ is the adjacency matrix of an undirected, connected, weighted
graph with positive weights.

  2. 2.

For $l\neq k$, $\bm{A}_{lk}=\epsilon_{lk}\bm{1}_{N}$ where $\bm{1}_{N}$ is a
matrix $N\times N$ with all entries equal to $1$.

As observed in the previous section, the (simpler) reduced system can be used
to study the dynamics of the multilayer network. In particular, equilibrium
points of the reduced system can be broadcast to equilibrium points of the
multilayer one.

However, for stability analysis, we cannot extend the idea in a
straightforward way. The reason is that a neighborhood of an equilibrium point
for the multilayer network $\bm{\theta}^{\ast}$ is $NM$-dimensional, while a
neighborhood of an equilibrium point for the reduced system
$\bm{\bar{\theta}}^{\ast}$ is $M$-dimensional. Broadcasting the neighborhood
of $\bar{\bm{\theta}}^{\ast}$

$B_{\varepsilon}(\bar{\bm{\theta}}^{\ast})=\left\\{\bar{\bm{\theta}}^{\ast}+\bm{\varepsilon}\leavevmode\nobreak\
|\leavevmode\nobreak\ \|\bm{\varepsilon}\|_{2}<\varepsilon\right\\}\,,$
(III.1)

to the multilayer system yields an $M$-dimensional subset of the
$NM$-dimensional neighborhood $B_{\varepsilon}(\bm{\theta})$.

In order to perform a stability analysis of the broadcasted fixed point, one
needs a more in-depth analysis of the Jacobian matrix. Fortunately, we will
show that such an analysis can still be performed, and leverages our knowledge
of the fixed point for the reduced system.

In order to study on the stability of these solution we can use the Jacobian
of each system. Here we denote the Jacobian for the multilayer system at the
equilibrium point $\bm{\theta}^{\ast}$ as
$J(\bm{\theta}^{\ast})=J_{\bm{A}}(\bm{\theta}^{\ast})$, and the Jacobian for
the reduced system at the equilibrium point $\bar{\bm{\theta}}^{\ast}$ as
$J(\bar{\bm{\theta}}^{\ast})=J_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast})$.

We first recall the matrix $\bar{\bm{A}}$, defined in Eq. (II.6). Then, we
notice that main diagonal of $\bar{\bm{A}}$ does not affect the dynamics of
the system (since $\sin(\bar{\theta}_{i}-\bar{\theta}_{i})=0$), so we can
write the matrix as

$\bar{\bm{A}}=\begin{pmatrix}0&r_{\bm{A}_{12}}&\cdots&r_{\bm{A}_{1M}}\\\
r_{\bm{A}_{21}}&0&\cdots&r_{\bm{A}_{2M}}\\\ \vdots&\vdots&\ddots&\vdots\\\
r_{\bm{A}_{M1}}&r_{\bm{A}_{M2}}&\cdots&0\end{pmatrix},$ (III.2)

where, by definition, $r_{\bm{A}_{ij}}=N\epsilon_{ij}$. Based on this matrix,
we can now write the Jacobian for the reduced system as

$J(\bar{\bm{\theta}}^{\ast})=\begin{pmatrix}-\lambda_{1}&r_{\bm{A}_{12}}\cos(\bar{\theta}^{\ast}_{2}-\bar{\theta}^{\ast}_{1})&\cdots&r_{\bm{A}_{1M}}\cos(\bar{\theta}^{\ast}_{M}-\bar{\theta}^{\ast}_{1})\\\
r_{\bm{A}_{21}}\cos(\bar{\theta}^{\ast}_{1}-\bar{\theta}^{\ast}_{2})&-\lambda_{2}&\cdots&r_{\bm{A}_{2M}}\cos(\bar{\theta}^{\ast}_{M}-\bar{\theta}^{\ast}_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{\bm{A}_{M1}}\cos(\bar{\theta}^{\ast}_{1}-\bar{\theta}^{\ast}_{M})&r_{\bm{A}_{M2}}\cos(\bar{\theta}^{\ast}_{2}-\bar{\theta}^{\ast}_{M})&\cdots&-\lambda_{M}\end{pmatrix},$
(III.3)

where

$\lambda_{i}=\sum_{j=1,j\neq
i}^{M}r_{\bm{A}_{ij}}\cos(\bar{\theta}^{\ast}_{j}-\bar{\theta}^{\ast}_{i}).$
(III.4)

In particular, the above Jacobians are semimagic square matrices with line sum
zero.

In order to compute the Jacobian for the multilayer system
$J(\bm{\theta}^{\ast})$, we first recall the definition of the Laplacian. For
an undirected weighted graph $H$ with adjacency matrix $\bm{B}_{H}$ and degree
matrix $\bm{D}_{H}$, the Laplacian is defined as:

$\bm{L}_{H}=\bm{D}_{H}-\bm{B}_{H}.$ (III.5)

We note that, as long as $H$ is undirected, $\bm{L}_{H}$ is always a symmetric
semi-magic square matrix (with row sum equals to zero). Regarding the spectrum
of $\bm{L}_{H}$, we have the following proposition.

######  Proposition III.1.

Let $\bm{L}_{H}$ be the weighted Laplacian matrix of an undirected, connected,
weighted graph $H$. Then $\bm{L}_{H}$ has a unique zero eigenvalue and all
other eigenvalues are strictly positive.

###### Proof.

By Gershgorin’s Discs Theorem [41, Theroem 6.1.1], all eigenvalues of
$\bm{L}_{H}$ are nonnegative. Note also that the eigenvectors of $\bm{L}_{H}$
form an orthonormal basis of $\mathbb{R}^{n}$, because $\bm{L}_{H}$ is
symmetric. Now assume that $\bm{v}$ is an eigenvector of $\bm{L}_{H}$
associated to the eigenvalue $0$. Then,

$\displaystyle 0$
$\displaystyle=\bm{v}^{\top}\bm{L}_{H}\bm{v}=\sum_{i<j}a_{ij}(v_{i}-v_{j})^{2}\,,$
(III.6)

where $a_{ij}$ is the weight of the edge linking nodes $i$ and $j$. Each term
in the sum above is necessarily zero. Therefore, for any pair of nodes $(i,j)$
that are connected, $v_{i}=v_{j}$. As the graph is connected, the eigenvector
$\bm{v}$ is necessarily a multiple of the vector $(1,\ldots,1)^{\top}$. Hence
$\bm{L}_{H}$ has a unique zero eigenvalue and all other eigenvalues are
strictly positive. ∎

We now compute the Jacobian $J(\bm{\theta}^{\ast})$. By definition, for two
indices $(i_{1},j_{1})$ such that $i_{1}\in V(G_{i}),j_{1}\in V(G_{j})$ where
$i\neq j$ then

$[J(\bm{\theta}^{\ast})]_{i_{1},j_{1}}=A_{i_{1},j_{1}}\cos(\bar{\theta}_{j}^{\ast}-\bar{\theta}_{i}^{\ast})=\epsilon_{ij}\cos(\bar{\theta}_{j}^{\ast}-\bar{\theta}_{i}^{\ast}).$
(III.7)

Here, $G_{i}$ is the graph representing the $i^{\mathrm{th}}$ layer, and
$\epsilon_{ij}$ is the coupling strength between nodes in the layer $i$ and
$j$. We emphasize that for the stability analysis, the coupling between layers
is considered uniform, i.e., oscillators in layer $i$ are connected to
oscillators in layer $j$ with the same coupling strength.

If $i_{1},i_{2}\in V(G_{i})$ and $i_{1}\neq i_{2},$ then

$[J(\bm{\theta}^{\ast})]_{i_{1},i_{2}}=(\bm{A}_{ii})_{i_{1},i_{2}}\cos(\bar{\theta}_{i}^{\ast}-\bar{\theta}_{i}^{\ast})=(\bm{A}_{ii})_{i_{1},i_{2}}.$
(III.8)

Finally, we need to consider the case $i_{1}=i_{2}$ and $i_{1}\in V(G_{i})$.
For this part, we observe that $J(\bm{\theta}{\ast})$ is a semimagic square
matrix with line sum equal to zero. We can then see that

$[J(\bm{\theta}^{\ast})]_{i_{1},i_{1}}=-\lambda_{i}-\deg(i_{1}).$ (III.9)

with $\deg(i_{1})$ denoting the weighed degree of node $i_{1}$.

By combining these facts, we can write the Jacobian for the multilayer system
at the equilibrium point $\bm{\theta}^{\ast}$ as:

$J(\bm{\theta}^{\ast})=J_{\bm{A}}(\bm{\theta}^{\ast})=\left(\begin{array}[]{c|c|c|c}-\lambda_{1}\bm{I}-\bm{L}_{\bm{A}_{11}}&\epsilon_{12}\cos(\bar{\theta}_{2}^{\ast}-\bar{\theta}_{1}^{\ast}))\bm{1}&\cdots&\epsilon_{1M}\cos(\bar{\theta}_{M}^{\ast}-\bar{\theta}_{1}^{\ast})\bm{1}\\\
\hline\cr\epsilon_{21}\cos(\bar{\theta}_{1}^{\ast}-\bar{\theta}_{2}^{\ast})\bm{1}&-\lambda_{2}\bm{I}-\bm{L}_{\bm{A}_{22}}&\cdots&\epsilon_{2M}\cos(\bar{\theta}_{M}^{\ast}-\bar{\theta}_{2}^{\ast})\bm{1}\\\
\hline\cr\vdots&\vdots&\ddots&\vdots\\\
\hline\cr\epsilon_{M1}\cos(\bar{\theta}_{1}^{\ast}-\bar{\theta}_{M}^{\ast})\bm{1}&\epsilon_{M2}\cos(\bar{\theta}_{2}^{\ast}-\bar{\theta}_{M}^{\ast})\bm{1}&\cdots&-\lambda_{M}\bm{I}-\bm{L}_{\bm{A}_{MM}}\end{array}\right).$
(III.10)

where, we recall that $\epsilon_{ij}$ is the coupling between oscillators in
layer $i$ and $j$.

One realizes, in particular, that $J(\bm{\theta}^{\ast})$ is a joined union of
semimagic square matrices (see [18, 19] for details). Furthermore, we have the
following equality

$\overline{J_{\bm{A}}(\bm{\theta}^{\ast})}=J_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast}).$
(III.11)

It is important to emphasize that Eq. (III.11) is valid if the equilibrium
point $\bm{\theta}^{\ast}$ of the multilayer network is obtained through the
broadcasting of the solution $\bar{\bm{\theta}}^{\ast}$ of the reduced system
- see Eqs. (II.9) and (II.10).

We summarize these ideas in the following proposition:

######  Proposition III.2.

Suppose $\bar{\bm{\theta}}^{\ast}$ is an equilibrium point of the Kuramoto
model associated with the reduced system ($\bar{\bm{A}}$). Then
$\bm{\theta}^{\ast}$ is an equilibrium point of the Kuramoto model on the
multilayer system ($\bm{A}$) following Eq. (II.10). In addition to that,
suppose that the matrices representing the inter-layer connections
($\bm{A}_{ij},\,i\neq j$) have all entries the same. The Jacobian
$J_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast})$ and
$J_{\bm{A}}(\bm{\theta}^{\ast})$ are given by Eq. (III.3) and Eq. (III.10),
respectively. Furthermore, we have the following relation

$\overline{J_{\bm{A}}(\bm{\theta}^{\ast})}=J_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast}).$

For the next proposition, we need to introduce a couple of notations. For a
set $S=\\{s_{1},s_{2},\ldots,s_{m}\\}\subset\operatorname{\mathbb{C}}^{m}$ and
$a,b\in\operatorname{\mathbb{C}}$, we denote

$aS-b=\\{as_{1}-b,as_{2}-b,\ldots,as_{m}-b\\}.$

We also denote by $S\setminus s$ the set difference of $S$ by $s.$ By [19,
Theorem 3.3], we have the following statement about the spectrum of
$J_{\bm{A}}(\bm{\theta}^{\ast})$:

######  Proposition III.3.

The spectrum of $J_{\bm{A}}(\bm{\theta}^{\ast})$ can be defined in terms of:

  * •

The spectra of the Laplacian matrices of each layer, $\bm{L}_{\bm{A}_{i}}$;

  * •

The Jacobian of the reduced system,
$J_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast})$.

Namely,

$\operatorname{{\rm
Spec}}(\bm{J}_{\bm{A}}(\bm{\theta}^{\ast}))=\bigcup_{i=1}^{d}\left\\{(-\operatorname{{\rm
Spec}}(\bm{L}_{\\\
bm{A}_{i}})\setminus\\{0\\})-\lambda_{i}\right\\}\cup\operatorname{{\rm
Spec}}(\bm{J}_{\bar{\bm{A}}}(\bar{\bm{\theta}}^{\ast}))\,.$ (III.12)

In Sec. II, we showed how we can use the reduced system to broadcast solutions
to the multilayer Kuramoto network. We can now study the stability of these
solutions.

######  Proposition III.4.

The fixed point $\bar{\bm{\theta}}^{\ast}$ for the reduced system is linearly
stable if and only if the corresponding broadcasted fixed point for the
multilayer system is linearly stable.

###### Proof.

Assume that $\bar{\bm{\theta}}^{\ast}$ is linearly stable, then
$J_{\bar{A}}(\bar{\bm{\theta}}^{\ast})$ is symmetric negative-semidefinite. By
Sylvester’s criterion [42, Theorem 3.3.12], we know that $\lambda_{i}\geq 0.$
Furthermore, by Proposition III.1, we know that all eigenvalues of
$\bm{L}_{\bm{A}_{i}}$, except 0, are positive. These two facts show that all
elements in

$\bigcup_{i=1}^{d}\left\\{(-\operatorname{{\rm
Spec}}(\bm{L}_{\bm{A}_{i}})\setminus\\{0\\})-\lambda_{i}\right\\},$

are negative. By Proposition III.3 and the assumption that
$\bar{\bm{\theta}}^{\ast}$ is linearly stable, we conclude that all
eigenvalues of $\bm{J}_{\bar{A}}(\bm{\theta}i^{\ast})$, other than $0$, must
be negative. Hence $\bm{\theta}^{\ast}$ is linearly stable.

The other direction is straightforward. Indeed, by Proposition III.3, the
spectrum of $\bm{J}_{\bar{A}}(\bar{\bm{\theta}}^{\ast})$ is a subset of the
spectrum of $\bm{J}_{\bm{A}}(\bm{\theta}^{\ast})$. Therefore, if
$\bm{\theta}^{\ast}$ is linearly stable, then $\bar{\bm{\theta}}^{\ast}$ is
linearly stable as well. ∎

##  IV Numerical simulations and examples

In this section, we show examples and numerical simulations of Kuramoto
networks on multilayer systems and their respective reduced version that we
introduce in this paper.

Figure 1 shows a graphic representation of this approach. We study Kuramoto
oscillators on multilayer networks, where each layer is composed of $N$ nodes
and has an internal connection structure and coupling (intra-coupling). These
layers ($M$ in total) are then coupled with an external connection structure
(inter-coupling), thus leading to a two-levels coupling scheme (left). In this
paper, we introduced a reduced “inter-layers” representation that allows us to
investigate on the dynamics of the multilayer system in a simplified way. In
this simpler representation, each layer is given by a node, which is connected
to other nodes following the inter-coupling scheme (right). As shown in Sec.
II, the reduced system and the multilayer network have equivalent dynamics and
we can broadcast solutions from the former to the latter.

Figure 1: Graphic representation of the system. The multilayer network is
given by $M$ layers with $N$ nodes each one. Here, we use a two-levels system,
where each layer has internal connections in addition to the inter-layer
connections. Our approach allows us to study the dynamics of these system
using a reduced $M\times M$ system. This figure depicts an example with $M=3$
layers.

###  IV.1 Phase synchronization

We first consider a mutilayer network with $M=3$, where each layer is composed
of $N=100$ nodes. The internal connection structure of each layer is given by
ring with periodic boundary conditions, described by a circulant matrix where
each node has $10$ connections to each direction. The adjacency matrix for the
internal coupling is given by $A_{ij,\,\text{intra}}=1$ if nodes $i$ and $j$
are connected, and by $A_{ij,\,\text{intra}}=0$ otherwise. The inter-layers
coupling scheme is given by a complete graph, where all nodes in one layer are
connected to all nodes in another layer with coupling strength given by
$\epsilon$.

The first example is depicted in Fig. 2. Here, the initial state is given by
randomly selected phases for each layer, which are represented in the unitary
circle in color-code (Fig. 2a). All oscillators in each layer have the same
phase, therefore each layer is phase synchronized, but the entire, multilayer
system is not – in the beginning of the simulation. As time evolve, the
multilayer system transitions to phase synchrony due to the inter-layer
coupling, which leads all oscillators in all layers to depict the same phase
in the final state (Fig. 2b).

Figure 2: Dynamics on a $3$ layers system and its reduced system with random
initial conditions. In the beginning of the simulation, each layer is phase
synchronized but in different phases (a), therefore the multilayer system is
nonsynchronized. As time evolves, due to the inter-layer coupling, the
internal dynamics changes and the multilayer system transition to phase
synchronization (b). Here, the dynamics of the multilayer system and the
reduced one are the same, which is expressed by the Kuramoto order parameter
$R(t)$ (c). The trajectories for both the reduced system and the multilayer
one are shown in (d), where we can observe that (i) these system have an
equivalent dynamics; and (ii) the network transitions to phase synchrony as
time evolves.

The change on the dynamics over time is represented by the Kuramoto order
parameter (Fig. 2c), which is evaluated through Eq. (II.11) for the multilayer
network and through Eq. (II.12) for the reduced system. Here, one can observe
that both systems have an equivalent dynamics: at $t=0$, $R$ has a low value
due to the random phases on the initial state; due to the coupling, $R(t)$
increases as time evolves and, at $t\approx 2$, the order parameter reaches
the unity once both systems reach phase synchronization.

In order to show details on the dynamical equivalence between the reduced
system and the multilayer network, we plot the trajectories given by the
phases of the oscillators in Fig. 2d. Here, the phases in the reduced systems
are represented by the solid-shaded lines, and the phases of the oscillators
in the multilayer network are given by the dashed lines. This result
emphasizes that: (i) these system have an equivalent dynamics, and also (ii)
the network transitions to phase synchrony as time evolves, since the
oscillators converge to the same phase.

###  IV.2 Unstable twisted states

We use the same system to analyze a different kind of dynamical behavior:
phase-locking states, where a constant phase difference is observed across
oscillators [31]. This kind of state is also known as “twisted states”, and,
for a network with $M$ units, the $p^{\mathrm{th}}$ twisted state is given by:

$\bm{\theta}^{(p)}=\left(0,\frac{-2\pi p}{M},\cdots,\frac{-2\pi
p(M-1)}{M}\right).$ (IV.1)

We use it equation to obtain the twisted states for the reduced system and
then use the broadcasting approach to extend this to the multilayer system.

We consider here the same network: a $3$ layered system, where each layer is
composed of $N=100$ nodes, and the internal connection structure of each layer
is given by ring with periodic boundary conditions. The reduced system is then
given by a $3$ nodes networks, which are coupled in a all-to-all scheme. A
possible phase-locking state for this network, following Eq. (IV.1) is
represented in Fig. 3a, where constant difference of $\nicefrac{{2\pi}}{{3}}$
is observed. As stated before, all oscillators within a layer have the same
phase, so each layer is phase synchronized, but the entire system, in this
particular case, has $R=0$, given the phase-locking solution.

Figure 3: Dynamics on a $3$ layers systems and its reduced system on a phase-
locking state. We consider a twisted state for the reduced system following
Eq. (IV.1) with $M=3$ and $p=1$. This leads to a initial state where a
constant phase difference is observed across nodes (a). However, this twisted
state is not stable, so as time evolves, a transition to phase synchrony is
observed, and all oscillators have the same phase in the final state (b). The
time to this transition occurs, however, depends on the perturbation applied
to the system. We considered different perturbations amplitudes (color-code)
and analyze the Kuramoto order parameter as a function of time for the reduced
system (solid-shaded lines) and for the multilayer network (dashed lines) (c).
The bigger the perturbation the faster the transition and our approach shows a
perfect match between the dynamics of the reduced system and the multilayer
network.

This state, however, is not stable, so small perturbation, e.g. due to the
coupling or perturbation on the initial state, can lead the system to
transition to phase synchronization (Fig. 3b). Our approach is able to capture
the details on this transition. To show this point, we perform a detailed
simulation protocol, where different perturbations are applied to the initial
state. This perturbation is given by a uniform, random distribution of phases,
multiplied by an amplitude factor. Mathematically, the perturbation can be
described by: $\mathcal{A}\,\mathcal{U}(-\pi,\pi)$, where $\mathcal{A}$ is the
perturbation amplitude. Figure 3c depicts the Kuramoto order parameter as a
function of time when different amplitude are considered, which are
represented in color-code. Here, the order parameter for the reduced system is
shown in the solid-shaded lines, and the order parameter for the multilayer
systems is shown in the dashed lines. We observe that the smaller the
perturbation, the longer it takes to the phase-locking state to transition to
phase synchrony. We emphasize that our reduced approach has an equivalent
dynamics to the multilayered system, so we observe a perfect match.

###  IV.3 Higher number of layers

The approach introduced here through the reduced system and the broadcasting
process is general to study multilayer system with an arbitrary number of
layers and internal connection scheme. To show this point, we then analyze a
multilayer system composed of $M=50$ layers, where each layer is given by a
random network with $N=100$ nodes. Here, all layers are connected, such that
the reduced system is given by a complete graph with $M=50$ nodes.

At first, we consider the case of random initial conditions for this system.
Figure 4a show the trajectories for the reduced (solid-shaded lines) and the
multilayer systems (dashed lines) based on the phases of the oscillators
($\bar{\theta}$ and $\theta$). Due to the random initial conditions (see the
inner panel for initial state), the systems are desynchronized in the
beginning of the simulation. However, as time evolves, phase synchronization
is reached and all oscillators have the same phase (see inner panel for final
state). This result also shows the equivalence in the dynamical behavior of
both systems.

Figure 4: Dynamics on a $50$ layers network and its reduced system. We first
consider the case where the inter-layer coupling scheme is given by a complete
graph (all-to-all). So, random initial conditions evolve to phase synchrony.
Here, solid-shaded lines represent the reduced system and dashed lines
represent the multilayer network. The trajectories of the systems emphasize
this point (a) and also show the equivalent dynamics between the systems. When
we consider higher coupling strengths, the systems transition to phase
synchrony quicker (b). We then consider the case where the inter-layer
connection scheme is given by a first-neighborhood architecture, where each
layer is connected to two layers (one in each direction) with periodic
boundary conditions. In this case, phase-locking solutions are stable. We then
show that the systems assume the form of a twisted state (c); and that this is
robust to different perturbations (d).

Moreover, if we consider different coupling strengths between the layers
(inter coupling), the transition to phase synchrony occurs at different times,
as described in Proposition (II.1). To show that our approach is able to
capture these details, we consider the same random initial conditions and we
change the inter coupling strength. Figure 4b shows the order parameter for
the reduced system (solid-shaded lines) and for the equivalent multilayer
network (dashed lines) as a function of time for different values of
$\epsilon_{\mathrm{inter}}$. We observe that the higher the coupling the
faster the transition to phase synchronization.

We also use the approach introduced in this paper to analyze the broadcasting
a stable phase-locking states. To do so, we consider the multilayer system
with $M=50$ layers, each one being described by a random network with $N=100$
nodes. Each layer, however, is connected to two layers only, in a first-
neighborhood fashion. This leads to a reduced system being described as a ring
graph with $N=50$ nodes, where each node has degree $2$ with periodic boundary
conditions. This kind of network is known to support stable twisted states
[31]. Particularly, we use Eq. (IV.1) with $M=50$ and $p=1$ to generate the
stable phase-locking solution [31] for the reduced system, which is broadcast
to the multilayer network following the approach described in this paper.
Moreover, to show the stability of this solution and that our approach is able
to capture these details, we apply perturbation to the this phase-locking
state using the same approach as before introduced:
$\mathcal{A}\,\mathcal{U}(-\pi,\pi)$.

The trajectories for both, the reduced and the multilayer systems are shown in
Fig. 4c. In this case, the initial state is given by a twisted state as shown
in the inner panel. We then apply a perturbation to this state, which leads to
a change in the phase configuration. However, given the stability of the
$1^{\mathrm{st}}$ ($p=1$) twisted state for this system, despite the
perturbation, the phase-locking is recovered and we can observed the constant
phase difference across the oscillators, which leads to the final state being
exact the same as the initial one without perturbation (Fig. 4c inner panels).
Again, the phases for the reduced system are represented by the solid-shaded
lines and the phases for the multilayer networks are represented by the dashed
lines.

We also consider different perturbation amplitudes applied on the initial
state. Figure 4d shows the Kuramoto order parameter for both systems (solid-
shaded lines represent the reduced system and dashed lines the multilayer
network) as a function of time. The perturbation amplitude is shown in color-
code. We observe that the stronger the amplitude the bigger the variation of
the order parameter, but given the stability of this state, in the end $R=0$
for all cases, which characterizes the solution for these systems.

##  V Discussions and conclusions

In this paper, we have introduced an alternative approach to study and predict
the dynamical behavior of Kuramoto oscillators on a wide class of multilayer
networks. Based on the results presented in [18, 19], we can represent a
multilayer system as a join of matrices, each one representing a different
connection between nodes either in the same layer or in different layers. From
this representation, we have developed a “reduced system” that holds the
important information regarding the original, multilayer network. While the
multilayer network is composed of $M$ layers with $N$ nodes, the reduced
systems is given by $M$ elements. We have shown that, when the initial state
of the multilayer system has a relation with the initial state of the reduced
one, both systems have an equivalent dynamics. Thus, it allows us to
investigate into the dynamical behavior of the multilayer network in a simpler
way.

Based on this result, we now can find solutions for the reduced system and
broadcast them to the multilayer network. This method offers an alternative
and simple way to find equilibrium point for Kuramoto oscillators on
multilayer connection scheme. This approach is general to arbitrary topologies
for intra-layer connections (within each layer). Furthermore, we can use a
similar approach to obtain information on the linear stability of these
equilibrium points. We can write the Jacobian matrix for both systems and, by
using the results from [18], we can obtain the spectrum of the Jacobian of the
multilayer network based on the spectrum of the Jacobian of the reduced
system. Therefore, we are now able to investigate into the linear stability of
equilibrium points of multilayer networks in a simple way.

Multilayer networks have been used to model a variety of systems. Extensive
numerical studies have shown a rich repertory of dynamical behavior [43, 44,
45]. Moreover, experimental analyses confirm this feature [46, 47].
Furthermore, the use of multilayer networks have been helpful in the
understanding of complex systems, e.g. neural systems [4, 48]. However,
analytical treatment for this kind of system is still an open problem in these
diverse fields. Our approach thus offers a novel path in the study of
multilayer network, opening the possibility of analytical and mechanistic
insights on the dynamics of these important systems.

###### Acknowledgements.

This work was supported by BrainsCAN at Western University through the Canada
First Research Excellence Fund (CFREF), the NSF through a NeuroNex award
(#2015276), the Natural Sciences and Engineering Research Council of Canada
(NSERC) grant R0370A01, ONR N00014-16-1-2829, SPIRITS 2020 of Kyoto
University, Compute Ontario (computeontario.ca), Compute Canada
(computecanada.ca), and by the Western Academy for Advanced Research. J.M.
gratefully acknowledges the Western University Faculty of Science
Distinguished Professorship in 2020-2021. R.C.B gratefully acknowledges the
Western Institute for Neuroscience Clinical Research Postdoctoral Fellowship.
R.D. was supported by the Swiss National Science Foundation, under grant
number P400P2_194359.

## References

  * Watts and Strogatz [1998] D. J. Watts and S. H. Strogatz, Collective dynamics of ‘small-world’networks, Nature 393, 440 (1998).
  * Strogatz [2001] S. H. Strogatz, Exploring complex networks, Nature 410, 268 (2001).
  * Boccara [2010] N. Boccara,  _Modeling complex systems_ (Springer Science & Business Media, 2010).
  * Bassett  _et al._ [2018] D. S. Bassett, P. Zurn, and J. I. Gold, On the nature and use of models in network neuroscience, Nature Reviews Neuroscience 19, 566 (2018).
  * Banerjee  _et al._ [2016] T. Banerjee, P. S. Dutta, A. Zakharova, and E. Schöll, Chimera patterns induced by distance-dependent power-law coupling in ecological networks, Physical Review E 94, 032206 (2016).
  * Kinney  _et al._ [2005] R. Kinney, P. Crucitti, R. Albert, and V. Latora, Modeling cascading failures in the north american power grid, The European Physical Journal B-Condensed Matter and Complex Systems 46, 101 (2005).
  * Motter  _et al._ [2013] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Spontaneous synchrony in power-grid networks, Nature Physics 9, 191 (2013).
  * Dörfler  _et al._ [2013] F. Dörfler, M. Chertkov, and F. Bullo, Synchronization in complex oscillator networks and smart grids, Proc. Natl. Acad. Sci. USA 110, 2005 (2013).
  * Battiston  _et al._ [2020] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Networks beyond pairwise interactions: structure and dynamics, Physics Reports 874, 1 (2020).
  * Battiston  _et al._ [2021] F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. Ferraz de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, _et al._ , The physics of higher-order interactions in complex systems, Nature Physics 17, 1093 (2021).
  * Boccaletti  _et al._ [2014] 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 (2014).
  * Kivelä  _et al._ [2014] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, Multilayer networks, Journal of Complex Networks 2, 203 (2014).
  * De Domenico  _et al._ [2016] M. De Domenico, C. Granell, M. A. Porter, and A. Arenas, The physics of spreading processes in multilayer networks, Nature Physics 12, 901 (2016).
  * Silk  _et al._ [2018] M. J. Silk, K. R. Finn, M. A. Porter, and N. Pinter-Wollman, Can multilayer networks advance animal behavior research?, Trends in ecology & evolution 33, 376 (2018).
  * Braun  _et al._ [2015] U. Braun, A. Schäfer, H. Walter, S. Erk, N. Romanczuk-Seiferth, L. Haddad, J. I. Schweiger, O. Grimm, A. Heinz, H. Tost, _et al._ , Dynamic reconfiguration of frontal brain networks during executive cognition in humans, Proceedings of the National Academy of Sciences 112, 11678 (2015).
  * Zhang  _et al._ [2015] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Explosive synchronization in adaptive and multilayer networks, Physical review letters 114, 038701 (2015).
  * Della Rossa  _et al._ [2020] F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino, Symmetries and cluster synchronization in multilayer networks, Nature communications 11, 1 (2020).
  * Đoàn  _et al._ [2022] J. Đoàn, J. Mináč, L. Muller, T. T. Nguyen, and F. W. Pasini, Joins of circulant matrices, Linear Algebra and its Applications 650, 190 (2022).
  * Doan  _et al._ [2022] J. Doan, J. Mináč, L. Muller, T. T. Nguyen, and F. W. Pasini, On the spectrum of the joins of normal matrices and applications, arXiv e-prints , arXiv (2022).
  * Chebolu  _et al._ [2022] S. K. Chebolu, J. L. Merzel, J. Mináč, L. Muller, T. T. Nguyen, F. W. Pasini, and N. D. Tân, On the joins of group rings, arXiv preprint arXiv:2208.07413 (2022).
  * Kuramoto [1975] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in  _International symposium on mathematical problems in theoretical physics_ (Springer, 1975) pp. 420–422.
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, The kuramoto model: A simple paradigm for synchronization phenomena, Reviews of Modern Physics 77, 137 (2005).
  * Rodrigues  _et al._ [2016] F. A. Rodrigues, T. K. D. M. Peron, P. Ji, and J. Kurths, The kuramoto model in complex networks, Physics Reports 610, 1 (2016).
  * Pluchino  _et al._ [2006] A. Pluchino, V. Latora, and A. Rapisarda, Compromise and synchronization in opinion dynamics, The European Physical Journal B-Condensed Matter and Complex Systems 50, 169 (2006).
  * Breakspear  _et al._ [2010] M. Breakspear, S. Heitmann, and A. Daffertshofer, Generative models of cortical oscillations: neurobiological implications of the kuramoto model, Frontiers in Human Neuroscience 4, 190 (2010).
  * Bick  _et al._ [2020] 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 (2020).
  * Arenas  _et al._ [2008] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Synchronization in complex networks, Physics Reports 469, 93 (2008).
  * Boccaletti  _et al._ [2002] S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, and C. Zhou, The synchronization of chaotic systems, Physics reports 366, 1 (2002).
  * Strogatz [2000] S. H. Strogatz, From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Nicosia  _et al._ [2013] V. Nicosia, M. Valencia, M. Chavez, A. Diaz-Guilera, and V. Latora, Remote synchronization reveals network symmetries and functional modules, Physical Review Letters 110, 174102 (2013).
  * Townsend  _et al._ [2020] 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, 083142 (2020).
  * Delabays  _et al._ [2016] R. Delabays, T. Coletta, and P. Jacquod, Multistability of phase-locking and topological winding numbers in locally coupled kuramoto models on single-loop networks, Journal of Mathematical Physics 57, 032701 (2016).
  * Abrams and Strogatz [2004] D. M. Abrams and S. H. Strogatz, Chimera states for coupled oscillators, Physical review letters 93, 174102 (2004).
  * Parastesh  _et al._ [2021] F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti, and M. Perc, Chimeras, Physics Reports 898, 1 (2021).
  * Xu  _et al._ [2018] C. Xu, S. Boccaletti, S. Guan, and Z. Zheng, Origin of bellerophon states in globally coupled phase oscillators, Physical Review E 98, 050202 (2018).
  * Bovier [2006] A. Bovier,  _Statistical Mechanics of Disordered Systems: A Mathematical Perspective_, 1st ed. (Cambridge University Press, 2006).
  * Chayes [2009] L. Chayes, Mean Field Analysis of Low–Dimensional Systems, Commun. Math. Phys. 292, 303 (2009).
  * Kadanoff [2009] L. P. Kadanoff, More is the Same; Phase Transitions and Mean Field Theories, J. Stat. Phys. 137, 777 (2009).
  * Strogatz and Mirollo [1988] S. H. Strogatz and R. E. Mirollo, Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random intrinsic frequencies, Physica D 31, 143 (1988).
  * Dörfler and Bullo [2014] F. Dörfler and F. Bullo, Synchronization in complex networks of phase oscillators: A survey, Automatica 50, 1539 (2014).
  * Horn and Johnson [2012] R. A. Horn and C. R. Johnson,  _Matrix analysis_ (Cambridge University Press, 2012).
  * Bazaraa  _et al._ [2013] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty,  _Nonlinear programming: theory and algorithms_ (John Wiley & Sons, 2013).
  * Kasatkin  _et al._ [2017] D. Kasatkin, S. Yanchuk, E. Schöll, and V. Nekorkin, Self-organized emergence of multilayer structure and chimera states in dynamical networks with adaptive couplings, Physical Review E 96, 062211 (2017).
  * Kachhvah and Jalan [2021] A. D. Kachhvah and S. Jalan, Explosive synchronization and chimera in interpinned multilayer networks, Physical Review E 104, L042301 (2021).
  * Majhi  _et al._ [2017] S. Majhi, M. Perc, and D. Ghosh, Chimera states in a multilayer network of coupled and uncoupled neurons, Chaos: an interdisciplinary journal of nonlinear science 27, 073109 (2017).
  * Blaha  _et al._ [2019] K. A. Blaha, K. Huang, F. Della Rossa, L. Pecora, M. Hossein-Zadeh, and F. Sorrentino, Cluster synchronization in multilayer networks: A fully analog experiment with l c oscillators with physically dissimilar coupling, Physical review letters 122, 014101 (2019).
  * Leyva  _et al._ [2017] I. Leyva, R. Sevilla-Escoboza, I. Sendiña-Nadal, R. Gutiérrez, J. Buldú, and S. Boccaletti, Inter-layer synchronization in non-identical multi-layer networks, Scientific Reports 7, 1 (2017).
  * Palmigiano  _et al._ [2017] A. Palmigiano, T. Geisel, F. Wolf, and D. Battaglia, Flexible information routing by transient synchrony, Nature neuroscience 20, 1014 (2017).