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. 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 asRecently, 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 ...