better than translation B. B: translation B is significantly better than
translation A. C: none of the translations are significantly better. NOTICE:
directly generate your decision in the format ‘Decision: A/B/C.‘
———————————————————————————————————————————————
# Generation: Decision: A.
Translation A is more accurate and faithful to the original sentence. It
correctly translates “jazz idioms” as “爵士乐的音调,” which conveys the intended
meaning of Milhaud’s use of jazz elements in his music. On the other hand,
translation B mistranslates “Suite for clarinet, violin and piano” as
“长号、小提琴和钢琴的四重奏,” which is not an accurate representation of the original work.
Additionally, the use of “米尔豪德” instead of “米哈伊尔德” may cause confusion for
readers who are familiar with the composer’s name. Overall, translation A is a
better quality translation. Figure 11: One case of comparison-based generation
on translation task. The preference label is generated by Qwen-72B-Chat, a
very powerful open-source LLM.

##  Appendix I Complete Results

###  I.1 Complete Results for Different Qualities

Table 9 and Table 10 show the overall evaluation results for different
qualities granularities on the subjective evaluation of the test and dev set
the feedback, comparison, and correction critique dimensions. Table 11 and
Table 12 show the overall evaluation results for different qualities
granularities on the objective evaluation of the test and dev set the
feedback, comparison, and correction critique dimensions.

###  I.2 Complete Results for Different Tasks

Table 13 and Table 14 show evaluation results on each tasks (test set) for the
feedback and comparison-based feedback dimensions. Table 15 and Table 16 show
complete evaluation results on each tasks (dev set) for the feedback and
comparison-based feedback dimensions. Table 17 and Table 18 show complete
evaluation results on each task (test and dev set) for the correction
dimension.


††thanks<EMAIL_ADDRESS>
# Group synchrony, parameter mismatches, and intragroup connections

Shirin Panahi  Francesco Sorrentino  University of New Mexico, Albuquerque New
Mexico 87106, USA

###### Abstract

Group synchronization arises when two or more synchronization patterns coexist
in a network formed of oscillators of different types, with the systems in
each group synchronizing on the same time-evolution, but systems in different
groups synchronizing on distinct time-evolutions. Group synchronization has
been observed and characterized when the systems in each group are identical
and the couplings between the systems satisfy specific conditions. By relaxing
these constraints and allowing them to be satisfied in an approximate rather
than exact way, we observe that stable group synchronization may still occur
in the presence of small deviations of the parameters of the individual
systems and of the couplings from their nominal values. We analyze this case
and provide necessary and sufficient conditions for stability through a master
stability function approach, which also allows us to quantify the
synchronization error. We also investigate the stability of group
synchronization in the presence of intra group connections and for this case,
extend some of the existing results in the literature. Our analysis points out
a broader class of matrices describing intra group connections for which the
stability problem can be reduced in a low-dimensional form.

##  I Introduction

A large literature has considered synchronization in networks of coupled
oscillators, see for example the review in Ref. [1]. In the case of identical
oscillators, stability of the synchronous solution for arbitrary network
topologies can be investigated through the master stability function approach,
introduced in Ref. [2]. More complex patterns of synchronous dynamics may
arise in networks formed of different types of oscillators. Group
synchronization in networks of coupled oscillators was first analyzed in [3],
followed by further analysis in Ref. [4] and an experimental realization of
group synchronization [5]. Moreover, group synchronization in networks of
fractional-order chaotic oscillators was investigated in Ref. [6].

According to the definition in Ref. [3], a _group_ is a set of oscillators
characterized by the same (uncoupled) dynamics, with oscillators of different
types belonging to different groups. Group synchronization is achieved when
the oscillators in each group synchronize on the same time evolution, with
these synchronized dynamics being different from group to group. Cluster
synchronization can be considered as a particular case of group
synchronization, but for which the oscillators in different clusters are of
the same type. Most recently, a canonical transformation for simultaneous
block diagonalization of matrices has been proposed to decouple the cluster
synchronization stability problem into subproblems of minimal dimensionality
[7]. The more general case of cluster synchronization in multilayer networks
formed of oscillators of different types has been studied in Ref. [8].

An experimental realization of group synchronization was performed in Ref.
[5], where stable group synchronization was observed under different
configurations for a small network of four coupled optoelectronic oscillators.
This observation suggests that group synchronization is robust with respect to
small parametric mismatches that are inevitable in experiments. The analysis
of stability of the group-synchronous solution in the presence of such small
parametric mismatches (affecting both the individual nodes being coupled and
the strengths of the couplings between them) is studied in this paper. The
effects of parameter mismatches in a network of _identical_ oscillators
achieving complete synchronization have been studied in Refs. [9, 10, 11] and
achieving cluster synchronization in Refs. [12, 13]. A second order expansion
to study complete synchronization in networks of identical nodes with
parameter mismatches has been studied in Refs. [14, 15]. However, to the best
of our knowledge, no paper has investigated the role of parameter mismatches
in networks formed of systems of different types, which is the main focus of
this paper. The other subject of this paper is the study of the conditions
under which group synchronization may arise in the presence of intragroup
connections, which was first considered in Ref. [4].

##  II Group synchronization in the presence of parametric mismatches

Though our results can be generalized to the case of an arbitrary number of
groups (for the case in which the stability of the group-synchronous solution
can be reduced in a low-dimensional form [4]), for simplicity in what follows
we will focus on the case of two groups. For the time being, we also maintain
the assumption introduced in [3] that the network topology is bipartite, i.e.,
connections exist only from group $X$ to group $Y$ and vice versa (this
assumption will be removed later in Sec.III).

We first present the ideal case that parameter mismatches are absent. For this
very special case, the dynamical equations are:

$\dot{\pmb{x}}_{i}=\pmb{F}\Bigl{(}\pmb{x}_{i},u_{i}^{x}\Bigr{)}=\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(\pmb{y}_{j})\Bigr{)},\quad
i=1,\cdots,N_{x}.$ (1a)
$\dot{\pmb{y}}_{j}=\pmb{G}\Bigl{(}\pmb{y}_{j},u_{j}^{y}\Bigr{)}=\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(\pmb{x}_{i})\Bigr{)},\quad
j=1,\cdots,N_{y}.$ (1b)

where $\pmb{x}_{i}$ ($\pmb{y}_{j}$) is an $n_{x}$-dimensional
($n_{y}$-dimensional) state vector of systems in group $\pmb{X}$ and in group
$\pmb{Y}$ respectively,
$\pmb{u}_{i}^{x}=\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(y_{j})$
$[\pmb{u}_{j}^{y}=\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(x_{i})]$ is
the input received by node $i$ in group $\pmb{X}$ ($j$ in group $\pmb{Y}$).
The functions
$\pmb{F}:\mathbf{R}^{n_{x}}\times\mathbf{R}^{n_{x}}\rightarrow\mathbf{R}^{n_{x}}$
and
$\pmb{G}:\mathbf{R}^{n_{y}}\times\mathbf{R}^{n_{y}}\rightarrow\mathbf{R}^{n_{y}}$
define the time evolution of the systems in group $\pmb{X}$ and group
$\pmb{Y}$ respectively, and the interaction functions
$\pmb{H}:\mathbf{R}^{n_{y}}\rightarrow\mathbf{R}^{n_{x}}$ and
$\pmb{L}:\mathbf{R}^{n_{x}}\rightarrow\mathbf{R}^{n_{y}}$ define the output of
systems in group $\pmb{Y}$ and group $\pmb{X}$, respectively. $N_{x}$ is the
number of systems in group $\pmb{X}$ and $N_{y}$ is the number of systems in
group $\pmb{Y}$. In what follows we assume without loss of generality that
$N_{x}\geq N_{y}$. The input $\pmb{u}_{i}^{x}$ received by node $i$ in group
$X$ ($\pmb{u}_{j}^{y}$ received by node $j$ in group $\pmb{Y}$) corresponds to
a superposition of the outputs from the nodes in group $\pmb{Y}$ (group
$\pmb{X}$) through the coefficients $A_{ij}^{NOM}$ ($B_{ji}^{NOM}$).
$A^{\mathrm{NOM}}$ is an $N_{x}\times N_{y}$ coupling matrix, whose entries
${A_{ij}^{\mathrm{NOM}}}$ represent the nominal strength of the direct
interaction form system $j$ in the $\pmb{Y}$-group to system $i$ in the
$\pmb{X}$-group. Analogously, $B_{ji}^{\mathrm{NOM}}$ is an $N_{y}\times
N_{x}$ coupling matrix, whose entries $B_{ji}^{\mathrm{NOM}}$ represent the
nominal strength of the direct interaction from system $i$ in the
$\pmb{X}$-group to $j$ in the $\pmb{Y}$-group. Note that Eqs. (1) are a
generalization of those considered in Refs. [3, 4].

In what follows, we consider the following two types of small parameter
mismatches from nominal conditions: (i) deviations from identicality of the
uncoupled dynamics of the individual systems in each group and (ii) deviations
from nominal conditions of the coupling strengths between systems.

If the following conditions hold,

$\displaystyle\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}$ $\displaystyle=a\neq
0,\quad i=1,\cdots,N_{x}.$ (2a)
$\displaystyle\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}$ $\displaystyle=b\neq
0,\quad j=1,\cdots,N_{y}.$ (2b)

then the synchronous solution exists,

$\displaystyle\dot{\pmb{x}}_{s}$
$\displaystyle=\pmb{F}\Bigl{(}\pmb{x}_{s},a\pmb{H}(\pmb{y}_{s})\Bigr{)},$ (3a)
$\displaystyle\dot{\pmb{y}}_{s}$
$\displaystyle=\pmb{G}\Bigl{(}\pmb{y}_{s},b\pmb{L}(\pmb{x}_{s})\Bigr{)}.$ (3b)

We note here that by replacing $a\pmb{H}\rightarrow\pmb{H}$ and
$A^{NOM}/a\rightarrow A^{NOM}$ ($b\pmb{L}\rightarrow\pmb{L}$ and
$B^{NOM}/b\rightarrow B^{NOM}$), it is always possible to set $a=b=1$. Hence,
without loss of generality, in what follows we set $a=b=1$. This has the
significant advantage that the synchronous solution, which obeys Eqs. (3) with
$a=b=1$, is independent of the particular choice of the bipartite network.

Stability of this synchronized solution can be reduced in a master stability
form [3]. Unfortunately, exact satisfaction of condition (2) is difficult to
implement in experiments and in real-world situations. The other requirement
whose exact satisfaction is hardly accomplished in experiments is that the
evolution of all the systems in the $\pmb{X}$-group (in the $\pmb{Y}$-group)
is exactly characterized by the same function $F$ ($G$). What can be
realistically achieved is that
$\pmb{F}(\pmb{x}_{i},\pmb{u}_{i}^{x})\rightarrow\pmb{F}(\pmb{x}_{i},\pmb{u}_{i}^{x},\mu_{i})$,
$i=1,...,N_{x}$ in Eq.(1a) and
$\pmb{G}(\pmb{y}_{j},\pmb{u}_{j}^{y})\rightarrow\pmb{G}(\pmb{y}_{j},\pmb{u}_{j}^{y},{\nu_{j}})$,
$j=1,...,N_{y}$ in Eq. (1b) where $\mu_{i}$ and $\nu_{j}$ are scalar
parameters that slightly vary from oscillator to oscillator 111Even though we
do not consider the case of two or more parameters that vary from system to
system, the extension of our results to this case is straightforward..
Analogously, it is realistic to assume that the couplings $A_{ij}$ and
$B_{ji}$ are affected by small mismatches with respect to their nominal
values,

$\displaystyle A_{ij}$ $\displaystyle=A_{ij}^{\mathrm{NOM}}+\delta A_{ij},$
(4a) $\displaystyle B_{ji}$ $\displaystyle=B_{ji}^{\mathrm{NOM}}+\delta
B_{ji},$ (4b)

where $\delta A_{ij}$ and $\delta B_{ji}$ are small deviations. Under these
assumptions, Eqs. (1) become

$\displaystyle\dot{\pmb{x}}_{i}$
$\displaystyle=\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}\pmb{H}(\pmb{y}_{j}),\mu_{i}\Bigr{)},\quad
i=1,\cdots,N_{x},$ (5a) $\displaystyle\dot{\pmb{y}}_{j}$
$\displaystyle=\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}\pmb{L}(\pmb{x}_{i}),\nu_{j}\Bigr{)},\quad
j=1,\cdots,N_{y},$ (5b)

which describe the dynamics of the two groups under realistic circumstances,
such as those that would be observed in any experimental setting.

We write $\mu_{i}=\bar{\mu}+\delta\mu_{i}$ where
$\bar{\mu}=(N_{x})^{-1}\sum_{i=1}^{N_{x}}\mu_{i}$ and $\delta\mu_{i}$ is a
small deviation. Similarly, we write $\nu_{j}=\bar{\nu}+\delta\nu_{j}$ where
$\bar{\nu}=(N_{y})^{-1}\sum_{j=1}^{N_{y}}\nu_{j}$ and $\delta\nu_{j}$ is a
small deviation. Note that by construction $\sum_{i=1}^{N_{x}}\delta\mu_{i}=0$
and $\sum_{j=1}^{N_{y}}\delta\nu_{j}=0$.

We can also write

$\displaystyle\sum_{j=1}^{N_{y}}A_{ij}$
$\displaystyle=a+\sum_{j=1}^{N_{y}}\delta A_{ij}=a+\delta\bar{a}+\delta
a_{i},$ (6a) $\displaystyle\sum_{i=1}^{N_{x}}B_{ji}$
$\displaystyle=b+\sum_{i=1}^{N_{x}}\delta B_{ji}=b+\delta\bar{b}+\delta
b_{j},$ (6b)

where

$\displaystyle\delta\bar{a}$
$\displaystyle=N_{x}^{-1}\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\delta A_{ij},$
(7a) $\displaystyle\delta\bar{b}$
$\displaystyle=N_{y}^{-1}\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\delta B_{ji},$
(7b)

are the average sums of the rows of the matrices $A$ and $B$ and

$\displaystyle\delta a_{i}$ $\displaystyle=\Bigg{(}\sum_{j=1}^{N_{y}}\delta
A_{ij}\Bigg{)}-\delta\bar{a}=\Bigg{(}\sum_{j=1}^{N_{y}}\delta
A_{ij}\Bigg{)}-N_{x}^{-1}\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\delta A_{ij},$
(8a) $\displaystyle\delta b_{j}$
$\displaystyle=\Bigg{(}\sum_{i=1}^{N_{x}}\delta
B_{ji}\Bigg{)}-\delta\bar{b}=\Bigg{(}\sum_{i=1}^{N_{x}}\delta
B_{ji}\Bigg{)}-N_{y}^{-1}\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\delta B_{ji},$
(8b)

are small deviations. The deviations $\delta a_{i}$ and $\delta b_{j}$ are
calculated with respect to the average row-sums $\delta\bar{a}$ and
$\delta\bar{b}$, hence they sum to zero, that is, $\sum_{i=1}^{N_{x}}\delta
a_{i}=0$ and $\sum_{j=1}^{N_{y}}\delta b_{j}=0$.

Unfortunately, for the case of equations (5), different from the case of
equations (1) an exact synchronous solution does not exist. It is possible,
however, that the individual trajectories stabilize in a nearly synchronous
state, where the trajectories of the nodes in the $\pmb{X}$-group (and in the
$\pmb{Y}$-group) remain close to each other, i.e., for which both

$\displaystyle\delta\pmb{x}_{i}(t)$
$\displaystyle=(\pmb{x}_{i}(t)-\bar{\pmb{x}}(t)),$ (9a)
$\displaystyle\delta\pmb{y}_{j}(t)$
$\displaystyle=(\pmb{y}_{j}(t)-\bar{\pmb{y}}(t)),$ (9b)

remain small in time, where $\bar{\pmb{x}}(t)$ and $\bar{\pmb{y}}(t)$ are the
average solutions,

$\displaystyle\bar{\pmb{x}}(t)$
$\displaystyle=(N_{x}^{-1})\sum_{i=1}^{N_{x}}\pmb{x}_{i}(t),$ (10a)
$\displaystyle\bar{\pmb{y}}(t)$
$\displaystyle=(N_{y}^{-1})\sum_{j=1}^{N_{y}}\pmb{y}_{j}(t),$ (10b)

obeying,

$\displaystyle\dot{\bar{\pmb{x}}}$
$\displaystyle=(N_{x})^{-1}\sum_{i=1}^{N_{x}}\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}\pmb{H}(y_{j}),\mu_{i}\Bigr{)},$
(11a) $\displaystyle\dot{\bar{\pmb{y}}}$
$\displaystyle=(N_{y})^{-1}\sum_{j=1}^{N_{y}}\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}\pmb{L}(\pmb{x}_{i}),\nu_{j}\Bigr{)}.$
(11b)

By differentiating (9) with respect to time, we obtain

$\displaystyle\delta\dot{\pmb{x}}_{i}$
$\displaystyle=\Bigl{[}\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}\pmb{H}(\pmb{y}_{j}),\mu_{i}\Bigr{)}-(N_{x})^{-1}\sum_{i=1}^{N_{x}}\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}\pmb{H}(\pmb{y}_{j}),\mu_{i}\Bigr{)}\Bigr{]},\quad
i=1,\cdots,N_{x},$ (12a) $\displaystyle\delta\dot{\pmb{y}}_{j}$
$\displaystyle=\Bigl{[}\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}\pmb{L}(\pmb{x}_{i}),\nu_{j}\Bigr{)}-(N_{y})^{-1}\sum_{j=1}^{N_{y}}\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}\pmb{L}(\pmb{x}_{i}),\nu_{j}\Bigr{)}\Bigr{]},\quad
j=1,\cdots,N_{y}.$ (12b)

By expanding to first order the function $\pmb{F}$ in (12a) about the point
$(\bar{\pmb{x}},\pmb{H}(\bar{\pmb{y}}),\bar{\mu})$ and the function $\pmb{G}$
in (12b) about the point $(\bar{\pmb{y}},\pmb{L}(\bar{\pmb{x}}),\bar{\nu})$,
we obtain,

$\displaystyle\pmb{F}\Bigl{(}\pmb{x}_{i},\sum_{j=1}^{N_{y}}A_{ij}\pmb{H}(\pmb{y}_{j}),\mu_{i}\Bigr{)}$
$\displaystyle\simeq\pmb{F}\Bigl{(}\bar{\pmb{x}},\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(\bar{\pmb{y}}),\bar{\mu}\Bigr{)}+\mathrm{D}\pmb{F}_{x}\Bigl{(}\bar{\pmb{x}},\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(\bar{\pmb{y}}),\bar{\mu}\Bigr{)}\delta{\pmb{x}}_{i}$
$\displaystyle+\mathrm{D}\pmb{F}_{u}\Bigl{(}\bar{\pmb{x}},\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(\bar{\pmb{y}}),\bar{\mu}\Bigr{)}\Bigl{[}\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\mathrm{D}\pmb{H}(\bar{\pmb{y}})\delta\pmb{y}_{j}+\sum_{j=1}^{N_{y}}\pmb{H}(\bar{\pmb{y}})\delta
A_{ij}\Bigr{]}$
$\displaystyle+\mathrm{D}F_{\mu}\Bigl{(}\bar{\pmb{x}},\sum_{j=1}^{N_{y}}A_{ij}^{\mathrm{NOM}}\pmb{H}(\bar{\pmb{y}}),\bar{\mu}\Bigr{)}\delta\mu_{i}$
$\displaystyle\pmb{G}\Bigl{(}\pmb{y}_{j},\sum_{i=1}^{N_{x}}B_{ji}\pmb{L}(\pmb{x}_{i}),\nu_{j}\Bigr{)}$
$\displaystyle\simeq\pmb{G}\Bigl{(}\bar{\pmb{y}},\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(\bar{\pmb{x}}),\bar{\nu}\Bigr{)}+\mathrm{D}\pmb{G}_{y}\Bigl{(}\bar{\pmb{y}},\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(\bar{\pmb{x}}),\bar{\nu}\Bigr{)}\delta{\pmb{y}}_{j}$
$\displaystyle+\mathrm{D}\pmb{G}_{u}\Bigl{(}\bar{\pmb{y}},\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(\bar{\pmb{x}}),\bar{\nu}\Bigr{)}\Bigl{[}\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\mathrm{D}\pmb{L}(\bar{\pmb{x}})\delta\pmb{x}_{i}+\sum_{i=1}^{N_{x}}\pmb{L}(\bar{\pmb{x}})\delta
B_{ji}\Bigr{]}$
$\displaystyle+\mathrm{D}\pmb{G}_{\nu}\Bigl{(}\bar{\pmb{y}},\sum_{i=1}^{N_{x}}B_{ji}^{\mathrm{NOM}}\pmb{L}(\bar{\pmb{x}}),\bar{\nu}\Bigr{)}\delta\nu_{i}$

Then Eqs. (12) can be rewritten,

$\displaystyle\delta\dot{\pmb{x}}_{i}(t)$
$\displaystyle=\mathrm{D}\pmb{F}_{x}^{*}\delta\pmb{x}_{i}(t)+\sum_{j=1}^{N_{y}}\tilde{A}_{ij}\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})\delta\pmb{y}_{j}+\pmb{U}^{x}_{i},\quad
i=1,...,N_{x},$ (13a) $\displaystyle\delta\dot{\pmb{y}}_{j}(t)$
$\displaystyle=\mathrm{D}\pmb{G}_{y}^{*}\delta\pmb{y}_{j}(t)+\sum_{i=1}^{N_{x}}\tilde{B}_{ji}\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})\delta\pmb{x}_{i}+\pmb{U}^{y}_{j},\quad
j=1,...,N_{y},$ (13b)

where $\pmb{U}^{x}_{i}=\mathrm{D}\pmb{F}_{u}^{*}\pmb{H}(\bar{\pmb{y}})\delta
a_{i}+\mathrm{D}\pmb{F}_{\mu}^{*}\delta\mu_{i}$ and
$\pmb{U}^{y}_{j}=\mathrm{D}\pmb{G}_{u}^{*}\pmb{L}(\bar{\pmb{x}})\delta
b_{j}+\mathrm{D}\pmb{G}_{\nu}^{*}\delta\nu_{j}$, and the superscript ∗ denotes
that the partial derivatives of the functions $\pmb{F}$ and $\pmb{G}$ are
evaluated about $(\bar{\pmb{x}},\pmb{H}(\bar{\pmb{y}}),\bar{\mu})$ and
$(\bar{\pmb{y}},\pmb{L}(\bar{\pmb{x}}),\bar{\nu})$, respectively. The matrices
$\tilde{A}=\\{\tilde{A}_{ij}\\}$ and $\tilde{B}=\\{\tilde{B}_{ij}\\}$ in (13)
are defined as follows,

$\displaystyle\tilde{A}_{ij}$ $\displaystyle=(A_{ij}^{NOM}-a_{j}),\qquad
i=1,\dots,N_{x},\quad j=1,\dots,N_{y},$ (14a) $\displaystyle\tilde{B}_{ji}$
$\displaystyle=(B_{ji}^{NOM}-b_{i}),\qquad i=1,\dots,N_{x},\quad
j=1,\dots,N_{y},$ (14b)

where $a_{j}=(N_{x})^{-1}\sum_{i=1}^{N_{x}}A_{ij}^{NOM}$ is the mean value
over the elements in the $j^{th}$ column of the matrix $A^{NOM}$ and
$b_{i}=(N_{y})^{-1}\sum_{j=1}^{N_{y}}B_{ji}^{NOM}$ is the mean value over the
elements in the $i^{th}$ column of the matrix $B^{NOM}$. Note that in order to
obtain (13) we have used Eqs. (8) and the following properties
$\sum_{i=1}^{N_{x}}\delta\mu_{i}=0$, $\sum_{j=1}^{N_{y}}\delta\nu_{j}=0$,
$\sum_{i=1}^{N_{x}}\delta a_{i}=0$, $\sum_{j=1}^{N_{y}}\delta b_{j}=0$,
$\sum_{i=1}^{N_{x}}\delta\pmb{x}_{i}=0$ and
$\sum_{j=1}^{N_{y}}\delta\pmb{y}_{j}=0$. Note also that the two matrices
$\tilde{A}$ and $\tilde{B}$ have sums over both their columns and rows equal
to zero.

In what follows, we will study stability of the system (13). If this system is
found to be stable, that implies that the trajectories (1) will tend to remain
close to the average solution (10) [see the definition of the variations in
(9)].

In what follows, we assume that for each connection from node $i$ in group
$\pmb{X}$ to node $j$ in group $\pmb{Y}$, there is a connection from node $j$
in group $\pmb{Y}$ to node $i$ in group $\pmb{X}$ and vice versa. We note that
we must still satisfy the assumption (which we imposed at the beginning of the
paper without loss of generality) that the sums over the entries in the rows
of both matrices $A^{NOM}$ and $B^{NOM}$ are equal to one. With these
considerations in mind, we take $B^{NOM}=e{A^{NOM}}^{T}$, where
$e=N_{y}/N_{x}$. It then follows $\tilde{B}=e\tilde{A}^{T}$.

We then rewrite Eq. (13) in vectorial form,

$\displaystyle\delta\dot{\pmb{x}}=[I_{N_{x}}\otimes\mathrm{D}\pmb{F}_{x}^{*}]\delta\pmb{x}+[\tilde{A}\otimes\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})]\delta\pmb{y}+\pmb{U}_{x}$
(15a)
$\displaystyle\delta\dot{\pmb{y}}=[I_{N_{y}}\otimes\mathrm{D}\pmb{G}_{y}^{*}]\delta\pmb{y}+[e\tilde{A}^{T}\otimes\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})]\delta\pmb{x}+\pmb{U}_{y}$
(15b)

By the substitution $\delta\pmb{y}\leftarrow\delta\pmb{y}/\sqrt{e}$, Eq. (15)
becomes,

$\displaystyle\delta\dot{\pmb{x}}=[I_{N_{x}}\otimes\mathrm{D}\pmb{F}_{x}^{*}]\delta\pmb{x}+[\sqrt{e}\tilde{A}\otimes\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})]\delta\pmb{y}+\pmb{U}_{x}$
(16a)
$\displaystyle\delta\dot{\pmb{y}}=[I_{N_{y}}\otimes\mathrm{D}\pmb{G}_{y}^{*}]\delta\pmb{y}+[\sqrt{e}\tilde{A}^{T}\otimes\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})]\delta\pmb{x}+\pmb{U}_{y}/\sqrt{e}$
(16b)

By stacking all perturbation vectors together in one vector
${\pmb{z}}=[\delta{\pmb{x}}^{T},\delta{\pmb{y}}^{T}]^{T}$ and by introducing
the vector $\pmb{U}=[\pmb{U}_{x}^{T},\pmb{U}_{y}^{T}/\sqrt{e}]^{T}$, Eq. (15)
can be rewritten as

$\begin{array}[]{c}\dot{\pmb{z}}=\begin{pmatrix}I_{N_{x}}\otimes\mathrm{D}\pmb{F}_{x}^{*}&0\\\
0&{0}_{N_{y}}\otimes\mathrm{D}\pmb{G}_{y}^{*}\end{pmatrix}\pmb{z}+\begin{pmatrix}0_{N_{x}}\otimes\mathrm{D}\pmb{F}_{x}^{*}&0\\\
0&I_{N_{y}}\otimes\mathrm{D}\pmb{G}_{y}^{*}\end{pmatrix}\pmb{z}+\\\ \\\
\sqrt{e}\begin{pmatrix}0&\tilde{A}\otimes\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})\\\
\tilde{A}^{T}\otimes\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})&0\end{pmatrix}\pmb{z}+\pmb{U}\par\end{array}$
(17)

We define the matrices,

$E^{x}=\begin{pmatrix}I_{N_{x}}&0\\\ 0&{0}_{N_{y}}\end{pmatrix}\quad
E^{y}=\begin{pmatrix}0_{N_{x}}&0\\\
0&{I}_{N_{y}}\end{pmatrix}\quad\tilde{M}=\begin{pmatrix}0_{n_{x}}&\tilde{A}\\\
\tilde{A}^{T}&{0}_{n_{y}}\end{pmatrix}$ (18)

and seek to compute the $N_{x}+N_{y}$-dimensional transformation matrix $T$
that simultaneously block diagonalizes the matrices in Eq. (18),

$T=\mathcal{SBD}(E^{x},E^{y},\tilde{M}),$ (19)

where here by $\mathcal{SBD}$ we indicate that the transformed matrices
$TE^{x}T^{T}$, $TE^{y}T^{T}$, and $T\tilde{M}T^{T}$ are in the same block-
diagonal form with blocks of minimal dimension [16, 17, 18, 19, 20]. To
compute the matrix $T$ we first need to find a matrix $P$ that commutes with
the symmetric matrices $E^{x}$, $E^{y}$, and $\tilde{M}$. The matrix $T$ has
for columns the eigenvectors of the matrix $P$ [17]. In order for the matrix
$P$ to commute with the diagonal matrices $E^{x}$ and $E^{y}$, it needs to be
in the following block-diagonal form,

$P=\begin{pmatrix}P_{1}&0\\\ 0&P_{2}\end{pmatrix}$ (20)

where the square block $P_{1}$ has dimension $N_{x}$ and the square block
$P_{2}$ has dimension $N_{y}$. Finally, the matrix $P$ in (20) needs also to
commute with $\tilde{M}$. Then the following two relations must be satisfied,

$\displaystyle\tilde{A}P_{2}=$ $\displaystyle P_{1}\tilde{A}$ (21a)
$\displaystyle\tilde{A}^{T}P_{1}=$ $\displaystyle P_{2}\tilde{A}^{T},$ (21b)

A simple solution to the above set of equations is given by

$\displaystyle P_{1}=$ $\displaystyle\tilde{A}\tilde{A}^{T}$ (22a)
$\displaystyle P_{2}=$ $\displaystyle\tilde{A}^{T}\tilde{A}.$ (22b)

From Eq. (20) it follows that the matrix $T$ has the following block-diagonal
structure,

$T=\begin{pmatrix}T^{x}&0\\\ 0&T^{y},\end{pmatrix}$ (23)

where the columns of the matrix $T^{x}$ are the eigenvectors of the matrix
$\tilde{A}\tilde{A}^{T}$ and the columns of the matrix $T^{y}$ are the
eigenvectors of the matrix $\tilde{A}^{T}\tilde{A}$. As both matrices
$\tilde{A}\tilde{A}^{T}$ and $\tilde{A}^{T}\tilde{A}$ are symmetric, the
matrices $T^{x}$ and $T^{y}$ are orthogonal. Hence, $T$ is also an orthogonal
matrix.

We can then define the $N_{x}n_{x}+N_{y}n_{y}$-dimensional orthogonal matrix

$\hat{T}=(T^{x}\otimes I_{n_{x}})\oplus(T^{y}\otimes I_{n_{y}})$ (24)

and use it to block-diagonalize Eq. (17). We introduce
$\pmb{w}=\hat{T}\pmb{z}$, then we have:

$\dot{\pmb{w}}=\begin{pmatrix}I_{N_{x}}\otimes\mathrm{D}\pmb{F}_{x}^{*}&0\\\
0&I_{N_{y}}\otimes\mathrm{D}\pmb{G}_{y}^{*}\end{pmatrix}\pmb{w}+\sqrt{e}\begin{pmatrix}0&S_{1}\otimes\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})\\\
S_{1}^{T}\otimes\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})&0\end{pmatrix}\pmb{w}+T\pmb{U}$
(25)

where

$S_{1}={T^{x}}\tilde{A}{T^{y}}^{T}.$ (26)

As mentioned before, both $T^{x}$ and $T^{y}$ are orthogonal matrices, hence
${T^{x}}^{T}{T^{x}}=I_{N_{x}}$ and ${T^{y}}^{T}T^{y}=I_{N_{y}}$. From (26), we
see that $S_{1}$ is the $N_{x}$ by $N_{y}$ matrix whose entries
${S_{1}}_{kk}=s_{\kappa}$ are equal to the singular values of the matrix
$\tilde{A}$ (all the other entries are zero).

There are many possible numerical methods to compute the singular values of a
matrix [21]. We note here that the eigenvalues of the matrix $\tilde{M}$ in
(18) are equal to $\\{\pm s_{1},\pm s_{2},...,\pm
s_{r}\\}\cup\\{\underbrace{0,0,...,0}_{\text{$N_{x}-N_{y}$ times}}\\}$ [3].
From a computational point of view the method of computing the singular values
from the eigenvalues of the matrix $\tilde{M}$ is more stable than direct
calculation of the matrix $S_{1}$ in (26).

We note that through the transformation $\pmb{w}=\hat{T}\pmb{z}$, each row of
the matrix $T^{x}$, say $T^{x\kappa}$, is associated with either a singular
value $s_{\kappa}$ or a zero value; analogously, each row of the matrix
$T^{y}$, say $T^{y\kappa}$, is associated with a singular value $s_{\kappa}$.

We then introduce the matrix:

$Q=\left(\begin{array}[]{cc}0&S_{1}\\\
S_{1}^{T}&0\end{array}\right)=\left(\begin{array}[]{ccccccc|cccc}&&&&&&&s_{1}&0&\cdots&0\\\
&&&&&&&0&s_{2}&\cdots&0\\\ &&&&&&&0&0&\ddots&0\\\
&&&\makebox(0.0,0.0)[]{\text{\huge 0}}&&&&0&0&\cdots&s_{r}\\\
&&&&&&&0&0&\cdots&0\\\ &&&&&&&\vdots&\vdots&\ddots&\vdots\\\
&&&&&&&0&0&\cdots&0\\\ \hline\cr s_{1}&0&0&0&0&\cdots&0&&&&\\\
0&s_{2}&0&0&0&\cdots&0&&&&\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&&&\makebox(0.0,0.0)[]{\text{\huge
0}}&\\\ 0&0&0&s_{r}&0&\cdots&0&&&&\\\ \end{array}\right)$ (27)

Based on our assumption that $N_{x}\geq N_{y}$, we have that $r=N_{y}$. We
call $d=(N_{x}-N_{y})$, $d+r=N_{x}$.

Now by permuting rows and columns of matrix $Q$, we obtain the block diagonal
matrix,

$\tilde{Q}=\left(\begin{array}[]{ccccc|ccc}\tilde{Q}_{1}&0&0&\cdots&0&0&\cdots&0\\\
0&\tilde{Q}_{2}&0&\cdots&0&0&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&\tilde{Q}_{r}&0&\cdots&0\\\ \hline\cr
0&0&0&\cdots&0&0&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&0&0&\cdots&0\\\ \end{array}\right)$ (28)

where each $\tilde{Q}_{i}\quad i=1,2,\cdots,r$ is a $2$ by $2$ block,

$\tilde{Q}_{\kappa}=\begin{pmatrix}0&s_{\kappa}\\\
s_{\kappa}&0\end{pmatrix}\quad\kappa=1,2,\cdots,r.$ (29)

One of the blocks $\kappa=1$ corresponds to motion parallel to synchronization
manifold and the other $r-1$ blocks, $i=2,...,r$ correspond to motion
transverse to the synchronization manifold [3]. By applying the same
permutation to the first term on the right hand side of Eq. (25), for each
$\kappa=1,..,r$ we can write:

$\dot{\pmb{p}}_{\kappa}=\Biggl{[}\begin{pmatrix}\mathrm{D}\pmb{F}_{x}^{*}&0\\\
0&\mathrm{D}\pmb{G}_{y}^{*}\end{pmatrix}+\begin{pmatrix}0&\lambda_{\kappa}\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}(\bar{\pmb{y}})\\\
\lambda_{\kappa}\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}(\bar{\pmb{x}})&0\end{pmatrix}\Biggr{]}\pmb{p}_{\kappa}+\Biggl{[}\begin{matrix}{\Big{(}\sum_{i}^{N_{x}}{T_{i}}^{x\kappa}\delta
a_{i}\Big{)}}\mathrm{D}\pmb{F}_{u}^{*}H(\bar{y})+{\Big{(}\sum_{i}^{N_{x}}{T_{i}}^{x\kappa}\delta\mu_{i}\Big{)}}\mathrm{D}\pmb{F}_{\mu}^{*}\\\
{\Big{(}\frac{\sum_{i}^{N_{y}}{T_{i}}^{y\kappa}\delta
b_{i}}{\sqrt{e}}\Big{)}}\mathrm{D}\pmb{G}_{u}^{*}L(\bar{x})+{\Big{(}\frac{\sum_{i}^{N_{y}}{T_{i}}^{y\kappa}\delta\nu_{i}}{\sqrt{e}}\Big{)}}\mathrm{D}\pmb{G}_{\nu}^{*}\end{matrix}\Biggl{]},$
(30)

where $\lambda_{\kappa}=\sqrt{e}s_{\kappa}$, $T^{x\kappa}_{i}$
($T^{x\kappa}_{i}$) is entry $i$ of row $T^{x\kappa}$ ($T^{x\kappa}$).
Analogously, to the right-lower block of the matrix $\tilde{Q}$ (Eq. (28)), we
can associate $d$ equations,

$\dot{\pmb{q}}_{\kappa}=\mathrm{D}\pmb{F}_{x}^{*}\pmb{q}_{\kappa}+{\Big{(}\sum_{i}T^{x\kappa}_{i}\delta
a_{i}\Big{)}}\mathrm{D}\pmb{F}_{u}^{*}H(\bar{y})+{\Big{(}\sum_{i}T^{x\kappa}_{i}\delta\mu_{i}\Big{)}}\mathrm{D}\pmb{F}_{\mu}^{*},\quad\kappa=1,2,\cdots,d.$
(31)

The homogeneous part of Eqs. (30) and (31) coincides with the reduction
obtained in Ref. [3].

The vector $\pmb{w}$ can be written as
$[{\pmb{p}}_{1}^{T},{\pmb{p}}_{2}^{T},...,{\pmb{p}}_{r}^{T},\pmb{q}_{1}^{T},\pmb{q}_{2}^{T},...,\pmb{q}_{p}^{T}]^{T}$,
where ${\pmb{p}}_{1}^{T}$ is a parallel perturbation and all the remaining
perturbations are transverse.

We now define the following master stability functions associated with Eqs.
(30) and (31),

$\mathcal{M}_{\kappa}(s_{\kappa},\alpha_{\kappa},\beta_{\kappa},\gamma_{k},\delta_{\kappa})=\lim_{\tau\to\infty}\sqrt{\tau^{-1}\int_{0}^{\tau}\|{\pmb{p}}_{\kappa}(t)\|^{2}dt},$
(32a)
$\mathcal{M}_{0}(0,\alpha_{\kappa},\beta_{\kappa})=\lim_{\tau\to\infty}\sqrt{\tau^{-1}\int_{0}^{\tau}\|\pmb{q}_{\kappa}(t)\|^{2}dt},$
(32b)

where $\alpha_{\kappa}={(\sum_{i}T^{x\kappa}_{i}\delta a_{i})}$,
$\beta_{\kappa}={(\sum_{i}T^{x\kappa}_{i}\delta\mu_{i})}$,
$\gamma_{\kappa}={(\sum_{i}T^{y\kappa}_{i}\delta b_{i})}/\sqrt{e}$,
$\delta_{\kappa}={(\sum_{i}T^{y\kappa}_{i}\delta\nu_{i})}/\sqrt{e}$.

Moreover, as the system of equations (30) and (31) is linear, in case of
stability, the master stability functions (32) can be approximated as follows:

$\displaystyle\mathcal{M}_{\kappa}=$ $\displaystyle c_{\alpha
1}(s)|\alpha_{\kappa}|+c_{\beta 1}(s)|\beta_{\kappa}|+c_{\gamma
1}(s)|\gamma_{\kappa}|+c_{\delta 1}(s)|\delta_{\kappa}|,$ (33a)
$\displaystyle\mathcal{M}_{0}=$ $\displaystyle c_{\alpha
2}|\alpha_{\kappa}|+c_{\beta 2}|\beta_{\kappa}|,$ (33b)

where $c_{\alpha 1},c_{\beta 1},c_{\gamma 1},c_{\delta 1},c_{\alpha
2},andc_{\beta 2}$ are appropriate coefficients (see also [11]). Note that
Eqs. (33) relate the master stability functions with the parameter mismatches
through the coefficients $\alpha_{\kappa},\beta_{\kappa},\gamma_{\kappa}$ and
$\delta_{\kappa}$ defined above.

Next, following Refs. [10, 11], we will show how the low-dimensional approach
[and the master stability functions defined in (32)] can be used to quantify
the synchronization error observed in the high-dimensional system. The main
motivation for performing this analysis is that typically dealing with the
low-dimensional systems (30) and (31) is computationally more convenient than
with the high-dimensional system (13) (see also Refs. [10, 11]).

We define the synchronization error,

$E(t)=\sqrt{\sum_{i=1}^{N_{x}}\delta\pmb{x}_{i}^{T}(t)\delta\pmb{x}_{i}(t)+\sum_{i=1}^{N_{y}}\delta\pmb{y}_{i}^{T}(t)\delta\pmb{y}_{i}(t)}=\sqrt{\sum_{\kappa=2}^{r}{\pmb{p}}_{\kappa}^{T}(t){\pmb{p}}_{\kappa}(t)+\sum_{\kappa=1}^{d}\pmb{q}_{\kappa}^{T}(t)\pmb{q}_{\kappa}(t)}$
(34)

Assuming stability, the time average of the synchronization error can be
computed in terms of the master stability functions,

$<E(t)>_{t}=\sqrt{\sum_{\kappa=2}^{r}{\mathcal{M}_{\kappa}}^{2}+\sum_{\kappa=1}^{d}{\mathcal{M}_{0}}^{2},}$
(35)

where with the symbol $<...>_{t}$, we indicate a time-average. Not that
$\mathcal{M}_{0}(0,\alpha_{\kappa},\beta_{\kappa})$ still depends on $\kappa$
through both arguments $\alpha_{\kappa}$ and $\beta_{\kappa}$, see Eq. (32b)

###  II.1 Spectrum of the matrix $\tilde{M}$

We have already stated that the singular values of the matrix $\tilde{A}$ can
be computed from the eigenvalues of the matrix $\tilde{M}$. We want to show
that there is a direct relation between the eigenvalues of the matrix
$\tilde{M}$ and those of the nominal matrix

$M^{NOM}=\left(\begin{array}[]{cc}0&A^{NOM}\\\
{A^{NOM}}^{T}&0\end{array}\right),$ (36)

where here we retain the assumption (with no loss of generality) that the sum
of the entries in the rows of the matrix $A^{NOM}$ is equal to $1$. This also
implies that the sum of the entries in the columns of the matrix $A^{NOM}$ is
equal to $1/e$. We know from Ref. [3] that the spectrum of the matrix
$M^{NOM}$ is characterized by the following properties: (i) the spectrum is
symmetric with respect to the real and imaginary axes, (ii) at least
$d=(N_{x}-N_{y})$ eigenvalues are equal to zero, (iii) two eigenvalues are
equal to $+\sqrt{e^{-1}}$ and $-\sqrt{e^{-1}}$ and (iv) the eigenvectors of
the matrix $M^{NOM}$,
$\pmb{v}^{\kappa}=[\pmb{v}^{x\kappa}\quad\pmb{v}^{y\kappa}]^{T}$ can be of
either one of two different types:

  * •

Type I eigenvectors associated with eigenvalues $\pm\sqrt{e^{-1}}$,
$\pmb{v}^{1}=[\underbrace{1,1,...,1}_{\text{$N_{x}$
times}},\underbrace{0,0,...,0}_{\text{$N_{y}$ times}}]$ and
$\pmb{v}^{2}=[\underbrace{0,0,...,0}_{\text{$N_{x}$
times}},\underbrace{1,1,...,1}_{\text{$N_{y}$ times}}]$. These correspond to
modes parallel to the synchronization manifold.

  * •

Type II eigenvectors associated with the remaining $(N_{x}+N_{y}-2)$
eigenvalues, for which $\sum_{l=1}^{N_{x}}v^{x\kappa}_{l}=0$ and
$\sum_{l=1}^{N_{y}}v^{y\kappa}_{l}=0$, $\kappa=3,...,N_{x}+N_{y}$. These
correspond to modes orthogonal to the synchronization manifold.

In the first part of this Section we have seen that, in case of parametric
mismatches, the master stability function (32) and the extended master
stability function (33) depend on the eigenvalues of the matrix $\tilde{M}$
rather than those of the nominal matrix $M^{NOM}$. We want to show that the
spectrum of the matrix $\tilde{M}$ is closely related to that of the matrix
$M^{NOM}$, namely,

$M^{NOM}\pmb{v}^{\kappa}=\varsigma^{\kappa}\pmb{v}^{\kappa}\Leftrightarrow\tilde{M}\pmb{v}^{\kappa}=\vartheta^{k}\pmb{v}^{\kappa},$
(37)

where (a) $\vartheta^{1}=\vartheta^{2}=0$ and (b)
$\vartheta^{\kappa}=\varsigma^{\kappa}$, $\kappa=3,..,N_{x}+N_{y}$.

To prove proberty (a) we just need to observe that from the structure of the
matrix $\tilde{M}$ in (18) and the property that the sums over the entries in
the rows and in the columns of the matrix $\tilde{A}$ is equal to zero, it
follows that $\tilde{M}\pmb{v}^{\kappa}=\pmb{0}$ for $\kappa=1,2$. Hence,
$\pmb{v}^{1}$ and $\pmb{v}^{2}$ are still eigenvectors for the matrix
$\tilde{M}$ but with associated eigenvalue $0$.

To prove property (b), we see that the left eigenvalue equation for the matrix
$\tilde{M}$ corresponds to the following equations:
$\pmb{v}^{x\kappa^{T}}\tilde{A}=\varsigma^{\kappa}\pmb{v}^{y\kappa^{T}}$ and
$\pmb{v}^{y\kappa^{T}}\tilde{A}^{T}=\varsigma^{\kappa}\pmb{v}^{x\kappa^{T}}$.
We then note that from Eq. (14) we can write that the matrix
$\tilde{A}=(A^{NOM}-\Delta)$, where the entries over column $j$ of the matrix
$\Delta$ are all the same and equal to $a_{j}$, $j=1,...,N_{y}$. Then for
$\kappa=3,...,N_{x}+N_{y}$,
$\pmb{v}^{x\kappa^{T}}\tilde{A}=\pmb{v}^{x\kappa^{T}}({A}^{NOM}-\Delta)=\varsigma^{\kappa}\pmb{v}^{y\kappa^{T}}$
since $\pmb{v}^{x\kappa^{T}}\Delta=\pmb{0}$ from the property that the sum
over the entries of the vector $\pmb{v}^{x\kappa}$ is equal to zero, for
$\kappa=3,...,N_{x}+N_{y}$. Hence, for $\kappa=3,...,N_{x}+N_{y}$,
$\pmb{v}^{\kappa}$ is still an eigenvector for the matrix $\tilde{M}$ with
associated eigenvalue $\varsigma^{\kappa}$.

###  II.2 Example

Consider the example of the bipartite network with two groups,
$N_{x}=12,N_{y}=6$ shown in Fig. 1.

Figure 1: A bipartite network with $N_{x}=12$ nodes in the $X$ group, shown in
light green, and $N_{y}=6$ nodes in the $Y$ group, shown in dark green. The
network is bipartite so it contains just inter group connections.

The matrix $A^{NOM}$ for this network is

$A^{NOM}=0.2\begin{pmatrix}0&1&1&1&1&1\\\ 0&1&1&1&1&1\\\ 1&0&1&1&1&1\\\
1&0&1&1&1&1\\\ 1&1&0&1&1&1\\\ 1&1&0&1&1&1\\\ 1&1&1&0&1&1\\\ 1&1&1&0&1&1\\\
1&1&1&1&0&1\\\ 1&1&1&1&0&1\\\ 1&1&1&1&1&0\\\ 1&1&1&1&1&0\end{pmatrix},$ (38)

$r=d=6$. By calculating the transformation matrix $T$ and following Eq. (25),
the matrix $\tilde{Q}$ for this network is

$\tilde{Q}=\frac{0.2}{\sqrt{1/2}}\left(\begin{array}[]{cccccccccccccccccc}\cline{1-2}\cr\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\vrule\lx@intercol\hfil 0\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\cline{1-4}\cr 0&0&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil
1\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&\vrule\lx@intercol\hfil 1\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\cline{3-6}\cr 0&0&0&0&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil
1\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\vrule\lx@intercol\hfil 1\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0\\\ \cline{5-8}\cr
0&0&0&0&0&0&\vrule\lx@intercol\hfil 0\hfil\lx@intercol&\lx@intercol\hfil
1\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&\vrule\lx@intercol\hfil 1\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0\\\ \cline{7-10}\cr
0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil 0\hfil\lx@intercol&\lx@intercol\hfil
1\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil 1\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0\\\ \cline{9-12}\cr
0&0&0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil
1\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil
1\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0\\\ \cline{11-12}\cr
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\end{array}\right)$ (39)

where each one of the $r$ two-dimensional $Q_{i}$ blocks is shown in a box.
The first box ($Q_{1}$) corresponds to the parallel motion, with associated
singular value $s_{1}=0$. The other blocks correspond to the transverse
motions with associated singular values
$s_{2}=s_{3}=...=s_{6}={0.2}/{\sqrt{1/2}}$. By applying the same permutation
applied to the matrix $Q$ to the first term on the right hand side of Eq.
(25), we have:

$\left(\begin{array}[]{cccccccccccccccccc}\cline{1-2}\cr\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil{0}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\cline{1-4}\cr
0&0&\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil{0}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
\cline{3-6}\cr
0&0&0&0&\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil{0}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0&0&0\\\
\cline{5-8}\cr
0&0&0&0&0&0&\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil{0}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0&0&0\\\
\cline{7-10}\cr
0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil{0}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0&0&0\\\
\cline{9-12}\cr
0&0&0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil\mathrm{D}\pmb{F}_{x}^{*}\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&\vrule\lx@intercol\hfil{0}\hfil\lx@intercol&\lx@intercol\hfil\mathrm{D}\pmb{G}_{y}^{*}\hfil\lx@intercol\vrule\lx@intercol&0&0&0&0&0&0\\\
\cline{11-12}\cr
0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\mathrm{D}\pmb{F}_{x}^{*}\\\
\end{array}\right)$ (40)

We further note that in case there are no mismatches, that is, $\delta
a_{i}=0$, $\delta b_{j}=0$, $\delta\mu_{i}=0$, and $\delta\nu_{j}=0$, the low-
dimensional system (30) reduces to the unforced system

$\displaystyle\dot{\tilde{\pmb{x}}}(t)$
$\displaystyle=\mathrm{D}\pmb{F}_{x}^{*}\tilde{\pmb{x}}(t)+\lambda_{\kappa}\mathrm{D}\pmb{F}_{u}^{*}\mathrm{D}\pmb{H}_{y}^{*}\tilde{\pmb{y}}(t),$
(41a) $\displaystyle\dot{\tilde{\pmb{y}}}(t)$
$\displaystyle=\mathrm{D}\pmb{G}_{y}^{*}\tilde{\pmb{y}}(t)+\lambda_{\kappa}\mathrm{D}\pmb{G}_{u}^{*}\mathrm{D}\pmb{L}_{x}^{*}\tilde{\pmb{x}}(t),$
(41b)

which is the low dimensional solution found in Ref. [3]. Note that Eqs. (25)
is a system of linear time-varying nonhomogeneous differential equations and
is characterized by the same stability range as the associated homogeneous
system (41) [22]. The condition for stability is that the maximum Lyapunov
exponents of (41) are negative for $\kappa=2,...,r$, excluding
$\lambda_{1}=0$, which is associated with perturbations tangent to the
synchronization manifold [3].

It is important to emphasize that Eqs. (41) depend on $\bar{x}(t)$ and
$\bar{y}(t)$, which are averaged trajectories over all the systems in the
network. In a large network, calculating $\bar{x}(t)$ and $\bar{y}(t)$ may be
computationally expensive, as it requires full integration of $(N_{x}+N_{y})$
individual systems, see Eq. (11). However, for practical purposes,
$\bar{x}(t)$ and $\bar{y}(t)$ in (41) can be replaced by $x_{s}(t)$ and
$y_{s}(t)$ obeying Eq. (3) (see also Ref. [11]).

We conclude that in the presence of small parametric mismatches, we can still
use (41) to determine whether the synchronous solution is stable. That also
explains the excellent agreement found between the stability region predicted
by the theory and experimentally observed in Ref. [5]. While this is not a
surprising result, as it is expected that linear stability of a smooth
dynamical system is not affected by small parametric variations, by using this
approach we will be able to quantify the synchronization error as a function
of the parameter mismatches, based on the low-dimensional reduction (30) and
(31).

We now consider the following set of equations,

$\displaystyle\dot{x}_{i1}=$
$\displaystyle-\sigma_{i}{x_{i1}}+\sigma_{i}\sum_{j=1}^{N_{y}}A_{ij}{y_{j}},$
(42a) $\displaystyle\dot{x}_{i2}=$ $\displaystyle
x_{i1}\sum_{j=1}^{N_{y}}A_{ij}{y_{j}}-\xi x_{i2}$ (42b)

$i=1,...,N_{x}$.

$\dot{y}_{j}=-y_{j}+\sum_{i=1}^{N_{x}}B_{ji}{x}_{i1}\Bigl{(}\rho-\sum_{i=1}^{N_{x}}B_{ji}x_{i2}\Bigr{)},$
(43)

$j=1,...,N_{y}$.

We set the nominal parameters as follows: $\sigma_{i}=\sigma=10$,
$i=1,...,N_{x}$, $\rho=28$, and $\xi=2$. If $\sum_{j}A_{ij}=1$ and
$\sum_{i}B_{ji}=1$, then the synchronous dynamics corresponds to that of the
chaotic Lorenz system [23]. Given the set of equations (42) and (43), Eqs.
(41) become

$\displaystyle\dot{\tilde{x}}_{1}=$
$\displaystyle-\sigma\tilde{x}_{1}+\sigma\lambda_{\kappa}\tilde{y},$
$\displaystyle\dot{\tilde{x}}_{2}=$
$\displaystyle\bar{y}\tilde{x}_{1}+\bar{x}_{1}\lambda_{\kappa}\tilde{y}-\xi{\tilde{x}}_{2},$
$\displaystyle\dot{\tilde{y}}=$
$\displaystyle-\tilde{y}+(\rho-\bar{x}_{2})\lambda_{\kappa}\tilde{x}_{1}-\bar{x}_{1}\lambda_{\kappa}\tilde{x}_{2}.$

The maximum Lyapunov exponent of this system is plotted versus
$\lambda_{\kappa}$ in Fig. 4, from which we see that the synchronous solution
is stable for $\lambda$ approximately less than $0.435$.

Figure 2: Maximum Lyapunov Exponent for the system (41) as a function of the
parameter $\lambda_{k}$. The dashed line is the $0$-ordinate line.

Figure 3 shows $c_{\alpha_{1}}$, $c_{\beta_{1}}$, $c_{\gamma}$, and
$c_{\delta}$ as functions of $\lambda_{\kappa}$. The two constant values
$c_{\alpha_{2}}=51.7524$ and $c_{\beta_{2}}=0.8002$. Using this information
and Eqs. (33), we can compute the master stability function as a function of
$\lambda_{\kappa}$. This in turn allows us to approximate the synchronization
error, using Eq. (35) from knowledge of $\alpha_{\kappa}$, $\beta_{\kappa}$,
$\gamma_{\kappa}$, and $\delta_{\kappa}$, as long as all the
$\lambda_{\kappa}$’s are in the region of stability,
$0\leq\lambda_{\kappa}\leq 0.435$. We wish to emphasize that this result is
general and applies to any bipartite network with bidirectional connections,
whose time-evolution obeys Eqs. (42) and (43).

Figure 3: $c_{\alpha_{1}}$, $c_{\beta_{1}}$, $c_{\gamma}$, and $c_{\delta}$
versus $\lambda_{\kappa}$. The vertical lines highlight the values of the
$c_{\alpha_{1}}$, $c_{\beta_{1}}$, $c_{\gamma}$ and $c_{\delta}$ for
$\lambda_{\kappa}={0.2}.$

In order to validate the theory, we consider the set of Eqs. (42) and (43) for
the network topology in Fig. 1 [$A^{NOM}$ in Eq. (38)] and study the following
two cases:

  1. (I)

Small mismatches affecting the couplings, namely
$A_{ij}=A_{ij}^{NOM}(1+5\times 10^{-5}\phi_{ij})$, where $\phi_{ij}$ is a
random number from a standard normal distribution.

  2. (II)

Small mismatches affecting both the couplings and the individual system
parameters, namely $A_{ij}=A_{ij}^{NOM}(1+5\times 10^{-5}\phi_{ij})$ and
$\sigma_{i}$ in Eq. (42a) is given by $\sigma_{i}=10(1+10^{-6}\psi_{i})$,
where both $\phi_{ij}$ and $\psi_{i}$ are random numbers from a standard
normal distribution.

All the transverse singular values of the matrix $A^{NOM}$ in Eq. (38) are
equal to $s_{\kappa}={0.2}/{\sqrt{1/2}}$. To these corresponds only one value
of $\lambda_{\kappa}=0.2$. Table 1 compares the synchronization error obtained
from integration of the full high-dimensional system (42) and (43) [left hand
side of Eq. (32)] and the estimated synchronization error using the master
stability function [right hand side of Eq. (32)], for both cases described
above.

Table 1: The synchronization error observed in the high-dimensional system and the estimated synchronization error using the master stability function defined in Eq. (32) for two different cases.  Case |  $<E(t)>_{t}$ |  $\sqrt{\sum_{\kappa=2}^{r}{\mathcal{M}_{\kappa}}^{2}+\sum_{\kappa=1}^{d}{\mathcal{M}_{0}}^{2}}$
---|---|---
(I) | $1.529\times 10^{-3}$ | $1.448\times 10^{-3}$
(II) | $1.609\times 10^{-3}$ | $1.614\times 10^{-3}$

##  III Couplings among members of the same group

So far we have only considered coupling from the $\pmb{X}$ group to the
$\pmb{Y}$ group and vice versa. While in general group synchronization is
possible in the presence of intra group connections, stability of the group
synchronous solution when also intra group connections are present has not
been fully elucidated. This problem was investigated in Ref. [4], where
conditions were presented for the dimensionality reduction of the stability
problem. Here we will extend some of the results in Ref. [4] and we will show
that the stability of group synchronization can be reduced in a low-
dimensional form for a broader class of networks.

Consider two groups of coupled oscillators, described by the following
equations:

$\displaystyle\dot{\pmb{x}}_{i}$
$\displaystyle=F(\pmb{x}_{i})+\sum_{j=1}^{N_{y}}A_{ij}H(\pmb{y}_{j})+\sum_{\ell=1}^{N_{x}}C_{i\ell}R(\pmb{x}_{\ell}),\quad
i=1,\ldots,N_{x},$ (44a) $\displaystyle\dot{\pmb{y}}_{j}$
$\displaystyle=G(\pmb{y}_{j})+\sum_{i=1}^{N_{x}}B_{ji}L(\pmb{y}_{i})+\sum_{\ell=1}^{N_{y}}D_{j\ell}S(\pmb{y}_{\ell}),\quad
j=1,\ldots,N_{y},$ (44b)

where the output functions $R:\mathbf{R}^{n_{x}}\rightarrow\mathbf{R}^{n_{x}}$
and $S:\mathbf{R}^{n_{y}}\rightarrow\mathbf{R}^{n_{y}}$ and the $N_{x}$-square
matrix $C$ and the $N_{y}$-square matrix $D$ define the intra group
connections of group $\pmb{X}$ and group $\pmb{Y}$, respectively.

If $\sum_{j=1}^{N_{y}}A_{ij}=a\neq 0$ and
$\sum_{\ell=1}^{N_{x}}C_{i\ell}=c\neq 0$ independent of $i$ and
$\sum_{i}^{N_{x}}B_{ji}=b\neq 0$ and $\sum_{\ell=1}^{N_{y}}D_{j\ell}=d\neq 0$
independent of $j$, a synchronous solution exists
$x_{1}=x_{2}=...=x_{N_{x}}=x_{s}$, $y_{1}=y_{2}=...=y_{N_{y}}=y_{s}$, obeying,

$\displaystyle\dot{\pmb{x}}_{s}$
$\displaystyle=F(\pmb{x}_{s})+aH(\pmb{y}_{s})+cR(\pmb{x}_{s})$ (45a)
$\displaystyle\dot{\pmb{y}}_{s}$
$\displaystyle=G(\pmb{y}_{s})+bL(\pmb{x}_{s})+dS(\pmb{y}_{s}).$ (45b)

Note that by appropriately rescaling the functions
$\pmb{H},\pmb{L},\pmb{R},\pmb{S}$ it is always possible to set $a=b=c=d=1$.
Therefore, without loss of generality, in what follows we will proceed under
this assumption.

In order to analyze stability, we linearize (44) about the synchronous
solution,

$\displaystyle\delta\dot{\pmb{x}}_{i}$
$\displaystyle=DF(\pmb{x}_{s})\delta\pmb{x}_{i}+\sum_{j=1}^{N_{y}}A_{ij}DH(\pmb{y}_{s})\delta\pmb{y}_{j}+\sum_{\ell=1}^{N_{x}}C_{i\ell}DR(\pmb{x}_{s})\delta\pmb{x}_{\ell},\quad
i=1,...,N_{x}$ (46a) $\displaystyle\delta\dot{\pmb{y}}_{j}$
$\displaystyle=DG(\pmb{y}_{s})\delta\pmb{y}_{j}+\sum_{i=1}^{N_{y}}B_{ji}DL(\pmb{x}_{s})\delta\pmb{x}_{i}+\sum_{\ell=1}^{N_{y}}D_{j\ell}DS(\pmb{y}_{s})\delta\pmb{y}_{\ell},\quad
j=1,...,N_{y}.$ (46b)

In [4] it was shown that Eq. (46) can be reduced in a low-dimensional form,
provided that the matrix $N=\begin{bmatrix}C&0\\\ 0&D\end{bmatrix}$ has the
same set of eigenvectors as the matrix $M=\begin{bmatrix}0&A\\\
B&0\end{bmatrix}$. It was thus concluded in [4] that if this condition is
satisfied, in a similar way to the case of bipartite topologies, stability of
the synchronous solution can be reduced in the following low dimensional form,

$\displaystyle\dot{\bar{\pmb{x}}}^{\kappa}$
$\displaystyle=DF(\pmb{x}_{s})\bar{\pmb{x}}^{\kappa}+{\lambda^{\kappa}}DH(\pmb{y}_{s})\bar{\pmb{y}}^{\kappa}+\eta^{\kappa}DR(\pmb{x}_{s})\bar{\pmb{x}}^{\kappa}$
(47a) $\displaystyle\dot{\bar{\pmb{y}}}^{\kappa}$
$\displaystyle=DG(\pmb{y}_{s})\bar{\pmb{y}}^{\kappa}+\lambda^{\kappa}DL(\pmb{x}_{s})\bar{\pmb{x}}^{\kappa}+\eta^{\kappa}DS(\pmb{y}_{s})\bar{\pmb{y}}^{\kappa},$
(47b)

$\kappa=1,...,N_{x}+N_{y}$, where $\lambda^{\kappa}$ and $\eta^{\kappa}$ are
the eigenvalues associated with the same eigenvector of the matrices $M$ and
$N$, respectively. By construction, there are always two $\kappa$ eigenvectors
which are associated with perturbations tangent to the synchronization
manifold, and the corresponding pairs of eigenvalues
$(\lambda^{\kappa},\eta^{\kappa})$ do not need to be considered in order to
assess stability. Without loss of generality, we label $\kappa=1$ and
$\kappa=2$ the two tangent eigenmodes and $\kappa=3,...,N_{x}+N_{y}$ the
remaining transverse eigenmodes.

In what follows we extend the results in [4] and show that a reduction in a
low dimensional form is possible for any pair of matrices

$\displaystyle C=C^{\prime}+J^{x},$ (48a) $\displaystyle D=D^{\prime}+J^{y}.$
(48b)

where (i) the matrix $\begin{bmatrix}C^{\prime}&0\\\
0&D^{\prime}\end{bmatrix}$ has the same set of eigenvectors as the matrix $M$
and (ii) $J^{x}$ and $J^{y}$ are any two $N_{x}$-square and $N_{y}$-square
matrices whose columns are composed of entries that are all the same. This
form of the matrices $J^{x}$ and $J^{y}$ corresponds to a particular coupling
configuration for which for any pair of nodes $(i,j)$ in either group
$\pmb{X}$ or group $\pmb{Y}$, the coupling strength from node $i$ to node $j$
is a function of $i$ but not of $j$.

Using the ansatz
$\delta\pmb{x}_{i}(t)=p_{i}^{\kappa}\hat{\pmb{x}}^{\kappa}(t)$ and
$\delta\pmb{y}_{j}(t)=q_{j}^{\kappa}\hat{\pmb{y}}^{\kappa}(t)$ in (46), we
obtain,

$\displaystyle{p_{i}^{\kappa}}\dot{\hat{\pmb{x}}}^{\kappa}$
$\displaystyle=DF(\pmb{x}_{s}){p_{i}^{\kappa}}\hat{\pmb{x}}^{\kappa}+\sum_{j=1}^{N_{y}}A_{ij}{q_{j}^{\kappa}}DH(\pmb{y}_{s})\hat{\pmb{y}}^{\kappa}+{\sum_{\ell=1}^{N_{x}}C_{i\ell}DR(\pmb{x}_{s})p_{\ell}^{\kappa}}\hat{\pmb{x}}^{\kappa}$
(49a) $\displaystyle{q_{j}^{\kappa}}\dot{\hat{\pmb{y}}}^{\kappa}$
$\displaystyle=DG(\pmb{y}_{s}){q_{j}^{\kappa}}\hat{\pmb{y}}^{\kappa}+\sum_{i=1}^{N_{x}}B_{ji}{p_{i}^{\kappa}}DL(\pmb{x}_{s})\hat{\pmb{x}}^{\kappa}+{\sum_{\ell=1}^{N_{y}}D_{j\ell}DS(\pmb{y}_{s})q_{\ell}^{\kappa}}\hat{\pmb{y}}^{\kappa}$
(49b)

We now consider the left eigenvalue equation,

$\pmb{\omega}^{\kappa^{T}}M=\lambda^{\kappa}\pmb{\omega}^{\kappa^{T}},$ (50)

where the eigenvectors
$\pmb{\omega}^{\kappa^{T}}=\Big{(}p_{1}^{\kappa^{\mathit{l}}},p_{2}^{\kappa^{\mathit{l}}},\dots,p_{N_{x}}^{\kappa^{\mathit{l}}},q_{1}^{\kappa^{\mathit{l}}},q_{2}^{\kappa^{\mathit{l}}},\dots,q_{N_{y}}^{\kappa^{\mathit{l}}}\Big{)}$.
By multiplying (49a) by $p_{i}^{\kappa^{\mathit{l}}}$, summing over $i$, and
dividing by $(\sum_{i}p_{i}^{\kappa^{\mathit{l}}}{p_{i}^{\kappa}})$ [by
multiplying (49b) by $q_{j}^{\kappa^{\mathit{l}}}$, summing over $j$, and
dividing by $(\sum_{j}q_{j}^{\kappa^{\mathit{l}}}{q_{j}^{\kappa}})$], we
obtain,

$\displaystyle\dot{\hat{\pmb{x}}}^{\kappa}$
$\displaystyle=DF(\pmb{x}_{s})\hat{\pmb{x}}^{\kappa}+{\lambda^{\kappa}}DH(\pmb{y}_{s})\hat{\pmb{y}}^{\kappa}+\eta^{\kappa}DR(\pmb{x}_{s})\hat{\pmb{x}}^{\kappa}+\frac{\Big{(}\sum_{\ell}\sum_{i}p_{i}^{\kappa^{\mathit{l}}}J^{x}_{i\ell}p_{\ell}^{\kappa}\Big{)}}{\Big{(}\sum_{i}{p_{i}^{\kappa^{\mathit{l}}}}{p_{i}^{\kappa}}\Big{)}}DR(\pmb{x}_{s})\hat{\pmb{x}}^{\kappa},$
(51a) $\displaystyle\dot{\hat{\pmb{y}}}^{\kappa}$
$\displaystyle=DG(\pmb{y}_{s})\hat{\pmb{y}}^{\kappa}+\lambda^{\kappa}DL(\pmb{x}_{s})\hat{\pmb{x}}^{\kappa}+\eta^{\kappa}DS(\pmb{y}_{s})\hat{\pmb{y}}^{\kappa}+\frac{\Big{(}\sum_{\ell}\sum_{j}q_{j}^{\kappa^{\mathit{l}}}J^{y}_{j\ell}q_{\ell}^{\kappa}\Big{)}}{\Big{(}\sum_{j}{q_{j}^{\kappa^{\mathit{l}}}}{q_{j}^{\kappa}}\Big{)}}DS(\pmb{y}_{s})\hat{\pmb{y}}^{\kappa},$
(51b)

$\kappa=1,...,N_{x}+N_{y}$.

We shall now show that $\sum_{i}p_{i}^{\kappa^{\mathit{l}}}J^{x}_{i\ell}=0$,
for $\kappa=3,...,(N_{x}+N_{y})$ in Eq. (51a) and
$\sum_{j}q_{j}^{\kappa^{\mathit{l}}}J^{y}_{j\ell}=0$, for
$\kappa=3,...,(N_{x}+N_{y})$ in Eq. (51b). Since by definition the matrices
$J^{x}$ and $J^{y}$ have columns whose entries are all the same, it suffices
to show that $\sum_{i}p_{i}^{\kappa^{\mathit{l}}}=0$, for
$\kappa=3,...,(N_{x}+N_{y})$ and $\sum_{j}q_{j}^{\kappa^{\mathit{l}}}=0$, for
$\kappa=3,...,(N_{x}+N_{y})$. We will prove that this latter property is
indeed satisfied by the $(N_{x}+N_{y}-2)$ eigenvectors that are associated
with the transverse eigenmodes of the matrix $M$.

Note that the left eigenvalue equation
$\pmb{\omega}^{\kappa^{T}}{M}=\lambda^{\kappa}\pmb{\omega}^{\kappa^{T}}$
implies the following two equations:
$\pmb{q}^{\kappa^{T}}{B}=\lambda^{\kappa}\pmb{p}^{\kappa^{T}}$ and
$\pmb{p}^{\kappa^{T}}{A}=\lambda^{\kappa}\pmb{q}^{\kappa^{T}}$ where the
vectors
$\pmb{p}^{\kappa^{T}}=\Big{(}p_{1}^{\kappa^{\mathit{l}}},p_{2}^{\kappa^{\mathit{l}}},\dots,p_{N_{x}}^{\kappa^{\mathit{l}}}\Big{)}$
and
$\pmb{q}^{\kappa^{T}}=\Big{(}q_{1}^{\kappa^{\mathit{l}}},q_{2}^{\kappa^{\mathit{l}}},\dots,q_{N_{y}}^{\kappa^{\mathit{l}}}\Big{)}$,
or equivalently,

$\displaystyle\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}{B}_{ji}$
$\displaystyle=\lambda^{\kappa}p_{i}^{\kappa^{\mathit{l}}}$ (52a)
$\displaystyle\sum_{i=1}^{N_{x}}p_{i}^{\kappa^{\mathit{l}}}{A}_{ij}$
$\displaystyle=\lambda^{\kappa}q_{j}^{\kappa^{\mathit{l}}}.$ (52b)

By summing Eq. (52a) over $i=1,...,N_{x}$ and Eq. (52b) over $j=1,...,N_{y}$,
we obtain

$\displaystyle\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}$
$\displaystyle=\lambda^{\kappa}\sum_{i=1}^{N_{x}}p_{i}^{\kappa^{\mathit{l}}}$
(53a) $\displaystyle\sum_{i=1}^{N_{x}}p_{i}^{\kappa^{\mathit{l}}}$
$\displaystyle=\lambda^{\kappa}\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}.$
(53b)

By plugging (53b) into (53a) we obtain
$\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}=(\lambda^{\kappa})^{2}\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}$,
which can be satisfied in two possible ways: either $\lambda^{\kappa}=\pm 1$,
which applies to the parallel eigenmodes $\kappa=1,2$ of the matrix $M$ or
$\sum_{j=1}^{N_{y}}q_{j}^{\kappa^{\mathit{l}}}=0$, which applies to the
transverse eigenmodes $\kappa=3,...,(N_{x}+N_{y})$ of the matrix $M$.
Analogously, by plugging (53a) into (53b), we can show that
$\sum_{i=1}^{N_{x}}p_{j}^{\kappa^{\mathit{l}}}=0$ for
$\kappa=3,...,(N_{x}+N_{y})$.

We can thus conclude that Eq. (51) reduces to (47) for
$\kappa=3,...,(N_{x}+N_{y})$ [while in general it does not reduce to (47) for
$\kappa=1,2$]. This implies that when intra group connections are in the
general form of Eq. (48), transverse stability is still determined by the low-
dimensional system (47).

Finally, we note that the analysis of parameter mismatches presented in Sec.II
can be extended to the case that small mismatches affect the intra group
connections. In particular, if the intra group connections are affected by
small mismatches but are in a form that approximately satisfies Eq. (48), then
the low-dimensional system (47) can still be used to predict stability of the
synchronous state.

###  III.1 Numerical Example

As a numerical example, we considered a small network for which $N_{x}=2$,
$N_{y}=3$, and

$A=\left(\begin{array}[]{ccc}w&1-w&0\\\ 0&1-w&w\\\ \end{array}\right),\quad
B=\left(\begin{array}[]{cc}1&0\\\ 1\over 2&1\over 2\\\ 0&1\\\
\end{array}\right).$ (54)

In Ref. [3] we showed that for this choice of the coupling matrices $A$ and
$B$, the eigenvalues $\lambda^{\kappa}$ were equal to $0,\pm\sqrt{w},\pm 1$.
Then, by setting ${w}=0.64$ stability is only affected by the pair of
eigenvalues $\lambda^{\kappa}$, $0$ and $\sqrt{w}=0.8$ [3]. In our numerical
experiments, we retained the set of equation (42) and replaced (43) by the
following:

$\dot{y}_{j}=-y_{j}+\sum_{i=1}^{N_{x}}B_{ji}{x}_{i1}\Bigl{(}\rho-\sum_{i=1}^{N_{x}}B_{ji}x_{i2}\Bigr{)}+\sum_{\ell=1}^{N_{y}}D_{j\ell}S(y_{\ell}),$
(55)

$j=1,...,N_{y}$, where the last term on the right-hand-side of (55) represents
coupling between systems of the same group $\pmb{Y}$ and we set the coupling
function $S(y_{\ell})=sy_{\ell}$, with $s$ being a variable scalar quantity.

Figure 4:  (a) Maximum Lyapunov Exponent (MLE) for the system (47) vs. $s$ for
$\lambda^{\kappa}=0$ and $\lambda^{\kappa}=0.8$. The dashed line corresponds
to the $0$-ordinate line. (b) Synchronization error $E_{y}$ vs. $s$ for the
two choices of the matrix $D$ in Eq. (56).

We considered the following two cases of intragroup coupling matrices $D$,

$\left(\begin{array}[]{ccc}1/3&1/3&1/3\\\ 1/3&1/3&1/3\\\
1/3&1/3&1/3\end{array}\right),\qquad\left(\begin{array}[]{ccc}0&1&0\\\
0&1&0\\\ 0&1&0\end{array}\right),$ (56)

which we label case 1 and case 2, respectively. For both cases 1 and 2, the
corresponding matrices $M$ and $N$ do not commute, so neither case can be
studied within the framework presented in Ref. [4]. Note that both choices of
$D$ above are matrices whose columns are composed of entries that are all the
same. Therefore, in both cases we expect stability to be described by the low-
dimensional system (47) with $\eta^{\kappa}=0$. Fig. 4(a) shows the Maximum
Lyapunov Exponent (MLE) for the system (47) as a function of the parameter $s$
for $\lambda^{\kappa}=0$ and $\lambda^{\kappa}=0.8$, from which we see that
stability is expected for $s\gtrsim 1.5$.

In Fig.4(b) we run full numerical simulations of the network dynamics [Eqs.
(42) and (55)] as we vary the parameter $s$. For each run, the network systems
are evolved from initial conditions that are close to the synchronization
manifold. The figure shows the final synchronization error,

$E_{y}=N_{y}^{-1}(\tau)^{-1}\int_{t_{1}}^{t_{1}+\tau}\sum_{i=1}^{N_{y}}\Bigl{(}y_{i}(t)-\bar{y}(t)\Bigr{)}^{2}dt,$
(57)

with $t_{1}$ a large-enough time past the initial transient dynamics and
$\tau$ a large-enough averaging window, versus $s$. As can be seen, perfect
agreement is attained with the low-dimensional predictions of Fig. 4(a).

##  IV Conclusions

The main contribution of this paper is to extend the analysis of the effects
of parameter mismatches on the stability of the network synchronous solution
[9, 10, 11, 13, 12, 14, 15] to the case of networks formed of systems of
different types. An important reference for our work is Ref. [10] which
introduced a master stability function to characterize the synchronization
error in networks with parameter mismatches.

Group synchronization was first studied in [3] and has been the subject of
both theoretical [4] and experimental [5] investigations. Group
synchronization has been observed and characterized when the systems in each
group are identical and the couplings between the systems satisfy specific
conditions. In this paper, we have defined a master stability function that
describes stability of the group synchronization solution in the presence of
mismatches on the individual parameters of the network oscillators and on the
couplings between the oscillators in the groups. This is relevant to
experimental realizations of group synchrony [5], for which parameter
mismatches are practically unavoidable. Our analysis applies to the case that
for each connection from node $i$ in group $\pmb{X}$ to node $j$ in group
$\pmb{Y}$, there is a connection from node $j$ in group $\pmb{Y}$ to node $i$
in group $\pmb{X}$ and vice versa. We have also extended the analysis
presented in Ref. [4] to study stability of the group synchronization solution
in the presence of intra group couplings. Our analysis has pointed out a
broader class of matrices describing intra group connectivity for which the
stability problem can be reduced in a low-dimensional form.

## Acknowledgement

The authors thank Prof. Kazuo Murota for insightful discussions on the subject
of $*$-algebra. This research is supported by NIH grant 1R21EB028489-01A1.

## References

  * Arenas et al. [2008] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, , and C. Zhou, Phys. Rep. 469, 93 (2008).
  * Pecora and Carroll [1998] L. Pecora and T. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
  * Sorrentino and Ott [2007] F. Sorrentino and E. Ott, Phys. Rev. E 76, 056114 (2007).
  * Dahms et al. [2012] T. Dahms, J. Lehnert, and E. Schöll, Physical Review E 86, 016202 (2012).
  * Williams et al. [2013] C. R. Williams, T. E. Murphy, R. Roy, F. Sorrentino, T. Dahms, and E. Schöll, Phys. Rev. Lett. 110, 064104 (2013).
  * Yang [2015] Y. F. Yang, in _Applied Mechanics and Materials_ (Trans Tech Publ, 2015), vol. 733, pp. 902–905.
  * Panahi et al. [2021] S. Panahi, I. Klickstein, and F. Sorrentino, arXiv preprint arXiv:2109.13792 (2021).
  * Della Rossa et al. [2020] F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino, Nature communications 11, 1 (2020).
  * Restrepo et al. [2004] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 69, 066215 (2004).
  * Sun et al. [2009] J. Sun, E. M. Bollt, and T. Nishikawa, Europhys. Lett. 85, 60011 (2009).
  * Sorrentino and Porfiri [2011] F. Sorrentino and M. Porfiri, Europhys. Lett. 93, 50002 (2011).
  * Sorrentino and Pecora [2016] F. Sorrentino and L. Pecora, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094823 (2016).
  * Cho [2019] Y. S. Cho, Physical Review E 99, 052215 (2019).
  * Acharyya and Amritkar [2012] S. Acharyya and R. Amritkar, EPL (Europhysics Letters) 99, 40005 (2012).
  * Acharyya and Amritkar [2015] S. Acharyya and R. Amritkar, Physical Review E 92, 052902 (2015).
  * Shapiro [1979] H. Shapiro, Linear Algebra and its applications 25, 129 (1979).
  * Maehara and Murota [2011] T. Maehara and K. Murota, SIAM J. Matrix Anal. Appl. 33, 605 (2011).
  * Murota et al. [2010] K. Murota, Y. Kanno, M. Kojima, and S. Kojima, Jpn. J. Ind. Appl. Math. 27, 125 (2010).
  * Irving and Sorrentino [2012] D. Irving and F. Sorrentino, Physical Review E 86, 056102 (2012).
  * Zhang and Motter [2020] Y. Zhang and A. E. Motter, SIAM Review 62, 817 (2020).
  * Golub and Kahan [1965] G. Golub and W. Kahan, Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis 2, 205 (1965).
  * Liao et al. [2007] X. Liao, L. Wang, and P. Yu, _Stability of dynamical systems_ (Elsevier, 2007).
  * Lorenz [1963] E. N. Lorenz, Journal of the atmospheric sciences 20, 130 (1963).