# Dynamics of multilayer networks with amplification

Thierry Njougouo  Research Unit Condensed Matter, Electronics and Signal
Processing, University of Dschang, P.O. Box 67 Dschang, Cameroon.  Victor
Camargo  Center for Interdisciplinary Research on Complex Systems,University
of Sao Paulo, Av. Arlindo Bettio 1000, 03828-000 São Paulo, Brazil.
Department of Physics-FFCLRP, University of São Paulo, Ribeirao Preto-SP,
14040-901, Brasil.  Patrick Louodop  Research Unit Condensed Matter,
Electronics and Signal Processing, University of Dschang, P.O. Box 67 Dschang,
Cameroon.  Fernando Fagundes Ferreira  Center for Interdisciplinary Research
on Complex Systems,University of Sao Paulo, Av. Arlindo Bettio 1000, 03828-000
São Paulo, Brazil.  Department of Physics-FFCLRP, University of São Paulo,
Ribeirao Preto-SP, 14040-901, Brasil.  Pierre K. Talla  L2MSP, University of
Dschang, P.O. Box 67 Dschang, Cameroon.  Hilda A. Cerdeira  São Paulo State
University (UNESP), Instituto de Física Teórica, Rua Dr. Bento Teobaldo Ferraz
271, Bloco II, Barra Funda, 01140-070 São Paulo, Brazil.

###### Abstract

We study the dynamics of a multilayer network of chaotic oscillators subject
to an amplification. Previous studies have proven that multilayer networks
present phenomena such as synchronization, cluster and chimera states. Here we
consider a network with two layers of Rössler chaotic oscillators as well as
applications to multilayer networks of chaotic jerk and Liénard oscillators.
Intra-layer coupling is considered to be all to all in the case of Rössler
oscillators, a ring for jerk oscillators and global mean field coupling in the
case of Liénard, the inter-layer coupling is unidirectional in all these three
cases. The second layer has an amplification coefficient. An in-depth study on
the case of a network of Rössler oscillators using master stability function
and order parameter leads to several phenomena such as complete
synchronization, generalized, cluster and phase synchronization with
amplification. For the case of Rössler oscillators, we note that there are
also certain values of coupling parameters and amplification where the
synchronization doesn’t exist or the synchronization can exist but without
amplification. Using other systems with different topologies, we obtain some
interesting results such as chimera state with amplification, cluster state
with amplification and complete synchronization with amplification.

††preprint: AIP/123-QED

> Research on multilayer networks has attracted a lot of attention in recent
> years in many areas of physics, engeneering, social sciences etcLeyva  _et
> al._ (2018); De Domenico  _et al._ (2013); Zhang  _et al._ (2015); Jalan,
> Kachhvah, and Jeong (2020); Belykh, Carter, and Jeter (2019); Blaha  _et
> al._ (2019). Some emergent behaviors in such systems due to interaction
> among the dynamical units reveal a variety of interesting phenomena, such as
> synchronization Arenas  _et al._ (2008); Leyva  _et al._ (2018); del Genio
> _et al._ (2016), cluster formation Zhang  _et al._ (2001), explosive
> synchronization Ji  _et al._ (2014), explosive desynchronization Lauter,
> Mitra, and Marquardt (2017), chimera Kuramoto and Battogtokh (2002); Schmidt
> _et al._ (2014); Schmidt and Krischer (2015); Mishra  _et al._ (2015);
> Dudkowski, Maistrenko, and Kapitaniak (2014) etc. Among these,
> synchronization and chimeras are the most widely studied. The notion of
> amplification is very important in science and technology. This work
> presents an investigation of different phenomena such as complete
> synchronization, cluster formation, phase synchronization and chimera states
> in a network with amplification. For an extended study we present three
> cases with three different topologies.

##  I Introduction

The structure of many real world problems in nature, engineering, science and
technology is defined as a set of entities interacting with each other in
complicated patterns that can produce multiple types of relationships which
change in time and exhibits a plethora of emergent patterns or behaviors as
synchronization, chimeras, chaos, consensus, cooperation, only to mention a
few. One very important ingredient behind most of those phenomena is the way
in which the particles or agents forming the systems interact, i.e., the
topology of the underlying network.

Network theory is used as an important tool for the modelling of dynamical
processes in complex systems Arenas  _et al._ (2008); Leyva  _et al._ (2018);
De Domenico  _et al._ (2013). It plays also a major role in the investigation
of collective behavior. It finds many applications in epidemiology where it is
used to investigate epidemic spreading, in industry where it is used in
control of behavior of machines, in dynamics of populations with the control
of displacement of the individuals, the cars, the drones or the airplanes.
According to the these applications, we can mention that the main objective is
the controllability of the network to lead to a certain state (it can be
synchronization, cluster state, phase-flip, chimera state, etc.) del Genio
_et al._ (2016); Castellano, Fortunato, and Loreto (2009); Pastor-Satorras
_et al._ (2015); Liu, Slotine, and Barabási (2011); Schmidt  _et al._ (2014);
Schmidt and Krischer (2015). Thus, the investigation of the dynamics of the
networks need the expertise of some mathematical tools such as the Master
Stability Function (MSF) developed by Pecora and Carroll, the transversal
Lyapunov exponent, the correlation between the oscillators of the same or of a
different layer, etc. One of the best methods to study the stability of the
synchronization in the network is the Master Stability Function Pecora and
Carroll (1998). This method is used for coupled identical oscillators.

Many researchers are studying several phenomena which take place in multiplex
networks. Such interest is motivated in understanding how the complete or
partial synchronization occurs in this type of systems and also because the
topologies of multiplex networks appear in several natural and tecnological
systems. Multiplex network may be described as being a collection of two or
more coupled networks where a set of networks is connected by links where the
interactions are of different types Boccaletti  _et al._ (2014); Mondragon and
Bianconi (2018). These links characterize the connections existing between any
node or network of the multiplex network. Many recent works addressing
multilayer structures and systems were summarized in Boccaletti  _et al._
(2016).

The study of inter-layer synchronization in non-identical multilayer networks
was address in Leyva  _et al._ . The authors were able to show an analytical
treatment for a two-layer multiplex using the Master Stability Function
method. One interesting outcome was to predict the effect that missing links
in one of the layers has on the inter-layer synchronization. Later, in Ghosh,
Zakharova, and Jalan (2018) it was found that a sparse inhomogeneous second
layer can promote chimera states in a sparse homogeneous first layer. The
study of synchronization of non-identical multilayers is very recent, thus
many collective behavior properties and patterns may unravel.

Here, we consider a network with two layers where we choose an all to all
coupling in the layer for the case of Rössler oscillators, a ring with
bidirectional coupling in each layer for the case of jerk oscillators and
aglobal mean field coupling for the last case mentioned above. The connection
between the systems of both layers (interlayer coupling) is unidirectional.
The main goal of this work is to investigate in each case the dynamics of each
layer as well as the whole network with amplification in the second layer. We
found the key values of the parameters to control synchronization with and
without amplification.

The remainder of this work is organized as follows. In section II, we present
our multilayer network with the mathematical description of the model and the
systems. The dynamics of the main case is presented in section III emphasizing
on the intralayer and interlayer coupling in which numerical simulations are
done. An application to another two systems has been studied in section IV and
finally, we present the conclusions in section V.

##  II Multilayer network

The model consists in a multilayer network constituted of $N$ nodes connected
in each layer, which can be represented by a $2N\times 2N$ adjacency matrix
$A_{ij}$ where the elements of this matrix are respectively $1$ if the nodes
$i$ and $j$ are connected and $0$ if not. Based on Ref.Ghosh, Zakharova, and
Jalan (2018), the adjacency matrix of the whole network consisting two layers
can be expressed as follow:

$\displaystyle\textbf{A}=\left({\begin{array}[]{*{4}{c}}A^{1}&0\\\ I&A^{2}\\\
\end{array}}\right)$ (3)

where $A^{1}$ and $A^{2}$ are the $N\times N$ adjacency matrix modelling the
intralayer connectivity in the first and second layer respectively. $I$ is an
$N\times N$ identity matrix representing the unidirectionnal interactions
(Layer 1 -> Layer 2) between the oscillators with the same index in both
layers. The use of the null matrix is justified by the non-existence of a
connection from the slave layer to the master layer. In the following, we
consider a model of multilayer network constituted of $N$ nodes in each layer
connected using an all-to-all coupling scheme in each layer (see Fig.1). To
each node corresponds a nonlinear autonomous Rössler oscillator as described
inRössler (1979). Notice that it is this combined oscillator which defines our
network as a two-layer system made up of a driving system and a slave one.

Figure 1: Schematic representation of a network with two layers of
interaction.

The dynamics of the first layer also considered as the driver for the network
is described by Eq.4 where $\epsilon_{1}$ is the intra-layer coupling strength
for the first layer.

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{x}_{i}^{1}=-x_{i}^{2}-x_{i}^{3}+{\epsilon_{1}}\sum\limits_{k=1}^{N}A_{ik}^{1}(x_{k}^{1}-x_{i}^{1}),\\\
\dot{x}_{i}^{2}=x_{i}^{1}+ax_{i}^{2},\\\
\dot{x}_{i}^{3}=bx_{i}^{1}+x_{i}^{3}(x_{i}^{1}-c).\end{array}\right.\par\end{array}$
(4)

The second layer has exactly the same intra-layer topology with similar
systems in the nodes. The choice of systems in the second layer follows
Louodop et al. Louodop  _et al._ (2017) where the authors show in Appendix A
that, this form of coupling between elements of different layers produces
generalized synchronizationLouodop  _et al._ (2017). Notice that it is this
system defined in Louodop et al. Louodop  _et al._ (2017), what defines which
oscillators in Layer 2 interact with those in Layer one, while keeping the
same all to all intra-layer topology. Therefore, the dynamics of the second or
slave layer is given by Eq.5:

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}{{\dot{y}}_{j}^{1}}=-{y_{j}^{2}}-{y_{j}^{3}}+{\epsilon_{2}}\sum\limits_{k=1}^{N}A_{jk}^{2}({y_{k}^{1}}-{y_{j}^{1}})+C_{0}(x_{j}^{1}-C_{2}y_{j}^{1}),\\\
{{\dot{y}}_{j}^{2}}=\frac{{x_{j}^{1}}+a{x_{j}^{2}}}{C_{2}},\\\
{{\dot{y}}_{j}^{3}}=b{y_{j}^{1}}+\frac{{{x_{j}^{3}}{x_{j}^{1}}}}{{{C_{2}}}}-c{y_{j}^{3}}.\end{array}\right.\end{array}$
(5)

Here, $\epsilon_{2}$ is the intra-layer coupling strength and $C_{0}$ is the
interlayer coupling strength. It is important to mention that, this inter-
layer coupling exists only between the oscillators with the same index (i.e.
for $i=j$ where the $j$ index runs from $1$ to $N$). $C_{2}$ is the parameter
of proportionality named amplification coefficient. This term or interaction
via a $conjugate$ variable has been used in the literature to model $revival$
Zou  _et al._ (2015); Ghosh, Banerjee, and Kurths (2015) as well as
$amplitude$ $death$ Karnatak, Ramaswamy, and Prasad (2007). For all these
layers we consider $a=0.36$, $b=0.4$ and $c=4.5$ and we note that at these
values of the parameters the systems operate in the chaotic regimeRössler
(1979). This topology of connectivity between the master and the slave layer
imposes generalized synchronization between both layers in the absence of
intralayer coupling because as it is conceived (see Appendix A), the slave
layer is supposed to function as an observer of the master layer with certain
conditions $C_{2}\neq 0$. In this work, we were interested in the notion of
synchronization with amplification (or reduction) depending on the value of
$C_{2}$. For $C_{2}>1$ we have an amplification of the systems of the master
layer (or a reduction of the systems of the slave layer) and conversely for
$C_{2}<1$. It should be noted that when $C_{2}$ takes negative values there is
an anti-synchronization between the systems of the master layer and those of
the slave layer having the same index. The amplification coefficient must be
different from zero ($C_{2}\neq 0$) and bounded, because if $C_{2}=0$ the
systems of the slave layer will diverge and if $C_{2}$ is too large the
slaving factor will tend to zero. Therefore we keep $C_{2}$ in the interval
between $0.005$ and $2$. This type of topology finds applications in many
domains such as aircraft control where recently they have proven that an
optimal control permitted to regulate the air traffic in the sky Cafieri  _et
al._ (2018). Considering the domains of application of this network topology,
it seems rather important to investigate the dynamics of this network with
different parameters, which is presented in the next section.

##  III Dynamics of networks

###  III.1 Dynamics of Different Layers

Considering the topology given in Fig.1 and the mathematical equations for
both layers Eq.4 and Eq.5, respectively, we are going to investigate
numerically the dynamics in the different layers using the MSF described in
Appendix B. A numerical calculation is done using Runge-Kutta fourth order for
a long time simulation and the permanent solutions are considered at
$t_{min}=0.6$ $t_{max}$.

Figure 2: Dynamics of the two layers of the network: (a) Master stability
function of the first layer as a function of the intralayer coupling
considering $\epsilon_{2}=10\epsilon_{1},C_{0}=1$ and $C_{2}=2$. (b) Order
parameter showing the dynamics of the slave layer for different values of the
amplification coefficient for $\epsilon_{2}=10\epsilon_{1}$ and $C_{0}=1$.
(c,d,e and f) Mean phase of the first and second layer respectively for
$\epsilon_{1}=0.005,0.007,0.009,0.04,C_{0}=1$ and $C_{2}=2$. (g) Attractors of
the oscillators $1$, $6$, $11$, $16$, $21$ and $26$ in the master and slave
layer for $\epsilon_{1}=0.007$. (h) Attractors of the oscillators $1$, $6$,
$11$, $16$, $21$ and $26$ in the master and slave layer $\epsilon_{1}=0.009$
respectively. $S_{i}L_{l}$ means system $i$ of the layer $l$ ($i=1,2,...,N$
and $l=1,2$).

To illustrate the behaviour of the oscillators of the first layer we present
in Fig.2(a) the Master Stability Function or Largest Lyapunov Exponent ($LLE$)
of the variational equation Eq. 29, used to characterize the stability of the
synchronization in the first layer. In this figure there are two important
regions in terms of characterization of the stability of the synchronization:
if the $LLE\leq 0$ the synchronization is stable and the $LLE>0$ the
synchronization is unstable. By varying smoothly the intralayer coupling
$\epsilon_{1}$, this figure shows that, when $\epsilon_{1}$ increases the
systems evolve to the synchronous state at $\epsilon_{1}=0.009$. This
synchronization is obtained at a minimal value of the coupling strength in the
first layer and we note that the synchronization in the slave layer is highly
influenced not only by the coupling but also by the amplification, this can be
seen on the behavior of the order parameter Chen  _et al._ (2017); Hong and
Choi (2000) which is plotted in Fig.2(b) as a function of $\epsilon_{1}$. The
procedure to compute the order parameter is explained in Appendix C. It should
be noted that the notion of amplification here is related to the amplitude of
the oscillations of the state variables of the systems. We are talking about
synchronization with amplification if and only if the ratio $X/Y=C_{2}$ is
respected, as shown in Appendix A. So, we investigate here the impact of the
amplification parameter on the dynamics of the network, where the order
parameter of the slave layer is plotted for different values of $C_{2}$.
According to this figure although interaction between the layers does not
impede synchronization in the slave layer, it becomes more effective for small
values of $C_{2}$. Based on the demonstration given in Appendix A
synchronization occurs at $Y=\frac{X}{C_{2}}$. Therefore, when $C_{2}<1$ we
obtain an amplification in the slave layer and respectively a reduction in
master layer and vice-versa when $C_{2}>1$. Thus, if we need to achieve
synchronization in both layers at the same value of the interlayer coupling,
$C_{2}$ must be very small, leading to a significant amplification at the
second layer. Also from Fig.2(b) we see that for the minimum value of $C_{2}$
= $0.005$, the synchronization in both layers happens at the corresponding
values of $\epsilon_{1}$ and $\epsilon_{2}$ ($\epsilon_{2}=10\epsilon_{1}$),
this is represented in Fig.2(b), where the order parameter $r_{1}$ and $r_{2}$
for both layers is seen to reach the value for complete synchronization at the
same point. This value of $C_{2}$ produces an amplification of around $200$ of
all the variables of the master layer $(Y\approx 200X)$ in the slave layer. In
Fig.2(c) we plot the mean phase of the driver and the slave layer (the first
30 systems are for the driver layer and the rest for the slave layer) for the
value $\epsilon_{1}=0.005$ of the intralayer coupling. The situation shown
here is confirmed by the Fig.2(a) ($LLE>0$). If we consider
$\epsilon_{1}=0.007$, Fig.2(a) shows that the Largest Lyapunov Exponent is non
negative but very close to zero, then in the Fig.2(d) we show the mean phase
where the master layer has a two cluster synchronization with equal phases
Zhang  _et al._ (2001); Hu  _et al._ (2000); Ji  _et al._ (2014); Maistrenko
_et al._ (2004); Allefeld, Müller, and Kurths (2007). In the slave layer,
while the clusters follow the systems with the same index as that of the
master layer, we see an oblique sliding of the systems reminding of a splay
state. To better appreciate the dynamics of the oscillators in this behaviour,
we show in Fig.2(g) the attractors of the oscillators labeled
$1,6,11,16,21,26$ in the master and slave layers for $C_{2}=2$ and
$\epsilon_{1}=0.007$ ( $LLE\geq 0$). For $\epsilon_{1}=0.009$ (with $LLE<0$),
we obtain Fig.2(e) which represents the synchronization in the first layer and
a coherent oblique sliding in the second layer. By computing the phase
difference between consecutive oscillators (here consecutive refers to the
indices of the oscillators) we verified that the phase distance between
oscillators of consecutive index in the slave layer is constant, therefore the
second layer presents indeed a phenomenon of splay Singha  _et al._ (2017);
Nichols and Wiesenfeld (1992). In Fig.2(f) we show for $\epsilon_{1}=0.04$ the
phase synchronization in both layers but not at the same value of the mean
phase. So, for the multilayer network, the dynamics is equivalent to that of
two clusters. To illustrate the dynamics of the system we show in Fig.2(h) the
attractors of some oscillators (labeled $1,6,11,16,21,26$) for
$\epsilon_{1}=0.04$ to appreciate the behaviour of the entire network.

Therefore, we can conclude that, according to the different values of the
intralayer coupling the network leads to different phenomena such as, cluster
synchronization, splay, synchronization and the stability of this
synchronization is confirmed by the MSF.

###  III.2 Impact of the interlayer coupling and amplification on the dynamics
of the network

In the previous Sections we have shown the influence of the intralayer
coupling on the dynamics of the network. We see that the network presents many
phenomena depending on the amplification coefficient $C_{2}$ and the
interlayer coupling $C_{0}$. In this section, our goal is twofold: first we
investigate the behaviour of the network under the impact of these two
parameters and then we show the effect of the amplification as well as its
effectiveness in the network. We keep the intra-coupling constants
($\epsilon_{1},\epsilon_{2}$) fixed, varying smoothly the amplification
coefficient ($C_{2}$) from 0.005 to 2 and the inter-coupling ($C_{0}$) from 0
to 20.

Figure 3: Network synchronization regions: (a) Order parameter of the
multilayer network for $\epsilon_{1}=0.03$ and $\epsilon_{2}=0.25$ as a
function $C_{0}$ and $C_{2}$. (b) Number of states in the slave layer for
$\epsilon_{1}=0.03$ and $\epsilon_{2}=0.25$ as function $C_{0}$ and $C_{2}$.
Two parameter phase diagram by simultaneously varying the intra-layer coupling
$\epsilon_{1}$ and inter-layer coupling $C_{0}$ with
$\epsilon_{2}=10\epsilon_{1}$ and for different values of the amplification
coefficient $C_{2}$: (c) $C_{2}=0.5$, (d) $C_{2}=1$ and (e) $C_{2}=2$.

As mentioned in Fig.2(b), the synchronization in the slave layer is imposed by
the synchronization in the master layer. According to the literature Fujisaka
and Yamada (1983); Pecora and Carroll (1990), in order to bring our multilayer
network towards a desired behaviour such as synchronization, cluster
formation, splay and so on, it suffices to adjust the coupling. Although this
is usually the case, in our system we have two important parameters acting as
interlayer couplings ($C_{0}$ and $C_{2}$) with the difference that one of
them ($C_{2})$ allows us to increase or decrease the amplitude of the signal
in one of these layers. To illustrate the evolution towards the
synchronization as a function of the amplification parameter $C_{2}$ and the
interlayer coupling $C_{0}$ we show in Fig.3(a) the order parameter of the
slave layer, since we are interested in synchronization throughout the
network. By considering the intralayer coupling $\epsilon_{1}=0.03$ and
$\epsilon_{2}=0.25$ all the oscillators of the first layer are synchronized
but due to the effects of the amplification parameter $C_{2}$ and the
interlayer coupling, the dynamics of the slave layer can not be the same.
Therefore, the dynamics of the multilayer network and particularly the slave
layer is greatly affected by $C_{0}$ and $C_{2}$. This can be observed in
Fig.3(a) where we plot the order parameter of the second layer. In red we
represent the region where synchronization and amplification is obtained. We
note that in this red zone the relation $X=C_{2}Y$ is verified with a
precision of $10^{-4}$, while in the blue domain it is possible to obtain
synchronization in the first layer only but not in both layers of the network.
This information is corroborated by the number of states of the second layer
shown in Fig.3(b). This Fig.3(b) confirms the order parameter by presenting a
number of states equal to one in the case of synchronisation (dark blue). When
the order parameter is different from one, the synchronisation of all the
systems in the network is not achieved. At the moment, we can find partial
synchronization called cluster which is characterized by a number of
independent states of the network much less than the total number of elements.
Based on Fig.3(b), we see a thin region with cluster formation in the area of
transition from synchronization to desynchronization and vice versa. To better
highlight the impact of the amplification ($C_{2}$) on the transition to
synchronization, we have presented the order parameter of the multilayer
network showing the transition to the synchronous state for three values of
the amplification $C_{2}=0.5$, $C_{2}=1$ and $C_{2}=2$, showed in Fig.3 (c, d
and e) respectively. Varying smoothly the intra and interlayer coupling for
these three fixed values of $C_{2}$, we have obtained different areas of
synchronization of the systems of the network. For these three figures the
first domain (D1) represents the region where there is no synchronization in
any layer and the second domain (D2) where the synchronization exists only in
the first layer. The third domain (D3) represents the zone of synchronization
of both layers and the last domain (D4) is showing the area where there is
divergence between the states of the oscillators of the slave layer. Here,
divergence means an infinite amplification of the amplitude of the
oscillations of the systems of the slave layer leading to an explosion
($y^{1},y^{2},y^{3}$ tend towards the inifinite). In summary, we notice that,
the area of synchronization of both layers increases when the amplification
coefficient decreases. So in Fig.3(c) this zone of synchronization of both
layers is the largest domain (D3) and the zone (D2) where only the master
layer synchronizes is almost nonexistent. Therefore when $C_{2}$ increases, it
takes stronger values of the coupling in the master layer for the slave layer
to become synchronized as can be seen in Fig.2(b). Fig.4 is used to confirm
and to present the transition to the synchronization using the Pearson
correlation Pearson (1920) (defined by Eq.6) between the first variables of
each oscillator of both layers, the time series of the first variables of the
oscillators label 1 and 15 of both layers and the mean phase of the
oscillators of the network.

$\rho(x^{1},y^{1})=\frac{\sum\limits_{i=1}^{N}(x_{i}^{1}-\bar{x}^{1})(y_{i}^{1}-\bar{y}^{1})}{\sqrt{\sum\limits_{i=1}^{N}(x_{i}^{1}-\bar{x}^{1})^{2}}\sqrt{\sum\limits_{i=1}^{N}(y_{i}^{1}-\bar{y}^{1})^{2}}}\\\
$ (6)

where $\bar{x}^{1}=\sum\limits_{i=1}^{N}x_{i}^{1}$ and
$\bar{y}^{1}=\sum\limits_{i=1}^{N}y_{i}^{1}$ are the mean of the states
variables $x_{i}^{1}$ and $y_{i}^{1}$. In this manuscript, the colour yellow
used in the correlation refers to the oscillators that synchronize and the
colour blue to those that do not synchronize.

Figure 4: Dynamics of the multi-layer network for $C_{0}$ and $C_{2}$: (a)
Correlation between $x_{i}^{1}$ and $y_{j}^{1}$ (Yellow color indicates the
oscillators in the synchronized state and blue the oscillators in an
unsynchronized state.) (b) Mean phase and (c) time series of some oscillators
of the network for $C_{2}=0.2$ and $C_{0}=10$. $S_{i}L_{l}$ means system $i$
of the layer $l$ (i = 1, 2, …, N and l = 1, 2). (d,e,f) Synchronization
between some oscillators in the first and second layer of the network for
$C_{2}=0.2$ and $C_{0}=10$. (g) Correlation between $x_{i}^{1}$ and
$y_{j}^{1}$. (h) Mean phase and (i) time series of some oscillators of the
network for $C_{2}=1.55$ and $C_{0}=10$. (j,k,l) Synchronization between some
oscillators in the first and second layer of the network for $C_{2}=1.55$ and
$C_{0}=10$. Figures (d,e,f) and (j,k,l) are special cases of the correlations
presented in Figures 4(b and g). However, these figures even if they may
confuse for being similar show that amplification clearly depends on the value
of $C_{2}$, and also show that in some cases the synchronisation of
oscillators of different index is possible (see Figs. 4d and 4k).

Fig.4(a) shows the correlation between the first variables of the oscillators
of the master and slave layer for $C_{0}=10$ and $C_{2}=0.2$ (these values are
taken almost at the border of the separation of the red and blue zone in
Fig.3(a)). This figure presents yellow and blue colors to identify those
oscillators of the multilayer network that synchronize and those who do not
synchronize, respectively. In this figure, the x-axis corresponds to the index
used to identify the oscillators of the master layer ($i$) and the y-axis
corresponds to the index used to identify the different oscillators of the
slave layer ($j$). According to this figure, the first 14 oscillators of both
layers synchronize. To present in detail the situation we plot in Fig.4(b) for
the same fixed values of the parameters the mean phase for oscillators in both
layers, where the labels 1 to 30 belong to the master layer and 30 to 60 to
the slave layer. Here we notice two different groups in each layer: the
synchronous group (first 14 oscillators of each layer) and the asynchronous
group, so the transition to the synchronization is done by cluster formation.
To show the effect of the amplification for these fixed values of parameters,
we have in Fig.4(c) the time series of the variables for oscillators labels
$1$,and $15$. This figure presents a chaotic evolution of these variables as a
function of time as well as the amplification, these effects of
synchronization with amplification are well appreciated in Fig.4(d,f) for the
oscillators indicated in the axes. According to Fig.4(a) the $15^{th}$
oscillator of the master layer and the $14^{th}$ oscillator of the slave layer
can not synchronize (see Fig.4(e)). At the end, to confirm the synchronization
with amplification in the multilayer network, we present in Fig.4(g) the
correlation between the elements of both layers for $C_{0}=10$ and
$C_{2}=1.55$. This correlation shows the synchronization of the multilayer
network which can be appreciated in Fig.4(h) where we plot the mean phase for
oscillators in both layers. The time series of the oscillators $1$ and $15$
(Fig.4(i)) show the amplification described in Eq.5. We can appreciate the
synchronization with amplification of the multilayer network in the
Fig.4(j,k,l) where complete synchronization of some chosen oscillators is
shown.

Figure 5: Road to divergence: (a) Variation of the phases of the slave layer
for $C_{2}=2$, $\epsilon_{1}=0.065721$ and $\epsilon_{2}=10\epsilon_{1}$. (b)
Zoom of the variation of the phases of the slave layer for $C_{2}=2$,
$\epsilon_{1}=0.065721$ and $\epsilon_{2}=10\epsilon_{1}$. (c)) Order
parameter of the slave layer for $C_{2}=2$, $\epsilon_{1}=0.065721$ and
$\epsilon_{2}=10\epsilon_{1}$. (d) Mean phase of the oscillators of the
network and (e) correlation between the oscillators of the slave layer for
$\epsilon_{1}=0.065721$ and $\epsilon_{2}=10\epsilon_{1}$, $C_{2}=2$ and
$C_{0}=18$. (The yellow colour indicates the oscillators that synchronize and
the blue colour those that do not synchronize).

To appreciate the dynamics at the border of the separation of the domain D3
and D4 of the Fig.3(e) we present in Fig.5(a) the variation of the phase of
all oscillators of the slave layer for $C_{2}=2$, $\epsilon_{1}=0.065721$,
$\epsilon_{2}=10\epsilon_{1}$ and interlayer coupling ($C_{0}$) varying from
17.8 to 18.2. At this range of values of the interlayer coupling, the slave
layer presents different dynamics. Before $C_{0}=18$ this slave layer shows a
synchronization of all oscillators of the slave layer. After this
synchronisation there follows a slight zone of disturbance (the zoom is given
at Fig.5(b)) before the division into two groups which drives the layer
towards the divergence. This abrupt change presented reminds of an explosive
desynchronization Lauter, Mitra, and Marquardt (2017). Given the form of the
connection (unidirectionnal coupling) between the first and second layer, the
oscillators of the slave layer can only remain synchronous and stable for a
certain range of values of $C_{0}$. This dynamic leads the slave layer to the
divergence that we observe in Fig.5(a) on the mean phase and in Fig.5(c) with
the order parameter. To understand more clearly the dynamics of the slave
layer at this value of interlayer coupling where we have destruction of the
synchronization, we show in Fig.5(d,e) the mean phase of the multilayer
network and the correlation between the oscillators of the slave layer
respectively. The mean phase presents a phase synchronization of all the
oscillators of the first layer but in the slave layer we have two clusters
formations. This cluster formation in the slave layer is confirmed using the
correlation between the oscillators of the slave layer and then we can
appreciate the formation of the these two clusters by the yellow colour.

##  IV Applications to other systems

The behavior of a multilayer network shown in the previous Sections is not
restricted to a system of Rössler oscillators, it can also be obtained with
other systems and topologies. In this Section we shall study a network of jerk
oscillators and another of Liénard oscillators with different topologies and
show that they reproduce the same behaviors.

###  IV.1 Synchronization with amplification in a multilayer network of jerk
oscillators

Here we investigate the dynamics of a multilayer network of jerk111Based on
Ref.Louodop  _et al._ (2017), a system is considered as jerk if the flow can
be rewritten as a third order differential equation in a single scalar
variable. For an isolated system, the jerk system can be defined as:
$\dddot{x}_{2}=\alpha\left(-\beta\left(-\dot{x}_{1}+\dot{x}_{2}-\gamma\ddot{x}_{1}+\frac{\gamma}{\alpha}\ddot{x}_{2}\right)+\dddot{x}_{1}\right)$.
where $\ddot{x}_{1}=I^{\prime}(x_{2}))\dot{x}_{2}$ and
$\dddot{x}_{1}=I^{\prime\prime}(x_{2}))\dot{x}_{2}^{2}+I^{\prime}(x_{2}))\ddot{x}_{2}.$
Then, $\dddot{x}_{1}$ is called jerk function. chaotic oscillators where the
first layer is described by Eq.7.

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{x}_{i}^{1}=I(x_{i}^{2})+{\epsilon_{1}}(x_{i+1}^{1}+x_{i-1}^{1}-2x_{i}^{1}),\\\
\dot{x}_{i}^{3}=\alpha(-x_{i}^{3}+I(x_{i}^{2})),\\\
\dot{x}_{i}^{3}=\beta(-x_{i}^{2}+x_{i}^{2}-\gamma
x_{i}^{3}).\end{array}\right.\end{array}$ (7)

The slave layer is described by Eq.8

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{y}_{j}^{1}=\frac{I(x_{j}^{2})}{C_{2}}+{\epsilon_{2}}(y_{j+1}^{1}+y_{j-1}^{1}-2y_{j}^{1})+C_{0}({x_{j}^{1}}-{C_{2}}{y_{j}^{1}}),\\\
\dot{y}_{j}^{2}=\alpha(-y_{j}^{3}+I(x_{j}^{2})/C_{2}),\\\
\dot{y}_{j}^{3}=\beta(-y_{j}^{1}+y_{j}^{2}-\gamma
y_{j}^{3}).\end{array}\right.\end{array}$ (8)

Where the piecewise linear function is

$I(x^{2})=\left\\{\begin{array}[]{ccc}-x^{2}\quad{if}\quad x^{2}\leq 1,\\\
-1\quad{otherwise}.\end{array}\right.$ (9)

$\alpha=0.025$, $\beta=0.765$, $\gamma=0.0938$ are the systems parameter.
$\epsilon_{1}$ and $\epsilon_{2}$ are the intralayer coupling strength of the
drive and slave layer of the network. The connection between the nodes of the
same layer is bidirectional and they are arranged on a ring. Moreover,
connection between the nodes of different layers is unidirectional and it only
concerns the oscillators with the same index in both layers.

Figure 6: Dynamics of the two layers of the network of jerk oscillators, where
the master layer is represented in red and the slave layer in blue. For
$\epsilon_{1}=\epsilon_{2}=2,C_{0}=5$ and $C_{2}=0.5$: (a) mean phase of the
first and second layer, (b) attractors of the oscillators in the master and
slave layer, (c) time series of the oscillators of the master and slave layer
and (d) synchronization between the oscillators of the first and second layer
in the network. For $\epsilon_{1}=\epsilon_{2}=2,C_{0}=5$ and $C_{2}=2$: (e)
mean phase of the first and second layer, (f) attractors of the oscillators in
the master and slave layer, (g) time series of the oscillators of the first
and second layer and (h) synchronization between the oscillators of the first
and second layer in the network.

Eq.7 and Eq.8 were solved numerically considering $N=30$ jerk oscillators per
layer with the systems parameter defined above. Fig.6 is obtained for
$\epsilon_{1}=2$, $\epsilon_{2}=2$ and $C_{0}=2$. Two values for the
amplification coefficient are considered: $C_{2}=0.5$ and $C_{2}=2$. For
$C_{2}=0.5$ we can clearly appreciate in Fig.6(a) the complete phase
synchronization of all oscillators in both layers. This phase synchronization
is a major condition to obtain amplification in the systems of the network. As
defined by our model, it emerges that for this values of $C_{2}$ the
oscillators of the slave network are supposed to be amplified compared to
oscillators of the master layer (see Fig.6(b,c) where in red we show the
oscillators of the master layer and in blue those of the slave layer). The
amplification between master and slave layer is perfectly observed in
Fig.6(d). In the second case we consider $C_{2}=2$. As in the previous case we
have phase synchronization between the oscillators of both layers (see
Fig.6(e)). For this value we have an amplification in the oscillators of the
master layer (see Fig.6(f,g)). In the same vein we present in Fig.6(h) the
synchronization between the oscillators in the first and second layer of the
network.

Therefore, interlayer synchronization can be obtained with amplification or
reduction depending on the value of the coefficient $C_{2}$.

###  IV.2 Chimera states with amplification in a multilayer network of Liénard
system

Let us consider a network of Liénard systems expressed as in reference Mishra
_et al._ (2015) where the authors chose to investigate the dynamics of the
oscillators basing themselves on an attractive and repulsive global coupling.
In the same vein, in the Refs.Hong and Strogatz (2012); Sathiyadevi  _et al._
(2018) the authors present some behaviour such as clusters, pattern formation,
synchronization and so on according to the attractive and repulsive coupling.

In this subsection, we consider the Liénard model with a intralayer topology
defined as in Mishra  _et al._ (2015).
First or master layer

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{x}_{i}^{1}=x_{i}^{2},\\\
\dot{x}_{i}^{2}=-\alpha x_{i}^{1}x_{i}^{2}-\beta(x_{i}^{1})^{3}-\gamma
x_{i}^{1}+\displaystyle
K\left[\left(\overline{x}^{2}-x_{i}^{2}\right)+\epsilon_{1}\left(\overline{x}^{1}-x_{i}^{1}\right)\right].\end{array}\right.\end{array}$
(10)

Where $\overline{x}^{1}=\frac{1}{N}\sum\limits_{i=1}^{N}x_{i}^{1}$ and
$\overline{x}^{2}=\frac{1}{N}\sum\limits_{i=1}^{N}x_{i}^{2}$.
Second or slave layer

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{y}_{j}^{1}=y_{j}^{2}+C_{0}({x_{i}^{1}}-{C_{2}}{y_{j}^{1}}),\\\
\dot{y}_{j}^{2}=-\alpha\frac{x_{j}^{1}x_{j}^{2}}{C_{2}}-\beta\frac{(x_{i}^{1})^{3}}{C_{2}}-\gamma
y_{j}^{1}+\displaystyle
K\left[\left(\overline{y}^{2}-y_{j}^{2}\right)+\epsilon_{2}\left(\overline{y}^{1}-y_{j}^{1}\right)\right]\\\
+C_{0}({x_{j}^{2}}-{C_{2}}{y_{j}^{2}}).\end{array}\right.\end{array}$ (11)

With $\overline{y}^{1}=\frac{1}{N}\sum\limits_{i=1}^{N}y_{i}^{1}$ and
$\overline{y}^{2}=\frac{1}{N}\sum\limits_{i=1}^{N}y_{i}^{2}$.

The system’s parameters $\alpha$, $\beta$ and $\gamma$ are selected exactly as
in Mishra  _et al._ (2015), $K$ is the strength of coupling, $\epsilon_{1}$
and $\epsilon_{2}$ are the intralayer global mean field coupling and $C_{2}$
is the amplification coefficient. The coupling between the nodes of different
layers is unidirectional and it only concerns the oscillators with the same
index in both layers and $N=100$ Liénard systems.

For the numerical simulation, we consider two cases, one where the
amplifcation is less than one (the systems of the slave layer are amplified
and the systems of master layer are reduced) and another when the
amplification is greather than one (the systems of the slave layer are reduced
and the systems of master layer are amplified).

Figure 7: Examples of behaviour of the multilayer network of Liénard systems.
Amplification for $K=0.9$, $\epsilon_{1}=\epsilon_{2}=-0.57$, $C_{2}=0.5$: (a)
Snapshot of all the oscillators of the multilayer network showing a
multichimera state, (b) phase portrait of the oscillators $j=1$ and $j=30$ of
both layers and (c,d) temporal dynamics of all the nodes in the first and
second layer in the network. Reduction for $K=0.9$,
$\epsilon_{1}=\epsilon_{2}=-0.57$, $C_{2}=2$: (e) Snapshot of all the
oscillators of the multilayer network showing a multichimera state, (f) phase
portrait of the oscillators $j=1$ and $j=30$ of both layers and (g,h) temporal
dynamics of all the nodes in the first and second layer in the network.

Fig. 7 elaborates the Chimera-I states based on Refs.Mishra  _et al._ (2015);
Dudkowski, Maistrenko, and Kapitaniak (2014). Fig. 7(a) presents a snapshot of
all the oscillators of both layers of the network (the first 100 systems
correspond to the master layer and the rest is for slave layer). This figure
presents a multichimera state in both master and slave layer (see Fig. 7(c,d))
and a phase synchronization of both layers for $K=0.9$,
$\epsilon_{1}=\epsilon_{2}=-0.57$, and $C_{2}=0.5$. For the same parameters we
show in Fig.7(b) the attractor of the system for oscillator $j=1$ in the
synchronization state and the attractor for $j=30$ in the incoherent state for
both layers of the network (attractor red corresponds to the master layer and
attractor blue for the slave layer). At $C_{2}=0.5$ the systems of the slave
layer are supposed to be amplified as in Fig.7(b). In the same vein we present
in Fig.7(e,f,g,h) for $C_{2}=2$ the same behaviours as in Fig.7(a,b,c,d). We
observe the same behaviours except that are the systems of the master layer
which are amplified as $X=2Y$.
Based on these results, we conclude that this form of coupling can lead to a
chimera or multichimera state with amplification or reduction depending on the
value of the amplification coeficient.

Figure 8: Examples of behaviour of the multilayer network of Liénard systems.
Amplification for $K=1.5$, $\epsilon_{1}=\epsilon_{2}=-0.1$, $C_{2}=0.5$: (a)
Snapshot of all the oscillators of the multilayer network showing a cluster,
(b) phase portrait of the oscillators $j=1$, $j=50$ and $j=90$ of both layers
and (c,d) temporal dynamics of all the nodes in the first and second layer in
the network. Reduction for $K=1.5$, $\epsilon_{1}=\epsilon_{2}=-0.17$,
$C_{2}=2$: (e) Snapshot of all the oscillators of the multilayer network
showing a cluster, (f) phase portrait of the oscillators $j=1$, $j=50$ and
$j=70$ of both layers and (g,h) temporal dynamics of all the nodes in the
first and second layer in the network.

Let us now consider the following parameters $K=1.5$
$\epsilon_{1}=\epsilon_{2}=-0.1$ that lead to a cluster formation and we
choose two values of amplification parameter: (a) $C_{2}=0.5$ (see
Fig.8(a,b,c,d)) and (b) $C_{2}=2$ (see Fig.7(e,f,g,h)). According to these two
cases we notice that this cluster state can be maintained with an
amplification or reduction based on the value of $C_{2}$.

##  V Conclusion

Summarizing, we have studied and characterized numerically the synchronization
and the amplification of signals in a multilayer network of Rössler, jerk or
Liénard oscillators. Using tools for studying synchronization in the network
such as master stability function, the order parameter, we have demonstrated
that, the existence of synchronization in the second layer is conditioned by
the first. To obtain synchronization of both layers at the same value of the
intralayer coupling, the amplification coefficient must be sufficiently low
(see Fig.2(b)). The key role of amplification is demonstrated by analyzing the
order parameter of the first and second layer. This parameter leads the
network to different dynamics such as cluster formation and synchronization.

From a theoretical point of view, this work contributes to the advancement in
the understanding of the phenomenon of synchronization and amplification of
signals between two coupled networks. From a practical point of view, the
results may be useful in many technological applications. For example, in
several mechanical systems, the transmission of movements or orders is done
through a driving belt or gear Denijs and Lindner (1997). Concerning this
process, one of the most important parameters is the transmission ratio, here
represented by $C_{2}$, the amplification parameter. The amplification
coefficient could be the transmission ratio between gears or pulley-belt
systems Farshidianfar and Saghafi (2014).

##  VI Data Availability Statement

The data that supports the findings of this study are available within this
article.

## Acknowledgements

HAC thanks ICTP-SAIFR and FAPESP grant 2016/01343-7 for partial support. VC
and FFF acknowledge that this study was financed in part by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
001.

## Appendix A

Synchronization conditions

We consider $X(x_{i}^{1},x_{i}^{2},x_{i}^{3})$ and
$Y(y_{i}^{1},y_{i}^{2},y_{i}^{3})$ the state vector of the $i^{th}$ oscillator
of the first and second layer respectively.
The dynamics of the first and second layers are given respectively by:

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{x}_{i}^{1}=-x_{i}^{2}-x_{i}^{3}+{\epsilon_{1}}\sum\limits_{k=1}^{N}A_{ik}^{1}(x_{k}^{1}-x_{i}^{1}),\\\
\dot{x}_{i}^{2}=x_{i}^{1}+ax_{i}^{2},\\\
\dot{x}_{i}^{3}=bx_{i}^{1}+x_{i}^{3}(x_{i}^{1}-c).\end{array}\right.\end{array}$
(12)

and

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}{{\dot{y}}_{j}^{1}}=-{y_{j}^{2}}-{y_{j}^{3}}+{\epsilon_{2}}\sum\limits_{k=1}^{N}A_{jk}^{2}({y_{k}^{1}}-{y_{j}^{1}})+C_{0}(x_{j}^{1}-C_{2}y_{j}^{1}),\\\
{{\dot{y}}_{j}^{2}}=\frac{{x_{j}^{1}}+a{x_{j}^{2}}}{C_{2}},\\\
{{\dot{y}}_{j}^{3}}=b{y_{j}^{1}}+\frac{{{x_{j}^{3}}{x_{j}^{1}}}}{{{C_{2}}}}-c{y_{j}^{3}}.\end{array}\right.\end{array}$
(13)

Without any intralayer coupling $(\epsilon_{1}=\epsilon_{2}=0)$ both layers
could synchronize with amplification depending on $C_{2}$.
We consider the synchronization error $e=X-C_{2}Y$ with
$\epsilon_{1}=\epsilon_{2}=0$. So, the error dynamical system which is
obtained only for the couple of systems with the same index in the first and
second layers is given by Eq.14:

$\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{e}_{i}^{1}=-e_{i}^{2}-e_{i}^{3}-C_{2}C_{0}e_{i}^{1},\\\
\dot{e}_{i}^{2}=0,\\\
\dot{e}_{i}^{3}=be_{i}^{1}-ce_{i}^{3}.\end{array}\right.\end{array}$ (14)

To simplify the demonstration, we consider the Lyapunov functions
$v_{i}(i=1,2,...,N)$ of system pairs $i$ and $j$ (with $i=j$) of the first and
second layer described as follow:

$\begin{array}[]{cc}v_{i}=\frac{1}{2}\left((e_{i}^{1})^{2}+(e_{i}^{2})^{2}+\frac{1}{b}(e_{i}^{3})^{2}\right)\end{array}$
(15)

For any couple $(X,Y)$ we choose the following Lyapunov function candidate:

$\begin{array}[]{cc}v=\frac{1}{2}\left((e^{1})^{2}+(e^{2})^{2}+\frac{1}{b}(e^{3})^{2}\right).\end{array}$
(16)

It is established (assuming $a>0,b>0$ and $c>0$) that the system defined by
Eq.14 is practically stable since the time derivative of the Lyapunov function
in Eq.16 is bounded by a positive constant. This also means that the error
between the driver and response systems is sufficiently small, but different
from zero and could be considered as a tolerance in the synchronization
condition Kakmeni  _et al._ (2010).

$\begin{array}[]{cc}\dot{v}\leq\frac{(e^{2})^{2}}{4c}.\end{array}$ (17)

For the whole network, the Lyapunov function candidate $V$ can be defined as a
sum of $v_{i}$:

$\begin{array}[]{cc}V=\sum\limits_{i=1}^{N}v_{i}.\end{array}$ (18)

And then for the whole network we can have:

$\begin{array}[]{cc}\dot{V}\leq\sum\limits_{i=1}^{N}\frac{(e_{i}^{2})^{2}}{4c}.\end{array}$
(19)

Based on Eq.14, this boundedness is ensured by the fact that $e_{i}^{2}(t)$ is
constant due to the fact that $\dot{e}_{i}^{2}(t)=0$. Thereby, from Eq.19
$e\rightarrow 0$, $X-C_{2}Y=0$ and induce $X=C_{2}Y$. we can obtain
amplification or reduction depending on the value of the coefficient $C_{2}$.

## Appendix B

Preliminaries to the Investigation of the Dynamics of the Network: the Master
Stability Function

A regular problem that arises when analyzing the dynamics of a network is to
find conditions that guarantee the synchronization of a system of coupled
identical nonlinear oscillators, so that all the oscillators converge
asymptotically towards the same state.

The Master Stability Function (MSF) developed by Pecora and Carroll Pecora and
Carroll (1998), constitutes one of the most useful tools to analyze the
synchronization stability of a system of coupled identical nonlinear
oscillators Leyva  _et al._ (2018); Pecora and Carroll (1998); Majhi, Perc,
and Ghosh (2016). We develop here only some points of the principal idea of
this method.

Considering a network of $N$ identical coupled chaotic oscillators (or nodes),
let $\mathbf{x}_{i}$ a vector with $m$ components necessary to describe the
state of the $i^{th}$ node. In general, in absence of any interaction between
the nodes of the network, the evolution of a node is given by Eq.20:

$\displaystyle\mathbf{\dot{x}}_{i}=\mathbf{F(x}_{i}),$ (20)

In this Eq.20, $\mathbf{F}$ is a function defined from $\mathbb{R}^{m}$ to
$\mathbb{R}^{m}$ and is used to define the local dynamics of the oscillators.
To describe how the oscillators evolve when they are connected in a network,
we need to consider not only the local dynamics presented at Eq.20, but also
how each node is affected by the ones to which it is connected. So, the law
governing the dynamical interaction of the $i^{th}$ node is defined as:

$\displaystyle\mathbf{\dot{x}}_{i}=\mathbf{F(x}_{i})+\sigma\sum\limits_{j=1}^{N}\mathbf{G}_{ij}\mathbf{H(x}_{j}),$
(21)

where $\sigma$ is a coupling strength,
$H:\mathbb{R}^{m}\longrightarrow\mathbb{R}^{m}$ is an arbitrary output
function of each node’s variables using in the coupling. If we put the network
in a synchronized state, we have $\mathbf{x}_{i}=\mathbf{s}$ for all nodes,
where $\mathbf{s}$ is any m-dimensional vector. The only way all nodes have
the same behavior is to have the sum, $\sum\limits_{j=1}^{N}\mathbf{G}_{ij}$,
be the same for all $i$. So, to obtain complete or identical synchronization
the row sums of the coupling matrix must be the same for all rows. According
to Pecora and Carrol Pecora and Carroll (1998), we can collect the node
dynamical variables, functions and coupling in:

$\displaystyle\mathbf{x}$
$\displaystyle=\left[\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{N}\right],$
(22a) $\displaystyle\mathbf{F(x)}$
$\displaystyle=\left[\mathbf{F}(\mathbf{x}_{1}),\mathbf{F}(\mathbf{x}_{2}),\ldots,\mathbf{F}(\mathbf{x}_{N})\right],$
(22b) $\displaystyle\mathbf{H(x)}$
$\displaystyle=\left[\mathbf{H}(\mathbf{x}_{1}),\mathbf{H}(\mathbf{x}_{2}),\ldots,\mathbf{H}(\mathbf{x}_{N})\right],$
(22c)

and $G$ be the matrix coupling coefficients $G_{ij}$, then based on Eq.22(a,b
and c), Eq.21 can be written in a compact form as follow:

$\displaystyle\mathbf{\dot{x}}=\mathbf{F(x)}+\sigma\mathbf{G}\otimes\mathbf{H(x)},$
(23)

where $\otimes$ is the Kronecker product. So, according to the form of Eq.21,
we can rewrite Eq.4 in the same form with:

$\mathbf{F(x}_{i})\begin{array}[]{cc}\left\\{\begin{array}[]{l}\dot{x}_{i}^{1}=-x_{i}^{2}-x_{i}^{3},\\\
\dot{x}_{i}^{2}=x_{i}^{1}+ax_{i}^{2},\\\
\dot{x}_{i}^{3}=bx_{i}^{1}+x_{i}^{3}(x_{i}^{1}-c).\end{array}\right.\end{array}$
(24)

According to Pecora and Carroll (1998), for an all-to-all coupling scheme, the
connectivity matrix can be defined as in matrix $G$. To couple the nodes of
the layer we choose the $\mathbf{x}^{1}$ component and then the matrix $H$ can
be defined as in Eq.25:

$\displaystyle\mathbf{H}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix},\;\;\mathbf{G}=\begin{pmatrix}1-N&1&\cdots&1&1&1\\\
1&1-N&\cdots&1&1&1\\\ \vdots&\vdots&\cdots&\vdots&\vdots&\vdots\\\
1&1&\cdots&1&1-N&1\\\ 1&1&\cdots&1&1&1-N\end{pmatrix},$ (25)

The Master Stability Function studies the stability of the global
synchronization in the network. Therefore, the synchronous state is obtained
when $\mathbf{x}_{1}=\mathbf{x}_{2}=\ldots=\mathbf{x}_{N}=\mathbf{s}$.

Suppose our system is synchronized and we perturb it so that each node is, in
general, slightly “away” from the synchronized motion. Let us consider
$\xi_{i}$ a small perturbation of the $i^{th}$ node of the network. so that
after the perturbation $\mathbf{x}_{i}=\mathbf{s}+\xi_{i}$. For N oscillators
of the first layer the collections of the variations can be expressed as
$\xi=(\xi_{1},\xi_{2},...,\xi_{N})$. Now we can derive an equation of motion
for the small perturbations that we will use to explore if the synchronized
state is unstable or stable. So, replacing the perturbation
$\mathbf{x}_{i}=\mathbf{s}+\xi_{i}$ in Eq.21 and using Taylor theorem expand
of $\mathbf{F(s+\xi}_{i})$ and $\mathbf{H(s+\xi}_{i})$ to first order (since
$\xi_{i}$ is small) we have the following variational equation:

$\displaystyle\mathbf{\dot{\xi}}_{i}=\mathbf{DF(s)\xi}_{i}+\sigma\sum\limits_{j=1}^{N}\mathbf{G}_{ij}\mathbf{DH(s)\xi}_{j},$
(26)

Using tensor notation, we can write Eq.26 in a more compact form:

$\displaystyle\mathbf{\dot{\xi}}=[\mathbf{1_{N}\otimes
D\textbf{F(s)}}+\sigma\mathbf{G}\otimes\mathbf{DH(s)}]\xi,$ (27)

where $\mathbf{1_{N}}$ is the identity matrix of order $N$, DF and DH are the
$N\times N$ Jacobian matrices of the corresponding vector functions.

The solution of Eq.27 can be in the form $\xi_{i}\sim\exp\mu_{i}t$. The
exponents $\mu$ tell us if the perturbation grows ($\mu>0$) or shrinks
($\mu<0$), the former indicating a direction that is unstable and the latter a
stable direction. After diagonalization of the second term of Eq.27 we obtain
the variational equations which are diagonal in the node coordinates and are
now uncoupled and individually given by:

$\displaystyle\mathbf{\dot{\xi}_{k}}=[\textbf{DF(s)}+\sigma\alpha_{k}D\mathbf{H(s)}]\xi_{k}.$
(28)

where $\alpha_{k}$ is an eigenvalue of G, $k=1,2,...N$. For each $k$, the form
of each block of the Eq.27 does not change, only the scalar multiplier
$\sigma\alpha_{k}$ differs for each block.

Therefore, these steps lead us to design the following master stability
equation:

$\displaystyle\mathbf{\dot{\xi}}=[\textbf{DF(s)}+\sigma\alpha\mathbf{DH(s)}]\xi.$
(29)

Computing the Largest Lyapunov exponent of this Master Stability Equation
Eq.29 we obtain what Pecora and Carroll called The Master Stability Function
and therefore we achieve a stable synchronization state if the MSF turns
negative Pecora and Carroll (1998); Tang  _et al._ (2019); Huang  _et al._
(2009).

## Appendix C

Calculation of the Order Parameter

Collective behavior of such an N-oscillator system is conveniently described
by the order parameter. The evaluation of this order parameter Kuramoto and
Battogtokh (2002) used the phase of each oscillator of the network. To define
the phase let us consider an arbitrary signal $s(\tau)$ with time $\tau$ and
its Hilbert transformation to be $\tilde{s}(\tau)$, we have:

$\psi(\tau)=s(\tau)+i\tilde{s}(\tau)=R(\tau)\exp^{i\phi(\tau)},$ (30)

where $R(\tau)$ is the amplitude and $\phi(\tau)$ the phase of the variable
$s(\tau)$. If the instantaneous phase is $\phi_{i}(\tau)$, it can be
determined through the following relation:

$\phi_{i}(\tau)=\tan^{-1}\left[\displaystyle\frac{\tilde{s}_{i}(\tau)}{s_{i}(\tau)}\right].$
(31)

In this paper the calculation of the phase was carried out using in each case
the variable that best describes the dynamics of the system.
Thus, from the expression of the phase $\phi_{i}$, the mean phase $\phi$ is an
algebraic average calculated on the $N$ oscillators of the layer. So, for a
network of N oscillators the order parameter can be expressed as:

$\begin{array}[]{cc}r=\frac{1}{N}\sum\limits_{i=1}^{N}e^{j\phi_{i}}\end{array}$
(32)

where $j^{2}=-1$, when $r\rightarrow 1$, phase synchronization is reached and
when $r\approx 0$, the network is desynchronized.

## References

## References

  * Leyva  _et al._ (2018) I. Leyva, I. Sendiña-Nadal, R. Sevilla-Escoboza, V. Vera-Avila, P. Chholak, and S. Boccaletti, “Relay synchronization in multiplex networks,” Sci. Rep. 8, 1 (2018).
  * De Domenico  _et al._ (2013) M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, “Mathematical formulation of multilayer networks,” Phys. Rev. X 3, 041022 (2013).
  * Zhang  _et al._ (2015) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, “Explosive synchronization in adaptive and multilayer networks,” Phys. Rev. Lett. 114, 038701 (2015).
  * Jalan, Kachhvah, and Jeong (2020) S. Jalan, A. D. Kachhvah, and H. Jeong, “Explosive synchronization in multilayer dynamically dissimilar networks,” Journal of Computational Science , 101177 (2020).
  * Belykh, Carter, and Jeter (2019) I. Belykh, D. Carter, and R. Jeter, “Synchronization in multilayer networks: when good links go bad,” SIAM Journal on Applied Dynamical Systems 18, 2267 (2019).
  * 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,” Phys. Rev. Lett. 122, 014101 (2019).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Phys. Rep. 469, 93 (2008).
  * del Genio  _et al._ (2016) C. del Genio, J. Gómez-Gardeñes, I. Bonamassa, and S. Boccaletti, “Synchronization in networks with multiple action layers,” Sci. Adv. 2, e1601679 (2016).
  * Zhang  _et al._ (2001) Y. Zhang, G. Hu, H. A. Cerdeira, S. Chen, T. Braun, and Y. Yao, “Partial synchronization and spontaneous spatial ordering in coupled chaotic systems,” Phys. Rev. E 63, 026211 (2001).
  * Ji  _et al._ (2014) P. Ji, T. K. D. Peron, F. A. Rodrigues, and J. Kurths, “Analysis of cluster explosive synchronization in complex networks,” Phys. Rev. E 90, 062810 (2014).
  * Lauter, Mitra, and Marquardt (2017) R. Lauter, A. Mitra, and F. Marquardt, “From kardar-parisi-zhang scaling to explosive desynchronization in arrays of limit-cycle oscillators,” Phys. Rev. E 96, 012220 (2017).
  * Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, “Coexistence of coherence and incoherence in nonlocally coupled phase oscillators,” Nonlinear Phenom. Complex Syst. 5, 380 (2002).
  * Schmidt  _et al._ (2014) L. Schmidt, K. Schönleber, K. Krischer, and V. García-Morales, “Coexistence of synchrony and incoherence in oscillatory media under nonlinear global coupling,” Chaos 24, 013102 (2014).
  * Schmidt and Krischer (2015) L. Schmidt and K. Krischer, “Clustering as a prerequisite for chimera states in globally coupled systems,” Phys. Rev. Lett. 114, 034101 (2015).
  * Mishra  _et al._ (2015) A. Mishra, C. Hens, M. Bose, P. K. Roy, and S. K. Dana, “Chimeralike states in a network of oscillators under attractive and repulsive global coupling,” Phys. Rev. E 92, 062920 (2015).
  * Dudkowski, Maistrenko, and Kapitaniak (2014) D. Dudkowski, Y. Maistrenko, and T. Kapitaniak, “Different types of chimera states: An interplay between spatial and dynamical chaos,” Phys. Rev. E 90, 032920 (2014).
  * Castellano, Fortunato, and Loreto (2009) C. Castellano, S. Fortunato, and V. Loreto, “Statistical physics of social dynamics,” Rev. Mod. Phys. 81, 591 (2009).
  * Pastor-Satorras  _et al._ (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, “Epidemic processes in complex networks,” Rev. Mod. Phys. 87, 925 (2015).
  * Liu, Slotine, and Barabási (2011) Y.-Y. Liu, J.-J. Slotine, and A.-L. Barabási, “Controllability of complex networks,” nature 473, 167 (2011).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett. 80, 2109 (1998).
  * 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,” Phys. Rep. 544, 1 (2014).
  * Mondragon and Bianconi (2018) J. I. Mondragon, Raul J. and G. Bianconi, “Multilink communities of multiplex networks,” PloS one 13, e0193821 (2018).
  * Boccaletti  _et al._ (2016) S. Boccaletti, R. Criado, M. Romance, and J. J. Torres, “Introduction to focus issue: Complex dynamics in networks, multilayered structures and systems,” Chaos 26, 065101 (2016).
  * (24) I. Leyva, R. Sevilla-Escoboza, I. Sendina-Nadal, R. Gutierrez, J. Buldu, and S. Boccaletti, “Inter-layer synchronization in nonidentical multi-layer networks,” .
  * Ghosh, Zakharova, and Jalan (2018) S. Ghosh, A. Zakharova, and S. Jalan, “Non-identical multiplexing promotes chimera states,” Chaos, Solitons & Fractals 106, 56 (2018).
  * Rössler (1979) O. E. Rössler, “Continuous chaos—four prototype equations,” Annals of the New York Academy of Sciences 316, 376 (1979).
  * Louodop  _et al._ (2017) P. Louodop, S. Saha, R. Tchitnga, P. Muruganandam, S. K. Dana, and H. A. Cerdeira, “Coherent motion of chaotic attractors,” Phys. Rev. E 96, 042210 (2017).
  * Zou  _et al._ (2015) W. Zou, D. Senthilkumar, R. Nagao, I. Z. Kiss, Y. Tang, A. Koseska, J. Duan, and J. Kurths, “Restoration of rhythmicity in diffusively coupled dynamical networks,” Nature communications 6, 1 (2015).
  * Ghosh, Banerjee, and Kurths (2015) D. Ghosh, T. Banerjee, and J. Kurths, “Revival of oscillation from mean-field-induced death: Theory and experiment,” Physical Review E 92, 052908 (2015).
  * Karnatak, Ramaswamy, and Prasad (2007) R. Karnatak, R. Ramaswamy, and A. Prasad, “Amplitude death in the absence of time delays in identical coupled oscillators,” Physical Review E 76, 035201 (2007).
  * Cafieri  _et al._ (2018) S. Cafieri, L. Cellier, F. Messine, and R. Omheni, “Combination of optimal control approaches for aircraft conflict avoidance via velocity regulation,” Optimal Control Applications and Methods 39, 181 (2018).
  * Chen  _et al._ (2017) H.-B. Chen, Y.-T. Sun, J. Gao, C. Xu, and Z.-G. Zheng, “Order parameter analysis of synchronization transitions on star networks,” Frontiers of Physics 12, 120504 (2017).
  * Hong and Choi (2000) H. Hong and M. Choi, “Phase synchronization and noise-induced resonance in systems of coupled oscillators,” Physical Review E 62, 6462 (2000).
  * Hu  _et al._ (2000) G. Hu, Y. Zhang, H. A. Cerdeira, and S. Chen, “From low-dimensional synchronous chaos to high-dimensional desynchronous spatiotemporal chaos in coupled systems,” Phys. Rev. Lett. 85, 3377 (2000).
  * Maistrenko  _et al._ (2004) Y. Maistrenko, O. Popovych, O. Burylko, and P. A. Tass, “Mechanism of desynchronization in the finite-dimensional kuramoto model,” Phys. Rev. lett. 93, 084102 (2004).
  * Allefeld, Müller, and Kurths (2007) C. Allefeld, M. Müller, and J. Kurths, “Eigenvalue decomposition as a generalized synchronization cluster analysis,” International Journal of Bifurcation and Chaos 17, 3493 (2007).
  * Singha  _et al._ (2017) J. Singha, V. Tchuyang, P. Louodop, R. Tchitnga, H. A. Cerdeira, and N. Gupte, “Spatial splay states in coupled map lattices and josephson junction arrays,” Indian Academy of Sciences 1, 3493 (2017).
  * Nichols and Wiesenfeld (1992) S. Nichols and K. Wiesenfeld, “Ubiquitous neutral stability of splay-phase states,” Phys. Rev. A 45, 8430 (1992).
  * Fujisaka and Yamada (1983) H. Fujisaka and T. Yamada, “Stability theory of synchronized motion in coupled-oscillator systems,” Prog. Theor. Phys. 69, 32 (1983).
  * Pecora and Carroll (1990) L. M. Pecora and T. L. Carroll, “Synchronization in chaotic systems,” Phys. Rev. Lett. 64, 821 (1990).
  * Pearson (1920) K. Pearson, “Notes on the history of correlation,” Biometrika 13, 25 (1920).
  * Hong and Strogatz (2012) H. Hong and S. H. Strogatz, “Mean-field behavior in coupled oscillators with attractive and repulsive interactions,” Phys. Rev. E 85, 056210 (2012).
  * Sathiyadevi  _et al._ (2018) K. Sathiyadevi, V. Chandrasekar, D. Senthilkumar, and M. Lakshmanan, “Distinct collective states due to trade-off between attractive and repulsive couplings,” Phys. Rev. E 97, 032207 (2018).
  * Denijs and Lindner (1997) J. F. Denijs and M. Lindner, “Transmission for converting rotary motion into linear motion,” (1997), uS Patent 5,690,567.
  * Farshidianfar and Saghafi (2014) A. Farshidianfar and A. Saghafi, “Identification and control of chaos in nonlinear gear dynamic systems using melnikov analysis,” Physics Letters A 378, 3457–3463 (2014).
  * Kakmeni  _et al._ (2010) F. M. Kakmeni, S. Bowong, D. Senthikumar, and J. Kurths, “Practical time-delay synchronization of a periodically modulated self-excited oscillators with uncertainties,” Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 043121 (2010).
  * Majhi, Perc, and Ghosh (2016) S. Majhi, M. Perc, and D. Ghosh, “Chimera states in uncoupled neurons induced by a multilayer structure,” Sci. Rep. 6, 39033 (2016).
  * Tang  _et al._ (2019) L. Tang, X. Wu, J. Lü, J.-a. Lu, and R. M. D’Souza, “Master stability functions for complete, intralayer, and interlayer synchronization in multiplex networks of coupled rössler oscillators,” Phys. Rev. E 99, 012304 (2019).
  * Huang  _et al._ (2009) L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, “Generic behavior of master-stability functions in coupled nonlinear dynamical systems,” Phys. Rev. E 80, 036204 (2009).