Bifurcation and scaling at the aging transition boundary in globally coupled excitable and oscillatory units

Following a previous paper [Phys. Rev. E 88, 052907 (2013)], we study in detail the mechanism of aging transition in globally and diffusively coupled excitable and oscillatory units. Here two of the three models taken up in the earlier work are used, each composed of a large number of units with their bifurcation parameters forming a uniform distribution. The control parameters are the coupling strength and the average of the bifurcation parameters. The present work is mostly devoted to a region of the phase diagram near the aging transition boundary with the coupling strength greater than its critical value for the onset of bistability and hysteresis. The bifurcation structure of each system at the aging transition boundary is investigated theoretically as well as numerically. Moreover, we show that critical scaling laws of order parameters $S$ and $M$ used in the previous paper are different depending on which region of the coupling strength to be chosen.

Figure 1

Phase diagrams of the phase oscillator model, Eq. (1). (a) $\gamma =0.5,N=2000$: DP1 (upper left side, red), DP2 (upper right side, green). The broken curve (blue) is a theoretical boundary at which the order parameter $R$ vanishes. (b) $\gamma =0.5,N=3200$: Relayed initial conditions are used, which means using the final point in phase space computed for the last value of $a_{\text{av}}$ as the initial condition for its next value with $K$ fixed. The real curve (label A, red) and the broken curve with shorter line segments (label B, blue) are the boundaries between SP and DP for growing and decreasing $a_{\text{av}}$, respectively. The broken curve with longer line segments (green) and the dotted curve (magenta) are the boundaries between DP1 and DP2 found for $a_{\text{av}}$ changing in the same directions as above. The boundaries A and B are drawn by detecting the start of oscillation of the $N\mathrm{th}$ unit, because it has the largest value of $a$ in the ensemble and the system is in SP if it does not oscillate. The criterion of oscillation is that the magnitude of its average frequency is larger than $\pi /t_{\text{cp}}$, where $t_{\text{cp}}$ is the computational time after a transient. The boundaries separating DP1 and DP2 are drawn under the same criterion as in Ref. [26]; i.e., if $M>\sqrt{2/N}$, then the system is considered to lie in DP2.

Figure 2

Phase diagram for the phase-oscillator model beyond and near $K=K_{\text{oc}}\sim 0.90 (\gamma =0.5,N=3200,6400)$. The regions labeled “SP, DP1” and “SP, DP2” are where both phases coexist. The squares (magenta) and asterisks (blue) locate the upper and lower ends of the floating branch. The solid squares (orange) show the boundary between DP1 and DP2, part of which for $K<K_{\text{f}}$ is the same as the dotted curve (magenta) in Fig. 1. The two black curves are theoretical boundaries below (above) which a solution belonging to DP1 (SP) emerges for the first time as $a_{\text{av}}$ is increased (decreased). The agreement between theory and simulation is excellent. Other details are the same as Fig. 1.

Figure 3

Behaviors of the order parameters $S, M$, and $R$ in the phase-oscillator model $(\gamma =0.5,N=6400)$. Relayed initial conditions are used (see the caption of Fig. 1). The data displayed (asterisks, red; squares, green) are those for increasing and decreasing $a_{\text{av}}$, respectively, while solid circles (blue) show the floating branch with its neighborhoods included. The vertical lines in each panel show where discontinuous transitions take place, except for the rightmost vertical line of panel (d). (a) $K=1.13$; (b) $K=1.148$ (immediately before $K=K_{\text{f}}$); (c) $K=1.149$ (immediately after $K=K_{\text{f}}$); (d) $K=1.16(>K_{\text{h}})$. The floating branches in (c) and (d) were found using the final values of all variables around a roughly middle point of the same branch for a previous value of $K$ as the initial condition.

Figure 4

Critical behavior of the order parameter $S$ in the phase-oscillator model $(\gamma =0.5,N=6400)$. (a) $K=0.85(<K_{\text{oc}})$. The symbols ($+$, red) are simulation results. The curves (blue, green) are due to theory. (b) $K=1.05(>K_{\text{oc}})$. The curves (magenta, blue) are due to theory. Panel (b) displays the simulation data obtained for changing $a_{\text{av}}$ with relayed initial conditions (upward, $+$, red; downward, $\times$, green) (see the caption of Fig. 1). See the point at which the lower theoretical curve with $S>0$ (blue) intersects the horizontal axis, whose abscissa is $a_{\text{c}}$. The right (red) and left (green) vertical lines correspond to boundaries A and B, respectively, in Fig. 1.

Figure 5

Behavior of $\left|Z\right|$ in the phase oscillator model $(\gamma =0.5,N=6400)$, where the simulation data show the moduli of time averages of $Z$. (a) $K=0.85(<K_{\text{oc}})$. The simulation results: $+$, red. The curves are theoretical: the left one (blue) and the right one (green) are subcritical $(R=0)$ and supercritical $(R>0)$ branches, respectively. (b) $K=1.05(>K_{\text{oc}})$. The simulation results are obtained by increasing ($+$, red) and decreasing ($\times$, green) $a_{\text{av}}$ with relayed initial conditions (see the caption of Fig. 1). The curves are theoretical: the magenta (mostly on the left side) and the blue one (mostly on the right side) correspond to the branches of $R=0$ and $R>0$, respectively.

Figure 6

Mechanism of the branching of dynamic solutions corresponding to DP2 in the phase-oscillator model $(\gamma =0.5,N=6400)$, where the ordinate is $\left|Z\right|$. Simulation results have been obtained using relayed initial conditions (see the caption of Fig. 1). The simulation data in panels (a) and (b) are the moduli of $Z$'s time averages. (a) $K=1.149$ (immediately after $K=K_{\text{f}}$). Simulation results: data for increasing $a_{\text{av}}$ ($+$, red) and for decreasing $a_{\text{av}}$ ($\times$, green) and the floating branch with its neighborhoods included (asterisk, blue). The curves are due to theory: The upper cyan and lower magenta curves correspond to SP and DP1, respectively. The discontinuous transitions from DP2 to SP occur at $a_{\text{av}}\sim 0.9604$ and $\sim 0.9631$. (b) $K=1.16(>K_{\text{h}})$. Simulation results: data for increasing $a_{\text{av}}$ ($+$, red) and the floating branch with its neighborhoods included (asterisk, blue). The curves are due to theory: The upper cyan and lower magenta curves correspond to SP and DP1, respectively. Note that in each of the above two panels, the point at which the floating branch deviates from the theoretical DP1 branch divides the former branch into DP1 and DP2. (c) Behavior of $Z$ on the complex plane: $K=1.149,a_{\text{av}}=0.9645$. Simulation results: the closed curve (red, DP2) and the point ($+$, green, SP). Theory: three points ($\times$, blue, SP, stable), (asterisk, magenta, SP, unstable), and (square, cyan, DP1, unstable). (d) Behavior of $Z$ on the complex plane: $K=1.149,a_{\text{av}}=0.9635$. Details are the same as panel (c).

Figure 7

Trajectories of $Z$ on the complex plane obtained from the phase-oscillator model $(\gamma =0.5,N=6400)$. The smallest (red), middle (green), and largest (blue) “closed” curves are for $(K,a_{\text{av}})=(1.05,0.958),(1.13,0.965)$, and $(1.2,0.97)$, respectively.

Figure 8

Average frequency vs $j$ for the phase oscillator model $(\gamma =0.5,N=6400)$. The parameters $K$ and $a_{\text{av}}$ of the data labeled “a”, “b”, and “c” are the same as those of the smallest, middle, and largest “curves” in Fig. 7, respectively.

Figure 9

Phase diagrams of the simplified oscillator model, Eq. (3). (a) $\gamma =0.02,N=1000$: DP1 (upper left side, red), DP2 (upper right side, green) [26]. The broken curve (blue) is a theoretical boundary at which the order parameter $R$ vanishes. (b) $\gamma =0.02,N=4000$: Relayed initial conditions are used (see the caption of Fig. 1). The meanings of the line styles (colors) and the labels are the same as in Fig. 1. The boundaries A and B are drawn using the criterion that the system is in DP if $S>0.02$. Moreover, the boundaries separating DP1 and DP2 are determined by the same criterion: The system is in DP2 if $M>0.02$. Note that the value 0.02 is not very different from $\sqrt{1/N}(=0.0158...\text{for}N=4000)$ used in drawing panel (a) [26].

Figure 10

Phase diagram near $K=K_{\text{f}}\sim 0.355$ for the simplified oscillator model $(\gamma =0.02,N=4000)$. Details are the same as in Fig. 2. The criteria to identify the boundaries are the same as in Fig. 9. The upper theoretical boundary (black curve) is drawn using Eq. (34).

Figure 11

Behaviors of the order parameters $S, M$, and $R$ in the simplified oscillator model $(\gamma =0.02,N=4000)$. Relayed initial conditions are used (see the caption of Fig. 1). (a) $K=0.354$; (b) $K=0.354934$ (immediately before $K=K_{\text{f}}$); (c) $K=0.36(<K_{\text{h}}\sim 0.374)$; (d) $K=0.377(>K_{\text{h}})$. Other details are the same as in Fig. 3.

Figure 12

Behavior of $X$ in the simplified oscillator model $(\gamma =0.02,N=4000)$. The simulation results here are time averages of $X$. (a) $K=0.15(<K_{\text{oc}}=0.2)$ (see the formula of $K_{\text{oc}}$ in the text and Ref. [26]). The simulation results ($+$, red) exist in the whole region of ${\mu }_{\text{av}}$, which are perfectly covered by the theoretical results that consist of two parts: ${\mu }_{\text{av}}<{\mu }_{\text{c}}\sim -0.005$, where $R=0$ (SP, green curve), and ${\mu }_{\text{av}}>{\mu }_{\text{c}}$, where $R>0$ (DP1, blue curve). (b) $K=0.35$. The branches connected with the two vertical lines are the simulation results which have been obtained by increasing ($+$, red) and decreasing ($\times$, green) ${\mu }_{\text{av}}$ with relayed initial conditions (see the caption of Fig. 1). There are also two theoretical branches: Part of the lower branch (SP, $R=0$, blue curve) with positive slopes perfectly overlaps with the simulation result, while the upper branch (DP1, $R>0$, magenta curve) initially agrees with the branch obtained by simulation ($\times$, green), but starts to disagree after that, which is because the system enters DP2 there.

Figure 13

Bifurcation diagram and trajectories in the simplified oscillator model $(\gamma =0.02,K=0.356>K_{\text{f}})$. Simulation results $(N=4000)$ have been obtained using relayed initial conditions (see the caption of Fig. 1). (a) Simulation results, which are time averages of $X$: data for increasing ${\mu }_{\text{av}}$ ($+$, red) and for decreasing ${\mu }_{\text{av}}$ ($\times$, green) and the floating branch with its neighborhoods included (asterisks, blue). The curves are based on theory: The upper black and lower cyan curves respectively correspond to DP1 and SP. The discontinuous transitions from DP2 to SP happen at ${\mu }_{\text{av}}\sim -0.00506$ and $\sim -0.00373$. (b) Trajectories on the $(X,Y)$ plane, where ${\mu }_{\text{av}}=-0.0051458333...$ and $Y$ is defined as ${\sum }_{j=1}^N{x}_j^2/N$. The orbit (red) is due to simulation, while the point ($+$, green) corresponds to DP1 (unstable, theory). Calculations of $Y$ in SP and DP1 are easily performed on the basis of the theory developed in the text (see the Appendix).

Figure 14

Time series of $X$ in the simplified oscillator model $(\gamma =0.02,N=4000)$. The smallest (red), middle (green), and largest (blue) amplitude curves are the data for $(K,{\mu }_{\text{av}})=(0.3,-0.001),(0.33,-0.0005)$, and $(0.377,-0.0005)$, respectively.

Figure 15

Average frequency vs $j$ for the simplified oscillator model $(\gamma =0.02,N=4000)$. The parameters $K$ and ${\mu }_{\text{av}}$ of the data labeled “a”, “b”, and “c” are the same as those of the smallest, middle, and largest amplitude curves in Fig. 14, respectively.

Figure 16

Critical scaling behavior of $S$ in the region $0\le K<K_{\text{oc}}$. The straight line in each panel shows the slope of $3/4$. (a) The phase-oscillator model with $\gamma =0.5$ and $K=0.5$. The system size $N$: 1600 ($+$, red), 6400 ($\times$, green), 25 600 (asterisks, blue), 102 400 (squares, magenta). (b) The simplified oscillator model with $\gamma =0.02$ and $K=0.12$. The system size $N$: 4000 ($+$, red), 16 000 ($\times$, green), 64 000 (asterisks, blue). The straight line here is drawn using the asymptotic formula of $S$, namely, Eq. (36) combined with Eq. (10) of Ref. [26].

Figure 17

Critical scaling behavior of $S$ at $K=K_{\text{oc}}$. The straight line in each panel shows the slope of $1/2$. (a) The phase-oscillator model with $\gamma =0.5$. The system size $N$: 4000 ($+$, red), 16 000 ($\times$, green), 64 000 (asterisks, blue), and 256 000 (squares, magenta). (b) The simplified oscillator model with $\gamma =0.02$. The system size $N$: 8000 ($+$, red), 32 000 ($\times$, green), 128 000 (asterisks, blue), and 512 000 (squares, magenta). The straight line here is drawn in the same way as in Fig. 16.

Figure 18

Critical scaling behaviors of $S$ and $M$ at $K=K_{\text{h}}$. The symbols: $S$ ($+$, red) and $M$ ($\times$, green). The straight lines in each panel show the slope of $1/5$. (a) The phase-oscillator model with $\gamma =0.5$ and $N=25600$. Here $K=K_{\text{h}}=1.156275$. (b) The simplified oscillator model with $\gamma =0.02$ and $N=16000$. In this case, $K=K_{\text{h}}=0.37375$.

Figure 19

Critical scaling behaviors of $S$ and $M$ for $K>K_{\text{h}}$. The symbols: $S$ ($+$, red) and $M$ ($\times$, green). The straight line in each panel shows the slope of $1/4$. (a) The phase-oscillator model with $\gamma =0.5,K=1.19$, and $N=12800$. (b) The simplified oscillator model with $\gamma =0.02,K=0.5$, and $N=16000$.

Figure 20

Behavior of $Z$ for $K$ beyond and near $K_{\text{h}}$ in the phase-oscillator model $(\gamma =0.5,N=6400,K=1.16>K_{\text{h}})$. (a) A trajectory and fixed points on the complex plane. Simulation results: the orbit (red, $a_{\text{av}}=0.96575$), stable SP points; ($+$, green, $a_{\text{av}}=0.96525$) and ($\times$, blue, $a_{\text{av}}=0.9655$). Note that the latter SP points overlap SP points obtained by theory for the same parameter values (see below). Theory: two pairs of SP points (asterisks, magenta, $a_{\text{av}}=0.96525$; squares, cyan, $a_{\text{av}}=0.9655$), where in each pair, the left (right) one is unstable (stable), and the other is an unstable DP1 point (solid square, orange, $a_{\text{av}}=0.96575$). The validity of theory is confirmed by the excellent agreement between its SP points and corresponding simulation results. (b) Time series of $\text{Re}Z$ (lower curve, red) and $\text{Im}Z$ (upper curve, green), where the same data as the closed orbit in panel (a) are used $(a_{\text{av}}=0.96575)$.