Collective dynamics of random Janus oscillator networks

PHYSICAL REVIEW RESEARCH 2, 013255 (2020) Collective dynamics of random Janus oscillator networks Thomas Peron,1,2,*Deniz Eroglu,3,†Francisco A. Rodrigues,1and Yamir Moreno 2,4,5 1Institute of Mathematics and Computer Science, University of São Paulo, São Carlos 13566-590, São Paulo, Brazil 2Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, Zaragoza 50018, Spain 3Department of Bioinformatics and Genetics, Kadir Has University, 34083 Istanbul, Turkey 4Department of Theoretical Physics, University of Zaragoza, Zaragoza 50009, Spain 5ISI Foundation, Torino 10126, Italy (Received 11 September 2019; accepted 29 January 2020; published 4 March 2020) 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 specifically 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 find 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. DOI: 10.1103/PhysRevResearch.2.013255 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]. *[email protected] †[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. 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˝ os-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 NJanus oscillators [6]on heterogeneous networks as ˙ θi=ωi+ 2N  j=1 Wij sin(θj−θi),(i=1,...,2N),(1) 2643-1564/2020/2(1)/013255(9) 013255-1 Published by the American Physical Society PERON, EROGLU, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013255 (2020) + + + + - - - - 12 FIG. 1. Illustration of heterogeneous networks of Janus oscillators. Janus faces in layer “+” (“–”) oscillate with frequency ˆω1(ˆω2). where the 2N×2Nmatrix Wis defined as W=0βI+σA βI+σAT0.(2) Natural frequencies are assigned as ωi=ω0+/2≡ˆω1, for i=1,...,N, and ωi=ω0−/2≡ˆω2,fori=N+ 1,...,2N, where is the frequency mismatch and ω0is 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. 1for 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”1defined in Eqs. (1) and (2). III. THEORY AND SIMULATIONS By defining the local order parameters Ri=2N j=1Wijeiθj, Eqs. (1) are then decoupled as ˙ θi=ωi+Im[e−iθiRi]. 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 tby 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 ωiare distributed according to a function g(ωi). By expanding ρiin Fourier series and applying the OA ansatz to its coefficients [7,10]wehave ρi(θi,ω i,t)=g(ωi) 2π1+∞  n=1 ˆαn i(ωi,t)einθi+c.c.. 1Phase oscillators coupled on bipartite networks have been addressed previously; see, e.g., [14–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: ˙ ˆαi+iˆαiωi+1 2ˆα2 iRi−R∗ i=0(i=1,...2N),(3) where the coefficients Riare calculated as Ri= 2N  j=1 Wij ∞ −∞ 2π 0 ρjeiθjdθjdωj. Inserting the expansion of ρjin the previous equation yields Ri= 2N  j=1 Wij ∞ −∞ ˆα∗ j(ωj,t)g(ωj)dωj. 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=/2is randomly connected to koscillators 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≡α1and ˆαN+1=···= ˆα2N≡α2. Hence, the local order parameters Riare reduced to Ri=α∗ 2(β+kσ)≡R1if i=1,...,N, and Ri=α∗ 1(β+ kσ)≡R2if i=N+1,...,2N. Inserting solutions α1,2and R1,2in Eq. (3), we obtain the reduced set of equations, ˙α1,2+iˆω1,2α1,2+1 2α2 1,2R1,2−R∗ 1,2=0,(4) which in terms of the coordinates α1,2=r1,2eiψ1,2are written as ˙r1=1 2(β+kσ)r21−r2 1cosδ ˙r2=1 2(β+kσ)r11−r2 2cosδ(5) ˙ δ=−−1 2(β+kσ)sinδ2r1r2+r1 r2+r2 r1, where δ=ψ1−ψ2is the phase lag between subpopulations [full derivation of Eq. (5) is given in Appendix A]. The variables r1and r2turn 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|r1eiψt(t)+r2eiψ2(t)|. States that emerge from system (5) are summarized as follows: (1) a partially synchronized state in which r1,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 (r1,2=1) rotate in opposite directions, yielding an incessantly rotating δ; and (3) a distinct form of SW emerges when 0 <r1,2<1: along with the increase or decrease in δ, the order parameters r1,2exhibit 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]. 013255-2 COLLECTIVE DYNAMICS OF RANDOM JANUS … PHYSICAL REVIEW RESEARCH 2, 013255 (2020) (a) (b) FIG. 2. Dynamics of a random regular network with k=40. Temporal evolution of (a) order parameters r1,2measuring the level of synchronization within subpopulations, and total order parameter R=1 2r2 1+r2 2+2r1r2cosδ; and (b) evolution of phase lag δ. Initial conditions: r1(0) =0.61, r2(0) =0.83, δ(0) =2π/3. Parameters: β=0.015, σ=0.005 and =1, N=2500. To uncover the conditions for the existence of the partially synchronized state, we set r1(t)=r2(t)≡r(t)inEqs.(5), leading to ˙r=(1/2)(β+kσ)r(1 −r2)cosδand ˙ δ=−− (β+kσ)(1 +r2)sinδ. Imposing ˙r=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 σc2=−/2+β kand σc1=/2−β k,(6) which determines the coupling range where the partially synchronized state exists. The total order parameter Ris then given by R=1 √2   1±1−2 4(β+kσ)2,(7) where the “–” branch is stable for σ⩽σc1, whereas the “+” branch is stable in the region σ⩾σc2.Forσc1<σ <σ c2,the limit cycle solution of ˙ δholds and SW states emerge with perfectly synchronized subpopulations (r1,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 −0.4−0.3−0.2−0.10.00.10.20.30.4 β −0.06 −0.04 −0.02 0.00 0.02 0.04 0.06 σ Partial sync Partial sync Partial sync/breathing SW Partial sync/breathing SW SW/breathing SW σc4 σc2 σc1 σc3 (a) −0.08 −0.06 −0.04 −0.02 0.00 0.02 0.04 0.06 0.08 σ 0.0 0.2 0.4 0.6 0.8 1.0 Order parameters σc3 σc1 σc2 σc4 (b) r1=r2 R FIG. 3. (a) Stability diagram of system (5)fork=40 and =1. “Partial sync” refers to the state in which r1,2=1andδ= 0. SW denotes parameter regions for each r1,2=1 and for which the phase lag δrotates with a nonzero frequency. The regions where a multitude of solutions with r1,2<1and˙ δ= 0 are found are labeled as “breathing SW.” (b) Order parameter for β=−0.3. Solid and dashed lines of Rfor σ∈(−∞,σ 2]∪[σ1,∞) are obtained with Eq. (7). Line for r1=r2denotes 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 σ∈[σc4,σ c3] 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 σc4≡−(+β) k<σ<−β k≡σc3.(8) Figure 3(a) outlines the critical conditions given by Eqs. (6) and (8) in the plane spanned by couplings βand σ.Asitis 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 (σ<σ c4), the order parameter Rcollapses in the solution given by Eq. (7). As σis further increased, at σ=σc2 the unstable and stable branches of Rmerge 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 013255-3 PERON, EROGLU, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013255 (2020) at σ=σc1and the system is brought back to the partially synchronized state. However, besides the branch of Robtained under the symmetry condition r1=r2=1, the numerical evolution of the reduced system predicts the existence several other curves, which are upper bounded by Rof the r1,2=1 solution for σ∈[σc4,σ c3]. 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 r1=r2 [11]. By defining the transversal and longitudinal coordinates r⊥=(r1−r2)/2 and r=(r1+r2)/2, we have ˙r⊥=1 2(β+kσ)(r2 ⊥−r2 −1)r⊥cosδ, which in terms of variables b=r2 and b⊥=r2 ⊥reads ˙ b⊥=(β+kσ)b⊥(b⊥−b−1)cosδ. Linearization at a point (b0,δ 0) lying on a limit cycle solution of the manifold r1=r2yields the variation equation ˙ δb⊥= λ⊥δb⊥, where λ⊥=−(β+kσ)(1 +b0)cosδ0. By averaging the previous equation over a period of oscillation and using the periodicity of the limit cycle ((d/dt)lnb0=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 σ∈[σc2,σ c1], 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 r1,2<1, explaining the origin of the numerous solutions encountered in the coupling region encompassed by σc4and σc3 in Fig. 3(b). Notice also that the lower branches in the region σ∈[σc4,σ c3] do not correspond precisely to the classical incoherent state but rather represent limit cycle solutions with small amplitudes r1,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 σc1,2again marks the limits of the partially synchronized branch of R; however, no state analogous to the perfectly symmetric SW (r1,2=1) is observed in σc2<σ<σ c1for ER networks. Let us take a step further in the analysis of heterogeneous structures and consider general uncorrelated networks with degree distribution pk. In this case, we assume that nodes with thesamedegreekadmit the same solution, i.e., αi=αk,if −0.06 −0.04 −0.02 0.00 0.02 0.04 0.06 0.00 0.25 0.50 0.75 1.00 Order parameters (a) R r1 r2 −0.04 −0.02 0.00 0.02 0.04 0.00 0.25 0.50 0.75 1.00 Order parameters σ(net) c2 σ(net) c1 (b) −0.04 −0.02 0.00 0.02 0.04 σ −0.02 0.00 0.02 ˙ Θ (c) FIG. 4. (a) Order parameters Rand r1,2for ER networks with k=40. Dots are numerical experiments with N=2500 oscillators. To highlight the dynamical similarity between random regular and dense ER networks, critical conditions and branches in this panel are obtained as in Fig. 3. (b) Synchronization curves of an uncorrelated SF network with pk∼k−γ,whereγ=2.25 and N= 104oscillators. Conditions σ(net) 1,2are given by Eq. (10). (c) Evolution of mean-field frequency ˙ associated with the order parameter R(t)ei(t)=(r1eiψ1+r2eiψ2)/2 for the same SF network of panel (b). Internal coupling is set to β=0.25 in all panels. ki=k. Thus, Eqs. (3) are reduced to ˙αk,1+iωαk,1+1 2α2 k,1βα∗ k,2+σk k k kpkα∗ k,2 −βαk,2+σk k k kpkαk,2=0, where αk,1describes the dynamics of oscillators with degree kand 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: ˙ δα =−iωδ ¯α+1 2β+σk2 kδα, ˙ δ¯α=−iωδα −1 2β+σk2 kδ¯α, (9) 013255-4 COLLECTIVE DYNAMICS OF RANDOM JANUS … PHYSICAL REVIEW RESEARCH 2, 013255 (2020) where δα is a small perturbation of the complex order parameter α=1 2kkkpk(α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kkpk(αk,1−αk,2). The eigenvalues of the Jacobian matrix of system (9) become purely imaginary at =|β+σk2/ k|; therefore, we predict the appearance of states with global oscillations in the range σ(net) c1≡−(+β)k k2<σ< k k2(−β)≡σ(net) c2.(10) 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) c1provides an inaccurate estimation of the region where the order parameters R and r1,2are expected to exhibit an erratic behavior, suggesting perhaps that finite-size effects could be behind the deviation between σ(net) c1and the point σ≃0.17 at which the branch of Rcollapses 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) c1,2actually define very accurately the boundaries of the states with multiple oscillating solutions (σ(net) c1<σ<σ(net) c2). Given the dependence on k/k2, one further envisions from Eq. (10)the absence of such oscillating states in the thermodynamic limit for SF networks with 2 <γ ⩽3, since σ(net) c1,2are 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 r1,2were 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. ACKNOWLEDGMENTS T.P. acknowledges FAPESP (Grants No. 2016/23827-6 and No. 2018/15589-3). D.E. acknowledges TÜB˙ ITAK (Grant No. 119F125) and a Kadir Has University internal Scientific Research Grant (BAF). F.A.R. acknowledges CNPq (Grant No. 305940/2010-4) and FAPESP (Grants No. 2016/256825 and No 2013/07375-0) for the financial support given to this research. Y.M. acknowledges partial support from the Government of Aragon, Spain, through Grant No. E36-17R (FENOL), by MINECO and FEDER funds (FIS2017-87519P), and by the Intesa Sanpaolo Innovation Center. This research was carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (Grant No. 2013/07375-0). 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 as ˙ θi=ωi+1 2i[e−iθiRi−eiθiR∗ i],(A1) where the oscillators are indexed by i=1,...,2N. The local order parameters Riare given by Ri= 2N  j=1 Wijeiθj,(A2) where the 2N×2Ncoupling matrix Wis defined as in the main text, that is, W=0βI+σA βI+σAT0, with Ibeing the identity matrix and Athe 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 tby 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 ρ2Nis then dictated by ∂ρ2N ∂t+ 2N  i=1 ∂ ∂θi(ρ2N˙ θi)=0.(A3) 013255-5 PERON, EROGLU, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013255 (2020) Multiplying the continuity equation of ρ2Nby j=idωjdθj and integrating, we obtain the evolution equation for the marginal oscillator density ρi(θi,ω i,t)=ρ2Nj=idωjθj, that is, ∂ρi ∂t+∂ ∂θi(ρi˙ θi)=0.(A4) Following the Ott and Antonsen (OA) ansatz [7], we expand the one-oscillator density ρi(θi,ω i,t)as ρi(θi,ω i,t)=g(ωi) 2π1+∞  n=1 ˆαn i(ωi,t)einθi+c.c..(A5) 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 ∂ρi ∂t+∂ ∂θi(ρi˙ θi)=0.(A6) By calculating the time derivative of Eq. (A5), we get ∂ρi ∂t=g(ωi) 2π∞  n=1 nˆαn−1 i˙ ˆαieinθi+n(ˆα∗ i)n−1˙ ˆα∗ ie−inθi.(A7) The product inside the second term of Eq. (A6) reads ρi(θi,ω i,t)˙ θi=g(ωi) 4πi2iωi+ˆαiRi−ˆα∗ iR∗ i +∞  n=1ˆαn+1 iRi−ˆαn−1 iR∗ i+2iωiˆαn ieinθi +∞  n=1ˆαn−1 iRi−ˆαn+1 iR∗ i+2iωiˆαn ie−inθi. Calculating then the derivative on θiyields ∂ ∂θi(ρi˙ θi)=g(ωi) 4πi ∞  n=1ˆαn+1 iRi−ˆαn−1 iR∗ i+2iωiˆαn i(in)einθi +g(ωi) 4πi ∞  n=1ˆαn−1 iRi−ˆαn+1 iR∗ i+2iωiˆαn i ×(−in)e−inθi. In order for Eq. (A6) be satisfied, all the Fourier coefficients must vanish; it follows then that ˙ ˆαi+iωiˆαi+1 2ˆα2 iRi−R∗ i=0(i=1,...,2N).(A8) For the case in which the adjacency matrix Aencodes the topology of a random regular graph characterized by degree k, i.e., in this coupling scheme, one Janus face rotating with frequency ˆω1is connected randomly to kfaces with frequency ˆω2and vice versa, we assume the following solution for coefficients ˆαi: ˆα1=ˆα2=···= ˆαN≡α1, ˆαN+1=···= ˆα2N≡α2.(A9) By substituting the solutions above in Eq. (A8), we get Ri=α∗ 2(β+kσ)≡R1if i=1,...,N; α∗ 1(β+kσ)≡R2if i=N+1,...,2N.(A10) Therefore, Eqs. (A8) are reduced to ˙α1,2+iωα1,2+1 2α2 1,2R1,2−R∗ 1,2=0.(A11) By writing in terms of the coordinates α1,2=r1,2eiψ1,2,we have ˙r1=−β+kσ 2r2 1r2cos(ψ1−ψ2)−r2cos(ψ1−ψ2), ˙ ψ1=−ω1+β+kσ 2r1r2sin(ψ1−ψ2)−r2 r1sin(ψ1−ψ2). (A12) Equations for r2and ψ2are obtained by interchanging labels “1” and “2” in the equation above. By defining the phase lag δ=ψ1−ψ2,Eqs.(A12) are further reduced to ˙r1=1 2(β+kσ)r21−r2 1cosδ ˙r2=1 2(β+kσ)r11−r2 2cosδ(A13) ˙ δ=−−1 2(β+kσ)sinδ2r1r2+r1 r2+r2 r1. APPENDIX B: ADDITIONAL DIAGRAMS In this section we present complementary diagrams related to the Janus dynamics on random regular networks. Besides the stability diagram [Fig. 5(a)] and the curves of the order parameters for β=−0.3 [Fig. 5(b)] present in the main text, here we also show information about collective frequencies. Specifically, we consider as the total average frequency the quantity given by =1 2N 2N  j=1dθj dt t ,(B1) where ·tis a temporal average. Analogously, we define the average frequency of each group of Janus faces as 1=1 N N  j=1dθj dt tand 2=1 N 2N  j=N+1dθj dt t .(B2) 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 Rei=(1/2)(r1eiψ1+r2eiψ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 1and 2in the region σ∈[σc2,σ 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 to ˙ ≈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 013255-6 COLLECTIVE DYNAMICS OF RANDOM JANUS … PHYSICAL REVIEW RESEARCH 2, 013255 (2020) FIG. 5. (a) Stability diagram of system (5) considering a random regular network of Janus oscillators with k=40 and frequency mismatch =1. Dashed lines depict the regions in the parameter space considered in the following panels. Descriptions of panels (b) and (c) follow as in Fig. 2of the main text. Panels (d) and (e) show the evolution of the total average frequency and average frequencies of each subpopulation 1,2for β=−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 σ. 013255-7 PERON, EROGLU, RODRIGUES, AND MORENO PHYSICAL REVIEW RESEARCH 2, 013255 (2020) FIG. 6. (a) Order parameters, (b) average frequencies, and (c) locking frequency ˙ of ER networks with N=2500 Janus oscillators, and average degree k=40. Panels (d), (e), and (f) display the corresponding results for SF networks with P(k)∼k−γ,whereγ=2.25, and N=104Janus oscillators. Other parameters used in all panels: β=0.25 and =1. phenomenon akin to the so-called traveling wave states, which are typically observed in ensembles of oscillators with positive and negative interactions [13]. Figure 6shows the corresponding results for ER and SF networks. [1] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469,93 (2008). [2] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, Phys. Rep. 610,1(2016). [3] M. J. Panaggio and D. M. Abrams, Nonlinearity 28,R67 (2015). [4] J. Gómez-Gardenes, S. Gómez, A. Arenas, and Y. Moreno, Phys. Rev. Lett. 106,128701 (2011). [5] T. Nishikawa and A. E. Motter, Phys. Rev. Lett. 117,114101 (2016); Y. Zhang, T. Nishikawa, and A. E. Motter, Phys. Rev. E 95,062215 (2017). [6] Z. G. Nicolaou, D. Eroglu, and A. E. Motter, Phys. Rev. X 9, 011017 (2019). [7] E. Ott and T. M. Antonsen, Chaos 18,037113 (2008). [8] A. Pikovsky and M. Rosenblum, Chaos 25,097616 (2015). [9] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, andM.A.Porter,J. Complex Networks 2,203 (2014). [10] G. Barlev, T. M. Antonsen, and E. Ott, Chaos 21,025103 (2011). [11] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. M. Antonsen, Phys.Rev.E79,026204 (2009). 013255-8 COLLECTIVE DYNAMICS OF RANDOM JANUS … PHYSICAL REVIEW RESEARCH 2, 013255 (2020) [12] T. Peron, B. Messias F. de Resende, A. S. Mata, F. A. Rodrigues, and Y. Moreno, Phys. Rev. E 100,042302 (2019). [13] H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106,054102 (2011); B. Sonnenschein, T. K. DM. Peron, F. A. Rodrigues, J. Kurths, and L. Schimansky-Geier, Phys.Rev.E91,062910 (2015); K. O’Keeffe and C. Bettstetter, in Micro-and Nanotechnology Sensors, Systems, and Applications XI (International Society for Optics and Photonics, Bellingham, WA, 2019), Vol. 10982, p. 109822E. [14] K. A. Samani and S. Ghanbarian, Phys.Rev.E77,036209 (2008). [15] N. Punetha, R. Ramaswamy, and F. M. Atay, Phys.Rev.E91, 042906 (2015). [16] E. Montbrió and D. Pazó, Phys.Rev.Lett.120,244101 (2018). 013255-9