Collective dynamics of random Janus oscillator networks

Janus oscillators have been recently introduced as a remarkably simple phase oscillator model that exhibits nontrivial dynamical patterns—such as chimeras, explosive transitions, and asymmetry-induced synchronization—that were once observed only in speciﬁcally tailored models. Here we study ensembles of Janus oscillators coupled on large homogeneous and heterogeneous networks. By virtue of the Ott-Antonsen reduction scheme, we ﬁnd that the rich dynamics of Janus oscillators persists in the thermodynamic limit of random regular, Erd˝os-Rényi, and scale-free random networks. We uncover for all these networks the coexistence between partially synchronized states and a multitude of solutions of a collective state we denominate as a breathing standing wave , which displays global oscillations. Furthermore, abrupt transitions of the global and local order parameters are observed for all topologies considered. Interestingly, only for scale-free networks, it is found that states displaying global oscillations vanish in the thermodynamic limit.


I. INTRODUCTION
Research on coupled oscillators in the past decade has been marked by the discovery of many intriguing patterns in the collective behavior of networks [1,2]. Notable examples of such patterns are chimeras [3], states in which populations of synchronous and asynchronous oscillators coexist; explosive synchronization transitions [2,4], which appear as a consequence of constraints in the natural frequency assignment; and asymmetry-induced synchronization [5], a state in which synchrony is counterintuitively favored by oscillator heterogeneity. In all these cases, phase oscillator models had to be specially designed so that those nontrivial states could be scrutinized. Very recently, however, Nicolaou et al. [6] defined an oscillator model coined as Janus oscillator. The name is inspired in the homonym two-faced god of Roman mythology and reflects the two-dimensional character of an isolated oscillator. Each "face" of a Janus unit consists of a Kuramoto oscillator, whose natural frequency has the same absolute value but opposite sign to the frequency of its counterface. When coupled on one-dimensional regular graphs, Janus oscillators have been found to exhibit a striking rich dynamical behavior that encompasses the co-occurrence of several dynamical patterns, in spite of the simplicity of the topology and the oscillator model itself [6].
Here we ask whether such a rich collective dynamics exists on more general networks made up of Janus oscillators. To address this issue, we employ the Ott-Antonsen (OA) ansatz [7,8] and obtain a reduced set of equations describing the system's evolution. From this reduced representation we find that, indeed, peculiar patterns of synchrony persist when Janus oscillators are placed on random regular, Erdős-Rényi (ER), and scale-free (SF) random networks. More specifically, we provide analytical and numerical evidence that the multitude of states in Janus dynamics on heterogeneous networks is a consequence of the coexistence of infinite neutrally stable limit-cycle trajectories, a state which we denominate "breathing standing-waves" and which has a distinct nature from the states reported in [6]. In fact, we find, unexpectedly, that heterogeneous connections do not yield chimeras and other nontrivial states that appear in rings of Janus oscillators [6]. In addition to breathing standing-waves, we report co-occurrence between classical partially synchronized and standing waves. We further show that for high average degrees the collective states of ER networks are accurately described by the reduced system obtained for random regular ones. Interestingly, we demonstrate that the coupling range in which global oscillations are possible vanishes in the thermodynamic limit of SF networks.

II. MODEL
We define the dynamics of N Janus oscillators [6] on heterogeneous networks aṡ where the 2N×2N matrix W is defined as Natural frequencies are assigned as ω i = ω 0 + /2 ≡ω 1 , for i = 1, . . . , N, and ω i = ω 0 − /2 ≡ω 2 , for i = N + 1, . . . , 2N, where is the frequency mismatch and ω 0 is the average frequency, which we set to ω 0 = 0. System (1) is analogous to a bipartite or multilayer network [9] in which oscillators belonging to the same group do not interact with one another, while connections between groups are encoded in matrix A (see Fig. 1 for illustration). Notice also that interactions between oscillators are weighted by coupling σ , except oscillators with indexes (i, i + N ), i ∈ [1, N]; these pairs of nodes interact with coupling β. We then refer to the bipartite network with bimodal frequency distribution as a "Janus oscillator network" 1 defined in Eqs. (1) and (2).

III. THEORY AND SIMULATIONS
By defining the local order parameters R i = 2N j=1 W i j e iθ j , Eqs. (1) are then decoupled aṡ Following [10], we consider an ensemble of systems defined by Eq. (1) with fixed coupling matrix W. In this formulation (see Appendix A), we describe the oscillators at time t by the marginal density ρ i (θ i , ω i , t ), which satisfies the continuity equation ∂ t ρ i + ∂ θ i (ρ iθi ) = 0. Let us suppose that at first frequencies ω i are distributed according to a function g(ω i ). By expanding ρ i in Fourier series and applying the OA ansatz to its coefficients [7,10] we have 1 Phase oscillators coupled on bipartite networks have been addressed previously; see, e.g., [14][15][16]. Of particular interest is the work by Samani and Ghanbarian [14], who investigated Kuramoto oscillators on fully connected bipartite networks in which the natural frequencies are bimodally distributed similarly as in the setup depicted in Fig. 1. However, despite the similarity of both models, the new breathing standing-waves and the routes to different states uncovered here remained unnoticed in [14].
Inserting the previous equation into the continuity equation of ρ i , we obtain the evolution for coefficientsα i : where the coefficients R i are calculated as Inserting the expansion of ρ j in the previous equation yields

A. Random regular networks
Let us now consider the case in which each subpopulation of the Janus coupling scheme is a random regular network with degree k. More precisely, each oscillator of the subpopulation that rotates with frequencyω 1 = /2 is randomly connected to k oscillators of subpopulation 2 (ω 2 = − /2), and vice versa. Since oscillators within each group are identical, i.e., g(ω i ) = δ(ω i − /2) for i = 1, . . . , N, and g(ω i ) = δ(ω i + /2) for i = N + 1, . . . , 2N, we then assume the solutionsα 1 = · · · =α N ≡ α 1 andα N+1 = · · · = α 2N ≡ α 2 . Hence, the local order parameters R i are reduced Inserting solutions α 1,2 and R 1,2 in Eq. (3), we obtain the reduced set of equations, which in terms of the coordinates α 1,2 = r 1,2 e iψ 1,2 are written asṙ where δ = ψ 1 − ψ 2 is the phase lag between subpopulations [full derivation of Eq. (5) is given in Appendix A]. The variables r 1 and r 2 turn out to be the order parameters measuring the level of synchronization within each subpopulation; the traditional Kuramoto order parameter evaluating the global synchrony is obtained through R(t ) = 1 2 |r 1 e iψ t (t ) + r 2 e iψ 2 (t ) |. States that emerge from system (5) are summarized as follows: (1) a partially synchronized state in which r 1,2 = 1, while the subpopulations remain separated by a constant δ (hence, R < 1); (2) a standing-wave (SW) state, where the bulks of the two fully synchronized populations (r 1,2 = 1) rotate in opposite directions, yielding an incessantly rotating δ; and (3) a distinct form of SW emerges when 0 < r 1,2 < 1: along with the increase or decrease in δ, the order parameters r 1,2 exhibit a breathing behavior, as depicted in the simulation shown Fig. 2. Henceforth we refer to this state as "breathing SW." Importantly, we emphasize that the breathing SW is not a form of chimera state and hence is a distinct collective state from the traveling chimeras reported in [6]. To uncover the conditions for the existence of the partially synchronized state, we set r 1 (t ) = r 2 (t ) ≡ r(t ) in Eqs. (5), leading toṙ = (1/2)(β + kσ )r(1 − r 2 ) cos δ andδ = − − (β + kσ )(1 + r 2 ) sin δ. Imposingṙ = 0, we notice that r = 1 is always a fixed point. Inserting the latter solution in the equation forδ = 0, we find that sin δ = − /2(β + kσ ). Thus, the partially synchronized state remains stable when − /2|β + kσ | 1 is satisfied. In terms of coupling σ , we then write the first critical conditions as which determines the coupling range where the partially synchronized state exists. The total order parameter R is then given by where the "-" branch is stable for σ σ c 1 , whereas the "+" branch is stable in the region σ σ c 2 . For σ c 1 < σ < σ c 2 , the limit cycle solution ofδ holds and SW states emerge with perfectly synchronized subpopulations (r 1,2 = 1). Linear stability analysis of the incoherent state (α 1,2 = 0) in Eqs. (4) reveals that the eigenvalues of the Jacobian matrix become purely imaginary at |β + kσ | = . Therefore, limit Order parameters 3. (a) Stability diagram of system (5) for k=40 and =1. "Partial sync" refers to the state in which r 1,2 = 1 and δ = 0. SW denotes parameter regions for each r 1,2 = 1 and for which the phase lag δ rotates with a nonzero frequency. The regions where a multitude of solutions with r 1,2 < 1 andδ = 0 are found are labeled as "breathing SW." (b) Order parameter for β = −0.3. Solid and dashed lines of R for σ ∈ (−∞, σ 2 ] ∪ [σ 1 , ∞) are obtained with Eq. (7). Line for r 1 = r 2 denotes the symmetric solution of Eq. (5). Dots are obtained by numerically evolving the original system (1) with N = 2500; each dot is an average over t ∈ [250, 500], with time step dt = 0.01. Gray lines for σ ∈ [σ c 4 , σ c 3 ] are generated numerically by evolving system (5) with random initial conditions. Gray dots depict the corresponding result yielded by the original system [Eq. (1)]. For the sake of clarity, we show only a small sample of the possible states attainable in the gray area. cycle solutions arise for Figure 3(a) outlines the critical conditions given by Eqs. (6) and (8) in the plane spanned by couplings β and σ . As it is seen, the partially synchronized state is favored by extreme values of both couplings, whereas states with oscillating synchrony appear for intermediate values in the parameter space. In Fig. 3(b) we visualize the evolution of the order parameters over a vertical section of the diagram in (a). Supposing we initiate the system with a negative σ in the "partial sync" region (σ < σ c 4 ), the order parameter R collapses in the solution given by Eq. (7). As σ is further increased, at σ = σ c 2 the unstable and stable branches of R merge via a saddle-node infinite-period bifurcation, whereby the limit cycle solution of the SW arises [see Fig. 3(b)]. A saddle node appears again at σ = σ c 1 and the system is brought back to the partially synchronized state.
However, besides the branch of R obtained under the symmetry condition r 1 = r 2 = 1, the numerical evolution of the reduced system predicts the existence several other curves, which are upper bounded by R of the r 1,2 = 1 solution for σ ∈ [σ c 4 , σ c 3 ]. Insights about the nature of such states can be gained by investigating the stability of the SW state under perturbations transversal to the symmetric manifold r 1 = r 2 [11]. By defining the transversal and longitudinal coordinates r ⊥ = (r 1 − r 2 )/2 and r = (r 1 + r 2 )/2, we havė Linearization at a point (b 0 , δ 0 ) lying on a limit cycle solution of the manifold r 1 = r 2 yields the variation equationδb ⊥ = λ ⊥ δb ⊥ , where λ ⊥ = −(β + kσ )(1 + b 0 ) cos δ 0 . By averaging the previous equation over a period of oscillation and using the periodicity of the limit cycle ( (d/dt ) ln b 0 = 0) we find λ ⊥ = −2(β + kσ ) cos δ 0 . Numerical calculations with the original system (1) for extensive parameter combinations show that cos δ 0 ≈ 0 (and, consequently, λ ⊥ ≈ 0) in the region σ ∈ [σ c 2 , σ c 1 ], suggesting then that the limit cycle solution of the SW is neutrally stable. Although our numerical estimate does not give us an exact proof of the stability of the limit cycle solution, it sheds light on the existence of the multitude of curves observed in Fig. 3(b). Essentially, since nearby trajectories are not attracted nor repelled, any perturbation of the SW state leads to a new limit cycle with r 1,2 < 1, explaining the origin of the numerous solutions encountered in the coupling region encompassed by σ c 4 and σ c 3 in Fig. 3(b). Notice also that the lower branches in the region σ ∈ [σ c 4 , σ c 3 ] do not correspond precisely to the classical incoherent state but rather represent limit cycle solutions with small amplitudes r 1,2 .

B. Heterogeneous uncorrelated networks
The theory developed for random regular networks can also give us insights into the dynamics of networks with mildly heterogeneous degree distributions. In Fig. 4(a) we superimpose numerical results for ER networks with the branches for the partially synchronized states [Eq. (7)] and critical conditions given by Eqs. (6) and (8). Interestingly, we see that the dependence of the order parameters is reproduced with good precision with the expressions derived for simpler networks. Boundaries enclosing the breathing SW states in the random regular network also delineate the region with global oscillations for the ER network. Notice also that σ c 1,2 again marks the limits of the partially synchronized branch of R; however, no state analogous to the perfectly symmetric SW (r 1,2 = 1) is observed in σ c 2 < σ < σ c 1 for ER networks.
Let us take a step further in the analysis of heterogeneous structures and consider general uncorrelated networks with degree distribution p k . In this case, we assume that nodes with the same degree k admit the same solution, i.e., α i = α k , if where α k,1 describes the dynamics of oscillators with degree k and frequency ω i = /2; equations for coefficients α k,2 standing for the second face of Janus oscillators (ω i = − /2) are obtained accordingly. Linearizing the system around α k,1 = δα k,1 1 yields the following variational system: where δα is a small perturbation of the complex order parameter α = 1 2 k k k p k (α k,1 + α k,2 ), and δᾱ is the analogous quantity for the parameter measuring the difference in the internal synchrony in Janus oscillators, i.e.,ᾱ = 1 2 k k k p k (α k,1 − α k,2 ). The eigenvalues of the Jacobian matrix of system (9) become purely imaginary at =|β+σ k 2 / k |; therefore, we predict the appearance of states with global oscillations in the range We compare the predictions of the equation above in Fig. 4(b) for SF networks with the power-law exponent γ = 2.25. At first sight, it seems that the condition σ (net) provides an inaccurate estimation of the region where the order parameters R and r 1,2 are expected to exhibit an erratic behavior, suggesting perhaps that finite-size effects could be behind the deviation between σ (net) c 1 and the point σ 0.17 at which the branch of R collapses to values R ≈ 0. However, Eq. (10) refers to the coupling range in which multiple oscillating states are expected to emerge. Visualizing in Fig. 4(c) the evolution of the mean-field frequency˙ , we observe that σ (net) c 1,2 actually define very accurately the boundaries of the states with multiple oscillating solutions (σ (net) c 2 ). Given the dependence on k / k 2 , one further envisions from Eq. (10) the absence of such oscillating states in the thermodynamic limit for SF networks with 2 < γ 3, since σ (net) c 1,2 are expected to vanish as N → ∞ in similar fashion as the classical coupling for the onset of synchronization in such structures [2,12].
It is noteworthy that although we have considered negative and positive couplings, we have not observed π states, which are collective phenomena that are characteristic of the interplay between attractive and repulsive interactions [13]. On the other hand, Fig. 4(c) (see also Appendix B) indicates the emergence of multiple traveling waves [13] in the partially synchronized state.

IV. CONCLUSION
In conclusion, we have explored the collective dynamics of Janus oscillators on large homogeneous and heterogeneous random networks. By employing the OA ansatz, we obtained for random regular networks a reduced set of equations whereby critical points of the dynamics were revealed. We found that several collective behaviors coexist for intermediate coupling values, elucidating the findings in [6]. Although initially obtained for homogeneous networks, we verified that the solutions of the reduced system fitted accurately numerical experiments for dense ER networks. By analyzing the stability of general uncorrelated networks, we further verified that the coupling range for which global oscillations are possible shrinks in the thermodynamic limit of SF networks. It is pertinent to mention that the accuracy of the OA ansatz in predicting the transition points is deteriorated for σ and β values beyond the region depicted in Fig. 3. Deviations from the temporal signature yielded by the reduced system for solutions of r 1,2 were also observed in simulations for some couplings (β, σ ) in the breathing SW area.
All in all, we provided theoretical and numerical analysis of ensembles of Janus oscillators on homogeneous and heterogeneous random networks. As such, our work raises further interesting questions about the study initiated by Nicolaou et al. [6]. For instance, future investigations should target the dynamics on sparse and correlated networks-situations in which the ensemble approach in [10] and mean-field techniques are inaccurate in predicting tipping points of the system-as well as limitations of the OA manifold in capturing the Janus dynamics. The funders had no role in study design, data collection and analysis, or preparation of the manuscript.

APPENDIX A: DERIVING THE REDUCED SYSTEM WITH THE OTT-ANTONSEN ANSATZ
This section offers further details of the derivation of the reduced system (5) in the main text. The equations that govern the Janus dynamics can be rewritten aṡ where the oscillators are indexed by i = 1, . . . , 2N. The local order parameters R i are given by where the 2N×2N coupling matrix W is defined as in the main text, that is, with I being the identity matrix and A the matrix accounting for the connections between the Janus oscillators. Following [10], we consider an ensemble of systems defined by Eq. (A1) with fixed coupling matrix W. In this formulation, we describe system (A1) at time t by the joint probability density ρ 2N (θ, ω, t ), where θ = (θ 1 , . . . , θ 2N ) contains the phases at time t, and ω = (ω 1 , . . . , ω 2N ) is the time-independent vector with the natural frequencies of the individual oscillators. The evolution of ρ 2N is then dictated by Multiplying the continuity equation of ρ 2N by j =i dω j dθ j and integrating, we obtain the evolution equation for the marginal oscillator density Following the Ott and Antonsen (OA) ansatz [7], we expand the one-oscillator density ρ i (θ i , ω i , t ) as Our goal now is to find the governing equations for coefficientsα i . In the continuum limit, the distribution ρ i (θ i , ω i , t ) satisfies the continuity equation By calculating the time derivative of Eq. (A5), we get The product inside the second term of Eq. (A6) reads Calculating then the derivative on θ i yields In order for Eq. (A6) be satisfied, all the Fourier coefficients must vanish; it follows then thaṫ 2N ). (A8) For the case in which the adjacency matrix A encodes the topology of a random regular graph characterized by degree k, i.e., in this coupling scheme, one Janus face rotating with frequencyω 1 is connected randomly to k faces with frequencŷ ω 2 and vice versa, we assume the following solution for coefficientsα i :α We also quantify spontaneous drifts in the oscillation frequencies by calculating the mean-field frequency˙ , which quantifies how fast the centroid given by the total order parameter Re i = (1/2)(r 1 e iψ 1 + r 2 e iψ 2 ) rotates. In Figs. 5(d) and 5(e) we clearly see the SW phenomenon, i.e., the separation of the population of oscillators into two counteracting clusters, quantified by the mirrored behavior of frequencies 1 and 2 in the region σ ∈ [σ c 2 , σ c1 ]. Notice also that the modules of these frequencies are slightly different; as it is seen in panels (d) and (e), exhibits in the SW state a modest trend as coupling σ is varied. Furthermore, in panels (f) and (g) we also see that multiple stationary solutions for frequency˙ emerge in the region designated as "Partial sync." At first, one could expect to observe values close tȯ ≈ 0, since the average natural frequency of the system is ω i = 0. However, as shown in Fig. 5, rhythms of oscillation different from the mean natural frequency emerge-a  Fig. 2 of the main text. Panels (d) and (e) show the evolution of the total average frequency and average frequencies of each subpopulation 1,2 for β = −0.3 and 0.25, respectively. For the total average frequency , we highlight two cases with different symbols, namely, " (SW)" denotes the average frequency calculated in the SW state, while " (bSW)" stands for the same quantity calculated now in the breathing-SW state. Panels (f) and (g) depict the locking frequency˙ as a function of coupling σ . phenomenon akin to the so-called traveling wave states, which are typically observed in ensembles of oscillators with positive and negative interactions [13]. Figure 6 shows the corresponding results for ER and SF networks.