# Small- and large-amplitude limit cycles in Kukles systems with algebraic
invariant curves

José Mujica  Department of Mathematics, Vrije Universiteit Amsterdam, The
Netherlands<EMAIL_ADDRESS>
(Date: January 2022.)

###### Abstract.

Limit cycles of planar polynomial vector fields have been an active area of
research for decades; the interest in periodic-orbit related dynamics comes
from Hilbert’s 16th problem and the fact that oscillatory states are often
found in applications. We study the existence of limit cycles and their
coexistence with invariant algebraic curves in two families of Kukles systems,
via Lyapunov quantities and Melnikov functions of first and second order. We
show center conditions, as well as a connection between small- and large-
amplitude limit cycles arising in one of the families, in which the first
coefficients of the Melnikov function correspond to the first Lyapunov
quantities. We also provide an example of a planar polynomial system in which
the cyclicity is not fully controlled by the first nonzero Melnikov function.

###### Key words and phrases:

Limit cycles, Kukles systems, Lyapunov quantities, Melnikov functions,
Invariant algebraic curves

###### 2000 Mathematics Subject Classification:

Primary 34C07, 34C05; Secondary 34C25

##  1\. Introduction

We consider planar systems of the form

(1.1) $\begin{array}[]{rcl}\dot{x}&=&-y\\\\[8.5359pt]
\dot{y}&=&f(x,y),\end{array}$

where $f(x,y)$ is a polynomial function with real coefficients and not
divisible by $y$. Systems like (1.1) are known as _Kukles_ systems, named
after _Isaak Solomonovich Kukles_ (1905-1977), who showed necessary and
sufficient center conditions for a cubic system of the form

(1.2) $\left\\{\begin{array}[]{rcl}\dot{x}&=&-y,\\\\[8.5359pt]
\dot{y}&=&x+a_{1}x^{2}+a_{2}xy+a_{3}y^{2}+a_{4}x^{3}+a_{5}x^{2}y+a_{6}xy^{2}+a_{7}y^{3},\end{array}\right.$

where $a_{1},\ldots,a_{7}$ are real coefficients. Finding necessary and
sufficient center conditions for a planar system means to find conditions on
the parameters so that the system has an equilibrium point (usually the
origin) with a Jacobian matrix having complex eigenvalues with zero real part
and all the orbits in a neighborhood of the equilibrium point are
topologically equivalent to a circle. Finding such conditions, specially the
second one, is far from being straight since a center is a codimension-
infinity phenomenom that depends highly on the system’s nature. The conditions
given by Kukles for system (1.1) were:

$\begin{array}[]{rcl}k_{\alpha}=k_{\beta}=k_{\gamma}=a_{7}=0&&(K1),\\\\[8.5359pt]
a_{7}=a_{2}=a_{5}=0&&(K2),\\\\[8.5359pt]
a_{7}=a_{5}=a_{3}=a_{1}=0&&(K3),\\\\[8.5359pt]
k_{\alpha}=k_{\beta}=k_{\gamma}=k_{\delta}=0&&(K4),\end{array}$

with

$\begin{array}[]{rcl}k_{\alpha}&=&a_{4}a_{2}^{2}+a_{5}\lambda\\\\[8.5359pt]
k_{\beta}&=&(3a_{7}\lambda+\lambda^{2}+a_{6}a_{2}^{2})a_{5}-3a_{7}\lambda^{2}-a_{6}a_{2}^{2}\lambda\\\\[8.5359pt]
k_{\gamma}&=&\lambda+a_{1}a_{2}+a_{5}\\\\[8.5359pt]
k_{\delta}&=&9a_{6}a_{2}^{2}+2a_{4}^{2}+9\lambda^{2}+27a_{7}\lambda\end{array}$

and

$\lambda=a_{2}a_{3}+3a_{7}.$

Conditions $(K1)$–$(K4)$ are known as Kukles conditions [29].

Since Kukles’ work more researchers became interested in Kukles systems, due
to its connection with Hilbert’s 16th problem and the fact that Kukles systems
are often found in applications via Liénard equations; they can model several
oscillatory phenomena [5, 29]. In particular, these systems appear in models
for heart rythm and breathing patterns [37, 38], neuronal activity [13, 30],
chemical reactions [12], microtube kinetics and self organization on cellular
level [18], as well as dynamics of populations of neurons, genetic models and
biological regulation [36, 41].

Using computational techniques, Jin y Wang [28] computed the focal quantities
of (1.2) and noticed that there was a condition missed by Kukles, which
provided a weak focus of much higher order than expected for these systems.
Such conditions were

$a_{2}=a_{6}=0;\ a_{3}=-2a_{1};\ a_{5}=-3a_{7};\ a_{7}^{2}=a_{4}^{2};\
-a_{4}=\dfrac{a_{1}^{2}}{3}.$

In [7], Cristopher and Lloyd proved that Jin and Wang’s conditiones were
indeed center conditions, claiming that Kukles conditions where incomplete.
They classify all center conditions when $a_{7}=0$, verifying the case
completely studied by Kukles. They also proved that for $a_{7}=0$ the maximum
order of a weak focus is five. In [9] the case $a_{7}\neq 0$ is studied,
together with the center conditions not considered previously by Kukles, Jin
and Wang.

As a remark, it is known that there exists a one-to-one correspondence between
Kukles systems and second-order ordinary differential equations of the form

$\ddot{x}+f(x,\dot{x})=0,$

where $f(x,y)$ is a polinomial function. In addition, all singularities of a
Kukles system lie on the $x$-axis.

###  1.1. Limit cycles in Kukles systems

One of the most relevant questions in planar polynomial systems and, in
particular for Kukles systems, is about the existence and number of limit
cycles, together to their stability and configuration on phase plane. A limit
cycle is an isolated periodic orbit, and it is one of the responsibles for the
organization of oscillatory patterns. The problem of finding limit cycles also
takes special relevance due to its connection with Hilbert’s 16th problem [8,
10, 11, 23, 26, 27, 33].

In order to obtain limit cycles in (1.1), it necessary that $f(x,y)$ is not
divisible by $y$ in (1.1). If $y$ divides $f(x,y)$, we can write (1.1) as

(1.3) $\left\\{\begin{array}[]{rcl}\dot{x}&=&-y\\\\[8.5359pt]
\dot{y}&=&yg(x,y),\end{array}\right.$

where $g(x,y)$ is a polynomial function with real coefficients. Here, all
points on the straight line $y=0$ are singularities. It is known that any
limit cycle of (1.3) must enclose a singularity, and since all singularities
lie on the $x$-axis, any limit cycle must cross the $x$-axis and contain a
singularity. This contradicts the uniqueness of solutions for system (1.3).
Therefore, system (1.3) does not have limit cycles. In what follows, we
consider (1.1) with $f(x,y)$ not divisible by $y$.

The number and configuration of limit cycles for Kukles systems have been
widely studied. In [34], Sadovskii proved via Hopf bifurcations that (1.2) has
at least 7 small-amplitude limit cycles that bifurcate from a weak focus. In
[15], the authors showed that the number of limit cycles in (1.2) is at least
3, while in [43] it is shown that (1.2) has at least 6 limit cycles enclosing
the origin. In [44] the authors study a reduced Kukles system with a cubic
perturbation:

$\left\\{\begin{array}[]{rcl}\dot{x}&=&-y\\\\[8.5359pt]
\dot{y}&=&x-\dfrac{a+1}{a}x^{2}+\dfrac{2-a}{1-a}y^{2}+\dfrac{1}{a}x^{3}+\epsilon
g(x,y)\end{array}\right.$

with $a>2$, $\epsilon>0$ and $g$ given by:

$g(x,y)=b_{01}y+b_{11}xy+b_{02}y^{2}+b_{12}xy^{2}+b_{21}x^{2}y+b_{03}y^{3}$

and $a_{ij},b_{ij}\in\mathbb{R}$.

The authors obtained 3 different configurations for 5 limit cycles in the
system, with 2 of them having global (large-amplitude) limit cycles. In [6],
the authors study all possible configurations of invariant straight lines in
Kukles systems, together with providing bounds for the number of limit cycles.
They also give necessary conditions for the existence of an invariant
algebraic curve of degree greater or equal to 2 that coexists with limit
cycles. Sáez and Szántó [35] found a family of Kukles systems of arbitrary
degree having an invariant non-degenerate ellipse $h(x,y)=0$:

$\left\\{\begin{array}[]{rcl}\dot{x}&=&-y\\\\[8.5359pt]
\dot{y}&=&-dc+bx+q_{n-2}(x,y)f(x,y)\end{array}\right.$

with $h(x,y)=d-2dcx+bx^{2}+y^{2},\ b>0,\ d>0$ and for $0\leq e<1,\ a>0$,

$\begin{array}[]{rcl}b&=&1-e^{2}\\\\[8.5359pt] c&=&-\dfrac{1}{a}\\\\[8.5359pt]
d&=&-e^{2}a^{2}\end{array}$

and

$q_{n-2}(x,y)=-\dfrac{1}{a}+\displaystyle\sum_{i=1}^{n-2}q_{i0}x^{i}+q_{0i}y^{i}.$

One can rewrite $h(x,y)=x^{2}+y^{2}-e^{2}(x+a)^{2}$, where $e$ corresponds to
the excentricity of the ellipse. They proved that the ellipse is an invariant
algebraic limit cycle. They also found bounds for the number of limit cycles
and showed that the coefficients of the first Melnikov function coincide with
the Lyapunov quantities. There are also studies on modified Kukles systems,
e.g. generalized Kukles systems [31] and extended Kukles systems [24], to name
a few.

In this paper, we study the existence of limit cycles for systems of the form
(1.1), and their coexistence with invariant algebraic curves. Recall that an
algebraic curve $C$ is invariant with respect to a system such as (1.1) if
there exists a polynomial $K$ such that $\dot{C}=CK$; here, $\dot{C}$ denotes
the rate of change along the orbits. We consider two leading examples: we
first consider a family of Kukles systems of degree four having an invariant
circle. We show center conditions and also prove the existence and uniqueness
of a limit cycle via Lyapunov quantities and Melnikov functions, it turns out
that the unique limit cycle is the invariant circle; the novelty here is the
use Melnikov functions of second order. Secondly, we consider a family of
Kukles systems of arbitrary odd degree with an invariant circle. We show an
upper bound for the number of limit cycles of large amplitude and provide a
concrete example when there is a connection between the small- and large-
amplitude limit cycles via a relation between the Lyapunov quantities and the
coefficients of the first Melnikov function. Here the total number of limit
cycles is not totally controled by the Melnikov function and there is an
interplay between limit cycles bifurcating from a weak focus and a center.

This paper is organized a follows. Sections 2 and 3 provide a review of the
techniques and methods used in this work. Even though they are not extensive
and far from being a complete survey of methods for studying limit cycles in
planar poltnomial systems, they intend to be a help for non-expert readers, so
that they can go through the rest of the paper. In particular, section 2
introduces Lyapunov quantities and show how they can be used for studying
limit cycles that bifurcate form a weak focus. Also, section 3 describes the
situation when the vector field is a perturbation of a Hamitonian system and
introduce Melnikov functions to study limit cycles that bifurcate from a
center; these are the main tools used throughout this paper. We provide the
corresponding references for the results exposed there. Sections 4 and 5
present the main results. In section 4 we study the coexistence of limit
cycles with an invariant circle in a Kukles family of degree four. Finally, in
section 5 we study a more general case, when limit cycles coexist with an
invariant circle in a Kukles family of arbitrary odd degree.

##  2\. Center conditions and limit cycles

By using Lyapunov functions [19, 40] one can determine the stability of a
singularity in a planar system. In this section, we will see that the good
choice of a Lyapunov function also gives information on the number of limit
cycles in planar polynomial systems and can be used to provide center
conditions.

In what follows we consider a planar polynomial vector field $X$ with an
isolated singularity at the origin. We assume that that the origin is a
center-focus of $X$, that is, the Jacobian matrix at the origin $DX(0,0)$ has
eigenvalues with zero real part. A normal form for such system is given by

(2.1) $X=\left\\{\begin{array}[]{llll}\dot{x}&=&P(x,y)=\lambda x-y+p(x,y)\\\
&\\\ \dot{y}&=&Q(x,y)=x+\lambda y+q(x,y)\end{array}\right.$

where $p,q$ are polynomials without linear terms and $\lambda\in{\mathbb{R}}$.

For the origin to be a center of (2.1), it is necessary that $\lambda=0$. If
$\lambda=0$ and the origin is not a center, we call it a _weak focus_. The
necessary conditions for a center can be obtained by calculating the _focal
quantities_ , which are polynomials in the coeficients of (2.1) that are
computed using a Lyapunov function of the form

(2.2)
$V(x,y)=\displaystyle\frac{1}{2}(x^{2}+y^{2})+V_{3}(x,y)+V_{4}(x,y)+\ldots,$

where $V_{j}(x,y)$ is a homogeneous polynomial of degree $j$ such that

$\dot{V}=\eta_{2}r^{2}+\eta_{4}r^{4}+\ldots+\eta_{2k}r^{2k}+\ldots$

with $r=x^{2}+y^{2}$.

To find $\eta_{2k},k\in\mathbb{N}$, one solves a sequence of linear equations
in the coefficients of $P$ and $Q$ as follows. If we denote by $D_{k}$ the
coefficients of degree $k$ in $\dot{V}$, the condition $D_{k}=0$ provides a
system of linear equations that allows to find the coefficients of $V_{k}$.
This linear system has a unique solution for $k$ even. In this case one can
find $\eta_{2k}$ by setting $D_{2k}=\eta_{2k}(x^{2}+y^{2})^{k}$. For more
details, see [3, 32].

The coefficients $\eta_{2k}$ are known as _focal values_. From Lyapunov
theorems [19, 40], the stability of the origin depends on the first nonzero
focal value, and the origin is a center if and only if all the focal values
are zero. Since the focal values $\eta_{2k}$ are all polynomials, Hilbert’s
Nullstellensatz [1, 22] implies that the ideal spanned by them has a finite
basis. Therefore, there exists $M>0$ such that $\eta_{2\ell}=0,$ for $\ell\leq
M$, implies $\eta_{2\ell}=0$ for all $\ell$. The value of $M$ is not known _a
priori_ , so the number of focal values one needs to calculate for finding the
basis is also unknown. Due to the extension of the calculations one needs to
perform in order to compute the focal values, the use of computational
packages becomes handy. For instance, the software _Mathematica_ [42] can be
used to calculate the first focal values, using the relations
$\eta_{2}=\eta_{4}=\dots=\eta_{2k}=0$ to eliminate some of the variables in
$\eta_{2k+2}$. Removing common factors from the reduced focal values one
proceeds until it is possible to prove that the reminding expressions do not
vanish simultaneously. The conditions over the parameters of $X$ used for
computing the focal values provide the complete necessary center conditions.
To prove that the conditions obtained are also sufficient one can use, for
instance, symmetries or integrating factors for the vector field $X$. The
reduced focal value $\eta_{2k+2}$, with the strictly positive factors removed,
is known as the _Lyapunov quantity_ $L(k)$. Note that $L(0)=\lambda$ and, with
this notation, the necessary center conditions are given by a basis for
$\\{L(k)=0,\ \forall k\geq 0\\}$. In this context, we say that the origin is a
_weak focus of order_ $k$ if $L(i)=0$, for $i=0,1,2,\ldots,k-1$ and $L(k)\neq
0$.

From Lyapunov’s theorems [19, 40] follows that a weak focus of order $k$ is
attracting if $L(k)<0$ and repelling if $L(k)>0$. In addition, if the origin
is a weak focus of order $k$, at most $k$ limit cycles can bifurcate from the
origin [3]. These limit cycles arise from Hopf bifurcations and are known as
_small-amplitude_ or _infinitesimal_ limit cycles.

##  3\. Perturbation of Hamiltonian systems and Melnikov functions

From section 2 we know that one way to obtain limit cycles in planar
polynomial systems and, in particular, in a Kukles system, is via changes of
stability of a saddle-focus equilibrium point via Hopf bifurcations; the use
of Lyapunov quantities plays a relevant role. In this section we consider
limit cycles that bifurcate from a center, via the perturbation of a
Hamiltonian system. We consider a planar vector field $X_{\epsilon}$ of the
form

(3.1)
$X_{\epsilon}=\left\\{\begin{array}[]{rcl}\dot{x}&=&\displaystyle\frac{\partial
H}{\partial
y}(x,y)+\displaystyle\sum_{k=1}^{\infty}\epsilon^{k}f_{k}(x,y,\delta)\\\\[8.5359pt]
\dot{y}&=&-\displaystyle\frac{\partial H}{\partial
x}(x,y)+\displaystyle\sum_{k=1}^{\infty}\epsilon^{k}g_{k}(x,y,\delta)\end{array}\right.,$

where $H,f_{k},g_{k}$ are analytic functions in an open set
$U\subset{\mathbb{R}}^{2},\ \epsilon$ is a small non-negative parameter and
$\delta\in D\subset{\mathbb{R}}^{n}$ is a vector of parameters in a compact
set $D$. We assume that the origin is a singularity of (3.1) for all
$\epsilon\in\mathbb{R}$ and the system has a center at the origin for
$\epsilon=0$. We call (3.1), with $\epsilon=0$ the unpertutbed system, and for
$\epsilon\neq 0$ perturbed one. The unperturbed system $X_{0}$ is Hamiltonian,
and has the form

(3.2) $X_{0}=\left\\{\begin{array}[]{rcl}\dot{x}&=&\displaystyle\frac{\partial
H}{\partial y}(x,y)\\\\[8.5359pt] \dot{y}&=&-\displaystyle\frac{\partial
H}{\partial x}(x,y)\end{array}\right.$

For $\epsilon_{0}>0$, we assume the existence of an open interval $J=(0,a)$ on
the $x$-axis, transverse to the flow of $X_{\epsilon}$ for
$|\epsilon|\leq\epsilon_{0}$. There exists a subinterval $I\subset J$ such
that the Poincaré return map

$\begin{array}[]{rcl}\pi:I\times(-\epsilon_{0},\epsilon_{0})&\mapsto&J\\\\[8.5359pt]
(x,\epsilon)&\mapsto&\pi(x,\epsilon)\end{array}$

is well defined, and maps $(x,\epsilon)$ to the $x$-coordinate
$\pi(x,\epsilon)$ of the first return of the point $(x,0)$ to $J$ by the flow
of $X_{\epsilon}$. We can define the displacement function

$\begin{array}[]{rcl}d:I\times(-\epsilon_{0},\epsilon_{0})&\mapsto&{\mathbb{R}}\\\\[8.5359pt]
(x,\epsilon)&\mapsto&d(x,\epsilon)\end{array}$

by $d(x,\epsilon)=\pi(x,\epsilon)-x$.

An orbit of $X_{\epsilon}$ through the point $(x,0)$, with $x\neq 0$ is
periodic if and only if $(x,\epsilon)$ is as zero of the displacement function
$d$. Since the origin is a center for $\epsilon=0$, we have $d(x,0)=0$ for
$x\in I$. Therefore, for $x>0$ and $\epsilon_{0}$ small enough, we have the
following series expansion for $d$:

(3.3) $d(x,\epsilon)\ =\ \displaystyle\sum_{k=1}^{\infty}M_{k}(x)\epsilon^{k}\
=\
\left.\displaystyle\sum_{k=1}^{\infty}\displaystyle\frac{1}{k!}\displaystyle\frac{\partial
d(x,\epsilon)}{\partial\epsilon^{k}}\right|_{\epsilon=0}\epsilon^{k}.$

The first nonzero function $M_{n}(x)$ in (3.3) is known as the _n-th Melnikov
function_. For $n=1$, $M_{1}$ is known as the _Melnivok integral_. Finding
simple zeroes of the first nonzero Melnikov function $M_{n}$ means finding
limit cycles that bifurcate from perturbing the Hamiltonian system (3.2). One
usually aims for finding $M_{1}$, but if $M_{1}(x)\equiv 0$ one needs to
calculate higher-order Melnikov functions. Unfortunately, there is no explicit
formula in the literature for the $n$-th Melnikov function, except for special
cases when the hamiltonian function asocciated to $X_{0}$ has a particular
form [14], and for some recursive formulas.

###  3.1. Representations for Melnikov functions and limit cycles

Suppose that system (3.2) has a family of periodic orbits, given by
$L_{h}:H(x,y)=h,\ h\in(0,\beta)$ and such that $L_{h}$ approaches the origin
when $h\to 0$. Let $h=h_{0}\in(0,\beta)$ and $A(h_{0})\in L_{h_{0}}$. Consider
the transverse section $l$ of (3.2) through $A(h_{0})$. For $h$ close to
$h_{0}$ the periodic orbit $L_{h}$ intersects $l$ in a unique point $A(h)$.
Consider the positive orbit $\gamma(h,\epsilon,\delta)$ of (3.1) with initial
point $A(h)$. Let $B(h,\epsilon,\delta)$ be the first intersection of $\gamma$
with $l$. Then,

$\begin{array}[]{rcl}H(B)-H(A)&=&\displaystyle\int_{AB}dH\\\\[8.5359pt]
&=&\epsilon[M(h,\delta)+\mathcal{O}(\epsilon)]\\\\[8.5359pt] &=&\epsilon
F(h,\epsilon,\delta)\end{array}$

with

$\begin{array}[]{rcl}M(h,\delta)&=&\displaystyle\oint_{L_{h}}(H_{y}q+H_{x}p)dt\\\\[8.5359pt]
&=&\displaystyle\oint_{L_{h}}(qdx-pdy)dt\\\\[8.5359pt]
&=&\displaystyle\iint_{H\leq h}(p_{x}+q_{y})dxdy.\end{array}$

For $\epsilon$ small enough, system (3.1) has a limit cycle near the origin if
and only if $F(h,\epsilon,\delta)$ has a positive zero for $h$ near $h=0$.
From here, $M(h,\delta)$ controls the number of limit cycles of (3.1) near the
origin. We can assume that the Hamiltonian function has the form

(3.4)
$H(x,y)=\displaystyle\frac{\omega}{2}(x^{2}+y^{2})+\displaystyle\sum_{i+j\geq
3}h_{ij}x^{i}y^{j},\ \omega>0.$

The following theorem holds [20]:

######  Theorem 3.1.

Let H be of the form (3.4). Then, $M(h,\delta)$ is $C^{\infty}$ in $0\leq h\ll
1$, with

$M(h,\delta)=h\displaystyle\sum_{s\geq 0}b_{s}(\delta)h^{s}$

If (3.2) is analytic, so is $M$. Moreover, if there exists $k\geq
1,\delta_{0}\in D$ such that $b_{k}(\delta_{0})\neq 0$ and

(3.5) $\begin{array}[]{ccc}&b_{j}(\delta_{0})=0,\qquad
j=0,1,\ldots,k-1,&\\\\[8.5359pt]
&det\left(\displaystyle\frac{\partial(b_{0},\ldots,b_{k-1})}{\partial(\delta_{1},\ldots,\delta_{k})}\right)(\delta_{0})\neq
0,&\end{array}$

where $\delta=(\delta_{1},\ldots,\delta_{m}),\ m\geq k,$ then there exist a
constant $\epsilon_{0}>0$ and a neighborhood $V$ of the origin such that for
all $0<|\epsilon|<\epsilon_{0}$ and $|\delta-\delta_{0}|<\epsilon_{0}$, (3.1)
has at most $k$ limit cycles in $V$. Moreover, for any neighborhood $V_{1}$ of
the origin there exists $(\epsilon,\delta)$ near $(0,\delta_{0})$ such that
system (3.1) has $k$ limit cycles in $V_{1}$.

This theorem provides an analytic expression for the first Melnikov function
and states that the number of limit cycles that bifurcate from a center is
related to the algebraic multiplicity of $0$ as a root of the first Melnikov
function.

Between 1996 and 2002 Giacomini, Llibre and Viano worked on a method for
computing analytically the shape of the limit cycles that bifurcate from a
center. The method is based on the following theorem [16]:

######  Theorem 3.2.

Let $X=(P,Q)$ be a $C^{1}$ vector field on an open set
$U\subset{\mathbb{R}}^{2}$. Let $V=V(x,y)$ be a $C^{1}$ solution of the
partial differential equation

(3.6) $P\displaystyle\frac{\partial V}{\partial
x}+Q\displaystyle\frac{\partial V}{\partial
y}=\left(\displaystyle\frac{\partial P}{\partial
x}+\displaystyle\frac{\partial Q}{\partial y}\right)V.$

If $\gamma$ is a limit cycle of $X$, then $\gamma\subset\left\\{(x,y)\in
U:V(x,y)=0\right\\}.$

The function $V$ is known as a _reciprocal integrating factor_ , name after
the fact that the function $R=\frac{1}{V}$ defined on $U-\\{V=0\\}$ is an
integrating factor of the vector field $X$. In [39], the authors use the
function $V$ in order to describe the limit cycles that bifurcate from a
Hamiltonian center. In [17], they look for a representation for $V$ when the
vector field $(P,Q)$ is of the form (3.1). It is proved that an analytic
solution for (3.6) has the form

$V=V(x,y,\epsilon)=\displaystyle\sum_{k=0}^{\infty}\epsilon^{k}V_{k}(x,y),$

and obtain a recursive formula for $V_{k}$ by solving a sequence of linear
partial differential equations. For $V_{0}$ the equation is

$\displaystyle\frac{\partial H}{\partial y}\frac{\partial V_{0}}{\partial
x}-\displaystyle\frac{\partial H}{\partial x}\frac{\partial V_{0}}{\partial
y}=0$

while for $V_{k}$ the equation is

$\displaystyle\frac{\partial H}{\partial y}\frac{\partial V_{k}}{\partial
x}-\displaystyle\frac{\partial H}{\partial x}\frac{\partial V_{k}}{\partial
y}=\displaystyle\sum_{j=0}^{k-1}V_{j}\left(\displaystyle\frac{\partial
f_{k-j}}{\partial x}+\displaystyle\frac{\partial g_{k-j}}{\partial
y}\right)-\displaystyle\sum_{j=0}^{k-1}\left(f_{k-j}\displaystyle\frac{\partial
V_{j}}{\partial x}+g_{k-j}\displaystyle\frac{\partial V_{j}}{\partial
y}\right).$

The solution for $V_{0}$ is used for finding $V_{1}$, and so on. One of the
main results of [17] is the following:

######  Theorem 3.3.

(Giacomini, Llibre, Viano, [17])

  1. (i)

If $T=T(h)$ is the period of the periodic orbit $H(x,y)=h$ of (3.2), then

$V_{0}(h)=\displaystyle\int_{0}^{T(h)}\left(f_{1}\displaystyle\frac{\partial
H}{\partial x}+g_{1}\displaystyle\frac{\partial H}{\partial y}\right)dt$

is the Poincaré-Melnikov integral of (3.1) associated to the orbit $H(x,y)=h$
when $V_{0}(h)\not\equiv 0$.

  2. (ii)

$V_{k}$, with $k\geq 1$ is calculated recursively as
$V_{k}(x,y)=V_{k}^{p}(x,y)+W_{k}(h)$, where

$\begin{array}[]{rcl}V_{k+1}^{p}(x,y)&=&V_{0}(h)\displaystyle\int_{0}^{t}\left(\displaystyle\frac{\partial
f_{k+1}}{\partial x}+\displaystyle\frac{\partial g_{k+1}}{\partial
y}\right)dt-
V_{0}^{\prime}(h)\displaystyle\int_{0}^{t}\left(f_{k+1}\displaystyle\frac{\partial
H}{\partial x}+g_{k+1}\displaystyle\frac{\partial H}{\partial
y}\right)dt\\\\[8.5359pt]
&&+\displaystyle\sum_{j=1}^{k}\displaystyle\int_{0}^{t}\left[V_{j}\left(\displaystyle\frac{\partial
f_{k+1-j}}{\partial x}+\displaystyle\frac{\partial g_{k+1-j}}{\partial
y}\right)-\left(f_{k+1-j}\displaystyle\frac{\partial V_{j}}{\partial
x}+g_{k+1-j}\displaystyle\frac{\partial V_{j}}{\partial
y}\right)\right]dt\end{array}$

and $V_{0}(h)W_{k}^{\prime}(h)-V_{0}^{\prime}(h)W_{k}(h)=\xi(h)$, with

$\begin{array}[]{rcl}\xi(h)&=&V_{0}(h)\varphi_{k+1}^{\prime}(h)-V_{0}^{\prime}(h)\varphi_{k+1}(h)+\displaystyle\sum_{j=1}^{k-1}\left[W_{j}(h)\varphi_{k+1-j}^{\prime}(h)-W_{j}^{\prime}(h)\varphi_{k+1-j}(h)\right]\\\\[8.5359pt]
&&+\displaystyle\int_{0}^{T(h)}\displaystyle\sum_{j=1}^{k}\left[V_{j}^{p}\left(\displaystyle\frac{\partial
f_{k+1-j}}{\partial x}+\displaystyle\frac{\partial g_{k+1-j}}{\partial
y}\right)-\left(f_{k+1-j}\displaystyle\frac{\partial V_{j}^{p}}{\partial
x}+g_{k+1-j}\displaystyle\frac{\partial V_{j}^{p}}{\partial
y}\right)\right]dt.\end{array}$

The symbol ′ denotes derivation with respect to $h$, and

$\varphi_{k}(h)=\displaystyle\int_{0}^{T(h)}\left(f_{k}\displaystyle\frac{\partial
H}{\partial x}+g_{k}\displaystyle\frac{\partial H}{\partial y}\right)dt.$

From part (i) of the theorem above one can compute the first Melnikov function
of a perturbed Hamiltonian system, while part (ii) is for computing the
following $V_{k}$, which determine the shape and position of the limit cycles
that bifurcate from a center. In [17] there are examples showing explicit
ecuations for the limit cycles that bifurcate in the Van der Pol and Liénard
systems.

Poincaré proved that the maximum number of simple zeroes of $V_{0}(h)$
coincide with the maximum number of limit cycles of system (3.1) that
bifurcate from periodic orbits $H(x,y)=h$ of (3.2) when $\epsilon$ is small.
This limit cycles are known as _global_ or _large-amplitude_ limit cycles.

When the first Melnikov function is identically null, one needs to compute
higher order Melnikov functions, which is a challenging task. There are some
results for second order Melnikov due to Iliev [25], and for higher order
Melnikov funcions via Abelian integrals [14]. Some explicit examples in wich
the vector field is written in complex coordinates can be found in [4] and
[14].

Overall, in one hand we have small-amplitude limit cycles bifurcating from an
equilibrium point, while on the other hand we have large-amplitude limit
cycles bifurcating from a Hamiltonian center. These two mechanisms for
obtainig limit cycles and the limit cycles themselves are different in
principle. However, they can be related. To our understanding, the relation
between these two type of limit cycles, conditions for which they can coincide
and what type of vector fields allow this to happen is far from being
completely understood and remains an interesting challenge beyond the scope of
this paper; we show an example of this relation in section 5.

##  4\. Coexistence of limit cycles with an invariant circle in Kukles systems

In this section we study a family of Kukles systems, in which there is an
invariant circle that coexists with an unique limit cycle. We prove this via
Lyapunov quantities and second-order Melnikov functions.

We consider the vector field family given by:

(4.1) $X_{\mu}:\left\\{\begin{array}[]{rcl}\dot{x}&=&-y\\\
\dot{y}&=&x+y(x^{2}+y^{2}-1)(ax+by+c),\end{array}\right.$

where $\mu=(a,b,c)\in{\mathbb{R}}^{3}$. We prove that the algebraic curve
$\mathcal{C}(x,y)=x^{2}+y^{2}-1$ is invariant for the system, as well as an
algebraic limit cycle. We prove that the limit cycle is unique and can be
obtained as a small-amplitude limit cycle that bifurcate from a nonhyperbolic
focus at the origin. We begin with the following:

######  Lemma 4.1.

For all $\mu\in{\mathbb{R}}^{3}$, the circle $\mathcal{C}=0$ is an algebraic
invariant curve of (4.1).

###### Proof.

We note that, for all $\mu\in{\mathbb{R}}^{3}$, the flow along the curve
satisfies

$\dot{\mathcal{C}}=\displaystyle\frac{\partial\mathcal{C}}{\partial
x}\dot{x}+\displaystyle\frac{\partial\mathcal{C}}{\partial
y}\dot{y}=\mathcal{C}(x,y)2y^{2}(ax+by+c).$

Therefore, $\mathcal{C}=0$ is an invariant circle. ∎

Note that the curve $\mathcal{C}=0$ encloses the equilibrium point at the
origin. In addition, there are conditions over the parameters of system (4.1)
for which there are no limit cycles inside the open region bounded by
$\mathcal{C}=0$. This is shown in the following result:

######  Lemma 4.2.

If $\mu=(0,0,c)$, with $c\neq 0$, system (4.1) does not have limit cycles in
the region $A=\\{(x,y)\in{\mathbb{R}}^{2}:x^{2}+y^{2}<1\\}$.

###### Proof.

Let $\tilde{\mu}=(0,0,c),\ c\neq 0$. The function

$\mathcal{D}(x,y)=\displaystyle\frac{1}{x^{2}+y^{2}-1}$

is a Dulac function of $X_{\tilde{\mu}}$ on $A$, which implies that $div(\psi
X_{\tilde{\mu}})=cte$. From Dulac’s criterion [2] there are no periodic orbits
of $X_{\tilde{\mu}}$ in $A$. ∎

In what follows we consider $a^{2}+b^{2}\neq 0$. In this case the circle
$\mathcal{C}=0$ is an invariant limit cycle:

######  Theorem 4.3.

Let $c\neq 0$. The curve $\mathcal{C}=0$ is an algebraic limit cycle of (4.1).

###### Proof.

We calculate the first Melnikov function for the vector field $X_{\mu}$.
Rescaling parameters by

$a\to\epsilon a;\ b\to\epsilon b;\ c\to\epsilon c,$

and reversing time $t\to-t$, system (4.1) can be written as

(4.2) $\begin{array}[]{rcl}\dot{x}&=&y\\\ \dot{y}&=&-x+\epsilon
y(1-x^{2}-y^{2})(ax+by+c).\end{array}$

System (4.2) is a perturbation of the Hamiltonian system

$\begin{array}[]{rcl}\dot{x}&=&y\\\ \dot{y}&=&-x,\end{array}$

with Hamiltonian associated function
$H(x,y)=\displaystyle\frac{x^{2}+y^{2}}{2}$. We parameterize the orbit of the
isocronus center $H=h$ by

$\begin{array}[]{rcl}x(t)&=&\sqrt{2h}\cos t\\\ y(t)&=&-\sqrt{2h}\sin
t\end{array};\ t\in[0,2\pi].$

From theorem 3.3, the first Melnikov function has the form

$M_{1}(h)=2h(2h-1)c\pi.$

Then, $M_{1}$ has as simple zero for $h=\displaystyle\frac{1}{2}$, which
corresponds to the curve $\mathcal{C}=0$. Therefore, the curve $\mathcal{C}=0$
in an algebraic limit cycle of (4.1) ∎

The limit cycle obtained above can bifurcate from the origin, as we can see in
the following theorems:

######  Theorem 4.4.

Let $a,b\in{\mathbb{R}}-\\{0\\}$. If $c=0$, system (4.1) has a weak focus of
order $1$ at the origin. This limit cycle is repelling for $ab<0$ and
attracting for $ab>0$.

###### Proof.

The linear part of (4.1) is given by:

$DX_{\mu}(0,0)\;=\;\left(\begin{array}[]{cc}0&-1\\\ 1&-c\\\
\end{array}\right).$

Note that the origin is a weak focus if $c=0$. To calculate its order, we
compute the Lyapunov quantities $L(k),\ k\geq 0$ in order. We have $L(0)=c=0$,
and

$L(1)=-\displaystyle\frac{ab}{8},$

which proves the theorem. ∎

######  Theorem 4.5.

In parameter space ${\mathbb{R}}^{3}$, there exists an open set $\mathcal{N}$
such that, for all $\mu\in\mathcal{N}$, system (4.1) has at least one small-
amplitude limit cycle.

###### Proof.

From theorem 4.4, system (4.1) has a weak focus of order 1 at the origin for
$\tilde{\mu}=(a,b,0)$, with $ab\neq 0$. If $ab<0$ the weak focus is repelling,
then we can perturb the system in such a way that $c<0$, changing the
stability of the origin and thus creating a new repelling limit cycle via a
Hopf bifurcation. Similarly, if $ab>0$, the weak focus is attracting and we
can perturb the system with $c>0$, creating an attracting limit cycle. ∎

We can also show center conditions in system (4.1):

######  Theorem 4.6.

System (4.1) has a center at the origin if and only if $c=a=0$ or $c=b=0$.

###### Proof.

From theorems 4.4 and 4.5, for $c=0$ the Lyapunov quatities are $\ L(0)=0$ and
$L(1)=-\displaystyle\frac{ab}{8}$. If $a=0$ or $b=0,\ L(1)=L(2)=L(k)=0,\
\forall k\geq 0$ and we have necessary center conditions. On the other hand,
for $c=a=0;\ b\neq 0$, system (4.1) has the form:

$\begin{array}[]{rcl}\dot{x}&=&p(x,y)\ =\ -y,\\\ \dot{y}&=&q(x,y)\ =\
x+b(-y^{2}+x^{2}y^{2}+y^{4}).\end{array}$

In this case we have the symmetries $p(x,-y)=-p(x,y)$ and $q(x,-y)=q(x,y)$,
which prove that $X_{\mu}$ has a center at the origin if $c=a=0$.
Analogously, if $c=b=0;\ a\neq 0$, the vector field has the form

$\begin{array}[]{rcl}\dot{x}&=&p(x,y)\ =\ -y,\\\ \dot{y}&=&q(x,y)\ =\
x+a(-xy+x^{3}y+xy^{3}),\end{array}$

and the symmetries $p(-x,y)=p(x,y);\ q(-x,y)=-q(x,y)$. Therefore, the origin
is a center of $X_{\mu}$ for $c=b=0$.
This proves that the conditions are also sufficient. ∎

When the parameter $c$ vanishes in (4.1), the information obtain from the
first Melnikov function is lost, since $M_{1}\equiv 0$ for $c=0$. However,
following the method from [14] and [25], it is possible to compute the second
Melnikov function. This is shown in the following result:

######  Theorem 4.7.

If $ab\neq 0$, system (4.1) has exactly one limit cycle.

###### Proof.

From theorem 4.4, for $c=0$ there exists a limit cycle inside the region
bounded by $\mathcal{C}=0$. For small values of $h$, the second Melnikov
function can be represented by:

$M_{2}(h)=\displaystyle\int_{H=h}\left[G_{1h}(x,y)P_{2}(x,h)-G_{1}(x,y)P_{2h}(x,y)\right]dx,$

where

$\begin{array}[]{rcl}G_{1h}(x,y)&=&\displaystyle\frac{ax-
ax^{3}-3axy^{2}}{y}\\\\[5.69046pt]
P_{2}(x,h)&=&b2hx-4bh^{2}x-b\displaystyle\frac{x^{3}}{3}+\displaystyle\frac{2}{3}bhx^{3}\\\\[5.69046pt]
G_{1}(x,y)&=&-axy-ax^{3}y-axy^{3}\\\\[8.5359pt]
P_{2h}(x,y)&=&\displaystyle\frac{6bx+2bx^{3}-24bxy}{3y}\end{array}$

then $M_{2}(h)=\displaystyle\frac{1}{3}abh(2h-1)(6-7h+28h^{2})\pi$. The closed
orbit associated to $h=\displaystyle\frac{1}{2}$ is still a limit cycle.
Furthermore, the quadratic polynomial $6-7h+28h^{2}$ does ot have real roots;
this holds for all $a,b\in{\mathbb{R}}-\\{0\\}$. Then, the limit cycle is
unique. ∎

######  Remark 4.8.

It is expected that the Melnikov function up to some order provides the total
number of limit cycles, but there are cases when the Melnikov function does
not control the limit cycles of small amplitude. It is worth asking about the
relation between small- and large-amplitude limit cycles, or equivalently
about the relation (if exists) between the coefficients of the Melnikov
functions and the Lyapunov quantities. There is no general answer to this
question, and we think an answer may be of help on tackling Hilbert’s 16th
problem. In [21] the authors give a partial answer to this question, about
which we mention in next section.

##  5\. Coexistence of limit cycles with an invariant circle in a family of
Kukles systems with arbitrary odd degree

In this section we study the coexistence of limit cycles with an invariant
circle in a family of Kukles systems of arbitrary odd degree. We study
Lyapunov quantities and the first Melnikov function, and show a new concrete
example, different from the one given in [20], in which theorem 3.1 does not
apply, i.e. it is not possible to stablish the cyclicity of the system with
only the coefficients of the first Melnikov function.

We consider the following family of Kukles systems:

(5.1) $X=\left\\{\begin{array}[]{rcl}\dot{x}&=&-y,\\\
\dot{y}&=&x+y(1-x^{2}-y^{2})(b_{00}+\displaystyle\sum_{i+j=1}^{n}b_{2i2j}x^{2i}y^{2j}),\end{array}\right.$

where $b_{2i2j}\in{\mathbb{R}}$ and $i,j$ are non-negative integer numbers.
This family of planar vector fields of degree $2n+3$ has the invariant circle
$C(x,y)=0$, with

$C(x,y)=1-x^{2}-y^{2},$

as we show in the following lemma:

######  Lemma 5.1.

The circle $C=0$ is an algebraic invariant curce of (5.1).

###### Proof.

A straightforward calculation shows that

$\displaystyle\frac{\partial C}{\partial x}\dot{x}+\displaystyle\frac{\partial
C}{\partial
y}\dot{y}=C(x,y)(-2y^{2})(b_{00}+\displaystyle\sum_{i+j=1}^{n}b_{2i\;2j}x^{2i}y^{2j}),$

which implies that $C$ is invariant. ∎

The origin is the only singularity of system (5.1), and it is a _center-focus_
for $b_{00}=0$. This can be easily checked by looking at the linearization of
the system, given by

$DX(0,0)\;=\;\left(\begin{array}[]{cc}0&-1\\\ 1&b_{00}\\\ \end{array}\right)$

We study the coexistence of limit cycles of (5.1) with the invariant curve
$C=0$. To this end, we calculate the first Melnikov function. We have the
following theorem:

######  Theorem 5.2.

For all $n\in{\mathbb{N}}$, system (5.1) has at most $n+1$ limit cycles of
large amplitude.

###### Proof.

We write system (5.1) as the perturbation of a Hamiltonian system. The
parameter rescaling $b_{2i2j}\to\epsilon\cdot b_{2i,2j}$, for $\epsilon>0$,
gives the vector field

(5.2) $X_{\epsilon}=\left\\{\begin{array}[]{rcl}\dot{x}&=&-y,\\\
\dot{y}&=&x+\epsilon
y(1-x^{2}-y^{2})(b_{00}+\displaystyle\sum_{i+j=1}^{n}b_{2i\;2j}x^{2i}y^{2j}).\end{array}\right.$

Note that the unperturbed system $(\epsilon=0)$ is Hamiltonian:

(5.3) $X_{0}=\left\\{\begin{array}[]{rcl}\dot{x}&=&-y,\\\
\dot{y}&=&x,\end{array}\right.$

with associated Hamiltonian function
$H(x,y)=\displaystyle\frac{x^{2}+y^{2}}{2}$. Reverting the direction of the
flow via the time rescaling $t\to-t$, we finally write the system as

$X_{\epsilon}=\left\\{\begin{array}[]{rcl}\dot{x}&=&\displaystyle\frac{\partial
H}{\partial y},\\\\[11.38092pt] \dot{y}&=&-\displaystyle\frac{\partial
H}{\partial x}+\epsilon g_{1}(x,y),\end{array}\right.$

where

$g_{1}(x,y)=y(x^{2}+y^{2}-1)(b_{00}+\displaystyle\sum_{i+j=0}^{n}b_{2i\;2j}x^{2i}y^{2j}).$

We now consider the period $T=2\pi$ of the orbit $H(x,y)=h$ of the unperturbed
system. We parameterize the orbit by

$\begin{array}[]{rcl}x(t)&=&\sqrt{2h}\cos(t),\\\\[8.5359pt]
y(t)&=&-\sqrt{2h}\sin(t),\end{array}$

with $h>0$. From theorem 3.3, the first Melnikov function has the form

$\begin{array}[]{rcl}M_{1}(h)&=&\displaystyle\int_{0}^{T}y(t)g_{1}(x(t),y(t))dt\\\\[14.22636pt]
&=&2h(2h-1)\left(\pi
b_{00}+\displaystyle\sum_{i+j=1}^{n}b_{2i,2j}(2h)^{i+j}\displaystyle\int_{0}^{2\pi}\cos^{2i}(t)\sin^{2j+2}(t)dt\right).\end{array}$

Note that

$\displaystyle\int_{0}^{2\pi}\cos^{2i}(t)\sin^{2j+2}(t)dt=4\displaystyle\int_{0}^{\frac{\pi}{2}}\cos^{2i}(t)\sin^{2j+2}(t)dt,$

and for $r,s>0$

$\displaystyle\int_{0}^{\frac{\pi}{2}}\cos^{2r-1}(t)\sin^{2s-1}(t)dt=\displaystyle\frac{\Gamma(r)\Gamma(s)}{2\Gamma(r+s)},$

where $\Gamma(\cdot)$ corresponds to the _Gamma_ function. Then

$\begin{array}[]{rcl}M_{1}(h)&=&h(2h-1)\left(2b_{00}+\displaystyle\sum_{k=1}^{n}\sum_{i+j=k}\displaystyle\frac{(2i)!(2j)!(2j+1)}{2^{i+j-1}i!j!(i+j+1)!}b_{2i2j}h^{i+j}\right)\pi.\end{array}$

This polynomial has degree $n+2$.
Now note that $h=\frac{1}{2}$ is a simple zero of $V_{0}(h)$, i.e. the orbit
$H(x,y)=\frac{1}{2}$ is a limit cycle or equivalently the curve $C=0$ is a
limit cycle. The simple zero $h=0$ does not provide any limit cycle since the
orbit $H=0$ reduces to the origin. Therefore, system (5.1) has at most $n+1$
large-amplitude limit cycles. ∎

######  Example 5.3.

System (5.1) with $n=3$ gives:

(5.4) $X=\left\\{\begin{array}[]{rcl}\dot{x}&=&-y,\\\\[2.84544pt]
\dot{y}&=&x+y(1-x^{2}-y^{2})(b_{00}+b_{20}x^{2}+b_{02}y^{2}+b_{40}x^{4}+b_{22}x^{2}y^{2}+\\\\[2.84544pt]
&&b_{04}y^{4}+b_{60}x^{6}+b_{42}x^{4}y^{2}+b_{24}x^{2}y^{4}+b_{06}y^{6}),\end{array}\right.$

which has degree 9. Using theorem 5.2, the first Melnikov function has the
form

$\begin{array}[]{rcl}M_{1}(h)&=&\displaystyle\frac{\pi}{4}h(2h-1)(8b_{00}+4(3b_{02}+b_{20})h+4(5b_{04}+b_{22}+b_{40})h^{2}+\\\\[8.5359pt]
&&(35b_{06}+5b_{24}+3b_{42}+5b_{60})h^{3}).\end{array}$

From here, at most 4 limit cycles can bifurcate from the Hamiltonian center,
including the invariant circle $C=0$.

Using the software _Mathematica_ we calculate the Lyapunov quantities $L(k)$,
$k\geq 0$. We have:

$L(0)=b_{00},$

and if $b_{00}=0,\ L(0)=0$. Then

$L(1)=\displaystyle\frac{3b_{02}+b_{20}}{8}.$

If $b_{20}=-3b_{02},\ L(1)=0$, then we have

$L(2)=\displaystyle\frac{5b_{04}+b_{22}+b_{40}}{16},$

and if $b_{22}=-5b_{04}-b_{40},\ L(2)=0$, then

$L(3)=\displaystyle\frac{-6b_{02}^{3}+35b_{06}+5b_{24}+3b_{42}+5b_{60}}{128}.$

We assume the condition

(5.5) $5b_{24}+3b_{42}+5b_{60}\neq 0,$

so that the higher order term in $M_{1}$ is not zero. If we set
$b_{42}=\displaystyle\frac{6b_{02}^{3}-35b_{06}-5b_{24}-5b_{60}}{3}$,
$L(3)=0$, we have

$L(4)=\displaystyle\frac{3}{128}b_{02}^{2}(2b_{02}-b_{04}+b_{40}).$

In addition, if $b_{04}=2b_{02}+b_{40}$ we have that $L(4)=0$. It follows that

$L(5)=\displaystyle\frac{b_{02}(144b_{02}^{2}+9b_{02}b_{24}+204b_{02}b_{40}+80b_{40}^{2}+21b_{02}b_{60})}{1536},$

Now, if $b_{40}=0$, we can write

$\ L(5)=\displaystyle\frac{3}{1536}b_{02}^{2}(48b_{02}+3b_{24}+7b_{60}).$

From here, setting $b_{24}=-\displaystyle\frac{1}{3}(48b_{02}+7b_{60})$, we
have $L(5)=0$ and then

$L(6)=-2b_{02}^{2}(246b_{02}+9b_{02}^{3}-183b_{06}-5b_{60}).$

For $b_{60}=\displaystyle\frac{3}{5}(82b_{02}+3b_{02}^{2}-61b_{06})$, it
follows that $L(6)=0$ and we can write

$L(7)=\displaystyle\frac{18}{25}b_{02}(749152b_{02}^{2}+77586b_{02}^{4}+1737b_{02}^{6}-785952b_{02}b_{06}-45918b_{02}^{3}b_{06}+171488b_{06}^{2})$

For this set of conditions upon the parameters, condition (5.5) is equivalent
to $b_{02}\neq 0$. If $b_{02}=0$, system (5.4) has a center at the origin.
Let’s consider $b_{02}\neq 0$ and define functions $\xi$ and $\zeta$ defined
over an open interval in $\mathbb{R}$ by

$\begin{array}[]{rcl}\xi(b_{02})&=&392976b_{02}+22959b_{02}^{3}-5\sqrt{1038382336b_{02}^{2}+819584160b_{02}^{4}+916941b_{02}^{6}},\\\\[8.5359pt]
\zeta(b_{02})&=&392976b_{02}+22959b_{02}^{3}+5\sqrt{1038382336b_{02}^{2}+819584160b_{02}^{4}+916941b_{02}^{6}}.\end{array}$

We can now write the following Lyapunov quantities by

$\begin{array}[]{rcl}L(7)&=&b_{02}(b_{06}-\xi)(b_{06}-\zeta),\\\\[8.5359pt]
L(8)&=&F(b_{02},b_{06}),\end{array}$

where $F:A\subset\mathbb{R}^{2}\to\mathbb{R}$ and $A$ is an open set. On the
curve $b_{06}=\xi(b_{02})$, we have that $L(7)=0$. The implicit function
theorem implies the existence of a point $P=(b_{02}^{\ast},b_{06}^{\ast})$ on
the graph of $\xi$ on the $(b_{02},b_{06})$-parameter plane, and a
neighborhood $\mathcal{U}$ of $P$ where it is possible to define a real
function $\gamma$ satisfying $b_{06}=\gamma(b_{02})$, with
$\gamma^{\prime}>0,\;\gamma^{\prime\prime}<0$ and with
$F(b_{02},\gamma(b_{02}))=0$, such that $L(9)>0$ in $\mathcal{U}$. Similarly,
on the curve $b_{06}=\zeta(b_{02})$, $L(7)$ also vanishes and in a similar way
we can use the implicit function theorem to find a point
$Q(\tilde{b}_{02},\tilde{b}_{06})$ on the graph of $\zeta$ and a neighborhood
$\mathcal{V}$ of $Q$ when it is possible to define a curve
$b_{06}=\phi(b_{02})$ with $\phi^{\prime}>0,\;\phi^{\prime\prime}>0$ y and
with $F(b_{02},\phi(b_{02}))=0$, such that $L(9)<0$ in $\mathcal{V}$; we omit
the exact expressions of $L(8)$ and $L(9)$ due to their extension. From here,
at least 9 limit cycles can bifurcate from the origin.

Note that the Lyapunov quantities $L(0),L(1)$ and $L(2)$ appear in the
coefficients of the first Melnikvo funciton $M_{1}(h)$. This suggests a
relation between the small- and large-amplitude limit cycles of system (5.4).
In fact, this also suggests that the limit cycles obtained perturbing the
Lyapunov quatities are related with the ones guaranteed by theorem 3.1, but
there are more limit cycles that are not controled by the first Melnikov
function. In this case the condition (3.5) over the determinant is not
satisfied. Example 5.3 provides a new system in which the number of limit
cycles that bifurcate from a weak focus at the origin is larger than the
number of limit cycles that can be obtained as simple zeros of the first
Melnikov function. The fact that the Lyapunov quantities appear in the
coefficients of the first Melnikov function gives a hint of that this relation
may imply that some of the limit cycles obtained by the two different methods
(perturbation from an equilibrium point vs perturbation from a family of
periodic orbits) can coincide. The actual nature of this relation and a
complete caracterization remains an interesting challenge for future research.

##  6\. Aknowledgements

The author would like to express his gratitude to Iván Szántó for his
guidance, support and patience.

## References

  * [1] J. M. Almira, Nullstellensatz revisited, Rend. Sem. Mat. Univ. Pol. Torino - Vol. 65 (3) (2007), 365–369.
  * [2] A. Andronov, Qualitative Theory of Second-Order Dynamic Systems, New York:Wiley, 1973.
  * [3] T. R. Blows & N. G. Lloyd, The number of limit cycles of certain polynomial differential equations, Proc. Roy. Soc. Edinburg Sect. A 98 (1984), 215-239
  * [4] A. Buica, A. Gasull & J. Yang, The third order Melnikov function of a quadratic center under quadratic perturbations, J. Math. Anal. Appl. 331 (2007) 443-454.
  * [5] T. Carletti & G. Villari, A note on the existence and uniqueness of limit cycles for Liénard systems, J. Math. Anal. Appl., 307, (2005), 763-773.
  * [6] J. Chavarriga, E. Sáez, I. Szantó & M. Grau, Coexistence of limit cycles and invariant algebraic curves for a Kukles system, Nonlinear Analysis 59 (2004), 673-693.
  * [7] C.J. Christopher & N.G. Lloyd, On the paper of Jin and Wang concerning the conditions for a centre in certain cubic systems, Bull. London Math. Soc. 22 (1990), 5-12.
  * [8] C.J. Christopher & N.G. Lloyd, Polynomial systems: a lower bound for the Hilbert numbers, T Proc. Roy. Soc. London Ser. A450, no. 1938 (1995), 219–224.
  * [9] C.J. Christopher,  Invariant algebraic curves and conditions for a centre, Proceedings of the Royal Society of Edinburgh, 124A, 1209-1229, 1994.
  * [10] F. Dumortier, M. ElMorsalami & C.Rousseau, Hilbert’s 16th problem for quadratic systems and cyclicity of elementary graphics, Nonlinearity 9 (1996), 1209–1261.
  * [11] F. Dumortier, R. Roussarie & C. Rousseau, Hilbert’s 16th problem for quadratic vector fields, Differential Equations 110, no. 1 (1994), 86–133.
  * [12] R. Field & R. Noyes, Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction, Journal of Chemical Physics (1974), 60(March), 1877–1884.
  * [13] R. FitzHugh, Impulses and Physiological States in Theoretical Models of Nerve Membrane, Biophysical Journal (1961), 1(6), 445–466.
  * [14] J.P. Francoise, Successive derivatives of a first return map, application to the study of quadratic vector fields, Ergod. Theory Dyn. Syst. 16 (1996) 87-96.
  * [15] V.A. Gaiko & W.T. van Horssen, Global bifurcation of limit and separatrix cycles in a generelized Lienard systems, Nonlinear Anal. 59 (2004), 189-198.
  * [16] H.Giacomini, J. Llibre & M. Viano, On the nonexistence, existence and uniqueness of limit cycles, Nonlinearity 9 (1996) 501-516.
  * [17] H.Giacomini, J. Llibre & M. Viano, On the shape of limit cycles that bifurcate from Hamiltonian centers, Nonlinear Analysis 41 (2000) 523-537.
  * [18] N. Glade, J. Demongeot, J. Tabony, Microtubule self-organisation by reaction-diffusion processes causes collective transport and organisation of cellular particles, BMC Cell Biology (2004), 523.
  * [19] W. Hahn, Stability of Motion, Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen, Springer-Verlag, 1967.
  * [20] M. Han,  Hopf Bifurcations for near-Hamiltonian Systems, International Journal of Bifurcation and Chaos, Vol. 19, No. 12 (2009) 4117-4130.
  * [21] M.Han, On Hopf Cyclicity of Planar Systems, J. of Math. Anal. and Appl. 245, 404-422 (2000).
  * [22] D. Hilbert, Ueber die Theorie der algebraischen Formen, Mathematische Annalen, 36 (4), 473–534.
  * [23] D. Hilbert,  Mathematical Problems, Bull. Amer. Math. Soc. 8, 437-479.
  * [24] J.M Hill, N. G. Lloyd & J. Pearson, Centres and limit cycles for an extended Kukles system, Electronic Journal of Differential Equations, Vol. 2007(2007), No. 119, pp. 1-23.
  * [25] I.D.Iliev, On second order bifurcations of limit cycles, J. London Math. Soc. 58 (1998) 353-366.
  * [26] Yu. Ilyashenko, Centennial History of Hilbert’s 16th problem, Bulletin (New Series) of the American Mathematical Society, Volume 39, Number 3, 301–354.
  * [27] Yu. Ilyashenko & S. Yakovenko, Concerning the Hilbert 16th problem, American Mathematical Society, Providence, RI, 1995.
  * [28] X. Jing & D. Wang,  On the conditions of Kukles for the existence of a centre, Bull. London Math. Soc. 22 (1990), 1-4.
  * [29] I.S. Kukles,  Sur quelques cas de distinction entre un foyer et un centre, Dokl. Akad. Nauk. SSSR 42 (1944), 208 - 211.
  * [30] J. Li, S. Inada, H. Dobrzynski & H. Zhang, A modified FitzHugh-Nagumo model that allows control of action potential duration and refractory period, Computing in Cardiology (2009), 65–68.
  * [31] J. Llibre & A. C. Mereu, Limit cycles for generalized Kukles polynomial differential systems, Nonlinear Annalysis 74 (2011), 1261-1271.
  * [32] N.G. Lloyd, Limit cycles of polynomial systems - some recent developments, In New Direction in Dynamical Systems. Cambridge University Press, 1988, pp. 192-234.
  * [33] R. Roussarie, Bifurcations of Planar Vector Fields and Hilbert’s Sixteenth Problem, Birkhauser, Basel, 1998.
  * [34] A. P. Sadovskii, Cubic systems of nonlinear oscillations with seven limi cycles, Fiffer. Uravne. 39 (2003), 472-481 (in Russian).
  * [35] E. Sáez & I. Szántó, Bifurcations of limit cycles in Kukles systems of arbitrary degree with invariant ellipse, Appl. Math. Lett. (2012), doi: 10.1016/j.aml.2012.01.039 (in press).
  * [36] R. Thomas, On the Relation Between the Logical Structure of Systems and Their Ability to Generate Multiple Steady States or Sustained Oscillations, in: J. Della Dora, J. Demongeot, B. Lacolle (Eds.), Numer. Methods Study Crit. Phenom., Springer Berlin Heidelberg (1981), 180–193.
  * [37] B. Van der Pol, Biological Rhythms Considered as Relaxation Oscillations, Acta Medica Scandinavica (1940), 103(S108), 76–88.
  * [38] B. Van der Pol & J. Van der Mark, The heartbeat considered as a relaxation oscillation, and an electrical model of the heart, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science (1928), 6(38), 763–775.
  * [39] M. Viano, J. Llibre & H.Giacomini, Arbitrary order bifurcations for perturbed Hamiltonian planar systems via the reciprocal of an integrating factor, Nonlinear Analysis 48 (2002) 117-136.
  * [40] S.W. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, Ney York, 1990.
  * [41] H.R Wilson & J.D. Cowan, Excitatory and inhibitory interactions in localized populations of model neurons, Biophysical Journal (1972), 12(1), 1–24.
  * [42] Wolfram Research Mathematica , A System for Doing Mathematics by Computer, Champaign, IL, (1988).
  * [43] W. Ye & Y. Ye, Qualitative theory of the Kukles system, I. Number of critical points, Ann. Differ. Equations 17 (2001), 275-286.
  * [44] H. Zhang, T. Zhang, Y. Tian & M. Tadé , Limit cycles for the Kukles system, Journal of Dynamical and Control Systems, Vol. 14, No. 2 (April 2008), 283-298.