Analysis of the Bogdanov-Takens bifurcation in a retarded oscillator with negative damping and double delay

1.   Introduction

A harmonic oscillator is a system that consists of a mass and a spring with a restoring force (proportional to the displacement from the equilibrium position). Harmonic oscillator has been discussed for a long time due to its vast applications, such as in physics and many other fields. If there is a frictional force (damping) proportionate to the velocity, the harmonic oscillator is called a damped oscillator. The amplitude of vibration in damped oscillators reduces over time. Damping is vital in actual oscillatory systems because practically all physical systems include factors like air resistance, friction, and intermolecular interactions, all of which cause energy to be wasted as heat or sound. Positive, zero, or negative damping can be used. Negative damping can be found in a variety of real-world situations. An aeroplane's nose wheel shimmy, for example. A damped harmonic oscillator is any genuine oscillatory system, such as a yo-yo, clock pendulum, or guitar string: after starting to vibrate, the yo-yo, clock, or guitar string eventually slows down and stops, corresponding to the decay of sound or amplitude in general. We also know that positive damping takes the energy out of the system and causes stability, but for negative damping adds energy to the system and causes instability rather oscillation with higher amplitude. When the coefficient of friction curve produce a large negative slope, then the friction induced system possesses negative damping, which can result in self-excited vibration instability. Negative damping also found in the mechanical linkage and power system with great potential of being applied in practical field of applications ^[1]. It is quite obvious in the oscillation caused by the string of violin. Therefore, negative damping is a reality, so we cannot ignore their presence in some nonlinear oscillators.

If we introduce the time delay into the ordinary differential equations (ODEs), we get the delay differential equations (DDEs). Clearly it is more realistic. Following Pyragas ^[2] pioneering work, time-delayed feedback control has been used to a variety of disciplines of inquiry, including chaos communication, optics, electronic systems, biology, and engineering. We know that location feedback cannot affect the amplitude of oscillations without delay ^[3]. As a result, the delay acts as a derivative feedback in modifying the amplitude. The following functional differential equation can be used to explain the damped harmonic oscillator with delayed feedback:

where $x(t)\in{\bf R}$ signifies the distance between the equilibrium position and the current position, $a > 0$ is the spring stiffness constant, $f(x(t-\tau))$ is a $C^r\; (r \geq3)$ function that describes delayed feedback, where $\tau$ is the time delay and $b(<0)$ is the negative damping coefficient. System (1.1) has been analyzed by many authors ^{[4,5,6,7,8,9,10]}. In ^[6] the authors have shown that the system (1.1) with negative damping, steady-state bifurcation, Bogdanov-Takens (B-T) bifurcation, triple-zero, and Hopf-zero bifurcations exist. The stability and steady-state bifurcation problems for the system at the simple zero eigenvalue singularity were solved using center manifold reduction and normal form theory ^[6,7]. Moreover, the second- and third-order normal forms at the origin, as well as the system's stability at the double-zero eigenvalue singularity, have been extracted ^[4,6,11].

So far, there are many studies in the literature on the Hopf/Hopf-zero/triple-zero bifurcations of van der-Pols' oscillators and some unusual neural network models with a single delay or that can be converted to the case with a single delay ^[12,13]. However, there are just a few studies that discuss the B-T bifurcation of harmonic oscillators with multiple delays and negative damping. A very few articles on nonlinear retarded oscillators have been studied by researchers McGahan ^[14], Wang et al ^[15].

Inspired from the discussions in ^{[16,17,18,19,20,21,22]}, when a time delay is inserted into the system (1.1), the result is a retarded damped oscillator with double delays as follows:

We can easily check that under the criteria specified below, the systems' origin (1.2) is a double-zero singularity. The primary goal of this study is to discuss the B-T singularity and extract the systems' corresponding normal form (1.2). The B-T singularity is a zero eigenvalue equilibrium with algebraic multiplicity two and geometric multiplicity one, as we know. We also know that if an ODE has a B-T point, its order is at least two, although this is not the case for DDEs. B-T singularity analysis is a constructive way can provide a plenty of information on the system's dynamics. The normal form computation is a powerful tool for analyzing local bifurcation and stability, and the results for ODEs have been researched for decades. At the B-T singularity, DDEs can be reduced to two-dimensional ODEs using center manifold reduction and normal form theory. Many models, such as predator-prey or neural network models, can describe the B-T bifurcation under certain crucial conditions. The B-T bifurcation of a retarded oscillator with negative damping will be studied in this work.

The following is a summary of the rest of the paper: The criteria under which the origin is a B-T singularity are stated in Section 2. The second and third-order normal forms, as well as the accompanying bifurcation curves at the B-T singularity are precisely described in Section 3. Section 4 presents various numerical simulations to demonstrate the obtained results. Conclusions in Section 5 bring the paper to an end.

2.   Analysis of B-T singularity

Retarded damped harmonic oscillator with two delays can be taken as follows:

Let $x_1(t) = x(t)$ and $x_2(t) = \dot{x}_1(t)$. Then the equation (2.1) is broken up to the following system:

Throughout this paper, it is assumed that $a > 0, \; b < 0;\; f(0) = 0,$ and $f'(0) = d.$ Henceforth, the characteristic equation of the linearized part of (2.2) at the origin is given by

Lemma 1. Let $d = a, \; b = -a\tau_2$ and $2+b(\tau_2-2\tau_1) > 0.$ Then $\lambda = 0$ is a double-zero root of Eq (2.3).

Proof. Clearly, $F(0) = a-d$. Since $d = a$, we have $F(0) = 0$. Also, we have $F'(0) = 2\lambda+be^{-\lambda\tau_1}-b\tau_1\lambda e^{-\lambda \tau_1}+d\tau_2e^{-\lambda\tau_2}.$ Then $F^{\prime}(0) = 0$ as $b = -a\tau_2$. Again, $F{''}(0) = 2+b(\tau_2-2\tau_1).$ Since, $2+b(\tau_2-2\tau_1) > 0$, we have $F^{\prime\prime}(0)\neq0$. Henceforth, we can conclude that $0$ is a double-zero root of the characteristic Eq (2.3).

Lemma 2. If $d = a$, $b = -a\tau_2$, $2+b(\tau_2-2\tau_1) > 0$, and $a\in (0, a_0^{+})$, then all roots of (2.3) except $\lambda = 0$, have non-zero real parts, i.e. origin of the system (2.2) is B-T singularity, where $a_0^{+} = \{\min a_j:\; a_j = \frac{w_j^2}{1-\tau_2w_j\sin{\tau_1w_j}-\cos{\tau_2w_j}} > 0, \; 1\leq j\leq m, \; and\; w_j\; are\; the\; roots\; of \; the\; equation$ $w^2-(b^2+2a)+\frac{2ab\sin(\tau_1-\tau_2)w}{w} = 0\}.$

Proof. Let us consider $iw$ be a root of the the Eq (2.3) if $-w^2+ibwe^{-i\tau_1 w}+a-ae^{-i\tau_2 w} = 0,$ $\text{i.e. if}\; -w^2+a+bw(\sin \tau_1w+i\cos\tau_1w)-a(\cos \tau_2w-i\sin\tau_2 w) = 0$, i.e., if

Given that $d = a, \; b = -a\tau_2$. Hence, under these conditions $-iw$ is also a root of the Eq (2.3).

Thus, we can assume $w > 0$.

Now (2.4) gives $w^4-(b^2+2a)w^2+2abw\sin(\tau_1-\tau_2)w = 0,$ i.e.,

We set $G(w) = w^2-(b^2+2a)+2ab\frac{\sin(\tau_1-\tau_2)w}{w}.$ We know that $\lim_{x\rightarrow 0}g(x) = t, \; g(x) = \frac{\sin tx}{x}.$ We can compute $G(0) = -a[2+b(\tau_2-2\tau_1)] < 0.$ Since, $G(w)\rightarrow \infty$ as $w\rightarrow \infty,$ there exists $w^* > 0$ such that $G(w^*) = 0.$ Here we observe that for large value of $w$, $G(w)\approx w^2$. Hence (2.6) have finite number of positive roots. Let the roots are $w_1, w_2, \cdots, w_m\; (m\geq1).$

Also we see that $w\leq A$, where $A = \sqrt{a[2-b(\tau_2+2|\tau_1-\tau_2|)]}.$

Now, we will differentiate the Eq (2.3) with respect to $a$. Then, we have

Thus, $\frac{d\lambda}{da}|_{\lambda = iw} = \frac{-w^2+ibw(\cos\tau_1w-i\sin\tau_1w)}{a(\Theta+i\chi)}$, where $\Theta = b(\cos\tau_1w-\tau_1w\sin\tau_1w-\cos\tau_2w)$ and $\chi = 2w-b(\sin\tau_1w+\tau_1w\cos\tau_1w-\sin\tau_2w).$

Consequently, $\Re\left(\frac{d\lambda}{da}\right)|_{\lambda = iw} = \frac{b[-b\tau_1\omega^2 +\omega^2(\cos \tau_1\omega+\cos\tau_2\omega)-b\omega\sin(\tau_1-\tau_2)w+\tau_1w^3\sin\tau_1w]}{a(\Theta^2+\chi^2)}\ne 0.$

3.   Bogdanov-Takens bifurcation analysis

By Lemma 2, it is clear that, if $d = a\in(0, a_0^{+})$ and $b = -a\tau_2$, then the system (2.2) at the origin experiences a B-T bifurcation. Thus, we can take d and b as bifurcation parameters. Hence, we can introduce two small parameters $\mu_1$ and $\mu_2$, which changing in a small neighborhood $V$ of $(0, 0)^T$, then discuss the effect of perturbation on the system (2.2)

where h.o.t. stands for higher order terms of $x_1(t-\tau_2).$ By simplifying, we can write the system (3.1) as the following retarded functional differential equation on the phase space $C$:

where $C$ is the Banach space of all continuous functions from $\phi:\, [-\tau_1, 0] \rightarrow {\bf R}^2$ with supremum norm $|\phi| = \sup_{[-\tau_1, 0]}|\phi(\theta)|,$ $X_t\in C$, defined by $X_t(\theta) = X(t+\theta)$, $\forall$ $\theta\in[-\tau_1, 0],$ $\mu = (\mu_1, \mu_2)^T\in V,$ and $L(\mu):C\rightarrow {\bf R}^2$ is a set of parameterized bounded linear operators defined as follows:

and $G:C\times V\rightarrow {\bf R}^2$ is a $C^m\; (m\geq2)$ function satisfying $G(0, 0), \; DG(0, 0) = 0$ with

Now, we have the linearization of the system (3.2) at the point $(X_t, \mu) = (0, 0)$ as follows:

The operator $L_0 = L(0)$ can be written as :

where

In the same way, we can represent the linear operator $L_1$ also.

Here $L_0$ be a bounded linear operator. Then by the Riesz representation theorem, we have a matrix $\eta(\cdot)$ of order two, defined on $[-\tau_1, 0]$ with components of bounded variation such that

From ^[23] and ^[24], we know very well that the fundamental solutions of system the (3.5) form a $C_0$-semigroup $\{T_0(t):t\geq0\}$ on $C$ with infinitesimal generator $A_0: C\rightarrow C$, defined as $A_0(\phi) = \dot{\phi}$ on the domain

Let us consider the adjoint space of $C$ as $C^* = C([0, \tau_1], {\bf R}^{2*})$, where ${\bf R}^{2*}$ be the 2-dimensional row vector space. The adjoint inner product on $C^*\times C$ is given by

Let $\Lambda_0$ be the set of eigenvalues of $A_0$ with zero real parts and counting multiplicity. Clearly, for the B-T bifurcation $\Lambda_0 = \{0, 0\}$. Let $P$ be the invariant space of $A_0$ associated with the zero real part eigenvalues and $P^*$ be the corresponding dual space. So, the dimension of $P$ is two. Also, the dimension of $P^*$ is two. Now, we will use the formal adjoint theory for a functional differential equation and express the phase space $C$ as $C = P\oplus Q$, where $Q = \{\phi\in C:\left < \psi, \phi\right > = 0\; \forall\psi\in P^*\}$. Also, it is know that $Q$ is invariant under both of $T_0(t)$ and $A_0$.

Let $\Phi$ and $\Psi$ be the respective bases of $P$ and $P^*$. We can chose $\Phi$ and $\Psi$ as follows: $\Phi = (\phi_1(\theta), \phi_2(\theta))$ for $-\tau_1\leq\theta\leq 0$, and $\Psi = (\psi_1(s), \psi_2(s))^T$ for $0\leq s\leq\tau_1.$ Then $\left < \Psi(s), \Phi(\theta)\right > = I_2$, and $\dot{\Phi} = \Phi\bar {B},$ $-\dot{\Psi} = \bar{B}\Psi,$ where $\bar{B} = \left(\begin{array}{cc} 0 & 1 \\ 0 & 0 \\ \end{array} \right).$ After applying the same method as in Lemma $3.1$ in ^[25], we get $\Phi(\theta) = \left(\begin{array}{cc} 1 & \theta \\ 0 & 1 \\ \end{array} \right), \; -\tau_1\leq\theta\leq 0$ and

$\Psi(0) = \left(\begin{array}{cc} \frac{2b^2(\tau_2^2-3\tau_1^2)}{3[2+b(\tau_2-2\tau_1)]^2}+\frac{2(1-b\tau_1)}{2+b(\tau_2-2\tau_1)} & \frac{2b(\tau_2^2-3\tau_1^2)}{3[2+b(\tau_2-2\tau_1)]^2} \\ \frac{2b}{2+b(\tau_2-2\tau_1)} & \frac{2}{2+b(\tau_2-2\tau_1)} \\ \end{array} \right) = \left(\begin{array}{cc} \psi_{11} & \psi_{12} \\ \psi_{21} & \psi_{22} \\ \end{array} \right).$

For discussion of B-T bifurcation of system (3.2), we enlarge the phase space $C$ by Banach space $BC = \{\phi:[-\tau_1, 0]\rightarrow {\bf R}^2:\phi\; \text{is continuous on}\; [-\tau_1, 0), \; \text{with aprobable}\; \text{jump discontinuity near zero}\}.$ Any element $\phi$ of $BC$ can be taken as $\phi = \varphi+X_0 c$. In $BC$, the norm is given by $||\phi|| = ||\varphi+X_0c|| = ||\varphi||_C+|c|$, where $X_0$ is a square matrix valued function of order $2$, satisfying

Thus, in $BC$, the system (3.2) can be rewrite as the following abstract ordinary differential equations:

where $F(u, \mu) = (L(\mu)-L_0)u+G(u, \mu) = L_1(\mu)u+G(u, \mu), \; \mu\in V$ and $\bar{A}$ is the extension of infinitesimal generator of $A_0$, defined by $\bar{A}:C^1\rightarrow BC$, satisfying the following equation:

The continuous projection operator $\pi:BC\rightarrow P$ is defined as follows:

Hence, by $A_0$, the phase space $BC$ can be decomposed as $BC = P\oplus \ker\pi$. As $\pi$ commutes with $\bar{A}$ in $C^1$. We can take $u = {\bf \phi}(\theta)z(t)+y$ and then from the result of ^[26], we may decomposed the abstract ODE (3.9) into the following system

where $\Psi(0) = \left < \Psi, X_0\right > , \; z = (z_1, z_2)^T\in{\bf R}^2 = P, \; y = (y_1, y_2)^T\in Q^1 = Q\cap C^1\subset \ker\pi,$ $A_{Q^1}$ is the restriction of $\bar{A}$ as an operator from $Q^1$ to the Banach space $\ker\pi$. Assume $f^1(z, y, \mu) = \Psi(0)F(\Phi z+y, \mu)$ and $f^2(z, y, \mu) = (I-\pi)X_0F(\Phi z+y, \mu)$. Expanding $F(\Phi z+y, \mu)$ at $(z, y, \mu) = (0, 0, 0),$ (3.12) according to Taylor series expression, we have

where $f^k_{j}(z, y, \mu)\; , k = 1, 2$ denote the homogeneous polynomials of $z, y, \mu$ of degree $j$. With the help of the Eq (3.1) and (3.13), we have the following relations:

$\frac{1}{2}f^1_2(z, y, \mu) = \Psi(0)[\frac{1}{2!}F_2(\Phi z+y, \mu)] = \Psi(0)\left(\begin{array}{c} \bar{F}^1_2(z, y, \mu) \\ \bar{F}^2_2(z, y, \mu) \\ \end{array} \right),$

$\frac{1}{2}f^2_2(z, y, \mu) = (I-\pi)X_0[\frac{1}{2!}F_2(\Phi z+y, \mu)] = (I-\pi)X_0\left(\begin{array}{c} \bar{F}^1_2(z, y, \mu) \\ \bar{F}^2_2(z, y, \mu) \\ \end{array} \right),$

$\frac{1}{3!}f^1_3(z, y, \mu) = \Psi(0)[\frac{1}{3!}F_3(\Phi z+y, \mu)] = \Psi(0)\left(\begin{array}{c} 0\\ \bar{F}^2_3(z, y, \mu) \\ \end{array} \right),$

$\frac{1}{3!}f^2_3(z, y, \mu) = (I-\pi)X_0[\frac{1}{3!}F_3(\Phi z+y, \mu)] = (I-\pi)X_0\left(\begin{array}{c} 0\\ \bar{F}^2_3(z, y, \mu) \\ \end{array} \right),$ where

Let $V_j^4({\bf R}^2)$ be the vector space of homogeneous polynomials of $(z, \mu) = (z_1, z_2;\mu_1, \mu_2)$ of degree $j$ having coefficients in ${\bf R}^2$. Then

where $(q, l) = (q_1, q_2;l_1, l_2)\in {\bf N}^4, \; z^q\mu^l = z_1^{q_1}z_2^{q_2}\mu_1^{l_1}\mu_2^{l_2}, \; q_1+q_2+l_1+l_2 = j.$ We can take the canonical basis for $V_2^4({\bf R}^2)$ as:

$\Bigg\{\left(\begin{array}{c} z_i^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} z_1z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i z_1 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_1\mu_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_i^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i^2 \\ \end{array} \right),$ $\left(\begin{array}{c} 0 \\ z_1z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_1 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_1\mu_2 \\ \end{array} \right); \, i = 1, 2\Bigg\}.$

Canonical basis for $V_3^4({\bf R}^2)$ can be taken as:

$\Bigg\{\left(\begin{array}{c} z_i^3 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i^3 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} z_1^2z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_iz_1^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i^2z_1 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} z_1z_2^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_1\mu_2 z_i \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i z_1z_2 \\ 0 \\ \end{array} \right),$ $\left(\begin{array}{c} \mu_i z_2^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i^2 z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_1^2\mu_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_1\mu_2^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_i^3 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i^3 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_1^2z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_1^2 \\ \end{array} \right),$ $\left(\begin{array}{c} 0 \\ \mu_i^2z_1 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_1z_2^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_1\mu_2 z_i \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_1z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_2^2\mu_i \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_2\mu_i^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_1^2\mu_2 \\ \end{array} \right),$ $\left(\begin{array}{c} 0 \\ \mu_1\mu_2^2 \\ \end{array} \right); \, i = 1, 2\Bigg\}.$

Define the operator $M_j^1$ on $V_j^4({\bf R}^2)$ by

From ^[25] and ^[27], $V_2^4({\bf R}^2)$ can be decomposed as follows:

Since $M_2^1\left(\begin{array}{c} z_1^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 2z_1z_2 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} 0 \\ z_1^2 \\ \end{array} \right) = \left(\begin{array}{c} -z_1^2 \\ 2z_1z_2 \\ \end{array} \right), M_2^1\left(\begin{array}{c} z_1z_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} z_2^2 \\ 0 \\ \end{array} \right),$ $M_2^1\left(\begin{array}{c} 0 \\ z_1z_2 \\ \end{array} \right) = \left(\begin{array}{c} -z_1z_2 \\ z_2^2 \\ \end{array} \right), M_2^1\left(\begin{array}{c} z_2^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} 0 \\ z_2^2 \\ \end{array} \right) = \left(\begin{array}{c} -z_2^2 \\ 0 \\ \end{array} \right),$ $M_2^1\left(\begin{array}{c} \mu_iz_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} \mu_iz_2 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} 0 \\ \mu_iz_1 \\ \end{array} \right) = \left(\begin{array}{c} -\mu_iz_1 \\ \mu_iz_2 \\ \end{array} \right), M_2^1\left(\begin{array}{c} \mu_iz_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right),$ $M_2^1\left(\begin{array}{c} 0 \\ \mu_iz_2 \\ \end{array} \right) = \left(\begin{array}{c} -\mu_iz_2 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} \mu_i^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} 0 \\ \mu_i^2 \\ \end{array} \right) = \left(\begin{array}{c} -\mu_i^2 \\ 0 \\ \end{array} \right),$ $M_2^1\left(\begin{array}{c} \mu_1\mu_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right), M_2^1\left(\begin{array}{c} 0 \\ \mu_1\mu_2 \\ \end{array} \right) = \left(\begin{array}{c} -\mu_1\mu_2 \\ 0 \\ \end{array} \right)$ for $i = 1, 2$; we have $I_m(M_2^1) = \Bigg < \Bigg\{\left(\begin{array}{c} -z_1^2 \\ 2z_1z_2 \\ \end{array} \right), \left(\begin{array}{c} z_1z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} z_2^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_2^2 \\ \end{array} \right), \left(\begin{array}{c} \mu_i z_1 \\ -\mu_i z_2 \\ \end{array} \right), \left(\begin{array}{c} \mu_i z_2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} \mu_i^2 \\ 0 \\ \end{array} \right),$ $\left(\begin{array}{c} \mu_1\mu_2 \\ 0 \\ \end{array} \right):i = 1, 2\Bigg\}\Bigg > .$

As $\left(\begin{array}{c} \mu_iz_1 \\ -\mu_iz_2 \\ \end{array} \right)+\left(\begin{array}{c} 0 \\ \mu_iz_2 \\ \end{array} \right) = \left(\begin{array}{c} \mu_iz_1 \\ 0 \\ \end{array} \right)$ for $i = 1, 2$; we have

$I_m(M_2^1)^c = \left < \left\{\left(\begin{array}{c} 0 \\ z_1^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_1z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_1 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_1\mu_2 \\ \end{array} \right):i = 1, 2\right\}\right > .$

Let $P_{I, j}^1$ be the mapping from $V_j^4({\bf R}^2)$ to $I_m(M_j^1)$ for $j = 2, 3$ which satisfies

For $j = 2$, we can see that

Also $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} z_1^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 2z_1z_2 \\ \end{array} \right)$, and $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} \mu_iz_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ \mu_iz_2 \\ \end{array} \right)\; \text{for}\; i = 1, 2.$

By ^[26,28], on the centre manifold connected to the space $P$, the normal form of (2.2) can be written as

If $f_2^1(z, 0, \mu) = \left(\begin{array}{c} a_1z_1^2+a_2z_1z_2+a_3z_2^2+a_4\mu_1z_1+a_5\mu_2z_1+a_6\mu_1z_2+a_7\mu_2z_2 \\ b_1z_1^2+b_2z_1z_2+b_3z_2^2+b_4\mu_1z_1+b_5\mu_2z_1+b_6\mu_1z_2+b_7\mu_2z_2 \\ \end{array} \right),$ then

$Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_1z_1^2 \\ 0 \\ \end{array} \right) =$ $\left(\begin{array}{c} 0 \\ 2a_1z_1z_2 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_2z_1z_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_3z_2^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_4\mu_1z_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ a_4\mu_1z_2 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_5\mu_2z_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ a_5\mu_2z_2 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_6\mu_1z_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} a_7\mu_2z_2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_1z_1^2 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ b_1z_1^2 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_2z_1z_2 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ b_2z_1z_2 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_3z_2^2 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_4\mu_1z_1 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ b_4\mu_1z_1 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_5\mu_2z_1 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ b_5\mu_2z_1 \\ \end{array} \right),$ $Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_6\mu_1z_2 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ b_6\mu_1z_2 \\ \end{array} \right), Pr_{I_m(M_2^1)^c}\left(\begin{array}{c} 0 \\ b_7\mu_2z_2 \\ \end{array} \right) =$ $\left(\begin{array}{c} 0 \\ b_7\mu_2z_2 \\ \end{array} \right).$

Hence, $g_2^1(z, 0, \mu)$ in the normal (3.17) is given by $g_2^1(z, 0, \mu) = Pr_{I_m(M_2^1)^c}f_2^1(z, 0, \mu)$$= \left(\begin{array}{c} 0 \\ 2\psi_{22}\mu_1z_1+2[(\psi_{12}-\tau_2\psi_{22})\mu_1-\psi_{22}\mu_2]z_2+f{''}(0)\psi_{22}z_1^2+2(\psi_{12}-\tau_2\psi_{22})f{''}(0)z_1z_2 \\ \end{array} \right).$

Now, on the center manifold, the system (2.2) can be transformed into the following normal form:

where $\delta_1 = \psi_{22}\mu_1, \; \delta_2 = (\psi_{12}-\tau_2\psi_{22})\mu_1-\psi_{22}\mu_2, a_{2} = \frac{f^{\prime\prime}(0)}{2}\psi_{22}, \; b_{2} = (\psi_{12}-\tau_2\psi_{22})f^{\prime\prime}(0).$

Since $a\in(0, a_0^+)$, we may get a small $a$ such that $\psi_{12}-\tau_2\psi_{22} < 0.$ In addition, it is obvious that $\psi_{22} > 0$ as we have assumed that $2+b(\tau_2-2\tau_1) > 0$. As a result, discovering the sign of $f^{\prime\prime}(0)$ can yield the signs of the coefficients $a_2$ and $b_2$.

Now we will discuss different cases depending on the sign of $a_2$ and $b_2$.

Case Ⅰ: If $f^{\prime\prime}(0) > 0,$ then $a_2 > 0$ and $b_2 < 0.$ Re-scaling the time parameter and transform the coordinates in the following way:

Then on the center manifold, the system (3.18) becomes (after dropping bars)

where $\nu_1 = \frac{4(\psi_{12}-\tau_2\psi_{22})^2}{\psi_{22}}\mu_1$ and $\nu_2 = -\frac{2(\psi_{12}-\tau_2\psi_{22})^2}{\psi_{22}}\mu_1+2(\psi_{12}-\tau_2\psi_{22})\mu_2.$

The bifurcation curves associated to the perturbation parameters $\mu_1, \; \mu_2$ are sum up as follows ^[4,7,9,27]:

$1)$ Transcritical bifurcation (TB) occurs when $\mu_1 = 0.$

$2)$ The system (3.19) experiencing Hopf bifurcation around the zero equilibrium point when $\mu_2 = -\frac{6\tau_2+b(3\tau_1^2+2\tau_2^2-6\tau_1\tau_2)}{3[2+b(\tau_2-2\tau_1)]}\mu_1$ and $\mu_1 < 0$.

$3)$ The system (3.19) experiencing Hopf bifurcation around the axial equilibrium point when $\mu_2 = \frac{6\tau_2+b(3\tau_1^2+2\tau_2^2-6\tau_1\tau_2)}{3[2+b(\tau_2-2\tau_1)]}\mu_1$ and $\mu_1 > 0$.

$4)$ The system (3.19) experiencing Homoclinic bifurcation at the trivial equilibrium when $\mu_2 = \frac{5[6\tau_2+b(3\tau_1^2+2\tau_2^2-6\tau_1\tau_2)]}{21[2+b(\tau_2-2\tau_1)]}\mu_1$ and $\mu_1 > 0$.

$5)$ The system (3.19) experiencing Hopf bifurcation around the axial equilibrium point when $\mu_2 = -\frac{5[6\tau_2+b(3\tau_1^2+2\tau_2^2-6\tau_1\tau_2)]}{21[2+b(\tau_2-2\tau_1)]}\mu_1$ and $\mu_1 > 0$.

Case Ⅱ: If $f^{\prime\prime}(0) < 0,$ then $a_2 < 0$ and $b_2 > 0.$ Now we re-scale the time and transform the coordinates in the following way:

Then on the center manifold, the system (3.18) becomes (after dropping bars)

where $\nu_1, \; \nu_2$, and the associated bifurcation curves are the same as those of Case Ⅰ.

Case Ⅲ: Here we will discuss the case for $f^{\prime\prime}(0) = 0.$ Then, clearly $a_2 = b_2 = 0,$ and hence system becomes degenerate. To discuss the dynamics near the B-T singularity we have to compute the higher order normal forms. From ^[26] and ^[27], we may get

where $\bar{f}_3^1(z, 0, \mu) = f_3^1(z, 0, \mu)+\frac{3}{2}[(D_zf_2^1)(z, 0, \mu)U_2^1(z, \mu)$$-(D_zU_2^1)(z, \mu)g_2^1(z, 0, \mu)+ (D_yf_2^1)(z, 0, \mu)U_2^2(z, \mu)].$ Now, we can easily obtain

To get $g_3^1(z, 0, \mu)$, we have to calculate $U_2(z, \mu) = (U_2^1(z, \mu), U_2^2(z, \mu))^T.$ From ^[4] and ^[27], one can get

where $U_2^1\in\ker(M_2^1)^c$ and $\ker(M_2^1)^c$ is spanned by

$\left\{\left(\begin{array}{c} z_1^2 \\ 0 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_i^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ z_1z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_1 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i z_2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_i^2 \\ \end{array} \right), \left(\begin{array}{c} 0 \\ \mu_1\mu_2 \\ \end{array} \right):i = 1, 2 \right\}$

Thus, $U_2^1(z, \mu) = 2\tau_2\psi_{12}\left(\begin{array}{c} 0 \\ \mu_1z_2 \\ \end{array} \right)+\psi_{12}\left(\begin{array}{c} 0 \\ \mu_2z_2 \\ \end{array} \right)-2\psi_{12}\left(\begin{array}{c} 0 \\ \mu_1z_1 \\ \end{array} \right).$

Hence, $(D_zf_2^1)(z, 0, \mu)U_2^1(z, \mu) = \left(\begin{array}{c} \psi_{12}\upsilon \\ \psi_{22}\upsilon \\ \end{array} \right),$

where $\upsilon = 4\psi_{12}[\tau_2\mu_1^2(z_1-\tau_2z_2)+\mu_1\mu_2(z_1-2\tau_2z_2)-\mu_2^2z_2].$

Now, $(D_zU_2^1)g_2^1(z, 0, \mu) = \left(\begin{array}{c} 0 \\ 4\psi_{12}(\tau_2\mu_1+\mu_2)[\psi_{22}\mu_1z_1+(\psi_{12}-\tau_2\psi_{22})\mu_1z_2-\psi_{22}\mu_2z_2] \\ \end{array} \right).$

On the other side, $U_2^2(z, \mu) = h_2(\theta)(z, \mu)\in V_2^4(Q^1)$, where $h_2(\theta) = (h_2^1(\theta), h_2^2(\theta))^T$ which satisfies $(M_2^2U_2^2)(z, \mu) = f_2^2(z, 0, \mu).$ Using the formula of $\bar{A}$, one can get

Consequently,

The expression of $h_2(\theta)(z, \mu) = \left(\begin{array}{c} h_2^1(z, \mu) \\ h_2^2(z, \mu) \\ \end{array} \right)$ with degree 2 can be evaluated as $h_2^i(\theta)(z, \mu) = \sum_{|(q, l)| = 2}h_{2(q, l)}^i(\theta)z^q\mu^l = h_{22000}^i(\theta)z_1^2+h_{20200}^i(\theta)z_2^2+h_{20020}^i(\theta)\mu_1^2+h_{20002}^i(\theta)\mu_2^2 +$$h_{21100}^i(\theta)z_1z_2+h_{21010}^i(\theta)\mu_1z_1+h_{21002}^i(\theta)\mu_2z_1+h_{20110}^i(\theta)\mu_1z_2+h_{20101}^i(\theta)\mu_2z_2+h_{20011}^i(\theta)\mu_1\mu_2$ for $i = 1, 2.$

Comparing the coefficients of $\mu_1z_1, \; \mu_1z_2, \; \mu_2z_2$ in (3.24) and (3.25), we have the following conditions:

From the above relations, we have

$h_{21010} = \left(\begin{array}{c} \psi_{22}\theta^2+2\psi_{12}\theta+d_1 \\ 2\psi_{22}\theta+2\psi_{12} \\ \end{array} \right),$

$h_{20110}(\theta) = \left(\begin{array}{c} \frac{1}{3}\psi_{22}\theta^3+(\psi_{12}-\tau_2\psi_{22})\theta^2+(d_1-2\tau_2\psi_{12})\theta+d_2\\ \psi_{22}\theta^2+2(\psi_{12}-\tau_2\psi_{22})\theta+d_1-2\tau_2\psi_{12} \\ \end{array} \right),$

$h_{20101}(\theta) = \left(\begin{array}{c} -\psi_{22}\theta^2-2\psi_{12}\theta+d_3 \\ -2\psi_{22}\theta-2\psi_{12} \\ \end{array} \right);$

and the other elements of $h_{2(q, l)}$ are all zero. Furthermore, $h_{2(q, l)}(\theta)\in Q^1 = Q\cap C^1$ and satisfies

Then $d_j$ can be determined from (3.29), as follows:

$d_1 = \frac{1}{2+b\tau_2}[\frac{\psi_{22}(4\tau_1^3-\tau_2^3)b}{6}+\frac{4\psi_{12}^2}{\psi_{22}^2}],$

$d_2 = \frac{2b}{2+b\tau_2}[\frac{\psi_{22}(6\tau_2^4-5\tau_1^4)}{60}+\frac{(\psi_{12}-\tau_2\psi_{22})\tau_1^3-2\psi_{12}\tau_2^3-(d_1-3\tau_2\psi_{12})\tau_1^2}{3}],$

$d_3 = -d_1.$

Thus, we get

To get the third order normal form, we need the relationship

By ^[27], we have the basis for $I_m(M_3^1)$ and $I_m(M_3^1)^c$. We know that

and for the other bases in $V_3^4({\bf R}^2)$, we have

$Pr_{I_m(M_3^1)^c}\left(\begin{array}{c} z_1^3 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 3z_1^2z_2 \\ \end{array} \right), \; Pr_{I_m(M_3^1)^c}\left(\begin{array}{c} \mu_i^2z_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ \mu_i^2z_2 \\ \end{array} \right),$

$Pr_{I_m(M_3^1)^c}\left(\begin{array}{c} \mu_1\mu_2z_1 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ \mu_1\mu_2z_2 \\ \end{array} \right), \; Pr_{I_m(M_3^1)^c}\left(\begin{array}{c} \mu_i z_1^2 \\ 0 \\ \end{array} \right) = \left(\begin{array}{c} 0 \\ 2\mu_iz_1z_2 \\ \end{array} \right), i = 1, 2.$

Now, from the above, we can see that

From Eqs (3.31)–(3.34), we have the normal form of the system (2.2) as

where $\lambda_1 = \psi_{22}\mu_1+\frac{1}{2}\psi_{22}h_{21010}^1(-\tau_2)\mu_1^2-\frac{1}{2}\psi_{22}h_{21010}^2(-\tau_1)\mu_1\mu_2$, $\lambda_2 = (\psi_{12}-\tau_2\psi_{22})\mu_1-\psi_{22}\mu_2+\frac{1}{2}[\psi_{22}h_{20110}^1(-\tau_2)+\psi_{12}h_{21010}^1(-\tau_2)]\mu_1^2$$+\frac{1}{2}[\psi_{22}h_{20101}^1(-\tau_2)-\psi_{22}h_{20110}^2(-\tau_1)-$$\psi_{12}h_{21010}^2(-\tau_1)] \mu_1\mu_2-3\psi_{22}h_{20101}^2(-\tau_1)\mu_2^2,$ $a_3 = \frac{1}{6}f^{'''}(0)\psi_{22}, \; b_3 = \frac{1}{2}f^{'''}(0)(\psi_{12}-\tau_2\psi_{22}).$

Now, we will use the following time re-scaling and co-ordinate transformation :

Then (3.35) becomes

where $\sigma_1 = (\frac{b_3}{a_3})^2\lambda_1, \; \sigma_2 = -\frac{b_3}{|a_3|}\lambda_2, \; s = sgn(a_3).$

From ^[29,30], we know that the bifurcations of the system (3.37) are linked with the sign of $s$. If $s = 1$, the bifurcation curves associated to the perturbation parameters $\mu_1, \; \mu_2$ are sum up as follows ^[4,9,27]:

$({\rm{a}})$ The system (3.37) attains a pitchfork bifurcation on the parametric curve

$({\rm{b}})$ The system (3.37) attains a Hopf bifurcation at the trivial equilibrium on the parametric curve

${\rm{(c}})$ The system (3.37) attains a Heteroclinic bifurcation at the trivial equilibrium on the parametric curve

If $s = -1$, the bifurcation curves associated to the perturbation parameters $\mu_1, \; \mu_2$ are sum up as follows ^[5,6,9]:

$({\rm{i}})$ The system (3.37) attains a pitchfork bifurcation on the parametric curve

$({\rm{ii}})$ The system (3.37) attains a Hopf bifurcation at the trivial equilibrium on the parametric curve

$({\rm{iii}})$ The system (3.37) attains a Hopf bifurcation at the non-trivial equilibrium on the parametric curve

$({\rm{iv}})$ The system (3.37) attains a Homoclinic bifurcation on the parametric curve

$({\rm{v}})$ The system (3.37) attains a fold bifurcation of the limit cycle on the parametric curve

4.   Numerical simulations

Here, numerical simulations of the delayed system (2.2) are performed to illustrate the results obtained above. To confirm our analytical and theoretical finding, we cite some numerical simulations with the help of computing softwares MATLAB-R2015a, Maple-18 and Mathematica7.0. Firstly, let $f(x) = \sin{x},$ we have $f (0) = 0$, $f^{\prime}(0) = 1$. To simulate the dynamics of the delayed system (2.2) near the B-T bifurcation, we assume $a = f^{'}(0) = 0.0222, \; \tau_1 = 20, \; \tau_2 = 1.$ By Lemma 2, the critical B-T bifurcation point is $(d, b) = (a, -a) = (0.0222, -0.0222).$ Now, we consider a small perturbation of the bifurcation parameters by letting $(d, b) = (0.0222+\mu_1, 0.0222+\mu_2).$ The following cases are discussed to test the bifurcation.

If the perturbation parameter $(d, b) = (-0.0122, 0.0322)$ corresponding to the perturbation $(\mu_1, \mu_2) = (0.01, 0.01)$, the trivial equilibrium of the system (2.2) is a saddle point and other two non-trivial equilibria are stable focus (cf. Figure 1 (a)). For the same but opposite perturbation results the merging of all three equilibria into one stable trivial equilibrium (cf. Figure 1 (b)). For the slight smaller perturbation $(\mu_1, \mu_2) = (0.001, 0.001)$ and slightly higher feedback delay $\tau_2 = 2$ makes the domain into four sub-domains by two sepatrices: separatrix curve shown by incoming thick arrows and the separatrix curve, shown by outgoing thick arrows (boundaries of saddle point $(0, 0)$. The separatrix $x_2\approx -x_1$ divide the domain into two sub-domains: the left one is the basin of attraction of $E_{-}(- 0.5118910020, 0)$ and right one is the basin of attraction of $E_{+}(+0.5118910020, 0)$ (cf. Figure 2). Next we find the dynamics near the B-T point $(-a\tau_2, a) = (-0.0444, 0.0222)$ for the set $a = 0.0222, \, b = -0.0222, \, d = 0.0222, \, \tau_1 = 20$, with the perturbation $(\mu_1, \mu_2) = (-0.001, 0.001)$, a saddle point origin and two stable focus placed symmetrically left and right of $(0.0)$ (cf. Figure 3 (a)) for the feedback delay parameter $\tau_2 = 2$, two small stable limit cycles around the axial equilibria, surrounded by a bigger stable limit cycle for$(\mu_1, \mu_2) = (0.02, 0.02)$ (cf. Figure 3 (b)). This higher amplitude limit cycles surrounding the stable equilibrium points is an important phenomenon from a practical point of view. From the Figure 4, it is observed that without delay $(\tau_1 = \tau_2 = 0)$, the system (2.2) will become unbounded, while it was confined within a bounded domain under both the effects of retarded delay $(\tau_1 = 20)$ and feedback delay $(\tau_2 = 2)$.

If $f(x) = \sin{x},$ and the fixed set $a = 0.0222, \, b = -0.0222, \, d = 0.0222, \, \tau_2 = 1$, with the perturbation $(\mu_1, \mu_2) = (-0.1, -0.1)$, the system becomes asymptotically stable around origin (cf. Figure 5) for the retarded delay parameter $\tau_1 = 10$, an unstable/semi-stable has been emerged from the origin for $\tau_1 = 15$ (cf. Figure 6), a chaotic attractor is emerged for $\tau_1 = 20$ (cf. Figure 7) and system becomes unbounded for all $\tau_1\geq 25$ (cf. Figure 8). For same set of parameters and without any perturbation around B-T critical value the system generates a double-homoclinic loop around origin for no feedback delay $\tau_2 = 0$ (cf. Figure 9), but it becomes an attractor around origin for $\tau_2 = 2$ and $(\mu_1, \mu_2) = (-0.02, \, 0.1)$ (cf. Figure 10). There exists an attractor covering two unstable axially symmetric equilibria and the saddle origin for for $f(x) = \sin{4x}$ and $d = 0.0222\times f'(0) = 0.0888$ (cf. Figure 11).

Next, we consider $f(x) = \tanh(x+1)-\tanh{(1)},$ then $f (0) = 0$, $f^{'}(0) = 1-\left(\frac{e^2-1}{e^2 +1}\right)^2 > 0$ and $f^{\prime\prime}(0) = -2\left(\frac{e^2-1}{e^2 +1}\right)\left(1-\left(\frac{e^2-1}{e^2 +1}\right)^2\right) < 0$. In this case, we find complicated dynamics near the B-T point. For the set $a = 0.0222, \, b = -0.0222, \, d = 0.0222, \, \tau_1 = 20$, with the perturbation $(\mu_1, \mu_2) = (-0.13, 0.02)$, a double-well chaotic attractor emerges from the origin (cf. Figure 12) for the feedback delay parameter $\tau_2 = 1$, a triple-well attractor is emerged from the origin for $\tau_2 = 2$ (cf. Figure 13) and a chaotic attractor is emerged for $\tau_2 = 3$ (cf. Figure 14). It is also observed that the system (2.2) becomes regular after a small transient period around origin for $\tau_2 = 7$, lastly it becomes unbounded for all $\tau_2 > 7$ (Figures are not reported here). Bifurcation diagram in Figure 15 shows that the aperiodic oscillation can be suppressed by increasing the coefficient stiffness of the spring. Is is seen from the Figure 15 that the the aperiodic solution converges to equilibrium state when the value of the parameter $a$ becomes larger than $a = 0.0222\in (0, a^{+}) = (0, 1.541191394)$. Hence, the feedback delay plays a crucial role in regulating the amplitude on the oscillations of the system (2.2) as we observed for the retarded delays $\tau_1 = 10, \, 15, \, 20, \, 25$ in Figures 5–8.

5.   Conclusions

We have studied a retarded oscillator with negative damping and two delays. It is found that the origin of the delayed system (2.2) is a B-T bifurcation point if $d = a\in (0, a_0^{+})$ and $b = -a\tau_2$. Utilizing the center manifold and normal form theories, we derived the canonical forms of B-T singularity. Moreover, the phase portraits and bifurcation diagrams along with associated criteria for several cases of the normal form have shown. Finally, numerical simulations are given which confirmed the obtained criteria and observed that both the retarded delay $\tau_1$ and feedback delay $\tau_2$ have crucial functioning on the qualitative change of the dynamics of the delayed system (2.2). From the numerical study, it is evident that the retarded delay plays a significant role in regulating the dynamics of the system under consideration. It is clear from the Figures 12–14 that the feedback delay is also an important factor in controlling the dynamics of the proposed system. It can be concluded that the the coefficient of stiffness has a crucial role in regulating the irregular dynamics of the system (2.2) (cf. Figure 15).

We know very well that a system with multiple delays makes it very hard to investigate the distribution of the eigenvalues of the characteristic equation and induce the system to show richer dynamical behaviors. There are hardly any universal unfolding results about the triple zero bifurcation for a system with retarded delay as well as feedback delay. As a result, research on the quadruple bifurcation for delayed systems is scarce. These results have guiding importance for the engineers to choose the values of the delay to acquire the desired dynamical effects. These facts will be discussed in the future studies.

Acknowledgements

The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the financial support for this research under the number 10367-cos-2020-1-3-I during the academic year 1442 AH / 2020 AD.

Conflict of interest

The authors declare no conflicts of interest.