Collective dynamics of heterogeneously and nonlinearly coupled phase oscillators

Coupled oscillators have been used to study synchronization in a wide range of social, biological, and physical systems, including pedestrian-induced bridge resonances, coordinated lighting up of firefly swarms, and enhanced output peak intensity in synchronizing laser arrays. Here we advance this subject by studying a variant of the Kuramoto model, where the coupling between the phase oscillators is heterogeneous and nonlinear. In particular, the quenched disorder in the coupling strength and the intrinsic frequencies are correlated, and the coupling itself depends on the amplitude of the mean-field of the system. We show that the interplay of these factors leads to a fascinatingly rich collective dynamics, including explosive synchronization transitions, hybrid transitions with hysteresis absence, abrupt irreversible desynchronization transitions, and tiered phase transitions with or without a vanishing onset. We develop an analytical treatment that enables us to determine the observed equilibrium states of the system, as well as to explore their asymptotic stability at various levels. Our research thus provides theoretical foundations for a number of self-organized phenomena that may be responsible for the emergence of collective rhythms in complex systems.


I. INTRODUCTION
The Kuramoto model is a well-known and widely applicable paradigm for theoretically describing and modeling collective behavior in large ensembles of interacting units [1]. Since its conception in 1975, the model has been used to explore synchronization dynamics in numerous social, biological, and physical systems, with examples ranging from arrays of Josephson junctions [2], flashing fireflies [3], cardiac pacemaker cells [4], to economic markets [5]. The model has been particularly popular in physics and mathematics, because it often affords analytical insights towards better understanding different types of collective behavior and synchronization transitions that occur in systems consisting of a large number of coupled nonlinear self-sustained oscillators [6][7][8][9].
In its original version, the Kuramoto model elucidates synchronization at the onset of a phase transition that emerges due to the interplay of the intrinsic frequencies of individual oscillators and the global coupling among them. Specifically, the heterogeneity of the natural frequencies of oscillators, which are distributed randomly across the population according to a prescribed probability density, tends to desynchronize the system, while the attracting global coupling between the oscillators opposes this disorder and tends towards synchrony. These two competing forces result in the rich dynamics of the Kuramoto model that has been explored in great detail during the past decades [10].
In addition to the heterogeneity in natural frequencies, recent work in physics and network science has highlighted the importance of heterogeneous coupling between dynamical units. Namely, the coupling strength between individual oscillators need not be uniform, but rather oscillator-dependent [11,12]. Such heterogeneities can be encoded by means of quenched random interactions, which often represent a realistic setup for a wide range of systems and applications that go beyond the traditional identical coupling [13][14][15][16]. For example, such heterogeneous interactions have been used to shed light on the collective dynamics of excitatory-inhibitory neurons, or to describe opinion formation amongst conformists and contrarians in sociology [17][18][19]. In terms of theoretical results obtained with heterogeneous coupling between phase oscillators, previous works have reported the emergence of traveling waves, standing waves, π-states, glass states, cluster states, as well as chimera states [20][21][22][23][24].
Almost in parallel to coherent phases observed in net-work dynamics that emerge from heterogeneous coupling, explosive synchronization transitions have also received a lot of attention [25][26][27][28]. The latter represent first-orderlike phase transitions that switch abruptly between incoherent and coherent states, and they are typically due to correlations in network topology and oscillator dynamics [29][30][31]. Since abrupt transitions between two stable states underlie many real-life systems, the phenomenon has received ample attention, for example to better understand bistable perception in the brain dynamics [32] and the robustness and optimal functioning of large-scale power grids [33]. More recently, it has also been shown that explosive synchronization transitions can be due to higher-order interactions, where a link connects more than two oscillators, or due to nonlinear coupling, where, for example, the coupling strength between oscillators is proportional to the coherence of the system [34][35][36][37][38]. Research has also shown that these two effects can actually be considered equivalent [39][40][41]. But despite the wealth of findings and insights related to the Kuramoto model and its variants, an unexplored question is what novel effects could different coupling patterns have on the macroscopic dynamics of the system. Particularly so if both heterogeneity and nonlinearity are present in the coupling among phase oscillators. In this case, indeed little is known about the possible collective behaviors and the phase transitions between them that may be observed as a result. In this paper, we therefore consider a new variant of the Kuramoto model, which incorporates heterogeneity and nonlinearity in the coupling. More specifically, we consider correlations between the quenched disorder in the coupling strength and the intrinsic frequencies, and we also consider the coupling to be dependent on the amplitude of the mean-filed of the population. As we will show, the interplay between these effects significantly shapes the overall collective dynamics of the system, giving rise to a number of fascinating phenomena of relevance for the emergence of synchronization. Importantly, the studied model is analytically tractable, so that our research provides theoretical foundations, and in fact a natural mechanism, for the spontaneous emergence of various rhythmic states and the bifurcations associated with them in complex systems, and it does so with relatively small changes to system parameters.
The remainder of this paper is organized as follows. In Sec. II, we present the new variant of the Kuramoto model, and we establish the self-consistency approach to determine the long-term macroscopic dynamics. In Sec. III, we carry out a linear stability analysis for the fully locked, two-cluster, states in a special onedimensional manifold, and we set up a necessary stability criterion for their occurrence. In Sec. IV, we extend this stability result to all perturbed directions in the phase space, and we explore the characteristic dynamics of the system due to perturbations. In Sec. V, we investigate the asymptotic stability of the steady states in the infinite limit, and we obtain the associated eigenspectrum by using the Ott-Antonsen reduction. In Sec. VI, we investigate the transitions by means of which the system can move between different coherent states. Various types of phase transition towards synchronization are identified, including the explosive synchronization transition, the abrupt irreversible desynchronization transition, and the tiered synchronization transition with or without a vanishing onset. In Sec. VII, we summarize our main conclusions and briefly discuss their possible wider implications.

II. MATHEMATICAL MODEL AND ITS STATIONARY SOLUTIONS
As noted, we consider an extension of the Kuramoto model, whose governing equations arė Here θ i is the phase of i−th oscillator with i = 1, ..., N . {ω i } are the natural frequencies selected from a probability density function g(ω), which is assumed to be symmetric throughout the paper, i.e., the mean frequency is set to zero and g(ω) = g(−ω). Eq. (1) differs from the classical Kuramoto model in that the uniform coupling strength has been replaced by the heterogeneous coupling K i with a feedback factor f (R), which is a generic function of the amplitude of the Kuramoto order parameter defined by R(t) measures the level of coherence of the system and Θ(t) denotes the average phase of the population. K i and f (R) account for the inhomogeneity and nonlinearity of the coupling, respectively. Without loss of generality, we set f (R) > 0 (which can be made by rescaling the time) and K i = K N −i . As we will show below, the long-term dynamics of this model is significantly richer than the dynamics of the traditional Kuramoto mode with linear homogeneous coupling. We emphasize that the inherent nonlinearity, quenched random interactions, and the large number of degrees of freedom make it difficult to understand the quantitative dynamics of Eq. (1) with arbitrary choices of K i and f (R). We therefore consider only the case where the coupling is symmetric around index i of each oscillator. Aside from this constraint, however, we are free to choose K i and f (R) arbitrarily.
Before proceeding the analysis, we make a few comments of the model. As mentioned above, the heterogeneous interactions have been used to shed light on the intrinsic dynamics of glass states, or to describe opinion formation amongst conformists and contrarians in sociology. On the other hand, the natural frequencies, as a typical heterogeneity of the population, are the intrinsic properties of oscillators themselves. Therefore, it is reasonable to establish correlations between the coupling strength and natural frequencies [12,46,58]. In all the quoted studies, the coupling strength depends on the natural frequencies implicitly, and the aim is to uncover the traveling waves and glass states in coupled oscillator populations.
In our model, instead, the coupling strength depends on the natural frequencies explicitly, i.e., we set K i = K|ω i | corresponding to a frequency-weighted coupling (FWC). The FWC was proposed by Ref. [42], and then it was generalized in Ref. [43,45]. In all those studies, it was shown that a FWC is equivalent to the frequency-degree correlations in complex networks reported in Ref. [29]. Both models turn out to give rise to Explosive synchronization phenomena [59].
Regarding the nonlinear coupling (NC), where f (R) is a nonlinear function with respect to R. Relevant examples include Ref. [51,60], where NC can be achieved by means of phase reduction in a system of Stuart-Landau oscillators. For convenience, one may set f (R) = R β−1 (a power-law function), which serves as a typical NC. Similar discussions can be found in Ref. [52,55]. At the same time, recent works in physics and neuroscience highlight the potential importance of higher-order (or group) interactions, e.g., three-way or more, connections that may be organized via higher-order simplexes or a simplicial complex. The effect of such interactions on the dynamics thus represents an important topic of research in the complex system community. More importantly, recent research has also shown that the higher-order interactions are equivalent to NC [41,61].
The definition of the order parameter Eq. (2) allows us to rewrite Eq. (1) aṡ By doing so, Eq. (3) retains the mean-field character, and K|ω i |R β displays an effective force acting on each phase oscillator θ i . Next, we briefly outline the self-consistency approach to analyze the stationary dynamics of Eq. (1). The key idea for the self-consistent method is to assume that, in the long-term evolution, the macroscopic quantity R(t) is a constant, and Θ(t) rotates uniformly, i.e., Θ(t) = Ωt + Θ 0 . Indeed, we may set Ω = Θ 0 = 0 by entering into an appropriate rotating frame and shifting initial conditions [62,63]. Then, Eq. (3) can be reformulated aṡ where the auxiliary parameter q = KR β is defined to simplify notation. We can see from Eq. (4) that each individual oscillator exhibits two types of long-term behavior depending on the magnitude of q. Specifically, if q > 1, all oscillators become phase-locked with and When q < 1, all oscillators are drifting, i.e., they rotate nonuniformly on the unit circle and form a stationary distribution given by The sign function sgn(ω) appearing in Eq. (7) is needed to ensure the positivity of the distribution. In contrast to the traditional Kuramoto model and its previously studied variants, we find that in our case the coexistence of the phase-locked and the drifting populations can never occur. In other words, the system remains either in the fully phase-locked or the totally drifting state for a fixed value of q. The splitting of the system into the fully phase-locked or the completely drifting state allows us to obtain an expression for the order parameter R self-consistently and where · denotes the average over the population.
Remarkably, the symmetry of the system implies that sin θ q>1 = e iθ q<1 = 0. Therefore, for q < 1, the only solution to Eq. (8) is R = q = 0, which corresponds to the incoherent state with the phases scattered uniformly around the unit circle. Conveniently, for q > 1, the expression for R becomes simply The derivations presented above reveal three important results. In the first place, the presence of correlations between frequency and coupling results in a universal phase transition towards synchronization, whose equilibrium is independent of the frequency distribution and system size. Secondly, there are no solutions to Eq. (11) for a sufficiently small value of K, since F (q) remains finite for q ∈ (1, +∞). And thirdly, there exists a critical coupling K c , beyond which the two coherent state solutions exist located at q < q c and q > q c , respectively. The critical parameter q c satisfies the condition F (q c ) = 0, i.e., which in turn determines the critical coupling K c and the critical order parameter R c , yielding and In what follows, we will prove that the solutions for R > R c and q > q c are attractive, and moreover, correspond exactly to the solutions along the F (q) < 0 branch.

III. STABILITY IN THE ONE-DIMENSIONAL MANIFOLD
We have shown by means of Eqs. (5) and (6) that, depending on the sign of the natural frequencies and depending on the intensity of q, the oscillator population splits into two populations in the long run at opposite locations. These observations lead us to consider the following question: What is the stability property of the steady states described by Eq. (11) in this manifold? As we will show, this consideration is quite convenient, because the low-dimensional description of the dynamics in Eq. (1) allows us to deduce the associated critical points and gain some intuition about the stability of phaselocked states in the presently studied modified Kuramoto model.
In the fully locked scenario, the system acts like two giant oscillators, labeled as θ ± , which correspond to natural frequencies ±ω, respectively. In such a manifold, Eq. (1) reduces tȯ Without loss of generality, we assume that ω > 0, thus having sin θ ± = ±q −1 and cos θ ± = 1 − q −2 , and the order parameter becoming Since θ + + θ − = 0, Eqs. (15) and (16) become onedimensional if we introduce the phase difference φ = θ + − θ − = 2θ + , which then yieldṡ with In this one-dimensional manifold, the dynamics is relatively simple to analyze. The steady state solution of Eq. (18) requires that which implies that the coupling strength K ≥ K c , such that We can then show that and that the critical order parameter R c is given by which are exactly the expressions obtained previously in Eqs. (13) and (14). In order to consider the stability of the solutions beyond K c , we impose a small perturbation δφ away from φ in Eq. (20) and neglect higher order terms of δφ. In doing so, we obtain from which it follows that the characteristic value of the perturbation in this direction is proportional to −Φ (φ). This means further that the stable (unstable) condition To corroborate this conclusion, we note that which can be simplified significantly by taking into account Eqs. (17) and (20), thus obtaining From this we find that R > R c leads to Φ (φ) > 0, thereby determining a stable branch of the solution, while R < R c leads to Φ (φ) < 0, corresponding to an unstable solution.
Based on the above treatment, we conclude that the steady state solution along the branch q > q c (or F (q) < 0) is stable at least in this special perturbed direction. However, the stationary solution along the branch q < q c (or F (q) > 0) is unstable in all perturbed directions, because it has a positive characteristic value of the perturbation even in the synchronized manifold.
In the subsequent two sections, we will show that these stability criteria can be extended to any perturbation direction, and this regardless of whether N is finite or infinite.

IV. STABILITY IN THE N -DIMENSIONAL MANIFOLD
It is worth pointing out that what is unusual about the currently studied system is that the dynamics of Eq. (4) naturally entrains the oscillators into two groups with equal numbers and with opposite positions. To better understand the nature of these fixed points, we restrict our attention to their asymptotic stability in the N -dimensional manifold. We will show that the solutions obtained with R > R c or q > q c are indeed stable in all the perturbed directions, while the solutions obtained with R < R c or q < q c are not.
For the locked states (θ 1 , θ 2 ...θ N ), let θ i (t) = θ i + x i (t) and x = (x 1 , x 2 ...x N ) be the small perturbations. The linearization of the dynamics described by Eq. (1) is given by the following equation [64] x = Jx, where J is the Jacobian matrix, whose entries are We note that the rotational invariance of the dynamics described by Eq. (1) implies that J always has a trivial eigenvalue λ = 0 corresponding to an eigenvector (1, 1...1). To verify that this indeed holds, we point out that the rows of J all sum to 0, i.e., where we have taken into account N j=1 sin θ j = 0 and N j=1 cos θ j = N R. Aside from λ = 0, the stable con-dition for the locked states requires that the remaining N − 1 eigenvalues are all negative. Before a more in-depth exploration of the eigenvalues of J, we here provide an intuitive argument. We first express where c = (c 1 , c 2 ...c N ) with c i = cos θ i and s = (s 1 , s 2 ...s N ) with s i = sin θ i . Hence, the stability condition is equivalent to Jxx < 0 for all x ∈ R N with x = 0, and x is orthogonal to (1, 1...1). That is, N i=1 x i = 0, and we thus have Following Ref. [64], it becomes apparent that J has at most one positive eigenvalue. If there are two or more positive eigenvalues for J, the span of the corresponding eigenvectors defines at least a two-dimensional subspace, on which Jxx > 0. On the other hand, one can find a nonzero vector x on a space that is orthogonal to s, but from the identity above, we have Jxx < 0, which thus leads to a contradiction.
To explore the eigenvalues of J more conveniently and in further detail, we introduce a change of coordinates in order to partially diagonalize the system. In particular, if y = (y 1 , y 2 ...y N ) with y i = |ω i |c i x i , we obtain The vectors u = (u 1 , , are introduced to ease notation, and their product is In this new coordinate system, we can express with a and b being the projections to the unit vector u u and v v , respectively. y ⊥ is the remaining vector, such that y ⊥ · u = y ⊥ · v = 0. We thus have In particular, for b = 0, the term around ab vanishes, and the identity above becomes a quadratic form. Since the coefficient a could be arbitrarily small, the stability condition Jxx < 0 requires that which again recovers the necessary stability condition obtained in Eqs. (22) and (23).
To proceed with extracting more details about the eigenvalues of J, we now compute its characteristic polynomial. We therefore express J in the matrix form where W and D are the diagonal matrices with the entries |ω i | and qc i along the diagonal, respectively. The matrix C is given by The characteristic polynomial is The key task from here on is to determine the second term above. We recall that the rank of C is only two, and det (D + λW −1 ) −1 = 0. The same is true for the matrix (D + λW −1 ) −1 C. Inspired by this property, we choose an orthogonal basis c c , s s and extend it into R N . The remaining N − 2 vectors span the kernel of (D+λW −1 ) −1 C. Therefore, the matrix (D+λW −1 ) −1 C if restricted to the orthogonal basis ( c c , s s ) reduces to a 2 × 2 matrix given by where its entries are defined as We note further that the terms Q 12 (λ) and Q 21 (λ) are in fact all zero, because the natural frequencies are grouped into plus-minus pairs, i.e., ω i = −ω N −i . This symmetrical case significantly simplifies the characteristic polynomial into The eigenvalues λ corresponding to P (λ) = 0 are therefore determined by and For convenience, we order the natural frequencies so that |ω 1 | ≤ |ω 2 | ≤ ... ≤ |ω N |. The obtained eigenvalues merit three important comments. In the first place, we note that λ = 0 is a trivial solution of Eq. (46) because K − 1 β = H c (0) is precisely the self-consistent equation described by Eq. (11). As already mentioned, this property stems from the rotational symmetry of the dynamical equation. Secondly, if λ < −|ω N | q 2 − 1, both functions H c (λ) and H s (λ) are negative. This implies that the characteristic Eqs. (46)- (47) have no roots in the region λ ∈ (−∞, −|ω N | q 2 − 1). And lastly, lim H c(s) (λ) = ±∞, which reflects that, away from their poles, the functions H c (λ) and H s (λ) are strictly decreasing. Accordingly, the characteristic equations must each have exactly one root between any two consecutive poles (−|ω i+1 | q 2 − 1, −|ω i | q 2 − 1).
Since functions H c (λ) and H s (λ) are decreasing for λ ∈ (−|ω 1 | q 2 − 1, +∞), and lim λ→+∞ H c(s) (λ) = 0, the remaining one eigenvalue is completely determined by the difference ∆ defined by Specifically, ∆ < 0 implies that there is a negative root to Eq. (47) in the region λ ∈ (−|ω 1 | q 2 − 1, 0). On the other hand, for ∆ > 0, a positive root to Eq. (47) exists in the region λ ∈ (0, +∞). By combining this with the self-consistency equation, we obtain the elegant identity ∆ = βqF (q), (49) which directly relates the stability criterion for the locked states to the shape of the self-consistent equation F (q). Namely, if F (q) > 0 the steady state solutions are unstable, whereas if F (q) < 0 the steady state solutions are attractive and stable, which completes the proof.

V. STABILITY IN THE INFINITE LIMIT
Next, to better understand the dynamics that emerges from the above considerations, we study the system in the thermodynamical limit. In particular, we investigate two aspects of the asymptotic stability of the equilibrium states, for which we reduce the system to a lowdimensional equation that governs the long-term macroscopic dynamics. Under the limit N → ∞, the dynamical system is equivalent to the continuity equation for the probability density function ρ(θ, ω, t) Here the velocity v(θ, ω, t) is given by and the complex order parameter Z(t) becomes By applying the Ott-Antonsen ansatz [65,66], the probability density function can be expressed as where the variable α(ω, t) is a complex valued function with |α(ω, t)| ≤ 1, and the bar denotes the complex conjugate. This consideration underlines the Poisson kernel that is dynamically invariant, as long as α(ω, t) evolves according to Meanwhile, Z(t) reduces to where we have defined the operatorĝ to denote the integral with respect to ω. In this setting, the quantity α(ω, t) can be interpreted as the local order parameter that is formed by the oscillators with frequencies near ω.
With respect to Eq. (54), we are interested in the steady statesα = 0, which yields two-composed solutions We note that the first branch of Eq. (56) corresponds to the drifting population, whereas the second branch corresponds to the phase-locked case. Other possible steady state solutions of α 0 (ω) are ruled out by the constraints |α 0 (ω)| ≤ 1 and R ≥ 0 [67].

VI. NUMERICAL RESULTS
To show how the heterogeneous and nonlinear coupling affects the phase coherence, we study the system numerically. Figs. 1 and 2 illustrate the macroscopical and microscopical synchronization properties, respectively. When the exponent β is fixed, and the natural frequencies are drawn independently from a symmetric distribution g(ω), we can indeed observe rich collective dynamics from the still relatively simple extension of the traditional Kuramoto model. In Fig. 1, we show the time-averaged amplitude of the order parameter, R = 1 T t+T t R(t)dt with T being the time-averaged window, as a function of the global coupling strength K. In the simulations, both directions for synchronization transitions are considered, i.e., K is first increased adiabatically from 0 and then decreased. The presented simulation results reveal that the inhomogeneity and nonlinearity in the coupling combines together to affect the overall collective dynamics of the system, leading to a number of interesting dynamical phenomena.
In particular, for β = 1 (left column), the system displays three types of phase transitions to synchronization depending on the difference between K f and K c . Specifically, if K f > K c [ Fig. 1(a1)], the system exhibits a firstorder (explosive) phase transition to synchronization. As the system transitions from the completely incoherent state ( R ∼ 0) to the fully locked state ( R ∼ 1) at the critical coupling K f under forward continuation, when K is gradually decreased, another abrupt transition from synchronization to incoherence occurs at K c . In the hysteretic region, K ∈ (K c , K f ), the system admits a bistability, where the stable incoherent state and synchronized state coexist. If K f = K c [ Fig. 1(a2)], the system undergoes a hybrid phase transition, which is characterized by a vanishing hysteresis loop. Interestingly, if K f < K c [ Fig. 1(a3)], a tiered phase transition to synchronization occurs, where the system converts from the incoherent state to the fully synchronized state, which is mediated by an oscillatory state that emerges when K ∈ (K f , K c ).
When the exponent β = 1, the phase diagrams change considerably. For β < 1 (middle column), we have shown that K f = 0 < K c , so that the system displays a tiered phase transition to synchronization with a vanishing onset, where the incoherent state is replaced by an oscillatory state the amplitude of the mean-field is time-varying that arises in the region K ∈ (0, K c ). However, for β > 1 (right column), K f = +∞ > K c , an irreversible abrupt desynchronization transition is observed, featuring a hysteresis area, whose width is infinite. The system experiences a discontinuous transition from the fully locked state to the incoherent state at K c in the backward continuation, where R jumps from R c to 0 directly. In contrast, there is no corresponding counterpart of this abrupt transition when going from the incoherent to the coherent state.
In Fig. 2, we highlight the oscillatory states that emerge for the intermediate region of coupling strengths [70]. We show the instantaneous phases, effective frequencies, and the time series of R(t). These results illustrate that, in the oscillatory state, the oscillators form coexisting synchronized and drifting groups, or spontaneously partition into different clusters according to their natural frequencies. Moreover, the effective frequencies remain entrained in each cluster, but they do not synchronize with one another. As a result, the instantaneous phases show a correlation at a fixed time, thereby causing the order parameter to behave periodically. For further details, we refer to the results shown in Fig. 2 and the corresponding description of the parameter values in the figure caption.

VII. DISCUSSION
In summary, we have studied an extension of the Kuramoto model, where the randomness of the coupling strength is related to the natural frequencies of the coupled oscillators, and the mean-field is modified to depend on the global coherence of the system. In doing so, we have thus jointly considered heterogeneity and nonlinearity in the coupling of the Kuramoto model, which has led us to uncover a number of fascinating dynamical states and transitions associated with their emergence. These include coherent and incoherent states, two-cluster states, and oscillatory states, which can emerge via various transitions among them, including explosive synchronization transitions, hybrid transitions with hysteresis absence, abrupt irreversible desynchronization transitions, and tiered phase transitions with or without a vanishing onset.
Despite the wealth of different dynamical states and the transitions between them, our model still proved to be analytically tractable, thus allowing us to provide firm theoretical foundations for the observed collective behavior. In particular, we have provided a rigorous stability analysis of the equilibrium states by studying their eigenspecturm properties in the finite and infinite size limit. Moreover, we have shown that the critical points that correspond to bifurcations of steady states and the stable conditions for their occurrence in the phase space can be obtained by means of the matrix theory and the Ott-Antonsen reduction.
Although the Kuramoto model and its many variants have been studied extensively in the past decades, we have shown that the joint consideration of heterogeneity and nonlinearity in the coupling among oscillators affords new insights into the emergent collective dynamics that could be of relevance for complex systems where such conditions apply. Previous research has already highlighted the importance of heterogeneous interactions for the dynamics of excitatory-inhibitory neurons, and for the emergence of consensus amongst conformists and contrarians in different social settings [17][18][19]. Heterogeneity and nonlinearity are also often related to plasticity, which leads to changes in the connection strength among individuals to achieve efficient global states. This may have implications in neuroscience as well as in the so-cial sciences, where such interaction scenarios often unfold. More generally, our research reveals how an alternative coupling scheme can lead to the emergence of diverse rhythmic states in complex systems that exhibit synchronization transitions with relatively small changes to system parameters. Since such complexity is rarely analytically tractable, we hope that our model will prove to be inspirational for future research along these lines and find applicability also in realistic complex systems.