Stability Switches and Double Hopf Bifurcation Analysis on Two-Degree-of-Freedom Coupled Delay van der Pol Oscillator

: In this paper, the normal form and central manifold theories are used to discuss the inﬂuence of two-degree-of-freedom coupled van der Pol oscillators with time delay feedback. Compared with the single-degree-of-freedom time delay van der Pol oscillator, the system studied in this paper has richer dynamical behavior. The results obtained include: the change of time delay causing the stability switching of the system, and the greater the time delay, the more complicated the stability switching. Near the double Hopf bifurcation point, the system is simpliﬁed by using the normal form and central manifold theories. The system is divided into six regions with different dynamical properties. With the above results, for practical engineering problems, we can perform time delay feedback adjustment to make the system show amplitude death, limit loop, and so on. It is worth noting that because of the existence of unstable limit cycles in the system, the limit cycle cannot be obtained by numerical solution. Therefore, we derive the approximate analytical solution of the system and simulate the time history of the interaction between two frequencies in Region IV.


Introduction
The van der Pol oscillator is a limit cycle oscillation of vacuum tube amplifiers discovered by Dutch scientists. The limit cycle oscillation can be expressed by the following nonlinear differential equations:ẍ The van der Pol oscillator exists in many aspects, such as image encryption [1] and signal detection [2]. It is worth noting that we can only solve the system by approximating analytical and numerical methods, including multiple time scales [3] and average methods [4]. Of course, in addition to the development of research methods, many scholars have studied the dynamic behavior of van der Pol oscillators under the influence of time delay [5,6] and non-smooth oscillators [7].
The coupling of oscillators usually generates many new phenomena, such as synchronization, phase locking, and amplitude death. Stankevich et al. [8] studied quasi-periodic bifurcations of five coupled van der Pol oscillators. Bukh and Anishchenko [9] studied the spiral wave, target wave, and chimeric wave in coupled van der Pol oscillators. Singh and Yadava [10] found transient chaos and stable chaotic dynamics in coupled autonomous van der Pol oscillations, but this is a rare case. Then they revealed that the nonlinear restoring forces in a pair of van der Pol oscillators can induce a transient chaotic route by a small disturbance to the amplitude of an oscillator. Algaba et al. [11] investigated canard explosion in van der Pol electronic oscillators. They developed a new effective program which can calculate the approximate value of critical parameters of any desired x 1 (t) = x 2 (t), The linear part of Equation (2) is as follows: The characteristic equation of the linearized System (3) reduces to When If the characteristic equation has zero real part eigenvalues at the equilibrium point, then the equilibrium point is hyperbolic. The Hartman-Grobman theorem tells us that the orbit topology near the hyperbolic equilibrium point is equivalent to the linearized system.
Eliminating τ 2 from Equation (7), we obtain Where there is fixed τ 2 , if G(ω 01 , τ 2 ) = 0 has one positive root, then there is a critical value τ determined by When the system passes through the critical value, it will lose stability. In order to make sure the Hopf bifurcation occurs, we also need to verify its transversality condition. Without loss of generality, here we calculate the derivative of λ to τ 1 .
Notice that if dλ dτ 1 = 0, then there is a Hopf bifurcation. Thus, we draw the following conclusion.
(iii) When G(ω 01 , τ 2 ) = 0 has several pairs of positive real roots, there exist several critical values. Then we divide the interval, and we can also find that trivial equilibriums of systems are asymptotically stable on finite intervals. When τ 2 = 1, G(ω 01 ) = 0 has two positive real roots, it means there will be an interval which makes the system stable. If τ 1 separates from this interval, the stability of the system will change and the positive real part will appear, which can be observed in Figure 2a. When τ 1 increases, the blue curve obviously crosses the zero line, then the blue curve falls from above the zero line. When τ 2 = 8, the situation becomes more complicated. G(ω 01 ) = 0 has several positive real roots, and the stability switching of the system will be more frequent. It can be seen from Figure 3a that multiple curves of the system move back and forth on the zero line, and it illustrates the complexity of positive real characteristic roots.
In order to explain the stability switching of the system more clearly. We take τ 2 = 1, τ 1 = 0.1, 1 and τ 2 = 8, τ 1 = 1, respectively. Comparing Figure 2b,c with Figure 3b, it is found that time delay can certainly make the system go from stable to unstable.

Analysis of Double Hopf Bifurcation
In this part, we will study the nature of double Hopf bifurcation about the direction and the stability of bifurcating periodic solutions with two time delays, τ 1 and τ 2 . Firstly, we study double Hopf bifurcation based on the normal form theory and center manifold theorem. Secondly, we carry out the numerical simulations.

Computation of Normal Form and Center-Manifold Reduction
For convenience, let β = β 1 = −β 2 ; then we derive the unfolding of the System (2) by the normal form theory and center manifold theorem, and then we are able to reach the stability near the critical value. Now, we need to regulate the time delay by changing the time. Let t → t τ 1 and r = τ 2 τ 1 , where the system becomes: Regarding the time delays τ 1 and r as bifurcation parameters, we set where εδ 1 and εδ 2 are unfolding parameters. It follows that We rewrite the system as follows: To apply the central manifold reduction, we need to rewrite it into functional differential form. C([−τ, 0], R 4 ) is the Banach space of continuous functions, where τ = max{1, r}. For any φ ∈ C, we define: , and the bilinear inner product. (The detailed process is in the Appendix A).
If System (11) has two pairs of purely imaginary eigenvalues Λ = {±iω 01 , ±iω 02 } and other eigenvalues are negative, the phase space C can be divided into two subspaces. That is, C = P Λ ⊕ Q Λ . P Λ is the central subspace obtained by extending the basis vector of linear operator A ε with respect to ±iω 01 , ±iω 02 . Q Λ is its complementary space.
We suppose φ j (θ) and ψ j (s) are the eigenvectors of A(0) and A * (0). They correspond to eigenvalue iω 0j , −iω 0j , j = 1, 2, respectively. By direct computations, we have where p j2 = iω 0j , The real basis of P Λ and its dual space can be expressed as follows: Therefore, it is easy to see thaṫ where Define Z = (z 1 , z 2 , z 3 , z 4 ) T = Ψ, x t , which represents the local coordinates on the central manifold caused by Ψ. We can decompose the phase space by C = P Λ + Q Λ . Then, where ΦZ is the projection of x i on the central manifold. Substituting Equation (21) intȯ x t = L(0)x t + L ε x t + R ε x t and expressing Ψ in its bilinear form, then we can obtain Then combining Equations (21)- (22), we obtain It meansŻ = BZ + Ψ(0)F(t, ΦZ). We know that and Following the computation of the normal forms introduced by [14], we can get the normal form as follows: Furthermore, we can derive the following results: We then apply double polar coordinate transformation by where r 1 , r 2 > 0. Then, Equation (26) becomeṡ where c 1 = Rem 11 εδ 1 + Rem 21 εδ 2 , c 2 = Rem 13 εδ 1 + Rem 23 εδ 2 , a 0 = Rem 2100 , b 0 = Rem 1011 , c 0 = Rem 0021 , d 0 = Rem 1110 . (30)

Classification of Dynamical Behaviours
We take k = −0.5, β = 0.4, l = 0.5, and let τ 1 and τ 2 be the bifurcation parameters. We can draw the curves of Hopf bifurcation when τ 1 , τ 2 vary. As shown clearly in Figure 4, two Hopf bifurcation curves intersect, and we call the intersections the double Hopf bifurcation point, denoted by HH. It follows that the (τ 2 , τ 1 ) plane is divided into different regions, in which the stability of trivial equilibria are different. When the system moves through these curves, the stability switches. Now we take Figure 4b as an example, when τ 1 = 1.693, τ 2 = 0.190, τ 0 01 , and τ 0 02 intersect. Figure 5 shows two pairs of pure virtual roots. Correspondingly, ω can be calculated, and we denote ω 11 = 0.5731, ω 12 = 0.9424. Using the central manifold method, we can get   The detailed process can be referred to in Chapter 7 of the Reference [14]. By classifying the bifurcation solutions, we can clearly analyze the dynamic system near the double Hopf bifurcation point. Letṙ 1 = 0,ṙ 2 = 0 of Equation (29), where we arrived at The value of the equilibrium point depends on c 1 , c 2 in Equation (30). That is when c 1 , c 2 changes near the critical value (τ c 2 , τ c 1 ), and the stability of the equilibrium will change as well. We need to pay attention to the straight lines c 2 = c 1 d 0 b 0 and c 2 = c 0 c 1 a 0 , because there are different dynamic behaviors on both sides of them. At this time, the neighborhood of (τ c 2 , τ c 1 ) is divided into several parts, including region I − V I : We study the normal form of System (29), which is reduced by the central manifold of Equation (11). It can reflect some properties of System (11). Now, Figure 6 shows the area division of the c 1 − c 2 and τ 2 − τ 1 planes at the critical point (τ c 2 , τ c 1 ). In addition, Figure 7 is a phase diagram of different regions. Hence, the phase diagrams of different regions are as follows: When (τ 2 , τ 1 ) is in Region I, Equation (29) shows a stable equilibrium point E 0 = (0, 0), and this area is called the amplitude death region.
When (τ 2 , τ 1 ) is in Region III, the equilibrium point E 0 , E 1 still exists, and a new equilibrium E 2 appears. At this time, E 0 is unstable, E 1 is stable, and E 2 is a saddle.
When (τ 2 , τ 1 ) enters Region IV, it is different from Region III. There appears a new equilibrium point, which is expressed as E * . E 0 , E 1 , and E 2 are unstable and E * is stable in this case.
When (τ 2 , τ 1 ) is in Region V, E * disappears and there are two other equilibria left. E 0 , E 1 remain unstable, and E 2 becomes stable.
When (τ 2 , τ 1 ) enters the last Region VI, the dynamical properties are similar to Region II. The equilibrium E 0 is unstable, but E 2 is stable.

Numerical Simulation
Due to the instability limit cycle generated by the system, the Runge Kutta method could not be calculated for numerical simulation in the case. Therefore, we could deduce the approximate analytical formula of the limit cycle of the System (2) by x = ΦZ. Figure 8 shows the comparison between the approximate analytical formula (red) and numerical solution (blue). It can be seen that when the parameter is closed to the bifurcations of double Hopf, the approximate solution is in good agreement with the numerical solution. It is worth noting that the accuracy of the approximate solution of the system will decrease rapidly when the parameter selection is far away from the double Hopf bifurcation point. When τ 1 = 1.52, τ 2 = 0.25, the approximate analytical solution of x 1 can be expressed as follows: where θ depends on the initial value. Through the phase diagram and time history diagram in Figure 9, it is found that the system is locally asymptotically stable in Region I. With the change of delay τ 1 , the equilibrium point loses its stability because of the supercritical Hopf bifurcation, and a stable limit cycle appeared in Region VI. Note that the system periodically oscillates with ω 01 = 0.5637. When the system is in Region II, the equilibrium point loses its stability again, and it was caused by a subcritical Hopf bifurcation and produced an unstable limit cycle. The system oscillates at ω 01 = 0.5871. However, in Area IV, the system is affected by two frequencies, ω 01 = 0.5695, ω 02 = 0.9602 (Figure 10b). Due to the unstable limit cycle, it is impossible to simulate the system numerically, so an approximate solution of the system is given by

Conclusions
Time delay is an inevitable factor in practice. However, in engineering, more attention is being paid to the phenomenon, and few people study the underlying mechanism theoretically. This paper discussed a type of two-degree-of-freedom van der Pol oscillator with time delay feedback. It was analyzed by using the normal form theory and central manifold theory. Research shows that the time delay transition can cause the stability of the system to switch. Moreover, the greater the time delay, the more complicated the stability switches. Secondly, there are abundant dynamic behaviors near the double Hopf bifurcation. According to the stability of the equilibria, it can be divided into six regions, including the amplitude dead zone and limit cycle zone. It should be noted that, due to the subcritical Hopf bifurcation in the system, it will lead to unstable limit cycles. The Runge Kutta method cannot directly calculate the numerical solution, and the approximate analytical solution of the system can be derived by using the central manifold method.

Data Availability Statement:
The calculation process can be found in [14], and the presentation results and source code of some programs can be found in [27].