On the stability analysis of a pair of van der Pol oscillators with delayed self-connection , position and velocity couplings

In this paper, we perform a stability analysis of a pair of van der Pol oscillators with delayed self-connection, position and velocity couplings. Bifurcation diagram of the damping, position and velocity coupling strengths is constructed, which gives insight into how stability boundary curves come into existence and how these curves evolve from small closed loops into open-ended curves. The van der Pol oscillator has been considered by many researchers as the nodes for various networks. It is inherently unstable at the zero equilibrium. Stability control of a network is always an important problem. Currently, the stabilization of the zero equilibrium of a pair of van der Pol oscillators can be achieved only for small damping strength by using delayed velocity coupling. An interesting question arises naturally: can the zero equilibrium be stabilized for an arbitrarily large value of the damping strength? We prove that it can be. In addition, a simple condition is given on how to choose the feedback parameters to achieve such goal. We further investigate how the in-phase mode or the out-of-phase mode of a periodic solution is related to the stability boundary curve that it emerges from a Hopf bifurcation. Analytical expression of a periodic solution is derived using an integration method. Some illustrative examples show that the theoretical prediction and numerical simulation are in good agreement.


I. INTRODUCTION
Many important physical, chemical and biological systems such as semiconductor lasers, 1,2 coupled Brusselator models 3,4 and neural networks for circadian pacemakers 5 are composed of coupled nonlinear oscillators.Ubiquitous in nature due to finite propagation speeds of signals, time delay may have profound effects on the collective dynamics of such systems.7][8][9] In particular, the van der Pol oscillator has been considered by many researchers as the nodes for various networks.Atay 10 investigated the effect of delayed feedback for the van der Pol oscillator on oscillatory behavior.Maccari 11 investigated the resonance of a parametrically excited van der Pol oscillator under state feedback control with a time delay.From the viewpoint of vibration control, they demonstrated that the time delay and the feedback pairs could enhance the control performance and reduce the amplitude peak.For a van der Pol-Duffing oscillator with delayed position feedback, Xu and Chung 12 showed that time delay might be used as a simple but efficient switch to control motions of a system: either from orderly motion to chaos or from chaotic motion to order for different applications.Wirkus and Rand 13 studied the dynamics of two weakly coupled van der Pol oscillators with delayed velocity coupling due to its relevance to coupled laser oscillators.They found that both the in-phase and out-of-phase modes were stable for delays of about a quarter of the uncoupled period of the oscillators.Li et al. 14 extended the above work by including both delayed position and velocity coupling.They showed that, for the case of 1:1 internal resonance, both the in-phase mode and out-of-phase mode existed when the two coupling coefficients were identical, and there were two death domains when these two modes did not exist.Zhang and Gu 15 considered the dynamics of a system of two van der Pol equations with delay position coupling.They showed the existence of stability switches and, as the delay is varied, a sequence of Hopf bifurcations occurred at the zero equilibrium.Song 16 investigated the stability switches of two van der Pol oscillators with delay velocity coupling, and obtained different in-phase and anti-phase patterns as the coupling delay was increased.In the above papers, only weakly nonlinear van der Pol oscillators were investigated.
One of the most interesting and important collective behaviors in coupled oscillators that have aroused much attention in recent years is the amplitude death which refers to the diffusivecoupling-induced stabilization of unstable fixed points in coupled oscillators. 17,18 he theoretical and practical meanings of the phenomenon of amplitude death in coupled systems are of great significance.For example, it is a desirable control mechanism in cases such as coupled lasers where it leads to stabilization 19,20 and a pathological case of oscillation suppression or disruption in cases like neuronal disorders such as Alzheimers disease, Parkinsons disease, etc. [21][22][23][24][25] For the occurrence of amplitude death, one of the following conditions is needed: the parameter mismatch, 18,26 the time-delayed coupling, 27 dynamical coupling 28 and conjugate coupling. 29Amplitude death by delay was first reported by Ramana Reddy et al. 27 in their study of a pair of limit cycle oscillators which were the normal form for the Hopf bifurcation.A novel result 30 that they found was the occurrence of amplitude death even in the absence of a frequency mismatch between the two oscillators.Based on the system studied in Ref. 30, Song et al. 31 gave more detailed and specific conditions on the existence of amplitude death for different delays.Since Pyragas 32 introduced a novel feedback control method of using a single time delay some two decades ago, the investigation has been extended to multiple time delays 33 for stabilizing steady states of various chaotic dynamical systems.For several well-known chaotic systems, Ahlborn and Parlitz 34 showed that multiple delay feedback control is more effective for fixed point stabilization in terms of stability and flexibility, in particular for large delay times.Blyuss et al. 35 applied delayed feedback control to stabilize an unstable steady state of a neutral delay differential equation.They showed that a number of amplitude death regions came into existence in the parameter space due to the interplay between the control strength and two time delays.For the van der Pol-Duffing system with delayed position feedback, Xu et al. 36 obtained an amplitude death region near a weak resonant double Hopf bifurcation point.A thorough review of the current works can be found in Ref. 37.
The effect of time delay on the amplitude death of the Stuart-Landau oscillators with linearly (diffusively) delay coupling has been well studied. 27,37 s for the van der Pol oscillators with delay coupling, investigations have been focused only on the weakly nonlinear situation.For the pair of van der Pol oscillators with delay velocity coupling studied in Ref. 16, amplitude death is possible only when the damping strength is less than 0.5.Complex dynamics such as periodic-doubling sequences leading to chaos occur for strongly nonlinear situation. 38An interesting question naturally arises: can the strongly nonlinear van der Pol oscillators be stabilized using delay coupling?In other words, is it possible to derive a delay feedback control strategy such that amplitude death exists for all positive values of the damping strength?Ahlborn and Parlitz 34 suggested that more delays entering into the control terms were more effective and flexible for fixed point stabilization.Due to the complicated analytical expressions, the analysis of systems with multiple delays is always based on numerical simulations.It would be invaluable to develop an efficient analytical method for problems with multiple delays.Delayed position and velocity feedbacks are two kinds of strategies commonly used for control purposes.However, there are very few investigations on using both kinds for the stability control of the van der Pol oscillators.These situations constitute the motivation of the present paper.
In this paper, our goal is to derive a delay feedback control strategy for the amplitude death of nonlinear van der Pol oscillators with arbitrary large damping strength and investigate periodic solutions of the in-phase and out-of-phase modes arising from Hopf bifurcation.In doing so, we introduce three kinds of feedbacks namely, position, velocity, self-connection and three time delays.The paper is organized as follows.Sec.II describes the model formulation and introduces a parameter γ for finding the solutions of the characteristic equation of the linearized system.Local stability analysis of the zero equilibrium is preformed in Sec.III.Bifurcation diagram of the system parameters is constructed and the properties of the stability boundary curves in the plane of time delays are discussed.In Sec.IV, we turn our attention to amplitude death region in the plane of time delays and prove a sufficient condition for its existence for arbitrary large damping strength.Sec.V is devoted to the investigation of the periodic solutions of the in-phase and the out-of-phase modes.We show how the type of mode is related to the stability boundary curves obtained in Sec.III.An integration method is employed to study the amplitude of periodic solutions arising from a Hopf bifurcation.As illustrated in Sec.VI, the analytical results are in good agreement with those obtained from numerical simulations.Finally, Sec.VII contains the conclusions.

II. MODEL FORMULATION
We study a pair of van der Pol oscillators in which there are distinct, discrete time delays in the signal transmission of self-connection, position and velocity couplings.The coupled system is shown schematically in Fig. 1 and expressed in two delay differential equations as ( where τ 1 , τ A and τ B are, respectively, the time delays of self-connection, from oscillator x 2 to x 1 , and from x 1 to x 2 ; α and β are, respectively, the position and velocity coupling strengths; μ is the damping strength which is always positive.Let y 1 (t) = ẋ1 (t) and y 2 (t) = ẋ2 (t), system (1) can be written as We note that system (2) is reduced to the system investigated in Ref. 15 if β = τ 1 = 0 and that of Ref. 16 if α = τ 1 = 0.The characteristic equation of the linearization of system (2) at the origin is given by where Equation (3a) determines the local stability of the origin in (1).When τ 1 = τ 2 = 0, (3a) has the following four roots: Since μ > 0, the origin is always unstable.To investigate the distribution of roots of (3a) for nonzero τ 1 and τ 2 , we have to consider the stability boundary curves of (3b) and (3c) by letting λ = iω for ω > 0. From (3b) and (3c), we obtain, respectively, and where Eliminating τ 2 in (4a) or (4b), we arrive at the same equation and Then, after simplification, (5) becomes It follows from (7) that 0 < γ ≤ 1. Substituting ( 6) and ( 7) into (4a), we obtain To simplify the subsequent calculations, we let α = kβ.It follows from ( 7) and ( 8) that 1 (iω, For the solutions of (3b), we may first choose γ ∈ (0, 1] and obtain ω if it exists from (9a).Then, τ 1 and τ 2 can be obtained from (9b) and (9c), respectively.

III. LOCAL STABILITY ANALYSIS
In this section, we investigate the number of positive solutions of ω in (3a) and construct bifurcation diagram in the parameter space of μ, β and k.The van der Pol oscillator is inherently unstable at the zero equilibrium.A study of positive solution of ω gives the condition of μ, β and k that a pair of eigenvalues of (3a) cross the imaginary axis, resulting in a change of stability at the zero equilibrium.For a given k, we will show that the bifurcation diagram (μ, β) is partitioned into three regions according to whether the stability boundary curves exist (closed or open-ended) or not.For a region with no positive solution, the zero equilibrium is always unstable even when delays exist.The existence of two positive solutions in a region may lead to the occurrence of amplitude death, i.e. the stabilization of the zero equilibrium. Let ) be the discriminant of (9a).Since (9a) is a quadratic equation in ω 2 , we have the following lemma regarding the number of positive roots of ω in (9a).
To investigate how varies with respect to μ > 0 and k, we consider the equation = 0 and obtain With some calculations on inequalities, we obtain the regions in the (k, μ) plane with 0, 1 and 2 positive real roots of γ in (11a) as follows: |k| .Then, (C 1 ) − (C 3 ) are the conditions for 0, 1 and 2 positive real roots, respectively, of γ in (11a) (see Fig. 2).
Next, we consider the bifurcation diagram of (9a) with respect to the number of positive roots in ω for γ ∈ (0, 1].Define so that 1 2|k| and together with Lemmas 1 and 2, we have the following conditions of γ 1 , γ 2 and γ + (which is defined in (11b)) on the number of positive roots of ω in (9a): For γ ∈ (0, 1], we define γ = sin c for 0 < c ≤ π 2 .Then, for a positive solution of (9a), (τ 1 , τ 2 ) in (9b) and (9c) can be expressed as where m, m 1 = m − n ≥ 0. Since (10a) is the same as (9a) and for a positive solution of (10a), (τ 1 , τ 2 ) in (10b) and (10c) is given by In the above consideration, we use γ as a parameter to find the purely imaginary roots of (3a).The following theorem gives the bifurcation diagram of (3a) with respect to the number of purely imaginary roots.i is greater/less than one.These curves intersect concurrently at Theorem 1.Let k, β ∈ R, μ ∈ R + and define the following regions in the (μ, β) plane as (c) If k, β and μ satisfy the condition of region III a , then (3a) has 0 or 1 pair of purely imaginary roots according to γ ≤ γ 1 , or γ 1 < γ ≤ 1, respectively.
Proof.See Appendix A.
Theorem 1 provides information of and construction method for the stability boundary curves in the (τ 1 , τ 2 ) plane.
Remark 2 For the values of k, β and μ in region II, stability boundary curves can be constructed using (13a)-(13d) for γ ∈ [γ + , 1].Since ω is always finite and so do τ 1 and τ 2 , the curves form closed loops which may intersect itself (see Fig. 4).When the values of μ, β and k tend to either OP + or OP − such that γ + = 1, i.e. c = π 2 , the loops shrink to the set of discrete points  where m, m 1 ≥ 0 and τ 1 , τ 2 > 0.
Remark 3 For the values of k, β and μ in regions III a or III b , one of the roots of (9a) vanishes at γ = γ 1 .Then, both τ 1 and τ 2 in (13a)-(13d) become infinite.Therefore, the stability boundary curves are open-ended (see Fig. 5).
For the parameters of μ, β and k in region II, stability boundary loops in the (τ 1 , τ 2 ) plane can be constructed in the following way: (a) Set γ = γ + in (9a).From (9a) and (11b), ω = ω + is given by From numerical simulation, we observe that two horizontally neighboring loops never intersect each other to form amplitude death region.However, it may occur to two vertically neighboring loops with small m and m 1 (see Fig. 4).For m = 0, the small loops are initially born in the left-half (τ 1 , τ 2 ) plane and, as |β| increases, part of their interior crosses the τ 2 -axis into the right-half (τ 1 , τ 2 ) plane.In fact, in Fig. 4 where μ = 0.05, (3a) has two pairs of eigenvalues in the right-half plane when (τ 1 , τ 2 ) is outside the loops.Inside each loop (but not inside an overlapping area), (3a) has a pair of eigenvalues in the right-half plane.Inside the overlapping areas, all the eigenvalues of (3a) are in the left-half plane and thus these areas are amplitude death regions.To find the possibility of amplitude death region for arbitrary large μ, we investigate the intersection of stability boundary curves with m = 0 and the τ 2 -axis.
We first find the value of β for which a stability boundary curve with m = 0 is tangential to τ 2 -axis.When τ 1 = 0, (5) is reduced to At a tangential point, ( 16) has a double root in ω 2 such that the discriminant vanishes, i.e.
The larger root β T of ( 17) and the double root of ω 2 in ( 16) are given by, respectively, and In the stability control of the coupled van der Pol system for arbitrary large value of μ, we may assume that k is positive.At β = β T , it follows from (4b) that and Let θ = θ T at β = β T .It follows from ( 6) that θ T is in the second quadrant.Therefore, the tangential point is generated from (13b) or (13d) since π /2 < θ T = π − c < π.For β > β T , ( 16) has two positive roots ω 1 and ω 2 where ω 1 < ω T < ω 2 and each stability boundary curve intersects the τ 2 -axis at two points.From (13a)-(13d) with τ 1 = m = 0, we denote the intersection points by and θ i is obtained from (6) with ω = ω i .For delay-coupled van der Pol oscillators with one delay, it was shown in Ref. 15 and 16 that stability switchings occur along the time-delay axis.For the generalized system (3a), we apply the results of Ref. 39 which studied the stability boundary curves of general linear systems with two delays.The following theorem states the directions of crossing the imaginary axis for solutions of i (λ, τ 1 , τ 2 ) = 0 (i = 1, 2) at these intersection points in the positive τ 2 direction.
all the roots of (λ, 0, τ 2 ) = 0 have negative real part for 2 ).Thus, the interior of the intersection of the two stability boundary curves is an amplitude death region (see the shaded areas in Figs. 4 and 5).

IV. EXISTENCE OF AMPLITUDE DEATH REGION FOR ARBITRARY DAMPING STRENGTH
The system investigated in Ref. 16 is a special case of system (2) in that α = τ 1 = 0 and amplitude death is possible only when μ < 0.5.By adding a position delay feedback (i.e.α = 0) to the system in Ref. 16 while keeping τ 1 = 0, is it possible to obtain amplitude death for arbitrary large μ > 0? We give an affirmative answer.
To find sufficient conditions in terms of μ, k, β which satisfy the second condition in (22) we note that Given μ > 0, the condition

2
. Substituting ( 6) into ( 23) and after simplification, we obtain the following sufficient condition for the existence of amplitude death region.
Theorem 3 provides a sufficient condition for the existence of amplitude death region based on (23).To obtain a more accurate region in Fig. 6 for the existence of amplitude death region, we may use the second condition of (22) as For a given value of k 1 , the minimum value of k 2 for which ( 25) is satisfied can be computed numerically.The region is also plotted in Fig. 6.The shaded region is inside the region to the right of the dashed curve since condition ( 23) is more conservative than (25).It is observed that the boundary curve from numerical result is asymptotically close to that of the second condition of (24) as k 1 → ∞.

V. IN-PHASE/OUT-OF-PHASE MODES AND AN INTEGRATION METHOD
In the study of two weakly coupled van der Pol oscillators with delayed velocity coupling, Wirkus and Rand 13 found that both the in-phase and out-of-phase modes were stable for delays of about a quarter of the uncoupled period of the oscillators.In the bifurcation analysis of system (1) with delayed position and velocity couplings, we are going to investigate how the in-phase and out-of-phase modes occur and derive an analytical expression for a periodic solution arising from Hopf bifurcation using an integration method.
In Ref. 13, an algebraic condition was given to distinguish an in-phase oscillation from an outof-phase oscillation.In our investigation, we describe how in-phase and out-of-phase oscillations are related to the stability boundary curves in the (τ 1 , τ 2 ) plane.Therefore, we consider the oscillations from a geometrical point of view which is different from that of Ref. 13.
Next, we derive an analytical expression for a periodic solution arising from Hopf bifurcation using the integration method described in Ref. 36.and w(t) = w(t + 2π/ω), then

Theorem 5. If w(t) is a periodic solution of the adjoint equation of (26) which is given by
Proof.Eq.( 28) can be obtained by multiplying both sides of (26) by w T , integrating with respective to t from 0 to the period 2π /ω and applying integration by part.
Based on the expression of (B1), a periodic solution of (26) for small can be considered to be a perturbation to (B1) as An analytical expression for the zero-order amplitude p 0 in (29b) can easily be obtained from Theorem 5. Now, a periodic solution of ( 27) can be expressed as w(t) = (q 1 cos(ω 0 t) + q 2 sin(ω 0 t))K , (30)   where q 1 and q 2 are independent constants and K is defined in (29b).Substituting (29a)-(29b) and ( 30) into (28), comparing coefficients of the term and noting the independence of q 1 and q 2 , we obtain two equations in p 0 and ω 1 .On solving the equations, we obtain the analytical expressions of p 0 and ω 1 as shown in Corollary 1.
The stability of a periodic solution arising from a Hopf bifurcation can be determined from the Floquet theory described in Ref. 36.

VI. NUMERICAL SIMULATION
In this section, we perform numerical simulations to illustrate the results obtained from previous sections.In particular, we consider the periodic solutions arising from Hopf bifurcation near the shaded death region shown in Fig. 4 where (μ, k, β) = (0.05, 0.05, 0.1).Let (τ 1 , τ 2 ) be the intersection point of the stability boundary loops (0, 0) and (0, 1).Since it satisfies both (9a)-(9c) and (10a)-(10c), we have (τ 1 , τ 2 ) = (0.8744, 1.4249).From Theorem 4, if τ 1 is increased from τ 10 = 0.8744 while keeping τ 2 unchanged at 1.4249, we would expect a stable periodic solution of the in-phase mode arising from a Hopf bifurcation out of the stability boundary loop (0, 0) and another stable periodic solution of the out-of-phase mode from the loop (0, 1).As predicted from theoretical consideration, we find two periodic solutions as τ 1 is increased from τ 10 .Figures 8(a) and 8(b) show the numerical simulation in time history for the in-phase mode and out-of-phase mode, respectively, at τ 1 = 1.A comparison of the phase portraits between the analytical solutions obtained from Corollary 1 and numerical simulations is depicted in Fig. 8(c).Figures 9(a)-9(c) show the periodic solutions at τ 1 = 1.5.In Fig. 10, the periodic solutions obtained from Corollary 1 are compared with those obtained from numerical simulation as τ 1 is increased.It can be seen that they are in good agreement.

VII. CONCLUSIONS
In this paper, we perform a stability analysis of a pair of van der Pol oscillators with delayed self-connection, position and velocity couplings.By introducing a parameter γ where 0 < γ ≤ 1 in the local stability analysis, stability boundary curves can easily be obtained if they exist.Bifurcation diagram on stability analysis of the damping, position and velocity coupling strengths is constructed, which consists of three types of regions: (i) absolutely unstable region, (ii) regions where closed stability boundary curves exist, and (iii) regions where the stability boundary curves are open-ended.Stability switching regions are observed near the τ 2 -axis.An interesting question arises naturally: for arbitrary large damping strength μ, can the zero equilibrium of a pair of van der Pol oscillators be stabilized using the delayed position and velocity couplings strategy?Theorem 3 gives an affirmative answer.A sufficient condition which requires only one time delay τ 2 is found for the stabilization of the zero equilibrium.With delayed velocity coupling alone, it was shown in Ref. 16 that stabilization of the zero equilibrium can be achieved when μ < 0.5.However, with both delayed position and velocity couplings, it is shown that stability of the zero equilibrium can be achieved for arbitrary large μ.A general question arises naturally, for a network of van der Pol oscillators with any topological connection, is it possible to derive a delayed couplings strategy to stabilize the zero equilibrium for arbitrary large damping strength μ?The question will be dealt with in future.
Periodic solutions of the in-phase and the out-of-phase modes coexist in system (1).We have shown how the type of mode is related to the stability boundary curve that a periodic solution emerges from a Hopf bifurcation.The amplitude of a periodic solution can easily be obtained using an integration method.The analytical solutions agree very well with those from numerical simulation.

ACKNOWLEDGMENT
This work was supported by the Hong Kong Research Grant Council under CERG Grant (CityU 1005/07E).The constructive comments by the anonymous reviewers are gratefully acknowledged.

APPENDIX A
Proof of Theorem 1: Since a positive solution of (9a) is 1 equivalent to a pair of purely imaginary roots of (3a), we only need to prove that the region specified by the inequalities in each part of Lemma 3 or Lemma 4 corresponds to one of the regions defined in this theorem.
Fig. 3(b) shows the curves of γ 2 i = 1 for i ∈ {+, 1, 2} and the regions where γ 2 i is greater or less than one.These curves intersect concurrently at two points .
The region specified by the inequalities of Lemma 3(a) is given by γ 1 > 1 and (γ 2 2 ≥ 1 or γ + > 1) and It can be shown easily that the points satisfying the above inequalities are those to the right of the curves OP + and OP − in Fig. 3(b) (see the shaded area), i.e. the points in region I defined in this theorem (see Fig. 3).From Lemma 3(a), since (9a) has no positive solution for this set of points, (3a) has no purely imaginary root.This proves part (a).Furthermore, the region specified by the inequalities of Lemma 3(b) corresponds to region II.Finally, since the regions specified by the inequalities of Lemma 4(a) and 4(b) correspond to regions III a and III b , respectively.This completes the proof.

FIG. 1 .
FIG.1.A pair of van der Pol oscillators with discrete time delays in the signal transmission of self-connection, position and velocity couplings.

FIG. 8 .
FIG.8.Time history of the oscillators x 1 (solid) and x 2 (dashed) in (a) the out-of-phase mode and (b) the in-phase mode at τ 1 = 1.(c) A comparison of the zero-order out-of-phase (cross) and in-phase (plus) analytical solutions of (29a) with those from numerical simulations in the phase plane of x 1 .

FIG. 9 . 1 )
FIG. 9. Time history of the oscillators x 1 (solid) and x 2 (dashed) in (a) the out-of-phase mode and (b) the in-phase mode at τ 1 = 1.5.(c) A comparison of the zero-order out-of-phase (cross) and in-phase (plus) analytical solutions of (29a) with those from numerical simulations in the phase plane of x 1 .