Eigenvalue spectrum and synchronizability of multiplex chain networks

Physica A 537 (2020) 122631

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Eigenvalue spectrum and synchronizability of multiplex chain networks ∗

Yang Deng, Zhen Jia , Guangming Deng, Qiongfen Zhang College of Science, Guilin University of Technology, Guilin 541004, China

article

info

Article history: Received 8 June 2019 Received in revised form 4 September 2019 Available online 14 September 2019 Keywords: Multiplex chain networks Synchronizability Eigenvalue spectrum

a b s t r a c t Synchronization phenomena are of broad interest across disciplines and increasingly of interest in a multiplex network setting. In this paper, the problem of synchronization of two multiplex chain networks is investigated, according to the master stability function method. We define two kinds of multiplex chain networks according to different coupling modes: one is a class of the multiplex chain networks with one-to-one undirected coupling between layers(Networks-A), and the other is a class of the multiplex chain networks with one-to-one unidirectional coupling between layers(Networks-B). The eigenvalue spectrum of the supra-Laplacian matrices of two kinds of the networks is strictly derived theoretically, and the relationships between the structural parameters and synchronizability of the networks are further revealed. The structural parameter values of the networks to achieve the optimal synchronizability are obtained. Numerical examples are also provided to verify the effectiveness of theoretical analysis. © 2019 Elsevier B.V. All rights reserved.

1. Introduction In recent years, the study of multiplex complex network has become a frontier and hot topic in network science. The existing researches on complex networks are plentiful [1–20], but most of them focus on the researches of single networks. With the deepening of researches and the arrival of the era of big data, people has realized that many real world networks do not exist in isolation, they are often interact with other networks of similar or different nature, forming what is known as Networks-of-Networks (NoN) [21–23]. Examples of such networks are widespread, for instance, in a social system, a set of individuals interact between each other in various modes of social interactions between the same people: an individual has interactions with others through online social systems (such as Facebook or Twitter) and off-line systems (such as professional or personal circles), which formed a two-layer network that is coupled between online and offline. In another example, there are many forms of interaction between digital rumors, such as blogs and emails, not only on their respective networks, but also across each other. Once the interactions on these networks are considered, the observations obtained in a single-layer network environment can change dramatically. Therefore, researches limited to single-layer networks are far from meeting the needs in realistic development, and the development of theories, methods and techniques in multiplex networks research is imminent. Synchronization in complex networks has been widely studied, very few works have investigated synchronization in NoNs. Gómez et al. developed a formalism to unveil the time scales of diffusive processes on multiplex networks [24], Wu et al. investigated generalized outer synchronization between two completely different complex dynamical networks with directional interlayer connection, but the interlayer couplings are directed and one-to-one [25,26]. Li et al. investigated ∗ Corresponding author. E-mail address: [email protected] (Z. Jia). https://doi.org/10.1016/j.physa.2019.122631 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

that, for duplex networks composed of two networks connecting by two links, the degree of connector nodes plays a crucial role in facilitating complete synchronization of the duplex [27]. The maximum eigenvalue of the Laplacian matrix of the coupled network has a decisive influence on the synchronization stability of the network [28,29]. In 2013, SoléRibalta et al. presented the asymptotic analysis of the spectrum of the Laplacian of multiplex networks and found analytical expressions that allow us to infer the behavior of dynamical processes [30]. In 2016, Xu Mingming et al. studied the eigenvalue spectrum and synchronizability of two-layer star network, and gave the analytic values of eigenvalues of twolayer star networks [31,32]. In 2017, Sun Juan et al. further explored the eigenvalue spectrum and synchronizability of multiplex unidirectional coupled star networks [33], Wei Juan et al. proposed the synchronizability of a two-layer rule network based on the master stability function method [34]. However, due to the complexity of multiplex networks, there are few rigorous theoretical deductions, most of the studies are based on the results of numerical simulation on synchronizability of two-layer networks. Based on this motivation, the aim of our work is to study the eigenvalue spectrum and synchronizability of more general multiplex chain networks (not limited to two layers), to strictly derive the eigenvalue spectrum of the supra-Laplacian matrices of two kinds of multiplex chain networks, and to analyze the relationships between the synchronizability and structural parameters. The theoretical results are verified by numerical simulation experiments. This paper is organized as follows. Section 2 introduces the dynamic equation, the supra-Laplacian matrices and the internal structure of the multiplex chain networks. Section 3 studies the eigenvalue spectrum of the multiplex chain networks with one-to-one undirected coupling between layers(Networks-A). Section 4 explores the eigenvalue spectrum of the multiplex chain networks with one-to-one unidirectional coupling between layers(Networks-B). Section 5 studies the synchronizability of two kinds of multiplex chain networks. Finally, concluding comments are given in Section 6. 2. Preliminaries 2.1. Multiplex network dynamics model We consider the multiplex network consisting of M layers each consisting of N nodes. The evolution of the full multiplex system can be written as [26,32] x˙ Ki = f xKi + a

( )

N ∑

M ∑ ( ) ( L) ωijK H xKj + d dKL i Γ xi , (i = 1, 2, . . . , N ; K = 1, 2, . . . , M )

j=1

(1)

L=1

∈ Rn is the state the ith node in the K th layer, and f (∗) : Rn → Rn is a well-defined vector function, → Rn and a are the inner coupling function and coupling strength for nodes within each layer, respectively, : Rn(→ )Rn and d are the inner coupling function and coupling strength for nodes across layers, respectively. = ωijK ∈ RN ×N is the coupling weight configuration matrix of the K th layer. Explicitly, if the ith node ∑N K is connected with the jth (j ̸ = i) node within the K th layer, ωijK = 1, otherwise ωijK = 0, and ωiiK = − j=1 ωij , for

where xKi H (∗) : Rn and Γ (∗) Here, W K

j̸ =i

i, j = 1, 2, . . . , N , K = 1, 2, . . . , M. Thus L(K ) = −aW K is a Laplacian matrix. When the connections between nodes across layers are undirected, if the ∑ith node in the K th layer is connected with KL KK its replica in the Lth (L ̸ = K ) layer, dKL = − ML=1 dKL i = 1, otherwise di = 0, and di i , for K , L = 1, 2, . . . , M. It is obvious L̸ =K

(

)

that D = dKL ∈ RM ×M is also a negative Laplacian matrix. i When the connections between nodes across layers are unidirectional, if the ith node in the K th layer is unidirectional KL connected with its replica in the Lth (L ̸ = K ) layer, dKL i = 1 (1 ≤ K ≤ M , (K <) L), or di = 1 (1 ≤ L ≤ M , K > L), otherwise ∑M KL KK KL KL M ×M di = 0, and di = − ∈R , and then D is an upper triangular L=1 di , for K , L = 1, 2, . . . , M. Similarly, D = di L̸ =K

matrix or a lower triangular matrix. Let L be the supra-Laplacian matrix of Eq. (1), LI be the supra-Laplacian matrix representing the interlayer topology, and LL be the supra-Laplacian matrix describing the intralayer topology. Then L can be written as: (2)

L = LI + LL

Taking LI to be the Laplacian matrix of the interlayer networks, we have L I = L I ⊗ IN

(3)

where ⊗ is the Kronecker product, LI = −dD, IN is the N × N identity matrix. As for LL , it can be represented by the direct sum of the Laplacian matrix L(K ) within each layer, namely, L(1) ⎜ 0 ⎜

⎛ LL = ⎜

⎝

.. .

0

0

L

(2)

··· ···

.. . 0

0 0

.. .

···

L(M )

⎞ ⎟ M ⎟ ⎟ = ⊕ L( K ) ⎠ K =1

(4)

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

3

Fig. 1. Schematic diagram of three-layer chain network structure, other layers can be seen as first layer replicas, and connecting each node with its counterpart in the other layer. (a) three-layer chain network structure with one-to-one undirected coupling between layers. (b) three-layer chain network structure with one-to-one unidirectional coupling between layers.

The eigenvalues of the supra-Laplacian matrix L are recorded as 0 =λ1 < λ2 ≤ λ3 ≤ · · · ≤ λmax . According to the master stability function method, the unbounded and bounded synchronous regions are more usually taken into consideration in real-world networks, the synchronizability of the network (1) is determined by the minimum non-zero eigenvalue λ2 or the ratio R = λmax /λ2 of the maximum eigenvalue to the minimum non-zero eigenvalue of supraLaplacian matrix L. Generally, for a network whose synchronous region is unbounded, λ2 determines the synchronizability. The larger the λ2 is, the better synchronizability the network has. For a network with a bounded synchronous region, the synchronizability is decided by R = λmax /λ2 . If R is small enough (1 is the best), the network has strong synchronizability. In order to meet the needs of the following theoretical derivation, a lemma is given here: Lemma 1 ([33]). let A, B be two square matrices, M be an integer, then

⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐

A B

.. .

B

B A

.. .

··· ··· .. .

B

···

B ⏐ ⏐ B ⏐

⏐

.. ⏐⏐ . ⏐ ⏐ A

= |A + (M − 1) B| · |A − B|(M −1)

(5)

M ×M

2.2. Multiplex chain network structure Generally, a chain network is a rule network, and many complex networks can be approximated as being formed by interconnecting and interacting with nodes of a plurality of chain networks. Multiplex chain networks M studied in this paper, with M layers each consisting of N nodes are a pair M = (G , C ) where G = {GK ; K ∈ {1, . . . , M }} is a family of undirected graphs GK = (XK , EK ) with XK = {x1 , . . . , xN }, in which X1 = X2 = · · · = XM = X . C = {EKL ⊆ XK × XL ; K , L ∈ {1, . . . , M }, K ̸ = L} is the set of connections between nodes of layers GK and GL with K ̸ = L. The elements of C are called crossed layers, and the elements of each EK are called intralayer connections of M in contrast with the elements of each EKL (K ̸ = L) that are called interlayer connections, and the only possible type of interlayer connections are those in which a given node is only connected to its counterpart nodes in the rest of layers, i.e., EKL = {(x, x); x ∈ X } for every 1 ≤ K ̸ = L ≤ M [22]. There are two kinds of connection modes between nodes of layers : one-to-one undirected coupling, as shown in Fig. 1(a), there are undirected connections between nodes across each two layers; one-to-one unidirectional coupling, as shown in Fig. 1 (b), there are unidirectional connections between nodes across each two layers. In other words, the given nodes are one-to-one mapped into the corresponding nodes in other layers, the one to one mapping between the nodes is performed by interlayer connection probability p (p = 1). For the multiplex chain network, the nodes xi (i = 1, 2, . . . , N) on each layer are the same, which has N ∗ M nodes in all, a is the intralayer coupling strength, d is the interlayer coupling strength. In order to analyze the synchronizability of two kinds of multiplex chain networks, we will calculate the eigenvalues of the networks in the next two sections and then analyze the synchronizability of the networks by using computational simulation.

4

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

3. The eigenvalues of the multiplex chain networks with one-to-one undirected coupling between layers For the sake of simplicity, we put a class of the multiplex chain networks with one-to-one undirected coupling between layers record as Networks-A. First, we explore the eigenvalue spectrum of Networks-A consisting of M layers each consisting of N nodes. The dynamics of the network nodes are as shown in (1). The supra-Laplacian matrix of Networks-A is: ⎡

a + (M − 1) d

⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ L=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

−a

···

0

−d

0

···

0

··· ··· ··· ···

−d

0

···

0

⎤

2a + (M − 1) d · · ·

0

0

−d

···

0

··· ··· ··· ···

0

−d

···

0

. . .

. . .

. . .

.

. . .

. . .

.

··· ··· ··· ···

. . .

. . .

0

0

···

−d

··· ··· ··· ···

0

0

···

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

−d

0

···

0

0

−d

···

0

. . .

. . .

.

. . .

0

···

−a

···

0

2a + (M − 1) d · · ·

0

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

−a

. . .

. . .

.

0

0

· · · a + (M − 1) d

.

.

.

.

.

−d

0

···

0

a + (M − 1 ) d

0

−d

···

0

−a

. . .

. . .

.

. . .

. . .

. . .

.

0

0

···

−d

0

0

· · · a + (M − 1) d

.

.

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

.

.

−a

···

0

2a + (M − 1) d · · ·

0

.

.

.

.

.

.

. .

.

.

. .

.

.

. .

. .

. .

.

. . .

.

. .

. .

. .

.

.

. .

.

.

.

0

···

0

0

−d

···

0

···

···

···

···

0

. . .

. . .

.

. . .

···

···

···

. . .

. . .

.

···

0

0

···

0

0

· · · −d

.

.

0

−d 0 · · · 0 a + (M − 1) d

−d

−d

−d · · · 0 .

.

−a

. . .

.

.

.

.

−d

−d

. . .

. . .

.

0

0

· · · a + (M − 1) d

.

.

. . .

According to Lemma 1, the characteristic polynomial of L can be transformed into:

⏐ ⏐ λ−a ⏐ ⏐ a ⏐ ⏐ ⏐⏐ 0 ⏐ ⏐ ⏐ ⏐λI − L˜ ⏐ = ⏐⏐ . ⏐ .. ⏐ ⏐ 0 ⏐ ⏐

a

0

λ − 2a

a

a

.. .

λ − 2a .. .

0

0

0

0

0

⏐ ⏐ λ − a − Md ⏐ ⏐ a ⏐ ⏐ ⏐ 0 · ⏐⏐ .. ⏐ . ⏐ ⏐ 0 ⏐ ⏐

··· ··· ···

0

0

0

0

0

0

a

λ−a

.. . λ − 2a

··· ··· a

.. .

a

0

λ − 2a − Md

a

a

.. .

λ − 2a − Md .. .

0

0

0

0

0

⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐

··· ··· ···

0

0

0

0

0

0

.. .

··· ···

.. .

λ − 2a − Md

a

a

λ − a − Md

⏐M −1 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐

the eigenvalues of L are: 0, 4asin2 (

where Md +

kπ 2N

), Md

,

Md + 4asin2



kπ 4asin2 ( 2N )



M −1

( kπ )

, k = 1, 2, . . . , N − 1.

2N



and Md are M − 1 multiple roots. λmax is the maximum eigenvalue and λ2 is the minimum ( )

non-zero eigenvalue of L, then λmax = Md + 4asin2

( N −1 ) { } π π , λ2 = min Md, 4asin2 ( 2N ) ,R = 2N

λmax λ2

λmax λ2

=

−1 π Md+4asin2 N2N

{ π }. min Md,4asin2 2N

( )

The relationships between λ2 , R = and the structural parameters of the networks are summarized in Table 1. According to the master stability{function method, } the synchronizability of Networks-A with unbounded synchronous π regions is determined by λ2 = min Md, 4asin(2 ( 2N )) ; if Networks-A have bounded synchronous regions, the synchronizability is determined by R =

λmax λ2

parameters are shown in Table 2.

=

−1 π Md+4asin2 N2N

{ π }. min Md,4asin2 2N

( )

The relationships between the synchronizability and the structural

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

5

Table 1 Changes of λ2 , R with N, M, a, d of Networks-A.

π λ2 = min Md, 4asin2 ( 2N )

{

}

(

R=

d> d<

)

−1 π Md+4asin2 N2N { π } min Md,4asin2 2N

(

d<

)

d>

(

π 4asin2 2N M 4asin2 π

)

λ2 = Md

( 2N )

λ2 = 4asin ( 2N )

M 4asin2

(

( 2N ) π

R=

M 4asin2

π

2

−1 π Md+4asin2 N2N

R=

M

−1 π Md+4asin2 N2N

Increase of a

Increase of d

↑

−

↑

↓

−

↑

−

↑

↓

↑

↓

↑

↑

↓

↑

)

)

(

π 4asin2 2N

Increase of M

−

)

Md

(

π ( 2N )

Increase of N

− : unchange; ↑ : increase; ↓ : decrease. Table 2 Changes of synchronizability with N, M, a, d of Networks-A. The case of synchronous region unbounded d<

(

π 4asin2 2N M

)

d>

(

π 4asin2 2N M

)

The case of synchronous region bounded d<

(

π 4asin2 2N M

)

d>

(

π 4asin2 2N M

)

Synchronizability With With With With

the the the the

increase increase increase increase

of of of of

N M a d

− ↑ − ↑

↓ − ↑ −

↓ ↑ ↓ ↑

↓ ↓ ↑ ↓

− : unchange; ↑ : strengthen; ↓ : weaken.

4. The eigenvalues of the multiplex chain networks with one-to-one unidirectional coupling between layers A class of the multiplex chain networks with one-to-one unidirectional coupling between layers are called Networks-B for short. The eigenvalue spectrum of Networks-B consisting of M layers each consisting of N nodes is explored in this section. The supra-Laplacian matrix of Networks-B is: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ L=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

a + (M − 1) d

−a . . . 0

−a ··· 0 2a + (M − 1) d · · · 0 . . . . .. . . . 0 · · · a + (M − 1) d

−d

0

0

. . .

−d . . .

··· ··· .. .

0

0

···

−a

···

0

0

···

0

a + (M − 2) d

0

0

. . .

··· .. .

0

. . .

−a . . .

0

0

···

0

0

. . . . . . . . . . . .

..

. . . . . . . . . . . . 0

0

0

0

. . .

. . .

··· ··· .. .

0

0

···

. . .

2a + (M − 2) d · · ·

. . .

..

0

···

.

. ..

. ..

.

0

· · · · · · · · · · · · −d 0 · · · 0 · · · · · · · · · · · · 0 −d · · · . . . . . . . .. . ··· ··· ··· ··· . . −d ··· ··· ··· ··· 0 0 ··· . . .. . . . 0 . . . . .. . . . 0 . . . . . . . .. . . . . . . . .. . . . a + (M − 2) d . . .. . −d 0 · · · .. . 0 −d · · · . . . .. . . .. . . . ..

..

.

0

0

0

···

···

···

···

0

. . .

···

···

···

···

. . .

0

0

0

.

0

0

· · · 0 a −a 0 · · · 0 −a 2a . . . . . . .. . . . . . . . 0 ··· 0 0 0

··· ··· ··· .. . ···

0

⎤

0 ⎥

⎥ ⎥ . ⎥ . ⎥ . ⎥ ⎥ −d⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ . ⎥ . ⎥ . ⎥ ⎥ ⎥ ⎥ −d⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ . ⎥ . ⎥ . ⎦ a

According to the Laplace Expansion Theorem, the characteristic polynomial of L can be transformed into: ⏐ ⏐ a 0 ··· 0 0 ⏐λ − a − (M − 1) d ⏐ ⏐ ⏐ ⏐ ⏐ a λ − 2a − (M − 1) d a ··· 0 0 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ 0 a λ − 2a − (M − 1) d · · · 0 0 ⏐ ⏐ ⏐ ⏐ ⏐λI − L˜ ⏐ = ⏐ ⏐· .. .. .. .. .. .. ⏐ ⏐ . . . . . . ⏐ ⏐ ⏐ ⏐ 0 0 0 · · · λ − 2a − (M − 1) d a ⏐ ⏐ ⏐ ⏐ 0 0 0 ··· a λ − a − (M − 1) d

6

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

⏐ ⏐λ − a − (M − 2) d a 0 ⏐ a λ − 2a − (M − 2) d a ⏐ ⏐ 0 a λ − 2a − (M − 2) d ⏐ · ⏐⏐ .. .. .. . . . ⏐ ⏐ 0 0 0 ⏐ ⏐ 0 0 0 ⏐ ⏐ ⏐ λ−a ⏐ a 0 ··· 0 0 ⏐ ⏐ λ − 2a a ··· 0 0 ⏐ a ⏐ ⏐ 0 ⏐ a λ − 2a · · · 0 0 ⏐ ⏐ ⏐ · · · ⏐⏐ .. .. .. . . .. ⏐ . . . . . . . ⏐ . ⏐ ⏐ 0 ⏐ 0 0 · · · λ − 2a a ⏐ ⏐ ⏐ 0 0 0 ··· a λ−a ⏐

··· ··· ··· .. . ··· ···

0 0 0

.. . λ − 2a − (M − 2) d a

⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ .. ⏐ . ⏐ ⏐ a ⏐ λ − a − (M − 2) d⏐ 0 0 0

the eigenvalues of L are: kπ kπ kπ kπ kπ ), d + 4asin2 ( 2N ), 2d + 4asin2 ( 2N ), 3d + 4asin2 ( 2N ), . . . , (M − 1)d + 4asin2 ( 2N ), 0, d, 2d, 3d, . . . , (M − 1) d, 4asin2 ( 2N k = 1, 2, . . . , N − 1. ( )

λmax = (M − 1) d + 4asin2

( N −1 ) { } π π , λ2 = min d, 4asin2 ( 2N ) , R = 2N

λmax λ2

λmax λ2

=

−1 π (M −1)d+4asin2 N2N } { π 2 min d,4asin ( 2N )

. The relationships

between λ2 , R = and the structural parameters are summarized in Table 3. According to the master stability function {method, the }synchronizability of Networks-B with unbounded synπ chronous regions is determined by λ2 = min d, 4asin2 ( 2N ) ; the synchronizability is determined by R = λλmax = ) ( −1 π (M −1)d+4asin2 N2N { π } , min d,4asin2 ( 2N )

2

with bounded synchronous regions. The relationships between the synchronizability and the struc-

tural parameters are shown in Table 4. 5. Numerical simulation Numerical examples are provided to illustrate the theoretical results. It can be seen from Tables 1–4 that the synchronizability of the networks is closely related to the size of the interlayer coupling strength. Our numerical simulation experiments are performed in three steps. Firstly, when the interlayer coupling strength is very small (take d = 0.0001), the synchronizability varies with respect to the number of nodes in each layer, intralayer coupling strength, and the number of layers. Secondly, when the interlayer coupling strength is slightly larger (take d = 1), the synchronizability varies with respect to the number of nodes in each layer, intralayer coupling strength, and the number of layers. Finally, the number of nodes in each layer, the intralayer coupling strength, and the number of layers are fixed, to analyze the synchronizability varies with respect to the interlayer coupling strength. 5.1. The synchronizability of Networks-A To begin with, we explored when the interlayer coupling strength is very small (take d = 0.0001), the synchronizability varies with respect to the number of nodes in each layer, intralayer coupling strength, and the number of layers. (1) Fig. 2 displays λ2 and R for Networks-A varying with respect to the number of nodes in each layer N. It can be seen from the left panel of Fig. 2 that λ2 remains invariant at a specific value λ2 = Md = 1.5 × 10−3 with increasing N (when π π N < 1 ), and then λ2 decreases with increasing N (when N > 1 ), with unbounded synchronous regions. 2 arcsin ( Md )2 4a

2 arcsin ( Md )2 4a

The observation reveals that synchronizability of Networks-A remains invariant at first, then gradually gets weakened with increasing the number of nodes in each layer. The right panel reveals that R increases gently at the beginning with increasing N, then increases sharply, and the synchronizability is weakened with increasing the number of nodes in each layer, with bounded synchronous regions. Fig. 2 reveals that increasing the number of nodes in each layer (increasing the network modularity) hinders the synchronizability of the networks. Both panels reveal that the synchronizability of Networks-A is weakened with increasing nodes size N, especially for networks belonging to the class with bounded synchronous regions. Md (2) From Fig. 3(a), it is clear that λ2 increases linearly from zero (when a < 2 π ) and then almost levels off at an upper bound value λ2 = Md = 0.0015 (when a >

Md π ) ). 4sin2 ( 2N

4sin ( 2N )

This means that the synchronizability is strengthened at

first and then remains invariant with increasing a (unbounded synchronous region). From panel (b), one can observe Md that for Networks-A, with increasing a, R first decreases slowly (when a < 2 π ), and then increases sharply (when a >

Md π ) ). 4sin2 ( 2N

4sin ( 2N )

The right panel reveals that the synchronizability of Networks-A is strengthened slowly at first

and then weakened sharply with increasing intralayer coupling strength (bounded synchronous region). That is to say, when the interlayer coupling strength is fixed and the intralayer coupling strength is increased to a certain extent, the

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

7

Fig. 2. The synchronizability of Networks-A v s. varying the number of nodes N (d = 0.0001, a = 1, M = 15 ). (a) λ2 v s. varying N (N varies from 10 to 600) (subgraph: λ2 v s. varying N (N varies from 10 to 80)); (b) R v s. varying N (N varies from 10 to 600) (subgraph: R v s. varying N (R varies from 10 to 80)).

Fig. 3. The synchronizability of Networks-A v s. varying intralayer coupling strength a (N = 100, d = 0.0001, M = 15). (a) λ2 v s. varying a (a varies from 0.01 to 2); (b) R v s. varying a (a varies from 0.01 to 2) (subgraph: R v s. varying a (a varies from 0.01 to 1.5)). Table 3 Changes of λ2 , R with N, M, a, d of Networks-B.

( ) } d < 4asin2 π π ( 2N ) λ2 = d λ2 = min d, 4asin2 ( 2N ) π π d > 4asin2 2N λ2 = 4asin2 ( 2N ) {

R=

( ) −1 π (M −1)d+4asin2 N2N { π } min d,4asin2 ( 2N )

d < 4asin2 d > 4asin

(

( 2

π

)

2N

π 2N

)

R= R=

(

)

(

)

−1 π (M −1)d+4asin2 N2N

d −1 π (M −1)d+4asin2 N2N

(

π 4asin2 2N

)

Increase of N

Increase of M

Increase of a

Increase of d

− ↓

− −

− ↑

↑ −

↑

↑

↑

↓

↑

↑

↓

↑

− : unchange; ↑ : increase; ↓ : decrease.

synchronizability of the networks will be weakened sharply. This phenomenon occurs because the intralayer coupling strength far exceeds the interlayer coupling strength, the structure of the multiplex network is similar to the community structure of multiple clusters, the connections within the cluster blocks are tight while the connections between the blocks Md are sparse and therefore the synchronizability is weakened, but as long as a ∈ (0, ], the synchronizability will 2 π 4sin

( 2N )

continue to strengthen with the intralayer coupling strength increasing. The synchronizability of Networks-A is optimal at

8

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

Fig. 4. The synchronizability of Networks-A v s. varying the number of layers M (N = 100, d = 0.0001, a = 1). (a) λ2 v s. varying M (M varies from 2 to 30); (b) R v s. varying M (M varies from 2 to 30) (subgraph: R v s. varying M (M varies from 10 to 30)). Table 4 Changes of synchronizability with N, M, a, d of Networks-B. The case of synchronous region unbounded d < 4asin

2

(

π

)

2N

d > 4asin

2

(

π 2N

)

The case of synchronous region bounded d < 4asin2

(

π 2N

)

d > 4asin2

(

π

)

2N

Synchronizability With With With With

the the the the

increase increase increase increase

of of of of

N M a d

− − − ↑

↓ − ↑ −

↓ ↓ ↓ ↑

↓ ↓ ↑ ↓

− : unchange; ↑ : strengthen; ↓ : weaken.

a∗A =

Md π ), 4sin2 ( 2N

and the larger the network scale and the greater the interlayer coupling strength, the greater the intralayer

coupling strength is needed to achieve optimal synchronization. (3) In Fig. 4, the left panel (a) plots λ2 varying with respect to the number of network layers, with increasing M, λ2 4asin2

(π)

4asin2

(π)

2N 2N first increases linearly (when M < ), and then remains invariant (when M > ). The synchronizability d d is strengthened at first and then remains invariant with increasing M, with unbounded synchronous region. R decreases

4asin2

(π)

4asin2

(π)

2N 2N sharply at first (when M < ), and then increases slowly with increasing M (when M > ), as shown d d in Fig. 4(b). That is to say, the synchronizability of the networks is strengthened at first and then weakened with increasing of layers, with bounded synchronous region. The synchronizability of Networks-A is optimal at ⌈ the number ⌉

MA∗ =

( )

π 4asin2 2N d

.

Secondly, we explored when the interlayer coupling strength of the networks is slightly larger (d = 1), the synchronizability varies with respect to the number of nodes in each layer, the intralayer coupling strength, and the number of network layers. (1) Fig. 5 shows that increasing the number of nodes in each layer N decreases λ2 (decreases gradually and approaches 0) and therefore the synchronizability of Networks-A. Increasing N increases R and the synchronizability will be weakened. That is to say, the synchronizability is weakened with increasing the number of nodes in each layer for the slightly larger interlayer coupling strength. (2) As is shown in Fig. 6, it is observed from panel (a) that λ2 increases linearly with increasing intralayer coupling strength a, and finally arrives at a certain value. This means that, if the considered Networks-A have unbounded synchronous regions, the synchronizability is first enhanced and then remains invariant with increasing intralayer coupling strength. Panel (b) displays that the ratio R decreases first and then increases sharply with increasing a. It implies that the synchronizability of Networks-A with bounded synchronous regions is enhanced firstly, then reaches Md ≈ 1.5 × 104 . the maximum, and finally gets weakened. The synchronizability will also be optimal at a∗A = 2 π 4sin ( 2N )

The greater the interlayer coupling strength, the greater the intralayer coupling strength be needed to achieve optimal synchronization. (3) Fig. 7 shows that, for d = 1, λ2 remains invariant with increasing the number of layers M, while it has little effect on synchronizability of Networks-A belong to the class of unbounded synchronous regions. From the right panel, it is

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

9

Fig. 5. The synchronizability of Networks-A v s. varying the number of nodes N (d = 1, a = 1, M = 15 ). (a) λ2 v s. varying N (N varies from 10 to 600) (subgraph: λ2 v s. varying N (N varies from 100 to 600)); (b) R v s. varying N (N varies from 10 to 600).

Fig. 6. The synchronizability of Networks-A v s. varying intralayer coupling strength a (N = 100, d = 1, M = 15). (a) λ2 v s. varying a (a varies from 0.01 to 2) (subgraph: λ2 v s. varying a (a varies from 20 to 18 000)); (b) R v s. varying a (a varies from 0.01 to 2) (subgraph: R v s. varying a (a varies from 20 to 18 000)).

Fig. 7. The synchronizability of Networks-A v s. varying the number of layers M (N = 100, d = 1, a = 1). (a) λ2 v s. varying M (M varies from 2 to 30); (b) R v s. varying M (M varies from 2 to 30).

10

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

Fig. 8. The synchronizability of Networks-A v s. varying interlayer coupling strength d (N = 100, a = 30, M = 15). (a) λ2 v s. varying d (d varies from 0.0001 to 0.02) (subgraph: λ2 v s. varying d (d varies from 0.002 to 5)); (b) R v s. varying d (d varies from 0.0001 to 0.02) (subgraph: R v s. varying d (d varies from 0.002 to 5)).

obtained that the eigenratio R grows with increasing M. That is, for Networks-A with unbounded synchronous regions, the number of layers M has no impact on the synchronizability, but results in worse synchronizability with bounded synchronous regions. Finally, we explored the synchronizability of the networks varies with respect to the interlayer coupling strength when the intralayer coupling strength, the number of network layers and the number of nodes in each layer are fixed. Fig. 8 shows, the synchronizability of Networks-A varying with respect to the interlayer coupling strength d. It is clear that λ2 increases quickly for relatively small values of d (when d < π

bound value λ2 = 4asin ( 2N ) (when d > 2

4asin2

( 2N )

M

π

( )

π 4asin2 2N M

) and finally almost levels off at an upper

). The left panel reveals that the synchronizability of Networks-A is

strengthened with increasing interlayer coupling strength when that strength is low and is invariant to the interlayer coupling strength when that strength is slight high. It is apparent that R decreases sharply with increasing d (when d <

( )

π 4asin2 2N M

). At slightly larger interlayer coupling strength, R arrives at its minimal value and then increases slowly 4asin2

(π)

2N with increasing d (when d > ), with bounded synchronous region. The right panel reveals that synchronizability M of the networks is strengthened at first and then weakened with increasing interlayer coupling strength. Therefore, there

is an optimum value d∗A =

( )

π 4asin2 2N M

for maximizing the synchronizability.

5.2. The synchronizability of Networks-B In the following, using a similar approach to that of Section 5.1, we explored the synchronizability of Networks-B. From Figs. 9–13, the influences of the number of nodes in each layer, the interlayer coupling strength, and the intralayer coupling strength on synchronizability for Networks-B are similar to Networks-A. Two types of the networks have different structural parameter values to achieve the optimal synchronizability. There are an intralayer coupling strength ) ( πoptimal to optimize the synchronizability value a∗B = 2d π and an optimal interlayer coupling strength value d∗B = 4asin2 2N 4sin ( 2N )

for Networks-B. (π) Figs. 14–15 show that λ2 = d = 1 × 10−4 (when d < 4asin2 2N ), the synchronizability ( π ) of Networks-B is only π determined by interlayer coupling strength d; λ2 = 4asin2 ( 2N ) = 0.0099 (when d > 4asin2 2N ), the synchronizability is determined by intralayer coupling strength a and the number of nodes (in each layer N, with unbounded synchronous ( ) ) π π region. Regardless of the area where d > 4asin2 2N or d < 4asin2 2N , λ2 is not affected by the change of the number of layers,( the )synchronizability is independent of the number of) layers. If the synchronous region is bounded, ( R =

(M −1)d+4asin2 d

N −1 π 2N

(when d < 4asin2

(

π 2N

)

); R =

(M −1)d+4asin2 4asin2

N −1 π 2N

( 2N ) π

(when d > 4asin2

(

π 2N

)

). That is, R increases

with increasing M, which indicates that the synchronizability of Networks-B is weakened with increasing the number of layers. All the above experimental results are completely consistent with the theoretical conclusions of Sections 3 and 4, thus verifying the correctness of the theoretical derivation.

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

11

Fig. 9. The synchronizability of Networks-B v s. varying the number of nodes N (d = 0.0001, a = 1, M = 15 ). (a) λ2 v s. varying N (N varies from 10 to 600) (b) R v s. varying N (N varies from 10 to 600) (subgraph: R v s. varying N (N varies from 10 to 300)).

Fig. 10. The synchronizability of Networks-B v s. varying the number of nodes N (d = 1, a = 1, M = 15 ). (a) λ2 v s. varying N (N varies from 10 to 600) (b) R v s. varying N (N varies from 10 to 600).

Fig. 11. The synchronizability of Networks-B v s. varying intralayer coupling strength a (N = 100, d = 0.0001, M = 15). (a) λ2 v s. varying a (a varies from 0.001 to 0.2); (b) R v s. varying a (a varies from 0.001 to 0.2).

12

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

Fig. 12. The synchronizability of Networks-B v s. varying intralayer coupling strength a (N = 100, d = 1, M = 15). (a) λ2 v s. varying a (a varies from 0.001 to 0.2) (subgraph: λ2 v s. varying a (a varies from 20 to 2000)); (b) R v s. varying a (a varies from 0.001 to 0.2)(subgraph: R v s. varying a (a varies from 20 to 2000)).

Fig. 13. The synchronizability of Networks-B v s. varying interlayer coupling strength d (N = 100, a = 30, M = 15). (a) λ2 v s. varying d (d varies from 0.0001 to 0.04) (subgraph: λ2 v s. varying d (d varies from 0.03 to 9)); (b) R v s. varying d (d varies from 0.0001 to 0.04) (subgraph: R v s. varying d (d varies from 0.03 to 9)).

Fig. 14. The synchronizability of Networks-B v s. varying the number of layers M (N = 100, d = 0.0001, a = 10, d < 4asin2 M (M varies from 2 to 30); (b) R v s. varying M (M varies from 2 to 30).

(

π 2N

)

). (a) λ2 v s. varying

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

13

Fig. 15. The synchronizability of Networks-B v s. varying the number of layers M (N = 100, d = 1, a = 10, d > 4asin2 (M varies from 2 to 30); (b) R v s. varying M (M varies from 2 to 30).

(

π

). (a) λ2 v s. varying M

)

2N

6. Discussion and conclusion In this paper, two types of multiplex chain networks are studied, one is a class of the multiplex chain networks with one-to-one undirected coupling between layers (Networks-A), and the other is a class of the multiplex chain networks with one-to-one unidirectional coupling between layers (Networks-B). We studied analytically, numerically and experimentally how the structural parameters influence the ability of the whole system to synchronize. We have strictly derived the eigenvalue spectrum of the supra-Laplacian matrices of two kinds of multiplex chain networks, and analyzed the relationships between the synchronizability and structural parameters. According to the theoretical derivation and numerical simulation results, we found that, for two types of networks studied in this paper, the changes in the number of nodes in each layer, the interlayer coupling strength and the intralayer coupling strength have a similar effect on their synchronizability, but the two have different structural parameter values to achieve the optimal synchronizability, the optimal values for interlayer coupling strength and the intralayer coupling strength to achieving the optimal synchronizability are different (for Networks-A: a∗A =

Md π ), 4sin2 ( 2N

d∗A =

( )

π 4asin2 2N M

; for Networks-B: a∗B =

d π ), 4sin2 ( 2N

π ). Except that, for the sufficiently small value of interlayer coupling strength d, the synchronizability d∗B = 4asin 2N of Networks-A with unbounded synchronous regions is positively correlated with the number of network layers (d ∈

( 2

4asin2

)

(π)

π 2N (0, ]), that of Networks-B is independent of the number of network layers (d ∈ (0, 4asin2 2N ]). For bounded M synchronous regions, the synchronizability of Networks-A is also positively correlated with the number of network layers 4asin2

(

(π)

)

π 2N (d ∈ (0, ]), that of Networks-B is negatively correlated with the number of network layers (d ∈ (0, 4asin2 2N ]). M For unbounded synchronous regions, the synchronizability of two types of the multiplex chain networks is strengthened with increasing the interlayer coupling strength d, when d increases to a certain extent (d∗ ), then d increases again and no longer affects the synchronizability. For bounded synchronous regions, the synchronizability is enhanced with increasing interlayer coupling strength d when that strength is low and is weakened to the interlayer coupling strength when that strength is slight high. Therefore, there is an optimal value d∗ for the interlayer coupling strength. The intralayer coupling strength is positively correlated with the interlayer coupling strength (a ∝ d), increase a while increasing d, only in this way can the multiplex chain network achieve optimal synchronization. This is because when a far exceeds d, increasing d is beneficial to reducing the clustering property of network, more conducive to synchronization. Therefore, enhancing the synchronizability does not mean increasing the interlayer coupling strength infinitely, but increasing the interlayer coupling strength to a certain extent and then increasing the intralayer coupling strength, so as the two achieve a certain degree of equilibrium. These results provide guidances for building an optimized synchronization network. Generally, the multiplex networks are large in scale and complicated in structure, and it is difficult to obtain an exact solution for their eigenvalues. In this paper, we only studied the multiplex chain network with the same number of network nodes in each layer, the same intralayer coupling strength, and one-to-one fully coupled between layers. Due to the simplification of the studied network structure, there is a certain gap with the complex structure of the actual network. For example, when the interlayer connection probability of the networks is less than 1 (incomplete connection), when the nodes are connected to the central node at each layer, how are the eigenvalues of these networks represented ? How does the synchronizability of the networks change ? For multiplex chain networks, how to change other parameters to keep the synchronizability unchanged while the intralayer coupling strength or interlayer coupling strength changes the synchronizability of the networks, these are waiting for us to continue to explore. A more systematic analysis approach

(

)

14

Y. Deng, Z. Jia, G. Deng et al. / Physica A 537 (2020) 122631

is a part of our future work as well. The results and understanding based on these limiting structures can be illustrative and instructive for building insight towards synchronizability of more realistic multiplex networks. Acknowledgments This project is supported by National Natural Science Foundation of China (Nos. 61563013, 61663006) and the Natural Science Foundation of Guangxi, China (No. 2018GXNSFAA138095). References [1] L. Tang, J. Lu, X. Wu, J. Lü, Impact of node dynamics parameters on topology identification of complex dynamical networks, Nonlinear Dynam. 73 (2013) 1081–1097. [2] X. Wang, X. Fan, Complex networks: topology, dynamics and synchronization, Int. J. Bifurcation Chaos 12 (2002) 885–916. [3] P. Wang, Y. Chen, J. Lü, Q. Wang, X. Yu, Graphical features of functional genes in human protein interaction network, IEEE Trans. Biomed. Circuits Syst. 10 (2016) 707–720. [4] W. Wang, B. Wang, B. Hu, et al., General dynamics of topology and traffic on weighted technological networks, Phys. Rev. Lett. 94 (2005) 188702. [5] K. Havva, T. Mahbubeh, D.P. Luisa, et al., The basic principles of topology-dynamics relations in networks: An empirical approach, Physica A 508 (2018) 584–594. [6] H. Liu, J. Lu, J. Lü, et al., Structure identification of uncertain general complex dynamical networks with time delay, Automatica 45 (2009) 1799–1807. [7] P. Wang, X. Yu, J. Lu, Identification and evolution of structurally dominant nodes in protein-protein interaction networks, IEEE Trans. Biomed. Circuits Syst. 8 (2014) 87–97. [8] X. Wu, W. Wei, L. Tang, et al., Coreness and h-index for weighted networks, IEEE Trans. Circuits Syst. I. Regul. Pap. 66 (2019) 3113–3122. [9] G. Mei, X. Wu, Y. Wang, et al., Compressive-sensing-based structure identification for multilayer networks, IEEE Trans. Cybern. 48 (2018) 754–764. [10] E. Boukobza, M. Chuchem, D. Cohen, et al., Phase-diffusion dynamics in weakly coupled bose–Einstein condensates, Phys. Rev. Lett. 102 (2009) 180403. [11] G. Mei, X. Wu, D. Ning, et al., Finite-time stabilization of complex dynamical networks via optimal control, Complexity 21 (2016) n/a-n/a. [12] G. Siudem, J.A. Hołyst, Diffusion on hierarchical systems of weakly-coupled networks, Physica A 513 (2019) 675–686. [13] X. Liu, T. Chen, Cluster synchronization in directed networks via intermittent pinning control, IEEE Trans. Neural Netw. 22 (2011) 1009–1020. [14] P. Wang, S. Xu, Spectral coarse grained controllability of complex networks, Physica A 478 (2017) 168–176. [15] P. Wang, C. Tian, J. Lu, Identifying influential spreaders in artificial complex networks, J. Syst. Sci. Complexity 27 (2014) 650–665. [16] M. Denker, M. Timme, M. Diesmann, et al., Breaking synchrony by heterogeneity in complex networks, Phys. Rev. Lett. 92 (2004) 074103. [17] A.E. Motter, C. Zhou, Jrgen Kurths, Network synchronization diffusion and the paradox of heterogeneity, Phys. Rev. E 71 (2005) 016116. [18] M. Xue, E. Yeung, A. Rai, et al., Initial-condition estimation in network synchronization processes: Algebraic and graphical characterizations of the estimator, Complex Syst. 21 (2013) 297–333. [19] Z. Jia, L. Zeng, Y. Wang, et al., Optimization algorithms for spectral coarse-graining of complex networks, Physica A 514 (2019) 925–935. [20] Y. Long, Z. Jia, Y. Wang, Coarse graining method based on generalized degree in complex network, Physica A 505 (2018) 655–665. [21] K. Lee, J. Kim, S. Lee, et al., Multiplex networks, Underst. Complex Syst. (2014) 53–72. [22] S. Boccaletti, G. Bianconi, R. Criado, et al., The structure and dynamics of multilayer networks, Phys. Rep. 544 (2014) 1–122. [23] L. Tang, X. Wu, J. Lü, J. Lu, et al., Master stability functions for complete intralayer, and interlayer synchronization in multiplex networks of coupled Rössler oscillators, Phys. Rev. E 99 (2019) 012304. [24] S. Gómez, A. Díaz-Guilera, J. Gómez-Gardeñes, et al., Diffusion dynamics on multiplex networks, Phys. Rev. Lett. 110 (2013) 028701. [25] X. Wu, W.X. Zheng, J. Zhou, Generalized outer synchronization between complex dynamical networks, Chaos 19 (2009) 013109. [26] S. Zhang, X. Wu, J. Lu, H. Feng, J. Lü, Recovering structures of complex dynamical networks based on generalized outer synchronization, IEEE Trans. Circuits Syst. I 61 (2014) 3216–3224. [27] Y. Li, X. Wu, J. Lu, et al., Synchronizability of duplex networks, IEEE Trans. Circuits Syst. II 63 (2017) 206–210. [28] L.M. Pecora, T.L. Carroll, Master stability functions for synchronized coupled systems, Int. J. Bifurcation Chaos 9 (1999) 2315–2320. [29] Y. Liang, X.Y. Wang, Pinning chaotic synchronization in complex networks on maximum eigenvalue of low order matrix, Acta Phys. Sin. 61 (2012) 038901, (in Chinese). [30] A. Solé-Ribalta, M. De Domenico, N.E. Kouvaris, et al., Spectral properties of the Laplacian of multiplex networks, Phys. Rev. E 88 (2013) 032807. [31] M. Xu, J. Lu, J. Zhou, Eigenvalue spectrum and synchronization ability of two-layer star networks, Acta Phys. Sin. 65 (2016) 028902, (in Chinese). [32] M. Xu, J.Zhou, J. Lu, et al., Synchronizability of two-layer networks, Eur. Phys. J. B 88 (2015) 240. [33] J. Sun, X. Li, J.H. Zhang, et al., Synchronizability and eigenvalues of multilayer star networks through unidirectionally coupling, Acta Phys. Sin. 66 (2017) 188901, (in Chinese). [34] J. Wei, X. Wu, J. Lu, et al., Synchronizability of duplex regular networks, Europhys. Lett. 120 (2017) 20005.