Cartesian product of synchronization transitions and hysteresis

We present theoretical results when applying the Cartesian product of two Kuramoto models on different network topologies. By a detailed mathematical analysis, we prove that the dynamics on the Cartesian product graph can be described by the canonical equations as the Kuramoto model. We show that the order parameter of the Cartesian product is the product of the order parameters of the factors. On the product graph, we observe either continuous or discontinuous synchronization transitions. In addition, under certain conditions, the transition from an initially incoherent state to a coherent one is discontinuous, while the transition from a coherent state to an incoherent one is continuous, presenting a mixture state of first and second order synchronization transitions. Our numerical results are in a good agreement with the theoretical predictions. These results provide new insight for network design and synchronization control.


Introduction
Synchronization in nonlinear sciences is an emergent property that occurs in a broad range of dynamical systems, including neural signaling, the beating of the heart, fire-fly light waves, or power grids [1]. Kuramoto phase model is a paradigmatic example in synchronization analysis [2,3]. This model consists of self-sustained phase oscillators rotating at heterogeneous intrinsic frequencies coupled through the sine of their phase differences [4,5]. A transition from an initially incoherent state to a fully coherent state takes place when the coupling strength exceeds a critical threshold, which explains many collective behaviors in nature, science, society and technology [1,3]. There have been rapid developments in the study of the Kuramoto model focusing on the effects of network structures on synchronization, e.g., from the traditional all-to-all coupling to heterogeneous complex network topologies [6,7]. In the case of complex network structures, several models have been proposed to study effects of small-world structures, communities, degree correlations and multi-layer network structures [8]. Recent findings of the Kuramoto model on different network topologies have been reviewed in [9].
It is important to emphasize that most likely continuous synchronization transitions can be observed in these Kuramoto models on top of networks. Namely, the order parameter which characterizes the degree of synchronization grows continuously when the coupling strength passes the critical threshold. In other words, when the coupling strength is increased progressively, the sizes of the synchronized clusters grow gradually [10]. The findings of discontinuous phase transition to synchronization (also known as abrupt explosive synchronization) have triggered several rapid investigations [11][12][13][14], which is a consequence of correlations between network structure and local dynamics. The most intriguing phenomenon is that hysteresis has been largely observed in explosive synchronization [11]. More specifically, the network shows an explosive jump from an incoherent state to a coherent one when the coupling strength is increased adiabatically, which is conveniently called forward continuation curve below. In addition, it shows a sudden drop from the coherent state to the incoherent state when the coupling strength is decreased progressively (backward continuation curve). There is a clear hysteresis area because these two curves (forward and backward) do not overlap. In the case of continuous synchronization transitions, the forward and backward curves are overlapped. Hysteresis is a fundamental property of a first order phase transition, which shares some analogies with explosive percolation [15]. The hysteresis behavior at the onset of synchronization has been widely observed in various situations, such as scale-free networks [11,16], electronic circuits [12,17], time delayed systems [18], and a second order Kuramoto model [19]. Recent studies also focused on the importance of the frequency distribution [20], noise effects [21], various generalizations of the coupling patterns [22,23], and multi-layer networks [14].
Usually, the continuous and discontinuous synchronization transitions are separately discussed in the literature [15], because of the different requirements of network configurations and local dynamics. Here, we propose an algorithm based on graph product operation, which allows for discussing continuous and discontinuous synchronization transitions in a unified framework. In particular, we construct the Kuramoto model by the Cartesian product, which is one of the basic operations on a graph [24]. This graph operation helps to find hysteretic synchronization transitions with more intriguing phenomena. In particular, we find that the forward transition curve is discontinuous while the backward transition is continuous. Note that the effects of graph operations on the synchronization behavior remained largely unclear, although some discussions on the synchronizabilities have been presented in terms eigenspectra [25,26]. To the best of our knowledge, our work is the first attempt to address the effects of the Cartesian product on the synchronization behavior of the Kuramoto model, resulting in a mixture state of both continuous and discontinuous transitions.
Let us start by considering the Kuramoto model on a complex network. We consider n phase oscillators on top of a complex network G in the following framework without loss of generality where i w are the natural frequencies taken from a certain distribution, λ is the coupling strength, and A ij is the adjacency matrix of the network. Oscillator i is coupled with j if there is a link between node i and j, namely, A 1 ij = . Otherwise, A 0 ij = means that i and j are not connected. The Kuramoto order parameter R is defined as where the notation t |·| denotes the time average of the absolute value over t 1  . Small values of R indicate incoherent behavior. In contrast, as R 1  we encounter a highly coherent state.

Product of Kuramoto models
We start by introducing the method of Cartesian product of graphs, which is an important way to construct a bigger graph and plays an important role in network design and analysis.
Given two nonempty graphs G V E ,  of the two graphs is a graph such that: (i) the vertex set of G G 1 2  is the Cartesian product V V 1 2 . For example, given vertices i of G 1 and k of G 2 , we denote the vertex of G G 1 2  as ik ; á ñ (ii) two vertices ik á ñ and jl á ñ are connected in G G 1 2  if and only if (a) i=j and k is adjacent with l in G 2 , or (b) k=l and i is adjacent with j in G 1 . The graphs G 1 and G 2 are called factors of the product G G 1 2  [24]. In figure 1(a), we give an example of the Cartesian product of two chains and the resulting graph is a regular lattice. Further examples of the Cartesian product graphs are illustrated by two factor subgraphs of stars (Figure 2(a)), two rings ( Figure 2(e)), and one star and one ring (figure 2(i)).
The Cartesian product of graphs is a commutative, associative binary operation on graphs [27]. It has many useful properties, most of which can be derived from the factors. 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, A ) ( ) ( )(see appendix A for details on the Kronecker product and sum of matrices).
Theorem. Given two independent Kuramoto models on factor graphs G 1 and G 2 (n 1 oscillators on G 1 and n 2 oscillators on G 2 , respectively), the definition of the phase of the node ik á ñ as ik Proof. An example of constructing Cartesian product from two chain models is shown in figure 1(b). For two independent graphs G 1 and G 2 , the adjacency matrices are A 1 ( ) and A 2 ( ) respectively. The Kuramoto models on the respective factors are written as   of (a): two stars; (e) two rings; and (i) G 1 is a star and G 2 is a ring. Panels (b), (f) and (j) are order parameters for synchronization transition curves on the factor G 1 , where the red curve is the forward curve and the blue is the backward curve. The similar curves for factor G 2 are shown in panels (c), (g) and (k). Panels (d), (h) and (l) are order parameters for synchronization transition curves on the product G G The result of equation (9) has exactly the same expression as the equation of the Kuramoto model. Furthermore, it is easy to recognize that the natural frequencies are ik i k There are no interactions between the two factors G 1 and G 2 and the dynamics of oscillators on G 1 evolves independently with oscillators on G 2 . The oscillator dynamics on the Cartesian product graph G G 1 2  is reconstructed by means of the summation of the respective phases on G 1 and G 2 . The advantage of defining the phase of the node on the product graph as the summation of the respective phases of the factor subgraphs yields the canonical equations of the Kuramoto model on G G 1 2  . In addition, the Cartesian product operation can be easily generalized to the case of n subgraphs because G G 1 2  and G G 2 1  are isomorphic (the operation is commutative). Furthermore this operation is associative. For instance, given three factor subgraphs G 1 , G 2 and G 3 , the phase summation of three factors results in the traditional equations as the Kuramoto model on the product graph as With the commutative and associative properties of the phase summation, we generalize the theorem to the case of n factor subgraphs straightforwardly.
Lemma. The order parameter R on the Cartesian product graph of G G 1 2  is the product of the order parameters R 1 and R 2 of two independent factors, namely, R R R 1 2 = .
Proof. The order parameter R on the Cartesian product graph is computed as follows: Therefore, we prove that the order parameter R is fully determined by the two factors R 1 and R 2 . Furthermore, R does not converge if any of the order parameters of the two factors (either R 1 or R 2 ) does not converge. In addition, for n factor subgraphs, it is straightforward to show that the order parameter R on the product graph is the order parameters of n factors. For instance, given three factors G 1 , G 2 and G 3 , the order parameter R of the product G G G

Kuramoto model on star and ring structures
Before showing the main results of this work, we discuss the synchronization transitions when implementing the Kuramoto models on single star and ring structures (see appendix E for further results on chain structures and appendix F for trees). Both stars, rings, chains and trees of oscillators have been studied extensively in synchronization analysis because they are considered to be motifs to build complex networks [28][29][30][31][32][33][34][35]. Here, we show the critical coupling threshold values for synchronization on these topologies. First, taking into account different network topologies, we rewrite the Kuramoto model (equation (1)) in terms of link representation in the following. Assume that there are n nodes and m links in a network which is represented by the adjacency matrix A n n  Î´. In general, the link incidence matrix B m n  Î´of a directed network is defined as is the starting node of link , 0 otherwise, 1 is the end node of link , where α is the link and i is the oscillator. Note that no difference appears when using the above matrix B to represent an undirected network. Therefore, the Kuramoto model (equation (1)) is rewritten in the following compact matrix form and T is the transpose operation to the corresponding vectors. Furthermore, the sine function acts on the phase vector one by one, e.g., we define sin , , sin , sin , sin ) . It is straightforward to show that the compact form representation of the Kuramoto model (equation (15)) is identical to the canonical equations (equation (1)) (in appendix B, we show the equivalence for the two expressions). We then compute the average phase of all oscillators as Synchronization is achieved when the phase difference between all oscillators is a constant. Therefore, the synchronization condition is equivalently represented by T T q w q q q w w w = =   (¯¯¯)¯(¯¯¯) representing each oscillator has the same frequency as the population averaged frequency. Thus, synchronization suggests that there is a solution q in the following equation: In appendix D, we provide detailed discussions on the existence and uniqueness of the solutions to the linear equation (18). In addition, we prove that a solution always exists for connected networks. We summarize the results for two typical networks in terms of the following corollaries: In addition, we have which is the same as reported in [13,36].
Remark 3. The Kuramoto model on either a star or a chain can be generalized to a more general framework of a tree structure. The necessary condition to synchronization of the Kuramoto model on a tree is max Proof. The incidence matrix is  permutation matrix P that shifts the last entry to the beginning, the second row is denoted as bP T . With this notation, we re-write the Moore-Penrose pseudo-inverse as In addition, we have Note that the summation of the phase difference between two neighboring oscillators along the ring structure is zero. We consider a special case that n is an even number and i i n2 w w = -+ , which is equivalent to that P n 2 w w = -. This means that summation of all terms of arcsin B T w l + is zero. Therefore, a necessary condition for synchronization is where the maximum function runs over all links α. é

Product of synchronization transitions
In the following, we first summarize the results of the synchronization transitions of the Kuramoto model on a single network (star and ring).
(i) For n oscillators coupled in a star structure. We consider the special case of explosive transitions to synchronization [13,36]. In particular, we assume that the hub node has the frequency n 1 ) and all leaf nodes have the same frequency In this work, we focus on the dynamics of the Kuramoto model when performing the Cartesian product from two factor graphs. It is certainly an interesting topic to discuss the structural properties of the resulting product graph, which however is outside the scope of the present work. Given two independent networks of phase oscillators on G 1 and G 2 (e.g., G 1 is a star and G 2 is a ring), the synchronization transitions on the Cartesian product G G 1 2  graph are summarized in the following: for G 2 . The choice of natural frequencies yields different critical coupling thresholds for G 1 and G 2 , which are listed in table 3 and are further illustrated in figure 2. In particular, we observe discontinuous transitions in both forward and backward directions if the product graph G G 1 2  is obtained from two star networks, where we find a clear hysteresis as shown in figures 2(a)-(d). If the product G G 1 2  is reconstructed from two ring networks, only Table 1. Hysteresis area l D on G G  if factor G 1 is a star and G 2 is a ring.

Sync type Only discontinuous transitions
Discontinuous transition for the forward curve and continuous transition for the backward curve Only continuous transitions continuous transitions are observed in figures 2(e)-(h). In addition, if G 1 is a star and G 2 is a ring, the product graph G G 1 2  presents a mixture state of continuous transition for the backward direction and a discontinuous one for the forward direction, showing hysteretic behavior as shown in figures 2(i)-(l). In all these cases, the theoretical predictions for synchronization transitions are perfectly validated by our numerical simulations, which are highlighted by vertical dashed lines.

Discussion
In summary, we have provided theoretical analysis on effects of graph operations on synchronization transitions in the Kuramoto models, in particular with the Cartesian product. Given two independent Kuramoto models on different topologies, where the individual systems present either continuous or discontinuous transitions to synchronization. Taking into account different network topologies (star, chain, ring, and tree), we proposed a unified framework in obtaining necessary conditions for synchronization in each network structure, which is based on a detailed analysis of the existence of solutions to linear equations (18). Under the phase summation assumption, the results as summarized in the theorem and lemma can be easily generalized to the case of n factor subgraphs because of the commutative and associative properties of the Cartesian product operation.
Rich synchronization transitions have been obtained for the Cartesian product graph. Depending on the relation between the critical coupling thresholds, there are three different synchronization scenarios on the Cartesian product graph: (i) both the forward and backward transitions are discontinuous with a clear hysteresis area, (ii) both the forward and backward transitions are continuous and these two curves are overlapped and, hence, without hysteresis, and (iii) the forward transition from an incoherent state to a coherent one is discontinuous and the backward transition from a coherent state to an incoherent one is continuous with a hysteretic behavior.
In this work, our theoretical approaches are mainly focused on star, chain and ring structures because the necessary conditions of synchronization transitions have been analytically obtained in the same framework. From the viewpoint of numerical simulations, the phase summation operation on the Cartesian product graph can be performed for general network settings as well. For instance, given two independent factors G 1 and G 2 presenting continuous synchronization transitions at respective critical couplings c1 l and c2 l , the synchronization on the product graph G G 1 2  is fully determined by the product R R R 1 2 = . Therefore, the critical coupling on G G 1 2  is max , The Cartesian product can be performed recursively, for instance, by Kronecker power, which is one of the graph operations in generating a big graph from two or more small factor graphs [37]. This generative model presents some properties that are often observed in real networks, e.g., in heavy-tailed degree distributions and small diameters etc. From the viewpoint of illustration, it shares some similarities with generating multi-layer networks [8]. Therefore, our results provide some novel insight in network design and synchronization control.  Appendix B. Equivalence between equations (15) and (1) Proof. We show the coupling terms of both representations are identical. Suppose two phase oscillators attached to the link α are denoted as a + and a -, where a + is the starting point while ais the end point respectively. Thus we have the notation Bq Therefore, the coupling term for the node i is expanded as If oscillator i is the starting point i a = + , and j is the end point j a = - On the other hand, if i is the end point i a = -, and j is the starting point j a = -, we have Therefore the coupling term of node i is simplified as This proves the statement that the coupling terms are identical in equations (15) and (1). é

Appendix C. Moore-Penrose pseudo-inverse
Here we give a brief introduction to the Moore-Penrose pseudo-inverse, a generalization of the inverse of a matrix. The Moore-Penrose pseudo-inverse is defined for any matrix and is unique, which brings conceptual clarity to study the solutions of arbitrary systems of linear equations.  Solutions of q do not always exist because B is not row full rank. In addition, the summation of phase difference between two neighboring oscillators along the loop structure is zero, which suggests that the choice of u can not be arbitrary.
If synchronization exists (