Fully solvable lower dimensional dynamics of Cartesian product of Kuramoto models

A broad range of example systems shows synchronization properties as the interaction between units change, e.g. birds flocking, male fireflies flashing together, heart beating of mother and fetal, neurons in the brain, cardiac and respiratory system [1, 2]. One of traditional models to simulate synchronization dynamics is the Kuramoto phase oscillator model, which still attracts much attention in the literature [3–6].

In the traditional Kuramoto model, as the coupling strength increases, a transition from an incoherent to a coherent state takes place generically at a critical value of λ_c, after which the interacting units follow the same dynamical behavior. This macroscopic appearance of synchronization is often described by an order parameter r which is normalized to r ∈ [0, 1]. Namely, a small value of r ≈ 0 corresponds to an incoherent state, while a large value of r close to 1 indicates a high degree of synchrony. Both the natural frequencies of oscillators and the coupling strength determine the synchronization transition properties. As first pointed in [7], in most cases the transition is second-order like, with the order parameter r growing continuously when the coupling strength passes the critical threshold λ_c of synchronization [8]. On the other hand, abrupt discontinuous transitions to synchronization have been recently reported in both all-to-all coupling [3, 9] and complex network topologies [10–15] where an infinitesimal variation of the coupling strength gives rise to a macroscopic explosive transition to synchronization. In the study of explosive sync, the transition from an initially incoherent to a coherent state is referred to as forward continuation curve when the coupling strength is progressively increased and the corresponding critical coupling is termed as ${\lambda }_{c}^{f}$. In contrast, the desynchronization transition from an initially coherent to an incoherent state is called backward continuation curve when the coupling is decreased and the critical coupling is denoted as ${\lambda }_{c}^{b}$. In addition, a clear hysteresis area has been observed between these synchronization and desynchronization processes, which is further denoted as $S={\lambda }_{c}^{f}-{\lambda }_{c}^{b}$.

Recently, star graphs have been applied to understand the global explosive behavior of the order parameter r which shed insights for other cases of more general heterogeneous network settings [10, 13, 16]. The two fundamental results of explosive synchronization (discontinuity and hysteresis) have been delineated by the Watanabe–Strogatz (WS) approach [17, 18]. More specifically, the exact nonlinear equation for the order parameter r of the high dimensional coupled system has been explicitly obtained and the different synchronized states correspond to different steady states of the equation. Furthermore, different stability conditions of coexisting fixed points in the parameter space lead to the hysteresis behavior and the discontinuous transitions in both the forward and backward continuation curves [16].

From the viewpoint of network topologies, star graphs are considered as building motifs to generate a big graph by several graph operations, e.g. Cartesian product, direct tensor product and strong product [19, 20]. For example, the Cartesian product of graphs is a commutative, associative binary operation on graphs [21]. It has many useful properties, most of which can be derived from the factors. Furthermore, several multilayer network properties have been obtained by graph product operations [20]. On the other hand, from the viewpoint of dynamics on top of networks, it remains largely undisclosed that the synchronization process is obtained by similar graph product operations except some discussions on eigenspectra [22, 23]. In our work [24], we have provided a novel framework to obtain a canonical Kuramoto model by the Cartesian product operation from two independent factor graphs. In this earlier work, we focused on the Cartesian product for two basic network graphs of star and ring where we found a mixture state of both an explosive transition to sync in the forward curve and a continuous desynchronization transition in the backward curve. This mixture state of synchronization transitions cannot be easily observed in a single factor graph. However, the lower dimensional equations for the order parameters of the Cartesian product model have not been discussed in the literature.

In this work, we provide a more general dimension reduction treatment of the WS ansatz to the Cartesian product model, obtaining fully solved lower dimensional dynamical solutions of multiple fixed points for the order parameters. The stability of each fixed point has been obtained by a linear analysis. Comparing to the case of a single star graph, the results are richer depending on the interaction between the hysteresis areas of the two independent factor stars. In addition, cluster synchronization solutions have been obtained for the product model. A cluster synchronized state represents that the network evolves into subsets of oscillators in which members of the same cluster are synchronized, but members of different clusters are not [25, 26]. Together with Chimera states [27–29], cluster synchronization is one of most interesting partial synchronization scenarios that has attracted both theoretical and experimental studies [30, 31]. Recently, a computational group theory has been developed to characterize the emergence and stability conditions of cluster synchronization [26]. More specifically, one has to identify the set of symmetries of the network of interest by discrete algebra routines [25] or by approximation techniques when there are system parameter mismatches [32, 33]. Then the nodes of the network are partitioned into M clusters, which yields disjoint sets of nodes when all of the symmetry operations are applied to permute one from the other. Importantly, the dynamics of oscillators in each disjoint set is essentially unchanged by the permutations, forming a cluster of synchronized oscillators. Once the clusters are identified, the stability of the clusters can be further analyzed by the corresponding variational equations of the system. Differently from the literature, we emphasize that the proposed cluster synchronized behavior in this work are synthetic states that are constructed from the dynamics of two independent subgraphs by means of Cartesian graph operation, which therefore provides a novel 'bottom-up' framework to generate more complicate dynamics encountered in complex systems.

The outline of this paper is as follows: in section 2, we introduce the Cartesian product model from two independent star graphs and provide the WS ansatz. In section 3, we show the steady state solutions of the ensemble order parameters and their respective stability conditions. Numerical simulation are presented in section 4. Finally, our main conclusions are summarized in section 5.

We consider the dynamics of two independent star networks of N_1,2 leaf nodes and a central hub. The degree of a node is the number of connections it receives. So, the degree of the leaf nodes equals one and the degree of the central hub equals N. The equations of motion are described by

$$\begin{eqnarray}{G}_{1}:\left\{\begin{array}{rcl}{\dot{{\theta }_{i}}}^{\left(1\right)} & = & {\omega }_{1}+\lambda \sin ({\theta }_{h}^{\left(1\right)}-{\theta }_{i}^{\left(1\right)}),\,\ \mathrm{for}\ \ i=1,2,\ \cdots ,\ {N}_{1},\\ {\dot{\theta }}_{h}^{\left(1\right)} & = & {\beta }_{1}{\omega }_{1}+{\beta }_{1}\displaystyle \frac{\lambda }{{N}_{1}}\displaystyle \sum _{i=1}^{{N}_{1}}\sin ({\theta }_{i}^{\left(1\right)}-{\theta }_{h}^{\left(1\right)}),\end{array}\right.\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}{G}_{2}:\left\{\begin{array}{rcl}{\dot{{\theta }_{k}}}^{\left(2\right)} & = & {\omega }_{2}+\lambda \sin ({\theta }_{h}^{\left(2\right)}-{\theta }_{k}^{\left(2\right)}),\,\ \mathrm{for}\ \ k=1,2,\ \cdots ,\ {N}_{2},\\ {\dot{\theta }}_{h}^{\left(2\right)} & = & {\beta }_{2}{\omega }_{2}+{\beta }_{2}\displaystyle \frac{\lambda }{{N}_{2}}\displaystyle \sum _{k=1}^{{N}_{2}}\sin ({\theta }_{k}^{\left(2\right)}-{\theta }_{h}^{\left(2\right)}),\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where ω_1,2 is the natural frequency of the leaves, λ is the overall coupling strength, and β_1,2 is a parameter controlling the frequency mismatch between the hub and the leaves [16]. In this model, β_1,2 > 1 mimics a positive correlation between the hub's natural frequency and its degree, i.e. the hub of larger degree has a larger frequency than that of a leaf node [10, 13]. The parameter β_1,2 helps to understand a more general effect besides the network degrees. In addition, we consider a lower dimensional dynamics in the thermodynamic limit ${N}_{\mathrm{1,2}}\to \infty$, so the normalization is necessary to make sense of the limit, otherwise the hub would rotate infinitely fast.

Following [16], we introduce the phase differences as

$$\begin{eqnarray}&&{\varphi }_{i}^{\left(1\right)}={\theta }_{i}^{\left(1\right)}-{\theta }_{h}^{\left(1\right)},{\varphi }_{k}^{\left(2\right)}={\theta }_{k}^{\left(2\right)}-{\theta }_{h}^{\left(2\right)}.\end{eqnarray} \tag{ 3 }$$

In addition, the order parameters of the two independent networks are defined as

$$\begin{eqnarray}&&{z}_{1}(t)={r}_{1}(t){{\rm{e}}}^{{\bf{i}}{{\rm{\Phi }}}_{1}(t)}=\frac{1}{{N}_{1}}\displaystyle \sum _{i=1}^{{N}_{1}}{{\rm{e}}}^{{\bf{i}}{\varphi }_{i}^{(1)}},{z}_{2}(t)={r}_{2}(t){{\rm{e}}}^{{\bf{i}}{{\rm{\Phi }}}_{2}(t)}=\frac{1}{{N}_{2}}\displaystyle \sum _{k=1}^{{N}_{2}}{{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}},\left[\right]\end{eqnarray} \tag{ 4 }$$

where is the bold font i is for the imaginary unit throughout this paper. Then the original models (equations (1), (2)) are rewritten as the following compact forms

$$\begin{eqnarray}&&{G}_{1}:{\dot{{\varphi }_{i}}}^{(1)}=(1-{\beta }_{1}){\omega }_{1}-{\beta }_{1}\lambda \mathrm{Im}({z}_{1}(t))-\lambda \sin ({\varphi }_{i}^{(1)}),\,{\rm{for}}\,i=1,2,\ \cdots ,\ {N}_{1}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{G}_{2}:{\dot{{\varphi }_{k}}}^{(2)}=(1-{\beta }_{2}){\omega }_{2}-{\beta }_{2}\lambda \mathrm{Im}({z}_{2}(t))-\lambda \sin ({\varphi }_{k}^{(2)}),\,{\rm{for}}\,k=1,2,\ \cdots ,\ {N}_{2}.\end{eqnarray} \tag{ 6 }$$

For individual factor graphs, the order parameters of equations (4) represent the mean-field coupling terms in the above compact equations (5), (6).

The Cartesian product of the two stars G₁ and G₂ are schematically shown in figure 1. Meanwhile, we use the notation $\left\langle {ik}\right\rangle$ to represent the index of node on the product ${G}_{1}\square {G}_{2}$. In addition, on the Cartesian product graph, the phase of the node $\left\langle {ik}\right\rangle$ is defined as ${\varphi }_{\left\langle {ik}\right\rangle }={\varphi }_{i}^{\left(1\right)}+{\varphi }_{k}^{\left(2\right)}$. Note that the definition of the phase of the node on the product graph as the summation of the respective phases on the factor subgraphs yields the canonical equations of the Kuramoto model on ${G}_{1}\square {G}_{2}$ [24]. Furthermore, with the commutative and associative properties of the phase summation of the Cartesian product operation, we easily generalize the present results to the case of n factor subgraphs [24].

Zoom In Zoom Out Reset image size

In addition, the order parameter of the Cartesian product model is defined as

$$\begin{eqnarray}&&Z(t)=r(t){{\rm{e}}}^{{\bf{i}}{\rm{\Phi }}(t)}=\displaystyle \frac{1}{{N}_{1}{N}_{2}}\displaystyle \sum _{{i}=1}^{{N}_{1}}\displaystyle \sum _{k=1}^{{N}_{2}}{{\rm{e}}}^{{\bf{i}}{\varphi }_{\langle {ik}\rangle }}.\end{eqnarray} \tag{ 7 }$$

In a full analogy, we note that the above definition of Z(t) preserves the mean-field coupling properties. In the present framework of phase summation, the right-hand side of equation (7) can be further expanded, which yields Z(t) = z₁(t) z₂(t). Namely, the order parameter on ${G}_{1}\square {G}_{2}$ is the product of two factor graphs. The summation rule of phases can be generalized to the case of more than two factor graphs straightforwardly, while preserving that the order parameter Z is the product of subgraphs. For instance, given three factors G₁, G₂ and G₃, the order parameter R of the product ${G}_{1}\square {G}_{2}\square {G}_{3}$ is R = R₁R₂R₃ [24].

With the above phase definition, the time derivatives of the phases ${\varphi }_{\left\langle {ik}\right\rangle }$ on the Cartesian product ${G}_{1}\square {G}_{2}$ are

$$\begin{eqnarray}&&\,\displaystyle \frac{{\rm{d}}{\varphi }_{\left\langle {ik}\right\rangle }}{{\rm{d}}{t}}=\displaystyle \frac{{\rm{d}}{\varphi }_{i}^{\left(1\right)}}{{\rm{d}}{t}}+\displaystyle \frac{{\rm{d}}{\varphi }_{k}^{\left(2\right)}}{{\rm{d}}{t}}\,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\,=\ (1-{\beta }_{1}){\omega }_{1}+(1-{\beta }_{2}){\omega }_{2}-{\beta }_{1}\lambda \mathrm{Im}({z}_{1}(t))-{\beta }_{2}\lambda \mathrm{Im}({z}_{2}(t))-\lambda \sin ({\varphi }_{i}^{\left(1\right)})-\lambda \sin ({\varphi }_{k}^{\left(2\right)})\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\,=\ (1-{\beta }_{1}){\omega }_{1}+(1-{\beta }_{2}){\omega }_{2}-{\beta }_{1}\lambda {r}_{1}\sin {{\rm{\Phi }}}_{1}-{\beta }_{2}\lambda {r}_{2}\sin {{\rm{\Phi }}}_{2}+{\bf{i}}\displaystyle \frac{\lambda }{2}({{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}}-{{\rm{e}}}^{-{\bf{i}}{\varphi }_{{i}}^{(2)}})+{\bf{i}}\displaystyle \frac{\lambda }{2}({{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}}-{{\rm{e}}}^{-{\bf{i}}{\varphi }_{k}^{(2)}})\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\,=\ {g}_{1}+{g}_{2}\,+\,{{f}{\rm{e}}}^{{\bf{i}}{\varphi }_{{\bf{i}}}^{\left(1\right)}}+{f}^{* }{{\rm{e}}}^{-{\bf{i}}{\varphi }_{{\bf{i}}}^{\left(1\right)}}+{{f}{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{\left(2\right)}}+{f}^{* }{{\rm{e}}}^{-{\bf{i}}{\varphi }_{k}^{\left(2\right)}},\,\end{eqnarray} \tag{ 11 }$$

where f^* are the complex conjugates of f. The f and g_1,2 are further defined as

$$\begin{eqnarray}&&f={\bf{i}}\displaystyle \frac{\lambda }{2},{g}_{1}={\rm{\Delta }}{\omega }_{1}-{\beta }_{1}\lambda {r}_{1}\sin {{\rm{\Phi }}}_{1},{g}_{2}={\rm{\Delta }}{\omega }_{2}-{\beta }_{2}\lambda {r}_{2}\sin {{\rm{\Phi }}}_{2},{\rm{\Delta }}{\omega }_{1}=(1-{\beta }_{1}){\omega }_{1},{\rm{\Delta }}{\omega }_{2}=(1-{\beta }_{2}){\omega }_{2}.\end{eqnarray} \tag{ 12 }$$

Note that, in equation (9), the term (1 − β₁)ω₁ + (1 − β₂)ω₂ represents the natural frequency of the oscillator on the product ${G}_{1}\square {G}_{2}$, while the rest terms are the coupling [24].

The above phase equation (11) has exactly the form such that the WS ansatz can be applied. The WS approach [34, 35] is applicable for systems of identical oscillators driven by a common force. More specifically, in the product model of equation (11), identical oscillators g₁ + g₂ are driven by the arbitrary complex ${{f}{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}}+{{f}{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}}+{\rm{c}}.{\rm{c}}.$, where c.c. are the complex conjugates.

Next, the basic idea is to expand the model system (equation (11)) in terms of the global variables of z₁, z₂ and Z. Therefore, we first consider the relationship

$$\begin{eqnarray}&&\,\displaystyle \frac{\partial }{\partial {t}}({{\rm{e}}}^{{\bf{i}}{\varphi }_{\langle {ik}\rangle }})={\bf{i}}{{\rm{e}}}^{{\bf{i}}({\varphi }_{{i}}^{(1)}+{\varphi }_{k}^{(2)})}{\dot{\varphi }}_{\langle {ik}\rangle }\,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\,=\ {\bf{i}}{{f}{\rm{e}}}^{{\bf{i}}(2{\varphi }_{{i}}^{(1)}+{\varphi }_{k}^{(2)})}+{\bf{i}}{{f}{\rm{e}}}^{{\bf{i}}({\varphi }_{{i}}^{(1)}+2{\varphi }_{k}^{(2)})}+{\bf{i}}({g}_{1}+{g}_{2}){{\rm{e}}}^{{\bf{i}}({\varphi }_{{i}}^{(1)}+{\varphi }_{k}^{(2)})}+{\bf{i}}{f}^{* }{{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}}+{\bf{i}}{f}^{* }{{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}}.\end{eqnarray} \tag{ 14 }$$

Furthermore, in the formulation presented in [17, 18], we need to introduce a series of Möbius transformations that expand the exponential functions in the above polynomial in terms of the order parameters. Importantly, we introduce

$$\begin{eqnarray}&&{{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}(t)}=\displaystyle \frac{{z}_{1}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}{1+{z}_{1}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}.\end{eqnarray} \tag{ 15 }$$

Similar Möbius transformations of other terms have been included in the appendix. Note that the variables of ξ_i,k are constants. Remarkably, in the case of a uniform distribution of ξ_i,k, the global variables α_1,2 and ξ_i,k do not enter the equation for z_1,2 and Z [16, 17]. Thus, the equation for Z is a closed equation that fully describes the dynamics of the system (equation (11)), and therefore in the following, we focus on the discussion on Z only. Furthermore, it has been shown in [36, 37], for ${N}_{\mathrm{1,2}}\to \infty$ and the uniform distribution of constants of motion ξ_i,k, the WS variables z₁, z₂ coincide with the local Kuramoto mean-field, namely yielding the Ott–Antonsen (OA) ansatz [38–40].

Putting these Möbius transformations into the two sides of equation (14), we obtain the following relationship by equalizing the non-exponential terms of the two sides of the equations (13), (14).

$$\begin{eqnarray}&&{\dot{z}}_{1}{z}_{2}+{\dot{z}}_{2}{z}_{1}={\bf{i}}f({z}_{1}^{2}{z}_{2}+{z}_{1}{z}_{2}^{2})+{\bf{i}}({g}_{1}+{g}_{2}){z}_{1}{z}_{2}+{\bf{i}}{f}^{* }({z}_{1}+{z}_{2}).\end{eqnarray} \tag{ 16 }$$

As it has been proved in [24], the order parameter Z of the Cartesian model is the product of the order parameters of the two independent factors, namely, Z = z₁z₂. Then, we have $\dot{Z}={\dot{z}}_{1}{z}_{2}+{z}_{1}{\dot{z}}_{2}$. Therefore, in terms of the global variables Z, the phase equation equation (11) is expressed as

$$\begin{eqnarray}&&\dot{Z}={\bf{i}}({z}_{1}+{z}_{2})({fZ}+{f}^{* })+{\bf{i}}({g}_{1}+{g}_{2})Z.\end{eqnarray} \tag{ 17 }$$

Considering the definition of the order parameter $Z={{r}{\rm{e}}}^{{\bf{i}}{\rm{\Phi }}}$ (equation (7)), we obtain the following two nonlinear coupled equations for the global order parameter in the complex plane

$$\begin{eqnarray}&&\,\dot{r}=f(r,{\rm{\Phi }})=\displaystyle \frac{\lambda }{2}{r}_{1}(1-{r}_{2}^{2})\cos {{\rm{\Phi }}}_{2}+\displaystyle \frac{\lambda }{2}{r}_{2}(1-{r}_{1}^{2})\cos {{\rm{\Phi }}}_{1},\,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{rcl}\dot{{\rm{\Phi }}} & = & g(r,{\rm{\Phi }})=(1-{\beta }_{1}){w}_{1}+(1-{\beta }_{2}){w}_{2}-{\beta }_{1}\lambda {r}_{1}\sin {{\rm{\Phi }}}_{1}-{\beta }_{2}\lambda {r}_{2}\sin {{\rm{\Phi }}}_{2}\\ & & \,-\ \displaystyle \frac{\lambda (1+{r}_{1}^{2})}{2{r}_{1}}\sin {{\rm{\Phi }}}_{1}-\displaystyle \frac{\lambda (1+{r}_{2}^{2})}{2{r}_{2}}\sin {{\rm{\Phi }}}_{2}.\end{array}\end{eqnarray} \tag{ 19 }$$

There are three steady states of the nonlinear coupled system (equations (18), (19)), namely,

• 1.  
  Full synchronization of r₁ = r₂ = 1,
• 2.  
  Non-synchronization of ${{\rm{\Phi }}}_{1}=\pm \tfrac{\pi }{2}$ and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$,
• 3.  
  Cluster synchronized state of r₁ = 1 and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$ $\left(\mathrm{or}\,\mathrm{by}\,\mathrm{symmetry}\,{r}_{2}\,=\,1\,\mathrm{and}{{\rm{\Phi }}}_{1}=\pm \tfrac{\pi }{2}\right)$.

In the following subsections, we get the explicit expressions for the global order parameter Z of the Cartesian product model in terms of Z = z₁z₂. To this end, we first denote the forward critical coupling threshold values for the two independent factors following the notations of [16]:

$$\begin{eqnarray}&&{\lambda }_{c}^{{f}_{1}}=\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}}{\sqrt{2{\beta }_{1}+1}},\ {\lambda }_{c}^{{f}_{2}}=\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}}{\sqrt{2{\beta }_{2}+1}},\end{eqnarray} \tag{ 20 }$$

and two backward critical coupling thresholds

$$\begin{eqnarray}&&{\lambda }_{c}^{{b}_{1}}=\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}}{{\beta }_{1}+1},\ {\lambda }_{c}^{{b}_{2}}=\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}}{{\beta }_{2}+1}.\end{eqnarray} \tag{ 21 }$$

In addition, we consider β_1,2 > 1 which implement positive correlations between the node frequency and its associated number of connections [16]. In other words, the following inequalities do hold always

$$\begin{eqnarray}&&{\lambda }_{c}^{{f}_{1}}\gt {\lambda }_{c}^{{b}_{1}},\,{\lambda }_{c}^{{f}_{2}}\gt {\lambda }_{c}^{{b}_{2}},\end{eqnarray} \tag{ 22 }$$

which further ensure that the hysteretic areas exist for both independent star graphs, namely, ${S}_{{G}_{1}}={\lambda }_{c}^{{f}_{1}}-{\lambda }_{c}^{{b}_{1}}$ and ${S}_{{G}_{2}}={\lambda }_{c}^{{f}_{2}}-{\lambda }_{c}^{{b}_{2}}$.

3.1. Full synchrony for the product model r₁ = r₂ = 1

In this case, equations (18), (19) are simplified as

$$\begin{eqnarray}&&(1-{\beta }_{1}){\omega }_{1}-(1+{\beta }_{1})\lambda \sin {{\rm{\Phi }}}_{1}\,=\,0,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&(1-{\beta }_{2}){\omega }_{2}-(1+{\beta }_{2})\lambda \sin {{\rm{\Phi }}}_{2}\,=\,0.\end{eqnarray} \tag{ 24 }$$

Since β_1,2, ω_1,2 and λ are positive, the solutions Φ_1,2 exist only if $\left|\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right|\leqslant 1$ and $\left|\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right|\leqslant 1$. Note that this corresponds to the existence of the phase locking manifold [13], which therefore yields the critical coupling for the backward continuation curves as ${\lambda }_{c}^{{b}_{1}}$ and ${\lambda }_{c}^{{b}_{2}}$, respectively. In this case, the solutions of z₁ and z₂ are easily obtained as

$$\begin{eqnarray}&&{z}_{1}={{\rm{e}}}^{{\bf{i}}(\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right])},\mathrm{for}\ \lambda \gt {\lambda }_{c}^{{b}_{1}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{z}_{2}={{\rm{e}}}^{{\bf{i}}(\arcsin \left[\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right])},\mathrm{for}\ \lambda \gt {\lambda }_{c}^{{b}_{2}}.\end{eqnarray} \tag{ 26 }$$

Furthermore, keeping in mind the condition of ${\rm{\Phi }}={{\rm{\Phi }}}_{1}+{{\rm{\Phi }}}_{2}$, and therefore the solution of the product model has the following expression

$$\begin{eqnarray}&&Z={z}_{1}{z}_{2}={{\rm{e}}}^{{\bf{i}}(\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right]+\arcsin \left[\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right])},\mathrm{for}\ \lambda \gt \max \{{\lambda }_{c}^{{b}_{1}},{\lambda }_{c}^{{b}_{2}}\}.\end{eqnarray} \tag{ 27 }$$

3.2. Non-synchrony for the product model ${{\rm{\Phi }}}_{1}=\pm \tfrac{\pi }{2}$ and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$

Inserting the steady states of ${{\rm{\Phi }}}_{1}=\pm \tfrac{\pi }{2}$ and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$ into equations (18), (19), we obtain

$$\begin{eqnarray}&&(1-{\beta }_{1}){\omega }_{1}\mp {\beta }_{1}\lambda {r}_{1}\mp \displaystyle \frac{\lambda (1+{r}_{1}^{2})}{2{r}_{1}}=0,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&(1-{\beta }_{2}){\omega }_{2}\mp {\beta }_{2}\lambda {r}_{2}\mp \displaystyle \frac{\lambda (1+{r}_{2}^{2})}{2{r}_{2}}=0.\end{eqnarray} \tag{ 29 }$$

Considering the fact that ${z}_{\mathrm{1,2}}={r}_{\mathrm{1,2}}{{\rm{e}}}^{{{\rm{\Phi }}}_{\mathrm{1,2}}}=\pm {\bf{i}}{r}_{\mathrm{1,2}}$ or ${r}_{\mathrm{1,2}}=\mp {\bf{i}}{z}_{\mathrm{1,2}}$, we get

$$\begin{eqnarray}&&{z}_{1}=-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}\pm \sqrt{{\left[({\beta }_{1}-1){\omega }_{1}\right]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda },\mathrm{for}\ ({\beta }_{1}-1){\omega }_{1}\geqslant \sqrt{2{\beta }_{1}\,+\,1}\lambda ,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{z}_{2}=-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}\pm \sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda },\mathrm{for}\ ({\beta }_{2}-1){\omega }_{2}\geqslant \sqrt{2{\beta }_{2}\,+\,1}\lambda .\end{eqnarray} \tag{ 31 }$$

Note that both z_1,2 exist only if $({\beta }_{\mathrm{1,2}}-1){\omega }_{\mathrm{1,2}}\geqslant \sqrt{2{\beta }_{\mathrm{1,2}}+1}\lambda$ hold. Furthermore, both '+' and '−' are possible if (β_1,2 − 1)ω_1,2 ≤ (β_1,2 + 1)λ, while '−' is possible if (β_1,2 − 1)ω_1,2 > (β_1,2 + 1)λ. In consequence, the expression of the global order parameter of the product model Z = z₁z₂ has several different cases while considering ${r}_{\mathrm{1,2}}=| {z}_{\mathrm{1,2}}| \leqslant 1$. Depending on the relationship between ${\lambda }_{c}^{{f}_{1}}\lt {\lambda }_{c}^{{f}_{2}}$ or ${\lambda }_{c}^{{f}_{1}}\gt {\lambda }_{c}^{{f}_{2}}$, we will have different solutions for Z and the details are omitted here for the simplification purpose, which will be summarized in section 3.4.

3.3. Cluster synchronized states of r₁ = 1 and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$

We focus on the discussion on the case of r₁ = 1 and ${{\rm{\Phi }}}_{2}=\pm \tfrac{\pi }{2}$, while the symmetric case of r₂ = 1 and ${{\rm{\Phi }}}_{1}=\pm \tfrac{\pi }{2}$ shows a full analogy. In this case, we have the following two equations in order to obtain the steady states, namely

$$\begin{eqnarray}&&(1-{\beta }_{1}){\omega }_{1}-(1+{\beta }_{1})\lambda \sin {{\rm{\Phi }}}_{1}=0,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&(1-{\beta }_{2}){\omega }_{2}\mp {\beta }_{2}\lambda {r}_{2}\mp \displaystyle \frac{\lambda (1+{r}_{2}^{2})}{2{r}_{2}}=0.\end{eqnarray} \tag{ 33 }$$

From equation (32), we get

$$\begin{eqnarray}&&{z}_{1}={{\rm{e}}}^{{\bf{i}}(\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right])},\mathrm{only}\ \mathrm{if}\ ({\beta }_{1}-1){\omega }_{1}\leqslant (1+{\beta }_{1})\lambda .\end{eqnarray} \tag{ 34 }$$

From equation (33), we obtain

$$\begin{eqnarray}&&{z}_{2}=-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}\pm \sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda },\mathrm{for}\ ({\beta }_{2}-1){\omega }_{2}\geqslant \sqrt{2{\beta }_{2}+1}\lambda .\end{eqnarray} \tag{ 35 }$$

Both '+' and '−' are possible for z₂ if (β₂ − 1)ω₂ ≤ (β₂ + 1)λ, while only '−' is possible for z₂ if (β₂ − 1)ω₂ > (β₂ + 1)λ. The order parameter Z for the product model is Z = z₁z₂ again with different expressions depending on β_1,2 and ω_1,2, which will be summarized in section 3.4.

3.4. Intermediate summary of steady states

Note that the critical coupling thresholds ${\lambda }_{c}^{{b}_{1}},{\lambda }_{c}^{{f}_{1}},{\lambda }_{c}^{{b}_{2}}$ and ${\lambda }_{c}^{{f}_{2}}$ subdivide the coupling strength into five subintervals. All solutions for Z of the product model have been summarized in the following when ${\lambda }_{c}^{{b}_{2}}\lt {\lambda }_{c}^{{b}_{1}}$. When switching the notations of the two factor subgraphs G₁ and G₂, we obtain the case of ${\lambda }_{c}^{{b}_{1}}\lt {\lambda }_{c}^{{b}_{2}}$ which is in a full analogy as we summarized below. Therefore, we only focus on ${\lambda }_{c}^{{b}_{2}}\lt {\lambda }_{c}^{{b}_{1}}\lt {\lambda }_{c}^{{f}_{1}}$. In this case, we further have three subcategories as illustrated in figure 2, depending on the critical value of ${\lambda }_{c}^{{f}_{2}}$. From the viewpoint of hysteresis areas of ${S}_{{G}_{1}}$ and ${S}_{{G}_{2}}$, these three subcategories correspond respectively to ${S}_{{G}_{1}}\cap {S}_{{G}_{2}}=\varnothing$, ${S}_{{G}_{1}}\cap {S}_{{G}_{2}}\ne \varnothing$, and ${S}_{{G}_{1}}\subset {S}_{{G}_{2}}$.

Zoom In Zoom Out Reset image size

All possible solutions are expressed in terms of the product of the two subgraphs, namely, Z = z₁z₂. The respective existence regimes of Z are illustrated in figures 2(a)–(c).

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{{a}_{1}} & = & \left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}-\sqrt{{\left[({\beta }_{1}-1){\omega }_{1}\right]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\\ & & \cdot \ \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}-\sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{array}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{{a}_{2}} & = & \left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}-\sqrt{{\left[({\beta }_{1}-1){\omega }_{1}\right]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\\ & & \cdot \ \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}+\sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{array}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{{a}_{3}} & = & \left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}+\sqrt{{\left[({\beta }_{1}-1){\omega }_{1}\right]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\\ & & \cdot \ \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}-\sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{array}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{{a}_{4}} & = & \left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}+\sqrt{{\left[({\beta }_{1}-1){\omega }_{1}\right]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\\ & & \cdot \ \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}+\sqrt{{\left[({\beta }_{2}-1){\omega }_{2}\right]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{Z}_{{a}_{5}}=\left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}-\sqrt{{[({\beta }_{1}-1){\omega }_{1}]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\cdot {{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right]}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{Z}_{{a}_{6}}=\left[-{\bf{i}}\displaystyle \frac{({\beta }_{1}-1){\omega }_{1}+\sqrt{{[({\beta }_{1}-1){\omega }_{1}]}^{2}-(2{\beta }_{1}+1){\lambda }^{2}}}{(2{\beta }_{1}+1)\lambda }\right]\cdot {{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right]}\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{Z}_{{a}_{7}}={{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right]}\cdot \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}+\sqrt{{[({\beta }_{2}-1){\omega }_{2}]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{Z}_{{a}_{8}}={{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right]}\cdot \left[-{\bf{i}}\displaystyle \frac{({\beta }_{2}-1){\omega }_{2}-\sqrt{{[({\beta }_{2}-1){\omega }_{2}]}^{2}-(2{\beta }_{2}+1){\lambda }^{2}}}{(2{\beta }_{2}+1)\lambda }\right]\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{Z}_{{s}_{\mathrm{1,2},\mathrm{3,4}}}={{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{1}){\omega }_{1}}{(1+{\beta }_{1})\lambda }\right]}\cdot {{\rm{e}}}^{{\bf{i}}\arcsin \left[\tfrac{(1-{\beta }_{2}){\omega }_{2}}{(1+{\beta }_{2})\lambda }\right]}\end{eqnarray} \tag{ 44 }$$

Note that in the case of figure 2(a), the solutions of ${Z}_{{a}_{1}}$, ${Z}_{{a}_{2}}$, ${Z}_{{a}_{3}}$, ${Z}_{{a}_{4}}$, ${Z}_{{a}_{7}}$ and ${Z}_{{a}_{8}}$ do not exist. Additionally, the solutions ${Z}_{{s}_{\mathrm{1,2},\mathrm{3,4}}}$ correspond to four synchronous states, showing different stability conditions as we summarize in the next section.

3.5. Linear stability analysis

Accordingly, we obtain the stability of each fixed point (equations (36)–(43)). Therefore, we first obtain the Jacobian matrix of the system (equations (18), (19)) which is expressed as

$$\begin{eqnarray}{\bf{J}}=\left(\begin{array}{cc}\displaystyle \frac{\partial f}{\partial r} & \displaystyle \frac{\partial f}{\partial {\rm{\Phi }}}\\ \displaystyle \frac{\partial g}{\partial r} & \displaystyle \frac{\partial g}{\partial {\rm{\Phi }}}\end{array}\right),\end{eqnarray} \tag{ 45 }$$

where f and g are the right-hand side of the system (equations (18), (19)). The elements of J are the following:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial r}=\displaystyle \frac{\left[\tfrac{\lambda }{2}(1-{r}_{1}^{2})\cos {{\rm{\Phi }}}_{1}-\lambda {r}_{1}{r}_{2}\cos {{\rm{\Phi }}}_{2}\right]\left[\tfrac{\lambda }{2}(1-{r}_{2}^{2})\cos {{\rm{\Phi }}}_{2}-\lambda {r}_{1}{r}_{2}\cos {{\rm{\Phi }}}_{1}\right]}{\tfrac{\lambda }{2}{r}_{1}(1-{r}_{2}^{2})\cos {{\rm{\Phi }}}_{2}-\lambda {r}_{1}^{2}{r}_{2}\cos {{\rm{\Phi }}}_{1}+\tfrac{\lambda }{2}{r}_{2}(1-{r}_{1}^{2})\cos {{\rm{\Phi }}}_{1}-\lambda {r}_{1}{r}_{2}^{2}\cos {{\rm{\Phi }}}_{2}},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial {\rm{\Phi }}}=-\displaystyle \frac{[\lambda {r}_{2}(1-{r}_{1}^{2})\sin {{\rm{\Phi }}}_{1}][\lambda {r}_{1}(1-{r}_{2}^{2})\sin {{\rm{\Phi }}}_{2}]}{2\lambda {r}_{1}(1-{r}_{2}^{2})\sin {{\rm{\Phi }}}_{2}\,+\,2\lambda {r}_{2}(1-{r}_{1}^{2})\sin {{\rm{\Phi }}}_{1}},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial g}{\partial r}=\displaystyle \frac{\left[{\beta }_{2}\lambda \sin {{\rm{\Phi }}}_{2}+\tfrac{({r}_{2}^{2}-1)\lambda \sin {{\rm{\Phi }}}_{2}}{2{r}_{2}^{2}}\right]\left[{\beta }_{1}\lambda \sin {{\rm{\Phi }}}_{1}+\tfrac{({r}_{1}^{2}-1)\lambda \sin {{\rm{\Phi }}}_{1}}{2{r}_{1}^{2}}\right]}{-{\beta }_{1}\lambda {r}_{1}\sin {{\rm{\Phi }}}_{1}-\tfrac{({r}_{1}^{2}-1)\lambda \sin {{\rm{\Phi }}}_{1}}{2{r}_{1}}-{\beta }_{2}\lambda {r}_{2}\sin {{\rm{\Phi }}}_{2}-\tfrac{({r}_{2}^{2}-1)\lambda \sin {{\rm{\Phi }}}_{2}}{2{r}_{2}}},\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial g}{\partial {\rm{\Phi }}}=\displaystyle \frac{\left[-{\beta }_{1}\lambda {r}_{1}\cos {{\rm{\Phi }}}_{1}-\tfrac{\lambda (1+{r}_{1}^{2})}{2{r}_{1}}\cos {{\rm{\Phi }}}_{1}\right]\left[-{\beta }_{2}\lambda {r}_{2}\cos {{\rm{\Phi }}}_{2}-\tfrac{\lambda (1+{r}_{2}^{2})}{2{r}_{2}}\cos {{\rm{\Phi }}}_{2}\right]}{-{\beta }_{2}\lambda {r}_{2}\cos {{\rm{\Phi }}}_{2}-\tfrac{\lambda (1+{r}_{2}^{2})}{2{r}_{2}}\cos {{\rm{\Phi }}}_{2}-{\beta }_{1}\lambda {r}_{1}\cos {{\rm{\Phi }}}_{1}-\tfrac{\lambda (1+{r}_{1}^{2})}{2{r}_{1}}\cos {{\rm{\Phi }}}_{1}}.\end{eqnarray} \tag{ 49 }$$

The stability is studied by inserting each steady state solution into the trace and determinant of the Jacobian matrix (equation (45)). Due to the lengthy of the derivation, we only summarize the stability of the fixed points and the corresponding physical meaning according to the three subcategories as illustrated in figures 2(a)–(c), which are respectively shown in tables 1–3.

Note that the fixed points of ${{Z}_{s}}_{1,2,3,4}$ correspond to four synchronous states of different stabilities which are determined by the products from ${z}_{{s}_{1}}$ and ${z}_{{s}_{2}}$. Taking figure 2(a) as an example (table 1), ${{Z}_{s}}_{1}(++)$ is a stable sink if both ${z}_{{s}_{1}}$ and ${z}_{{s}_{2}}$ are positive real values, and ${{Z}_{s}}_{2}(--)$ is an unstable source if both ${z}_{{s}_{1}}$ and ${z}_{{s}_{2}}$ are negative real values. On the other hand, ${{Z}_{s}}_{3}(+-)$ is either a stable sink or an unstable source if ${z}_{{s}_{1}}$ is a positive real value and ${z}_{{s}_{2}}$ is a negative real value. Respectively, ${{Z}_{s}}_{4}(-+)$ is either a stable sink or an unstable source if ${z}_{{s}_{1}}$ is a negative real value and ${z}_{{s}_{2}}$ is a positive real value. The annotations of other two cases of tables 2 and 3 have similar stability conditions.

Table 1.  Fixed points of the order parameters with their stability and meaning for figure 2(a).

Fixed point	Stability	Existence region	Physical meaning
${{Z}_{s}}_{1}(++)$	Sink	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{2}(--)$	Source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{3}(+-)$	Sink or source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{4}(-+)$	Sink or source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{a}}_{1}$	Center	$\lambda \lt {\lambda }_{c}^{{f}_{2}}$	Asynchronous state
${{Z}_{a}}_{2}$	Saddle	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{5}$	Center	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Cluster sync
${{Z}_{a}}_{6}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix

Table 2.  Fixed points of the order parameters with their stability and meaning for figure 2(b).

Fixed point	Stability	Existence region	Physical meaning
${{Z}_{s}}_{1}(++)$	Sink	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{2}(--)$	Source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{3}(+-)$	Sink	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{4}(-+)$	Source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{a}}_{1}$	Center	$\lambda \lt {\lambda }_{c}^{{f}_{2}}$	Asynchronous state
${{Z}_{a}}_{2}$	Saddle	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{3}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{4}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{5}$	Center	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Cluster sync
${{Z}_{a}}_{6}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix
${{Z}_{a}}_{7}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{8}$	Center	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Cluster sync

Table 3.  Fixed points of the order parameters with their stability and meaning for figure 2(c).

Fixed point	Stability	Existence region	Physical meaning
${{Z}_{s}}_{1}(++)$	Sink	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{2}(--)$	Source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{3}(+-)$	Sink	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{s}}_{4}(-+)$	Source	$\lambda \gt {\lambda }_{c}^{{b}_{1}}$	Coherent state
${{Z}_{a}}_{1}$	Center	$\lambda \lt {\lambda }_{c}^{{f}_{1}}$	Asynchronous state
${{Z}_{a}}_{2}$	Saddle	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix
${{Z}_{a}}_{3}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix
${{Z}_{a}}_{4}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix
${{Z}_{a}}_{5}$	Center	${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Cluster sync
${{Z}_{a}}_{6}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$	Separatrix
${{Z}_{a}}_{7}$	Saddle	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Separatrix
${{Z}_{a}}_{8}$	Center	${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$	Cluster sync

Now, we follow mainly the simulation routines as presented in [16, 24] and have numerically simulated the model equations (1), (2) by using a fourth order Runge–Kutta integrator with the integration step h = 0.01. When the coupling λ = 0, the initial conditions (ICs) are uniformly distributed over the interval [−π, π]. Then the coupling is increased by a step size Δλ = 0.02 and the ICs for the coupling λ + Δλ are the final states when coupling equals to λ as suggested in [10, 13]. The first T = 10⁵ steps are discarded as transients and the next T iterations are used to estimate the order parameter. We consider the time average of the order parameter $\parallel z\parallel =\tfrac{1}{T}{\sum }_{i=1}^{T}| z(i)|$. Note that this average is useful as the asynchronous fixed points of the order parameter are centers, and so the order parameter in the asynchronous regime shows oscillatory behavior.

There are two equivalent ways to implement the dynamics of the Cartesian product model: (i) We simulate G₁ and G₂ (equations (1), (2)) simultaneously, while phase dynamics of the product G₁$\square$G₂ simply follows the Cartesian product summation rule of the corresponding phases. (ii) We simulate the Cartesian model equation (11), directly. The ODE integrator is performed for N₁ + N₂ phase oscillators in the former case, yielding a better computation efficiency than that of the latter case that is integrated for N₁N₂ oscillators. The additional requirement for the second simulation method is that we have to obtain the connectivity matrix of the product model, in particular, the adjacency matrix of ${G}_{1}\square {G}_{2}$ is the Kronecker sum of the adjacency matrices of G₁ and G₂, namely, ${\boldsymbol{A}}({G}_{1}\square {G}_{2})={\boldsymbol{A}}({G}_{1})\oplus {\boldsymbol{A}}({G}_{2})$. Throughout this work, we have obtained the same results for these two slightly different ways of numerical simulations. In the examples below, we choose N₁ = N₂ = 100.

For a better understanding of the product effects on synchronization transitions, we choose the parameters such that both the forward and backward critical coupling thresholds of G₁ do not vary among the three cases of figures 2(a)–(c), namely, β₁ = 9, ω₁ = 1.2, which yield ${\lambda }_{c}^{{b}_{1}}=0.96$, ${\lambda }_{c}^{{f}_{1}}=2.202$, and ${S}_{{G}_{1}}=1.242$. In the case of G₂, we choose β₂ and ω₂ such that ${\lambda }_{c}^{{b}_{2}}$ is fixed as 0.48, but ${\lambda }_{c}^{{f}_{2}}$ is changed in such a way that, respectively, represents the three different regimes that are illustrated in figures 2(a)–(c). Note that our theoretical predictions obtained by the lower dimensional dynamics of the product model agree very well with the numerical simulations.

4.1. Case A of figure 2(a)

The first case of figure 2(a) is implemented by choosing parameters β₂ = 3 and ω₂ = 0.96, which yield ${\lambda }_{c}^{{b}_{2}}=0.48$, ${\lambda }_{c}^{{f}_{2}}=0.726$, and ${S}_{{G}_{2}}=0.246$, namely ${\lambda }_{c}^{{b}_{2}}\lt {\lambda }_{c}^{{f}_{2}}\lt {\lambda }_{c}^{{b}_{1}}\lt {\lambda }_{c}^{{f}_{1}}$. In this case, the two hysteresis areas ${S}_{{G}_{1}}$ and ${S}_{{G}_{2}}$ are well separated, leading to ${S}_{{G}_{1}}\cap {S}_{{G}_{2}}=\varnothing$ as shown in figures 3(a), (b). Therefore, on the product model ${G}_{1}\square {G}_{2}$, we observe two separated hysteresis areas as well, which are shown in figures 3(c), (d). In particular, starting from the incoherent state on G₁ and G₂ when λ = 0, the product model remains incoherent until ${\lambda }_{c}^{{f}_{2}}=0.726$. Note that when G₂ experiences a jump at ${\lambda }_{c}^{{f}_{2}}$ to the coherent state as zoomed in figure 3(d), however, the product model keeps incoherent since G₁ is in an incoherent state. The forward continuation curve of the product model follows the incoherent solution of G₁ until the critical coupling value of ${\lambda }_{c}^{{f}_{1}}=2.202$ showing a jump to full synchronization. The backward continuation curve of the product model drops to an incoherent state at ${\lambda }_{c}^{{b}_{1}}=0.96$ (figure 3(c)) because G₁ loses synchronization first, while G₂ keeps sync. Reducing the coupling λ further, G₂ loses synchronization at ${\lambda }_{c}^{{b}_{2}}=0.48$, which yields the second drop to an even lower order parameter as zoomed in figure 3(d). From the viewpoint of global behavior, the hysteresis area of the product ${G}_{1}\square {G}_{2}$ is determined by G₁, namely, ${S}_{{G}_{1}\square {G}_{2}}={S}_{{G}_{1}}$.

Zoom In Zoom Out Reset image size

In this case, the cluster synchronized states on the product model are found in both the forward and backward transition processes. In particular, in the coupling regime ${\lambda }_{c}^{{f}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$ of the forward transition, G₁ is in an incoherent state, while G₂ is in a coherent one, which yields the cluster sync states on the product model. In a full analogy in the coupling regime ${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{b}_{1}}$ of the backward transition process, G₁ is again incoherent while G₂ is coherent, showing cluster sync dynamics. The microscopic details of these cluster synchronized states will be further illustrated in section 4.4.

4.2. Case B of figure 2(b)

The second case of figure 2(b) is implemented by choosing the parameters ${\beta }_{2}=13$ and ω₂ = 0.56, which yield ${\lambda }_{c}^{{b}_{2}}=0.48$, ${\lambda }_{c}^{{f}_{2}}=1.293$ and ${S}_{{G}_{2}}=0.813$, namely ${\lambda }_{c}^{{b}_{2}}\lt {\lambda }_{c}^{{b}_{1}}\lt {\lambda }_{c}^{{f}_{2}}\lt {\lambda }_{c}^{{f}_{1}}$. In this case, ${S}_{{G}_{1}}\cap {S}_{{G}_{2}}\ne \varnothing$ as shown in figures 4(a), (b), namely, ${\lambda }_{c}^{{f}_{2}}$ is inside the hysteresis area of G₁. In the forward transition process from an initially incoherent state at λ = 0, the product model ${G}_{1}\square {G}_{2}$ remains an incoherent state until ${\lambda }_{c}^{{f}_{2}}=1.293$ when G₂ experiences the first jump at to a coherent state of $\parallel {z}_{2}\parallel =1$ (as shown in figures 4(c), (d)). Note, however, that when the coupling is further increased in the interval ${\lambda }_{c}^{{f}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$, the ${G}_{1}\square {G}_{2}$ keeps the incoherent state since G₁ is incoherent. The product model shows the second jump at ${\lambda }_{c}^{{f}_{1}}=2.202$ when G₁ is synchronized as well. The backward transition is the same as the case of figure 2(a), namely, the product model loses synchronization by the first jump down at ${\lambda }_{c}^{{b}_{1}}=0.96$, after which G₁ is incoherent but G₂ is coherent. Reducing the coupling further to ${\lambda }_{c}^{{b}_{2}}=0.48$, G₂ loses synchronization leading to the second jump down to an even lower level of the incoherent state of the product model. The global hysteresis area of ${G}_{1}\square {G}_{2}$ is determined by the subgraph G₁, ${S}_{{G}_{1}}=[{\lambda }_{c}^{{b}_{1}},{\lambda }_{c}^{{f}_{1}}]$.

Zoom In Zoom Out Reset image size

In the forward transition process on the product model, the cluster sync states are observed in the coupling regime ${\lambda }_{c}^{{f}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$ since G₁ is incoherent while G₂ is coherent. The cluster sync scenario for the backward transition process is the same as figure 2(a) in the coupling regime of ${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{b}_{1}}$, where G₁ is incoherent but G₂ is coherent.

The stability analysis of steady states of the product model suggests that there are further stable cluster synchronized solutions of ${z}_{{a}_{5}}$ and ${z}_{{a}_{8}}$ in the coupling regime ${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$, which are not easily observed by the previous traditional ways of implementing the forward transition process on both subgraphs G₁ and G₂ simultaneously (or the backward transition process respectively). In contrast, these clustered states are obtained by the product of the forward curve of G₁ with the backward curve of G₂ transition processes, which leads to the solution ${z}_{{a}_{5}}$, or the product of the backward of G₁ with the forward of G₂ transition processes, which leads to ${z}_{{a}_{8}}$. Both cluster sync solutions are, respectively, highlighted by filled diamond and triangle dots in figures 4(c), (d).

4.3. Case C of figure 2(c)

The third case of figure 2(c) is implemented by choosing the parameters ${\beta }_{2}=60$ and ${\omega }_{2}=0.496$, which yield ${\lambda }_{c}^{{b}_{2}}=0.48$, ${\lambda }_{c}^{{f}_{2}}=2.66$ and ${S}_{{G}_{2}}=2.18$, namely ${\lambda }_{c}^{{b}_{2}}\lt {\lambda }_{c}^{{b}_{1}}\lt {\lambda }_{c}^{{f}_{1}}\lt {\lambda }_{c}^{{f}_{2}}$. In this case, ${S}_{{G}_{1}}\subseteq {S}_{{G}_{2}}$ as shown in figures 5(a), (b).

Zoom In Zoom Out Reset image size

This case shows a different forward transition process compared to the previous two cases. In particular, the product model ${G}_{1}\square {G}_{2}$ is incoherent until G₁ experiences a jump at ${\lambda }_{c}^{{f}_{1}}=2.202$ to a higher level of incoherent state that G₁ is coherent but G₂ remains incoherent. When the coupling is further increased to ${\lambda }_{c}^{{f}_{2}}=2.66$, ${G}_{1}\square {G}_{2}$ is fully synchronized because G₂ undergoes the second jump to synchronization. The backward transition process is the same as in the previous two cases of figures 2(a), (b) since both critical values ${\lambda }_{c}^{{b}_{1}}$ and ${\lambda }_{c}^{{b}_{2}}$ are not affected. The product model loses sync by the first jump down at ${\lambda }_{c}^{{b}_{1}}=0.96$ and then by the second jump at ${\lambda }_{c}^{{b}_{2}}=0.48$ to a lower level of incoherent state. The global hysteresis area of ${G}_{1}\square {G}_{2}$ is $[{\lambda }_{c}^{{b}_{1}},{\lambda }_{c}^{{f}_{2}}]$. This size is larger than that of the single G₁ subgraph because of the product effect from G₂.

In the forward transition process, the cluster synchronized states are observed in the coupling regime ${\lambda }_{c}^{{f}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{2}}$ since G₁ is synchronized but G₂ is not. Again because there are no changes for ${\lambda }_{c}^{{b}_{1}}$ nor ${\lambda }_{c}^{{b}_{2}}$, the cluster sync scenario is observed in the coupling regime of ${\lambda }_{c}^{{b}_{2}}\lt \lambda \lt {\lambda }_{c}^{{b}_{1}}$ of the backward transition process. Note again that the product model is obtained by performing the Cartesian operation on the forward continuation curves of G₁ and G₂ simultaneously (or the backward curves).

Furthermore, the stability analysis suggests further cluster synchronized solutions in the coupling regime of ${\lambda }_{c}^{{b}_{1}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$ which are similar to the results as presented in figure 4(d). In this particular coupling interval, we obtain the steady state of ${z}_{{a}_{5}}$ when performing the Cartesian product on the forward continuation process of G₁ with the backward process of G₂. On the other hand, ${z}_{{a}_{8}}$ is achieved by the product of the backward process of G₁ with the forward process of G₂, which are highlighted in figures 5(c), (d).

4.4. Microscopic views of cluster synchronized behavior

We have obtained cluster synchronized states in all three cases above. These states are generated in the coupling regime where one subgraph is synchronized, while the other subgraph is not. In this subsection, we numerically show the microscopic details of these states on the product model.

First, on the product model each node is denoted by ${\left\langle {ik}\right\rangle }_{i=1,k=1}^{{N}_{1},{N}_{2}}$ (equations (5), (6)). For an illustration purpose, we relabel the oscillators as $m=\{[(i-1){N}_{2}+k{]}_{k=1}^{{N}_{2}}\}{}_{i=1}^{{N}_{1}}$, which leads to $m\in [1,{N}_{1}{N}_{2}]$, for instance, m ∈ [1, N₂] corresponds to indices ${\left\langle 1k\right\rangle }_{k=1}^{{N}_{2}}$ and m ∈ [N₂ + 1, 2N₂] is for indices ${\left\langle 2k\right\rangle }_{k=1}^{{N}_{2}}$, etc. In the following numerical example, we choose N₁ = N₂ = 9. Hence there are N₁N₂ = 81 oscillators on the product model. In addition, we report only the case of ${\lambda }_{c}^{{f}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$ of the forward transition (in figure 3) and other cases of figures 4 and 5 show the same clustered states.

We focus on the coupling regime when G₁ is not synchronized while G₂ is synchronized, which is implemented by random ICs for G₁ while identical ICs for G₂. Namely, the oscillators ${\varphi }_{i}^{\left(1\right)}(i=1,\cdots ,{N}_{1})$ are not synchronized, but the oscillators ${\varphi }_{k}^{\left(2\right)}(k=1,\cdots ,{N}_{2})$ are synchronized. When implementing the Cartesian product operation, the phase differences between the ith oscillator ${\varphi }_{i}^{\left(1\right)}$ on G₁ and all other oscillators ${\varphi }_{k}^{\left(2\right)}(k=1,\cdots ,{N}_{2})$ on G₂ are fixed to the same value, all of them are locked to the phase locking manifold [13]. Therefore, N₁ different ${\varphi }_{i}^{\left(1\right)}$ oscillators on G₁ lead to N₁ clusters of synchronized states, which are shown by the temporal phase profile in figure 6(a). The N₁ cluster synchronized state has been observed for the coupling regime ${\lambda }_{c}^{{f}_{2}}\lt \lambda \lt {\lambda }_{c}^{{f}_{1}}$ as shown in figure 6(b).

Zoom In Zoom Out Reset image size

If the coupling strength is in a regime that G₁ is synchronized while G₂ is not, we obtain an N₂ clustered synchronization state. Again for illustration purpose, we relabel the oscillators as $n=\{[(k-1){N}_{1}+i{]}_{i=1}^{{N}_{1}}\}{}_{k=1}^{{N}_{2}}$ which helps to visualize the Cartesian product operation in the following. On the product model, the phase differences between the kth oscillator ${\varphi }_{k}^{\left(2\right)}$ on G₂ and all other oscillators ${\varphi }_{i}^{\left(1\right)}(i=1,\cdots ,{N}_{1})$ on G₁ are constants, forming one clustered state. Therefore, N₂ asynchronous phase oscillators yield N₂ clustered synchronization states on the product model.

In summary, we have provided the WS ansatz to the lower dimensional dynamic equations which describe the high dimensional Cartesian product model that is reconstructed from two independent Kuramoto models on star subgraphs of G₁ and G₂ by graph operation. The order parameters describing different synchronization states of the product model are expressed by fixed points of the lower dimensional equations. Furthermore, the steady states of the order parameter Z and their respective stabilities have been delineated theoretically. Our numerical simulations agree very well with the theoretical results.

In the case of a single star graph, there is only one discontinuous explosive transition to synchronization in the forward process. In contrast, two explosive jumps in the forward curve have been observed in the product model. The first jump corresponds to a local scale of synchronization of one subgraph only while the other subgraph is not synchronized. The second jump is for the global synchronization for both subgraphs. Between these two jumps, cluster synchronized behavior has been widely obtained, which provides complementary insights for the understanding of cluster synchronization. In the literature, many versions of cluster synchronization scenarios have been reported in various settings, for instance, for unidirectional coupling, time delays and some special network structures [41, 42]. Some numerical algorithms are required to identify synchronized clusters [43] or using graph partitions [44]. Recently, computational group theory has been proposed to characterize cluster synchronization, which hinges on the decomposition of the networked nodes into clusters with the help of network symmetries [25, 26].

In contrast to the literature, cluster synchronized states in this work are reconstructed by the Cartesian product operation which is performed from two independent star networks of phase oscillators. In the product model ${G}_{1}\square {G}_{2}$, the clustered states are widely observed, especially in the coupling regimes where one factor graph is synchronized while the other factor graph is not. Note that such clustered states are realized for the case that both G₁ and G₂ are in a forward continuation transition processes (or both are on the backward processes). Furthermore, the linear stability analysis of fixed points identifies further cluster synchronized states that are realized by the Cartesian product of the forward transition process of G₁ and the backward process of G₂ (or the vise versa). We emphasize that the cluster synchronization solutions are synthetic states which are obtained by graph product operations. Furthermore, we easily get the number of clusters which is determined by the number of asynchronous oscillators of one factor graphs. As it has been demonstrated in [26], there are six symmetries in a star of identical oscillators. One interesting but maybe challenging task is to study how these group symmetries change when performing Cartesian operations from two subgraphs of stars.

In our earlier work of the Cartesian product model on a star and a ring subgraphs [24] (e.g. G₁ is a star and G₂ is a ring), we focused on disclosing a hybrid state of an explosive forward synchronization transition and a continuous backward desynchronization transition. Furthermore, the critical coupling thresholds for synchronization transitions are obtained by the necessary conditions of synchronized solutions in the linearized equations. In contrast, the WS method provides lower dimensional nonlinear equations for the ensemble order parameter Z. However, in the product of a star and a ring subgraphs, the WS dimension reduction technique can not be applied straightforwardly since the complex common driving force of equation (11) can not be written down explicitly. Therefore, it remains to be challenging to obtain lower dimensional equations for such a case of Cartesian product of arbitrary subgraph structures.

The Cartesian product of graphs is a commutative, associative binary operation on graphs [21] and most of properties can be derived from the factors. Furthermore, in this work a single star graph is modeled by a population of identical units when introducing the phase difference between the hub and leaf nodes. Therefore, dimension reduction techniques like WS and OA can be applied straightforwardly. The Cartesian product model of two star graphs provides an easy way to build a big graph of the Kuramoto phase dynamics. The product operation can be further performed recursively for G₁ and G₂ or based on more than two subgraphs (i.e. G_n, n ≥ 3), which is one of the interesting topics for future work. In such cases, we expect that cluster synchronized states on the product can be easily constructed since more combinations of synchronization and desynchronization transition processes of subgraphs are involved. Therefore, graph product operations provide a 'bottom-up' framework which may generate much complicate dynamics that are observed in complex systems.

Parts of this work have been financially supported by the National Natural Science Foundation of China (Grant Nos. 11872182, 11875132, and 11835003) and the Natural Science Foundation of Shanghai (Grant Nos. 17ZR1444800 and 18ZR1411800).

All Möbius transformations are provided below that are necessary for equations (13), (14)

$$\begin{eqnarray}&&{{\rm{e}}}^{{\bf{i}}(2{\varphi }_{{i}}^{(1)}(t)+{\varphi }_{k}^{(2)}(t))}={({{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}(t)})}^{2}{{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}(t)}=\displaystyle \frac{{[{z}_{1}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}]}^{2}}{{[1+{z}_{1}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}]}^{2}}\displaystyle \frac{{z}_{2}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}}{1+{z}_{2}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}},\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{{\rm{e}}}^{{\bf{i}}({\varphi }_{{i}}^{(1)}(t)+2{\varphi }_{k}^{(2)}(t))}={{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}(t)}{({{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}})}^{2}=\displaystyle \frac{{z}_{1}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}{1+{z}_{1}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}\displaystyle \frac{{[{z}_{2}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}]}^{2}}{{[1+{z}_{2}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}]}^{2}},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{{\rm{e}}}^{{\bf{i}}({\varphi }_{{i}}^{(1)}(t)+{\varphi }_{k}^{(2)}(t))}={{\rm{e}}}^{{\bf{i}}{\varphi }_{{i}}^{(1)}(t)}{{\rm{e}}}^{{\bf{i}}{\varphi }_{k}^{(2)}(t)}=\displaystyle \frac{{z}_{1}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}{1+{z}_{1}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}}\displaystyle \frac{{z}_{2}(t)+{{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}}{1+{z}_{2}^{* }(t){{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}}.\end{eqnarray} \tag{ A.3 }$$

Note that the variables of ξ_i,k are constants. In addition, in the case of uniform distribution of ξ_i,k, the global variables α_1,2 and ξ_i,k do not enter the equation for z_1,2 and Z [16, 17]. Using these transformations, we can obtain the equation (16) by equalizing the non-exponential terms of the two sides of the equations (13), (14). More specifically, the non-exponential term of the left-hand side of equation (13) reads

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{z}_{1}}{z}_{2}+\dot{{z}_{2}}{z}_{1}}{{[1+{z}_{1}^{* }{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}]}^{2}{[1+{z}_{2}^{* }{{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}]}^{2}},\end{eqnarray} \tag{ A.4 }$$

and, respectively, that of the right-hand side of equation (14) is

$$\begin{eqnarray}&&\displaystyle \frac{{\bf{i}}f({z}_{1}^{2}{z}_{2}+{z}_{1}{z}_{2}^{2})+{\bf{i}}({g}_{1}+{g}_{2}){z}_{1}{z}_{2}+{\bf{i}}{f}^{* }({z}_{1}+{z}_{2})}{{[1+{z}_{1}^{* }{{\rm{e}}}^{{\bf{i}}({\xi }_{{i}}^{(1)}+{\alpha }_{1}(t))}]}^{2}{[1+{z}_{2}^{* }{{\rm{e}}}^{{\bf{i}}({\xi }_{k}^{(2)}+{\alpha }_{2}(t))}]}^{2}}.\end{eqnarray} \tag{ A.5 }$$

By equalizing the above two terms, we obtain

$$\begin{eqnarray}&&\dot{{z}_{1}}{z}_{2}+\dot{{z}_{2}}{z}_{1}={\bf{i}}f({z}_{1}^{2}{z}_{2}+{z}_{1}{z}_{2}^{2})+{\bf{i}}({g}_{1}+{g}_{2}){z}_{1}{z}_{2}+{\bf{i}}{f}^{* }({z}_{1}+{z}_{2}).\end{eqnarray} \tag{ A.6 }$$