%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% %% This is the mc_template.tex file for the mc document class. %% %% It is used to prepare s manuscript for Mathematical %% %% Communications journal. %% %% %% %% The mc.cls class works only with a pdflatex engine. %% %% The file newmc.cls should be placed where LaTeX %% %% can find it, e.g. in the current working directory. %% %% %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentclass{mc} %%===============================================================%% %% Please add here your own packages, macros and enviroments. %% %% It is not necessary to include ams* and graphicx packages %% %% since they are automatically included by the mc class. %% %% Avoid defining your own environments and use the already %% %% defined ones (e.g.~theorem, lemma etc.) %% %%===============================================================%% %\usepackage{enumerate} % uncomment to use this package %\newcommand{\E}{\mathbb{E}} % example of a macro %%===============================================================%% %%===============================================================%% %% Journal info will be edited by the typesetter %% %% DO NOT CHANGE THIS PART %% %%===============================================================%% \setcounter{page}{283} \renewcommand\thisnumber{x} \renewcommand\thisyear {2025} \renewcommand\thismonth{xxx} \renewcommand\thisvolume{30} \renewcommand\datereceived{October 30, 2024} \renewcommand\dateaccepted{July 9, 2025} \renewcommand\doinum{10.1000/100} %%===============================================================%% \begin{document} % \linenumbers %%===============================================================%% %% TITLE %% %% Please add the title with \title[Short title]{Title} %% %% Short title is a running head apearing in the header. %% %%===============================================================%% \title[On the dynamics of discrete-time Rayleigh-Duffing oscillators]{On the dynamics of discrete-time Rayleigh-Duffing oscillators} %%===============================================================%% %%===============================================================%% %% AUTHOR(S) %% %% %% %% Add author's details in the following format. For each author %% %% provide the affiliation, address and Orcid identifier. %% %% Mark the corresponding author with \comma\corrauth. %% %%===============================================================%% \author[M.~Fatehi Nia, Z.~Nurkanovi\' c, and N.~Khajoei] {Mehdi Fatehi Nia\affil1\orcidnumber{0000-0002-8965-9359}, Zehra Nurkanovi\' c\affil2\comma\corrauth\orcidnumber{0000-0002-0917-8433}, and Najmeh Khajoei\affil3\orcidnumber{0000-0002-1252-4210}} \address{\affilnum1 Department of Mathematics, Yazd University, P.\,O.\,Box 89195-741, Yazd, Iran\\ \affilnum2 Department of Mathematics, Faculty of Natural Sciences and Mathematics, University of Tuzla, Univerzitetska 4, 75\,000 Tuzla, Bosnia and Herzegovina\\ \affilnum3 Department of Mathematics, Shahid Bahonar University of Kerman, 22 Bahman Boulevard, Kerman 76169-14111, Iran } \emails{% \email{fatehiniam@yazd.ac.ir} (M.~Fatehi Nia), \email{zehra.nurkanovic@untz.ba} (Z.~Nurkanovi\' c), \email{khajuee.najmeh@yahoo.com} (N.~Khajoei) } %%===============================================================%% %% ABSTRACT %% %%===============================================================%% \begin{abstract} This paper investigates the dynamics of a discrete-time Rayleigh-Duffing oscillator exhibiting chaos in the Li-Yorke sense. We use Marotto's theorem to prove the existence of chaos by finding a snap-back repeller. It is shown that the system undergoes Neimark-Sacker and a period-doubling bifurcations. We compute the Lyapunov exponents numerically to show sensitive dependence on initial conditions and chaotic behavior. The results are illustrated using numerical simulations. \end{abstract} %%===============================================================%% %%===============================================================%% %% KEYWORDS %% %%===============================================================%% \keywords{period-doubling bifurcation; Neimark-Sacker bifurcation; Lyapunov exponent; Marotto's method; Li-Yorke chaos} %%===============================================================%% %%===============================================================%% %% AMS subject classification %% %%===============================================================%% \ams{39A10, 39A28, 39A30} %%===============================================================%% %%===============================================================%% \maketitle %%===============================================================%% %%===============================================================%% %% MAIN BODY %% %%===============================================================%% \section{Introduction} \label{Sec:1} Discretization of ordinary differential equations is an approach to approximate and analyze their dynamical behavior. Some discrete schemes are e.g. the trapezoidal rule method, midpoint, forward, and backward Euler's methods \cite{zbMATH06445212}. For an autonomous differential equation \begin{equation} \label{Eqn1} y^{\prime}=f(y), \end{equation} where $y$, $f\in\mathbb{R}^{n}$ the flow map is defined by the function $% \phi_{\tau}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, where $\phi_{\tau }(y_{0})=y(\tau)$ and $y(t)$ is the solution of the associated differential equation \eqref{Eqn1} with the initial condition $y(0)=y_{0}$ for $t\in \lbrack0,\tau]$. The flow map of a differential equation defines a discrete dynamical system (difference equation) when the vector field $f$ is globally Lipschitz, since there is a unique solution to any initial condition $% y_{0}\in\mathbb{R}^{n}$ for all time. Among the discretization methods of ordinary differential equations, forward Euler's method approximates the flow map. \newline This paper considers the Rayleigh-Duffing oscillator, which is demonstrated as $\ddot{x}+\beta\dot{x}+\dot{x}^{3}+\alpha x+x^{3}=0$. Duffing's function $% \alpha x+x^{3}$ shows a nonlinear stiffness in this system, and Rayleigh's function $\beta\dot{x}+\dot{x}^{3}$ shows nonlinear damping. Also, $\alpha$ and $\beta$ describe the linear damping coefficient and linear stiffness, respectively. This system can be shown equivalently in the form of the planar system \begin{equation} \label{n1} \begin{array}{rl} \dot{x}=y,\quad \dot{y}=-\alpha x-\beta y-x^{3}-y^{3}, & \end{array}% \end{equation} where $\alpha$, $\beta\in\mathbb{R}$. This system is a mathematical model in many physical and mechanical problems; for instance, see \cite{bikdash1994melnikov,herfort1987robust,zbMATH03651791}.% It behaves similarly to the Van der Pol-Duffing and Van der Pol oscillator, see \cite{YAMAPI2006187}. In \cite{kanai2012creation}, the authors studied qualitative analyses with small positive linear and cubic damping. In \cite{freire1998study}, Freire et al. studied the bifurcation diagram when $|\alpha|$ and $|\beta|$ are small. In \cite{chen2016global}% , Chen and Zou studied the global dynamics of Rayleigh-Duffing oscillators. Also, the dynamic behavior of the general form of Rayleigh-Duffing oscillators is investigated by Chen et al. in \cite{chen2018saddle}. By choosing various values for parameters $\alpha$ and $\beta$, system (\ref{n1}) has different behaviors. So, the forward Euler method applied to system \eqref{n1} is an approach to investigate the approximate dynamical behavior of its flow map. Hence, by applying the forward Euler scheme to system \eqref{n1}, we obtain the following discrete-time model: \begin{equation} \label{n2} \begin{array}{rl} x_{n+1}=x_{n}+hy_{n},\quad y_{n+1}=y_{n}-h(\alpha x_{n}+\beta y_{n}+x_{n}^{3}+y_{n}^{3}), & \end{array}% \end{equation} where the step size is $h>0$. The main goal of this paper is to study chaotic behaviors and bifurcation conditions of the given discrete system. This paper is organized as follows: The following section will briefly discuss the existence of fixed points and describe their local stability conditions. In Section \ref{Sec:3}, sufficient conditions for the existence of codimension-one bifurcations, including Neimark-Sacker and period-doubling bifurcations, are investigated. In Section \ref{Sec:4}, it is shown that system \eqref{n2} is chaotic in the sense of Li-Yorke using Marotto's method. Finally, we used Maple and Mathematica to conduct numerical simulations, confirming our theoretical results. \section{Local stability of fixed points} \label{Sec:2} In this section, we determine the existence of the fixed points of the discrete-time system and analyze their stability. It is clear that system % \eqref{n2} has $E_{0}=(0,0)$ as a fixed point for all $\alpha$, $\beta\in% \mathbb{R}$. Furthermore, it has the origin and $E_{\pm }=(\pm\sqrt{-\alpha}% ,0)$ as fixed points for $\alpha<0$. Let $T$ denote the map associated with system \eqref{n2}, i.e. \begin{equation*} T\left( \begin{array}{c} x \\ y% \end{array} \right) =\left( \begin{array}{c} x+hy \\ y-h\left( \alpha x+\beta y+x^{3}+y^{3}\right)% \end{array} %\right), \label{n2a} \right), \label{n2a} \end{equation*} the Jacobian matrix of $T$ at the point $\left( x,y \right) $ is \begin{equation*} J_{T}\left( x,y \right) =\left[ \begin{array}{cc} 1 & h \\ -\alpha h-3hx^{2} & 1-\beta h-3hy^{2}% \end{array}% \right] , %\right] , \label{JT} \end{equation*}% and characteristic polynomial \begin{equation} F(\lambda)=\lambda ^{2}+\left( 3hy^{2}+h\beta -2\right) \lambda +3h^{2}% x^{2}+h^{2}\alpha -3hy^{2}-h\beta +1. \label{CPT} \end{equation}% To study the stability of $E_{0}=(0,0)$ and $E_{\pm }=(\pm\sqrt{-\alpha}% ,0)$, we use \cite[Lemma 2.2]{LIU200780}. In the following proposition, the stability of system \eqref{n2} at the fixed point $% E_{0}=(0,0)$ will be considered and used in the rest of the paper to find some conditions that eventuate bifurcations at the origin. \begin{proposition} \label{n4} Assume that $\alpha$, $\beta\in\mathbb{R}$, and $h>0$. \newline Then the equilibrium point $E_{0}=(0,0)$ of system \eqref{n2} is: \begin{itemize} \item[$(a)$] a sink for $\max \left\{ 0,\frac{2\beta h-4}{h^{2}}\right\} <\alpha <\frac{\beta }{h}$ and $0<\beta <\frac{4}{h}$, \item[$(b)$] a saddle for $\frac{2\beta h-4}{h^{2}}<\alpha<0$ and $\beta<% \frac {2}{h}$ or $0<\alpha<\frac{2\beta h-4}{h^{2}}$ and $\beta>\frac{2}{h}$, \item[$(c)$] a source for $\beta\in\mathbb{R}$ and $\alpha<\min\left\{ 0,\frac{2\beta h-4}{h^{2}}\right\} $ or $% \alpha>\max\left\{ 0,\frac{2\beta h-4}{h^{2}},\frac{\beta}{h}\right\},$ \item[$(d)$] a non-hyperbolic with \begin{itemize} \item[$1)$] $\lambda_{1,2}$ are the conjugate roots and $\left\vert \lambda_{1,2}\right\vert =1$ for $\alpha=\frac{\beta}{h}\neq0$ and $0<\beta<% \frac{4}{h}$, \item[$2)$] $\lambda _{1}=-1$, $\left\vert \lambda _{2}\right\vert \neq 1$ or $% \left\vert \lambda _{1}\right\vert \neq 1$, $\lambda _{2}=-1$ for $\alpha =% \frac{2\beta h-4}{h^{2}}$, $\beta \neq \frac{2}{h}$, $% \beta \neq \frac{4}{h}$ and $\beta \neq 0$, \item[$3)$] $\lambda _{1}=1$, $-1\neq \lambda _{2}<1$ for $\alpha =0$ and $% \frac{2}{h}\neq \beta >0$, \item[$4)$] $\lambda _{1}=1$, $\lambda _{2}>1$ for $\alpha =0,\beta <0$, \item[$5)$] $\lambda _{1}=\lambda _{2}=1$\ for $M_{1}=\left( \beta ,\alpha \right) =\left( 0,0\right) $, \item[$6)$] $\lambda _{1}=-1$, $\lambda _{2}=1$\ for $M_{2}=\left( \beta ,\alpha \right) =\left( \frac{2}{h},0\right) $, \item[$7)$] $\lambda _{1}=-1$, $\lambda _{2}=-1$\ for $M_{3}=\left( \beta ,\alpha \right) =\left( \frac{4}{h},\frac{4}{h^{2}}\right) $. \end{itemize} \end{itemize} \end{proposition} See Figure \ref{P1_0}. \newpage \noindent \begin{figure}[h!] \centering \includegraphics[width=11.5cm]{fig1} \caption{The stability of $E_{0}=(0,0)$ in the $(\protect\beta,\protect\alpha)$% -plane, $M_{1}=\left( 0,0\right)$, $M_{2}=\left( \frac{2}{h},0\right)$ and $% M_{3}=\left( \frac{4}{h},\frac{4}{h^{2}}\right).$} \label{P1_0} \end{figure} \begin{proof} For $h>0$ and $\alpha $, $\beta \in %TCIMACRO{\U{211d} }% %BeginExpansion \mathbb{R} %EndExpansion $, we have \begin{equation*} J_{T}\left( E_{0}\right) =J_{T}\left( 0,0\right) =\left[ \begin{array}{cc} 1 & \quad h \\ -\alpha h & \quad 1-\beta h% \end{array}% \right] , \end{equation*}% and the corresponding characteristic equation is \begin{equation} F\left( \lambda \right) =\lambda ^{2}+\left( \beta h-2\right) \lambda +\alpha h^{2}-\beta h+1=\lambda ^{2}+B\lambda +C, \label{FE_0} \end{equation}% from which \begin{equation} \lambda _{1}=\frac{2-\beta h-h\sqrt{\beta ^{2}-4\alpha }}{2},\quad \lambda _{2}=\frac{2-\beta h+h\sqrt{\beta ^{2}-4\alpha }}{2}. \label{sv_vr} \end{equation}% Since \begin{equation*} F\left( 1\right) =\alpha h^{2}\left\{ \begin{array}{c} <0\text{ for }\alpha <0, \\ =0\text{ for }\alpha =0, \\ >0\text{ for }\alpha >0,% \end{array}% \right. \end{equation*}% we have three different cases. \begin{itemize} \item[i)] %\end{description} For $\alpha =0$: \begin{equation*} \lambda _{1}=\dfrac{2-\beta h-h\left\vert \beta \right\vert }{2}\left\{ \begin{array}{c} =1\text{ for }\beta \leqslant 0, \\ <1\text{ for }\beta >0,% \end{array}% \right. \end{equation*} and \begin{equation*} \lambda _{2}=\dfrac{2-\beta h+h\left\vert \beta \right\vert }{2}\left\{ \begin{array}{c} >1\text{ for }\beta <0, \\ =1\text{ for }\beta \geqslant 0.% \end{array}% \right. \end{equation*} \end{itemize} \begin{itemize} \item[ii)] %\end{description} For $\alpha <0$: \begin{equation*} \lambda _{1}<1\Leftrightarrow -\beta h\frac{2\beta h-4}{h^{2}}. \end{align*}% It is clear that $\lambda _{2}>1$ for $\alpha <0$ and $\beta \in %TCIMACRO{\U{211d} }% %BeginExpansion \mathbb{R} %EndExpansion $. \end{itemize} \begin{itemize} \item[iii)] If $\alpha >0$, we have that $B^{2}-4C=\beta ^{2}-4\alpha \geq 0$ for $0<\alpha \leq \frac{\beta ^{2}}{4}$ and \begin{align*} B&=\beta h-2\neq \left\{ \begin{array}{l} 0\text{ for }\beta \neq \frac{2}{h},\smallskip \\ 2\text{ for }\beta \neq \frac{4}{h},% \end{array}% \right. \\[1ex] C&=\left( \alpha h-\beta \right) h+1\left\{ \begin{array}{l} =1\text{ for }\alpha =\frac{\beta }{h},\smallskip \\ >1\text{ for }\alpha >\frac{\beta }{h},\smallskip \\ <1\text{ for }\alpha <\frac{\beta }{h},% \end{array}% \right. \\[1ex] C&=1\text{ and }B^{2}-4C=\beta ^{2}-4\alpha <0\text{ for }\alpha =\frac{\beta }{h}\text{ and }0<\beta <\frac{4}{h}, \\[1ex] F\left( -1\right) &=\alpha h^{2}-2\beta h+4\left\{ \begin{array}{c} >0\text{ for }\alpha >\frac{2\beta h-4}{h^{2}},\smallskip \\ =0\text{ for }\alpha =\frac{2\beta h-4}{h^{2}},\smallskip \\ <0\text{ for }\alpha <\frac{2\beta h-4}{h^{2}},% \end{array}% \right. \end{align*} \end{itemize} and by using Lemma 2.2 \cite{LIU200780}, there follows the conclusion. \end{proof} \begin{proposition}\label{np2} Assume $\alpha<0$, $\beta\in% %TCIMACRO{\U{211d} }% %BeginExpansion \mathbb{R} %EndExpansion $, and $h>0$. Then equilibrium points $E_{\pm}=\left( \pm\sqrt{-\alpha },0\right) $ of system \eqref{n2} are: \begin{itemize} \item[$(a)$] a source for $\alpha<0$, $\beta\leqslant0$ and $\alpha<-\frac{% \beta }{2h}$ and $\alpha<\frac{2-h\beta}{h^{2}}$, \item[$(b)$] a saddle for $\alpha<0$ and $\alpha>\frac{2-h\beta}{h^{2}}$, $% \beta>\frac{2}{h}$, \item[$(c)$] a non-hyperbolic with \begin{itemize} \item[1)] $\lambda_{1}=-1$, $\left\vert \lambda_{2}\right\vert \neq1$ for $% \alpha<0$ and $\alpha=\frac{2-h\beta}{h^{2}}$, $\beta>\frac{2}{h}$, $% \beta\neq\frac{4}{h}$, \item[2)] $\lambda_{1}=-1$, $\lambda_{2}=1$ for $M_{4}=\left( \beta ,\alpha\right) =\left( \frac{4}{h},-\frac{2}{h^{2}}\right) $, \item[3)] $\lambda_{1,2}$ are the conjugate roots and $\left\vert \lambda_{1,2}\right\vert =1$ for $\alpha=-\frac{\beta}{2h}$ and $0<\beta<% \frac{4}{h}$. \end{itemize} \end{itemize} \end{proposition} See Figure \ref{P2_E}. \newpage \noindent \begin{figure}[h!] \centering \includegraphics[width=11.5cm]{fig2} \caption{The stability of $E_{\pm}$ in the $(\protect\beta,\protect\alpha )$% -plane, $\protect\alpha<0$, $M_{1}=\left( 0,0\right)$, $M_{2}=\left( \frac{2% }{h},0\right)$ and $M_{4}=\left( \frac{4}{h},-\frac{2}{h^{2}}\right)$} \label{P2_E} \end{figure} \begin{proof} For $\alpha<0$, $\beta\in% %TCIMACRO{\U{211d} }% %BeginExpansion \mathbb{R} %EndExpansion $ and $h>0$, we have $J_{T}\left( E_{\pm}\right) =\left[ \begin{array}{cc} 1 & \quad h \\ 2\alpha h & \quad1-\beta h% \end{array} \right]$, and% \begin{equation*} \label{FE_1_2} F\left( \lambda\right) =\lambda^{2}+\left( \beta h-2\right) \lambda -2\alpha h^{2}-\beta h+1=\lambda^{2}+B \lambda+C. \end{equation*} For $\alpha<0$ and $\beta\leq0$, equilibriums $E_{\pm}$ are the source because \begin{equation*} C=-2\alpha h^{2}-\beta h+1>1, \end{equation*} \begin{equation*} F\left( 1\right) =-2h^{2}\alpha>0,\quad F\left( -1\right) =2\left( 2-h\beta-h^{2}\alpha\right) >0. \end{equation*} For $\beta>0$, we have \begin{align*} B&=\beta h-2\left\{ \begin{array}{c} \neq0\text{ for }\beta\neq\frac{2}{h},\smallskip \\ \neq2\text{ for }\beta\neq\frac{4}{h},% \end{array} \right. \\[1ex] C&=1-h\left( \beta+2h\alpha\right) \left\{ \begin{array}{c} <1\text{ for }\alpha>-\frac{\beta}{2h},\smallskip \\ =1\text{ for }\alpha=-\frac{\beta}{2h},\smallskip \\ >1\text{ for }\alpha<-\frac{\beta}{2h},% \end{array} \right. \\[1ex] C&=1\text{ and }B^{2}-4C=h\beta\left( h\beta-4\right) <0\text{ for }0<\beta<% \frac{4}{h},\alpha=-\frac{\beta}{2h}, \\[1ex] F\left( 1\right) &=-2h^{2}\alpha>0, \\[1ex] F\left( -1\right) &=2\left( 2-h\beta-h^{2}\alpha\right) \left\{ \begin{array}{l} <0\text{ for }\alpha>\frac{2-h\beta}{h^{2}},\beta>\frac{2}{h},\smallskip \\ =0\text{ for }\alpha=\frac{2-h\beta}{h^{2}},\beta>\frac{2}{h},\smallskip \\ >0\text{ for }\alpha<\frac{2-h\beta}{h^{2}}.% \end{array} \right. \end{align*} Using \cite[Lemma 2.2]{LIU200780}, the above proves all the claims. \end{proof} Also, for the local stability of all equilibrium points, see Figure \ref{P3}. \begin{figure}[h!] \centering \includegraphics[width=12.0cm]{fig3} \caption{The stability of $E_{0}$ and $E_{\pm}$ in the $(\protect\beta,\protect% \alpha)$-plane, $M_{1}=\left( 0,0\right)$, $M_{2}=\left( \frac{2}{h}% ,0\right) $, $M_{4}=\left( \frac{4}{h},-\frac{2}{h^{2}}\right)$} \label{P3} \end{figure} \newpage \noindent \begin{remark} Propositions \ref{n4} and \ref{np2} give us a clear idea of the local stability of the equilibrium points of system (\ref{n2}). From this, it can be concluded that system (\ref{n2}) dynamics are very complicated. Namely, in addition to Neimark-Sacker and period-doubling bifurcations, the equilibrium $E_0$ also has codimension-two bifurcations (1:1, 1:2, possibly 1:3 and 1:4, as well as the bifurcation that occurs in (d) 6) of Proposition \ref{n4}). On the other hand, $E_\pm$ equilibria are unstable, except possibly in cases (c) 2), and 3) of Proposition \ref{np2}). Examining all these cases requires much more time and would take up much more space in the manuscript so that it can be the subject of future research. It is very important in situations where a codimension-two bifurcation occurs and when the shape of the trajectory of the orbits that move between all equilibrium points needs to be precisely determined. For this reason, we decided only to examine Neimark-Sacker and period-doubling bifurcations of the equilibrium point $E_0$ and theoretical proof of the appearance of chaos in these situations. \end{remark} \section{Bifurcation analysis} \label{Sec:3} In this section, we study the conditions for the existence of Neimark-Sacker and period-doubling bifurcations of system \eqref{n2} at the fixed point $E_{0}=(0,0)$, see \cite{garic2017stability,zbMATH07937099,zbMATH07026661,zbMATH06603714,zbMATH07047385,kulenovic2002discrete,zbMATH07391783}. We choose the step size $h$ as a bifurcation parameter to study these bifurcations. \subsection{ Neimark-Sacker bifurcation} Based on statement $(d)$, $1)$ of Proposition \ref{n4}, system \eqref{n2} has a pair of complex multipliers on the unit circle. We use the standard version of the Naimark-Sacker result, see \cite[Theorems 15 and 16]{zbMATH07026661} and \cite{hale2012dynamics,zbMATH07047385}. \begin{theorem} Assume that {$h>0,h_{0}:=\frac{\beta }{\alpha }<\frac{4}{\beta }$, }$a\left( h_{0}\right) =\frac{3\beta ^{3}\left( \alpha ^{2}-\beta \right) }{8\alpha ^{4}}$, {$\alpha >0$, $\beta >0$}. %\begin{description} \begin{itemize} \item[$(i)$] For $\alpha \in \left( 0,\sqrt[3]{4}\right) $ and $\beta \in \left( 0,\alpha ^{2}\right) $ or $\alpha \geq \sqrt[3]{4}$ and $\beta \in \left( 0,2\sqrt{\alpha }\right) $, i.e., $a\left( h_{0}\right) >0$, there is a neighborhood $\mathcal{U}$ of the origin and a $\delta >0$ such that for $|h-h_{0}|<\delta $ and $\left( x_{0},y_{0}\right) \in \mathcal{U}$% , then the $\omega $-limit set of $\left( x_{0},y_{0}\right) $ is the origin when $hh_{0}$. Furthermore, $\Gamma (h_{0})=0.$ \item[$(ii)$] For $\alpha \in \left( 0,\sqrt[3]{4}\right) $ and $\beta \in \left( \alpha ^{2},2\sqrt{\alpha }\right) $, i.e., $a\left( h_{0}\right) <0$% , there is a neighborhood $\mathcal{U}$ of the origin and a $\delta >0$ such that for $|h-h_{0}|<\delta $ and $\left( x_{0},y_{0}\right) \in \mathcal{U}$, then the $\alpha $-limit set of $\left( x_{0},y_{0}\right) $ is the origin when $h>h_{0}$ and belongs to a closed invariant $C^{1}$ curve $\Gamma (\lambda )$ encircling the origin if $hh_{0}$) and \ref{f1NH} ($h=h_{0}$) for visual illustration. We rewrite system \eqref{n2} in the following form: \begin{equation} \begin{bmatrix} x_{n+1} \\ y_{n+1}% \end{bmatrix}% =% \begin{bmatrix} 1\quad & h \\ -h\alpha \quad & 1-h\beta \end{bmatrix}% \begin{bmatrix} x_{n} \\ y_{n}% \end{bmatrix}% +% \begin{bmatrix} 0 \\ -h(x_{n}^{3}+y_{n}^{3})% \end{bmatrix}% . \label{n6} \end{equation}% The Jacobian $J_{T}\left( 0,0\right) $ has two non-real eigenvalues: $\lambda \left( h\right) =\frac{2-h\beta +ih\sqrt{4\alpha -\beta ^{2}}}{% 2}$ and $\overline{\lambda \left( h\right) }$, for $4\alpha -\beta ^{2}>0$, which implies that $|\lambda \left( h\right) |_{h=h_{0}}=1$. Since% \begin{equation*} \frac{d|\lambda \left( h\right) |}{dh}=\frac{1}{2}\left( 2\alpha h-\beta \right) \left( 1-h\beta +\alpha h^{2}\right) ^{-\frac{1}{2}}, \end{equation*}% we obtain $d\left( h_{0}\right) :=\left. \frac{d|\lambda \left( h\right) |}{% dh}\right\vert _{h=h_{0}}=\frac{\beta }{2}>0$. Also, we have that \begin{align*} \left. \lambda ^{k}\left( h\right) \right\vert _{h=h_{0}}& =\left. exp\left( ik\arctan \left( \tfrac{h\sqrt{4\alpha -\beta ^{2}}}{2-h\beta }\right) \right) \right\vert _{h=h_{0}} \\ & =exp\left( ik\arctan \left( \tfrac{\beta \sqrt{4\alpha -\beta ^{2}}}{% 2\alpha -\beta ^{2}}\right) \right) , \end{align*}% and \begin{equation*} \left. \lambda ^{k}\left( h\right) \right\vert _{h=h_{0}}\neq 1\Leftrightarrow \arctan \left( \tfrac{\beta \sqrt{4\alpha -\beta ^{2}}}{% 2\alpha -\beta ^{2}}\right) \neq 0, \quad k=1,2,3,4, \end{equation*}% (since $4\alpha -\beta ^{2}>0$) which implies that the third condition of Poincare-Andronov-Hopf bifurcation for maps (see \cite[p.~474, Theorem 15.31]{hale2012dynamics}) is valid. Now, our goal is to obtain the normal form of system \eqref{n6}, when $h=h_{0}$. First, denote $\lambda \left( h_{0}\right) =k+ie$ (with $k=\frac{2-h_{0}\beta }{2} $, $e=\frac{h_{0}\sqrt{4\alpha -\beta ^{2}}}{2}$), and let \begin{equation} \label{sistem2} \begin{bmatrix} x_{n} \\ y_{n}% \end{bmatrix}% =P% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix},% \end{equation} where $\tilde{x}_{n}$, $\tilde{y}_{n}$ are new variables, and% \begin{equation*} P=\left[ \begin{array}{cc} h_{0}\quad & 0 \\ k-1\quad & -e% \end{array}% \right], \quad P^{-1}=\left[ \begin{array}{cc} \frac{1}{h_{0}}\quad & 0 \\ \frac{k-1}{h_{0}e}\quad & -\frac{1}{e}% \end{array}% \right] , \end{equation*} such that $P^{-1}J_{T}\left( 0,0\right) P=\left[ \begin{array}{cc} k & -e \\ e & k% \end{array}% \right] $, with $\left\vert \begin{array}{cc} k & -e \\ e & k% \end{array}% \right\vert =1$. Let us notice that from (\ref{sistem2}) we have that \begin{equation*} x_{n}=h_{0}\tilde{x}_{n},\quad y_{n}=\left( k-1\right) \tilde{x}_{n}-e\tilde{% y}_{n}. \end{equation*}% Then, using (\ref{sistem2}) and $h=h_{0}$, system \eqref{n6} has the form% \begin{equation*} P% \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =% \begin{bmatrix} 1\quad & h_{0} \\ -h_{0}\alpha \quad & 1-h_{0}\beta \end{bmatrix}% P% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +% \begin{bmatrix} \phi (\tilde{x}_{n},\tilde{y}_{n}) \\ \psi (\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% , \end{equation*}% where \begin{eqnarray*} \phi (\tilde{x}_{n},\tilde{y}_{n}) &=&0, \\ \psi (\tilde{x}_{n},\tilde{y}_{n}) &=&-h_{0}(\left( h_{0}\tilde{x}% _{n}\right) ^{3}+\left( \left( k-1\right) \tilde{x}_{n}-e\tilde{y}% _{n}\right) ^{3}), \end{eqnarray*}% from which we obtain its normal Jordan form as% \begin{equation*} \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =P^{-1}% \begin{bmatrix} 1\quad & h_{0} \\ -h_{0}\alpha \quad & 1-h_{0}\beta \end{bmatrix}% P% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +P^{-1}% \begin{bmatrix} \phi (\tilde{x}_{n},\tilde{y}_{n}) \\ \psi (\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% , \end{equation*}% i.e., \begin{equation*} \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =% \begin{bmatrix} k & -e \\ e & k% \end{bmatrix}% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +% \begin{bmatrix} f(\tilde{x}_{n},\tilde{y}_{n}) \\ g(\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% , \end{equation*}% where% \begin{equation*} \begin{bmatrix} f(\tilde{x}_{n},\tilde{y}_{n}) \\ g(\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% =P^{-1}% \begin{bmatrix} \phi (\tilde{x}_{n},\tilde{y}_{n}) \\ \psi (\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% , \end{equation*}% and \begin{align*} f(\tilde{x},\tilde{y}) &=0 \\ g(\tilde{x},\tilde{y}) &=h_{0}\dfrac{\left( \left( k-1\right) ^{3}+h_{0}^{3}\right) \tilde{x}^{3}-3e\left( k-1\right) ^{2}\tilde{x}^{2}% \tilde{y}+3e^{2}\left( k-1\right) \tilde{x}\tilde{y}^{2}-e^{3}\tilde{y}^{3}}{% e}. \end{align*}% Now, for $(\tilde{x},\tilde{y})=\left( 0,0\right) $, we get \begin{align*} g_{\tilde{x}\tilde{x}}&=g_{\tilde{x}\tilde{y}}=g_{\tilde{y}\tilde{y}}=0,\text{ \ \ }g_{\tilde{x}\tilde{y}\tilde{y}}=6h_{0}e\left( k-1\right) ,\text{ \ \ }% g_{\tilde{y}\tilde{y}\tilde{y}}=-6h_{0}e^{2}, \\[1ex] g_{\tilde{x}\tilde{x}\tilde{x}}&=\frac{6h_{0}\left( h_{0}^{3}+\left( k-1\right) ^{3}\right) }{e}, \\[1ex] \xi _{20}&=\frac{1}{8}\left\{ (f_{\tilde{x}\tilde{x}}-f_{\tilde{y}\tilde{y}% }+2g_{\tilde{x}\tilde{y}})+i(g_{\tilde{x}\tilde{x}}-g_{\tilde{y}\tilde{y}% }-2f_{\tilde{x}\tilde{y}})\right\} =0, \\[1ex] \xi _{11}&=\frac{1}{4}\left\{ (f_{\tilde{x}\tilde{x}}+f_{\tilde{y}\tilde{y}% })+i(g_{\tilde{x}\tilde{x}}-g_{\tilde{y}\tilde{y}})\right\} =0, \\[1ex] \xi _{02}&=(f_{\tilde{x}\tilde{x}}-f_{\tilde{y}\tilde{y}}-2g_{\tilde{x}\tilde{% y}})+i(g_{\tilde{x}\tilde{x}}-g_{\tilde{y}\tilde{y}}+2f_{\tilde{x}\tilde{y}% })=0, \\[1ex] \xi _{21}&=\frac{1}{16}\left\{ \left( f_{\tilde{x}\tilde{x}\tilde{x}}+f_{% \tilde{x}\tilde{y}\tilde{y}}+g_{\tilde{x}\tilde{x}\tilde{y}}+g_{\tilde{y}% \tilde{y}\tilde{y}}\right) +i\left( g_{\tilde{x}\tilde{x}\tilde{x}}+g_{% \tilde{x}\tilde{y}\tilde{y}}-f_{\tilde{x}\tilde{x}\tilde{y}}-f_{\tilde{y}% \tilde{y}\tilde{y}}\right) \right\} , \\ \xi _{21}&=-\frac{3}{8}h_{0}\left( e^{2}+\left( k-1\right) ^{2}+i\left( \left( 1-k\right) e-\frac{\left( k-1\right) ^{3}+h_{0}^{3}}{e}\right) \right) , \\[1ex] \operatorname{Re}\left\{ \overline{\lambda }\xi _{21}\right\} &=\frac{3}{8}% h_{0}\left( h_{0}^{3}-\left( k-1\right) ^{2}-e^{2}\right) . \end{align*}% Finally, it is necessary to calculate the following coefficient: \begin{equation*} a\left( h_{0}\right) =\operatorname{Re}\left\{ \frac{(1-2\lambda )\overline{\lambda }% ^{2}}{1-\lambda }\xi _{11}\xi _{20}\right\} +\frac{1}{2}|\xi _{11}|^{2}+|\xi _{02}|^{2}-\operatorname{Re}\left\{ \overline{\lambda }\xi _{21}\right\} , \end{equation*} i.e.,% \begin{equation*} a\left( h_{0}\right) =a\left( \frac{\beta }{\alpha }\right) =\frac{3\beta ^{3}\left( \alpha ^{2}-\beta \right) }{8\alpha ^{4}},\quad 4\alpha -\beta ^{2}>0,\alpha >0,\beta >0, \end{equation*}% which is the coefficient of the cubic term in \cite{garic2017stability} (4.2) in polar coordinates. It is easy to see that% \begin{equation*} a\left( h_{0}\right) =a\left( \frac{\beta }{\alpha }\right) \left\{ \begin{tabular}{lll} $>0$, & for & $\left\{ \begin{array}{l} \alpha \in \left( 0,\sqrt[3]{4}\right) \text{\ and\ }\beta \in \left( 0,\alpha ^{2}\right) \\ \text{or} \\ \alpha \geq \sqrt[3]{4}\text{ and }\beta \in \left( 0,2\sqrt{\alpha }\right) ,% \end{array}% \right. $ \\ $=0,$ & for & $\alpha \in \left( 0,\sqrt[3]{4}\right) $ and $\beta =\alpha ^{2}$, \\ $<0$, & for & $\alpha \in \left( 0,\sqrt[3]{4}\right) $ and $\beta \in \left( \alpha ^{2},2\sqrt{\alpha }\right) $.% \end{tabular}% \right. \end{equation*}% This completes the proof of the theorem. \end{proof} \subsection{Period-doubling bifurcation} The period-doubling bifurcation happens when $\lambda _{i}=-1$ and $|\lambda _{j}|\neq 1$ for $i,j\in \left\{ 1,2\right\} ,i\neq j$. Then, from \eqref{sv_vr} we have $\lambda _{1}=-1$ and $\lambda _{2}=\frac{2-{\alpha h}_{1}^{2}}{2}=3-\beta h_{1}$, for $% {h_{1}:=\frac{\beta -\sqrt{\beta ^{2}-4\alpha }}{\alpha }}$ (and for ${\alpha h}_{1}^{2}=2\beta h_{1}-4$, ${% h}_{1}{\neq \frac{2}{\beta },}${\ }${h}_{1}{\neq \frac{4}{\beta },}${\ }${% \beta \neq 0,}$ ${\beta ^{2}>4\alpha }$, see $(d)$ $2)$ of Proposition \ref% {n4}). Now, we investigate the period-doubling bifurcation of system \eqref{n2} at the origin when parameter $h$ changes in the small neighborhood of $h_{1}$. Therefore, consider system \eqref{n2} by choosing $\tilde{h}$ as a small bifurcation parameter ($h=h_{1}+\tilde{h}$): \begin{equation} x_{n+1}=x_{n}+(h_{1}+\tilde{h})y_{n},\quad y_{n+1}=y_{n}-(h_{1}+\tilde{h})(\alpha x_{n}+\beta y_{n}+x_{n}^{3}+y_{n}^{3}), \label{n18} \end{equation}% Then, system \eqref{n18} can be written in the following form: \begin{equation} \begin{bmatrix} x_{n+1} \\ y_{n+1}% \end{bmatrix}% =% \begin{bmatrix} 1\quad & h_{1} \\ -h_{1}\alpha \quad & 1-h_{1}\beta \end{bmatrix}% \begin{bmatrix} x_{n} \\ y_{n}% \end{bmatrix}% +% \begin{bmatrix} \Phi \left( x_{n},y_{n}\right) \\ \Psi \left( x_{n},y_{n}\right) \end{bmatrix}% , \label{n9} \end{equation}% where \begin{eqnarray*} \Phi \left( x_{n},y_{n}\right) &=&\tilde{h}y_{n} \\ \Psi \left( x_{n},y_{n}\right) &=&-\tilde{h}\left( \alpha x_{n}+\beta y_{n}\right) -\left( h_{1}+\tilde{h}\right) \left( x_{n}^{3}+y_{n}^{3}\right) , \end{eqnarray*}% and $|\tilde{h}|<<1$ is a small perturbation parameter. To make the normal form of system \eqref{n9}, we need an invertible matrix $B$ whose columns are eigenvectors of the corresponding $\lambda _{1}=-1$ and $\lambda _{2}=3-\beta h_{1}$: \begin{equation*} B=% \begin{bmatrix} h_{1} & h_{1} \\ -2 & \lambda _{2}-1% \end{bmatrix},% \end{equation*} where \begin{equation*} B^{-1}=\frac{1}{h_{1}\left( \lambda _{2}+1\right) }% \begin{bmatrix} \lambda _{2}-1 & -h_{1} \\ 2 & h_{1}% \end{bmatrix},% \quad B^{-1}% \begin{bmatrix} 1\quad & h_{1} \\ -h_{1}\alpha \quad & 1-h_{1}\beta \end{bmatrix}% B=% \begin{bmatrix} -1 & 0 \\ 0 & \lambda _{2}% \end{bmatrix}% . \end{equation*}% Now, let \begin{equation*} \begin{bmatrix} x_{n} \\ y_{n}% \end{bmatrix}% =B% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% =% \begin{bmatrix} h_{1} & h_{1} \\ -2 & \lambda _{2}-1% \end{bmatrix}% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% =% \begin{bmatrix} h_{1}\left( \tilde{x}_{n}+\tilde{y}_{n}\right) \\ \tilde{y}_{n}\left( \lambda _{2}-1\right) -2\tilde{x}_{n}% \end{bmatrix}% . \end{equation*}% Then, system \eqref{n9} changes into% \begin{equation*} B% \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =% \begin{bmatrix} 1\quad & h_{1} \\ -h_{1}\alpha \quad & 1-h_{1}\beta \end{bmatrix}% B% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +% \begin{bmatrix} \Phi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) \\ \Psi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) \end{bmatrix}% , \end{equation*} where% \begin{align*} \Phi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) =&\tilde{h}\left( \tilde{y}% _{n}\left( \lambda _{2}-1\right) -2\tilde{x}_{n}\right) , \\ \Psi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) =&-\tilde{h}\left( \alpha h_{1}\left( \tilde{x}_{n}+\tilde{y}_{n}\right) +\beta \left( \tilde{y}% _{n}\left( \lambda _{2}-1\right) -2\tilde{x}_{n}\right) \right) \\ &-\left( h_{1}+\tilde{h}\right) \left( \left( h_{1}\left( \tilde{x}_{n}+% \tilde{y}_{n}\right) \right) ^{3}+\left( \tilde{y}_{n}\left( \lambda _{2}-1\right) -2\tilde{x}_{n}\right) ^{3}\right) , \end{align*}% which implies that \begin{equation*} \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =B^{-1}% \begin{bmatrix} 1\quad & h_{1} \\ -h_{1}\alpha \quad & 1-h_{1}\beta \end{bmatrix}% B% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +B^{-1}% \begin{bmatrix} \Phi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) \\ \Psi _1 \left( \tilde{x}_{n},\tilde{y}_{n}\right) \end{bmatrix},% \end{equation*}% i.e., \begin{equation} \begin{bmatrix} \tilde{x}_{n+1} \\ \tilde{y}_{n+1}% \end{bmatrix}% =% \begin{bmatrix} -1 & 0 \\ 0 & \lambda _{2}% \end{bmatrix}% \begin{bmatrix} \tilde{x}_{n} \\ \tilde{y}_{n}% \end{bmatrix}% +% \begin{bmatrix} f(\tilde{x}_{n},\tilde{y}_{n}) \\ g(\tilde{x}_{n},\tilde{y}_{n})% \end{bmatrix}% , \label{n10} \end{equation}% where% %\begin{equation*} %\begin{bmatrix} %f(\tilde{x}_{n},\tilde{y}_{n}) \\ %g(\tilde{x}_{n},\tilde{y}_{n})% %\end{bmatrix}% %=B^{-1}% %\begin{bmatrix} %\Phi \left( \tilde{x}_{n},\tilde{y}_{n}\right) \\ %\Psi \left( \tilde{x}_{n},\tilde{y}_{n}\right) %\end{bmatrix}% %, %\end{equation*}% %and% \begin{align*} f(\tilde{x},\tilde{y}) =&\dfrac{1}{\lambda _{2}+1}\Big[\left( \dfrac{-2% \tilde{h}\left( \lambda _{2}+1\right) }{h_{1}}\right) \tilde{x} \\ & +\left( h_{1}+\tilde{h}\right) \left( \left( h_{1}^{3}-8\right) \tilde{x}% ^{3}+\left( \left( \lambda _{2}-1\right) ^{3}+h_{1}^{3}\right) \tilde{y}% ^{3}\right) \\ & +3\left( h_{1}+\tilde{h}\right) \left( \left( h_{1}^{3}-2\left( \lambda _{2}-1\right) ^{2}\right) \tilde{x}\tilde{y}^{2}+\left( h_{1}^{3}+4\left( \lambda _{2}-1\right) \right) \tilde{x}^{2}\tilde{y}\right) \Big], \end{align*}% \begin{align*} g(\tilde{x},\tilde{y}) =&\dfrac{1}{\lambda _{2}+1}\Big[\left( \dfrac{\left( 2-\beta h_{1}\right) \tilde{h}\left( \lambda _{2}+1\right) }{h_{1}}\right) \tilde{y} \\ & -\left( h_{1}+\tilde{h}\right) \left( \left( h_{1}^{3}-8\right) \tilde{x}% ^{3}+\left( \left( \lambda _{2}-1\right) ^{3}+h_{1}^{3}\right) \tilde{y}% ^{3}\right) \\ & +3\left( h_{1}+\tilde{h}\right) \left( \left( 2\left( \lambda _{2}-1\right) ^{2}-h_{1}^{3}\right) \tilde{x}\tilde{y}^{2}-\left( h_{1}^{3}+4\lambda _{2}-4\right) \tilde{x}^{2}\tilde{y}\right) \Big]. \end{align*}% Now, we compute the center manifold $W^{c}(0,0)$ at the fixed point $(0,0)$ in a small neighborhood of $\tilde{h}=0$ by using the center manifold theorem in \cite{gukenheimer1983nonlinear,zbMATH02120913}. To compute the center manifold, one can consider it by a polynomial of $\tilde{x}$ and $\tilde{h}$ of degree four $ \tilde{y}=\mathcal{C}\left( \tilde{x}\right)$, where \begin{align*} \mathcal{C}\left( \tilde{x}\right)& =A\tilde{x}^{2}+B\tilde{x}\tilde{h}+C\tilde{h}^{2}+D\tilde{x}^{3}+E\tilde{x}% ^{2}\tilde{h}+F\tilde{x}\tilde{h}^{2} \\ & +\sigma \tilde{h}^{3}+H\tilde{x}% ^{4}+\omega \tilde{x}^{3}\tilde{h}+J\tilde{x}^{2}\tilde{h}^{2}+N\tilde{x}% \tilde{h}^{3}+M\tilde{h}^{4}, \end{align*}% with the unknown coefficients $A$, $B$, $C$, $D$, $E$, $F$, $\sigma $, $H $, $\omega $, $J$, $N$, $M$. Then, it must satisfy the following condition: \begin{equation*} \mathcal{C}\left( -\tilde{x}+f\left( \tilde{x},\mathcal{C}\left( \tilde{x}% \right) \right) \right) -\lambda _{2}\mathcal{C}\left( \tilde{x}\right) -g\left( \tilde{x},\mathcal{C}\left( \tilde{x}\right) \right) =0. \end{equation*}% By approximate computation, for the center manifold, we obtain% \begin{equation*} A=B=C=E=F=G=H=J=N=M=0, \end{equation*}% \begin{equation} \begin{array}{l} \sigma =\dfrac{h_{1}\left( 2-h_{1}\right) \left( 2h_{1}+h_{1}^{2}+4\right) }{% \left( \lambda _{2}+1\right) ^{2}},\quad \omega =\dfrac{2\left( h_{1}-2\right) \left( \lambda _{2}+3\right) \left( 2h_{1}+h_{1}^{2}+4\right) }{\left( \lambda _{2}+1\right) ^{3}}.% \end{array} \label{CenMan} \end{equation}% Therefore, the center manifold% \begin{equation*} W^{c}\left( 0,0\right) =\left\{ \left. \left( \tilde{x},\tilde{y}\right) \right\vert \;\tilde{y}=\mathcal{C}\left( \tilde{x}\right) ,\;\mathcal{C}% \left( 0\right) =\mathcal{C}^{\prime }(0)=0\right\}, \end{equation*}% near $\tilde{h}=0$, can be approximately written in the following form: \begin{equation*} W^{c}\left( 0,0\right) =\left\{ \left. \left( \tilde{x},\tilde{y}\right) \right\vert \;\tilde{y}=\sigma \tilde{x}^{3}+\omega \tilde{x}^{3}\tilde{h}+% \mathcal{O}\left( \left( |\tilde{x}|+|\tilde{h}|\right) ^{5}\right) \right\} , \end{equation*}% where $\sigma $ and $\omega $ are given by \eqref{CenMan}. Now, we restrict system (\ref{n10}) (the first equation) to the center manifold: \begin{equation*} \widetilde{F}\left( \tilde{x}\right) =-\tilde{x}+f\left( \tilde{x},\mathcal{C% }(\tilde{x})\right) =-\tilde{x}+f\left( \tilde{x},\sigma \tilde{x}% ^{3}+\omega \tilde{x}^{3}\tilde{h}\right) , \end{equation*}% i.e., \begin{equation*} \widetilde{F}\left( \tilde{x}\right) =-\tilde{x}+\delta \tilde{x}\tilde{h}% +\eta \tilde{x}^{3}+\mathcal{O}\left( \left( |\tilde{x}|+|\tilde{h}|\right) ^{4}\right), \end{equation*}% where $\delta =-\frac{2}{h_{1}}$ and $\eta =\frac{\left( h_{1}+\tilde{h}% \right) \left( h_{1}-2\right) \left( h_{1}^{2}+2h_{1}+4\right) }{\lambda _{2}+1}$. \bigskip \\ The following three conditions are satisfied (see \cite[Theorem 3.1, p.~246]{zbMATH00918596}): \begin{itemize} \item $\widetilde{F}(0)=0$, \item $\widetilde{F}^{\prime}(0)=-1$ and \item $\kappa _{1}=\left. \left( \frac{\partial ^{2}\widetilde{F}}{\partial \tilde{x}\partial \tilde{h}}+\frac{1}{2}\frac{\partial \widetilde{F}}{% \partial \tilde{h}}\frac{\partial ^{2}\widetilde{F}}{\partial \tilde{x}^{2}}% \right) \right\vert _{\left( 0,0\right) }=\delta =-\tfrac{2}{h_{1}}\neq 0$, \\ $\kappa _{2}=\left. \left( \frac{1}{6}\frac{\partial ^{3}\widetilde{F}}{\partial \tilde{x}^{3}}+\left( \frac{1}{2}\frac{\partial ^{2}\widetilde{F}}{\partial \tilde{x}^{2}}\right) ^{2}\right) \right\vert _{\left( 0,0\right) }=\frac{h_{1}\left( h_{1}-2\right) \left( h_{1}^{2}+2h_{1}+4\right) }{\lambda _{2}+1}\neq 0$, when $h_{1}\neq 2$. \end{itemize} %when $h_{1}\neq 2$. With the above arguments, we have the following theorem, which shows that system \eqref{n2} experiences a period-doubling bifurcation at the origin, as the parameter $\tilde{h}$ varies through the bifurcation value $\tilde{h}% =0$. \begin{theorem} System \eqref{n2} undergoes a period-doubling bifurcation at the fixed point $E_0=(0,0)$ if $h=h_{1}:=\frac{\beta-\sqrt{\beta ^{2}-4\alpha}}{\alpha}\neq2 $ and $\beta^{2}-4\alpha>0 $. Furthermore, the orbits of period-$2$ that bifurcate from the origin are stable if $h_{1}<2$, and unstable if $h_{1}>2$ for $\beta ^{2}-4\alpha>0$. \end{theorem} \section{The existence of Li-Yorke chaos } \label{Sec:4} In this section, by using Marotto's theorem \cite{zbMATH07626819,axioms13040280,MAROTTO200525,MAROTTO1978199}, we prove that system \eqref{n2} has chaotic behavior for particular values of coefficients. In \cite{MAROTTO1978199}, Marotto extended the definition of chaos in the sense of Li-Yorke from one dimension to multiple dimensions by introducing the concept of a snap-back repeller. In \cite{MAROTTO200525}, Marotto redefined the definition of a snap-back repeller. He proved that the existence of a snap-back repeller is an indicator of chaos. We present the concept of a snap-back repeller and Marotto's theorem. We will prove that system \eqref{n2} possesses chaotic behavior in the sense of Marotto's definition whenever the origin is a source and has a pair of complex multipliers. From (\ref{FE_0}) we see that it is $\left\vert \lambda _{1,2}\right\vert >1$ if $C>1$ and $B^{2}-4C<0$, i.e., if $h\alpha -\beta >0\text{ and }% 4\alpha -\beta ^{2}>0$. Our next step is to determine a neighborhood $% U_{r}\left( 0,0\right) $ in which the norms of eigenvalues exceed 1 for all $% \left( x,y\right) \in U_{r}\left( 0,0\right) $. Then from (\ref{CPT}) conditions $C>1$ and $B^{2}-4C<0$ are equivalent to the conditions % (see \eqref{JT})% \begin{align*} \left( 3hy^{2}+h\beta -2\right) ^{2}-4\left( 3h^{2}x^{2}+\alpha h^{2}-3hy^{2}-\beta h+1\right) & <0, \\ 3h^{2}x^{2}+\alpha h^{2}-3hy^{2}-\beta h+1& >1, \end{align*}% or \begin{align*} \delta \left( x\right) &=12x^{2}-9y^{4}-6y^{2}\beta +4\alpha -\beta ^{2}>0, \\[1ex] \gamma \left( x\right) &=3hx^{2}-3y^{2}+h\alpha -\beta >0. \end{align*} These are two hyperbolas that have local extremes at zero.% \begin{align*} \delta _{\min }\left( x\right) & =\delta \left( 0\right) =-9y^{4}-6y^{2}\beta +4\alpha -\beta ^{2}>0, \\ \gamma _{\min }\left( x\right) & =\gamma \left( 0\right) =-3y^{2}+h\alpha -\beta >0. \end{align*}% For $4\alpha -\beta ^{2}>0$ and $h\alpha -\beta >0$, we have% \begin{align*} \delta _{\min }\left( x\right) & >0\Leftrightarrow y\in \left( -\sqrt{\tfrac{% 2\sqrt{\alpha }-\beta }{3}},\sqrt{\tfrac{2\sqrt{\alpha }-\beta }{3}}\right) , \\ \gamma _{\min }\left( x\right) & >0\Leftrightarrow y\in \left( -\sqrt{\tfrac{% h\alpha -\beta }{3}},\sqrt{\tfrac{h\alpha -\beta }{3}}\right) , \end{align*}% and $r=\min \left\{ \sqrt{\tfrac{2\sqrt{\alpha }-\beta }{3}},\sqrt{\tfrac{% h\alpha -\beta }{3}}\right\} $, $U_{r}\left( 0,0\right) =\left\{ \left. \left( x,y\right) \right\vert \text{ }x^{2}+y^{2}0 $, we have $x_{1}=\pm \sqrt{\tfrac{\left( h\left( h\alpha -\beta \right) +1\right) h}{1-h^{3}}}\neq 0$. For% \begin{equation*} \left( x_{1},y_{1}\right) =\left( \sqrt{\tfrac{h^{2}\left( h\alpha -\beta \right) +h}{1-h^{3}}},-\frac{x_{1}}{h}\right) \neq \left( 0,0\right) , \end{equation*}% we have $|J_{T}\left( x_{1},y_{1}\right) |\neq 0$ because it is% \begin{equation*} |J_{T}\left( x_{1},y_{1}\right) |=\left. -h\beta -3hy^{2}+h^{2}\alpha +3h^{2}x^{2}+1\right\vert _{\left( x_{1},y_{1}\right) }=-2\left( h\left( h\alpha -\beta \right) +1\right) <0. \end{equation*}% Using $y_{1}=-\frac{x_{1}}{h}$ from (\ref{T1}) we get that $x_{0}$ is the solution of the equation% \begin{equation} Ax_{0}^{3}+Bx_{0}^{2}+Cx_{0}+D=0, \label{Snap} \end{equation}% where $A=\left( 1-h^{3}\right) >0$, $B=-3x_{1}<0$, \begin{align*} C& =-\alpha h^{3}+\beta h^{2}-h+3x_{1}^{2}=\dfrac{h\left( h^{3}+2\right) \left( 1+h\left( h\alpha -\beta \right) \right) }{\left( 1-h\right) \left( h+h^{2}+1\right) }>0, \\ D& =-\beta h^{2}x_{1}+2hx_{1}-x_{1}^{3}=\dfrac{\left( h^{4}\beta -2h^{3}-h^{2}\alpha +1\right) x_{1}h}{\left( 1-h\right) \left( h+h^{2}+1\right) }, \end{align*}% for $x_{1}=\sqrt{\frac{h^{2}\left( h\alpha -\beta \right) +h}{1-h^{3}}}$ and a small enough $h$. \begin{example} For $\alpha =20$, $\beta =2$ and $h=0.2$, we find a neighborhood% \begin{equation*} U_{r}\left( 0,0\right) =\left\{ \left( x,y\right) :x^{2}+y^{2}1$, $\left\vert \lambda _{2}\right\vert =0.81246<1$. \newline Equation \eqref{Snap} has three real solutions: \begin{equation*} -0.03237199632171373`,0.6022113936915484`,1.0368530868914105`. \end{equation*}% Only for $x_{0}=0.6022113936915484`$ it is true that% \begin{equation*} \left( x_{0},y_{0}\right) =\left( 0.6022113936915484`,-0.3546587278124836`\right) \in U_{r}\left( 0,0\right) , \end{equation*}% \begin{equation*} T^{2}\left( x_{0},y_{0}\right) =\left( 0.,-8.881784197001252`\cdot 10^{-16}\right) \approx \left( 0,0\right) , \end{equation*} \begin{equation*} |J_{T}\left( x_{1},y_{1}\right) |\approx 1.368\,1\neq 0 \end{equation*}% and \begin{equation*} \left\vert \lambda _{1,2}\right\vert =\left\vert 0.76227\pm 0.88713i\right\vert =1.1696>1. \end{equation*}% Then $\left( x_{0},y_{0}\right) =\left( 0.6022113936915484`,-0.3546587278124836`\right) $ is a snap-back repeller and in this case, system \eqref{n2} is chaotic in the sense of Li-Yorke. \end{example} \section{Numerical simulations} To justify the theoretical analysis for system \eqref{n2}, we present the phase portraits for $\alpha =10$, $\beta =0.1$, and initial conditions $% (10^{-3},10^{-3})$ by choosing a different value of parameter $h$ ($h=0.001$% , $h=0.1$, and $h=0.01$). Based on Proposition \ref{n4}, when $\alpha <\frac{% \beta }{h}$, $\alpha >\frac{\beta }{h}$ and $\alpha =\frac{\beta }{h}$, the fixed point $(0,0)$ is a sink, source, and non-hyperbolic (NH), respectively. See Theorem \ref{Thh1} (i), for $\alpha =10> \sqrt[3]{4}$, $\beta =0.1\in \left( 0,2% \sqrt{\alpha }\right) $, $h_{0}=\frac{\beta }{\alpha }=0.01$ and $a\left( h_{0}\right) =\frac{3\beta ^{3}\left( \alpha ^{2}-\beta \right) }{8\alpha ^{4}}=\frac{2997}{800\,000\,000}>0.$\bigskip For illustration see: \begin{itemize} \item Figure \ref{f1s}: 10000 and 100000 iterations for $% (x_{0},y_{0})=(0.001,0.001)$ and $hh_{0}$; \item Figure \ref{f1NH}: 350 iterations for $(x_{0},y_{0})=(0.001,0.001)$ and ellipse $a^{2}x^{2}+b^{2}y^{2}=a^{2}b^{2}$ with $a=0.0033232476573834087$ and $b=0.0010486001027132572$ for $h=h_{0}$. \end{itemize} \begin{figure}[h!] \centering \includegraphics[width=6cm]{fig4a} % \includegraphics[width=6cm]{fig4b} \caption{The phase portrait for $\protect\alpha =10<\frac{\protect\beta }{h}=% \frac{0.1}{0.001}=100$, (a sink for $h=0.001$) for $h\frac{\protect\beta }{h}=% \frac{0.1}{0.1}=1$, (a source for $h=0.1$) for $h>h_{0}$} \label{f1r} \end{figure} \newpage \noindent \begin{figure}[h!] \centering \includegraphics[width=6.5cm]{fig6a} % \includegraphics[width=6.5cm]{fig6b} \caption{The phase portrait for $\protect\alpha =\frac{\protect\beta }{h}=% \frac{0.1}{0.01}=10$, (NH) and ellipse $% a^{2}x^{2}+b^{2}y^{2}=a^{2}b^{2}$ with $a=0.0033232476573834087$ and $% b=0.0010486001027132572$ for $h=h_{0}$} \label{f1NH} \end{figure} \begin{figure}[h!] \centering \includegraphics[width=8.5cm]{fig7} \caption{Bifurcation diagram generated by Code Bif2D \protect\cite{ufuktepe2014applications}} \label{f2xy} \end{figure} \begin{figure}[h!] \centering \includegraphics[width=8.5cm]{fig8} \caption{The Lyapunov exponent} \label{f2lce} \end{figure} \begin{figure}[h!] \centering \includegraphics[width=6.5cm]{fig9a} % \includegraphics[width=6.5cm]{fig9b} \caption{The snap-back repeller and chaotic attractor} \label{spr} \end{figure} \newpage Equation \eqref{Snap} has at least one real solution (it can have three real solutions). If there is such $x_{0}$ that $\left( x_{0},y_{0}\right) \in U_{r}\left( 0,0\right) $, $|J_{T}\left( x_{1},y_{1}\right) |\neq 0$ and $\left( x_{1},y_{1}\right) \notin U_{r}\left( 0,0\right) $, then $\left( x_{0},y_{0}\right) $ is a snap-back repeller. \noindent In Figure \ref{f2xy}, we plot the bifurcation diagram of system % \ref{n2} in $(\alpha,x)$, $(\alpha,y)$ plane for $h=0.5$, $\beta=0.1$, and initial conditions $(10^{-3},10^{-3})$ are generated by Code Bif2D \cite{ufuktepe2014applications}. Finally, the diagram of the Lyapunov exponent for $% h=0.5$, $\alpha<2$, $\beta=0.1$, and initial conditions $(10^{-3},10^{-3})$ are shown in Figure \ref{f2lce}. The Lyapunov exponents (based on the computational algorithm in \cite{he2016numerical,sandri1996numerical}) corresponding to Figure \ref{f2xy} are shown in Figure \ref{f2lce}, verifying the existence of chaos. In Section $4$, we proved that system % \eqref{n2} has chaotic behavior in the sense of Marotto's definition. So, we guess it has chaotic behavior in the other sense. Namely, the point \begin{equation*} \left( x_{0},y_{0}\right) =\left(0.6022113936915484`,-0.3546587278124836`\right) \end{equation*} is a snap-back repeller for $\alpha=20$, $\beta=2$ and $h=0.2$. In this case, system \eqref{n2} is chaotic in the sense of Li-Yorke. For $\left( x_{0},y_{0}\right) =\left( 0.6,-0.3\right) $, we have a chaotic attractor (see Figure \ref{spr}). \section{Discussion} The general form of the Rayleigh-Duffing oscillator, \begin{equation} \label{z1} \ddot{x} + a \dot{x} + b \dot{x}^3 + c x + d x^3 = 0, \end{equation} was investigated by Chen et al. in \cite{chen2018saddle}. This model incorporates both nonlinear damping and external periodic forcing, and its complex dynamical behavior, including bifurcations (pitchfork, Hopf, and heteroclinic), chaotic attractors, and multistability, is analyzed using both analytical and numerical techniques. The continuous system \eqref{n1} is obtained from system \eqref{z1} for $\dot{x}=y$, $c=\alpha$, $a=\beta$, $d=b=1$, and it exhibits Hopf bifurcation, a homoclinic bifurcation, and a double-limit cycle bifurcation. In particular, the homoclinic bifurcation, also known as a gluing bifurcation, is notable for marking a critical point where two symmetric periodic orbits merge into a single, more complex orbit. This process leads to a sudden qualitative change in the system's dynamics, often associated with the emergence or disappearance of chaotic attractors, and plays a crucial role in shaping the system's global bifurcation structure (see \cite{chen2016global}). By comparing the continuous system \eqref{n1} with its discrete counterpart \eqref{n2}, we show that the fixed point at the origin $(0,0)$ undergoes both a Neimark–Sacker bifurcation and a period-doubling bifurcation. Moreover, when the origin of system \eqref{n2} is a source with a pair of complex multipliers, the system exhibits chaotic behavior in the sense of Marotto's definition, an extension of Li–Yorke chaos from one-dimensional to multidimensional systems via the concept of a snap-back repeller. \section*{Acknowledgements} The authors sincerely thank the anonymous referees for their valuable comments and suggestions, which significantly improved this paper. %%===============================================================%% %% REFERENCES %% %% References should be provided in bibtex file. %% %% We suggest using MR Lookup for finding bibtex entries. %% %%===============================================================%% \bibliography{references} \end{document} \begin{thebibliography}{99} \bibitem {1} {\sc M. Bikdash, B. Balachandran, A.H. Nayfeh}, {\em Melnikov analysis for a ship with a general Roll-damping model}, Nonlinear Dyn. {\bf 6}(1994), 101--124. https://doi.org/10.1007/BF00045435 \bibitem {2} {\sc H. Chen, L. Zou}, {\em Global study of Rayleigh-Duffing oscillators}, J. Phys. A: Math. Theor. {\bf 49}(2016), 1--25. % https://ui.adsabs.harvard.edu/link_gateway/2016JPhA...49p5202C/doi:10.1088/1751-8113/49/16/165202 \bibitem {3} {\sc H. Chen, D. Huang, Y. Jian}, {\em The saddle case of Rayleigh-Duffing oscillators}, J. Nonlinear Dyn. 93, 2283–2300 (2018). https://doi.org/10.1007/s11071-018-4325-8 \bibitem {4} {\sc E. Freire, E. Gamero, A. Rodr\u{A}-guez-Luis}, {\em Study of a degenerate Bogdanov-Takens bifurcation in a family of mechanical oscillator}, Mech. Res. Commun. {\bf 25}(1998), 287--297. \bibitem {5} {\sc M. Gari\'{c}-Demirovi\'{c}, S. Hrusti\'{c}, S. Moranjki\'{c}, M. Nurkanovi\'{c}, Z. Nurkanovi\'{c}}, {\em The existence of Li-Yorke chaos in certain predator-prey system of difference equations}, Sarajevo Journal of Mathematics {\bf 18(31)(1)}(2022), 45--62. {https://doi.org/10.5644/SJM.18.01.04} \bibitem {5a} {\sc M. Gari\'{c}-Demirovi\'{c}, M.R.S. Kulenovi\'{c}, M. Nurkanovi\'{c}, Z. Nurkanovi\'{c}}, {\em The Existence of Li–Yorke Chaos in a Discrete-Time Glycolytic Oscillator Model.} Axioms {\bf 13(4)}(2024), 1--17. https://doi.org/10.3390/axioms13040280 \bibitem {6} {\sc M. Gari\'{c}-Demirovi\'{c}, M. Nurkanovi\'{c}, Z. Nurkanovi\'{c}}, {\em Stability, periodicity and Neimark-Sacker bifurcation of certain homogeneous fractional difference equations}, International Journal of Difference Equations. {\bf 12(1)}(2017), 27--53. \bibitem {7} {\sc J. Guckenheimer, P. Holmes}, {\em Nonlinear oscillations, dynamical system and bifurcation of vector fields}, New York: Springer-Verlag 1983. \bibitem {8} {\sc J.K. Hale, H. Buttanri, H. Kocak}, {\em Dynamics and Bifurcations. Texts in Applied Mathematics}, Springer, New York 2012. \bibitem {9} {\sc W. Herfort, H. Troger}, {\em Robust modelling of flow induced oscillations of bluff bodies}, Math. Modelling Sci. Tech. {\bf 8}(1987), 251--255. \bibitem {10} {\sc J. He, S. Yu, J. Cai}, {\em Numerical Analysis and Improved Algorithms for Lyapunov-Exponent Calculation of Discrete-Time Chaotic Systems}, International Journal of Bifurcation and Chaos {\bf 26}(2016), 1650219-1-11. \bibitem {10a} {\sc S. Hrusti\'{c}, S. Moranjki\'{c}, Z. Nurkanovi\'{c}}, {\em Local Stability, Period-Doubling and 1:2 Resonance Bifurcation for a Discretized Chemical Reaction System }, MATCH Commun. Math. Comput. Chem. {\bf 93} (2025) 349–378. doi: 10.46793/match.93-2.349H \bibitem {11} {\sc Y. Kanai, H. Yabuno}, {\em Creation-annihilation process of limit cycles in the Rayleigh-Duffing oscillator}, Nonlinear Dyn. {\bf 70(2)}(2012), 1007--1016. \bibitem{12} {\sc M.R.S. Kulenovi\'{c}, O. Merino}, {\em Discrete Dynamical Systems and Difference Equations with Mathematica}, Chapman and Hall/CRC: Boca Raton, FL, USA; London, UK, 2002. https://doi.org/10.1201/9781420035353 \bibitem{13} {\sc M.R.S. Kulenovi\'{c}, S. Moranjki\'{c}, M. Nurkanovi\'{c}, Z. Nurkanovi\'{c}}, {\em Global asymptotic stability and Naimark-Sacker bifurcation of certain mix monotone difference equation}, Discrete Dyn. Nat. Soc. (2018), 1--22. https://doi.org/10.1155/2018/7052935 \bibitem{14} {\sc M.R.S. Kulenovi\'{c}, S. Moranjki\'{c}, Z. Nurkanovi\'{c}}, {\em Naimark-Sacker bifurcation of second order rational difference equation with quadratic terms}, J. Nonlinear Sci. Appl. {\bf 10(7)}(2017), 3477--3489. http://dx.doi.org/10.22436/jnsa.010.07.11 \bibitem{15} {\sc M.R.S. Kulenovi\'{c}, S. Moranjki\'{c}, Z. Nurkanovi\'{c}}, {\em Global dynamics and bifurcation of a perturbed sigmoid Beverton-Holt difference equation}, Math. Methods Appl. Sci. {\bf 39(10)}(2015), 2696--2715. https://doi.org/10.1002/mma.3722 \bibitem {16} {\sc Y. Kuznetsov}, {\em Elements of Applied Bifurcation Theory}, Springer, 2004. \bibitem {17} {\sc X. Liu, X. Dongmei}, {\em Complex dynamic behaviors of a discrete-time predatorprey system}, Chaos Solitons Fractals {\bf 32(1)}(2007), 80--94. \bibitem {18} {\sc A.C.J. Luo}, {\em Discretization and implicit mapping dynamics}, Nonlinear Phys. Sci. Higher Education Press; Berlin: Springer, Beijing, 2015. http://dx.doi.org/10.1007/978-3-662-47275-0 \bibitem {19} {\sc F.R. Marotto}, {\em Snap-Back Repellers Imply Chaos in $\mathbb{R}^{n} $}, J. Math. Anal. and Appl. {\bf 63}(1978), 199--223. {http://doi.org/10.1016/0022-247X(78)90115-4} \bibitem {20} {\sc F.R. Marotto}, {\em On redefining a snap-back repeller}, Chaos Solitons Fractals {\bf 25(1)}(2005), 25--28. http://dx.doi.org/10.1016/j.chaos.2004.10.003 \bibitem {21} {\sc A.H. Nayfeh, D.T. Mook}, {\em Nonlinear Oscillations, Pure and Applied Mathematics}, A Wiley-Interscience Publication. New York etc.: John Wiley \& Sons. XIV, {\bf 704}(1979), 21--35. \bibitem {22} {\sc Z. Nurkanovi\'{c}, M. Nurkanovi\'{c}, M. Gari\'{c}-Demirovi\'{c}}, {\em Stability and Neimark-Sacker Bifurcation of Certain Mixed Monotone Rational Second-Order Difference Equation}, Qual. Theory Dyn. Syst. {\bf 20(3)}(2021), 1--41. https://doi.org/10.1007/s12346-021-00515-4 \bibitem {23} {\sc C. Robinson}, {\em Dynamical systems, stability, symbolic dynamics and chaos}, FL: CRC Press, Boca Raton, 1995. \bibitem {24} {\sc M. Sandri}, {\em Numerical calculation of Lyapunov exponents}, The Mathematica Journal {\bf 6(3)}(1996), 78--84. https://www.mathematica-journal.com/issue/v6i3/article/sandri/contents/63sandri.pdf \bibitem {25} {\sc \"{U}. Ufuktepe, S. Kap\c{c}ak}, {\em Applications of Discrete Dynamical Systems with Mathematica}, Kurenai {\bf 1909}(2014), 207--216. https://hdl.handle.net/2433/223175 \bibitem {26} {\sc R. Yamapi}, {\em Synchronization dynamics in a ring of four mutually inertia coupled self-sustained electrical systems}, Physica A-statistical Mechanics and Its Applications {\bf 366}(2006), 187--196. \end{thebibliography}