Frustration induced transient chaos, fractal and riddled basins in coupled limit cycle oscillators

TL;DR

This paper demonstrates that frustration in coupled limit cycle oscillators induces transient chaos and complex basin structures like fractal and riddled basins, revealing new insights into multistability and chaotic dynamics in simple oscillator arrays.

Contribution

It uncovers the role of frustration as a source of transient chaos and complex basin structures in coupled limit cycle oscillators, including phenomena like super persistent transient chaos.

Findings

Transient chaotic dynamics occur near chaotic itinerancy regions.

Complex basin structures such as fractal and riddled basins are observed.

Transient chaos persists even in oscillation death states.

Abstract

This report unravels frustration as a source of transient chaotic dynamics even in a simple array of coupled limit cycle oscillators. The transient chaotic dynamics along with the multistable nature of frustrated systems facilitates the existence of complex basin structures such as fractal and riddled basins. In particular, we report the emergence of transient chaotic dynamics around the chaotic itinerancy region, where the basin of attraction near the chaotic region is riddled and it becomes fractal basin when we move away from this region. Normally the complex basin structures are observed only for oscillatory states but surprisingly it is also observed even for oscillation death states. Interestingly super persistent transient chaos is also shown to emerge in the coupled limit cycle oscillators, which manifests as stable chaos for larger networks.

Equations

$${\dot{x}}_i$$

$${\dot{y}}_i$$

$$x_1^{\ast }=-x_2^{\ast }=\pm \sqrt{\frac{\alpha {ϵ}_1-1+{ϵ}_1{ϵ}_2}{\alpha {ϵ}_1}},y_1^{\ast }=-y_2^{\ast }={ϵ}_1{x}_1^{\ast }.$$

$${\mathrm{B}}_1:$$

$${\mathrm{C}}_1:$$

$${\dot{\eta }}_{1j}$$

$${\dot{\eta }}_{2j}$$

$$x^{\ast }=\pm \frac{\sqrt{1+2(p-1)\alpha {ϵ}_1-4(p-1{)}^2{ϵ}_1{ϵ}_2}}{\sqrt{2}\sqrt{(p-1)\alpha {ϵ}_1}},y^{\ast }=2{ϵ}_1(p-1)x^{\ast }$$

$$\displaystyle\lambda_{1,2}=\frac{1}{2}\bigg{(}d_{1}\pm\sqrt{d_{2}}\bigg{)},\quad\lambda_{3,4}=\frac{1}{2}\bigg{(}d_{1}+2(p-1)(\epsilon_{1}-\epsilon_{2})\pm\sqrt{d_{3}}\bigg{)}$$

$$p(n)=\sum _{j=1}^{n}x(j)\cos (jc),q(n)=\sum _{j=1}^{n}x(j)\sin (jc),n=1,2,3,...,N.$$

$$M(n)=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{k=1}^N{\left[p(k+n)-p(k)\right]}^2+{\left[q(k+n)-q(k)\right]}^2.$$

$$⟨\tau ⟩=\nu (|{ϵ}_1-{ϵ}_1^{\ast }|{)}^{-\gamma }.$$

Full text

††journal: Communications in Nonlinear Science and Numerical Simulation

1 Introduction

Frustrated systems have been receiving a great deal of interest over several years as they have become a source of multiple ground states and/or multistability because of their inability to simultaneously minimize the competing interaction energies between the interacting constituents [1]. As a classical example, one can consider the case of three Ising spins with anti-ferromagnetic interaction placed at the corner of an equilateral triangle, where the inability to align each spin anti-parallel with each other results in frustration and gives rise to multiple ground states [1, 2, 3]. Further, the frustrated systems facilitate the onset of metastable states when they tend to relax themselves from frustration [4, 5]. Such frustration induced multistable and metastable states are the pivotal factors for the underlying rich dynamics of biological networks such as neural networks, genetic networks, slime models [2, 4, 5, 6, 7, 8] and so on. For instance, one of the metastable states, namely, frustrated chaos (or chaotic itinerancy) observed in such systems can explain stimulus independent switching between the various cognitive states observed in the brain and is central to the information processing ability of the latter [4, 5, 9, 10]. The multistable nature results in multiple operating regimes in biological systems facilitating functional flexibility in response to external stimuli [6].

Despite the emergence of a rich variety of dynamical behaviors due to multistable and metastable dynamics observed in frustrated systems, frustration has not yet been reported as a source of finite time chaotic behavior, namely the transient chaotic states. Such a transient chaotic behavior characterized by the combined features of chaotic and convergent dynamics produces complex basin structures, namely fractal and riddled basins due to the interplay with multistable nature of the system [11, 12, 13, 14, 15]. It was shown that the rapid transitions among different ordered states observed during cognitive tasks or memory recall are mediated by chaotic transients [16, 18, 19, 20]. Due to these reasons, transient chaotic associative neural networks are found to be useful in the field of neurocomputing and others. For instance, these networks are useful for solving combinatorial optimization problems [21, 22], progressive memory recalling, temporal pattern associations [23], maintenance-scheduling problem of generators in a practical power system and assigning channels in a cellular mobile network [24, 25]. In view of these applications, including combinatorial optimization problem, it has been shown that the transient chaotic behavior prevents the system to be trapped from a local minimum and provide a quick search over the fractal space [21, 22]. The information about the input conditions should fade after converging to a particular ordered state in order to process the next input stimulus and hence the converging dynamics of transient chaos facilitates fading memory. Furthermore, the searching space is reduced due to fractal basin structure which favors an efficient searching for the required ordered or memory state [21, 22]. Due to these reasons, this type of transient chaotic networks combining chaotic and convergent dynamics have high ability to realize an efficient search for a variety of optimization problems concerned with neuro-computing. The efficiency of such networks has been studied from both theoretical and experimental levels. Experimental implementation of the scheme can be seen in [26] where neuro-computing system is based on amoeboid unicellular organism, namely the true slime mold Physarum. This system is found to show chaotic itinerant or wandering dynamics where spontaneous transition among the multiple quasi-stable modes even without any external perturbation can be seen. It allows one to reach a desirable orbit or attractor quickly by suitable external forcing. The results show the high capability of the system to search the optimal solution of the four-city traveling salesman problem with a high probability. Particularly, information processing based on edge of chaos and nonlinear transient computing are promising contemporary research topics both from computational neuroscience and dynamical systems points of view [19, 27, 28, 29, 30, 31], where the edge of chaos that lies between ordered and chaotic regimes is proved to be the optimal region for computations in the cutting edge computing paradigm called liquid state machines (reservoir computing) [32].

The transient chaotic behavior has been mainly achieved in the literature by controlling chaos [33, 34] or by chaotic anealing [20, 21, 22, 35, 36], while the associated complex basin structures have been largely reported in driven dissipative systems, chaotic systems, delay coupled systems and maps [37, 38, 39, 40, 41]. From the experimental realization and application point of view, simple non-chaotic systems exhibiting transient chaos along with the associated complex basin structures are preferable as any experimental error and associated noise will not get amplified across the network, due to non-chaoticity, which may result in more accurate and robust performance. Hence, in this report we explore the possibility of generating transient chaos and complex basin structures (including fractal and riddled basins) in a simple system of coupled limit cycle oscillators by exploiting frustration dynamics. In addition, this report will also elucidate the possibility to inherit the chaotic wandering dynamics in a simple system of frustrated limit cycle oscillators, where such a chaotic wandering or itinerant dynamics is found to facilitate the spontaneous recollection of associated memory states in accomplishing a task.

In particular, in this report, we demonstrate frustration induced transient chaos and multistability in a paradigmatic model of coupled van der Pol oscillators exhibiting limit cycle oscillations with both attractive and repulsive couplings. Importantly, the competing attractive and repulsive couplings plays a crucial role in many biological processes. For instance, coexistence of excitatory and inhibitory synaptic couplings can be found in a pair of neurons. In the gene regulatory network the counteracting positive and negative feedback loops are used to operate various functions including bistable switches, oscillators and excitable devices [42, 43, 44]. Hence, we also aim towards understanding dynamical behavior due to the competing effect of attractive and repulsive couplings. In the non-frustrated case of two coupled systems ($N=2$), the system exhibits interesting multistabilities and the frustration comes into play when $N$ is increased to $3$ which results in rich dynamics in the multistable region of anti-phase oscillations with oscillation death observed for $N=2$. Interestingly, the frustration that exists among inhomogeneous oscillations give rise to chaotic itinerancy and transient chaotic behaviors. The transient chaotic behavior in combination with the multistability that arises due to the inhomogeneous oscillatory and oscillation death states leads to riddled and fractal basin structures. In addition, we also observe super persisting transient chaos (i.e. the chaotic behavior persists for very long time) in the region where the attractive and repulsive couplings are almost of equal strengths. This type of super persisting transient chaos has been first reported theoretically by C. Grebogi et al. in [45]. With the coupled chaotic electronic circuits, the first experimental evidence of long transient chaos is reported in [46] and their existence has also been reported in a practical example of Tilt-A-Whirl [47]. The super persistent transient chaotic state is stabilized to stable chaos with the increase of network size ($N>9$) giving rise to interesting multistability between completely ordered state and disordered state. Further, we report the emergence of fractal basins in the oscillation death region for the first time in the literature. We have also extended our analysis to a network of $N=100$ oscillators in search of extreme multistability with rich dynamical behavior.

The rest of the paper is organized as follows: In Sec. 2, we introduce our model of globally coupled van der Pol oscillators and the corresponding dynamical behavior will be discussed from smaller network to larger networks. The observed results are summarized in Sec. 3. The additional information about transient chaotic nature, Hopf bifurcation curves, mean exit time of transient chaos induced crisis, dynamics in riddled and fractal basin regions, uncertainty exponents and the transition from transient chaotic nature to stable chaos while increasing size of the networks will be discussed in the Appendix A.

2 Dynamical behavior in coupled vdP oscillators

2.1 Model

We consider a simple paradigmatic model of the van der Pol oscillator (vdP) which is known to be a primary model for describing self-excited oscillations observed extensively in biological, physical and engineering systems. The FitzHugh-Nagumo model is a modified form of the vdP oscillator [48]. From a practical point of view, van der Pol oscillators can be implemented using appropriate electronic circuit [49]. The set of coupled van der Pol oscillators is represented as

[TABLE]

where, $x_{i}$, $y_{i}$ are the state variables of the $i^{\mathrm{th}}$ system and $\alpha$ is the strength of the linear (negative) and the nonlinear (positive) damping terms, while $\epsilon_{1}$ and $\epsilon_{2}$ denote the strengths of the attractive and the repulsive couplings, respectively. Throughout this report, $\alpha$ is taken to be $0.5$. In Eq. (1), ${\overline{x}}=\displaystyle{\frac{1}{N}\sum\limits_{i=1}^{N}{x}_{i}}$, ${\overline{y}}=\displaystyle{\frac{1}{N}\sum\limits_{i=1}^{N}{y}_{i}}$. The dynamical behavior are obtained by solving the system of equation (1) using Runge Kutta $4^{\mathrm{th}}$ order integration method with step size $h=0.01$. We also wish to note that there is no need to restrict $h=0.01$, and one may consider other suitable step sizes, for example $h=0.005$ or $0.05$, for the numerical integration and the results are reproducible for those values of $h$ too.

2.2 Non-frustrated case (N=2)

A minimal network that can facilitate the frustration dynamics is a network of three coupled oscillators. Prior to unraveling the dynamical transitions in a three coupled network of vdP oscillators, we will briefly summarize the dynamics of two coupled vdP oscillators. When $N=2$, the attractive interaction tends to align the coupled oscillators to exhibit in-phase oscillations, while the repulsive coupling has the tendency to align them to evolve in anti-phase with each other. The in-phase synchronized state (IPS) is stable for appreciably large attractive coupling ($\epsilon_{1}$) whereas the anti-phase synchronized (APS) state is predominant for a large repulsive coupling ($\epsilon_{2}$) as is evident from Fig. 1(a). Multistability between APS and IPS states is also observed in the region $R_{1}$ (see Fig. 1(a)) for smaller values of $\epsilon_{1}$ and $\epsilon_{2}$ as a result of the trade-off between the attractive and repulsive couplings. Oscillation death is also found to occur in the two coupled system for higher values of $\epsilon_{1}$ and $\epsilon_{2}$ where we have observed the stabilization of the system towards one of the inhomogeneous steady states,

[TABLE]

Linear stability analysis elucidates that this anti-symmetric state is stable in the region $\epsilon_{1}>\frac{1}{\epsilon_{2}}$ when $\epsilon_{2}<1.0$ and in the region $\epsilon_{1}>1.0$ when $\epsilon_{2}\geq 1.0$ while the stabilization of the inhomogeneous (OD) states occur through a sub-critical Hopf bifurcation. Also the OD state coexists with the APS and IPS states to give rise to bistability in the regions $R_{2}$ and $R_{3}$, respectively, while multistability among them arises in $R_{4}$. The analytical stability curve for OD state is shown in Fig. 1(a). Also, the boundary of in-phase synchronized state is found through master stability analysis (see Fig. 1(a)).

2.3 Frustrated case (N=3)

On adding one more element to the system of two coupled vdP oscillators, We find the attractive coupling has the tendency to force all the three oscillators to exhibit in-phase oscillations whereas the phase-repulsive coupling facilitates the three oscillators to align in anti-phase with each other in order to minimize the energy. But this sort of simultaneous alignment among the three oscillators showing anti-phase nature to each other is practically impossible and as a consequence geometric frustration comes into play. Thus in the presence of the repulsive coupling, instead of anti-phase alignment (with $\pi$ phase difference), the three oscillators are aligned in an out-of-phase synchronized state (OPS) with $\frac{2\pi}{3}$ phase difference which is also called as the rotating wave [50]. The different phase synchronization including tri-phase synchronization are also found in a frustrated system of calling frogs [51]. The stable regions of the IPS and OPS states in a system of three coupled vdP oscillators are depicted in Fig. 1(b). Further, the OD state observed in the three coupled oscillators is no more anti-symmetric ($x_{i}\neq-x_{j}$, $i,j\in\{1,2,3\}$ and $i\neq j$) and it takes either of the forms given by $B$ and $C$ below.

OD State - $B$:

[TABLE]

OD State - $C$:

[TABLE]

where, $z=\frac{5}{3}-\frac{2^{\frac{1}{3}}p}{3(q+\sqrt{q^{2}+4p^{3}})^{\frac{1}{3}}}+\frac{(q+\sqrt{q^{2}+4p^{3}})^{\frac{1}{3}}}{3\times 2^{\frac{1}{3}}}$ with $p=-7+12a+12b$, $q=-20-18a-18b$, $a=\epsilon_{1}\alpha$ and $b=\epsilon_{1}\epsilon_{2}$. Depending on the initial conditions, the system takes one of the steady state configurations given by $B$ or $C$. The subscripts of $B$ and $C$ in Eqs. (3)-(4) denote the number of oscillators that can populate the corresponding steady state configuration. In particular, if $B$ is the asymptotic state of the system then as $t\rightarrow\infty$, one among the three oscillators reach the state $B_{1}$ while the other two oscillators populate the state $B_{2}$. Considering the commutation symmetry ($x_{i}\leftrightarrow x_{j}$), of the system, totally six possible stable steady state configurations can be observed in the gray region (marked as OD) in Fig. 1(b).

Comparing the stable regions of IPS, OPS and OD states in the case of three coupled oscillators with that of the two coupled oscillators case in Fig. 1(a), significant changes and novel dynamics can be observed in the multistable regions of OD and APS states (that is, corresponding to regions $R_{3}$ and $R_{4}$ in Fig. 1(a)). Figure 1(b) shows that in the latter parametric range, the multistability between OPS and OD states is largely suppressed with the emergence of rich dynamical states including three periodic OPS (OPS(3P)), inhomogeneous oscillatory, and chaotic states. The stable regions of these states are also shown clearly in Fig. 1(c). The dynamical regimes between the two coupled and three coupled oscillators are almost identical for weaker attractive and repulsive couplings. In such parametric ranges, we observe transition from OPS state to IPS state through a multistable region of OPS and IPS (region $S_{0}$) as a function of $\epsilon_{1}$. Increasing the strength of $\epsilon_{1}$ further leads to the stabilization of OD state and a bifurcation analysis elucidates that it is stabilized via a subcritical Hopf bifurcation (as in the case of $N=2$) for all values of $\epsilon_{2}$ that lie below the line $L_{1}$ in Fig. 1(b).

2.4 Super persistent transient chaos

Surprisingly stronger repulsive coupling gives rise to the onset of chaotic transients and other non-trivial dynamics in specific parametric ranges (in the regions where APS-OD coexists in the case of two coupled vdP oscillators) in between the lines $L_{2}$ and $L_{3}$ in Fig. 1(b) for $\epsilon_{2}=1.0$ and $1.7$, respectively. First we will discuss the dynamical transitions of the system as a function of $\epsilon_{1}$ for $\epsilon_{2}=1.0$, which are illustrated in the bifurcation diagram in Fig. 2 (using XPPAUT software). It elucidates that for smaller values of $\epsilon_{1}$ ($\epsilon_{1}<0.85$), OPS state is the only stable state (stable nature of this state is represented by lines connected by filled (brown) squares). By increasing the value of $\epsilon_{1}$, the IPS state (whose stable nature is represented by lines connected by filled (green) circles) also becomes stable at the left boundary of region $S_{0}$. Thus in the region $S_{0}$, both OPS and IPS states are stable. On further increasing $\epsilon_{1}$ the OPS state gets destabilized ( which is indicated by lines connecting unfilled squares) through a fold (or saddle-node like) limit cycle bifurcation (the bifurcation point is represented as SN in Fig. 2) where the stable limit cycle collides with the unstable ones and they get annihilated via fold limit cycle. The occurrence of fold limit cycle bifurcation is further illustrated via the subset of Fig. 2. After the destabilization of the OPS state, the same oscillatory (unstable) branch, undergoes folding and stabilizes as a three periodic OPS (OPS(3P)) state (stable nature of it is denoted by filled brown colored triangles) in the regions $S_{2}$ and $S_{3}$. In this OPS(3P) state, the oscillators are aligned with a $\frac{2\pi}{3}$ phase difference (and its temporal behavior looks like the one given in the subset of Fig. 3).

In the above case, the OPS(3P) state is not stabilized soon after the destabilization of OPS state and so in the intermediate region (namely, $S_{1}$), IPS state is the only stable state. It is to be noted that these initial conditions lead to stable OPS states in the region $S_{0}$ which will also lead to stable IPS state asymptotically in the parametric range of $S_{1}$, while the transient dynamics in this region non-trivially shows longer time transient chaotic behavior or super persistent transient chaos as shown in Fig. 3(a). The unstable orbits present in the region $S_{1}$ constitute a chaotic saddle which results in the long transient chaotic evolution. The fold or saddle-node type limit cycle bifurcations that occur in the $S_{0}$ and $S_{3}$ regions may be responsible for the emergence of such chaotic saddle. One can also observe the short intervals of OPS(3P) states among the super persistent chaotic dynamics as illustrated in the inset of Fig. 3(a) and it mimics intermittent like behaviors in transient chaotic evolutions. Such an intermittent transient chaotic nature is more pronounced if one examines the dynamics near the boundary of the region $S_{2}$. The chaotic nature of the transient dynamics is confirmed using $0-1$ test and the details of which are illustrated in Appendix A.1. The transient chaotic nature can also be corroborated from the transient time taken by the system corresponding to different initial conditions of depicted in Fig. 3(b). In this figure, we have fixed $y_{1}(0)=0.4$, $x_{2}(0)=1.8$, $y_{2}(0)=0.2$, $x_{3}(0)=2.0$, $y_{3}(0)=0.1$ and calculated the transient times $t_{tr}$ for different values of $x_{1}(0)$, where the non-smooth variation of $t_{tr}$ with respect to various initial conditions elucidate the presence of chaotic saddles. Further, it is also evident from the sensitive dependence of $t_{tr}$ with respect to the initial conditions as illustrated in the subset of Fig. 3(b). This transient chaotic state is not only observed in the region between the stable arena of OPS and OPS(3P) states (region $S_{1}$) in Fig. 1(b) but also near the boundaries of stable regions of the OPS state (near the boundary of $S_{0}$) and OPS(3P) state (near the boundary of $S_{2}$). But in the latter regions, the transient chaotic dynamics exists only for a finite time and super persisting chaotic transients are observed only in the region $S_{1}$.

After the parametric regions with transient chaos and stable OPS(3P) states, Fig. 2 displays the existence of stable OD state for larger values of $\epsilon_{1}$ (regions $S_{3}$ and $S_{4}$). In Fig. 2, we observe four branches of lines ($B_{1}$, $B_{2}$, $C_{1}$ and $C_{2}$) which correspond to the two steady state configurations of $B$ and $C$ in Eqs. (3)-(4). Depending on the initial condition, either of the steady state configurations $B$ or $C$ will be stabilized. If $B$ ($C$) is the asymptotic state, two oscillators will settle at the value given by $B_{2}$ ($C_{2}$) while the other takes up the value $B_{1}$ ($C_{1}$). Figure 2 shows that the stabilization of these inhomogeneous steady states occur via the supercritical Hopf bifurcation in contrast to the one observed for lower values of $\epsilon_{2}$ (and also in the case of two coupled system) where OD occurs via a sub-critical Hopf bifurcation. The parametric ranges over which the OD state is stabilized through sub-critical and super-critical Hopf bifurcations are shown in the Appendix A.2. (see Fig. 8(c)). The inhomogeneous limit cycles resulting from this supercritical Hopf bifurcation are found in a very narrow range of the parametric space as shown in Fig. 2 (also in the Figs. 1(b) and 1(c)).

2.5 Chaotic itinerancy and fractal basins

In addition to the super persistent transient chaos, other interesting dynamics including inhomogeneous oscillations and stable and finite time transient chaos can be observed in the region between OPS and OD states. To illustrate these, we have fixed the strength of the repulsive coupling as $\epsilon_{2}=1.7$ and the corresponding bifurcation diagram along the line $L_{3}$ in Fig. 1(b) is depicted in Figs. 4(a) and 4(b) ( which is a closer view of Fig. 4(a)). As discussed for $\epsilon_{2}=1.0$, Fig. 4(a) depicts the fact that the OPS state, which is stable for lower values of $\epsilon_{1}$, becomes unstable through a fold (saddle-node like) limit cycle bifurcation as a function of $\epsilon_{1}$. The OPS(3P) state is stabilized even before the destabilization of OPS state and so there is no void between the out-of phase synchronized states (that is, OPS and OPS(3P)) as a function of $\epsilon_{1}$ in contrast to the one observed in Fig. 2. The phase space trajectory of OPS(3P) state is shown in Fig. 4(c). Increasing $\epsilon_{1}$, two inhomogeneous oscillatory branches (shown by green and pink curves with filled diamonds) and a homogeneous oscillatory branch (shown by orange colored lines with unfilled diamonds) are found to be stabilized in the region $I_{1}$. The dynamics of the system in this region of homogeneous and inhomogeneous (H-IN) oscillations is shown in Fig. 4(d) which elucidates that the three oscillators are distributed in different orbits. Even though the three oscillators are in different periodic orbits, they oscillate with the same frequency and are phase synchronized. Taking into account the commutation symmetry of the system, six different configurations are possible in this particular state.

After the destabilization of H-IN state, there are no stable states in the region $I_{2}$ as evident from Fig. 4(a). Before getting into the details of dynamics in the region $I_{2}$, we first discuss the dynamics observed in the region $I_{3}$ where the inhomogeneous (IN) oscillations emerge from the supercritical Hopf bifurcation. For a better visualization of this bifurcation and the oscillatory branches ($BO_{1}$, $BO_{2}$, $CO_{1}$ and $CO_{2}$ ) that emerge from the steady state branches, a portion of Fig. 4(a) is enlarged in Fig. 4(b). A closer look at the dynamics in this range of $\epsilon_{1}$ reveals that the inhomogeneous periodic orbits follow the same scenario followed by the inhomogeneous steady states from which they are bifurcated. In other words, similar to the OD states, the inhomogeneous periodic orbits can exist in two possible configurations, namely, $BO$ or $CO$. If $BO$ ($CO$) is asymptotically stable then any two of the oscillators converge to the orbit $BO_{2}$ ($CO_{2}$) (see Fig. 4(b)) while the third oscillator will be designated to the orbit $BO_{1}$ ($CO_{1}$) as shown in Fig. 4(b). Six different configurations are possible in this IN state giving rise to multistability in the region $I_{3}$. Figure 4(e) illustrates the dynamics of the three coupled vdP oscillators in the region $I_{3}$ for a particular initial condition. It shows that two of the oscillators evolve in the same periodic orbit (namely, $BO_{2}$) with a $\pi$ phase (anti-phase) difference while the third oscillator occupies a separate orbit (namely, $BO_{1}$). The frequencies of all the three oscillators are found to be the same in this IN state despite their amplitudes and shapes of the periodic orbits are not symmetric. It is interesting to note that in order to accommodate two oscillators, the orbits $BO_{2}$ and $CO_{2}$ are having larger amplitudes while that of $BO_{1}$ and $CO_{1}$ are smaller. This strongly corroborates an interesting form of spontaneous symmetry breaking where the identical set of oscillators having commutation symmetry exhibit heterogeneous dynamics. It is to be noted that this symmetry breaking is also evident in homogeneous and inhomogeneous oscillating states (see Fig. 4(d)).

The frustration is also prominent in the regions of OPS, OPS(3P), H-IN, IN and OD states where the system cannot accommodate all the three oscillators in anti-phase with each other. In the IN state, two of the oscillators evolve in anti-phase with each other while the third oscillator cannot do so and so the frustration effects are more prominent here. The frustrated system normally has the tendency to relax itself into a non-frustrated state and as a consequence the inhomogeneous oscillations lose their stability in the region $I_{2}$ and show switching between different possible states in an erratic way[4]. For instance, we have depicted the above mentioned metastable dynamics in Fig. 5(a) and the peculiar dynamical states observed in different time windows are displayed in Figs. 5(b)-5(e). The oscillations around the unstable OPS(3P) state can be seen in the time interval $r_{1}$ (see Fig. 5(b)) and in the time window $r_{2}$ two of the oscillators are in anti-phase state while the third oscillator is moving independently as illustrated in Fig. 5(c). In the intervals $r_{3}$ and $r_{4}$, the first and third oscillators, and the second and third oscillators, respectively, are in anti-phase with the other two as shown in Figs. 5(d) and 5(e), respectively. Such a type of erratic wandering among different unstable periodic orbits is known as chaotic itinerancy or frustrated chaos and this type of metastable dynamics is of great interest in neurodynamics where the spontaneous switching between different states enables sequential memory recalling. The region with chaotic itinerant behavior is denoted as $I_{2}$ in Figs. 1(b) and 1(c) in a certain range of $\epsilon_{2}$ which is suppressed largely for large $\epsilon_{2}$. This is because the stabilization of the inhomogeneous steady states occur via a sub-critical Hopf bifurcation (shown in Fig. 8(c) of Appendix A.2.) for large $\epsilon_{2}$ by destabilizing the stable inhomogeneous limit cycles in the corresponding ranges of $\epsilon_{2}$. Due to the sub-critical Hopf bifurcation, the emerging unstable periodic orbits coexist with the stable inhomogeneous steady states and so there will not be any chaotic states for large $\epsilon_{2}$.

Even though the above mentioned stable chaotic behavior in coupled limit cycle oscillators exist for a shorter range of parameters, its impact in the form of transient chaos is wide over the parametric ranges. One can find the existence of transient chaos widely around the stable chaotic region (region - $I_{2}$), that is in the regions $I_{1}$, $I_{3}$ and near the boundaries of $S_{5}$ and OD regions. Even in the present case, we observe the sudden destabilization of chaotic orbit and it does not follow any of the routes other than the existence of transient chaos, which is found to be short lived in time. This type of transient chaos around the chaotic regions may arise due to crisis [52]. In Appendix A.3., we have deduced the mean exit time scales of the transient chaos.

The observed transient chaotic regions ($I_{1}$, $I_{3}$, $S_{5}$ and OD regions) are multistable. The erratic transient feature along with the multistable nature results in interesting intertwining between the basins of attraction of the different states as can be seen in the form of fractal basins or in the form of extreme type of fractal basin, namely, the riddled basin [37, 52, 53]. For example, the basin of attraction for $\epsilon_{1}=1.16$ and $\epsilon_{2}=1.7$ is depicted in Fig. 6(a) and its subsequent enlarged figures show a closer view of the basin. Figure 6(a) and its subsets corroborate the riddled form of the basin where the strong intertwining among the attractors is evident. We have also displayed the basin of attraction of the OD states for $\epsilon_{1}=1.125$ and $\epsilon_{2}=1.7$ in Fig. 6(b) which shows the fractal nature of the basin. The dynamics of the system in these fractal and riddled basins are also illustrated Appendix A.4.. Further in order to corroborate and quantify the fractal nature of the basin, the uncertainty exponent ($\tilde{\alpha}$) is calculated for different vales of $\epsilon_{1}$ and is shown in Fig. 6(c) for a fixed $\epsilon_{2}$ (the procedure for calculating $\tilde{\alpha}$ is discussed in Appendix A.5.). Note that $\tilde{\alpha}$ is related to the probability that two nearby initial conditions lead to different attractors and $\tilde{\alpha}<1$ denotes the fractal nature of the basin while $\tilde{\alpha}\approx 0$ denotes the riddled basins [52, 53]. The basins are riddled very near to the chaotic itinerancy region as confirmed by the values of the uncertainty exponent $\tilde{\alpha}$ in Fig. 6(c), revealing the high sensitivity of final asymptotic states with respect to the initial conditions. Further, the sensitivity to the initial conditions reduces as $\epsilon_{1}$ deviates from the chaotic region ($I_{2}$) as is evident from the value of $\tilde{\alpha}$ with their corresponding basins as fractal basins (see Fig. 6(b)). Thus it is strongly corroborated that even coupled limit cycle oscillators can show transient chaotic behaviors that results in fractal or riddled basins. The emergence of transient chaos around chaotic itinerant region has not yet been reported in the literature to the best of our knowledge. Such a type of transient chaotic states around the chaotic itinerancy regions allow us to inherit both the converging dynamics and chaotic wandering dynamics that are essential for achieving cognitive tasks and brain inspired computing paradigms as pointed out in the introduction. This is because that in certain situations it is required that the system to remain in a particular memory state represented by a limit cycle or steady state (until the arrival of external stimuli) and in certain situations, the system is required to transit between different associated memory states to do sequential tasks or sequential memory recalling. Hence, it is highly desirable and important to inherit both the type of dynamics in a simple system of coupled limit cycle oscillators.

2.6 Dynamics in large networks

Systems with increasingly large number of dynamical states are more desirable than that with a few number of states in view of the larger multistability required to represent different operations or memory states. Thus the search for the onset of transient chaos and fractal basins with more number of dynamical states is extended to a network of coupled vdP oscillators. Stable regions of different dynamical states are depicted in Fig. 7(a) for globally coupled vdP oscillators with $N=100$. The OPS, IPS and OD states exist even in the case of $N=100$ oscillators. Segregation of oscillators into three clusters is observed in the OPS state (represented as OPS(3C)) and the three periodic OPS state in three clusters is represented as 3C(3P) in Fig. 7(a). Stable chaos is observed in the parametric range where super persistent transient chaos exists for $N=3$. We found that an increase in the size of the network leads to its stabilization which is illustrated in Appendix A.6. It is also to be noted that the stabilization of chaotic attractor in this region does not affect the stability of synchronized periodic solution. Thus interesting multistability between the completely ordered and completely disordered or chaotic states is observed around equal attractive and repulsive couplings $\epsilon_{1}\approx\epsilon_{2}\approx 1$. Existence of chaos due to the balanced excitatory and inhibitory couplings was one of the interesting research topics in the past [54].

Now considering the chaotic itinerant region and the transient chaotic dynamics, one can observe that the chaotic itinerancy is absent for $N>3$ but the transient chaotic nature prevails in all the cases which intertwines the homogeneous and inhomogeneous periodic orbits. In the region denoted by H-IN in Fig. 7(a), we observe multi-clusters with a few of the oscillators distributed in the homogeneous oscillatory branch (with the center at zero) and the others at the inhomogeneous oscillatory branches. For instance, Fig. 7(b) shows a five cluster state for $\epsilon_{1}=1.3$ and $\epsilon_{2}=3.0$ where the three groups of oscillators evolve in the homogeneous periodic orbit with $\frac{2\pi}{3}$ phase difference while the other two groups are distributed in the inhomogeneous orbits. Similarly, for certain parametric ranges, seven clusters can be observed where two more groups occupy the inhomogeneous orbits and so evolution of the two groups of oscillators with $\pi$ phase difference can be observed in each of the inhomogeneous orbits. In these regions, the transient chaotic behavior intertwines the basin. When we consider the H-IN region, the basins are found to be riddled and show more sensitivity towards initial conditions. However the sensitivity towards the initial condition is reduced when we consider the H-IN state that emerge for higher values of $\epsilon_{2}$. For example, in Fig. 7(c), we have shown the distribution of oscillators among the homogeneous and inhomogeneous limit cycles for four different initial conditions. The initial conditions are chosen as $x_{i}(0)=\zeta_{i}+\eta\tilde{\zeta}_{i}$, where $\zeta_{i}$’s are distributed randomly between $-1.0$ to $1.0$ and $\tilde{\zeta}_{i}$’s between [math] and $1$. $\eta$ decides the deviation of the initial conditions and we have chosen $\eta=0,10^{-4},10^{-3},10^{-2}$. The figure shows that when $\eta$ is of the order of $10^{-4}$ (or even lesser), there are no variations in the separation of oscillators into different orbits or no variation in the asymptotic state. But when the perturbation is increased to the order of $\eta=10^{-3}$ and $10^{-2}$, there are deviations in the asymptotic state mimicking the fractal nature of the basins. Such fractal nature may not present in the H-IN state that exist for lower values of $\epsilon_{2}$ where the basins are found to be riddled. But considering the OD states that are present near the boundary of H-IN state, their basins are found to be of fractal nature. Thus it is interesting to note that we have not only shown the possibility to have fractal basin structures in oscillatory states but also in the OD state.

2.6.1 Stability analysis

The stable boundary for synchronized region obtained through master stability analysis [55, 56, 57]. In the synchronized manifold, $x_{i}=x$, and $y_{i}=y$, $\forall i$, then the variational equation of system (1) can be written as,

[TABLE]

where $x(t)$ and $y(t)$ are the solutions for the uncoupled system (1). ${\eta_{ij}}$ $(i=1,2)$ is the perturbation from the synchronized manifold, where $\eta_{ij}={\boldsymbol{\beta_{ij}K_{j}}}$. Here $\boldsymbol{{\beta_{ij}}}=(\beta_{i1},\beta_{i2},...,\beta_{iN})$ and $\boldsymbol{{(\beta_{1j},\beta_{2j})}}$ are the deviation of ($x_{j}$, $y_{j}$) from the synchronized solution $(x,y)$. $\boldsymbol{K_{j}}$ is the eigenvector of the coupling matrix corresponding to the eigenvalue $\lambda_{j}$. The coupling matrix $\boldsymbol{K_{j}}$ is having the elements $a_{ij}={-(N-1)}/{N}$ for $i=j$ and $a_{ij}=a_{ji}={1}/{N}$ for $i\neq j$, $(j=1,2,...N)$. To find the stable region, we obtain the eigenvalues from the coupling matrix ${\boldsymbol{K}_{j}}$, which is $\lambda_{0}=0$ and $\lambda_{k}=-1$, $k=1,2,...N-1$. The eigenvalue $\lambda_{0}=0$ implies the perturbation parallel to the synchronization manifold and the remaining $N-1$ eigenvalues designates the perturbations transverse to the synchronization manifold. In the stable synchronization manifold the transverse eigenmodes should be damped out. By substituting the eigenvalues in Eq. 5 and by finding the largest Lyapunov exponents the stable boundary of the synchronized region is obtained and is shown in Fig. 7.

Further, using linear stability analysis, we can find the stable boundary of the oscillation death state. In this case, the oscillators in the network tend to an anti-symmetric inhomogeneous steady states which are of the form

[TABLE]

The corresponding eigen values are,

[TABLE]

where, $d_{1}=(1-x_{1}^{2})\alpha$, $d_{2}=\alpha^{2}-4(1+2x_{1}y_{1})+x_{1}^{2}\alpha^{2}(x_{1}^{2}-2)$ and $d_{3}=d_{1}+2(p-1)(\epsilon_{1}-\epsilon_{2})^{2}+4(2(p-1)(x_{1}^{2}-1)\alpha\epsilon_{1}-1-2x_{1}y_{1}\alpha+4(p-1)^{2}\epsilon_{1}\epsilon_{2})$. The stable region of OD state emerges when $\epsilon_{1}\geq 1/(4(p-1)^{2}\epsilon_{2})$ when $\epsilon_{2}<1$ and $\epsilon_{1}\geq 1/(\sqrt{4(1-2p-p^{2})})$ when $\epsilon_{2}>1$. The stable boundary for OD region is shown in Fig. 7 for $p=0.5$. The analytical regions of synchronized state and oscillation death states match well with the numerical boundaries.

3 Conclusion

The present report unraveled the possibility of onset of transient chaos and chaotic wandering even in a simple network of autonomous limit cycle oscillators. We have elucidated that the transient chaotic nature that stems from frustration acts as a source for this type of dynamical behavior. The transient chaotic nature in combination with the frustration induced multistability gives rise to fractal and riddled basin structures. It is interesting to note that the complex basin structures have not only been observed in the case of oscillatory states but also in the case of oscillation death states. Further, we have shown the existence of finite time transient chaotic dynamics in parametric regions around the chaotic itinerancy region including the existence of fractal and riddled basins for the first time in the literature using coupled limit cycle oscillators. The observed transient chaotic and chaotic itinerant region may find applications in the field of neuro-computation. However, the chaotic itinerancy behavior can be seen only in the three coupled system and it can not be observed in the case of larger networks with global coupling. But the transient chaotic dynamics and riddled or fractal basin structures are observed even in larger networks. In addition, we have also shown the existence of super persistent transient chaos in the region $\epsilon_{1}\approx\epsilon_{2}\approx 1$ in smaller networks and it turns into stable chaos with an increase in the size of the network. We firmly believe that large number of states with extreme multistability can be observed even in such a simple system of limit cycle oscillators with suitable coupling which we plan to extend in future in view of its applications pointed out in the introduction.

Acknowledgments

KS sincerely thanks the CSIR for fellowship under SRF Scheme (09/1095(0037)/18-EMR-I). SK thanks the Department of Science and Technology (DST), Government of India, for providing a INSPIRE Fellowship. The work of VKC forms part of a research project sponsored by INSA Young Scientist Project under Grant No. SP/YSP/96/2014 and SERB-DST Fast Track scheme for young scientists under Grant No. YSS/2014/000175. DVS is supported by the CSIR EMR Grant No. 03(1400)/17/EMR-II. The work of ML is supported by DST-SERB Distinguished Fellowship.

Appendix A

A.1 $0-1$ test to check transient chaotic nature

Regular and chaotic nature of a deterministic system can be tested using the $0-1$ test [58]. We consider a time series of length $N$ corresponding to the transient region to test chaotic nature of the transient dynamics using $0-1$ test. The first step along this direction is to compute the translation variables, $p(n)$ and $q(n)$, given as

[TABLE]

In the above, $x(j)$ represents the time series taken during the transient time. The parameter $c$ in Eq. (8) can have any value between [math] to $2\pi$. We have taken $c=0.7$ and $N=50,000$. We have plotted the dynamics of the translational variables in Fig. 8(a), where we have ensured that the transient dynamics prevails up to $N=50,000$. Brownian like motion observed in Fig. 8(a) corroborates that the transient dynamics is chaotic. Further, the chaotic nature of the transient dynamics can also be confirmed by studying the diffusive or non-diffusive nature of the translational variables where the latter can be identified using mean square displacement $M(n)$,

[TABLE]

The quantity $M(n)$ becomes bounded when the dynamics is regular and it becomes unbounded when the dynamics is irregular. The behavior of $M(n)$ corresponding to the transient dynamics is given in Fig. 8(b) and the unbounded nature observed in the figure clearly confirms that the transient dynamics is chaotic.

A.2 Hopf bifurcation curves

Figure 8(c) illustrates the Hopf bifurcation points in the $\epsilon_{1}$-$\epsilon_{2}$ space. The figure shows the three different Hopf bifurcation points, HB1, HB2, HB3. The inhomogeneous steady states are stabilized after one of these Hopf bifurcation points. For values of $\epsilon_{2}$ lying the under line $C_{1}$, the stabilization of inhomogeneous steady states is found to occur via sub-critical Hopf bifurcation HB1 (as observed in the case of two coupled systems). Increasing $\epsilon_{2}$, the OD transition occurs via the super-critical Hopf bifurcation (HB3). The latter bifurcation is responsible for the inhomogeneous oscillations that are observed in the IN region and the frustrated chaos in the region $I_{2}$. As shown in Fig. 8(c), above the line $C_{2}$, the stabilization of the OD state occurs via the sub-critical Hopf bifurcation so that the IN state and frustration induced chaos are suppressed with the increase of $\epsilon_{2}$ beyond $C_{2}$.

A.3 Transient chaos induced by crisis: Mean exit time scale

As mentioned in the main text, destabilization of the chaotic orbits in the region $I_{2}$ does not seem to follow any particular bifurcation route other than self manifestation of it as transient chaos. Such a transient chaotic state preceded by stable chaos is often referred(occur) as(via) boundary crisis in which the collision of chaotic attractor and unstable periodic attractor or its stable manifold occurs [59, 60, 61]. To elucidate the above, we have calculated the mean exit time scales in the transient chaotic region. Suppose that $\epsilon_{1}^{*}$ is the point at which the destabilization of the chaotic attractor occurs and the chaotic attractor is found to be stable in the region $\epsilon_{1}<\epsilon_{1}^{*}$ (or equivalently $\epsilon_{1}>\epsilon_{1}^{*}$). In the region $\epsilon_{1}>\epsilon_{1}^{*}$ (or equivalently $\epsilon_{1}<\epsilon_{1}^{*}$), transient chaotic evolution prevails where the system behaves chaotically upto a time $\tau$ and then it reemerges from the chaotic attractor to a stable periodic or point attractor. The mean exit time scale ($\langle\tau\rangle$) can be obtained by averaging $\tau$ for different initial conditions. In a variety of dynamical systems, it has been illustrated that this mean exit time scale follows a power law relationship for the values of $\epsilon_{1}$ near to the critical value $\epsilon_{1}^{*}$ as

[TABLE]

In the above, $\gamma$ is often known as critical exponent and it depends on the stability properties of the basic periodic orbit [62, 63].

To understand the destabilization of chaos due to boundary crisis, the mean exit time is plotted for different values of $\epsilon_{1}$ in Fig. 9(a) for $\epsilon_{2}=1.7$. The figure shows that the chaotic dynamics prevails in the region $IV$, $1.128\leq\epsilon_{1}\leq 1.159$, and so the mean exit time blows up to infinity at the two critical points $\epsilon_{1}^{*}=1.128$ and $\epsilon_{1}^{*}=1.159$. Before ($\epsilon_{1}<1.128$) and after ($\epsilon_{1}>1.159$) this chaotic region, the transient chaotic behavior is found to exist. First, let us consider the chaotic destabilization at the boundary between regions $IV$ and $V$ or at $\epsilon_{1}^{*}=1.159$. At this critical point $\epsilon_{1}^{*}=1.159$, the collision of chaotic attractor with the inhomogeneous periodic orbits corresponding to IN state gives rise to chaotic saddles and so we observe the finite time chaotic evolutions. When the IN state becomes more stable, the widening of the IN state basin occurs and it results in the decrease of mean exit time for the values of $\epsilon_{1}>1.159$. By fitting the data to power law form as given in Eq. (10), we have found the critical exponent as $\gamma=0.367$.

Similarly, when considering the transient chaotic regime present below $\epsilon_{1}^{*}=1.128$, we found that the mean exit time given in Fig. 9(a) shows power law for the values of $\epsilon_{1}$ near to $\epsilon_{1}^{*}=1.128$ and not for the ones lying away from $\epsilon_{1}^{*}$. When $\epsilon_{1}$ deviates from $\epsilon_{1}^{*}=1.128$, the mean exit time shows sudden increase in its value and attains a maximum near the boundary between regions $I$ and $II$. However, in the region $I$, the mean exit time again decreases. Such an absence of monotonic decrease in the mean exit time scale has also been reported in [62]. Here, at the critical value $\epsilon_{1}^{*}=1.128$, the collision between the chaotic attractor with the periodic orbit corresponding to HIN state results in the transient chaotic behavior. As $\epsilon_{1}$ deviates from $\epsilon_{1}^{*}$, the HIN state become more stabilized resulting in widening of their basin so that the decrease in the mean exit time occurs for the values of $\epsilon_{1}$ near to $\epsilon_{1}^{*}$. But for the values of $\epsilon_{1}$ lying near to the boundary between region $III$ and $II$ and also in the region $II$, the HIN state may become weak and paves way for the stable OPS and OPS(3P) states. Due to the above fact, we may observe an increase in the value of mean exit time. However after the stabilization of OPS(3P) state in the region $I$ and due to the increase in the stable nature of the OPS and OPS(3P) states, the mean exit time starts decreasing in region $I$. Due to this fact, the mean exit time follows power law relationship only near $\epsilon_{1}^{*}=1.128$ as shown in Fig. 9(a) corroborating the destabilization of chaotic orbit due to the boundary crisis at both the ends of the chaotic region (region - $IV$).

A.4 Dynamics in riddled and fractal basin regions

The sensitivity of the system towards initial conditions in the riddled and fractal basin regions is illustrated in Figs. 9(b) and 9(c) which are plotted for $\epsilon_{1}=1.16$ and $\epsilon_{1}=1.19$, respectively. By considering two different initial conditions that differ in $y_{1}(0)$ by $10^{-6}$ (and keeping other initial condition values fixed), we have plotted their dynamics in Fig. 9(b). The latter shows that the trajectories corresponding to two initial conditions are close only up to a finite time. After that, the two trajectories rapidly deviate from each other and evolve towards different attractors. But in the fractal basin, the system does not show extreme sensitivity which is obvious from Fig. 9(c). The closer initial conditions lead to the same attractor state whereas well separated initial conditions lead to different attractors.

A.5 Uncertainty exponent

Both fractal and riddled basins correspond to non-smooth basin boundaries that make sensitive dependence of final states towards initial conditions. Compared to fractal basins, the riddled ones show extreme sensitivity towards initial conditions. Considering the case of fractal basins, the probability to predict the asymptotic state is improved when the initial condition is more precisely defined. In the case of riddled basins, the final state is found to be uncertain even if the initial condition is precisely defined. These two types of basins can be classified by calculating the uncertainty exponent [64]. The latter is related to the probability that the two nearby initial conditions lead to same asymptotic state or not. For this purpose, the attractor reached by a set of randomly distributed or equally spaced initial conditions $\vec{X}_{i}(0)(=(x_{i}(0),y_{i}(0))$) has to be determined and then one has to determine the asymptotic state of another set of initial condition which is deviated from the earlier set by a finite value $\epsilon_{i}$ as $\vec{X}_{i}(0)+\vec{\epsilon}_{i}$. The probability ($P(\epsilon)$) that this pair of nearby initial conditions ($\vec{X}_{i}(0)$ and $\vec{X}_{i}(0)+\vec{\epsilon}_{i}$) leading to different asymptotic states is computed for different values of separation (${\epsilon}$). The uncertainty exponent is determined by fitting the scaling law $P(\epsilon)\sim\epsilon^{\tilde{\alpha}}$ where $\tilde{\alpha}$ represents the uncertainty exponent. For riddled basins $\tilde{\alpha}\approx 0$.

The probability of the system to reach the same attractor for a set of initial conditions that differ by length $\epsilon$ is shown Fig. 10(a). The slope of the straight line shown in Fig. 10(a) gives the value of the uncertainty exponent $\tilde{\alpha}$.

A.6 Transient chaos to stable chaos

As seen in the main article, the system exhibits long persisting transient chaos for $N=3$ in the region where the value of the attractive and repulsive couplings are approximately unity. The transient chaos turns out to be stable chaos when the size of the network is increased. To illustrate the above, we have plotted the probability of the system to reach the periodic state ($P_{tr}$) at a particular time in Fig. 10(b). For different values of $N$, $P_{tr}$ is calculated at times $t=10^{4},10^{5}$ and $10^{6}$ by considering different initial conditions. For $N=3$ and $t=10^{4}$, we observe that the system remains in the chaotic state for many of the initial conditions and so $P_{tr}$ takes finite values whereas for $t=10^{6}$, the system reaches the periodic state for almost all the initial conditions so that $P_{tr}\approx 1$. By increasing $N$ (upto $N=9$), the lifetime of transient chaos is also increased so that $P_{tr}$ takes up finite values even at $t=10^{6}$. $P_{tr}$ becomes zero for $N>9$ in Fig. 10(b) which indicates the stabilization of chaotic state.