Fully solvable lower dimensional dynamics of Cartesian product of Kuramoto models

Implementing a positive correlation between the natural frequencies of nodes and their connectivity on a single star graph leads to a pronounced explosive transition to synchronization, additionally presenting hysteresis behavior. From the viewpoint of network connectivity, a star has been considered as a building motif to generate a big graph by graph operations. On the other hand, we propose to construct complex synchronization dynamics by applying the Cartesian product of two Kuramoto models on two star networks. On the product model, the lower dimensional equations describing the ensemble dynamics in terms of collective order parameters are fully solved by the Watanabe–Strogatz method. Different graph parameter choices lead to three different interacting scenarios of the hysteresis areas of two individual factor graphs, which further change the basins of attraction of multiple fixed points. Furthermore, we obtain coupling regimes where cluster synchronization states are often present on the product graph and the number of clusters is fully controlled. More specifically, oscillators on one star graph are synchronized while those on the other star are not synchronized, which induces clustered state on the product model. The numerical results agree perfectly with the theoretic predictions.


Introduction
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][4][5][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][11][12][13][14][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 l 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 l c b . In addition, a clear hysteresis area has been observed between these synchronization and desynchronization processes, which is further denoted as l l = -S c f 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][28][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.

Cartesian product of two Kuramoto models on stars
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 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 1,2 , 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 In addition, the order parameters of the two independent networks are defined as å å 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

 
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 1 and G 2 are schematically shown in figure 1. Meanwhile, we use the notation á ñ ik to represent the index of node on the product  G G 1 2 . In addition, on the Cartesian product graph, the phase of the node á ñ ik is defined as j j j = + á ñ ik i k 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 G 1 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]. In addition, the order parameter of the Cartesian product model is defined as 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 1 (t) z 2 (t). Namely, the order parameter on  G G 1 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 1 , G 2 and G 3 , the order parameter R of the product . With the above phase definition, the time derivatives of the phases j á ñ ik on the Cartesian product where f * are the complex conjugates of f. The f and g 1,2 are further defined as Note that, in equation (9), the term (1−β 1 )ω 1 +(1−β 2 )ω 2 represents the natural frequency of the oscillator on the product  G G 1 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 1 +g 2 are driven by the arbitrary complex , 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 1 , z 2 and Z. Therefore, we first consider the relationship 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 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 1,2 and the uniform distribution of constants of motion ξ i,k , the WS variables z 1 , z 2 coincide with the local Kuramoto mean-field, namely yielding the Ott-Antonsen (OA) ansatz [38][39][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).
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 1 z 2 . Then, we have = + Z z z z z Therefore, in terms of the global variables Z, the phase equation equation (11) is expressed as Considering the definition of the order parameter = F Z re i (equation (7)), we obtain the following two nonlinear coupled equations for the global order parameter in the complex plane

Steady states of the product model
There are three steady states of the nonlinear coupled system (equations (18), (19)), namely, 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 1 z 2 . To this end, we first denote the forward critical coupling threshold values for the two independent factors following the notations of [16]: and two backward critical coupling thresholds 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 which further ensure that the hysteretic areas exist for both independent star graphs, namely, 3.1. Full synchrony for the product model r 1 =r 2 =1 In this case, equations (18), (19) are simplified as . 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 l c b 1 and l c b 2 , respectively. In this case, the solutions of z 1 and z 2 are easily obtained as   We focus on the discussion on the case of r 1 =1 and F =  p 2 2 , while the symmetric case of r 2 =1 and F =  p 1 2 shows a full analogy. In this case, we have the following two equations in order to obtain the steady states, namely The order parameter Z for the product model is Z=z 1 z 2 again with different expressions depending on β 1,2 and ω 1,2 , which will be summarized in section 3.4.

Intermediate summary of steady states
Note that the critical coupling thresholds l l l , , 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 1,2,3,4 correspond to four synchronous states, showing different stability conditions as we summarize in the next section.

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 where f and g are the right-hand side of the system (equations (18), (19)). The elements of J are the following: 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 s1,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), ) is a stable sink if both z s 1 and z s 2 are positive real values, and --Z s2 ( )is an unstable source if both z s 1 and z s 2 are negative real Table 1. Fixed points of the order parameters with their stability and meaning for figure 2(a).

Fixed point Stability
Existence region Physical meaning Coherent state Cluster sync Za 6 Saddle 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 Cluster sync Za 6 Saddle Cluster sync values. On the other hand, )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 s4 ( ) 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.

Numerical results
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 5 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 = å = z zi . 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 1 and G 2 (equations (1), (2)) simultaneously, while phase dynamics of the product G 1 G 2 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 1 +N 2 phase oscillators in the former case, yielding a better computation efficiency than that of the latter case that is integrated for N 1 N 2 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 G 1 2 is the Kronecker sum of the adjacency matrices of G 1 and G 2 , namely, 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 1 =N 2 =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 1 do not vary among the three cases of figures 2(a)-(c), namely, β 1 =9, ω 1 =1.2, which yield l = 0.96 . In the case of G 2 , we choose β 2 and ω 2 such that l c b 2 is fixed as 0.48, but l 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. showing a jump to full synchronization. The backward continuation curve of the product model drops to an Table 3. Fixed points of the order parameters with their stability and meaning for figure 2(c).

Fixed point Stability
Existence region Physical meaning Note again that the product model is obtained by performing the Cartesian operation on the forward continuation curves of G 1 and G 2 simultaneously (or the backward curves). Furthermore, the stability analysis suggests further cluster synchronized solutions in the coupling regime of l l l < < 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 1 with the backward process of G 2 . On the other hand, z a 8 is achieved by the product of the backward process of G 1 with the forward process of G 2 , which are highlighted in figures 5(c), (d).

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 á ñ = = ik i k (equations (5), (6)). For an illustration purpose, we relabel the oscillators as = -  We focus on the coupling regime when G 1 is not synchronized while G 2 is synchronized, which is implemented by random ICs for G 1 while identical ICs for G 2 . Namely, the oscillators j = i N 1, ,  If the coupling strength is in a regime that G 1 is synchronized while G 2 is not, we obtain an N 2 clustered synchronization state. Again for illustration purpose, we relabel the oscillators as = - ] } which helps to visualize the Cartesian product operation in the following. On the product model, the phase differences between the kth oscillator j k 2 ( ) on G 2 and all other oscillators on G 1 are constants, forming one clustered state. Therefore, N 2 asynchronous phase oscillators yield N 2 clustered synchronization states on the product model.

Discussion
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 1 and G 2 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 G 1 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 1 and G 2 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 1 and the backward process of G 2 (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 1 is a star and G 2 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 1 and G 2 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.