Multifaceted dynamics of Janus oscillator networks

Recent research has led to the discovery of fundamental new phenomena in network synchronization, including chimera states, explosive synchronization, and asymmetry-induced synchronization. Each of these phenomena has thus far been observed only in systems designed to exhibit that one phenomenon, which raises the questions of whether they are mutually compatible and, if so, under what conditions they co-occur. Here, we introduce a class of remarkably simple oscillator networks that concurrently exhibit all of these phenomena. The dynamical units consist of pairs of nonidentical phase oscillators, which we refer to as Janus oscillators by analogy with Janus particles and the mythological figure from which their name is derived. In contrast to previous studies, these networks exhibit (i) explosive synchronization with identical oscillators; (ii) extreme multistability of chimera states, including traveling, intermittent, and bouncing chimeras; and (iii) asymmetry-induced synchronization in which synchronization is promoted by random oscillator heterogeneity. These networks also exhibit the previously unobserved possibility of inverted synchronization transitions, in which a transition to a more synchronous state is induced by a reduction rather than an increase in the coupling strength. These various phenomena are shown to emerge under rather parsimonious conditions, and even in locally connected ring topologies, which has the potential to facilitate their use to control and manipulate synchronization in experiments.


I. INTRODUCTION
It has been the tradition of physics to describe complex behavior using simple mathematical models. In such a description, of which the phenomenon of chaos offers many compelling examples [1], the behavioral complexity is emergent rather than explicitly coded in the model. In recent studies of network dynamics, the possibility of devising a complex network structure, and thus a complex model, adds a new dimension to this tradition and allows for easy departure from parsimonious models. This has lead to the discovery of fascinating new phenomena but has also created difficulties in isolating the minimal requirements for the observed dynamical behavior. This trend is partially countered by the discovery of chimera states [2,3], characterized by coexisting incoherence and synchrony in identically coupled identical oscillators; explosive synchronization transitions [4], in which transition to synchronization becomes abrupt; and asymmetryinduced synchronization (AIS) [5][6][7], a partial converse to the symmetry breaking exhibited by chimera states, in which either the oscillators or their couplings need to be non-identical for synchronization to prevail. These behaviors are unequivocally emergent, since they are not manifestly forged into the model. Yet, previous demonstrations of these various phenomena required the design of specific systems, in which the occurrence of the different types of behavior seemed to require different and sometimes complicated types of intrinsic dynamics, interaction structures, and coupling schemes.
Here, we demonstrate the co-occurrence of chimera states, explosive synchronization, and a new form of AIS in a class of surprisingly simple oscillator networks. The dynamical units in these networks are two-dimensional phase-phase oscillators, which we term Janus oscillators by analogy with the homonymous two-faced particles (and the two-faced ancient Roman deity on which that name was based) [8], since the components of the phasephase pair are taken to have different natural frequencies. Figure 1 schematically shows an especially simple network we consider and the various dynamical behaviors it exhibits as a function of coupling strength and oscillator heterogeneity. Importantly, our analysis of this new class of systems demonstrates for the first time the occurrence of (i) inverted synchronization transitions, (ii) a plurality (a) Model network of phase-phase oscillators, in which phase oscillators with alternating frequencies are sinecoupled to their nearest neighbors on a ring. Regarding each pair of plus-minus oscillators as a Janus oscillator results in a symmetric system of oscillators. In such a ring (shown here with N = 6 oscillators), any two Janus oscillators can be swapped by rotations without changing the overall configuration. (b) The various phenomena exhibited by a ring of N = 50 Janus oscillators as the coupling strength as well as the heterogeneity across the Janus oscillators are varied. Quantitative details about these regimes are presented in Section IV A. of chimera-like states, (iii) explosive synchronization in the absence of correlations between oscillator frequency and network structure, and (iv) synchronization induced by random oscillator heterogeneity. In particular, adding small oscillator heterogeneity destroys partially coherent states and makes fully phase-locked states more attractive, which represents a new form of AIS; and, for a certain range of coupling strength, further increase in heterogeneity stabilizes a new type of dynamical behavior, which we call chimera intermittency.
Our model is partially inspired by the antiferromagnetic order first characterized in certain spin systems described below, but that has since repeatedly reappeared in different physical systems as one of several universality classes of order. Janus oscillator models are an appropriate description of the oscillatory dynamics that can emerge in any driven system that exhibits antiferromagnetic-like order in equilibrium.
The paper proceeds with a description of the model in Sec. II, the presentation of the results in Secs. III and IV, and a discussion of further implications in Sec. V. Our presentation is complemented by an animated visualization of the main findings, which is included as Supplemental Material [9] and can be consulted before or after the text.

A. The model
We define a Janus oscillator as a 2-dimensional phasephase oscillator in which each component has a distinct natural frequency. We first consider ring networks of n such oscillators, whose dynamics are governed bẏ for i = 0, . . . , n − 1, where subscripts indicate the Janus oscillator index while superscripts indicate the variable index. The parameters ω 1 i and ω 2 i are the natural frequencies of oscillator i, and the periodic boundary conditions are assured through the index convention i = i mod n, where mod is the modulo operation. Further, the frequencies are assumed to be ω 1 i =ω − ω/2 and ω 2 i =ω + ω/2, whereω and ω are constants. This corresponds to the nearest-neighbor rotationally-symmetric network of identical Janus oscillators (θ 1 i , θ 2 i ) illustrated in Fig. 1(a). Finally, β is the internal coupling strength between the phase-oscillator components of each Janus oscillator (i.e., in the same node) and σ is the external coupling strength between oscillators in different nodes [10].
The full parameter space of Eqs. (1)-(2) is fourdimensional, with coordinate axesω, ω, σ, and β. However, this space can be substantially reduced without loss of generality. First, by changing to the co-rotating reference frame, we can eliminateω so that each pair of oscillators has opposite natural frequencies. Second, by scaling time by 1/ω and σ and β by ω, it is possible to set ω = 1. Lastly, if either the σ or β are too large (i.e., greater than the critical value ω/2 at which the isolated pair of oscillators would synchronize), the two phase components of each oscillator phase-lock and the dynamics can be reduced to a ring of single-phase oscillator nodes. Thus the dynamics of interest on the ring can be captured by fixingω = 0, ω = 1, and β = ω/4, and varying the external coupling strength between 0 ≤ σ ≤ 0.6 [11].
B. Comparison with the existing models As noted above, our model is inspired in part by the image of an antiferromagnet driven by an external magnetic field. The ground state of an antiferromagnet consists of an alternating configuration of magnetic dipoles arranged on a lattice. When an external magnetic field is applied to this lattice, the dipoles precess with alternating angular frequencies ±ω/2, as illustrated in Fig. 2(a). The dynamics and synchronization of such coupled dipoles have recently found applications in spintronics. For example, the synchronization of arrays of spin torque and spin Hall nano-oscillators through coupling currents has been successfully described with the Kuramoto model (and with more complex models) and experimentally realized [12][13][14][15][16]. Our model could thus be applicable to an antiferromagnetic array of spin torque or spin Hall nano-oscillators. While still in development, the possibility of designing spintronic devices from antiferromagnets rather than ferromagnets has also recently attracted increasingly detailed dynamical modeling [17,18]. Another notable system that may be modeled with Janus oscillators is the counter-rotating flagella in the cells of certain communities of algae called Chlamydomonas. The geometry of Chlamydomonas cells is shown schematically in Fig. 2(b). In communities of Chlamydomonas, both internal cellular interactions and hydrodynamic interactions in the cellular environment are thought to induce local coupling between the phases describing each flagellum's position [19,20], which rotate with opposite natural frequencies, much like our Janus oscillators. Striking patterns of synchronization in communities of cells bearing single flagella have been noted [21], and we propose that communities of cells like Chlamydomonas may rely on even more complex synchronization dynamics similar to the Janus oscillator dynamics we show here.
Previous studies on chimera states have also considered minimal models in order to gain a more complete understanding of this complex phenomenon. One such example is the solvable model of two populations of coupled oscillators [22]. Another example is the chimera identified in a globally coupled population of oscillators with delay feedback coupling [23]. Moreover, while initially found in nonlocally-coupled oscillators, chimeras have also been noted in strictly locally coupled-models of phase-amplitude oscillators in both networks [27][28][29][30][31] and continuous systems [32]. Our phase-phase oscillator model constitutes a particularly simple totally symmetric system (i.e., both vertex and edge transitive in the language of graph theory) which exhibits chimeras.
Likewise, discontinuous synchronization transitions were already known to occur in fully connected networks of inertial oscillators [33] and of phase oscillators with equidistantly distributed natural frequencies [34]. However, abrupt phase transitions in complex networks attracted a great deal of attention mainly after the discovery of explosive percolation on random [35] and scale-free graphs [36]. Explosive synchronization was first observed numerically in scale-free networks [4] and the discontinuity of the transition was later proven analytically for star networks [37]. Networks of globally coupled oscillators with bimodal frequency distributions have been shown to exhibit explosive, subcritical synchronization transitions as well [38,39]. Here, in contrast with previous studies focused on explosive transitions from an asynchronous state to a single synchronous state, we explore an entire spectrum of transitions to a multitude of partially synchronous states. Furthermore, the loss of this massive multistability with increasing heterogeneity leads to a new mechanism for asymmetry-induced synchronization.
In previous models, the traditional way to measure the degree of global synchronization of the whole population of oscillators has been the magnitude of Kuramoto order parameter which ranges from 0 to 1 with larger values typically corresponding to more synchronous states [40,41]. However, the order parameter is not a perfect measure of the degree of synchronization. For example, the order parameter is r = 0 for equidistantly distributed and phase-locked oscillators, and thus the order parameter of a partially coherent state may be higher than a fully phase-locked state. Therefore, the order parameter can be misleading for some cases and it is important to quantify such solu-tions with extra care. Nevertheless, r is a useful quantity to distinguish dynamical behaviors that we consider here. In addition, to distinguish cases in which r does not adequately quantify synchronization, we introduce additional metrics below.

III. THE SYMMETRIC CASE
We start with the scenario described above in which all Janus oscillators are identical and they are identically coupled through a ring network topology. Under such conditions, the system is symmetric in the sense that any two nodes can be swapped while leaving the equations of motion invariant, and fully synchronous states are the solutions that inherit that symmetry. This is therefore a suitable model to investigate the interaction between symmetry breaking and synchronization phenomena.

A. Numerical observations
First we show a global depiction of the solutions branches and the transitions between those branches in the case of a ring of n = 50 Janus oscillators with no heterogeneity in Fig. 3. To identify all the branches of solution, simulations with 10 4 random initial conditions were performed for each value of σ = 0.30, σ = 0.35, and σ = 0.40. These simulations waited for transients to die out before averaging the order parameter r and the mean number of phase-locked oscillators N locked . If all oscillators become phase-locked, the state was deemed phase-locked; if some but not all of oscillators were phaselocked, the state was deemed partially locked; otherwise asynchronous. The final values of r and N locked were then binned to identify initial conditions that resulted in the same final state, and the number of initial conditions that ended in each state was counted to estimate the size of the basin of attraction for each state. To map out the solution branches as a function of the coupling constant σ, simulations were performed that quasistatically vary σ for each state identified from the random initial conditions. The network diagrams on the left and right in Fig. 3 shows, respectively, the probability of transitions between various states under quasistatic changes in σ (as determined through 10 3 simulated transitions from each solution branch) in the decreasing and increasing directions, respectively. The nodes in these network diagrams represent the dynamical states with identical timeaveraged values of both r and N locked , corresponding to each solution branch we identified. They are arranged in space according to their critical coupling constants and order parameters.

Solution branches
The time evolution of several states is shown in Fig. 4. After reaching a critical coupling, explosive synchronization occurs when the branch of the asynchronous state, shown in Fig. 4(a), disappears. The order parameter varies discontinuously at this critical coupling constant as the system jumps into another solution branch. The loss of stability for the asynchronous branch occurs first through chimera intermittency, shown in Fig. 4(b). This intermittency is transient and eventually the system settles into either a chimera state like the one in Fig. 4(c) or a fully phase-locked state like the twisted state shown in Fig. 4(d). As shown below, the transient chimera intermittency can become stable and thus persistent when the system is not constrained to be symmetric.
The main mechanism generating the abrupt transition in explosive synchronization is the following. In scale-free networks, some oscillators in the system partially synchronize in clusters according to the similarities in their degrees and natural frequencies. Similarly, in symmetric networks, oscillators whose frequencies are symmetrically assigned can synchronize even when there is no direct interaction between them. This kind of synchronization was first described as "indirect synchronization" [42], but later dubbed remote synchronization [43,44]. In our model, there exist such a remotely synchronized solution branch for arbitrarily small coupling constants, but this solution is neutrally stable in the linear stability analysis and is not attractive. This remotely synchronized state coexists with the asynchronous state, where no such clustering occurs. However, the synchronized clusters in this state become phase-locked with each other when the coupling increases past a critical value σ sync c = 0.25, and at this point this synchronous state becomes stable and attractive. This phase-locking bifurcation is an example of a saddle node bifurcation on the invariant circle [45], in which the center manifold of the saddle-node is the limit cycle corresponding to the remotely synchronized state. As the coupling constant quasistatically increases, eventually the asynchronous solution ceases to be stable at the critical coupling σ async c ≈ 0.34, and the system must move to either phase-locked cluster solutions or to other partially phase-locked solutions (see Fig. 3) during this explosive synchronization. The intermittent transient behavior shown in Fig. 4(b) occurs as σ increases quasistatically past σ async c . This behavior is reminiscent of spatiotemporal intermittency, which has been investigated, e.g., in the transition to turbulence [46].
In our simulations, we observed two kinds of phaselocked states: the twisted state with small order parameter where the phase increases through 2π around the ring depicted in Fig. 4 (d), and the uniform phase state with large order parameter, where all the oscillators have similar phase, which is qualitatively similar. For intermediate coupling constants, several partially locked solution branches were observed with chimera-like collective behaviors, like that shown in Fig. 4 (c).  shows a selection of these chimera states in more detail. These chimeras exist over differing intervals between 0.28 σ 0.50. The chimeras in our system exhibit a symmetric behavior in which a domain of neighboring oscillators are phase-locked and coexist with an asynchronous domain, as in standard chimeras, but these domains propagate to the left or right and visit each oscillator equally in time. Different chimera states with dynamical boundaries have been observed in several other networks of coupled oscillators. For example, the random meandering of chimera states has been previously noted as undesirable, and a control scheme has been proposed to affix them [47]. Another example is the breathing chimeras [22] in the solvable two-population model. Traveling chimera states have also been identified in phase and limit cycle oscillators with various non-local coupling schemes [24][25][26]. The local coupling scheme in the ring of Janus oscillators consistitutes a particularly simple model in which to study travelling chimera states.

Discontinuous and inverted transitions
A fundamental question associated with the observed multiplicity of stable solution branches concerns the nature of the transitions between branches as the coupling strength σ is varied. Among all the quasistatic transitions, those out of the asynchronous branch were the richest, with many possible final states. Figure 6(a) shows the results for detailed simulations of 10 4 transitions from the asynchronous branch (dashed line) as the coupling strength is quasistatically varied in the range 0.24 ≤ σ ≤ 0.35. It follows that all solutions transition directly back to the asynchronous branch when the corresponding branches come to an end, forming the various hysteresis loops shown in the figure. A snapshot of the phases before and after one such transition is shown in Fig. 6(b). As this example explicitly shows, the tran-sitions from the asynchronous branch to phase-locked branches are discontinuous and constitute genuine manifestations of explosive synchronization even though the network structure is a regular graph.
Closer inspection of the simulations used to generate Fig. 3 shows that not all transitions occur in the expected direction, namely from less to more synchronous states as σ is increased and vice versa. As shown in Fig. 7(a) for simulations in the range 0.3 ≤ σ ≤ 0.4, several transitions are inverted: these are transitions between partially phase-locked solution branches in which an increase (decrease) in the coupling strength results in a decrease (or increase) in synchronization. Explicit examples of inverted transitions are shown in Fig. 7(b)-(c), where states for larger coupling strengths are visually less synchronous than the those for smaller coupling strengths. All inverted transitions shown enjoy this defining characteristic and correspond to branches whose relative degree of synchronization is suitably measured by both the mean number of phase-locked oscillators N locked and the order parameter r. We excluded from Fig. 7(a) the transitions between asynchronous branch and the bottom two twisted branches, which appear inverted with respect to r but this perception is in fact an artifact of the defi- nition of r noted above, which does not properly reflect the synchrony of the twisted state. For the forward inverted transitions in Fig. 7(b), a partially locked state becomes unstable under quasistatic increase in the coupling constant, and chimera intermittency follows for a brief period. This effectively randomizes the state, and there is some probability that once the system settles into an attractor, the final state is less synchronous than the initial state. On the other hand, in the backward inverted transition in Fig. 7(c), a small portion of the incoherent domain detaches from the rest of the incoherent domain on a quasistatic decrease in the couple constant and joins the coherent domain. This process is not random, but happens consistently under quasistatic decrease in the coupling constant from this initial state. The animation in the Supplemental Material is especially informative for understanding the dynamics of these inverted transitions [9]. The inverted transitions characterized here are analogous to negative compressibility transitions identified in mechanical metamaterials [48].

B. Analytical results
These rigorous numerical observations motivated us to pursue an analytical characterization of the observed patterns, which is possible by taking advantage of the symmetry in the model to identify a large class of exact solution reductions. These reductions are parameterized by the number of clusters, 2q, representing the number of phase-locked groups of oscillators, and the twisting wavenumber, ν, representing the number of 2π phase windings each cluster undergoes around the ring. We propose the following twisted cluster ansatz: where mod is, as before, the modulo operation. When Eqs. (4)-(5) are substituted into Eqs. (1)-(2), the following reduction is derived: where j = 1, . . . , q and φ 1,2 0,q+1 ≡ φ 1,2 q,1 . While this may not appear any simpler on first glance, there are, in fact, only 2q equations in Eqs. (6)-(7) as opposed to the N equations in Eqs. (1)-(2). Accordingly, this represents a significant dimension reduction for small q.
As an application of these exact reductions, consider the two-cluster solution for q = 1. For completely phaselocked two-cluster solutions, the difference η ≡ φ 1 j − φ 2 j is a constant. It follows from Eqs. (6)-(7) that which has real solutions only when |σ| ≥ (ω/2 − β) cos ν. The bifurcation that occurs as σ is decreased below this value is the saddle node bifurcation on the invariant circle discussed above. For ν = 0, 2π/N , these phase-locking critical coupling constants are apparent in the left-most transitions in Fig. 6(a). Similar bifurcations in the cases with higher q and ν can be derived. However, these higher order twisted cluster solutions are not attractors or repellors but rather appear to be non-attracting invariant sets. The asynchronous domains in the chimera solutions in Fig. 5 seem to exhibit an approximate twisted cluster symmetry, as apparent in the visible pattern in the asynchronous domain in Fig. 4(c). Figure 8(a) shows the evolution of a ring Janus oscillators starting from an initial condition near one such unstable limit cycle solution to Eq. (6)-(7) with q = 2 and ν = 2π/12. The ring was chosen with N = 48 oscillators so an integer number of twists fits in the domain, and the initial condition was perturbed at the center to accelerate the decay of the unstable solution. The evolution consists of propagating fronts that closely resemble the patterns of the incoherent domains of a chimera, which is shown in Fig 8(b) for comparison. While outside the scope of this work, it may be possible to interpret the chimera solutions as heteroclinic cycles connecting various cluster twisted solutions to phase-locked solutions [49].
Another possible analytic direction would be to employ the Ott-Antonsen reduction technique, which has been successful at describing many systems of globally coupled oscillators [50]. For example, subcritical explosive synchronization transitions in globally coupled networks with bimodal frequency distributions have been clearly characterized [38,39]. Extensions to nonlocally coupled networks have been proposed [51], but challenges remain in applications to strictly local coupling schemes like the ring of Janus oscillators.

C. General networks
We conclude this section with a comment on more general symmetric networks of Janus oscillators, described by where A = (A ij ) is an adjacency matrix, of which the system (1)-(2) is a special case. The entry A ij of this matrix is assumed to be one when the first oscillator in node i is coupled to the second oscillator of node j, and zero otherwise. Each link represents a bidirectional coupling, but the adjacency matrix is not symmetric in general (and hence can be visualized as describing a di- FIG. 9. Impact of oscillator heterogeneity for one profile. Order parameter r vs. coupling strength σ for (a) small heterogeneity (δ = 0.005), (b) moderate heterogeneity (δ = 0.020), and (c) large heterogeneity (δ = 0.040). The lines indicate the asynchronous branch (green, dotted), fully phase-locked branches (blue, continuous), and partially phase-locked branches (red, dashed). As the heterogeneity increases, the phase-locked solutions move to the right and the partially locked solutions disappear. The blue shade marks regions where only fully phase-locked solutions are present. The red shaded area in (a) marks the region where partially locked states exist and in (c) marks the region of stabilized chimera intermittency. rected interaction network) as it also encodes which of the two oscillators in each node are connected. In this way, we consider networks with k nearest neighbor links on a D-dimensional lattice. Our numerical simulations reveal that D = 1 rings with k > 1 nearest neighbors also exhibit large degree of multistability, including propagating chimeras. Thus, our core results appear to generalize to such networks. On the other hand, for a square lattice topologies in D = 2 dimensions with k = 1 nearestneighbors, numerical simulations have shown that while the explosive synchronization transitions persists, only the asynchronous solution and the fully phase-locked solutions are attractors.

IV. EFFECTS OF HETEROGENEITY
We have seen that the simple system considered above exhibits rich dynamics, with many different solutions branches. This was shown to be the case when the Janus oscillators are all identical and identically coupled, resulting in a globally symmetric system. However, in nature the interacting elements are generally nonidentical and/or nonidentically coupled. Next we consider the counterintuitive implications of breaking the system symmetry by introducing oscillator heterogeneity or disorder in the network structure.

A. Oscillator heterogeneity and AIS
We consider ring networks as above but now for randomly perturbed oscillator frequencies: ω 1,2 i = ±ω/2 + δp 1,2 i , where p 1,2 i are i.i.d. random variables drawn from a uniform distribution in [−1/2, 1/2] and the parameter δ defines the level of heterogeneity. As illustrated in Fig.  9 for one realization of the heterogeneity profile p 1,2 i , the solution branches of the order parameter change in form as δ increases. One might have expected that introducing heterogeneity would decrease the degree of coherence ob-served in this system. However, this is not the case, as the partially phased-locked solutions in Fig. 9(a) disappear for increased heterogeneity, as shown in in Fig. 9(b) for δ = 0.020, leaving the system with only coherent, phaselocked stable solutions for a range of coupling strength. This is evidence of AIS, a recently discovered effect [5] in which oscillator heterogeneity can strengthen synchronization even when the oscillators are identically coupled. However, different from all previously reported cases of AIS [5][6][7], in this scenario the phenomenon is determined by random heterogeneity.
We have seen that for a single, fixed heterogeneity profile, increasing the heterogeneity parameter δ can destroy the partially phase-locked states and cause the totally phase-locked states to dominate, but does this trend hold statistically for random heterogeneity profiles as well? To further establish the occurrence of AIS in our system mediated by this mechanism, we analyze the general effect of varying the heterogeneity profile. Figure 10 shows results for 100 realizations of p 1,2 i as δ is increased from 0 to 0.05. We first examined the impact of heterogeneity on the size of the attraction basins of the various branches. As shown in Fig. 10(a), a moderate level of heterogeneity can significantly increase the percentage of initial conditions that become completely (100%) phase-locked. This is observed for intermediate coupling strengths since the partially phase-locked solutions are destroyed in that region. Since the AIS effect is quite prominent in this system, we can be rather strict and qualify the system as exhibiting AIS only when at least 95% of states (with differing random initial conditions and heterogeneity profiles) become completely phase-locked. In addition, we require that in the absense of heterogeneity (i.e. when δ = 0), a minority of states become phase-locked for these coupling strengths, so that states which are already mostly synchronous in the absense of heterogeneity are excluded. That is, for that range of coupling strengths, the introduction of random heterogeneity does not promote disorder. Instead, heterogeneity makes it overwhelmingly likely that the system will evolve into a phase-locked state. This constitutes a striking example of AIS.
The emergence of AIS in this system is intimately related with explosive synchronization and chimera intermittency. This connection is best understood by considering the critical coupling strengths, as shown in Fig.  10(b). For a small level of heterogeneity, the discontinuous transitions and resulting in explosive synchronization persist for all profile realizations, with the asynchronous branch disappearing at a larger coupling strength σ = σ async c than the coupling strength σ sync c at which the uniform, totally phase-locked synchronous branch appears. The other phase-locked solution branches appear at critical coupling constants that are slightly larger than σ sync c , but all the critical coupling constants follow the trend of σ sync c as δ increases. After increasing the level of heterogeneity past δ ≈ 0.025, where the average critical coupling strengths σ sync c and σ async c coincide, the explosive synchronization disappears and the transition becomes continuous (i.e., supercritical). As δ is increased further, an irregular behavior sets in for a range of σ in the asynchronous branch when all phase-locked branches exist only coupling strengths larger than σ async c . This asynchronous state corresponds to the red shaded region in Fig. 9(c) and is in fact a stabilized and persistent form of the previously transient chimera intermittency shown in Fig. 4(b). Being the only stable solution in this range of coupling parameters, this new chimera state is more accessible than the chimera states in Fig. 5. Thus, the occurrence of AIS in this system can be interpreted as follows: introducing a small amount of heterogeneity into the model inhibits the explosive synchronization and eliminates partially phase-locked solutions, which in turn results in a greater tendency toward phase-locking.
To summarize the effects that oscillator heterogeneity has on the ring, we now describe the details about the various parameter regimes synthesized in Fig. 1(b). The explosive synchronization regime (green) and intermittency regime (red) were determined from the critical coupling constants in Fig. 10(b). To quantify the chimera parameter regime (orange), simulations were performed to quasistatically increase δ and σ starting from chimera initial conditions for 50 different heterogeneity profiles. These simulations then detected when the chimeras ceased existing in order to map the boundary. The exact boundary where chimeras cease to exist depends on the profile, and the region in Fig. 1(b) shows the average δ over the 50 different profiles. To quantify the AIS parameter regime, 1000 simulations with different random initial conditions and heterogeneity profiles were performed for each point in the σ-δ plane. As in Fig. 10(a), the AIS region shows where the vast majority (≥ 95%) of random initial conditions and heterogeneity profiles evolve into completely (100%) phaselocked states; the AIS region is delineated under the additional constraint that < 50% of the states are phaselocked for δ = 0 in order to exclude parameters for which the homogeneous system is already mostly synchronized. The white region between the explosive synchronization, AIS, and chimera regions in Fig. 1(b) represents areas where some heterogeneity profiles still exhibit attractive chimera states while the majority of heterogeneity profiles do not. This is a finite-size effect in which different behaviors will be observed depending on the exact profile of the oscillator heterogeneities.

B. Network disorder and bouncing chimeras
Network heterogeneity in the form of random link addition and/or removal was also studied on a ring with k = 1 and k = 2 nearest-neighbor links. Figure 11 shows typical results on such networks. When random links are added to the network, most chimeras are destroyed but we observe that some chimera solutions persist and bounce off the network defects [ Fig. 11(a)-(c)]. That is, these solutions are defined by a moving group of asynchronous oscillators that remain trapped between two nodes with added links, and constitute what we term bouncing chimera states. On the other hand, when links are randomly removed, the nodes with removed links do not phase-lock with their neighbors and remain unsynchronized for moderate coupling strengths [ Fig. 11(b)-(d)]. Such nodes act as sources for traveling chimeras, resulting in chimera intermittency.
It follows that heterogeneity in the network struc- ture can either inhibit synchronization-as one might expect-or reduce the degree of multistability by destroying most of the chimera states (leaving only the relatively rare bouncing chimeras) and thus, counterintuitively, promote synchronization. The latter can be interpreted as yet another form of AIS, where the the symmetry of the system is broken by heterogeneity in the network structure rather than among the oscillators.

V. CONCLUDING REMARKS
The significance of our study of Janus oscillator networks is multifold. First, for demonstrating the cooccurrence of an array of behaviors previously observed in disparate systems, it shows how these different behaviors are related and helps to facilitate manipulation of synchronization in applications. Second, it shows that these behaviors can be far more common than anticipated since previously assumed conditions, such as correlations between oscillator frequency and network structure for the occurrence of explosive synchronization, are in fact not restrictive. Third, it allows the systematic characterization of transitions between an unprecedented number of stable branch states representing different levels of coherence, including various types of chimera states. Fourth, it reveals new kinds of behaviors, such as inverted synchronization transitions, characterized by switch to a more (less) synchronous state in response to a decrease (increase) in coupling strength. Fifth, it demonstrates the generic occurrence of AIS in a remarkably simple system, which can facilitate experimental demonstration of this phenomenon. The latter is especially important given that this counterintuitive phenomenon was discovered recently, is only starting to be understood, and is yet to be realized experimentally.
Here we introduced Janus oscillator networks as a new class of systems that can exhibit surprising behaviors even when the network structure itself is extremely simple. Our formulation builds on the principle that it is easier to understand complexity in simpler systems, as long as they can still include the relevant features. This principle has been successfully explored in previous work that applied phase or dimension reductions [40,52] and group theory [53] to the analysis coupled oscillators. The simplicity of our model allowed us to derive analytical results and we expect that further analytic progress will be possible in future studies. Moreover, owing to this simplicity, the prospects for experiments are very positive given that coupled phase oscillators, which are the components of Janus oscillator networks, have been implemented experimentally in electric [54,55], chemical [56,57], and mechanical [58,59] networks, among others [60].
This work can be expanded in many directions. On the one hand, while here we purposely focused on simple network structures, we speculate that an even richer range of behaviors might be possible for Janus oscillators in complex network structures, which we expect to be explored in future work. On the other hand, given that we focused on Janus oscillators with only two frequencies, it would be natural to also consider higher-dimensional Janus oscillators with three or more frequencies. Another extension would be to characterize Janus oscillators whose "faces" are distinguished in terms of parameters other than the natural frequency (e.g., frustration parameter, delay, or oscillator type). In particular, experimentally accessible chemical oscillators [56,57], which are limit cycles oscillators of reactions in fluid cells, could be paired using cells with distinct geometries and chemical conditions. Furthermore, it would be natural to consider Janus oscillators in time-varying networks to study oscillators that can move in space or co-evolve with the network structure, and analyze the resulting interactions between swarming, self-organization, and synchronization [61]. Finally, we expect that this line of theoretical work will stimulate the experimental study of Janus oscillator networks and lead to yet new discoveries in the lab.