Double explosive transitions to synchronization and cooperation in intertwined dynamics and evolutionary games

Collective behavior, from murmurations to synchronized beating of heart cells, governs some of the most beautiful and important aspects of nature. Likewise, cooperation—the act of sacrificing personal benefits for the common good—is one of the pillars of social evolution, and it is the basis for the emergence of collective organized actions from single-cell organisms to modern human societies. Here we merge these two phenomena into a single model, considering an ensemble of networked oscillators, where each oscillator can be either a cooperator or a defector, and with only cooperators contributing to synchrony. At the same time, the value of the order parameter in the neighborhood of each oscillator is considered as an effective local temperature which determines the strategy updating procedure in the evolutionary game. The emergence of cooperation is thus intertwined with that of synchronization, producing a novel and fascinating dynamics which includes a double explosive transition.


Introduction
The extensive cooperation in nature, ranging from microbiology [1] to human societies [2], has profoundly changed the history of life. Nevertheless, cooperation is also known as one of the great unsolved natural mysteries-being that it is inconsistent with the fundamental Darwinian principle of evolution that organisms should act so as to maximize their own fitness-and a renewed attention to it started after Nowak and May first placed the prisoner's dilemma on a square lattice and opened up a completely new field of research [3]. In their seminal work, indeed, cooperators, despite sacrificing personal benefits for a common good, can actually survive by forming compact clusters. This phenomenon is known as network reciprocity [4], and it spawned many relevant studies in self-organization [5], dilemma strength [6,7], spatial games [8], complex networks [9][10][11][12], coevolution [13,14], to name but a few. Several recent reviews cover the latest advances in detail [2,6,7,15].
The path of studies in synchronization is quite similar. Synchronization in large groups of oscillators is ubiquitous in nature, spanning from cellular to societal levels, and with applications ranging from ensuring proper organ health and functioning to mitigating unrest [16]. In the age of network science [17][18][19][20][21][22], it became clear that both the structure and function of a network play a key role in synchronization [22][23][24][25][26][27][28][29]. Studies concentrated on the paradigmatic Kuramoto model, where oscillators are rotating on a unit circle with an arbitrary natural frequency and are coupled by means of a harmonic function of their phase differences [30][31][32][33][34]. Generally, a continuous transition from disorder to coherence occurs as the coupling strength increases, but recent studies unveiled that altering the network structure or considering adaptation can lead to explosive synchronization [25,[35][36][37][38], where the transition is instead discontinuous and may be even irreversible.
Three years ago, Antonioni and Cardillo [26] studied the coevolution of synchronization and cooperation in costly networked interactions, where a dichotomous scenario is considered in which oscillators may decide to cooperate (and pay the cost in order to be synchronized with the rest of the population), or to free-ride (without incurring in any cost, but waiting that others synchronize to its state). In this letter, we further explore this family of models, and merge synchronization and cooperation processes in a single setting in which only cooperative oscillators aim to synchronize with each other. At the same time, the local order induced by synchronization is seen as an effective local temperature which determines the strategy updating process of oscillators' strategies. Such a double feedback leads to novel and fascinating dynamics, including a double explosive transition to synchronization and cooperation as the coupling strength increases, as well as bistable steady states.

Model
Let us start by considering an ensemble of N networked oscillators, which evolve using the following equations: here, θ i and ω i are the instantaneous phase and the natural frequency of the ith oscillator, respectively, (i = 1, 2, . . . , N), and λ is the coupling strength. Furthermore, {A ij } are the elements of the network's adjacency matrix (with A ij = 1 if oscillators i and j are connected, and A ij = 0 otherwise). The degree of the oscillator i is defined as the number of its neighbors, i.e., k i = N j=1 A ij . Finally, α i is a time-dependent binary variable which assumes the value of 1 when the oscillator is in the cooperative state C, and the value of 0 when the oscillator is in the defective state D. In other words, a cooperating oscillator stays attached to the network structure and in this way it cooperates to the formation of a synchronous state, whereas a defective oscillator detaches from the network as it does not want to participate to the setting of any collective state. Notice that equation (1) reduce to the classical Kuramoto model when, at all times, α i = 1∀i, and when the network is a clique.
Each oscillator can be seen as a point in the complex plane rotating on the unit circle, and the extent of synchronization is quantified by the order parameter R, defined by R e iΨ = 1 N N j=1 e iθ j , with Ψ being the average phase. R = 0 corresponds to a fully incoherent state, where all oscillators rotate independently from each other and with different frequencies. On the other hand, R = 1 indicates a perfectly synchronized state, where all oscillators are locked to a common frequency, and all phases evolve in unison. As the connectivity structure is that of a network, for each oscillator i one can also define a local order parameter r i (measuring the extent of synchronization in node i's neighborhood) by means of the following equation: Furthermore we consider that at discrete times t 1 , t 2 , . . . , t n (with t j − t j−1 = Δt for all integers 1 < j n) each oscillator i plays an evolutionary game with pairwise interactions with all of its neighbors in the network. At each round of the game and in each one of the existing links i − j, oscillator i accumulates therefore a payoff π ij according to the following matrix: In expression (3) C and D stand, respectively, for cooperator and defector. For instance, when a C oscillator encounters a D oscillator, the former receives the quantity 0 while the latter receives b. This framework can be seen as that of a prototypical prisoner's dilemma game. In particular, b > 1 is called the temptation to defect, in that it quantifies the extra payoff that a defector obtains when meeting a cooperator.
Each time the game is updated, the oscillator i calculates its total payoff Π i = N j=1 A ij π ij , and updates its strategy by means of the Fermi rule. Namely, the oscillator i's strategy α i imitates the strategy α j (with j being a randomly chosen index among those labeling the nodes that form the neighborhood of node i) with probability: In equation (4), K stands for an overall rationality coefficient (small K values indicate rational choices, where node i tends to imitate the strategy with higher payoff, whereas large K values stand for irrational choices), while each oscillator can behave differently because the local rationality, K i = (1 − r i )K is related to the local order parameter r i . As a more physical picture, here (1 − r i ) is interpreted as a local temperature, and payoffs are seen as energy levels. Such a physical interpretation is easy to understand: if r i = 1 is in-situ equal to one, then the system is fully ordered on the spot, i.e. the local level of entropy will be the minimum and we can associate to such a frozen state a zero local temperature. On the contrary, r i = 0 would imply maximum local disorder and entropy, and one can therefore associate such a tangled state to the maximum value of a local temperature.
As we are mixing processes which are essentially different in time (the synchronization dynamics is continuous, while the evolutionary game is discrete), it is then crucial to discuss (and ultimately fix) the time scale at which the two processes need to be compared. If one takes as the unit reference for time the timescale of equation (1)

Results
Let us first report the full phase diagram in the relevant parameter space, the plane of λ and b, which are the coupling strength of synchronization and the temptation to defect, respectively, and which respectively rule the evolution of the dynamics and of the game. The results are presented in figure 1, and are obtained on an Erdös-Rényi (ER) random network [39] of size N = 8, 192 and with an average degree of k = 12. Initially, the frequencies ω i are taken from a homogenous random distribution in interval U(−1, 1). Moreover, initial phases are taken randomly, and each oscillator is initially randomly chosen as C or D. In figure 1 one notices immediately a high correlation between the patterns of R and R L . However, for λ < 0.1 (where the system does not synchronize) and for b < 1.04, one sees that R is nearly zero whereas R L is finite, and as a result ρ C is almost unchanged in this region. After the critical coupling strength value at which the synchronization emerges, it is also visible that cooperation can be maintained with a large temptation value b-in other words, cooperation is promoted by synchronization. Finally, the novel result is that the induced double transition to synchronization and cooperation is sharp and abrupt, resembling, though in a completely different context, the features of explosive synchronization [25]. This is seen in figure 2 which indeed shows that for moderate b or λ, the transition shows an explosive character.
In order to gather more information on the observed double explosive transitions, we need to investigate the role of the overall rationality coefficient K. For this purpose, we gradually vary K, and monitor in parallel the two evolutions of the system generated by two distinct sets of initial conditions. In the first set, the system is initialized in its fully incoherent state, i.e. with R(t = 0) ≈ 0, while in the second set the initial phases are fully synchronized, i.e. one has R(t = 0) ≈ 1. The asymptotic values of R and ρ C in the two simulations are denoted as R 0 , R 1 and ρ C0 , ρ C1 , respectively. Figure 3 reports, for each value of K, the average results from 50 different independent trials. Note that both ρ C and R are always decreasing as K increases. This implies that rationality helps to maintain cooperation, and therefore leads indirectly to synchronization. The remarkable result here is the presence of a large hysteresis region (the green shaded  , where ρ C0 (R 0 ) and ρ C1 (R 1 ) coexist as bi-stable solutions, and which then implies irreversibility [24,40]. The bistability is easily understandable from figure 3, where for either large (>5) or small (<1) coefficient K, the effect of the local order parameter r i on the local rationality K i is neglectable. The group rationality is then maintained in a lower or high level, and thus totally favoring or suppressing the cooperation. For a moderate K, however, if set K i to a large value initially, the cooperation is inhibited, and indirectly suppresses the synchronization. As a feedback, the rationality K i is kept high and continue to suppress cooperation. If we set a small K i initially, the mechanism is similar, but the cooperation and synchronization are supported. We also obtain the hysteresis region on scale-free network [41] and the results are presented in details in our SI.
The landscape of bistability in the space of the parameters {λ, b, K} provides a more complete view of the coevolutionary dynamics. If we determine the existence of bistability at one point with boundary d = 0.1 (there are many other possible choices as the criteria, but the selection of criteria would not affect our main results), then three regions are found from the simulation results: where only cooperation shows bistability.
• Reg.II: where both cooperation and synchronization show bistability. • Reg.III: Otherwise, no bistability is observed. The Reg.I and Reg.II in the full space are depicted in figure 4 (see the caption for the color code). It is noteworthy that the bistability of synchronization is possible only if the bistability of cooperation exists. In Reg.I, the coupling strength λ is too small that R(t = ∞) = 0 irrespectively of the initial condition, i.e., there is no synchronization. However, when the temptation b is moderate, the bistable phenomenon of  cooperation influenced by the initial group rationality can be observed. Specifically, strong rationality stimulates a high level of cooperation, while weak rationality leads to a low level of cooperation. The simultaneous bistability of cooperation and synchronization emerges together in Reg.II, where the coupling strength λ is large enough and the temptation is moderate. Obviously, equation (4) shows that for an extremely small rationality coefficient K ≈ 0, local order parameter r i almost does not change the local rationality K i . At this time, the two dynamics are dominated by the evolutionary game process, and there is no bistability. While K is larger than a threshold (e.g., K ≈ 1.2, in figure 3), then the co-evolution actually begins and shows the bistability. To sum up, the bistability is induced by the interplay of the two dynamics.

Conclusion
In summary, we have studied the emergence of collective dynamics in a general model which intertwines cooperation and synchronization. In particular, we have observed a double explosive transition for synchronization and cooperation. Namely, we have shown that by either increasing the coupling or the dilemma strength, the transition from incoherence to coherence (as well as that from full defection to full cooperation) become abrupt, irreversible, discontinuous and explosive. By further considering the effect of rationality, the transition implies a bistable state, where either a high or a low fraction of cooperators can be maintained in a stable way. The exploration in the parameter space shows the dynamics demonstrating co-evolution induces the bistability.
It is natural to ask whether the model used for the game (synchronization) could have broader implications in the formation of the phenomena reported here. We argue that it does, e.g., the payoff matrix equation (3) used in this paper is the so-called boundary game, while the pairwise game can be classified according to the strength of two dilemmas, namely GID and GAD (see [42][43][44][45]). In future work, we will expand more universal games.
Further extensions of the model remain to be considered, especially at the interface of physics and society [46], whereby synchronization represents physics, and cooperation represents society. For instance, since multilayer networks go a step further in accurately reflecting the reality of interactions in nature and society [36], it would be very important to confirm these results with more realistic network's structures. It would also be interesting to go beyond pairwise interactions, and to study whether group interactions [2,8] have any role in shaping the dynamics of these family of models where dynamics and games co-evolve.