Frequency-amplitude correlation inducing first-order phase transition in coupled oscillators

The first-order phase transitions in coupled oscillators have been widely studied because of their discontinuity and irreversibility. In previous research, the designed coupling mechanisms between each pair of oscillators can cause the first-order phase transitions occur stably. In the present study, we propose a new mechanism which requires the existence of an inversely proportional relationship between the natural frequencies and the intrinsic amplitudes in the homogeneously coupled oscillators. Based on two classical oscillator models, i.e., the Poincaré model and the Stuart–Landau model, the emergence of explosive oscillation death is independent of the frequency distributions. Our findings indicate that the first-order phase transitions can be induced by the frequency-amplitude correlation for the first time. Therefore, it provides a novel perspective to understand explosive phenomena in coupled oscillators.


Introduction
Phase transitions occur every day in nature, e.g., ice melting into water, the ferromagnetic transition [1]. They refer to the processes where the system states change apparently when one of their control variables varies minimally. One of the well accepted modern classification schemes is that the phase transitions are divided into two broad categories: the first-order ones and the second-order ones [2]. In the former category, the phase transitions involve latent heat during which the system either absorbs or releases a typically large amount of energy per volume. They manifest as discontinuous and irreversible processes, such as the processes of water boiling and the super-cooling/heating transitions. In the latter category, by contrast, the processes of the phase transitions are continuous and reversible, such as the ferromagnetic transitions and the superconducting transitions.
In the last decade, the first-order phase transitions have attracted extensive attention in the fields of nonlinear dynamics and complex networks, including three remarkable phenomena, namely, the explosive percolation (EP) [3][4][5][6], the explosive synchronization (ES) [7][8][9][10][11], and the explosive death (ED) [12][13][14]. During EP, the sudden formation of an inter-connected giant cluster in a growing network which satisfies the Achlioptas-like growth rules, i.e., the formation of the small clusters is enhanced and the formation of the giant clusters is suppressed before the transition. ES refers to an abrupt and irreversible synchronization process, i.e., the transition from the incoherent state to the synchronized state, characterized by the existence of a hysteresis loop in the evolution of the order parameters with the increasing/decreasing coupling strength subsequently between coupled oscillators. ES can emerge stably, based on some designed systems, e.g., the networks of frequency-degree correlated oscillators [7], the frequency-weighted networks of oscillators [15][16][17], and the coupled oscillators via the adaptive control of local order parameters [18,19]. ED represents a first-order transition from an oscillatory state to a dead state. During the transition, the solutions of two states in the system, i.e., the stable oscillatory state and the stable steady state, co-exist over a range of coupling strength. Taken together, these three phenomena are related to each other and have been hypothesized to play a key role in the onset of the acoustical signal transduction in the cochlea [20], the cascading failure of the power grids [21], and the information spreading in the social network [22]. This article focus on the ED phenomenon in coupled oscillators.
As a type of ED, the explosive oscillation death (EOD) was first observed in coupled frequency-weighted Stuart-Landau (SL) oscillators [12], where all oscillators settle on different fixed points abruptly from disordered oscillatory states. Through numerical simulations, three types of EOD phenomena in the system emerge, corresponding to the different frequency distributions (FDs). Subsequently, the explosive amplitude death, as another type of ED, emerged in an ensemble of N indirectly coupled and identical SL oscillators through a common dynamic environment [13], indicating that all oscillators settle on the same fixed point suddenly. Moreover, ED also occurs under the other coupling mechanisms, such as mean-field diffusion [23], conjugate variables coupling [24], quorum sensing [25], and state-space-dependent dynamic interactions [26].
In the human body, the vital movements are usually regulated by the single or multiple biological oscillators, with varied frequencies and amplitudes [27,28]. The different frequencies indicate different functional consequences, while the different amplitudes abstract the differences in the expression levels of proteins, hormones, etc. A good example is the suprachiasmatic nucleus (SCN) in the mammalian brain which regulates circadian rhythms. It has been found that the neurons differ in both the intrinsic periods (frequencies) and the intrinsic amplitudes due to their different locations [29,30]. These differences have also been reported in the SCN neurons of different age groups. Specially, both the periods and the amplitudes are reduced in the SCN neurons of older mice [31,32]. In addition, the marine mammals and birds have the capacity to sleep unihemispherically [33,34]. During the unihemispheric slow-wave sleep, the sleeping hemisphere exhibits the synchronized high-amplitude and low-frequency rhythm, while the other hemisphere (i.e., awake hemi-sphere) exhibits the incoherent low-amplitude and high-frequency rhythm. This phenomenon can be explained by the chimera states in nonlinear dynamics, which have also been shown to be associated with the explosive phenomena [35,36]. In the cerebral cortex, the amplitudes of faster rhythm activities are usually lower than those of the slower components, in conformance with an empirical inverse power relation in normal electrocorticograms [37][38][39]. In this article, we introduce a correlation between the intrinsic amplitudes and the natural frequencies in an ensemble of coupled limit-cycle oscillators, and seek the possibility of the emergence of explosive transition from the oscillatory state to the dead state. Through the numerical simulation and the theoretical analysis, EOD emerges when an inversely proportional relationship between the intrinsic amplitudes and the natural frequencies exists in the homogeneously coupled oscillators. Therefore, our work provides a new perspective to understand the phenomenon of hemi-brain sleep in bio-rhythm.

Model description
The Poincaré model is a classical oscillator model which can mimic circadian rhythm neurons, due to its parameters containing the intrinsic amplitude, the natural frequency, and the stability with respect to amplitude perturbations (amplitude relaxation rate) [40][41][42]. Here, we consider the coupled Poincaré model in a fully connected network, as follows, x k , where A j is the intrinsic amplitude of the jth node, ω j is the natural frequency, r j is the coupled amplitude with r j = x 2 j + y 2 j , γ is relaxation parameters with γ = 1, and g is the global coupling strength, respectively.
In order to construct the amplitude-frequency correlation, there are two typical cases, i.e., A j = ω j and A j = 1 w j , as follows, (2) The absolute of ω j are necessary to ensure the oscillation of each Poincaré oscillator. In the next section, we show a comparison between these two cases. In addition, to avoid A j approaching infinity especially when ω j is close to zero, the upper bound of A i are set to be 100 for the convenience of calculation.
To characterize the collective behavior of coupled oscillators, two types of order parameter are defined, where θ j represents the phase of the jth oscillator, which is defined as θ j = arctan y i x i . The order parameter R (0 R 1) characterizes the phase coherence of the system, which does not involve any information of amplitude. Another order parameter R r (R r 0) characterizes the coherence of the dynamics with both the amplitude and the phase.
In this article, the size of network is N = 500 and the fourth-order Runge-Kutta method is employed for numerical simulations with time increments of 0.01 h. In order to avoid the effect of transients on the coherence quantification, the first 2000 000 time steps, i.e., 20 000 h, are neglected. The coupling strength g changes adiabatically by Δg = 0.01 and the order parameter R for each value of g is calculated. Without loss of generality, we also consider the SL oscillator model and present the results in the supplementary materials (https://stacks.iop.org/NJP/24/073038/mmedia) (figure S1).

Numerical results
The case of A j = 1 |ω j | is presented at first. Figure 1 shows the simulation results in the coupled Poincaré oscillators where FDs satisfy U(−0.5, 0.5) and N(0, 0.5), respectively. For U(−0.5, 0.5), the phase transition in the forward direction, i.e., from the incoherent state to the coherent state, occurs when g is equal to 0.3 (g = g 1 = 0.24), while the phase transition in the backward direction, i.e., from the coherent state to the incoherent state, occurs when g = g 2 = 0.18. For N(0, 0.5), the forward phase transition occurs when g = g 1 = 0.3, while the backward phase transition occurs when g = g 2 = 0.23. The similar phase transition processes occur for both R and R r . Thereby, for the FDs satisfying ω j = 0, the first-order phase transitions occur with a hysteresis loop during the evolution of the order parameters R and R r , respectively. In figure 2, the simulation results are given in the oscillators with FDs satisfying ω j = 0, including N(0.5, 0.5) and R(0.2), respectively. The left corresponds to N(0.5, 0.5), and the phase transition occurs when g = g 1 = 0.38 and g = g 2 = 0.49 in the forward and backward directions, respectively. The right column of figure 2 corresponds to R(0.2), and the two critical points do not coincide (g 1 = 0.29 and g 2 = 0.34). Clearly, the explosive phase transitions also emerge in the case of ω j = 0. We repeat the calculation 10 times for each FD to verify the stability of the numerical simulation. Although the width of hysteresis loop and the position of transition points may vary slightly, the explosive transition occurs for each simulation. These results show that our mechanism can stably generate the explosion phenomenon.
From above, we find that the emergence of explosive transitions are independent of the symmetry of FDs in our mechanism. However, there is a question whether the inversely proportional relationship between A j and ω j , i.e., A j = 1 |ω j | is necessary. Instead of the inversely proportional case, we choose the positive proportional case, i.e., A j = |ω j |, to simulate. Figure 3 shows the transitions of R in the same four FDs, respectively. Apparently, the system exhibits the continuous phase transitions instead of the first-order phase transitions. It is difficult for all oscillators to be coherent, especially for the Gaussian FDs in figures 3(b) and (c). These results indicate that the inversely proportional relationship between the intrinsic amplitudes and the natural frequencies of networked oscillators plays a significant role in inducing the emergence of explosive phase transition.
Subsequently, the microstates of the oscillators in the case of A j = 1 |ω j | are observed. In figure 4, the orbits of all oscillators, near the critical point during the forward and backward transitions, are depicted. The oscillators of FD obeying N(0, 0.5) are shown as an example. The panels (a) and (b) corresponds to the  processes of forward transition. Each oscillator rotates on its natural frequency before the critical point. However, when g transits the forward critical point, all oscillations suddenly cease at different fixed points. The same process also occurs in the backward transition as shown in (c) and (d), but the critical point is different. These phenomena indicate that EOD appears in the system. As shown in (b) and (c), all oscillators form two phase-symmetric clusters on the coordinate plane when they are OD, even though their  radii of rotation are different. Few oscillators do not belong to any clusters because their intrinsic amplitudes are cut off during the calculation as described in section 2.
With varying g, the frequencies and amplitudes of coupled oscillators will change due to the coupled interactions. To quantify how the coupled amplitudes and frequencies change respectively, we define two quantitative measures, (1) effective frequency ratio, P j ω , (2) effective amplitude ratio, P j a , where j = 1, 2, . . . , N, as follows, whereθ j represents the instantaneous frequency of the jth oscillator. P j ω is taken as the ratio of the averaged instantaneous frequency to the natural frequency. Similarly, P j a can be regarded as the ratio of the averaged coupled amplitude to the intrinsic amplitude. Here, the parameter T means the calculation steps, which is defined as 200 000. Figures 5 and 6 show the evolution of P j ω and P j a of all oscillators, respectively. Taking the FD obeying U(−0.5, 0.5) as an example, all P j ω retain 1 before the transition and then they suddenly jump to 0, as shown in figures 5(a) and (c). A similar phenomenon also appears for all P j a , they also retain 1 before transition, and then suddenly enlarge, as shown in figures 6(a) and (c). These results indicate that before all oscillators are dead, they still maintain their natural frequencies and intrinsic amplitudes.
The evolution processes of P j ω and P j a of all oscillators corresponding to the case of A j = |ω j | are presented in the supplementary materials (figures S2 and S3). When A j = ω j , P j ω and P j a will continuously change as g increases, that means the coupled frequency and amplitude of jth oscillator differs for the different values of g.

Theoretical analysis
By polar coordinate transformation: x j = r j cos θ j , and y j = r j sin θ j , we can rewrite equation (1) aṡ r n cos θ n cos θ j , r n sin θ n sin θ j .
We set a reference rotating frame, i.e., ψ(t) = ψ(0) + Ω t, where Ω is the average effective frequency of oscillators in the system. According to the definition of the order parameter R r e iψ 2 = 1 N N j=1 r j e iθ j ,  (7) can be rewritten after employing the mean field theorẏ Define η j = θ j − ψ. Thus equation (8) can be reduced tȯ r j = r j (A j − r j ) + gR r cos η j , The necessary condition of phase-locking between two frequency-weighted Kuramoto oscillators i and j is described as [11], Equation (10) is taken as a suppressive rule for the formation of synchronized local clusters in the explosive transition, which can be considered as a dynamical counterpart of the Achlioptas process. When all oscillators are attracted by the synchronized cluster, the further increase of coupling strength induces the explosive transition (i.e., the significant jumping of the order parameter).
Motivated by this idea, we propose the suppressive rule for the coupled Poincaré oscillators. For any two oscillators i and j, substituting the phase difference α ij = θ j − θ i = η j − η i into equation (7), we obtain When the two oscillators are phase-locked, letα ij = 0. Hence, equation (11) becomes Because of sin η i ∈ [−1, 1] and sin η j ∈ [−1, 1], we obtain Equation (13) is the necessary condition of phase-locking between any two Poincaré oscillators. Equation (13) agrees well with equation (10), if r i,j = 1 |ω i,j | . In the following analysis, we suppose r j = 1 |ω j | for any j. When all oscillators are phase-locked, we haveη j = 0. The stationary state η * j = 0 can be solved from equation (9) as After solving equation (14), we obtain two symmetrical solutions, i.e., η * a = arcsin( 1 gR r ) and η * b = arcsin( −1 gR r ), in the cases of U(−0.5, 0.5), N(0, 0.5) and N(0.5, 0.5), respectively. However, for the case of R (0.2), there is only one solution, i.e., η * a = arcsin( 1 gR r ), because all ω j are positive. Based on this analysis, we conclude that the effect of amplitudes prevents the dispersion of synchronized clusters if the coupled oscillators satisfy r i,j = 1 |ω i,j | . From the perspective of the phase plane, excluding the case of Rayleigh distribution, all oscillators will evenly split into two clusters that are symmetric with respect to the real axis in their phase-locked state (shown in figure 4(b)). In particular, the radii r i of all oscillators in each cluster are different. With further increasing g, they both approach the real axis gradually. Only in the limit case of g, i.e., g → ∞, the two clusters combine into one. Therefore, we can regard the whole system as two clusters a and b that are located at (r p a , η * a ) and (r q b , η * b ), rotating with frequency ω > 0 and −ω < 0, respectively. Here, the cluster a is composed of N 1 oscillators with amplitude r p a , p = 1, 2, . . . , N 1 , and cluster b is composed of N 2 oscillators with amplitude r q b , q = 1, 2, . . . , N 2 . We let N 1 + N 2 = N. According to the definition of order parameter R r , we have Not only for Rayleigh FD, but also for other FDs, cos η * = 1 − 1 (g 2 R 2 ) always exists due to η * = arcsin( ±1 gR r ), we have Therefore, two symmetrical solutions with respect to the order parameter R r can be obtained from equation (16), as Since R r is real, equation (17) remains if and only if: then we obtain the critical point g c = . In other words, when the coupling strength adiabatically reduces to g c , all oscillators will suddenly start to rotate. Since the amplitude must be real values, when g < g c , all amplitudes have steady solutions, i.e., r j = A j , which have been verified by numerical simulation that shown in figure 6. Hence, we can get the solution of critical point for the backward transition, Accordingly, we obtain the analytical solutions, i.e., g = 0.1987, 0.2322, 0.3864, 0.2878, corresponding to the above four FDs, respectively, which are quantitatively consistent with figures 1 and 2. Moreover, the equidistant distribution is taken into consideration to further verify our analytical results. Let ω j be in an ascending order from −δ to δ with N = 200 without ω j = 0. Namely, the difference between each frequency is fixed at δ 100 . The simulation results are deterministic since that the interference from random frequencies is excluded. Figure 7(a) shows the numerical results of the backward transition process when δ changes from 0.2 to 0.9. Figure 7(b) shows the comparison between numerical simulation and analytical solution. As the width of the FD δ increases regularly, the average value of A j = 1 ω j decreases regularly, so the backward phase transition point, g c =

Discussion
In this article, we have proposed a new mechanism to induce the first-order phase transition in coupled oscillators through the frequency-amplitude correlation. Based on an inversely proportional relationship between the intrinsic amplitudes and the natural frequencies of networked oscillators, EOD can emerge stably. We have obtained four highlighted results. First, the mechanism designs the correlation between the natural frequencies and the intrinsic amplitudes for oscillators, which has not been considered in the previous models. Second, the mechanism requires the heterogeneity of the intrinsic amplitudes but not the heterogeneity of the couplings. The first-order phase transitions have been found in the homogeneous coupled Kuramoto oscillators of a limited system size based on particular FDs [43]. However, it is the first time to study EOD in the homogeneous coupled Poincaré oscillators or SL oscillators. Hence, our mechanism relaxes the conditions for the emergence of EOD. Third, the transition point based our mechanism can vary depending on the average amplitude of the system. We have shown another case in the supplementary materials (figure S4), in which A j = −αω j + c. By controlling the slope coefficient α, the type of phase transition changes from discontinuity to continuity. It means that we can control the emergence of explosive transition by modulating the intrinsic amplitudes of oscillators. Last but not least, our mechanism for the frequency-amplitude correlation can induce EOD for the uniform and Rayleigh FDs, which cannot be induced by the frequency-weighted model [44][45][46]. The comparison diagrams are shown in the supplementary materials (figures S5 and S6).
The amplitude effect of coupled oscillators was a broad topic, which induced various phenomena observed in typical theoretical models, such as oscillation quenching [47,48], aging transition [49,50] and amplitude chimera states [51,52]. As an important property of oscillation, the amplitudes that were thought to be related to the frequencies, strongly affect the collective behaviors of coupled oscillators. Ermentrout studied the behavior for the rings of coupled oscillators, and found that the amplitudes of the system would be changed drastically, when the waves are present [53]. Miroll and Strogatz analyzed a large system of the mean-field coupled limit-cycle oscillators with random FDs, and determined the stable region of the coupling variance space for the amplitude death [54]. Based on two coupled Rossler units, Li and Zheng found that the amplitude played a significant role in the phase synchronous dynamics [55]. Here, we have proposed that the heterogeneous intrinsic amplitudes that are related to the natural frequencies can affect the type of phase transition. This finding have advanced the understanding of amplitude effects in coupled oscillators [56][57][58][59]. In addition, the amplitude modulation of biological or chemical oscillators around the bifurcation points was considered to be significant [27,60,61]. For the Brusselator model, a forcing term was added into each oscillator equation for their amplitudes and frequencies modulated in the periodic oscillation state [62][63][64]. In order to ensure the existence of an approximately inverse ratio between the amplitude and the frequency within each individual oscillator, we have designed the special forcing terms. We show an explosive phase transition in the forced Brusselators in the supplementary materials (figures S7 and S8).
The application of neurotoxin (TTX) results in the disperse phases of the SCN cells that show autonomous oscillations through the inhibition of intercellular coupling [65]. After washout of the TTX, the network of SCN neurons was firstly desynchronized and then synchronized. In other words, the coupling between SCN neurons was first inhibited (g decreasing) until there was no coupling (g = 0), and then the coupling was restored (g increasing). We analyzed the data, which the authors provided, and indeed observed the first-order phase transition. Since the application of TTX affects intracellular coupling, we infer the existence of an explosive phenomenon (figures S9 and S10, see the supplementary materials in details). Considering that the intrinsic oscillation amplitude and oscillation period of SCN neurons are highly heterogeneous [66,67], we infer that the correlation of neuronal intrinsic amplitude and frequency may play an important role in this first-order phase transition.
Explosive transitions have been mainly investigated in pair-wised coupling, and the recent studies have shown that it is possible to achieve explosive transitions in higher-order networks, such as hypergraphs and simplicial complexes [19,[68][69][70][71]. The first-order phase transitions have also been discovered in experiments [72,73]. Thus we look forward that our works can be studied in higher-order systems and verified by experiments. In summary, we provide a new mechanism for inducing the first-order phase transitions from the perspective of intrinsic properties' correlation.