Oscillation death in coupled oscillators

Abstract

We study dynamical behaviors in coupled nonlinear oscillators and find that under certain conditions, a whole coupled oscillator system can cease oscillation and transfer to a globally nonuniform stationary state [i.e., the so-called oscillation death (OD) state], and this phenomenon can be generally observed. This OD state depends on coupling strengths and is clearly different from previously studied amplitude death (AD) state, which refers to the phenomenon where the whole system is trapped into homogeneously steady state of a fixed point, which already exists but is unstable in the absence of coupling. For larger systems, very rich pattern structures of global death states are observed. These Turing-like patterns may share some essential features with the classical Turing pattern.

Graphical abstract

Keywords

coupled oscillators / oscillation death (OD) / amplitude death (AD) / synchronization

1 Introduction

Collective behaviors in coupled nonlinear oscillators [1–5] such as frequency locking [1], complete synchronization (CS) [6–11], anti-phase synchronization (AS) [12, 13], and oscillation death (OD) [14–41] have attracted much attention in the past several decades. These phenomena are ubiquitous and fundamental in physical, chemical, biological, and even social sciences, and recently remarkable achievements have been made in both theoretic developments and experimental studies. For the so-called frequency locking (or phase synchronization), under certain conditions the oscillators can be entrained to a common frequency when they are coupled. In contrast, complete synchronization means that all trajectories of coupled oscillators go to the same trajectory as the system evolves, and anti-phase synchronization means they have the same amplitude but differ in sign. For the other particularly interesting amplitude effect, oscillation death refers to the phenomenon where the coupled oscillator system cease oscillation, become dead in time, and exhibit a stationary state. Generally speaking, two types of OD exist for the different relations between an isolated oscillator without coupling and an oscillator under coupling. Namely, the end stationary state can be the same, which already exists and is an unstable fixed point in the absence of coupling, and the coupling simply changes its stability and makes it observable. As fixed points with zero amplitude are usually considered, this type of oscillation death is generally termed amplitude death (AD) for its null amplitude. Conversely, the eventual stationary state can be entirely newly created by the coupling. For this type of oscillation death, oscillators show a none-zero amplitude, and thus, it is simply called OD, to be opposite to AD. As the two types of death state in coupled oscillators, both AD and OD are important for the changes of system dynamics. Moreover, from the point of view of controlling oscillatory dynamics, their studies are clearly of great significance for potential applications. In theoretic studies, AD has already been investigated by many researchers, but obviously studies of OD are still few [27]. In experimental studies, however, OD was extensively reported. See, e.g., Refs. [22–24, 37, 38]. The distinction between AD and OD is not so great as in modelling studies, and usually we do not know the exact information for unstable states of experimental systems.

Coupled limit cycles provide a simple but powerful mathematical model for simulating collective behaviors of coupled nonlinear oscillators [1, 2]. One of the classical models for single limit cycle, the Landau-Stuart oscillator,

\begin{array}{}\text{(1)}& \dot{z}=(1+\mathrm{i}w-{|z|}^2)z\end{array}

represents a weakly nonlinear system near a supercritical Hopf bifurcation. Here, w denotes the intrinsic frequency and z(t)=r(t){\mathrm{e}}^{\mathrm{i}\theta (t)}=x(t)+\mathrm{i}y(t). The model describes a stable limit cycle of unit amplitude (|z|=1) on which it moves at its angular frequency w, and apart from this, an unstable fixed point at the origin (z=0) exists.

A simple model for a pair of Landau-Stuart oscillators, in which two of them are coupled linearly to each other, is

\begin{array}{}\text{(2)}& \dot{z_i}=(1+\mathrm{i}{w}_i-{|z_i|}^2)z_i+ϵ(z_j-z_i)\end{array}

where i,j=1 or 2, and ε is the coupling intensity. It has been well known that with the change of the coupling strength ε and the frequency mismatch between the oscillators \Delta =|w_1-w_2|, the system can exhibit very rich phenomena [21]: (1) phase drifting under small coupling (), where the oscillators behave independently and the relative phase of the two oscillators drifts through all phases, (2) frequency locking, where the two oscillators transit to a common frequency, and (3) amplitude death, in which the oscillators collapse into the origin (r_1=r_2=0) as time reaches infinity. As a result, the occurrence of AD needs both the coupling strength and the frequency mismatch to be sufficiently large [\Delta >2 and (1+{\Delta }^2/4)/2>ϵ>1], implying that AD is not available in coupled identical Landau-Stuart oscillators.

Considering time delay of information propagation, Reddy et al. [33, 34] studied the following coupled time-delay equations,

\begin{array}{}\text{(3)}& \dot{z_i}=(1+\mathrm{i}{w}_i-{|z_i|}^2)z_i+ϵ(z_j(t-\tau )-z_i(t))\end{array}

where i,j=1 or 2, and \tau ⩾0 is a measure of the time delay. Obviously, Eq. (2) is a special case of Eq. (3) for \tau =0. They found that for appropriate delay times (\tau \ne 0) AD can be generally observed even in coupled identical systems (w_1=w_2). This finding has been immediately verified by several experiments, such as a pair of electronic circuits [35], coupled living oscillators [36], thermo-optical oscillators [37], chemical oscillators [38], laser systems [39], etc.

Thus, time-delayed coupling and frequency mismatch of oscillators seem to have become an alternative necessary condition for the occurrence of AD, namely, to observe AD, we have to construct a coupling by tuning a proper delayed time, or set a pair of oscillators having a sufficiently large parameter mismatch. One may ask: can we release these two restrictions at the same time, or can we find death state in coupled identical oscillators in the absence of parameter mismatch and only through the interaction of simultaneous information of coupled oscillators? It seems that the only way to do so is to construct other couplings, which are different from that in Eq. (2). Several researchers have already tried to answer this interesting question. For instance, Konishi [40] tried a dynamic coupling. Different from the usual static coupling, the proposed dynamic coupling is the coupling with the combination of several variables, which depend on time. As already stated in the paper, however, it is a rough approximation of the time-delayed coupling for low-frequency oscillators and/or short-time delay. Thus, basically it is a time delay coupling. Very recently, Karnatak et al. [41] proposed another novel coupling scheme, dissimilar (or “conjugate”) variable coupling. As z(t)=x(t)+\mathrm{i}y(t), and x and y are two conjugate variables, the Landau-Stuart oscillator can be specified in Cartesian coordinates, \dot{x}=px-wy,\dot{y}=py+wx (p=1-x^2-y^2). They considered the couplings, using the variable x of one oscillator to the \dot{y} term of the other, or the variable difference y_2-x_1 to the \dot{x} term of one oscillator. See Ref. [41] for more details. Although this conjugate variable coupling is natural in a variety of experimental situations, such as coupled semiconductor laser systems investigated in Ref. [42], where subsystems are coupled by feeding the output of one into the other, the underlying mechanism still comes from the time-delayed signal propagation of one single variable. The reason is simple: From Takens’s embedding theorem [43], we know that nonlinear systems can be reconstructed by the so-called delay coordinate vector with appropriate choices of embedding dimension and time delay. (This technique has been extensively used in time series analysis.) Thus, the conjugate variable is essentially equal to the time-delayed variable with a certain delayed time. It is not truly devoid of the effect of time delay. Note that the achieved end state by using such a conjugate variable coupling can be either an OD or AD state. It is only natural to obtain AD (or OD) for these two types of coupling, as there is not much difference between them and the time-delayed coupling considered by Reddy et al. [33, 34].

Consequently, the question for the possibility of death state in coupled identical oscillators is still not answered. In experiments, however, the transitions from oscillation to death in coupled oscillator systems have already been extensively observed, for instance, in the symmetrically coupled two [15] or three [16] identical (or nearly identical) chemical oscillators. On the other hand, many efforts have been devoted recently to the studies of death state in coupled systems with a large number of oscillators [25–32]. Both global [25–27] and local [28–32] couplings have been studied. For a local coupling, a specific case where the natural frequencies of oscillators are distributed in a regular monotonic trend was mostly treated and only partial death state with some oscillators being dead and the other still oscillating has been observed [28–30]. Novel effects were studied even on complex networks [31]. Because a large mismatch of frequencies of neighboring oscillators is necessary for the occurrence of AD, it is natural to believe that to obtain a whole AD pattern in populations of oscillators coupled locally with randomly chosen natural frequencies, an extremely large coupling is needed [32]. For an intermediate coupling strength under usual conditions, AD is impossible. This may even get worse for larger systems.

To fill the gap between theory and experiments and search for death states in coupled oscillators, in the present work we intend to study these problems in coupled identical (or nearly identical) nonlinear oscillators with another type of coupling-regular coupling. Compared with the dissimilar (conjugate) coupling, which is believed as an irregular coupling, the regular coupling indicates that the units are coupled by feeding the output of the same variables isochronously. The usual diffusive coupling will mainly be investigated. Different from the mostly studied scalar coupling [i.e., the same coupling strength for x and y with Cartesian coordinates in Eq. (2)] with only one parameter changeable, in this work we will study a non-scalar coupling with two changeable coupling parameters. We find that OD exists under broad conditions. For larger-size coupled systems, with the same type of diffusive coupling, we find that OD state is generic and robust. Very rich pattern structures of nonuniform death states are observed and analyzed. Relevant to these stationary patterns are the so-called Turing patterns [5, 44], which show nonuniform stationary patterns, appearing from an initial homogeneous steady state by a novel mechanism found by Turing in 1952, the diffusion-induced instability. Thus, we may call our nonuniform frozen patterns Turing-like patterns. Interestingly, these Turing-like patterns may share some essential features with the classical Turing pattern, although they have very different natures. We hope that all these findings are of significance and help explain the experimental phenomena, possible techniques for controlling of oscillatory systems, and pattern formations in non-equilibrium state systems.

The paper is organized as follows. In Section 2, we study and analyze the dynamics of a pair of regularly coupled oscillators. Both diffusive and direct couplings are investigated. A phase diagram as a function of the couplings on x and y is obtained for each coupling form. Sections 3 and 4 are devoted to the studies on death states in small-size systems and large-size systems, respectively. The properties of the structures of observed patterns are studied. Finally, discussions and conclusion are presented in Section 5.

2 Oscillation death in a pair of coupled oscillators

2.1 Diffusive coupling

First, let us consider the following two linearly coupled Landau-Stuart oscillators,

\begin{array}{}\text{(4)}& \begin{array}{c}\dot{x_i}=p_i{x}_i-w_i{y}_i+{ϵ}_{11}(x_j-x_i)+{ϵ}_{12}(y_j-y_i)\hfill \\ \dot{y_i}=p_i{y}_i+w_i{x}_i+{ϵ}_{21}(x_j-x_i)+{ϵ}_{22}(y_j-y_i)\hfill \end{array}\end{array}

where p_i=1-x_i^2-y_i^2, and \mathit{ϵ} is the coupling matrix with {ϵ}_{11}, {ϵ}_{12}, {ϵ}_{21}, and {ϵ}_{22} being the matrix elements. This equation can be regarded as a canonical model of linearly coupled oscillators and has been studied in Ref. [20]. To be completely different from the dissimilar coupling considered by Karnatak et al., we delete the two crossing terms and simply keep the diffusive coupling terms on the same variables. Namely, we set {ϵ}_{12}={ϵ}_{21}=0, {ϵ}_1={ϵ}_{11}, and {ϵ}_2={ϵ}_{22}. w_1=w_2=w. Without losing generality, w=2 is chosen. Obviously, if {ϵ}_1={ϵ}_2, the above equation degenerates to Eq. (2), and we cannot find AD as w_1=w_2. Therefore, death state is only possible for {ϵ}_1\ne {ϵ}_2.

As a first step, we numerically computed the largest Lyapunov exponent (LLE) of this system for different {ϵ}_1 and {ϵ}_2. Generally speaking, the periodicity of the system (limit cycle attractor) can be recognized by a zero-value of LLE for a typical trajectory from random initial conditions and steady state behavior (AD or OD) can be characterized by a negative LLE. As illustrated in Fig. 1(a), where we plot black points for zero LLE and leave it blank for negative LLE, the stationary state can really exist within very large blank regions. It is notable that in the whole coupling space LLE is equal to or smaller than zero, and no chaos can be found. In this case, there are two sets of fixed point. One is the origin (0,0,0,0), which is unstable for all {ϵ}_1 and {ϵ}_2; this can be easily analyzed. The other is the anti-phase synchronous oscillation death (ASOD) state (x_{a\ast },y_{a\ast },-x_{a\ast },-y_{a\ast }) given by

\begin{array}{}\text{(5)}& \begin{array}{c}{x}_{a\ast }=\pm w\sqrt{\frac{1-p_a}{{(p_a-2{ϵ}_1)}^2+w^2}}\hfill \\ y_{a\ast }=\pm \sqrt{\frac{(1-p_a){(p_a-2{ϵ}_1)}^2}{{(p_a-2{ϵ}_1)}^2+w^2}}\hfill \end{array}\end{array}

where p_1=p_2=p_a=({ϵ}_1+{ϵ}_2)-\sqrt{{({ϵ}_1-{ϵ}_2)}^2-w^2}. Below, we will see that this nonuniform OD state plays a constructive role in both small systems and large systems. The existence conditions for this ASOD are |{ϵ}_1-{ϵ}_2|⩾w and p_a⩽1, which are indicated in Fig. 1(b) by the two dashed lines. For the stability of the ASOD, the characteristic eigenvalue equation is

\begin{array}{}\text{(6)}& {\lambda }^2+[-(a+d)\mp ({ϵ}_1+{ϵ}_2)]\lambda +ad-bc+{ϵ}_1{ϵ}_2\pm (a{ϵ}_2+d{ϵ}_1)=0\end{array}

where a=p_a-2x_{a\ast }^2-{ϵ}_1, b=-w-2x_{a\ast }{y}_{a\ast }, c=w-2x_{a\ast }{y}_{a\ast }, and d=p_a-2y_{a\ast }^2-{ϵ}_2. We get: \mathrm{Re}({\lambda }_{\max })=2p_a-1=2({ϵ}_1+{ϵ}_2)-2\sqrt{{({ϵ}_1-{ϵ}_2)}^2-w^2}-1 within the stable region in Fig. 1(b) marked by 1, \mathrm{Re}({\lambda }_{\max })=2p_a-1+\sqrt{{(1-p_a)}^2-w^2} marked by 2, and \mathrm{Re}({\lambda }_{\max })=-2\sqrt{{({ϵ}_1-{ϵ}_2)}^2-w^2} marked by 3.

It is interesting to study the stability of possible dynamical behaviors (including the above ASOD) from the analyses of complete synchronous manifold and anti-phase synchronous manifold, whose theories have been developed by several researchers [6–11]. The equations of the complete synchronous manifold and anti-phase synchronous manifold are {\dot{x}}_s=p{x}_s-w{y}_s, {\dot{y}}_s=p{y}_s+w{x}_s, and {\dot{x}}_a=p{x}_a-w{y}_a-2{ϵ}_1{x}_a, {\dot{y}}_a=p{y}_a+w{x}_s-2{ϵ}_2{y}_a, respectively. Let {\xi }_1=x_1-x_2,{\xi }_2=y_1-y_2, and {\eta }_1=x_1+x_2,{\eta }_2=y_1+y_2, then the equations determining the evolution of the perturbation that is transverse to the complete synchronous manifold and the anti-phase synchronous manifold are \dot{{\xi }_1}=(p-2x_s^2){\xi }_1-(w+2x_s{y}_s){\xi }_2-2{ϵ}_1{\xi }_1, \dot{{\xi }_2}=(w-2x_s{y}_s){\xi }_1+(p-2y_s^2){\xi }_2-2{ϵ}_2{\xi }_2, and \dot{{\eta }_1}=(p-2x_a^2){\eta }_1-(w+2x_a{y}_a){\eta }_2, \dot{{\eta }_2}=(w-2x_a{y}_a){\eta }_1+(p-2y_a^2){\eta }_2, respectively. We numerically computed the largest transverse Lyapunov exponents (LTLE) of these two manifolds. Only if LTLE is negative, can we observe stable synchronous (inphase or anti-phase) state. The results are illustrated in Fig. 1(c), where the region, enclosed by the solid (dashed) line and marked by the letter S (A), represents the stable complete (anti-phase) synchronous state. There are overlaps in-between, which clearly indicates that in these overlap regions bistable behaviors with both synchronous and anti-phase synchronous behaviors for different initial conditions are possible.

As a result, the parameter space can be divided into several sub-regions for the corresponding different behaviors [Fig. 1(d)]: the ASOD in region I, complete synchronous limit cycle (CSLC) in region II, and anti-phase synchronous limit cycle (ASLC) in region III. Apart from these mono-stable regions, several bistable regions exist in the overlap regions: the coexistence region of ASOD and CSLC (region IV) and the coexistence region of CSLC and ASLC (region V), which are clearly indicated in Fig. 1(d). Figure 2 gives the corresponding attractors for the above five different dynamical behaviors. In particular, we use solid circles to show the attractor of the death state and open circles to show the snapshot of limit cycles. These figures well verify the results in Fig. 1. Further simulations found that the basins of the above coexisting states are not riddled and the division curves are smooth (the numerical results are not shown here).

2.2 Direct coupling

To show the generality of death states in regularly coupled systems, we consider another coupling form, the direct coupling,

\begin{array}{}\text{(7)}& \begin{array}{c}\dot{x_i}=p_i{x}_i-w_i{y}_i+{ϵ}_1{x}_j\hfill \\ \dot{y_i}=p_i{y}_i+w_i{x}_i+{ϵ}_2{y}_j\hfill \end{array}\end{array}

where i,j=1,2 (i\ne j), w_1=w_2=w, and {ϵ}_1 and {ϵ}_2 are coupling strengths, which can be different. Similar to the diffusive coupling discussed above, direct coupling is also very common in nature and engineering [21, 45], such as in electronics. Direct coupling is the connection of the output of one circuit directly to the input of the next, and in theoretical biology, two directly coupled neural oscillators as a model of circadian pacemaker have been well studied [45].

Similar to the calculations in the diffusive coupling, we computed the value of the LLE. The result is shown in Fig. 3(a), which indicates that stationary state can also exist in large blank regions. Now there are three sets of fixed point: the origin (0,0,0,0), which is again unstable for all couplings, the complete synchronous oscillation death (CSOD) state (x_{s\ast },y_{s\ast },x_{s\ast },y_{s\ast }) given by

\begin{array}{}\text{(8)}& \begin{array}{c}{x}_{s\ast }=\pm w\sqrt{\frac{1-p_s}{{(p_s+{ϵ}_1)}^2+w^2}}\hfill \\ y_{s\ast }=\pm \sqrt{\frac{(1-p_s){(p_s+{ϵ}_1)}^2}{{(p_s+{ϵ}_1)}^2+w^2}}\hfill \end{array}\end{array}

where p_1=p_2=p_s=\frac{-({ϵ}_1+{ϵ}_2)-\sqrt{{({ϵ}_1-{ϵ}_2)}^2-4w^2}}{2}, and the anti-phase synchronous oscillation death state (x_{a\ast },y_{a\ast },-x_{a\ast },-y_{a\ast }) given by

\begin{array}{}\text{(9)}& \begin{array}{c}{x}_{a\ast }=\pm w\sqrt{\frac{1-p_a}{{(p_a-{ϵ}_1)}^2+w^2}}\hfill \\ y_{a\ast }=\pm \sqrt{\frac{(1-p_a){(p_a-{ϵ}_1)}^2}{{(p_a-{ϵ}_1)}^2+w^2}}\hfill \end{array}\end{array}

where p_1=p_2=p_a=\frac{({ϵ}_1+{ϵ}_2)-\sqrt{{({ϵ}_1-{ϵ}_2)}^2-4w^2}}{2}. Clearly, the CSOD exists only if |{ϵ}_1-{ϵ}_2|⩾2w and p_s⩽1, and the ASOD exists only if |{ϵ}_1-{ϵ}_2|⩾2w and p_a⩽1. We show their existence regions, which are surrounded by the dashed lines in Fig. 3(b) and (c) for the CSOD and ASOD, respectively.

For the linear stability analysis, the characteristic eigenvalue equation of the CSOD is

\begin{array}{}\text{(10)}& {\lambda }^2+[-(a+d)\mp ({ϵ}_1+{ϵ}_2)]\lambda +ad-bc+{ϵ}_1{ϵ}_2\pm (a{ϵ}_2+d{ϵ}_1)=0\end{array}

where a=p_s-2x_{s\ast }^2, b=-w-2x_{s\ast }{y}_{s\ast }, c=w-2x_{s\ast }{y}_{s\ast }, and d=p_s-2y_{s\ast }^2. Denote the four eigenvalues as {\lambda }_1, {\lambda }_2, {\lambda }_3, and {\lambda }_4, and let \mathrm{Re}({\lambda }_{\max })=\max \{\mathrm{Re}({\lambda }_1),\mathrm{Re}({\lambda }_2), \mathrm{Re}({\lambda }_3),\mathrm{Re}({\lambda }_4)\}. CSOD is stable if and only if . Within the stable (gray) regions in Fig. 3(b), the analytical expressions for \mathrm{Re}({\lambda }_{\max }) are different: \mathrm{Re}({\lambda }_{\max })=-\sqrt{{({ϵ}_1-{ϵ}_2)}^2-4w^2} in the region marked by the number 1, and \mathrm{Re}({\lambda }_{\max })=-({ϵ}_1+{ϵ}_2) in the region marked by 2. Out of the gray region, the CSOD is unstable. A similar analysis is conducted for the stability of the ASOD: \mathrm{Re}({\lambda }_{\max })={ϵ}_1+{ϵ}_2 in the region marked by 3, and \mathrm{Re}({\lambda }_{\max })=-\sqrt{{({ϵ}_1-{ϵ}_2)}^2-4w^2} in the region marked by 4, as shown in Fig. 3(c).

We also numerically calculated LTLE to determine the stability of the synchronization manifold. The equations of the complete synchronous manifold and the anti-phase synchronous manifold are {\dot{x}}_s=p{x}_s-w{y}_s+{ϵ}_1{x}_s, {\dot{y}}_s=p{y}_s+w{x}_s+{ϵ}_2{y}_s, and {\dot{x}}_a=p{x}_a-w{y}_a-{ϵ}_1{x}_a, {\dot{y}}_a=p{y}_a+w{x}_a-{ϵ}_2{y}_a, respectively. The equations for the transverse stability of the complete synchronous manifold and anti-phase synchronous manifold are {\dot{\xi }}_1=(p-2x_s^2){\xi }_1-(w+2x_s{y}_s){\xi }_2-{\xi }_1, {\dot{\xi }}_2=(w-2x_s{y}_s){\xi }_1+(p-2y_s^2){\xi }_2-{\xi }_2, and {\dot{\eta }}_1=(p-2x_a^2){\eta }_1-(w+2x_a{y}_a){\eta }_2+{\eta }_1, {\dot{\eta }}_2=(w-2x_a{y}_a){\eta }_1+(p-2y_a^2){\eta }_2+{\eta }_2, respectively. The results are shown in Fig. 3(d). The complete synchronous manifold is stable within the up-right region ({ϵ}_1+{ϵ}_2>0) and, conversely, the anti-phase synchronous manifold is stable within the down-left region (), which are indicated in the figure by the letters “S” and “A”, respectively. No overlap exists in-between. Note that on the stable (complete or anti-phase) synchronous manifold, the dynamics can be a fixed point or a limit cycle.

So far, the dynamical behaviors have been well classified by the above theoretical analyses, supplemented by the numerical simulations. We obtain four main regions in the phase diagram for CSOD, ASOD, CSLC, and ASLC. It remains uncertain for the dynamics on the anti-diagonal, which is exactly the division line of the synchronization and the anti-phase synchronization regions. Detailed observations find that on the two line segments in Fig. 3(e) ({ϵ}_1+{ϵ}_2=0, , or {ϵ}_1>w) denoted by “III”, a new OD state appears, which can be termed meta-stable oscillation death (MSOD); for different initial conditions the two oscillators can be trapped into different fixed points (locating on one limit cycle), which are neither completely synchronous nor anti-phase synchronous. Conversely, the dynamics on the parameter line ({ϵ}_1+{ϵ}_2=0, w>{ϵ}_1>-w) denoted by “VI” shows a novel periodic behavior, meta-stable limit cycle (MSLC), where the two oscillators rotate around their own periodic trajectories. Finally, we classify the phase diagram in Fig. 3(e) into six different regions I, II, III, IV, V, and VI, for the corresponding dynamical behaviors, CSOD, ASOD, MSOD, CSLC, ASLC, and MSLC, respectively. Their corresponding fixed point behaviors (black points) and the snapshots of periodic trajectories at one time instance (open points on the curves) are exemplified in Fig. 4, which well confirm our theoretical analyses.

2.3 A pair of coupled nonidentical oscillators

Based on the above observations, AD state really does not exist in these coupled identical oscillator systems, but oppositely a novel OD state begins to appear for unequal couplings {ϵ}_1 and {ϵ}_2. This point is clearly distinct from what was expected previously. Compared with the AD state that traps the oscillators into the origin, which is unstable in the absence of coupling and becomes stable due to coupling, the novel OD state is coupling-induced and changes with the coupling (see the explicit formula above). On the contrary, the origin now is unstable, and thus, is unobservable.

Because the system studied in the paper is only a special case for w_1=w_2, it is valuable to investigate how the system changes with the degree of parameter mismatch. Again we calculated the LLE for diffusively coupled nonidentical Landau-Stuart oscillators for unequal w₁ and w₂. The results are shown in Fig. 5(a) and (b) for w_1=6 and w_1=8, respectively, with w_2=2 fixed. Now we can find AD in a window near the diagonal and OD far from the diagonal. With the increase of the frequency mismatch \Delta (\Delta =|w_1-w_2|), the window gets wider and wider. It can break the limit cycle (black) region and connect with the outer OD region if \Delta is larger than a threshold 4.7 [comparing Fig. 5(a) and (b)]. Two dashed lines in Fig. 5(b) are used to divide the AD and OD regions. We may also compare these numerical results with the theoretic result on AD in the earlier work for {ϵ}_1={ϵ}_2=ϵ. As the condition for the AD needs \Delta >2 and (1+{\Delta }^2/4)/2>ϵ>1, the window appears only if {\omega }_1 is larger than 4 (w_2=2), and after that the size of the window on the diagonal direction is (-1+{\Delta }^2/4)/2. All these predictions have indeed been verified by our simulations.

Correspondingly, Fig. 5(c) and (d) are the results for the direct coupling. Now under the nonidentical parameter condition, AD always does not exist. By contrast, the OD persists, but obviously its region decreases with the increase of \Delta. Further study shows that this OD is neither complete synchronous nor anti-phase synchronous, as the system symmetry is immediately broken with the input of the parameter mismatch under the current situation.

3 Stationary patterns in small-size systems

So far, we have studied a model of a pair of linearly coupled oscillators. Next, let us investigate a ring of linearly coupled identical Landau-Stuart oscillators, which can be written as

\begin{array}{}\text{(11)}& \begin{array}{c}{\dot{x}}_i=p_i{x}_i-\omega y_i+{ϵ}_1(x_{i+1}+x_{i-1}-2x_i)\hfill \\ {\dot{y}}_i=p_i{y}_i+\omega x_i+{ϵ}_2(y_{i+1}+y_{i-1}-2y_i)\hfill \end{array}\end{array}

with i=1,2,\mathrm{\cdots },N being the oscillator index. Periodic boundary conditions are considered. Without losing generality, we choose {ϵ}_1=-3and {ϵ}_2=0. \omega =2 is fixed again.

For N=2 [then x_{i+1}=x_{i-1} in Eq. (11)], the only stationary state is the anti-phase oscillation death state, described in Eq. (5). If we denote it as \mathit{a}\equiv (x_{a\ast },y_{a\ast }), then the patterns can be written compactly as (\mathit{a},-\mathit{a}) or, equivalently, (-\mathit{a},\mathit{a}). To better characterize the observed patterns, we may neglect the difference between these interchangeable patterns and redefine a new variable by the death amplitude d_i, d_i^2=x_i^2+y_i^2. Under such a new definition, the inphase and anti-phase states are unified, the exchange of oscillator indices gives no extra pattern configuration, and then the complexity of the death patterns can be conveniently quantified and compared by checking the distribution of d. Under the present parameters ({ϵ}_1=-3 and \omega =2), d=d_0=\sqrt{1-p}=\sqrt{1-2{ϵ}_1+\sqrt{4{ϵ}_1^2-{\omega }^2}}\approx 3.558; this value is consistent with that of d on the plateau in Fig. 6(a) for N=2.

If the system size is larger than two, then no explicit form of death state is available and we have to rely on numerical simulations. Computer simulations show that the patterns of LLE for other N’s do not have much difference with Fig. 1(a) for N=2. The data are not shown here. In Fig. 6(b)–(e), we plot the typical distributions of d for N=3, 4, 5, and 6, respectively. Figure 6(e1) and (e2) are for N=6 and two different initial conditions instead. For N=3, remarkably the distribution begins to have a nonuniform form; one oscillator displays a different death amplitude to the others [Fig. 6(b)]. For N=4, two different patterns [one, (\mathit{a},\mathbf{-}\mathit{a},\mathit{a},\mathbf{-}\mathit{a}), and the other, (\mathit{a},\mathbf{-}\mathit{a},\mathbf{-}\mathit{a},\mathit{a})] with the same uniform configurations can be found, as shown in Fig. 6(c). The value of d of the upper pattern is the same as that in Fig. 6(a) for N=2, as the pattern can be simply viewed as consisting of two N=2 patterns. For the lower pattern, the value of d can be easily calculated; d\approx 2.497. We also found that these two different patterns actually have extremely different generating probabilities: among {10}^6 times of observation, the upper pattern was generally found except only about 100 times for the lower one; thus, their appearing probabilities have {10}^4 times difference. For the pattern of five coupled oscillators in Fig. 6(d), the same as N=3 in Fig. 6(b), a nonuniform pattern is found again with three different values of d, and only one frozen pattern is observed. In contrast to this, for N=6 in Fig. 6(e1) and (e2), two coexisting patterns can be found, and both are predictable from the patterns of smaller-size systems. In particular, the pattern in Fig. 6(e1) is a combination of three N=2 patterns in Fig. 6(a), and the pattern in Fig. 6(e2) is a combination of two N=3 patterns in Fig. 6(b).

From these observations, we may summarize the properties of these frozen patterns as follows: First, these patterns are well organized with rich but well-defined structures. For example, all the uniform patterns, which are usually found for even N, possess the rotation symmetry, as shown in Fig. 6(a), (c) and (e1), while all the nonuniform patterns, possess the left-right mirror symmetry, as shown in Fig. 6(b), (d) and (e2). This may be solely due to the ring-structure network topology employed in our model. Therefore, the pattern of large systems can be viewed as being divided into several disconnected pattern pieces and each of the pattern pieces is a copy of pattern of smaller systems. For instance, see the pattern in Fig. 6(e2). Second, the final steady states depend on initial conditions, and for different initial conditions their probabilities can be extremely different. Third, as the system size increases, the pattern structure may become more complicated and the number of the coexisting patterns may become larger. We have studied the systems of size up to N=20 and well verified these essential characters. All these observations are expected to be helpful for much larger systems, as we will see in the next section.

4 Stationary patterns in large-size systems

Now we investigate the dynamics in large-size systems. As an example, we consider a system of size N=2000. The qualitative behavior does not change for any other N if it is large. As a first step, it is convenient to numerically compute the LLE with the change of {ϵ}_1 and {ϵ}_2. The result is shown in Fig. 1(a), with a quite similar pattern to that in Fig. 7(a) for N=2. This finding indicates that the appearance of OD patterns in coupled identical oscillators [Eq. (11)] is generic indeed and it is even independent of the system size.

To see the detailed structure of frozen patterns, we choose the same parameters ({ϵ}_1=-3, {ϵ}_2=0, \omega =2) and show a segment of typical patterns in Fig. 8. Compared to the patterns in small-size systems (Fig. 6), the pattern obviously shows much more complicated structures. Roughly, it consists of some uniform flat fragments, which have a constant value and are randomly distributed, and some embedded “bursting” states, which seemly scatter into the flat patterns in a random way. Taking a closer look at Fig. 8(a), (b) and (c), we can find that the flat state is exactly constructed by the small pattern (\mathit{a},-\mathit{a}) of N=2 system. To be clear, we use several vertical lines to mark and separate these uniform plateaus and nonuniform bursting structures. A further study shows that they are very sensitive to the choices of initial conditions and a small change of initial states of the oscillators could lead to a complete change of pattern structure.

To view and characterize these rich pattern structures in large-size systems, we randomly selected {10}^6 initial conditions from the region x_i,y_i\in [-1,1] and studied the final states after a long lapse. Altogether, six types of bursting patterns (highlighted by two red lines in each panel) are found in order of decreasing probability of appearance, as illustrated in Fig. 9. Roughly, the probability of each type of bursting patterns is determined by the width of these burstings. As the width increases, the probability quickly decreases. For instance, the probability of the pattern plotted in Fig. 9(a) is about {10}^3 times of that of Fig. 9(e2). Some other patterns are not presented here due to their much small probability and our limited simulation time. For the structures of these patterns, we find that all of them possess a left-right mirror symmetry. More specifically, if the bursting pattern has an odd size, then the ceasing states will be symmetrically distributed around the middle oscillator, e.g., the pattern in Fig. 9(b), whereas if the size is even, the ceasing states will be symmetrized around the middle link, e.g., the patterns in the others.

For the uniform parts in the patterns of d, we may characterize them in a statistical way. In Fig. 10(a), P_l, which represents the probability of finding uniform structures of a width of l, is plotted against l. In simulations, the system size is fixed at N=2000 again and each point is averaged over {10}^6 random initial conditions. A semi-logarithmic fit shows a clear exponential scaling between P_l and l, namely, P_l\propto {\mathrm{e}}^{-l/⟨l⟩}, with ⟨l⟩ defined as the average length of l. The probability P_l quickly decreases with the increase of l. Thus, even the uniform patten for whole oscillators (l=N) is stable, it cannot be observed from randomly chosen initial conditions due to its extremely small probability. This point can be believed as one of the key differences between the death patterns in large and small systems (comparing Figs. 8 and 6). From this fit, we obtain ⟨l⟩\approx 6.7. Obviously, this average length depends only on system parameters, and thus, it can be viewed as a characteristic length of the uniform patterns and further that of the whole patterns of the coupled systems. Furthermore, we scan one coupling parameter {ϵ}_1 with the other fixed ({ϵ}_2=0) and plot ⟨l⟩ versus {ϵ}_1 in Fig. 10(b).

To proceed further, we briefly study death patterns in two-dimensional coupled systems, which can be described as

\begin{array}{}\text{(12)}& \begin{array}{c}{\dot{x}}_{i,j}=p_{i,j}{x}_{i,j}-\omega y_{i,j}+{ϵ}_1(x_{i+1,j}+x_{i-1,j}+x_{i,j+1}+x_{i,j-1}-4x_{i,j})\hfill \\ {\dot{y}}_{i,j}=p_{i,j}{y}_{i,j}+\omega x_{i,j}+{ϵ}_2(y_{i+1,j}+y_{i-1,j}+y_{i,j+1}+y_{i,j-1}-4y_{i,j})\hfill \end{array}\end{array}

with i (i=1,2,\mathrm{\cdots },N) and j (j=1,2,\mathrm{\cdots },N) being the oscillator indexes on the horizontal and vertical directions, respectively. Periodic boundary conditions are considered.

Searching the extended two-dimensional parameter space ({ϵ}_1,{ϵ}_2), we find that OD state can also be generically found here; the phase diagram from the computation of the LLE with blank regions gives not much difference for the two-dimensional and one-dimensional systems [comparing Fig. 7(b) to (a)]. The region with several scattered black points in Fig. 7(b) simply exhibits bistable behaviors. As an example, a stationary pattern is given in Fig. 11(a), where points of amplitude d>3.3 are marked by black, and otherwise left for blank. A size of N\times N (N=100) is selected with the parameters {ϵ}_1=-1.5, {ϵ}_2=0.0, and \omega =2.0 within the OD regime. Now the pattern is well characterized by randomly distributed spots and it seems to be divided into a number of loosely connected modular patterns, forming various “pattern islands”. This phenomenon is robust for a small change of {\omega }_{i,j}. In Fig. 11(b), we randomly choose {\omega }_{i,j} within (1⩽{\omega }_{i,j}⩽3) and find a similar pattern. Note that all these patterns including the black and white regions are frozen.

5 Discussion and conclusion

In conclusion, we have comprehensively studied the dynamical behaviors in coupled limit-cycle oscillators in the extended two-dimensional coupling parameter space. Our finding uncovers a new regime for death states in coupled oscillator systems. Compared with the usual diagonal coupling form with equal strengths on the two Cartesian coordinates x and y, the only difference is that now we consider the case of unequal coupling strength. Thus, the usual diagonal in parameter space is extended to a whole two-dimensional parameter space, and the dynamical behaviors can become richer. The coupling form we consider is also very important, as it can be thought of as a simplified canonical model for any linearly coupled oscillators. We prove that in this case the usual AD really does not exist, and it may only occur in coupled nonidentical systems with a sufficiently large system mismatch. On the contrary, the OD states can exist under very broad conditions and survive even in coupled nonidentical systems. It is our major contribution to give a novel coupling and find death states in coupled oscillators in the absence of any effects of parameter mismatch and delayed time.

For a pair of coupled oscillators, as two examples of regular couplings, both the diffusive and direct couplings are investigated theoretically and numerically in detail. Apart from the death state, we also find some other interesting dynamical behaviors, for example, the periodic behaviors of oscillators with the same phase or anti-phase relation, and the coexistences. We expect that these theoretic results are not only valuable for model study, but also give a suggestive clue for the experimental observations in coupled identical oscillators [15, 16].

For the phenomenon of death state in large systems, the pattern of the new-defined pattern amplitude d is constructed by two separated modular structures: one is the plateau, which is characterized by a constant d and an exponential distribution of width, and the other is the bursting, which is characterized by nonuniform values of d and rich spatial structures. These two modular structures are randomly distributed and intermingled with each other. Based on the facts that the global AD pattern cannot be observed in identical (or slightly non-identical) systems under the usual conditions (except that an extremely large coupling strength is maintained) and the global OD pattern can be observed generically in our systems, we may conclude that OD is more robust, compared to AD. In Ref. [14], on the study of coupled Brusselators, the author conjectured that the pattern of death state in large-size systems could be nonuniform, but no evidence was given. In the current study, we have not only provided solid evidences to support this conjecture, but also analyzed the general properties of the nonuniform patterns. Other than the studies on the coupled Landau-Stuart oscillators in the paper, we also investigated some other models, and even found OD in diffusively coupled systems with some small-world network connections and in diffusively coupled chaotic oscillators as well.

Finally, it is worthwhile to compare our OD patterns with the classical Turing patterns. In his famous paper in Ref. [44], Turing gave an explanation of biological morphogenesis, on the basis of rigorous mathematical analysis of the instability of steady states. The novel phenomenon predicted in the paper, the Turing pattern, has been proved to be very general and capable of being found in many disciplines. Although both of these two types of patterns show a coupling-induced steady state, they have quite different dynamical natures: the OD state comes from an originally oscillatory system and the oscillators exhibit oscillation in the absence of coupling, but the Turing pattern comes from an originally homogeneous steady state and without coupling the system is a stable fixed point. On the other hand, they really share some basic characteristics, both in their pattern structures and in their mechanisms. First, the Turing patterns show a nonuniform stationary structure, whose spatial period is determined by an intrinsic wavelength selected by the most unstable mode. Here, the width of uniform parts of OD patterns shows an exponential distribution and is characterized by a well-defined average width (or wavelength), which is also intrinsic and only determined by system parameters (not initial conditions). Thus, it is quite natural to find these Turing-like patterns in one-dimensional systems (Fig. 6) and two-dimensional systems (Fig. 11). Second, both of these patterns can be found in coupled identical (or nearly identical) systems. The condition of large parameter mismatches is not needed. Third, for the existence of Turing patterns it is necessary that the inhibitor must diffuse faster than the activator. Note that in a two-variable reaction-diffusion equation (or its discrete form, two-variable diffusively coupled oscillators), one variable, which may activate its own production, is called the activator, and the other variable, which may inhibit its own production, is called the inhibitor. Now in our coupled oscillator systems, the two x and y variables of Cartesian coordinates do not show this character; they play equal contributions to the coupled systems. However, it seems that unequal diffusive coefficients for x and y ({ϵ}_1\ne {ϵ}_2) are also needed for the existence of OD. See the plots in Figs. 1(a) and 7, where the regions of OD state are always away from the central forbidden zone (oscillation regions) near the diagonal. This point is also true for the observed OD in coupled chaotic oscillators. Fourth, it is well-known that the underlying mechanism of Turing stationary patterns comes from diffusion-induced instability (Turing bifurcation), which means that under some particular conditions, diffusion may even show a destructive effect. Here in our work, in diffusively coupled Landau-Stuart oscillators, the choice of coupling is also crucial: the appearance of OD needs at least one negative coupling strength ( or , see the phase diagrams). Usually, a positive coupling strength in diffusive coupling represents negative-feedback coupling, whereas a negative coupling strength denotes positive-feedback coupling. The same as negative-feedback coupling, positive-feedback coupling is also very general, such as in electrical engineering and biological systems [46, 47]. Therefore, similar to the activator-inhibitor mechanism in Turing patterns, here for the OD states, the coupled systems may need a sufficiently large positive-feedback coupling (i.e., a negative diffusive coupling strength) to induce the instability of the periodic oscillator and make the new-born steady state stable. It is also notable that for a positive equal diffusive coupling strength in Eq. (11), i.e., {ϵ}_1={ϵ}_2>0, only synchronization or splay states of oscillators can be found [48]. Nevertheless, this condition is not always necessary; in the coupled chaotic oscillators, we do observe OD under the usual condition of positive diffusive coupling. Last, we expect that all these findings in the paper are of significance and importance for nonlinear physics and possible applications, and can be justified in experiments in the future.

References

[1] A. T. Winfree, The Geometry of Biological Time, New York: Springer-Verlag, 1980

[2] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Berlin: Springer, 1984

[3] I. R. Epstein and J. A. Pojman, An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos, New York: Oxford University Press, 1998

[4] A. Pikovsky, M. Rosenblum, and J. Kurths , Synchronization: A Universal Concept in Nonlinear Dynamics, Cambridge: Cambridge University Press, 2001

[5] J. D. Murray, Mathematical Biology, 3rd Ed., Berlin: Springer, 2003

[6] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett., 1990, 64: 821

[7] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett., 1998, 80: 2109

[8] J. Z. Yang, G. Hu, and J. H. Xiao, Phys. Rev. Lett., 1998, 80: 496

[9] M. Zhan, G. Hu, and J. Z. Yang, Phys. Rev. E, 2000, 62: 2963

[10] M. Zhan, J. H. Gao, Y. Wu, and J. H. Xiao, Phys. Rev. E, 2007, 76: 036203

[11] W. Zou and M. Zhan, Europhys. Lett., 2008, 81: 10006

[12] M. Bennett, M. Schatz, H. Rockwood, and K. Wiesenfeld, Proc. Roy. Soc. London Ser. A, 2002, 458: 563

[13] W. Q. Liu, J. H. Xiao, X. Qian, and J. Z. Yang, Phys. Rev. E, 2006, 73: 057203

[14] K. Bar-Eli, Physica D, 1985, 14: 242

[15] K. Bar-Eli, Reac. Kinet. Catal. Lett., 1990, 42: 435

[16] M. Yoshimoto, K. Yoshikawa, and Y. Mori, Phys. Rev. E, 1993, 47: 864

[17] Y. Yamaguchi and H. Shimizu, Physica D, 1984, 11: 212

[18] R. E. Mirollo and S. H. Stogatz, J. Stat. Phys., 1989, 60: 245

[19] G. B. Ermentrout and W. C. Troy, SIAM J. Math. Anal., 1989, 20: 1436

[20] D. G. Aronson, E. J. Doedel, and H. G. Othmer, Physica D, 1987, 25: 20

[21] D. G. Aronson, G. B. Ermentrout, and N. Kopell, Physica D, 1990, 41: 403

[22] M. F. Crowley and I. R. Epstein, J. Phys. Chem., 1989, 93: 2496

[23] D.-S. Lee and J.-W. Ryu, Appl. Phys. Lett., 2005, 86: 181104

[24] M.-D. Wei and J.-C. Lun, Appl. Phys. Lett., 2007, 91: 061121

[25] G. G. Ermentrout, Physica D, 1990, 41: 219

[26] P. C. Matthews and S. H. Stogatz, Phys. Rev. Lett., 1990, 65: 1701

[27] E. Ullner, A. Zaikin, E. I. Volkov, and J. Garcia-Ojalvo, Phys. Rev. Lett., 2007, 99: 148103

[28] L. Rubchinsky and M. Sushchik, Phys. Rev. E, 2000, 62: 6440

[29] L. Rubchinsky, M. Sushchik, and G. V. Osipov, Math. Comput. Simul., 2002, 58: 443

[30] F. M. Atay, Physica D, 2003, 183: 1

[31] Z. Hou and H. Xin, Phys. Rev. E, 2003, 68: 055103(R)

[32] J. Yang, Phys. Rev. E, 2007, 76: 016204

[33] D. V. R. Reddy, A. Sen, and G. L. Johnston, Phys. Rev. Lett., 1998, 80: 5109

[34] D. V. R. Reddy, A. Sen, and G. L. Johnston, Physica D, 1999, 129: 15

[35] D. V. R. Reddy, A. Sen, and G. L. Johnston, Phys. Rev. Lett., 2000, 85: 3381

[36] A. Takamatsu, T. Fujii, and I. Endo, Phys. Rev. Lett., 2000, 85: 2026

[37] R. Herrero, M. Figueras, J. Rius, F. Pi, and G. Orriols, Phys. Rev. Lett., 2000, 84: 5312

[38] Y. Zhui, I. Z. Kiss, and J. L. Hudson, Phys. Rev. E, 2004, 69: 026208

[39] A. Prasad, Y. C. Lai, A. Gavrielides, and V. Kovanis, Phys. Lett. A, 2003, 318: 71

[40] K. Konishi, Phys. Rev. E, 2003, 68: 067202

[41] R. Karnatak, R. Ramaswamy, and A. Prasad, Phys. Rev. E, 2007, 76: 035201

[42] M. Y. Kim, R. Roy, J. L. Aron, T. W. Carr, and I. B. Schwartz, Phys. Rev. Lett., 2005, 94: 088101

[43] F. Takens, Detecting Strange Attractors in Turbulence, Lecture Notes in Mathematics, New York: Springer, 1981: 366–381

[44] A. M. Turing, Philos. Trans. Roy. Soc. London Ser. B, 1952, 237: 37

[45] M. Kawato and R. Suzuki, J. Theor. Biol., 1980, 86: 547

[46] K. J. Astrom and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers, Princetion: Princeton University Press, 2008

[47] J. Bechhoefer, Rev. Mod. Phys., 2005, 77: 783

[48] W. Zou and M. Zhan, unpublished

Rights & permissions

© Higher Education Press and Springer-Verlag 2009