Multistability in nonlinearly coupled ring of Duffing systems

In this paper we consider dynamics of three unidirectionally coupled Duffing oscillators with nonlinear coupling function in the form of third degree polynomial. We focus on the influence of the coupling on the occurrence of different bifurcation’s scenarios. The stability of equilibria, using Routh-Hurwitz criterion, is investigated. Moreover, we check how coefficients of the nonlinear coupling influence an appearance of different types of periodic solutions. The stable periodic solutions are computed using path-following. Finally, we show the two parameters’ bifurcation diagrams with marked areas where one can observe the coexistence of solutions.


Introduction
In recent years, the dynamics of coupled nonlinear oscillators become one of widely investigated topic with a lot of applications in different field of science. In networks of coupled systems we can distinguish two main coupling schemes. The first one is a global coupling, where all-to-all coupling is present, hence the oscillators are connected with each other by direct links or by the mean field [1][2][3][4]. The second type is a local coupling, where the single oscillator is connected to nodes in its nearest neighborhood [5][6][7][8][9][10]. In this work we consider the local coupling. Let us now present now the classification of possible local coupling schemes. The first indicator is the direction of signal. The interchange of signal can be bidirectional (both coupled nodes exchange informations). This type of coupling occurs in mechanical and structural systems, where the action always precedes reaction. Nevertheless, it is also present in electrical, biological and social systems [11][12][13][14][15]. Contrary to aforementioned scheme, the signal can be transmitted from one system to another without receiving the feedback. This type of coupling is called unidirectional. It is common for lasers, neurons and electrical systems [16][17][18][19].
In [20] authors take into account three coupled Lorenz systems. They show the existence of rotating waves, visible as a phase relation between coupled systems. In case of simplest periodic rotating waves all coupled systems follow same trajectory with constant, uniformly spaced phase shifts. This relation preserves even in chaotic regime (chaotic rotating waves). Afterwards, a detailed study of unidirectional coupled Stuart-Landau oscillators was provided in [21,22]. The analysis describe an analytical method to investigate the behavior of equilibrium and periodic solutions. In our recent paper we consider dynamics of non-linearly coupled Toda oscillators [23]. The destabilization of steady state may occur by sub-or supercritical Hopf bifurcation, while periodic solutions lose stability via period doubling, Neimark-Saker or saddlenode bifurcations. Hence, based on type of coupling function we can enforce different types of systems' dynamics.
In this paper we show the dependence of bifurcation scenarios on the coefficients of coupling. As a node system we take a Duffing systems and we create the ring consists of the single, identical systems. The Duffing oscillators are coupled via nonlinear function. We show that proper choice of coupling coefficients may lead to given bifurcations scenario.
The paper is organized as follows: in Sect. 2 we introduce the investigated system. Section 3 is devoted to stability analysis of steady states. The path following of periodic solutions is shown in Sect. 4. The results presented in paper are concluded in Sect. 5.

Model of the system
We consider three autonomous Duffing oscillators. The subsystems are unidirectionally coupled with left neighbor of each. The system is given by following first order ODEs: where j = 1, 2, 3, α > 0, d > 0 and γ ≥ 0 is a coupling coefficient. As the oscillators are linked in a ring, x j−1 should be thought of as x 3 for j = 1. Coupling function is nonlinear and is given by following third degree polynomial: In our work we investigate the influence of parameters k 1 , k 2 and k 3 on bifurcations' scenarios present in the system.

Analysis of steady states
In this section we investigate the existence and stability of equilibria. We show the complete analytical analysis of scenarios of stabilization and destabilization of steady states. We present how parameters influence the appearance of different types of bifurcations. Moreover, we use l'Hopital's rule to find asymptotes for equilibria and we apply Routh-Hurwitz criterion to investigate stability.
Eigenvalue analysis for x 0 . In the case of equilibrium x 0 we perform the stability analysis analytically. Eigenvalues of this steady state are given by: where j = 1, 2, 3. This formula remains true for every other coupling function, for which C (0) = k 1 = 0. Coefficients k 2 and k 3 do not influence stability of x 0 equilibrium, because they do not appear in the linearized system. From λ + j (γ) only λ + 3 (γ) can be equal to zero, because only then exp 2πji for which stationary bifurcation occurs. Now, we will show that the real parts of the rest of eigenvalues are negative. The characteristic polynomial for Jacobian matrix of system (1) for equilibrium x 0 at point γ S is given by: Equating P (η) to zero, we get that η 1 = 0, η 2 = −d and η 3,4,5,6 are the roots of fourth degree polynomial: In order to determine if polynomial (3) is stable (i.e. it's roots have negative real parts), we use Routh-Hurwitz stability criterion, which implies that polynomial: For polynomial (3): Δ 1 = 3αd > 0, Δ 2 = 3αd(d 2 + α) > 0, Δ 3 = Δ 4 = 3d 2 α(2d 2 − α). Both Δ 3 and Δ 4 are positive real numbers only if α < 2d 2 . In conclusion, we showed that all real parts of eigenvalues for equilibrium x 0 at point γ S are negative (despite the one equal zero responsible for bifurcation) unless α > 2d 2 . Let us assume that α < 2d 2 . Then, as λ + 3 (γ) is real valued function, we examine its derivative obtaining that λ + 3 (γ) is strictly increasing function. In a consequence at point γ = γ S λ + 3 (γ) change its sign from negative to positive. Moreover as all the other λ ± j (γ) are continuous functions, we know that in a certain neighborhood of γ = γ S their real part is negative. Hence x 0 looses its stability at γ = γ S .
The other occurring bifurcation for steady state x 0 is the Hopf bifurcation, in which the small-amplitude limit cycle is born. Conditions Re(λ ± j ) = 0 and Im(λ ± j ) = ω > 0 for Hopf bifurcation are fulfilled only for λ + 1 and γ = γ H . We search the value of γ H by solving equation λ + 1 = iω and obtaining that γ Conjugate condition Re(λ ± j ) = 0 and Im(λ ± j ) = −ω is fulfilled by λ + 2 and γ = γ H . Defining a function B(d, α) = γ H − γ S and analyzing its sign, one may find which bifurcation, stationary or Hopf, occurs as first one. If α < 2d 2 the bifurcation in which x 0 loses its stability is stationary bifurcation. If α > 2d 2 than the destabilization is caused by the Hopf bifurcation. This confirms the earlier results obtained from Routh-Hurwitz stability criterion.  (1): (4): In the paragraphs: "Existence of x 01 and x 02 ", "Asymptote" and "Condition for pitchfork bifurcation" we divide the space of parameters into following areas (cases): Existence of x 01 and x 02 . Equilibrium points x 01 and x 02 occur only if the function , is greater than zero, because only then x 01 and x 02 are real numbers.
We will investigate the behavior of f (γ), which is quadratic function of the variable γ. This way we will find out for which parameters γ equilibrium points x 01 and x 02 exist.
Let us notice that: γ 1,2 are always real numbers, because the term under the root is always positive: The fact that γ 1,2 are always real numbers means that there always exists a bifurcation parameter for which equilibrium points occur.
Investigating the behavior of function f (γ) we consider three cases concerning the sign of denominator of expressions for γ 1,2 (Case (1-3)) and the case when γ 1 = γ 2 (Case (4)). The shape of function f (γ) depending on aforementioned cases is presented in Fig. 1.
α and γ = α k1 are trivial and equal Asymptote. For equilibria x 01 and x 02 there might occur an asymptote, what means that one or two equilibria escape to infinity (blow up). In this section we present under what conditions and for which equilibrium asymptote occur. If asymptote is present, it exists for γ = 1/k 3 . Let us assume that k 2 > 0. Calculating the limit of x 01,02 with γ tending to 1/k 3 we use l'Hopital's rule, obtaining: Assuming k 2 < 0. Then If k 2 = 0, then Now, we examine under which conditions asymptote occurs. The behavior of the system is investigated for γ ∈ [0, ∞), that is why for k 3 ≤ 0 we do not obtain asymptote. Assuming that k 3 > 0, we examine three cases: We will show that asymptote always occurs in that case, what means that 1 k3 > γ 2 .
Both sides are positive real numbers, hence we may square both sides of inequality and after simple transformations we obtain: In this case the existence of asymptote is equivalent to fulfilling condition γ 2 < 1 k3 < γ 1 . Hence: k3 > 0. Hence in this case asymptote occurs for γ ∈ (γ 2 , γ 1 ) unless k 2 = 0. For k 2 = 0 asymptote occurs for γ 1 or γ 2 .
In order to square both sides of equation we need to assume that (k 3 α − k 1 ) > αk 2 2 2k 1 .

Periodic solutions
In the numerical simulations we choose α = 0.1, d = 0.3, so α < 2d 2 and according to Sect. 3 "Eigenvalue analysis for x 0 " the first bifurcation destabilizing equilibrium x 0 is stationary one. The stability of equilibria x 01 and x 02 is analyzed numerically by investigation of the signs of maximum real parts of eigenvalues. We present them in Fig. 4(a-d) on two-parameter diagrams. The black area indicates coexistence of stable x 01 and x 02 , the coral stands for stable x 01 and yellow one for stable x 02 . When the region is grey both equilibria are unstable then. Finally, white area means that equilibria do not occur. The regions, where stable periodic orbits occur are indicated by hatched areas. The bifurcation curves presented in this section are obtained using Auto-07p: Continuation and bifurcation software [24,25].
Equilibria calculated in Sect. 3 are used as the base to the calculation of periodic orbits, as all periodic orbits are created via Hopf bifurcations from equilibria. Figures 4(a-d) connects equilibria analysis, i.e. Sect. 3 with periodic orbit analysis (Sect. 4).
In each Figs. 4(a-d) the stationary bifurcation, where both x 01,02 are born from x 0 is indicated by the line l S . All of obtained periodic orbits are created via supercritical Hopf bifurcations l HB , in which x 01 or x 02 lose stability. But not only supercritical Hopf occurs in investigated system, but also subcritical one. In Fig. 4(a) in subcritical Hopf bifurcation (l SUB−HB line) stable x 02 appears. In Fig. 4(c) we may see that x 01 loses its stability in the subcritical Hopf bifurcation, since no stable periodic orbit is born. Now, we know that each equilibria may undergo Hopf bifurcation. In analytical part we obtained that x 0 undergoes a Hopf bifurcation at γ H (for the parameters α and d we chose this bifurcation occurs for unstable x 0 ). Hopf bifurcation of x 01 and x 02 is obtained in numerical studies. Finally, whole area where equilibria x 01,02 exist is bounded by l LP line, standing for fold bifurcation.
In Fig. 4(a) parameter k 1 = 1 and k 2 = 3 are fixed and k 3 is varied. Increasing coupling parameter γ we first meet Hopf bifurcation line, which coincides with the value of γ 2 . Created orbit loses it's stability at period doubling bifurcation for k 3 ∈ (−2, −0.534). For k 3 = −0.534 codimension-two bifurcation (the Neimark-Sacker-period-doubling) occurs. For k 3 ∈ (−0.534, 1.77606) periodic orbit bifurcates via the Neimark-Sacker bifurcation. This scenario is different from k 3 = 1.77606 (fold-Neimark-Sacker bifurcation) to k 3 = 2, where the periodic orbit loses it's stability in saddle-node bifurcation. Then the stable periodic solution appears, when x 01 becomes unstable in the next Hopf bifurcation. This periodic orbit loses its stability in period doubling or Neimark-Sacker bifurcation. The Neimark-Sacker-period-doubling bifurcation at k 3 = −1.753503 is the border between those two bifurcations. Another periodic orbit is created from x 02 and is stable until it does not bifurcates via the Neimark-Sacker bifurcation.
The last Fig. 4(d) shows behavior of the system in the case of pitchfork bifurcation of equilibria. For k 3 ∈ (0, 4.44924) the symmetrical limit cycles , created from symmetrical equilibria x 01 and x 02 , destabilize via fold bifurcation. For k 3 ∈ (4.44924, 10) loss of stability occurs due to the Neimark-Sacker bifurcation. When k 3 ∈ (10, 20) equilibria are unstable and none stable periodic orbit is born. Figures 5, 6 and 7 present bifurcation scenarios of equilibria and limit cycles in one-parameter plots, which are intersections of diagrams in Figs. 4(a-d).
The continuation in Figs. 5(a-d) corresponds to the line k 3 = −2 in Fig. 4(a) and the bifurcations are confirmed to occur in that particular order. The interesting detail in this scenario is heteroclinic bifurcation. In Fig. 5(b) one may see that bifurcation curve undergo many fold bifurcations before ending at a heteroclinic bifurcation. For γ = 0.80852 period of the limit cycle tends to infinity and that is the point where bifurcation occurs (Fig. 5(c)). The phase portrait of considered orbit is presented in Fig. 5(d), we marked by dots unstable equilibrium (x 1,2,3 = −0.2795, −0.0115, −0.0839) which is connected by heteroclinic orbit. Branch of steady states emerging from this equilibrium is mainly unstable. With increasing of γ the stability is not changed, but with decreasing of γ we observe a saddle-node bifurcation for γ = 0.6436, where the stable part appears but it destabilizes in Hopf bifurcation for γ = 0.6441 (it is stable in very narrow range). The Hopf bifurcation is supercritical, hence the stable periodic appears. For slightly larger value of γ it loses stability in Neimark-Sacker bifurcation (γ = 0.6443). We can reach those solutions only with very precise initial conditions, hence from practical point of view they importance is negligible and we decided to do not present them in our analysis.
indent Figure 6(a) is the intersection of Fig. 4(a) for k 3 = 2. Equilibrium x 0 is stable until γ = 0.1, where in transcritical bifurcation it loses stability and stable x 01 is created. With increase of bifurcation parameter, we reach point γ = 0.69335, where in Hopf bifurcation stable limit cycle is created, which becomes unstable in Neimark-Sacker bifurcation at γ = 0.96989. Following unstable x 01 from γ = 0.1 with decreasing γ, we reach limit point at γ = 0.081935, where stable x 02 appears. Created x 02 loses stability at γ = 0.082415 at Hopf bifurcation. Appeared stable limit cycle becomes unstable at the limit point at γ = 0.12235.
Dynamics of system for parameters k 1 = 1, k 2 = 1, k 3 = 2 is presented in Fig. 6(b), which is intersection of Fig. 4(b) for k 2 = 1. Solution x 0 loses stability  at γ = 0.1 in transcritical bifurcation. Newly created x 01 is stable with increasing γ and unstable with decreasing γ. The stable branch of x 01 becomes unstable at γ = 0.25875, where in Hopf bifurcation stable periodic orbit appears. Following period orbit we reach the limit point at γ = 0.26117, where it lose stability. Following unstable periodic solution we find another Hopf bifurcation at γ = 0.14358, in which periodic orbit is born from equilibrium x 0 . The unstable branch of x 01 (i.e. x 01 for γ < 0.1) exists until limit point is reached at γ = 0.097076. Then, stable x 02 is created, which preserves it's stability till γ = 0.11676, where in Hopf bifurcation stable periodic orbit is born. The stable periodic solution exists till γ = 0.15321, where another, unstable, solution appears.
at the pitchfork bifurcation and two stable equilibria x 01 and x 02 are born. Both x 01 and x 02 destabilize at γ = 0.11896 in Hopf bifurcation, where the stable limit cycles are created. The both limit cycles lose stability in the Neimark-Sacker bifurcations at γ = 0.13394. Following unstable branches we reach limit point at γ = 0.13466. With further continuation we reach symmetry breaking at γ = 0.13411, where two unstable branches connect.

Conclusions
In this paper we show the variety of bifurcation scenarios, which depend on coefficients in nonlinear coupling. Using analytical approach, we manage to investigate the existence and stability of equilibria. We show that in considered system we can observe changes in stability of steady states via the Hopf, pitchfork and transcritical bifurcations. Appearance of each type of bifurcation is dependent on the parameters of the coupling function. The stable equilibria coexist in the wide range of coupling's parameters.
With increasing the coupling parameter γ we observe destabilization of equilibria and appearance of periodic solutions via sub-or supercritical Hopf bifurcations. We merge results obtained in analytical and numerical calculations to show ranges of existence of stable and unstable steady states and stable periodic solutions. The destabilization of periodic orbits is also strongly dependent on parameters of the coupling function. With changes of coefficients k 1 , k 2 , k 3 we can observe following types of bifurcations: Hopf, Neimark-Sacker and saddle-node.
In summary, we show that with changing parameters of non-linear coupling function we can strongly affect the bifurcation scenario of the system and reach given sequence of bifurcations. Additionally, one can move between regions in parameter space where the multi-stability exists or does not exist. This work has been supported by Lodz University of Technology own Scholarship Fund (PJ).