Explosive synchronization in coupled stars

Ruby Varshney¹, Kaustubh Manchanda^2,3 and Haider Hasan Jafri¹ ¹Department of Physics, Aligarh Muslim University, Aligarh 202 002, India
²School of Arts and Sciences, Azim Premji University, Bengaluru, Karnataka 562 125, India
³College of Arts and Sciences, Abu Dhabi University, Abu Dhabi 59911, UAE.

Abstract

We study the effect of network topology on the collective dynamics of an oscillator ensemble. Specifically, we explore explosive synchronization in a system of interacting star networks. Explosive synchronization is characterized by an abrupt transition from an incoherent state to a coherent state. In this study, we couple multiple star networks through their hubs and study the emergent dynamics as a function of coupling strength. The dynamics of each node satisfies the equation of a Kuramoto oscillator. We observe that for a small inter-star coupling strength, the hysteresis width between the forward and backward transition point is minimal, which increases with an increase in the inter-star coupling strength. This observation is independent of the size of the network. Further, we find that the backward transition point is independent of the number of stars coupled together and the inter-star coupling strength, which is also verified using the Watanabe and Strogatz (WS) theory.

I Introduction

Synchronization is an important phenomenon observed in the study of networked dynamical systems Pikovsky-Book ; Boccaletti-Book . In the context of an ensemble of interacting units, synchronization refers to a transition from an incoherent state to a coherent state. Such an onset of coherence in real systems has been reported in power grids dorfler2012 ; motter2013 , neuronal networks cobb1995 , communication networks ling , and circadian rhythms winfree1987 . Recently, special attention has been paid to understand the properties of complex networks consisting of dynamical units on the nodes Watts-1998 ; Barabasi-1999 ; Boccaletti-2006 . These studies have revealed that the underlying topology plays a crucial role in the onset of synchronization Boccaletti-2006 ; Moreno-2004 ; Zhou-2006 ; arenas . Most studies have reported that the transition from desynchrony to synchrony is second-order in nature. This is characterized by an order parameter showing a smooth transition with changing coupling strength.

With the work of Gardenes et al. gardenes , it was observed that in the case of a scale-free network, the order parameter might change discontinuously, resulting in the so-called explosive synchronization (ES). This transition is characterized by a well-defined hysteresis associated with the backward and forward transitions. It has been found to occur as a result of the correlation between the natural frequency of the dynamical units and the network topology in the case of a scale-free network. Many networks in nature are heterogeneous and have scale-free topology where a few nodes, namely the hubs, hold most connections. In such a scenario, a typical motif is a star network. There have been several studies to understand the collective phenomenon in the case of star networks, namely, the chimera state meena ; muni , pattern formation nikos ; shena and explosive synchronization gardenes ; coutinho ; zou ; pavel ; vlasov ; yusra2024 . In Ref. gardenes ; vlasov ; yusra2024 , ES in a single star has been studied by introducing the degree-frequency correlation.

However, network synchronization can be influenced by the presence of other networks. In many realistic situations, it was observed that the network may display several levels of topological organization, which is not exactly described by typical network models Song-2005 . In the case of the cerebral cortex of mammalian brains, the connectivity is organized as subnetworks of interacting neurons Zhou-2006 ; Gardenes-2010 . In transportation and telecommunication systems, the units that serve as transshipment and switching points may be regarded as a system of interacting hubs Campbell-1994 ; Kim-2008 . Thus, it will be interesting to understand the collective dynamics in a situation wherein the underlying topology consists of a few central elements, namely the hubs that are interacting with each other. To simplify the situation, one can consider a star network that evolves under the influence of other star networks xu ; zhu2021 . Collective dynamics in the case of multilayer networks, where one layer might affect the dynamics of the other layer, has also been explored in Refs. boccaletti2014 ; leyva ; jalan23 ; zhang .

In this study, we consider a network of two or more coupled star networks and study the phenomenon of explosive synchronization in the presence of degree-frequency correlation. We derive analytically a criterion for the synchronization of two identical star networks. This analysis is further extended to study the transitions in case of $N$ coupled star networks. The analytical results are then compared with the numerical findings and we observe that the numerical and analytical results are in good agreement with each other. We note that the backward transition point is independent of the coupling strength and the size of the network; however, we predict an upper bound for the forward transition points. Throughout the investigation, we make use of the fourth-order Runge-Kutta method to solve the differential equations.

This report is organized as follows: In Sec. II, we study the dynamics in the case of two coupled star networks. We also give analytical support to the numerical findings for two coupled stars. In Sec. III, we describe the collective dynamics in situations where we have three or more star networks interacting through their hubs. We give an analytical prescription to find the backward transition in the case of the $N$ -coupled star networks. Finally, we summarize the findings in Sec. IV.

Refer to caption — Figure 1: Two-star networks coupled via hubs with the inter-star coupling strength ( $\varepsilon$ ).

II Two coupled star networks

II.1 Numerical Results

Consider a system of two coupled star networks (Fig. 1) interacting through their hubs. The dynamics of each node is modelled via a Kuramoto oscillator. The equations of motion for the system are given by

$\displaystyle\dot{\phi}_{j}$	$\displaystyle=\omega+\lambda_{1}\sin(\phi_{h}-\phi_{j}),\qquad 1\leq j\leq N_{1}$	(1)
$\displaystyle\frac{1}{\beta}\dot{\phi}_{h}$	$\displaystyle=\omega+\frac{\lambda_{1}}{N_{1}}\sum_{j=1}^{N_{1}}\sin(\phi_{j}-\phi_{h})+\varepsilon\sin(\phi^{\prime}_{h}-\phi_{h}),$
$\displaystyle\dot{\phi}^{\prime}_{j}$	$\displaystyle=\omega+\lambda_{2}\sin(\phi^{\prime}_{h}-\phi^{\prime}_{j}),\qquad 1\leq j\leq N_{2}$
$\displaystyle\frac{1}{\beta}\dot{\phi}^{\prime}_{h}$	$\displaystyle=\omega+\frac{\lambda_{2}}{N_{2}}\sum_{j=1}^{N_{2}}\sin(\phi^{\prime}_{j}-\phi^{\prime}_{h})+\varepsilon\sin(\phi_{h}-\phi^{\prime}_{h}).$

where $\phi_{j}$ ( $\phi^{\prime}_{j}$ ) and $\phi_{h}$ ( $\phi^{\prime}_{h}$ ) are the phases of the oscillator at the nodes and at the hub, respectively, of the first (second) star. $\lambda_{i}$ ’s are the intra-star coupling strengths and $\varepsilon$ denotes the inter-star coupling strength. All the oscillators have the same frequency i.e. $\omega=1$ , and the mismatch parameter is set to $\beta=10$ . $N_{1}$ and $N_{2}$ represent the number of nodes in each star. To understand the collective dynamics of the system, we have calculated the global order parameter $R$ , given by

\displaystyle\centering R\@add@centering

\displaystyle=

\displaystyle\frac{1}{N}\left\lvert\sum_{j=1}^{N}e^{i\phi_{j}}\right\rvert,

(2)

where $N$ ( $=N_{1}+N_{2}+2$ ) is the total number of oscillators in the system. We consider a situation where both the stars have identical intra-star coupling strengths, i.e. $\lambda_{1}=\lambda_{2}=\lambda$ and study the variation of the global order parameter $R$ as a function of $\lambda$ . In Fig. 2, we plot the global order parameter for different values of the inter-star coupling parameter ( $\varepsilon$ ). For any given value of $\varepsilon$ , we observe that the system shows a first-order transition to synchrony. As shown in Fig. 2, the order parameter changes discontinuously for both the forward and backward continuations as $\lambda$ varies. Since the backward and forward transition points are distinct, the system shows a well-defined hysteresis. The hysteresis width increases as the coupling strength $\varepsilon$ varies. As the inter-star coupling is switched on, we observe a transition in the order parameter with a non-zero hysteresis width. We plot one such value in Fig. 2 for $\varepsilon=0.01$ , for which the hysteresis width is very small. If $\varepsilon$ is increased further ( $\varepsilon=0.1,0.3,0.5$ ), the hysteresis width increases as shown in Fig. 2. Note that the increase in hysteresis width occurs as a result of the shift in the forward transition point while the backward transition point is fixed for different values of $\varepsilon$ .

II.2 Analytical Results

To understand the dynamics, we make use of the Watanabe-Strogatz (WS) approach Watanabe-1993 ; Watanabe-1994 by introducing phase differences $\theta_{j},\theta^{\prime}_{j}$ and $\theta_{h}$ as follows

\displaystyle\theta_{j}=\phi_{j}-\phi_{h},\quad\theta^{\prime}_{j}=\phi^{\prime}_{j}-\phi^{\prime}_{h},\quad\theta_{h}=\phi^{\prime}_{h}-\phi_{h}.

(3)

In terms of new variables (Eq. (3)) along with $N_{1}=N_{2}=N^{\prime}$ , the system Eqs. (1) can be written as,

$\displaystyle\dot{\theta}_{j}$	$\displaystyle=\omega(1-\beta)-\lambda\sin{\theta_{j}}-\frac{\beta\lambda}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta_{l}}-\beta\varepsilon\sin{\theta_{h}},$	(4)
$\displaystyle\dot{\theta}^{\prime}_{j}$	$\displaystyle=\omega(1-\beta)-\lambda\sin{\theta^{\prime}_{j}}-\frac{\beta\lambda}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta^{\prime}_{l}}+\beta\varepsilon\sin{\theta_{h}},$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\frac{\beta\lambda}{N^{\prime}}\left[\sum_{l=1}^{N^{\prime}}\sin{\theta^{\prime}_{l}}-\sum_{l=1}^{N^{\prime}}\sin{\theta_{l}}\right]-2\beta\varepsilon\sin{\theta_{h}}$

Note that the terms $\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta_{l}}$ and $\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta^{\prime}_{l}}$ can be written as

	$\displaystyle\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta_{l}}=\operatorname{Im}\left[\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}e^{\mathrm{i}\theta_{l}}\right]=\operatorname{Im}(R_{1}),$		(5)
	$\displaystyle\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}\sin{\theta^{\prime}_{l}}=\operatorname{Im}\left[\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}e^{\mathrm{i}\theta^{\prime}_{l}}\right]=\operatorname{Im}(R_{2}).$		(5)

where $R_{1}$ and $R_{2}$ are the local order parameters. Substituting Eq. (5) in Eq. (4) with $\Delta\omega=\omega(1-\beta)$ , we get

$\displaystyle\dot{\theta}_{j}$	$\displaystyle=\Delta\omega+\lambda\operatorname{Im}(e^{-\mathrm{i}\theta_{j}})-\beta\lambda\operatorname{Im}(R_{1})-\beta\varepsilon\sin{\theta_{h}},$	(6)
$\displaystyle\dot{\theta}^{\prime}_{j}$	$\displaystyle=\Delta\omega+\lambda\operatorname{Im}(e^{-\mathrm{i}\theta^{\prime}_{j}})-\beta\lambda\operatorname{Im}(R_{2})+\beta\varepsilon\sin{\theta_{h}},$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\operatorname{Im}(R_{2})-\operatorname{Im}(R_{1})\right]-2\beta\varepsilon\sin{\theta_{h}}$

This higher dimensional equation can be reduced to a lower dimensional equation using WS ansatz Watanabe-1993 ; Watanabe-1994 for the general form

\dot{\theta}_{k}=g(t)+\operatorname{Im}(G(t)e^{-\mathrm{i}\theta_{k}}).

(7)

By comparing Eq. (6) with Eq. (7), we get

	$\displaystyle g(t)$	$\displaystyle=\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(R_{1});\quad G(t)=\lambda,$		(8)
	$\displaystyle g^{\prime}(t)$	$\displaystyle=\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(R_{2});\quad G^{\prime}(t)=\lambda$		(8)

The WS transformation from the phase variables $\theta_{k}$ to a set of global variables $z,\alpha$ ( $z$ being a complex and $\alpha$ being a real variable) changes Eq. (7) to

	$\displaystyle\dot{z}$	$\displaystyle=\mathrm{i}g(t)z+\frac{G(t)}{2}-\frac{G(t)^{*}}{2}z^{2},$		(9)
	$\displaystyle\dot{\alpha}$	$\displaystyle=g(t)+\operatorname{Im}\left[z^{*}G(t)\right]$		(9)

With WS transformation applied for the two-star network, we are left with variables $z,\alpha,z^{\prime},\alpha^{\prime}$ and $\theta_{h}$ whose variation is given as

$\displaystyle\dot{z}$	$\displaystyle=\mathrm{i}\left[\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z)\right]z+\frac{\lambda}{2}\left(1-z^{2}\right),$	(10)
$\displaystyle\dot{\alpha}$	$\displaystyle=\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z)+\lambda\operatorname{Im}(z^{*}),$
$\displaystyle\dot{z^{\prime}}$	$\displaystyle=\mathrm{i}\left[\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z^{\prime})\right]z^{\prime}+\frac{\lambda}{2}\left(1-{z^{\prime}}^{2}\right),$
$\displaystyle\dot{\alpha}^{\prime}$	$\displaystyle=\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z^{\prime})+\lambda\operatorname{Im}({z^{\prime}}^{*}),$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\operatorname{Im}(z^{\prime})-\operatorname{Im}(z)\right]-2\beta\varepsilon\sin{\theta_{h}}$

where $z,z^{\prime}$ are the modified order parameters corresponding to $R_{1},R_{2}$ respectively. Since the equations for $z,z^{\prime}$ do not contain any $\alpha,\alpha^{\prime}$ terms, we can solve them independently along with the equation for $\theta_{h}$ as follows

$\displaystyle\dot{z}$	$\displaystyle=\mathrm{i}\left[\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z)\right]z+\frac{\lambda}{2}\left(1-z^{2}\right),$	(11)
$\displaystyle\dot{z^{\prime}}$	$\displaystyle=\mathrm{i}\left[\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\operatorname{Im}(z^{\prime})\right]z^{\prime}+\frac{\lambda}{2}\left(1-{z^{\prime}}^{2}\right),$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\operatorname{Im}(z^{\prime})-\operatorname{Im}(z)\right]-2\beta\varepsilon\sin{\theta_{h}}$

The complex variables $z$ and $z^{\prime}$ can be written as,

	$\displaystyle z=\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}e^{\mathrm{i}\theta_{l}}=\rho_{1}e^{\mathrm{i}\psi_{1}},$		(12)
	$\displaystyle z^{\prime}=\frac{1}{N^{\prime}}\sum_{l=1}^{N^{\prime}}e^{\mathrm{i}\theta^{\prime}_{l}}=\rho_{2}e^{\mathrm{i}\psi_{2}}$		(12)

Substituting Eq. (12) in Eq. (11), we get

$\displaystyle\dot{\rho}_{1}$	$\displaystyle=\frac{\lambda}{2}\left(1-\rho_{1}^{2}\right)\cos{\psi_{1}},$	(13)
$\displaystyle\dot{\psi}_{1}$	$\displaystyle=\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\rho_{1}\sin{\psi_{1}}-\frac{\lambda}{2\rho_{1}}\left(1+\rho_{1}^{2}\right)\sin{\psi_{1}},$
$\displaystyle\dot{\rho_{2}}$	$\displaystyle=\frac{\lambda}{2}\left(1-\rho_{2}^{2}\right)\cos{\psi_{2}},$
$\displaystyle\dot{\psi}_{2}$	$\displaystyle=\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\rho_{2}\sin{\psi_{2}}-\frac{\lambda}{2\rho_{2}}\left(1+\rho_{2}^{2}\right)\sin{\psi_{2}},$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\rho_{2}\sin{\psi_{2}}-\rho_{1}\sin{\psi_{1}}\right]-2\beta\varepsilon\sin{\theta_{h}}$

For steady state,

$\displaystyle\dot{\rho}_{1}$	$\displaystyle=\dot{\rho}_{2}=0,$	(14)
$\displaystyle\implies\rho_{1}$	$\displaystyle=\rho_{2}=1\quad\qquad\text{(synchronous branch),}$
$\displaystyle\text{or}\quad\cos{\psi_{1}}$	$\displaystyle=\cos{\psi_{2}}=0\quad\text{(asynchronous branch).}$

We shall now consider the cases to find the solutions for synchronous and asynchronous branches separately.

II.2.1 Synchronous branch

The stability of the synchronous branch gives the backward transition point of the explosive synchronization, so for $\rho_{1}=\rho_{2}=1$ , Eq. (13) becomes,

$\displaystyle\dot{\psi}_{1}$	$\displaystyle=\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\sin{\psi_{1}}-\lambda\sin{\psi_{1}},$	(15)
$\displaystyle\dot{\psi}_{2}$	$\displaystyle=\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\sin{\psi_{2}}-\lambda\sin{\psi_{2}},$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\sin{\psi_{2}}-\sin{\psi_{1}}\right]-2\beta\varepsilon\sin{\theta_{h}}$

Since for the synchronized state, all the hubs are in the same state (i.e. $\theta_{h}=0$ ) and putting $\dot{\psi}_{1}=0$ , we obtain

	$\displaystyle\sin{\psi_{1}}=\frac{\Delta\omega}{\lambda(\beta+1)},$		(16)
	$\displaystyle\left\|\frac{-\omega(\beta-1)}{\lambda(\beta+1)}\right\|\leq 1.$		(16)

Thus, the synchronous branch is stable for

	$\displaystyle\lambda\geq\frac{\omega(\beta-1)}{(\beta+1)}.$		(17)
	$\displaystyle\implies\lambda_{c}^{b}=\frac{\omega(\beta-1)}{(\beta+1)}$		(17)

From Eq. (17), we note that the backward transition point is independent of the inter-star coupling parameter ( $\varepsilon$ ). For $\omega=1$ and $\beta=10$ , this transition point is found to occur at $\lambda^{b}_{c}=0.8$ , as shown by the black dashed line in Fig. 2. This is found to be in good agreement with the numerical results (Fig. 2).

II.2.2 Asynchronous branch

The stability of the asynchronous branch gives the forward point of the explosive transition, so we put $\psi_{1}=\psi_{2}=\frac{\pi}{2}$ in Eq. (13),

$\displaystyle\dot{\psi}_{1}$	$\displaystyle=\Delta\omega-\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\rho_{1}-\frac{\lambda}{2\rho_{1}}\left(1+\rho_{1}^{2}\right),$	(18)
$\displaystyle\dot{\psi}_{2}$	$\displaystyle=\Delta\omega+\beta\varepsilon\sin{\theta_{h}}-\beta\lambda\rho_{2}-\frac{\lambda}{2\rho_{2}}\left(1+\rho_{2}^{2}\right),$
$\displaystyle\dot{\theta}_{h}$	$\displaystyle=\beta\lambda\left[\rho_{2}-\rho_{1}\right]-2\beta\varepsilon\sin{\theta_{h}}$

Put $\dot{\psi}_{1}=\dot{\psi}_{2}=0$ and calculate the value of $\rho_{1}$ from Eq. (18) as,

\rho_{1}^{2}(2\beta+1)\lambda+2\left[\omega(\beta-1)+\beta\varepsilon\sin{\theta_{h}}\right]\rho_{1}+\lambda=0

(19)

Now, $\rho_{1}$ has to be real, so the determinant should be positive.

	$\displaystyle\left[\omega(\beta-1)+\beta\varepsilon\sin{\theta_{h}}\right]^{2}-\lambda^{2}(2\beta+1)\geq 0,$
	$\displaystyle\sin{\theta_{h}}\geq\frac{\lambda\sqrt{2\beta+1}-\omega(\beta-1)}{\beta\varepsilon},$
	$\displaystyle 1\geq\sin{\theta_{h}}\geq\frac{\lambda\sqrt{2\beta+1}-\omega(\beta-1)}{\beta\varepsilon},$
	$\displaystyle 1\geq\frac{\lambda\sqrt{2\beta+1}-\omega(\beta-1)}{\beta\varepsilon},$
	$\displaystyle\lambda\leq\frac{\beta\varepsilon+\omega(\beta-1)}{\sqrt{2\beta+1}}$
	$\displaystyle\implies\lambda_{c}^{f}=\frac{\beta\varepsilon\pm\omega(\beta-1)}{\sqrt{2\beta+1}}$		(20)

This gives the upper bound for the forward transition point. The forward transition point depends on the coupling parameter $\varepsilon$ . Numerically, the forward transition point for each $\varepsilon$ is found to be below $\lambda_{c}^{f}$ given by Eq. (II.2.2) (c.f. Fig. 2).

III N-coupled Star Network (N $>$ 2)

In this section, we extend our analysis to a more general scenario. We consider a system of $N$ coupled star networks. The connection topology is that of a ring, where we place the hub of a given star network on each node of a ring. The ring has $N$ such nodes. Thus, the hub of each star network is connected to the hubs of the two other star networks. The schematic diagram in Fig. 3 describes a system of $10$ star networks interacting on a ring topology. The equations of motion of $N$ such coupled star networks are given by

$\displaystyle\dot{\phi}_{j}^{k}$	$\displaystyle=$	$\displaystyle\omega+\lambda_{j}\sin(\phi_{h}^{k}-\phi_{j}^{k}),\hskip 28.45274pt1\leq j\leq N_{j}$
$\displaystyle\frac{1}{\beta}\dot{\phi}_{h}^{k}$	$\displaystyle=$	$\displaystyle\omega+\frac{\lambda_{j}}{N_{j}}\sum_{j=1}^{N_{j}}\sin(\phi_{j}^{k}-\phi_{h}^{k})$	(21)
		$\displaystyle+\frac{\varepsilon}{2}\sum_{i=1}^{N}\left\{A(i,k)\sin(\phi_{h}^{i}-\phi_{h}^{k})\right\}$	(21)

where the superscripts $k,i=1,2\dots N~{}(k\neq i)$ represent each star network and the subscript $j=1,2\dots N_{j}$ represents the nodes of an individual star network. Thus, ${\phi_{j}^{k}}$ is the phase of the $j^{th}$ oscillator in the $k^{th}$ star network and ${\phi_{h}^{k}}$ is the phase of the hub in the $k^{th}$ star network. $A(i,k)$ is the adjacency matrix that describes the coupling between the hubs of the star networks $i$ and $k$ . $A(i,k)=1$ when there is a connection between the hub of the $i^{th}$ star and the hub of the $k^{th}$ star, otherwise it is equal to zero. Now, one can generalize the analysis done in Sec. II.2, by introducing the phase differences given by

	$\displaystyle\theta_{j}^{k}$	$\displaystyle=$	$\displaystyle\phi_{j}^{k}-\phi_{h}^{k},$
	$\displaystyle\theta_{h}^{k}$	$\displaystyle=$	$\displaystyle\phi_{h}^{k+1}-\phi_{h}^{k},$		(22)

For the star network $k$ , the order parameter $R_{k}$ may be defined as (c.f. Eq. (5))

\displaystyle\frac{1}{N_{k}}\sum_{l=1}^{N_{k}}\sin{\theta_{l}^{k}}=\operatorname{Im}\left[\frac{1}{N_{k}}\sum_{l=1}^{N_{k}}e^{\mathrm{i}\theta_{l}^{k}}\right]=\operatorname{Im}(R_{k})

(23)

Rewriting Eq. (III) in terms of new phase variables (Eq. (III)) and the order parameter $R_{k}$ (Eq. (23)), we obtain

$\displaystyle\dot{\theta}^{k}_{j}$	$\displaystyle=$	$\displaystyle\Delta\omega+\lambda\operatorname{Im}(e^{-\mathrm{i}\theta_{j}^{k}})-\beta\lambda\operatorname{Im}(R_{k})$
		$\displaystyle-\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k}}-\sin{\theta_{h}^{k-1}}\right],$
$\displaystyle\dot{\theta}^{k}_{h}$	$\displaystyle=$	$\displaystyle\beta\lambda\left[\operatorname{Im}(R_{k+1})-\operatorname{Im}(R_{k})\right]$
		$\displaystyle+\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k+1}}+\sin{\theta_{h}^{k-1}}-2\sin{\theta_{h}^{k}}\right],$

Comparing Eq. (III) with the standard equation

\dot{\theta_{j}^{k}}=g_{k}(t)+\operatorname{Im}\left[G_{k}(t)e^{-\mathrm{i}\theta_{j}^{k}}\right],

(25)

we get,

	$\displaystyle g_{k}$	$\displaystyle=$	$\displaystyle\Delta\omega-\beta\lambda\operatorname{Im}(R_{k})-\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k}}-\sin{\theta_{h}^{k-1}}\right],$
	$\displaystyle G_{k}$	$\displaystyle=$	$\displaystyle\lambda$

The variation of the complex global variables ( $z_{k}$ ) can be written as

\displaystyle\dot{z_{k}}=\mathrm{i}g_{k}G_{k}+\frac{G_{k}}{2}-\frac{G_{k}^{*}z_{k}^{2}}{2};\quad k=1,2...N

(27)

Put

z_{k}=\rho_{k}e^{\mathrm{i}\psi_{k}}

(28)

In terms of $\rho_{k}$ and $\psi_{k}$ , Eqs. (III) are given by

$\displaystyle\dot{\rho}_{k}$	$\displaystyle=$	$\displaystyle\frac{\lambda}{2}\left(1-\rho_{k}^{2}\right)\cos{\psi_{k}},$
$\displaystyle\dot{\psi}_{k}$	$\displaystyle=$	$\displaystyle\Delta\omega-\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k+1}}-\sin{\theta_{h}^{k}}\right]-\beta\lambda\rho_{k}\sin{\psi_{k}}$
		$\displaystyle-\frac{\lambda}{2\rho_{k}}\left(1+\rho_{k}^{2}\right)\sin{\psi_{k}},$
$\displaystyle\dot{\theta}^{k}_{h}$	$\displaystyle=$	$\displaystyle\beta\lambda\left[\rho_{k+1}\sin{\psi_{k+1}}-\rho_{k}\sin{\psi_{k}}\right]$	(29)
		$\displaystyle+\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k+1}}+\sin{\theta_{h}^{k-1}}-2\sin{\theta_{h}^{k}}\right]$	(29)

For steady state,

\displaystyle\dot{\rho}_{k}=\dot{\psi}_{k}=0

(30)

From Eq. (III), the solution $\rho_{k}=1$ corresponding to $\dot{\rho}_{k}=0$ accounts for the stability of the synchronized state, and hence gives the backward transition point. Put $\rho_{k}=1$ in Eq. (III) and solve for $\psi_{k}$ as follows

\displaystyle\Delta\omega-\frac{\beta\varepsilon}{2}\left[\sin{\theta_{h}^{k+1}}-\sin{\theta_{h}^{k}}\right]-\beta\lambda\sin{\psi_{k}}-\lambda\sin{\psi_{k}}=0

In the synchronized state, all the hubs are in the same state, so the phase difference between any two hubs:

\theta_{h}^{k+1}=\theta_{h}^{k}=0

(32)

Substituting Eq. (32) in Eq. (III),

	$\displaystyle\sin{\psi_{k}}=\frac{-\omega(\beta-1)}{\lambda(\beta+1)},$
	$\displaystyle\left\lvert\frac{-\omega(\beta-1)}{\lambda(\beta+1)}\right\rvert\leq 1,$
	$\displaystyle\lambda_{c}^{b}\geq\frac{\omega(\beta-1)}{(\beta+1)}.$		(33)

The minimum coupling strength for which the synchronized branch remains stable can be obtained by using the analysis done for a two-coupled star network. The synchronous branch is stable after $\lambda=0.8$ and is independent of the inter-star coupling ( $\epsilon$ ) strength and the number of stars ( $N$ ) coupled together.

For the asynchronous branch, we again put $\dot{\psi_{k}}=0$ (Eq. (30)) with $\psi_{k}=\pi/2$ in Eq. (III), we get

\rho_{k}^{2}(2\beta+1)\lambda+\left[2\omega(\beta-1)+\beta\varepsilon\left(\sin{\theta_{h}^{k+1}}-\sin{\theta_{h}^{k}}\right)\right]\rho_{k}+\lambda=0

(34)

The determinant of the roots of the above quadratic equation should be positive, and it will lead to

\displaystyle\sin{\theta_{h}^{k+1}}-\sin{\theta_{h}^{k}}\geq\frac{2\lambda\sqrt{2\beta+1}-2\omega(\beta-1)}{\beta\epsilon}

(35)

The difference of the sine terms can be a maximum of $2$ . This way, the stability of the asynchronous branch is achieved in the case of the $N$ -coupled star networks. This gives an upper bound for the forward transition points.

Now we shall numerically explore a situation when there are more star networks coupled together $(N>2)$ and study the order parameter as a function of the intra-star coupling strength for a fixed inter-star coupling. We consider a network of three coupled stars $(N=3)$ , for which the equations of motion are given by Eq. (III) by setting $N=3$ . The variation of the global order parameter $R$ as a function of $\lambda$ is shown in Fig. 4. We observe that the order parameter changes abruptly as the intra-star coupling ( $\lambda$ ) is varied in the forward and the backward directions. Since the forward and backward transition points are distinct, the system exhibits a well-defined hysteresis. We also note that as the value of $\varepsilon$ increases, the hysteresis width increases. At $\varepsilon=0.02$ , the hysteresis width is very small, which increases as the value of $\varepsilon$ is increased, as shown by the red line ( $\varepsilon=0.4$ ), green line ( $\varepsilon=1$ ) and the blue line ( $\varepsilon=2.4$ ). Moreover, the backward transition point is exactly the same as predicted by Eq. (33). This is shown by a dashed line in Fig. 4.

Further, we repeated the analysis for a bigger network taking $N=10$ (Eq. (III)) and found that the results are qualitatively the same as shown in Fig. 5. Again, we note that the forward and the backward transition points are distinct, resulting in well-defined hysteresis. We show this by plotting the order parameter for $\varepsilon=0.02,0.3,0.4,0.5$ as shown in Fig. 5. The backward transition is the same and matches well with the analytically predicted point given by Eq. (33).

IV Discussion

In this work, we study the collective dynamics of multiple star networks that are coupled through their hubs. We observe that the system exhibits a global first-order transition as the intra-star coupling ( $\lambda$ ) increases in the forward and backward directions for different values of the inter-star coupling ( $\varepsilon$ ). As we increase $\varepsilon$ , the hysteresis width increases. The forward transition point (FTP) shifts towards the higher $\lambda$ value while the backward transition point (BTP) remains the same, thereby changing the hysteresis width. The BTP is found to be independent of the number of stars coupled together and the strength of the inter-star coupling. We have shown results for two, three and ten-coupled stars. These findings have been verified through the analytical treatment of the system. The dependency of the forward transition point on the inter-star coupling strength is obtained by the WS analysis. We analytically calculate the upper bound of the forward transition point and note that the numerically observed transition points are below the predicted value. In the case of BTP, the analytically predicted value is found to be in good agreement with the numerical findings.

Our work highlights the role of topological connectivity in achieving synchronization in networks of heterogeneous degree distributions, namely star networks. Such network configurations have been known to form the backbone of several real world complex systems Bonifazi-2009 ; Bertolero-2018 . Unravelling the relationship between network architecture and emergent collective dynamics has been the focus of several studies found in the literature Singh-2011 . This work, in particular, elucidates the role of coupling via hubs in observing global synchrony. Similar studies have been performed on networks of chemical oscillators, where the authors have investigated the modes and mechanisms of achieving synchronization in systems coupled via hubs Sawicki-2021 . A system with interacting hubs, namely the central elements, appears in various fields, namely the communication systems, social networks and mammalian brain Zhou-2006 ; Gardenes-2010 ; Campbell-1994 ; Kim-2008 . The importance of hubs in relaying information and optimising network operation costs has also been discussed in several works Aykin-1995 ; Sawicki-2021 ; Campbell-1994 ; Kim-2008 . We, therefore, believe our results have a many potential practical applications which could be a subject of future investigations.

Further, in the current study, we consider the Kuramoto model, which is analytically tractable and allows for fast numerical simulations. However, in more realistic situations, it will be interesting to consider more complicated dynamics to explore other dynamical features such as cluster synchronization, chimera states, relay synchronization etc.

ACKNOWLEDGMENTS

RV wants to acknowledge CSIR, India, for the SRF fellowship under file no. 09/112(0601)/2018-EMR-I.

References

(1) A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Science (Cambridge University Press, Cambridge, 2001).
(2) S. Boccaletti, A. N. Pisarchik, C. I. D. Genio, and A. Amann, Synchronization: From Coupled Systems to Complex Networks (Cambridge University Press, Cambridge, 2018).
(3) F. Dorfler, and F. Bullo, SIAM J. Control Optim. 50, 3 (2012).
(4) A. E. Motter, S. A. Myers, M. Anghel, and T.Nishikawa, Nat. Phy. 9, 191 (2013).
(5) S. R. Cobb, E. H. Buhl, K. Halasy, O. Paulsen, and P. Somogyi, Nat. 378, 75 (1995).
(6) F. Ling, Synchronization in Digital Communication Systems (Cambridge University Press, Cambridge, 2017).
(7) A. T. Winfree, When Time Breaks Down (Princeton Univ Press, Princeton, 1987).
(8) D. J. Watts, and S. H. Strogatz, Nat. (London) 393, 440 (1998).
(9) A. L. Barabási, and R. Albert, Science 286, 509 (1999).
(10) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D. U. Hwang, Phys. Rep. 424, 175 (2006).
(11) Y. Moreno, and A. F. Pacheco, Europhys. Lett. 68, 603 (2004).
(12) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93 (2008).
(13) C. Zhou, L. Zemanová, G. Zamora, C. C. Hilgetag, and J. Kurths, Phys. Rev. Lett. 97, 238103 (2006).
(14) J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, Phys. Rev. Lett. 106, 128701 (2011).
(15) C. Meena, K. Murali and S. Sinha, Int J. Bifurcation and Chaos 26, 1630023 (2016).
(16) S. S. Muni, and A. Provata, Nonlinear Dyn. 101, 2509 (2020).
(17) N. E. Kouvaris, M. Sebek, A. Iribarne, A. Díaz-Guilera, and I. Z. Kiss, Phys. Rev. E 95, 042203 (2017).
(18) J. Shena, J. Hizanidis, N. E. Kouvaris, and G. P. Tsironis, Phys. Rev. A 98, 053817 (2018).
(19) V. Vlasov, Y. Zou, and T. Pereira, Phys. Rev. E 92, 012904 (2015).
(20) Y. A. Muthanna, and H. H. Jafri, Phys. Rev. E 109, 054206 (2024).
(21) B. C. Coutinho, A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes, Phys. Rev. E 87, 032106 (2013).
(22) Y. Zou, T. Pereira, M. Small, Z. Liu, and J. Kurths, Phys. Rev. Lett. 112, 114102 (2014).
(23) P. V. Kuptsov, and A. V. Kuptsova, Phys. Rev. E 92, 042912 (2015).
(24) C. Song, S. Havlin, and H. A. Makse, Nat. 433, 392 (2005).
(25) J. Gómez-Gardeñes, G. Zamora-López, Y. Moreno, and A. Arenas, PLoS ONE 8(5), e12313 (2010).
(26) J. F. Campbell, Eur. J. Oper. Res. 72, 387 (1994).
(27) H. Kim, and M. E. O’Kelly, Geographical Anal., 41, 283 (2009).
(28) C. Xu, Y. Sun, J. Gao, W. Jia, Z. Zheng, Nonlinear Dyn. 94, 1267 (2018).
(29) J. Zhu, D. Huang, Z. Yu, and P. Pei, Front. Phys. 9, 782607 (2021).
(30) S. Boccaletti, G. Bianconi, R. Criado, C. I. del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wang, M. Zaninu, Phys. Rep. 544, 1 (2014).
(31) I. Leyva, I. Sendiña-Nadal, R. Sevilla-Escoboza, V. P. Vera-Avila, P. Chholak, and S. Boccaletti, Sci. Rep. 8, 8629 (2018).
(32) V. Rathore, A. Suman, and S. Jalan, Chaos 33, 091105 (2023).
(33) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev. Lett. 114, 038701 (2015).
(34) S. Watanabe and S. H. Strogatz, Phys. Rev. Lett. 70, 2391 (1993).
(35) S. Watanabe and S. H. Strogatz, Physica D 74, 197 (1994).
(36) P. Bonifazi, M. Goldin, M. A. Picardo, I. Jorquera, A. Cattani, G. Bianconi, A. Represa, Y. Ben-Ari, and R. Cossart, Science 326, 1419 (2009).
(37) M. A. Bertolero, B. T. T. Yeo, D. S. Bassett, and M. D’Esposito, Nat. Human Behav. 2, 765 (2018).
(38) T. U. Singh, K. Manchanda, R. Ramaswamy, and A. Bose, SIAM J. Appl. Dyn. Syst. 10, 987 (2011).
(39) J. Sawicki, J. M. Loulen, and E. Schöll, Chaos 31, 073131 (2021).
(40) T. Aykin, Eur. J. Oper. Res. 83, 200 (1995).