Synchronization of interconnected heterogeneous networks: The role of network sizes

Increasing evidence shows that real networks interact with each other, forming a network of networks (NONs). Synchronization, a ubiquitous process in natural and engineering systems, has fascinatingly gained rising attentions in the context of NONs. Despite efforts to study the synchronization of NONs, it is still a challenge to understand how do the network sizes affect the synchronization and its phase diagram of NONs coupled with nonlinear dynamics. Here, we model such NONs as star-like motifs to analytically derive the critical values of both the internal and the external coupling strengths, at which a phase transition from synchronization to incoherence occurs. Our results show that the critical values strongly depend on the network sizes. Reducing the difference between network sizes will enhance the synchronization of the whole system, which indicates the irrationality of previous studies that assume the network sizes to be the same. The optimal connection strategy also changes as the network sizes change, a discovery contradicting to the previous conclusion that connecting the high-degree nodes of each network is always the most effective strategy to achieve synchronization unchangeably. This finding emphasizes the crucial role of network sizes which has been neglected in the previous studies and could contribute to the design of a global synchronized system.


S1 Dynamics in the integrated grid
We consider an integrated grid composed of a traditional power grid and a microgrid, with N I and N II nodes respectively. The dynamics governing the integrated grid can be written aṡ where θ ∈ [0, 2π) denotes the phase of a node as a time dependent activity, I (II) means the node belongs to Network I (Network II) with i = 1, ..., N I (m = 1, ..., N II ), ω stands for the natural frequency, λ captures the interaction between nodes that have one link in the same network as internal coupling strength, and β encapsulates the interaction between nodes connecting the two networks as external coupling strength. The grid frequency ω is a fixed value, 50Hz or 60Hz, and ω is proportional to its degree k 1 , i.e. ω I i = k I i ω, ω II m = k II m ω. For internal coupling strength, λ I i j = V I i V I j Q red i, j = λ I ji in the traditional power grid 2 and λ II mn = V II m V II n |Y mn | = λ II nm in the microgrid with the frequency-droop controller 3 , where V is the voltage, Y accounts for the admittance and Q red refers to the reduced admittance matrix. As for external coupling strength, in the integrated grid. Here, we study the uniform networks with the dynamical equationṡ (S2) The application of the result derived by the uniform networks (S2) can be well applied to the non-uniform networks (S1), seen in the manuscript Fig. 4i.

S2 One-link connection
Before the construction of a microgrid, the experts need to determine how the microgrid will be integrated into the traditional power grid. We suggest that the topology between the microgrid and the traditional power grid is the factor that needs to be considered. The topology between them is called the connection strategy, which is found to be very important according to the previous research 4,5 . The connection strategy for two networks coupled through one link has three cases: H-H, H-L and L-L 4 .

S2.1 Case H-H
When the system is synchronized, the two networks are in the state of locking manifolds 6 . The locking manifolds in H-H case can be expressed as where a and b are the constants relating with the phase difference between the hub and the leaf. Since all leaves in the star motif are exactly the same, so a and b are the scalar 6 .Here the node N I is the hub in Network I while the node N II is the hub in Network II, same for the other cases. Using locking manifolds, (S2) for the synchronized system here can be rewritten aṡ Since the phase difference between θ hub and θ leaf does not change with time, we haveθ I Since θ I N I − θ II N II is a constant, θ I N I − θ I j and θ I N II − θ II j are constants too. Thus, (S4) yields As sin a and sin θ I N I − θ II N II in the l.h.s of (S6) and (S7) are limited in [−1, 1], the r.h.s of (S6) and (S7) are limited in [−1, 1] too. Hence, the range of λ and β can be obtained as namely,

S2.2 Case H-L
When the two networks are coupled through H-L connection, there exits three locking manifolds for the synchronized system: where the connectors coupling two networks are the hub node N I and the leaf node 1 II , a, b and c are the constants relating with the phase difference between the hub and the leaf. With the locking manifolds (S12), the dynamics of the synchronized system iṡ (S13) According to the locking manifolds (S12),θ I N I −θ I i = 0,θ II N II −θ II m = 0 andθ II N I −θ II 1 = 0. Thus, from (S13) we haveθ I i =θ I N I =θ II 1 =θ II m =θ II N II and The ranges of the sine values in (S14-S16) determine the ranges of the coupling strengths leading to

S2.3 Case L-L
The locking manifolds for the synchronized system in L-L case can be expressed as where the leaf node 1 I and the leaf node 1 II are the connectors linking two networks, a, b, c and d are the constant phase differences. Substituting the locking manifolds (S22) into (S2) for L-L case, we have Thus, the critical coupling strengths are Through the comparison of λ c in the three cases above, we find that λ H-H when N I = N II , which is verified in S4. As for external coupling strength, β c for all the modes in one-link connection are the same. That means whether hub acts as the connector or not doesn't influence external coupling strength. In sum, hub plays a significant role in internal coupling strength, but has no effect on external coupling strength.
We ran (S2) with N I = 11, N II = 16, ω = 1 for case H-H (    Figure S3: Critical coupling strengths in L-L case. Using (S2) on two networks in case L-L with N I = 11, N II = 16, ω = 1, we get how the phase differences change with time. Phase differences include the one between the hub and the leaf (not the leaf connector, blue curves for Network I and red curves for Network II) and the one between two hubs (green curves). a, With λ = 1.22, β = 0.36, the phase differences change with time. Similar results are found in b (λ = 1.22, β = 0.38) and c (λ = 1.24, β = 0.36). d, Phase differences are invariant with λ = 1.24, β = 0.38.

S3 Two-links connection
In two-links connection, the connection strategies investigated for the integrated grid are H-H L-L, H-L L-H, L-L L-L and H-L L-L. To get the expression of the critical coupling strengths here, Taylor approximation is used. The use of it is based on the fact that both the phases of the converter and the generator need to be very close to the phase of the injection points, suggesting that all the phase differences are close to zero 7 .

S3.1 Case H-H L-L
The locking manifolds in case H-H L-L are as same as (S22) in case L-L. Using the locking manifolds we can write (S2) for the synchronized system aṡ where the hub node N I connects to the hub node N II and the leaf node 1 I connects to the leaf node 1 II . Since the phase differences are constants,θ I (S34) By the setting of x = sin a − sin c, (S34) can be rearranged as Obviously, if y 1 and y 2 in (S35) have intersection during the common domain of the two functions, (S34) has solution and the system can be synchronized. The common domain is determined by arcsin in (S35) since its domain is limited in [−1, 1]. Through the observation of simulation, we find that if y 1 y 2 at the left endpoint of the domain and y 1 y 2 at the right endpoint of the domain, y 1 and y 2 have intersection. Without loss of generation, we set N I N II . And we record the left endpoint as x 1 , while the right endpoint as x 2 , same for case H-L L-H and case H-L L-L. In sum, there are two conditions for the synchronization: The endpoints are based on the values of the coupling strengths and node numbers, seen as Table S1. Since the endpoints are divided into three groups, there are three groups of discussions as below.
With the left endpoint N I +N II ω, after the linearization of y 1 (x 1 ) y 2 (x 1 ), we have which yields As for the right endpoint Thus, the second condition y 1 (x 2 ) y 2 (x 2 ) is satisfied without the limitation on the coupling strengths.
Here, since the left endpoint is as same as above, we get the same result as (S38) and only need to deal with the right endpoint. For the right endpoint The second condition y 1 When λ is very large, the l.h.s. of (S41) equals to 0. Thus, the r.h.s. 0, which leads to When λ is not very large, (S41) needs to be linearized and is reduced as Since (S43) is a quadratic equation with 2 λ > 0 and 0, we obtain the constrain on β as Here, the right endpoint is the same as above and we can get the same result as (S44). For the left endpoint x 1 = − β λ , we have there is no limit on β for the second condition y 1 (x 1 ) y 2 (x 1 ). Unlike the case in one-link connection, the two critical coupling strengths are inter-dependence. For union, we set the critical external coupling strength as a function of internal coupling strength, i.e. β c = f (λ). In sum, we have
The endpoints of the common domain of y 1 and y 2 can be seen in Table. S2. According to the different values of endpoints, there are three groups of discussions in the following part. Without loss of generation, we set N I II .
For the left endpoint Since ) ω , which means the first condition y 1 (x 1 ) y 2 (x 1 ) is satisfied without limitation.
For the right endpoint x 2 = 1, we have N I +N II ω, the linearization of y 1 (1) y 2 (1) from (S63) is needed and leads to Namely, (S66) The first condition (S68) The linearization of (S68) yields Since (S69) is a quadratic equation with 2 λ > 0 and 0, we have The requirement brought by the second condition here is as same as (S65).
3) β < λ − N I +N II −4 N I +N II ω Owing to the same left endpoint, the requirement brought by the first condition y 1 (x 1 ) y 2 (x 1 ) is the same as (S70). As for the right endpoint x 2 = β λ + N I +N II −4 λ(N I +N II ) ω, we have ) ω , indicating that the second condition y 1 (x 2 ) y 2 (x 2 ) is satisfied.
Similar to case H-H L-L, we set β c as a function of λ and get

S3.3 Case L-L L-L
The locking manifolds for two synchronized networks in L-L L-L case are Here, the leaf node 1 I connects to the leaf node 1 II and the leaf node 2 I connects to the leaf node 2 II . The dynamical equations for the synchronized system arė According to the locking manifolds, we getθ I 1 =θ I 2 =θ I 3 = ... =θ I N I =θ II 1 =θ II 2 =θ II 3 = ... =θ II N II straightly, accompanying with which (S76) leads to sin c + sin d = 2 sin a, (S80) sin e + sin f = 2 sin b, (S81) we apply Taylor approximation sin x = x + O(x) to (S78) with, and get With the application of Taylor approximation sin x = x + O(x) to (S78) -(S81), the conclusion can be drawn as c = e, d = f.
It is noticed that the conclusion c = θ I N I − θ I 1 = θ I N I − θ I 2 = e can be seen straightly if the initial phases differences of θ I N I − θ I 1 and θ I N I − θ I 2 are the same, since the dynamics of θ I N I − θ I 1 and θ I N I − θ I Thus, (S80) and (S81) lead to Since the sine value is limited in [−1, 1], the r.h.s. of (S77) and (S86-S89) result in Consequently, the critical coupling strengths are:

S3.4 Case H-L L-L
The locking manifolds for two synchronized networks in H-L L-L case are Here, the hub node N I connects to the leaf node 1 II and the leaf node 1 I connects to the leaf node θ II 2 . For the synchronized system, the dynamical equations (S2) can be rewritten aṡ Seen from (S97), we can draw the conclusionθ I straightly. The substitution of this conclusion back into (S97) yields Substituting (S100) into (S101), we get The combination of (S102), (S103) and (S104) yields , we obtain the following equation with (S99) and (S100): Using (S103), (S104) and (S106) jointly, the transcendental equations become sin f = 2 sin a − sin c, The transcendental equations can be separated into two parts: (S112) The system can be synchronized when (S107) has a solution, which means y 1 and y 2 in (S112) have intersection during the common domain. Through the observation of simulation, there are two conditions for the intersection: 1. y 1 (x 1 ) y 2 (x 1 ); 2. y 1 (x 2 ) y 2 (x 2 ).
The endpoints of the common domain of y 1 and y 2 are based on the values of the coupling strengths and node numbers, seen in Table. S3. There are three groups of discussions: parameter endpoint left endpoint x 1 right endpoint x 2 For the left endpoint (S113) The linearization of the first condition y 1 (x 1 ) y 2 (x 1 ) is (S114) Since λ N I +N II −4 N I +N II and N I N II , the r.h.s. of (S114) is larger than zero. Thus, if the l.h.s. of (S114) is less than zero, (S114) is satisfied. Setting the l.h.s. of (S114) as zero, we have Hence, when λ < 2 (S116) As for the right endpoint x 2 = 1, (S112) becomes (S117) To find out the requirement that satisfies the second condition y 1 (1) y 2 (1), the linearization of y 1 (1) y 2 (1) is needed and provides (S118) It can be expressed as If the left side is less than 0 while the right side is larger than 0, than (S118) can never be found.
(S128) 4) N I < N II , β λ − N I +N II −4 N I +N II ω Since the left endpoint is as same as 1), the first condition y 1 (x 1 ) y 2 (x 1 ) also leads to (S114). It is noticed that 2 ω, there is no constraint on β for the first condition. When 2(N II −2) ω, (S116) is needed to satisfy the first condition. And λ < 2(N II −2) N I +N II ω is not permitted here, seen in (S113). Thus For the right endpoint which requires λ 2(N II −2) ) ω , indicating that the second condition is satisfied.
The left endpoint x 1 = − β λ + N I +N II −4 λ(N I +N II ) ω leads to (S131) (S132) Using (S129), the linearization of y 1 (x 1 ) y 2 (x 1 ) yields leading to For the right endpoint β λ + 3N I −N II −4 λ(N I +N II ) ω, (S112) becomes (S135) The linearization of the second condition y 1 which means In sum, For N I N II , , (S140) For N I < N II , Moreover, the maximal external coupling strength exits: • for λ < 2 ; (S142) • for λ < 4 (S143) In two-links connection, we derive the critical values for two kinds of coupling strengths, λ and β, the validations about which are shown in Fig. S4 and Fig. S5 respectively. In Fig. S4 and Fig. S5, there are 8 sub-figures that separate into two columns: 1) the first column includes a, c, e and g, describing the non-synchronized phase; 2) the second column includes b, d, f and h, capturing the synchronized state. There are four rows representing four cases respectively. For example, sub-figures a and b in the first row are for the same H-H L-L case.
In Fig. S4, all the plots are about the simulation with N I = 16, N II = 11, ω = 1, β = 10. But the first column is with λ = λ c − 0.01, while the second column is with λ = λ c + 0.01. N I +N II ω = 1.11, which is close to 1.10 found in the simulation.
In Fig. S5, we ran the simulation with N I = 16, N II = 11, ω = 1, λ = 10 for all plots. There are two columns: the first column is about β c − 0.01; the second column is about β c + 0.01. For all the connection modes in two-links connection, β c found in the simulations are 0.19, which are in consistent with the values derived from (S48), (S74), (S95) and (S137).
In Figs. S4 and S5, the values derived are very close to the values obtained from the simulations, which work as the verification.
In all connection modes, the values of sin a are the same. Since sin a dominates the order parameter, the synchronization levels for different connection modes are nearly the same. According to the result derived in S2 and S3, β c = 0 and λ c are the same for all connection modes when N I = N II . This special case is validated through the simulations using (S2) with N I = N II = 20, ω = 1, shown in Table. S4. It can be seen that for all the connection modes, even if β is as small as 0.01, the system can still be synchronized when λ is larger than the critical value. This result indicates that two coupled networks seem to be uncoupled when their network sizes are equal.
• shape: in the synchronization area of case H-H, the approximate appearance of the boundary. For Groups 1 and 3, the shape of their synchronization areas is a equilateral triangle. For Group 2, the shape of its synchronization area is a skinny triangle.
In one hand, Group 1 and Group 2 share the same λ and size, while they have different β and shape.
In the other hand, Group 1 and Group 3 have different λ and size, while they share the same β and shape. Thus, we can draw a conclusion that λ controls the size of the synchronization area while β decides its shape.