\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 89, pp. 1--17.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2017/89\hfil Second-order bifurcation of limit cycles]
{Second-order bifurcation of limit cycles from a quadratic reversible center}

\author[L. Peng, B. Huang \hfil EJDE-2017/89\hfilneg]
{Linping Peng, Bo Huang}

\address{Linping Peng \newline
School of Mathematics and System Sciences,
Beihang University,
LIMB of the Ministry of Education,
 Beijing 100191, China}
\email{penglp@buaa.edu.cn}

\address{Bo Huang \newline
School of Mathematics and System Sciences,
Beihang University,
LIMB of the Ministry of Education,
 Beijing, 100191, China}
\email{bohuang0407@buaa.edu.cn}

\dedicatory{Communicated by Zhaosheng Feng}

\thanks{Submitted March 2, 2016. Published March 28, 2017.}
\subjclass[2010]{34C07, 37G15, 34C05}
\keywords{Hamiltonian system; bifurcation; limit cycles; perturbation;
\hfill\break\indent  averaging method; quadratic center}

\begin{abstract}
 This article concerns the bifurcation of limit cycles from a quadratic
 integrable and non-Hamiltonian system. By using the averaging theory,
 we show that under any small quadratic homogeneous perturbation,
 there is at most one limit cycle  for the first order bifurcation
 and two for the second-order bifurcation arising from the period annulus
 of the unperturbed system, respectively. Moreover, in each case the upper
 bound is sharp.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\section{Introduction} \label{sect1}

In the qualitative theory of real planar differential systems, one
of the important problems is to determine the number of limit
cycles. To solve this problem, some innovative methods have been
proposed based on the Poincar\'e map \cite{BP, CJ, LLZ}, the
Poincar\'e-Pontryagin-Melnikov integrals or the Abelian integrals
\cite{AI, AN, CGP, XH}, the inverse integrating factor \cite{GLV1,
GLV2, GLV3, VLG}, and the averaging method \cite{ BL1, CLP, GL, L, LL,
LPR} which is actually equivalent to the Abelian integrals in the
plane.

The averaging method serves as one of powerful tools for studying limit
cycles, which can reduce the problem regarding the number of limit
cycles of some differential systems to the exploration of the number
of hyperbolic equilibrium points of their averaged differential
equations. By using the averaging theory, some elegant results on
the number of limit cycles of the differential systems have been
obtained, such as by Buic\u a and Llibre \cite{BL3}, by Gine and
Llibre \cite{GL}, by Li and Llibre \cite{LL} and so on.

In this article, we start with the  quadratic system
\begin{equation}\label{eq1}
\begin{gathered}
\dot{x}=-y+xy,\\
\dot{y}=x+y^2,
\end{gathered}
\end{equation}
and investigate the second-order bifurcation of limit cycles under
any small quadratic homogeneous perturbations. Obviously, system
\eqref{eq1} has
$$
H(x,y)=\frac{1-x}{\sqrt{x^2+y^2}}=c
$$
as its first integral with the integrating factor $1/(x^2+y^2)^{3/2}$,
and has the unique finite singularity $(0, 0)$ as its isochronous center.
The period annulus, denoted by
$$
\{(x,y)|H(x,y)=c,\quad c\in(1, +\infty)\},
$$
starts at the center $(0, 0)$ and
terminates at the separatrix passing the infinite degenerate
singularity on the equator. The phase portrait of system \eqref{eq1}
is shown in Figure \ref{fig1}.

\begin{figure}[htp]
\begin{center}
 \includegraphics[width=0.4\textwidth]{fig1}
\end{center}
  \caption{Phase portrait of system \eqref{eq1} in the Poincar\'{e} disk }
\label{fig1}
 \end{figure}


By using the averaging theory, we study
the bifurcation of limit cycles for system \eqref{eq1} under any
small quadratic homogeneous perturbations. Our main result is as follows.

\begin{theorem}\label{thm1}
For any sufficiently small parameter $|\varepsilon|$, and any real
constants $a_{ij}^{(k)}$ and $b_{ij}^{(k)}$  $(i,j=0,1,2,k=1,2)$,
considering the quadratic homogeneous perturbed system
\begin{equation}\label{eq2}
\begin{gathered}
\dot{x}=-y+xy+\sum_{k=1}^2\varepsilon^{k}\sum_{i+j=2}a_{ij}^{(k)}x^{i}y^{j},\\
\dot{y}=x+y^2+\sum_{k=1}^2\varepsilon^{k}\sum_{i+j=2}b_{ij}^{(k)}x^{i}y^{j},
\end{gathered}
\end{equation}
we have
\begin{enumerate}
\item By using the averaging theory of first order, system \eqref{eq2}
has at most one limit cycle bifurcating from
the periodic orbits of the unperturbed one, and this upper bound is
sharp.
\item By using the averaging theory of second order, system \eqref{eq2}
has at most two limit cycles bifurcating
 from the periodic orbits of the unperturbed one, and this upper bound is sharp.
\end{enumerate}
\end{theorem}



The rest of this article is organized as follows. In Section $2$, we
give an introduction on the averaging theory of first and second
order, including some technical lemmas and methods employed in the
averaging theory. Sections $3$ and $4$ are dedicated to the study of
the bifurcation of limit cycles by computing the first and second
order averaged functions related to  system
\eqref{eq2} and exploring the number of the simple zeros.
In addition, some examples are given to illustrate the
established results.



\section{Preliminary results} \label{sect2}


In this section, we briefly introduce the averaging theory of first
and second order, and some technical lemmas which will be used in
the proof of our main results.


\begin{lemma}[\cite{BL1}] \label{lem1}
Consider the differential system
\begin{equation}\label{eq3}
\dot{x}(t)=\varepsilon F_1(t, x)+\varepsilon^2F_2(t,
x)+\varepsilon^3F_3(t, x)+\varepsilon^{4}W(t, x, \varepsilon),
\end{equation}
where $F_1, F_2, F_3: \mathbb{R}\times D\to \mathbb{R}$,
$W: \mathbb{R}\times D\times (-\varepsilon_0,
\varepsilon_0)\to \mathbb{R}$ $(\varepsilon_0>0)$ are
continuous functions and $T$-periodic in the first variable, and
$D$ is an open subset of $\mathbb{R}$. Assume that the following two
hypotheses hold:

(i) $F_1(t, \cdot)\in C^2(D)$, $F_2(t, \cdot)\in C^{1}(D)$ for
all $t\in\mathbb{R}, F_1, F_2, F_3, W, D_{x}^2F_1, D_{x}
F_2$ are locally Lipschitz with respect to $x$, and $W$ is twice
differentiable with respect to $\varepsilon$.

Define $F_k^0: D\to\mathbb{R}$ for $k=1, 2, 3$ as
\begin{gather*}
F_1^0(x)=\frac{1}{T}\int_0^{T}F_1(s, x)ds,\\
F_2^0(x)=\frac{1}{T}\int_0^{T}\Big[\frac{\partial F_1(s,
x)}{\partial x}y_1(s, x)+F_2(s, x)\Big]ds,\\
\begin{aligned}
F_3^0(x)&=\frac{1}{T}\int_0^{T}\Big[\frac{1}{2}\frac{\partial^2
F_1(s, x)}{\partial x^2}y_1^2(s,x) 
 +\frac{1}{2}\frac{\partial F_1(s, x)}{\partial x}y_2(s,x)\\
&\quad +\frac{\partial F_2(s, x)}{\partial x}y_1(s, x)+F_3(s,x)\Big]ds,
\end{aligned}
\end{gather*}
where
\begin{gather*}
y_1(s, x)=\int_0^{s}F_1(t, x)dt,\\
y_2(s, x)=2\int_0^{s}\left[\frac{\partial F_1(t, x)}{\partial
x}y_1(t, x)+F_2(t, x)\right]dt.
\end{gather*}

(ii) For an open and bounded set $V\subset D$ and for each
$\varepsilon\in (-\varepsilon_0, \varepsilon_0)\backslash
\{0\}$, there exists $a\in V$ such that
$(F_1^0+\varepsilon F_2^0+\varepsilon^2 F_3^0)(a)=0$ and
$$
\frac{d}{dx} (F_1^0+\varepsilon F_2^0+\varepsilon^2
F_3^0)(a) \neq 0.
$$

Then for sufficiently small $|\varepsilon|>0$, there exists a
$T$-periodic solution $x(t, \varepsilon)$ of system \eqref{eq3} such
that $x(0, \varepsilon)\to a$ as $\varepsilon\to 0$.
\end{lemma}

\begin{corollary}\cite{BL1}\label{coro1}
Under the hypotheses of Lemma \ref{lem1}, if $F_1^0(x)$ is not
identically zero, then the zeros of $(F_1^0+\varepsilon
F_2^0+\varepsilon^2 F_3^0)(x)$ are mainly the zeros of
$F_1^0(x)$ for sufficiently small $|\varepsilon|$. In this case,
conclusions in Lemma \ref{lem1} are true.

If $F_1^0(x)$ is identically zero and $F_2^0(x)$ is not
identically zero, then the zeros of $(F_1^0+\varepsilon
F_2^0+\varepsilon^2 F_3^0)(x)$ are mainly the zeros of
$F_2^0(x)$ for sufficiently small $|\varepsilon|$. In this case,
conclusions in Lemma \ref{lem1} are true too.
\end{corollary}

For convenience,  we call the functions $F_k^0(x)~ (k=1, 2)$,
defined in Lemma \ref{lem1}, the first and second order averaged
functions associated with system \eqref{eq3}, respectively.

Consider the planar integrable system of the form
\begin{equation}\label{eq4}
\begin{gathered}
\dot{x}=P(x, y),\\
\dot{y}=Q(x, y),
\end{gathered}
\end{equation}
where $P(x, y), Q(x, y): \mathbb{R}^2\to \mathbb{R}$
are continuous functions such that \eqref{eq4} has a first integral
$H$ with the integrating factor $\mu(x, y)\neq 0$, and has a
continuous family of ovals
$$
\{\gamma_{h}\}\subset\{(x, y):H(x, y)=h,\, h_{c}<h<h_{s}\},
$$
around the center $(0, 0)$. Here $h_{c}$ is the critical level of
$H(x,y)$ corresponding to the center $(0, 0)$ and $h_{s}$ denotes
the value of $H(x,y)$ for which the period annulus terminates at a
separatrix polycycle. Without loss of generality, assume
$h_{s}>h_{c}\geq 0$. We perturb this system as follows
\begin{equation}\label{eq5}
\begin{gathered}
\dot{x}=P(x, y)+\varepsilon p(x, y, \varepsilon),\\
\dot{y}=Q(x, y)+\varepsilon q(x, y, \varepsilon),
\end{gathered}
\end{equation}
where $\varepsilon$ is a small parameter and $p(x, y, \varepsilon),
q(x, y, \varepsilon):  \mathbb{R}^2\times \mathbb{R}\to
\mathbb{R}$ are continuous functions. To study the number
of limit cycles for sufficiently small $|\varepsilon|$ by using the
above averaging theory, we first need to transform system
\eqref{eq5} into the canonical form described in Lemma \ref{lem1}. The
following result developed from \cite{BL2} provides a way for such
transformations.

\begin{lemma}[\cite{BL2}]\label{lem2}
For system \eqref{eq4}, assume $xQ(x, y)-yP(x, y)\neq 0$ for all
$(x, y)$ in the period annulus formed by the ovals $\gamma_{h}$. Let
$\rho:(\sqrt{h_{c}}, \sqrt{h_{s}})\times [0, 2\pi)\to [0,+\infty)$
 be a continuous function such that
$$
H\Big(\rho(R, \varphi)\cos\varphi, \rho(R,
\varphi)\sin\varphi\Big)=R^2
$$
for all $R\in(\sqrt{h_{c}},\sqrt{h_{s}})$ and $\varphi\in[0, 2\pi)$.
Then the differential equation which describes the dependence between
the square root of energy $R=\sqrt{h}$ and the angle $\varphi$
for system \eqref{eq5} is
\begin{equation}
\frac{dR}{d\varphi}
=\varepsilon\frac{\mu(x^2+y^2)(Qp-Pq)}{2R(Qx-Py)
+2R\varepsilon(qx-py)}\Big|_{x=\rho(R,
\varphi)\cos\varphi,\, y=\rho(R, \varphi)\sin\varphi},
\end{equation}
which is equivalent to
\begin{align*}
\frac{dR}{d\varphi}
&=\Big[\varepsilon\frac{\mu(x^2+y^2)(Qp-Pq)}{2R(Qx-Py)}\\
&\quad -\varepsilon^2\frac{\mu(x^2+y^2)(Qp-Pq)(qx-py)}{2R(Qx-Py)^2}
\Big]\Big|_{x=\rho(R, \varphi)\cos\varphi,\, y=\rho(R,
\varphi)\sin\varphi}+O(\varepsilon^3),
\end{align*}
where $P, Q, p$ and $q$ are defined as before.
\end{lemma}

\section{First-order limit cycle bifurcation} \label{sect3}

It is notable that for integrable and non-Hamiltonian systems,
it is generally difficult to find the suitable transformations as
described in Lemma \ref{lem2}.

For the first integral of system \eqref{eq1},
$$
H(x,y)=\frac{1-x}{\sqrt{x^2+y^2}},
$$
we choose the function $\rho=\rho(R, \varphi)$ as follows
\begin{equation}\label{eq6}
\rho(R, \varphi)=\frac{1}{R^2+\cos\varphi}
\end{equation}
such that
$$
H\Big(\rho(R, \varphi)\cos\varphi, \rho(R, \varphi)\sin\varphi\Big)=R^2.
$$
Applying Lemma \ref{lem2} to system
\eqref{eq2}, we obtain the following result.

\begin{lemma}\label{lem3}
With the transformation $x=\rho(R, \varphi)\cos\varphi$ and $y=\rho(R,
\varphi)\sin\varphi$ for $\varphi\in[0, 2\pi)$, system \eqref{eq2}
can be reduced to
\begin{equation}\label{eq7}
\begin{split}
\frac{dR}{d\varphi}
&=\frac{1}{2R}\Big\{\varepsilon\frac{(Qp_1-Pq_1)}{(x^2+y^2)^{3/2}}
 +\varepsilon^2\Big[\frac{Qp_2-Pq_2}{(x^2+y^2)^{3/2}}\\
&\quad -\frac{(Qp_1-Pq_1)(xq_1-yp_1)}{(x^2+y^2)^{5/2}}\Big]
\Big\}\Big|_{x=\rho(R, \varphi)\cos\varphi,\,  y=\rho(R,
\varphi)\sin\varphi}+O(\varepsilon^3),
\end{split}
\end{equation}
where
\begin{gather*}
\begin{split}
Qp_k-Pq_k&=-b_{20}^{(k)}x^3y+\Big(a_{20}^{(k)}-b_{11}^{(k)}\Big)x^2y^2
+\Big(a_{11}^{(k)}-b_{02}^{(k)}\Big)xy^3+a_{02}^{(k)}y^{4}\\
&\quad +a_{20}^{(k)}x^3+\Big(a_{11}^{(k)}+b_{20}^{(k)}\Big)x^2y
+\Big(a_{02}^{(k)}+b_{11}^{(k)}\Big)xy^2 +b_{02}^{(k)}y^3,
\end{split}\\
xq_k-yp_k=b_{20}^{(k)}x^3+\Big(b_{11}^{(k)}-a_{20}^{(k)}\Big)x^2y
+\Big(b_{02}^{(k)}-a_{11}^{(k)}\Big)xy^2-a_{02}^{(k)}y^3,
\end{gather*}
and $a_{ij}^{(k)}$ and $b_{ij}^{(k)}$ $(i, j=0, 1, 2, k=1, 2)$ are
real, and $\rho(R, \varphi)$ is given by \eqref{eq6}.
\end{lemma}

Now we begin with the computation of the first-order averaged
function of system \eqref{eq7}.
A straightforward calculation gives the following lemma.

\begin{lemma}\label{lem4}
The following integral equalities hold:
\begin{gather*}
\int_0^{2\pi}\frac{\cos^2\varphi\sin^2\varphi}{R^2
 +\cos\varphi}d\varphi=\pi\Big(2R^{6}-R^2
 -2\frac{R^{8}-R^{4}}{\sqrt{R^{4}-1}}\Big),\\
\int_0^{2\pi}\frac{\sin^{4}\varphi}{R^2
 +\cos\varphi}d\varphi=\pi\Big(-2R^{6}+3R^2
 +2\frac{R^{8}-2R^{4}+1}{\sqrt{R^{4}-1}}\Big).
\end{gather*}
\end{lemma}

\begin{proposition}\label{prop1}
The first order averaged function associated with system \eqref{eq7}
has at most one simple zero, and this upper bound can be reached.
\end{proposition}

\begin{proof}
 The first-order averaged equation corresponding to system
\eqref{eq7} is
\begin{equation}\label{eq8}
\dot{R}=\varepsilon F_1^0(R),
\end{equation}
where
\begin{equation}\label{eq9}
\begin{split}
F_1^0(R)
&=\frac{1}{2\pi}\int_0^{2\pi}
 \Big[\frac{Qp_1-Pq_1}{2R(x^2+y^2)^{3/2}}\Big]
\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi\\
&=\frac{1}{4\pi R}\int_0^{2\pi}\Big\{\frac{1}{R^2+\cos\varphi}
\Big[\Big(a_{20}^{(1)}-b_{11}^{(1)}\Big)\cos^2\varphi\sin^2
\varphi \\
&\quad +a_{02}^{(1)}\sin^{4}\varphi \Big]\Big\}d\varphi.
\end{split}
\end{equation}
Using Lemma \ref{lem4} in \eqref{eq9}, we obtain
\begin{equation}\label{eq10}
\begin{split}
F_1^0(R)&=\frac{1}{4R}\Big\{\Big(2a_{20}^{(1)}-2b_{11}^{(1)}-2a_{02}^{(1)}
 \Big)R^{6}+\Big(-a_{20}^{(1)}+b_{11}^{(1)}+3a_{02}^{(1)}\Big)R^2\\
&\quad +\Big[\Big(-2a_{20}^{(1)}+2b_{11}^{(1)}+2a_{02}^{(1)}\Big)R^{8}\\
&\quad +\Big(2a_{20}^{(1)}-2b_{11}^{(1)}-4a_{02}^{(1)}\Big)R^{4}
+2a_{02}^{(1)}\Big]\frac{1}{\sqrt{R^{4}-1}}\Big\}.
\end{split}
\end{equation}
Recall that $R>1$, and let
$$
R^2=\frac{1+w^2}{1-w^2}
$$
for $0<w<1$. Then formula \eqref{eq10} becomes
\begin{equation}\label{eq11}
\begin{split}
g(w):=&F_1^0(R)\Big|_{R^2=(1+w^2)/(1-w^2)}\\
=&\frac{\sqrt{1-w^2}}{4\sqrt{1+w^2}}
 \Big\{\Big(2a_{20}^{(1)}-2b_{11}^{(1)}-2a_{02}^{(1)}\Big)
 \frac{(1+w^2)^3}{(1-w^2)^3} \\
&\quad +\Big(-a_{20}^{(1)}+b_{11}^{(1)}+3a_{02}^{(1)}\Big)\frac{1+w^2}{1-w^2}\\
&+\Big(-a_{20}^{(1)}+b_{11}^{(1)}+a_{02}^{(1)}\Big)
 \frac{(1+w^2)^{4}}{w(1-w^2)^3} \\
&\quad +\Big(a_{20}^{(1)}-b_{11}^{(1)}-2a_{02}^{(1)}\Big)
 \frac{(1+w^2)^2}{w(1-w^2)}+a_{02}^{(1)}\frac{1-w^2}{w}\Big\}\\
&=\frac{(1-w)^{2/3}}{4(1+w^2)^{1/2}(1+w)^{5/2}}
\Big[\Big(a_{20}^{(1)}-b_{11}^{(1)}+a_{02}^{(1)}\Big)w^2 \\
&\quad +4a_{02}^{(1)}w +a_{20}^{(1)}-b_{11}^{(1)}+a_{02}^{(1)}\Big].
\end{split}
\end{equation}
For the function $g(w)$, we know that if $w_0\neq0$ is one root of
$g(w)=0$, so is $1/w_0$. Thus $g(w)$ has at most one zero in
$w\in(0,1)$ , which implies that there exists at most one zero for
$F_1^0(R)$
in $R\in(1,\infty)$. We will show that this upper bound can be reached
by illustrating an example.
Consider a family of systems
\begin{equation}\label{eq12}
\begin{gathered}
\dot{x}=-y+xy+\varepsilon\Big[\Big(b_{11}^{(1)}+\frac{13}{8}\Big)x^2
+a_{11}^{(1)}xy -\frac{5}{8}y^2\Big],\\
\dot{y}=x+y^2+\varepsilon\Big(b_{20}^{(1)}x^2+b_{11}^{(1)}xy+b_{02}^{(1)}
y^2\Big),
\end{gathered}
\end{equation}
where $a_{11}^{(1)}, b_{20}^{(1)}, b_{11}^{(1)}$ and $ b_{02}^{(1)}$
are real. In the polar coordinates $x=\rho(R, \varphi)\cos\varphi$
and $y=\rho(R, \varphi)\sin\varphi$, system \eqref{eq12} can be
rewritten as
\begin{equation}\label{eq13}
\frac{dR}{d\varphi}=\varepsilon G(R, \varphi)+O(\varepsilon^2),
\end{equation}
where
\begin{align*}
G(R, \varphi)
&=\Big[\frac{Qp_1-Pq_1}{2R(x^2+y^2)^{3/2}}\Big]
\Big|_{x=\rho(R,\varphi)\cos\varphi,\, y=\rho(R,\varphi)\sin\varphi}\\
&=\frac{1}{2R}\Big\{\frac{1}{R^2+cos\varphi}
 \Big[-b_{20}^{(1)}\cos^3\varphi\sin\varphi+\frac{13}{8}
\cos^2\varphi\sin\varphi \\
&\quad +\big(a_{11}^{(1)}-b_{02}^{(1)}\big)
\cos\varphi\sin^3\varphi-\frac{5}{8}\cos^{4}\varphi\Big]\\
&\quad +\Big[\big(b_{11}^{(1)}+\frac{13}{8}\big)\cos^3\varphi
+\big(a_{11}^{(1)}+b_{20}^{(1)}\big)\cos^2\varphi\sin\varphi\\
&\quad +\big(b_{11}^{(1)}-\frac{5}{8}\big)\cos\varphi\sin^2\varphi
 +b_{02}^{(1)}\sin^3\varphi\Big]\Big\}.
\end{align*}
So the first order averaged equation of system \eqref{eq13} is
\begin{equation}
\frac{dR}{d\varphi}=\varepsilon G_1^0(R),
\end{equation}
where \begin{equation}
\begin{split}
G_1^0(R)
&=\frac{1}{2\pi}\int_0^{2\pi}G_1(R, \varphi)d\varphi\\
&=\frac{1}{4\pi R}\int_0^{2\pi}\frac{1}{R^2
 +\cos\varphi}\Big(\frac{13}{8}\cos^2\varphi\sin^2\varphi
 -\frac{5}{8}\sin^{4}\varphi\Big)d\varphi\\
&=\frac{1}{16R}\Big[18R^{6}-14R^2+\left(-18R^{8}+23R^{4}-5\right)
 \frac{1}{\sqrt{R^{4}-1}}\Big]\\
&=\frac{(1-w)^{3/2}}{4(1+w^2)^{1/2}
 (1+w)^{5/2}}\big(w-\frac{1}{2}\big)(w-2),
\end{split}
\end{equation}
where $R$ and $w$ are defined as before. Apparently,
$G_1^0(R)$ has exactly one positive zero, denoted by
\begin{equation}
R_0^{(1)}=\frac{\sqrt{15}}{3},
\end{equation}
corresponding to $w_0^{(1)}=1/2$ in $R\in(1, +\infty)$. Moreover,
we have
\begin{equation}
\frac{d}{dR} G_1^0\Big(R_0^{(1)}\Big)=-\frac{1}{32}<0.
\end{equation}
This completes the proof of Proposition \ref{prop1}.
\end{proof}

On the basis of Lemma \ref{lem1}, Corollary \ref{coro1} and Proposition
\ref{prop1}, we have the following proposition.

\begin{proposition}\label{prop2}
For $|\varepsilon|\neq 0$ sufficiently small, system \eqref{eq2} has
at most one limit cycle for the first order bifurcation arising from
the period annulus around the center of the unperturbed system
\eqref{eq2} with $\varepsilon=0$, and this upper bound is sharp.
\end{proposition}


\section{Second-order limit cycle bifurcation} \label{sect4}

In this section, we study the number of the zeros of second-order
averaged function associated with system \eqref{eq7}, in the case
where the first order averaged function $F_1^0(R)\equiv 0$
holds.
On the basis of formula \eqref{eq11}, we obtain

\begin{lemma}\label{l5}
For system \eqref{eq7}, the first-order averaged function
$F_1^0(R)\equiv 0$ holds if and only if
\begin{equation}\label{eq14}
a_{20}^{(1)}=b_{11}^{(1)},\quad a_{02}^{(1)}=0.
\end{equation}
\end{lemma}

When condition \eqref{eq14} holds, the
second-order averaged function associated with system \eqref{eq7}
takes the form
\begin{equation}\label{eq15}
F_2^0(R)=\frac{1}{2\pi}\int_0^{2\pi}\Big[ \frac{\partial
F_1(R, \varphi)}{\partial R} y_1(R, \varphi)+F_2(R,
\varphi)\Big]d\varphi,
\end{equation}
where
\begin{gather*}
\begin{split}
&F_1(R, \varphi)\\
&=\frac{(Qp_1-Pq_1)}{2R(x^2+y^2)^{3/2}}\Big|_{x=\rho\cos\varphi,\,
  y=\rho\sin\varphi}\\
&=\frac{1}{2R}\Big[-b_{20}^{(1)}\frac{\cos^3\varphi\sin\varphi}{R^2+\cos\varphi}
+\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\frac{\cos\varphi\sin^3\varphi}{R^2
 +\cos\varphi}+\Big(a_{11}^{(1)}
+b_{20}^{(1)}\Big)\cos^2\varphi\sin\varphi\\
&\quad +b_{02}^{(1)}\sin^3\varphi+b_{11}^{(1)}\cos\varphi\Big],
\end{split} \\
\begin{split}
&F_2(R,\varphi) \\
&=\Big[\frac{Qp_2-Pq_2}{2R(x^2+y^2)^{3/2}}-\frac{(Qp_1-Pq_1)
 (xq_1-yp_1)}{2R(x^2+y^2)^{5/2}}\Big]
\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\
&=\Big[\frac{Qp_2-Pq_2}{2R(x^2+y^2)^{3/2}}\Big]
 \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\
&\quad -\frac{1}{2R\Big(R^2+\cos\varphi\Big)^2}
 \Big[-\Big(b_{20}^{(1)}\Big)^2\cos^{6}\varphi\sin\varphi
+2b_{20}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\cos^{4}\varphi\sin^3\varphi\\
&\quad -\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)^2\cos^2\varphi\sin^{5}\varphi\Big]\\
&\quad -\frac{1}{2R\Big(R^2+\cos\varphi\Big)}\Big\{b_{11}^{(1)}b_{20}^{(1)}
 \cos^{6}\varphi
+b_{11}^{(1)}\Big(b_{20}^{(1)}+b_{02}^{(1)}-a_{11}^{(1)}\Big)
 \cos^{4}\varphi\sin^2\varphi\\
&\quad +b_{20}^{(1)}\Big(a_{11}^{(1)}+b_{20}^{(1)}\Big)\cos^{5}\varphi\sin\varphi
+\Big[\Big(a_{11}^{(1)}+b_{20}^{(1)}\Big)\Big(b_{02}^{(1)}-a_{11}^{(1)}\Big) \\
&\quad +b_{02}^{(1)}b_{20}^{(1)}\Big]\cos^3\varphi\sin^3\varphi
 -b_{11}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\cos^2\varphi\sin^{4}\varphi \\
&\quad -b_{02}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\cos\varphi\sin^{5}\varphi\Big\},
\end{split} \\
\begin{split}
&y_1(R, \varphi)\\
&=\int_0^{\varphi}F_1(R, \theta)d\theta\\
&=\int_0^{\varphi}\frac{1}{2R}\Big[b_{11}^{(1)}\cos^3\theta
 +\Big(a_{11}^{(1)}+b_{20}^{(1)}\Big)\cos^2\theta\sin\theta
+b_{11}^{(1)}\cos\theta\sin^2\theta+b_{02}^{(1)}\sin^3\theta \Big]d\theta\\
&\quad +\frac{1}{2R}\Big[-b_{20}^{(1)}\int_0^{\varphi}
 \frac{\cos^3\theta\sin\theta}{R^2+\cos\theta}d\theta
+\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)
 \int_0^{\varphi}\frac{\cos\theta\sin^3\theta}{R^2+\cos\theta}d\theta\Big],
\end{split}
\end{gather*}
and $P, Q, p_k$ and $q_k$ $(k=1, 2)$ are defined as before.

To compute the function $y_1(R, \varphi)$, in the following we first need to 
figure out some integral equalities.

\begin{lemma}\label{lem6} The
following integral equalities hold:
\begin{gather*}
\int_0^{\varphi}\frac{\cos\theta}{R^2+\cos\theta}d\cos\theta
=-1+R^2\ln\left(R^2+1\right)+\cos\varphi-R^2\ln\left(R^2+\cos\varphi\right),\\
\begin{aligned}
\int_0^{\varphi}\frac{\cos^3\theta}{R^2+\cos\theta}d\cos\theta
&=-R^{4}+\frac{1}{2}R^2-\frac{1}{3}+R^{6}\ln\left(R^2+1\right)
 +R^{4}\cos\varphi-\frac{1}{2}R^2\cos^2\varphi\\
&\quad +\frac{1}{3}\cos^3\varphi-R^{6}\ln\big(R^2+\cos\varphi\big).
\end{aligned}
\end{gather*}
\end{lemma}

Based on Lemma \ref{lem6}, we obtain the following lemma.

\begin{lemma}\label{lem7} 
The following integral equalities hold:
\begin{gather*}
\begin{split}
\int_0^{\varphi}\frac{\cos^3\theta\sin\theta}{R^2+\cos\theta}d\theta
&=R^{4}-\frac{1}{2}R^2+\frac{1}{3}-R^{6}\ln\left(R^2+1\right)
 -R^{4}\cos\varphi \\
&\quad +\frac{1}{2}R^2\cos^2\varphi
 -\frac{1}{3}\cos^3\varphi+R^{6}\ln\left(R^2+\cos\varphi\right),
\end{split}\\
\begin{split}
\int_0^{\varphi}\frac{\cos\theta\sin^3\theta}{1+R^2\cos^2\theta}d\theta
&=-R^{4}+\frac{1}{2}R^2+\frac{2}{3}+\left(R^{6}-R^2\right)\ln\left(R^2+1\right)
 +\left(R^{4}-1\right)\cos\varphi \\
&\quad -\frac{1}{2}R^2\cos^2\varphi
 +\frac{1}{3}\cos^3\varphi+\left(R^2-R^{6}\right)
 \ln\left(R^2+\cos\varphi\right).
\end{split}
\end{gather*}
\end{lemma}

By using Lemmas \ref{lem6} and \ref{lem7}, a straightforward computation
yields
\begin{equation}
\begin{split}
&y_1(R,\varphi)\\
&=-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{4}R\cos^2\varphi
+\Big[\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^3
 -\frac{a_{11}^{(1)}}{2R}\Big]\cos\varphi \\
&\quad +\frac{b_{11}^{(1)}}{2R}\sin\varphi
 +\Big[-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^{5}
 +\frac{a_{11}^{(1)}-b_{02}^{(1)}}{2}R\Big]\ln\left(R^2+\cos\varphi\right) \\
&\quad +\frac{a_{11}^{(1)}}{2R}+\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{4}R
-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^3 \\
&\quad +\Big[\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^{5}
 -\frac{a_{11}^{(1)}-b_{02}^{(1)}}{2}R\Big]\ln\left(R^2+1\right).
\end{split}
\end{equation}

\begin{lemma}\label{lem8}
 The following integral equalities are true:
\begin{gather*}
\int_0^{2\pi}\frac{1}{R^2+\cos\varphi}d\varphi=\frac{2\pi}{\sqrt{R^{4}-1}},\\
\int_0^{2\pi}\frac{\cos\varphi}{R^2+\cos\varphi}d\varphi
=2\pi\Big[-\frac{R^2}{\sqrt{R^{4}-1}}+1\Big],\\
\int_0^{2\pi}\frac{\cos^2\varphi}{R^2+\cos\varphi}d\varphi
=2\pi\Big[\frac{R^{4}}{\sqrt{R^{4}-1}}-R^2\Big],\\
\int_0^{2\pi}\frac{\cos^3\varphi}{R^2+\cos\varphi}d\varphi
=\pi\Big[-\frac{2R^{6}}{\sqrt{R^{4}-1}}+2R^{4}+1\Big],\\
\int_0^{2\pi}\frac{\cos^{4}\varphi}{R^2+\cos\varphi}d\varphi
=\pi\Big[\frac{2R^{8}}{\sqrt{R^{4}-1}}-2R^{6}-R^2\Big],\\
\int_0^{2\pi}\frac{\cos^{5}\varphi}{R^2+\cos\varphi}d\varphi
=\frac{\pi}{4}\Big[-\frac{8R^{10}}{\sqrt{R^{4}-1}}+8R^{8}+4R^{4}+3\Big],\\
\int_0^{2\pi}\frac{\cos^{6}\varphi}{R^2+\cos\varphi}d\varphi
=-\frac{\pi}{4}\Big[-\frac{8R^{12}}{\sqrt{R^{4}-1}}+8R^{10}+4R^{6}+3R^2\Big],\\
\int_0^{2\pi}\cos\varphi\ln\Big(R^2+\cos\varphi\Big)d\varphi
=2\pi\big[R^2-\sqrt{R^{4}-1}\big].
\end{gather*}
\end{lemma}

\begin{proof}
Most of integral equalities can be obtained by a direct
computation. Here we only show the derivation of the last integral formula.
Let
\begin{equation}
N(r)=\int_0^{2\pi}\cos\varphi\ln\big(1+r\cos\varphi\big)d\varphi,
\end{equation}
where $r=1/R^2$.
 Since
\begin{equation}
N'(r)=\int_0^{2\pi}\frac{\cos^2\varphi}{1+r\cos\varphi}d\varphi
=\frac{2\pi\Big(1-\sqrt{1-r^2}\Big)}{r^2\sqrt{1-r^2}},
\end{equation}
and $N(0)=0$, we obtain
\begin{equation}
N(r)=\int_0^{r}N'(s)ds=2\pi\big(R^2-\sqrt{R^{4}-1}\big),
\end{equation}
which implies the last formula in Lemma \ref{lem8}.
\end{proof}

By  Lemma \ref{lem8}, a straightforward calculation yields the following
lemma.

\begin{lemma}\label{lem9} The following integral equalities are true:
\begin{gather*}
\int_0^{2\pi}\frac{\cos\varphi\sin^{4}\varphi}{R^2+cos\varphi}d\varphi
=\pi\Big[-2(R^{6}-R^2)\sqrt{R^{4}-1}+2R^{8}-3R^{4}+\frac{3}{4}\Big],\\
\int_0^{2\pi}\frac{\cos^2\varphi\sin^{4}\varphi}{R^2+cos\varphi}d\varphi
=-\frac{\pi}{4}\Big[\frac{-8R^{12}+16R^{8}-8R^{4}}{\sqrt{R^{4}-1}}
 +8R^{10}-12R^{6}+3R^2\Big],\\
\int_0^{2\pi}\frac{\cos^3\varphi\sin^2\varphi}{R^2+cos\varphi}d\varphi
=\pi\Big[2R^{6}\sqrt{R^{4}-1}-2R^{8}+R^{4}+\frac{1}{4}\Big],\\
\int_0^{2\pi}\frac{\cos^{4}\varphi\sin^2\varphi}{R^2+cos\varphi}d\varphi
=\frac{\pi}{4}\Big[\frac{-8R^{12}+8R^{8}}{\sqrt{R^{4}-1}}+8R^{10}-4R^{6}-R^2\Big],\\
\int_0^{2\pi}\frac{\cos\varphi\sin^{4}\varphi}{(R^2+cos\varphi)^2}d\varphi
=-2\pi\Big[\frac{-4R^{8}+5R^{4}-1}{\sqrt{R^{4}-1}}+4R^{6}-3R^2\Big],\\
\int_0^{2\pi}\frac{\cos^3\varphi\sin^2\varphi}{(R^2+cos\varphi)^2}d\varphi
=2\pi\Big[\frac{-4R^{8}+3R^{4}}{\sqrt{R^{4}-1}}+4R^{6}-R^2\Big].
\end{gather*}
\end{lemma}

\begin{proposition}\label{prop3}
Under condition \eqref{eq14}, the second order averaged function
associated with system \eqref{eq7} has at most two simple zeros, and
this upper bound can be reached.
\end{proposition}

\begin{proof}
Define
\[
F_{21}^0(R)=\frac{1}{2\pi}\int_0^{2\pi}\frac{\partial F_1(R,
\varphi)}{\partial R} y_1(R, \varphi)d\varphi,\quad
F_{22}^0(R)=\frac{1}{2\pi}\int_0^{2\pi}F_2(R,
\varphi)d\varphi.
\]
Then \eqref{eq15} becomes
\begin{equation}
F_2^0(R)=F_{21}^0(R)+F_{22}^0(R).
\end{equation}
\smallskip

\noindent\textbf{Step 1.}  Computation of the function $F_{21}^0(R)$.
Let
\begin{gather*}
\begin{split}
A_1&=-\frac{1}{2R^2}\Big[-b_{20}^{(1)}\frac{\cos^3\varphi\sin\varphi}
 {R^2+\cos\varphi}
+\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\frac{\cos\varphi\sin^3\varphi}
 {R^2+\cos\varphi}\\
&\quad +\Big(a_{11}^{(1)}+b_{20}^{(1)}\Big)\cos^2\varphi\sin\varphi
 +b_{02}^{(1)}\sin^3\varphi\Big]+b_{20}^{(1)}
 \frac{\cos^3\varphi\sin\varphi}{(R^2+\cos\varphi)^2} \\
&\quad-\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)
 \frac{\cos\varphi\sin^3\varphi}{(R^2+\cos\varphi)^2},
\end{split}\\
A_2=-\frac{b_{11}^{(1)}}{2R^2}\cos\varphi,\\
\begin{split}
B_1=&-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{4}R\cos^2\varphi
+\Big[\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^3
 -\frac{a_{11}^{(1)}}{2R}\Big]\cos\varphi\\
&+\Big[-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^{5}
 +\frac{a_{11}^{(1)}-b_{02}^{(1)}}{2}R\Big]\ln\left(R^2+\cos\varphi\right),
\end{split}\\
\begin{split}
B_2=&\frac{b_{11}^{(1)}}{2R}\sin\varphi+\frac{a_{11}^{(1)}}{2R}
 +\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{4}R
-\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^3\\
&+\Big[\frac{b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}}{2}R^{5}
 -\frac{a_{11}^{(1)}-b_{02}^{(1)}}{2}R\Big]\ln\left(R^2+1\right).
\end{split}
\end{gather*}
Then
\[
\frac{\partial F_1(R, \varphi)}{\partial R}=A_1+A_2,\quad  y_1(R, \varphi)=B_1+B_2,
\]
and
\begin{equation} \label{eq16}
F_{21}^0(R)=\frac{1}{2\pi}\int_0^{2\pi}(A_1B_1+A_1B_2+A_2B_1+A_2B_2)d\varphi.
\end{equation}
A straightforward calculation shows
\begin{equation}\label{eq17}
\int_0^{2\pi}A_2B_2d\varphi=0.
\end{equation}
Recalling that the function $A_1B_1$ is odd with respect
to $\varphi$, we obtain
\begin{equation}\label{eq18}
\int_0^{2\pi}A_1B_1d\varphi=0.
\end{equation}
In addition, it is not difficult to verify that
\begin{equation}\label{eq19}
\begin{split}
&\frac{1}{2\pi}\int_0^{2\pi}A_1B_2d\varphi \\
&=\frac{1}{2\pi}\int_0^{2\pi}\Big\{-\frac{b_{11}^{(1)}}{4R^3}
 \Big[-b_{20}^{(1)}\frac{\cos^3\varphi\sin^2\varphi}{R^2+\cos\varphi}
+\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\frac{\cos\varphi\sin^{4}\varphi}
 {R^2+\cos\varphi} \\
&\quad +\Big(a_{11}^{(1)}+b_{20}^{(1)}\Big)\cos^2\varphi\sin^2\varphi
 +b_{02}^{(1)}\sin^{4}\varphi\Big]+\frac{b_{20}^{(1)}b_{11}^{(1)}}
 {2R}\frac{\cos^3\varphi\sin^2\varphi}{(R^2+\cos\varphi)^2}\\
&\quad -\frac{b_{11}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)}{2R}
 \frac{\cos\varphi\sin^{4}\varphi}{(R^2+\cos\varphi)^2}\Big\}d\varphi,
\end{split}
\end{equation}
and
\begin{equation}\label{eq20}
\begin{split}
&\frac{1}{2\pi}\int_0^{2\pi}A_2B_1d\varphi \\
&=\frac{1}{2\pi}\int_0^{2\pi}\Big\{\Big[\frac{-b_{11}^{(1)}
 \Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)}{4}R
+\frac{b_{11}^{(1)}a_{11}^{(1)}}{4R^3}\Big]\cos^2\varphi \\
&\quad  +\Big[\frac{b_{11}^{(1)}\Big(b_{20}^{(1)}+a_{11}^{(1)}
 -b_{02}^{(1)}\Big)}{4}R^3 \\
&\quad -\frac{b_{11}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)}{4R}\Big]
 \cos\varphi\ln\left(R^2+\cos\varphi\right)\Big\}d\varphi.
\end{split}
\end{equation}
Applying Lemmas \ref{lem8} and \ref{lem9} to \eqref{eq19} and
\eqref{eq20} gives
\begin{equation}\label{eq21}
\begin{split}
&\frac{1}{2\pi}\int_0^{2\pi}A_1B_2d\varphi\\
&=-\frac{b_{11}^{(1)}}{8R^3}\Big\{\Big[-2\Big(b_{20}^{(1)}
+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{6}
+2\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)R^2\Big]\sqrt{R^{4}-1}
\\
&\quad -14\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{8}
 +\Big(3b_{20}^{(1)}+9a_{11}^{(1)}-9b_{02}^{(1)}\Big)R^{4}+a_{11}^{(1)}\\
&\quad +\Big(16\big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\big)
 R^{10}-\big(12b_{20}^{(1)}+20a_{11}^{(1)}-20b_{02}^{(1)}\big)R^{6} \\
&\quad +4\big(a_{11}^{(1)}-b_{02}^{(1)}\big)R^2\Big)\big/ \sqrt{R^{4}-1} \Big\},
\end{split}
\end{equation}
and 
\begin{equation}\label{eq22}
\begin{split}
&\frac{1}{2\pi}\int_0^{2\pi}A_2B_1d\varphi\\
&=-\frac{b_{11}^{(1)}}{8R^3}\Big\{\Big[2\Big(b_{20}^{(1)}
 +a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{6}
-2\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)R^2\Big]\sqrt{R^{4}-1}\\
&\quad -2\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{8}
 +\Big(b_{20}^{(1)}+3a_{11}^{(1)}-3b_{02}^{(1)}\Big)R^{4}-a_{11}^{(1)}\Big\}.
\end{split}
\end{equation}
Substituting \eqref{eq17}, \eqref{eq18}, \eqref{eq21} and
\eqref{eq22} in \eqref{eq16} yields
\begin{equation}\label{eq23}
\begin{split}
F_{21}^0(R) 
&=-\frac{b_{11}^{(1)}}{2R^3}
 \Big\{-4\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{8}
+\Big(b_{20}^{(1)}+3a_{11}^{(1)}-3b_{02}^{(1)}\Big)R^{4}\\
&\quad +\Big(4\big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\big)R^{10}
 -\big(3b_{20}^{(1)}+5a_{11}^{(1)}-5b_{02}^{(1)}\big)R^{6} \\
&\quad +\big(a_{11}^{(1)}-b_{02}^{(1)}\big)R^2\Big)\big/ \sqrt{R^{4}-1}\Big\}.
\end{split}
\end{equation}
\smallskip

\noindent\textbf{Step 2.} Computation of the function $F_{22}^0(R)$.
As above, we have
\begin{equation}\label{eq24}
\begin{split}
&\frac{1}{2\pi}\int_0^{2\pi}\Big[\frac{Qp_2-Pq_2}{2R(x^2+y^2)^{3/2}}
\Big]\Big|_{x=\rho\cos\varphi,\,
y=\rho\sin\varphi}d\varphi\\
&=\frac{1}{4\pi R}\Big[\Big(a_{20}^{(2)}-b_{11}^{(2)}\Big)
 \int_0^{2\pi}\frac{\cos^2\varphi\sin^2\varphi}{R^2+\cos\varphi}d\varphi
+a_{02}^{(2)}\int_0^{2\pi}\frac{\sin^{4}\varphi}{R^2+\cos\varphi}d\varphi\Big]\\
&=\frac{1}{4R}\Big\{\Big(2a_{20}^{(2)}-2b_{11}^{(2)}-2a_{02}^{(2)}\Big)R^{6}
 +\Big(-a_{20}^{(2)}+b_{11}^{(2)}+3a_{02}^{(2)}\Big)R^2\\
&\quad +\Big[\Big(-2a_{20}^{(2)}+2b_{11}^{(2)}+2a_{02}^{(2)}\Big)R^{8}
 +\Big(2a_{20}^{(2)}-2b_{11}^{(2)}-4a_{02}^{(2)}\Big)R^{4} \\
&\quad +2a_{02}^{(2)}\Big]\frac{1}{\sqrt{R^{4}-1}}\Big\}.
\end{split}
\end{equation}
Using Lemmas \ref{lem8} and \ref{lem9}, we obtain
\begin{equation}\label{eq25}
\begin{split}
&-\frac{1}{2\pi}\int_0^{2\pi}\Big[\frac{(Qp_1-Pq_1)(xq_1-yp_1)}
 {2R(x^2+y^2)^{5/2}}\Big]
\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi \\
&=-\frac{1}{4\pi R}\Big[b_{11}^{(1)}b_{20}^{(1)}\int_0^{2\pi}
 \frac{\cos^{6}\varphi}{R^2+\cos\varphi}d\varphi \\
&\quad +b_{11}^{(1)}\Big(b_{20}^{(1)}+b_{02}^{(1)}-a_{11}^{(1)}\Big)
 \int_0^{2\pi}\frac{\cos^{4}\varphi\sin^2\varphi}{R^2+\cos\varphi}d\varphi\\
&\quad -b_{11}^{(1)}\left(a_{11}^{(1)}-b_{02}^{(1)}\right)\int_0^{2\pi}
\frac{\cos^2\varphi\sin^{4}\varphi}{R^2+\cos\varphi}d\varphi\Big]\\
&=-\frac{1}{4R}\Big\{-2b_{11}^{(1)}\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{6}
+b_{11}^{(1)}\Big(-b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^2\\
&\quad +\Big[2b_{11}^{(1)}\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{8}
-2b_{11}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)R^{4}\Big]
 \frac{1}{\sqrt{R^{4}-1}}\Big\}.
\end{split}
\end{equation}
It follows from \eqref{eq24} and \eqref{eq25} that
\begin{equation}\label{eq26}
F_{22}^0(R)=\frac{1}{4R}\Big(l_1R^{6}+l_2R^2+\frac{l_3}{\sqrt{R^{4}-1}}\Big).
\end{equation}
where
\begin{gather*}
l_1=2a_{20}^{(2)}-2b_{11}^{(2)}-2a_{02}^{(2)}+2b_{11}^{(1)}\Big(b_{20}^{(1)}
 +a_{11}^{(1)}-b_{02}^{(1)}\Big),\\
l_2=-a_{20}^{(2)}+b_{11}^{(2)}+3a_{02}^{(2)}+b_{11}^{(1)}\Big(b_{20}^{(1)}
 -a_{11}^{(1)}+b_{02}^{(1)}\Big),\\
\begin{aligned}
l_3&=\Big[-2a_{20}^{(2)}+2b_{11}^{(2)}+2a_{02}^{(2)}
 -2b_{11}^{(1)}\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big)\Big]R^{8}\\
&\quad +\Big[2a_{20}^{(2)}-2b_{11}^{(2)}-4a_{02}^{(2)}+2b_{11}^{(1)}
 \Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)\Big]R^{4}+2a_{02}^{(2)}.
\end{aligned}
\end{gather*}
Based on \eqref{eq23} and \eqref{eq26}, $F_2^0(R)$ becomes
\begin{equation}\label{eq27}
F_2^0(R)
=-\frac{1}{4R^3}\Big(l_{4}R^{8}+l_{5}R^{4}+\frac{l_{6}}{\sqrt{R^{4}-1}}\Big).
\end{equation}
where
\begin{gather*}
l_{4}=-2a_{20}^{(2)}+2b_{11}^{(2)}+2a_{02}^{(2)}-10b_{11}^{(1)}
 \Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big),\\
l_{5}=a_{20}^{(2)}-b_{11}^{(2)}-3a_{02}^{(2)}+b_{11}^{(1)}
 \Big(b_{20}^{(1)}+7a_{11}^{(1)}-7b_{02}^{(1)}\Big),\\
l_{6}=l_{6,1}R^{10}+l_{6,2}R^{6}+l_{6,3}R^2,\\
l_{6,1}=2a_{20}^{(2)}-2b_{11}^{(2)}-2a_{02}^{(2)}+10b_{11}^{(1)}
 \Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big),\\
l_{6,2}=-2a_{20}^{(2)}+2b_{11}^{(2)}+4a_{02}^{(2)}-6b_{11}^{(1)}
 \Big(b_{20}^{(1)}+2a_{11}^{(1)}-2b_{02}^{(1)}\Big),\\
l_{6,3}=2b_{11}^{(1)}\Big(a_{11}^{(1)}-b_{02}^{(1)}\Big)-2a_{02}^{(2)}.
\end{gather*}

After making the same transformations as before, \eqref{eq27}
becomes
\begin{equation}\label{eq28}
\begin{split}
&F_2^0(R)\\
&=-\frac{1}{8w(1+w^2)^{1/2}(1-w^2)^{5/2}}\Big[2l_{4}w(1+w^2)^3+2l_{5}w(1+w^2)
(1-w^2)^2\\
&\quad +l_{6,1}(1+w^2)^4+l_{6,2}(1+w^2)^2(1-w^2)^2+l_{6,3}(1-w^2)^4\Big]\\
&=-\frac{(1-w)^{3/2}}{4w(1+w^2)^{1/2}(1+w)^{5/2}}
\Big[N_{2,1}w^{4}+N_{2,2}w^3+N_{2,3}w^2 \\
&\quad +N_{2,2}w+N_{2,1}\Big]\\
\end{split}
\end{equation}
where
\begin{equation}
\begin{split}
&N_{2,1}=2b_{11}^{(1)}b_{20}^{(1)},\\
&N_{2,2}=-a_{20}^{(2)}+b_{11}^{(2)}-a_{02}^{(2)}-b_{11}^{(1)}\Big(b_{20}^{(1)}+3a_{11}^{(1)}-3b_{02}^{(1)}\Big),\\
&N_{2,3}=-4a_{02}^{(2)}+4b_{11}^{(1)}\Big(b_{20}^{(1)}+a_{11}^{(1)}-b_{02}^{(1)}\Big).
\end{split}
\end{equation}
This shows that the second order averaged function $F_2^0(R)$
associated with system \eqref{eq7} has at most two zeros in $R\in(1,
+\infty)$, by taking into account the multiplicity.

Next we will provide
an example to demonstrate that this upper bound can be reached.
Consider the system
\begin{equation}\label{eq29}
\begin{gathered}
\dot{x}=-y+xy+\varepsilon\Big[x^2+\Big(b_{02}^{(1)}+\frac{25}{12}\Big)xy\Big]
+\varepsilon^2\Big[\Big(b_{11}^{(2)}-\frac{11}{12}\Big)x^2+a_{11}^{(2)}xy\Big],\\
\dot{y}=x+y^2+\varepsilon\Big[\frac{1}{2}x^2+xy+b_{02}^{(1)}y^2\Big]
+\varepsilon^2\Big[b_{20}^{(2)}x^2+b_{11}^{(2)}xy+b_{02}^{(2)}y^2\Big],
\end{gathered}
\end{equation}
where $b_{02}^{(1)}$, $a_{11}^{(2)}$ and $b_{ij}^{(2)}(i,j=0,1,2)$ are
real. In the polar coordinates $x=\rho(R, \varphi)\cos\varphi$ and
$y=\rho(R, \varphi)\sin\varphi$, system \eqref{eq29} becomes
\begin{equation}\label{eq30}
\frac{dR}{d\varphi}=\varepsilon M_1(R,
\varphi)+\varepsilon^2M_2(R, \varphi)+O(\varepsilon^3),
\end{equation}
where
\begin{gather*}
\begin{split}
M_1(R,\varphi)
&=\frac{1}{2R}\Big[-\frac{\cos^3\varphi\sin\varphi}{2(R^2+\cos\varphi)}
+\frac{25\cos\varphi\sin^3\varphi}{12(R^2+\cos\varphi)}+\cos^3\varphi\\
&\quad +\Big(b_{02}^{(1)}+\frac{31}{12}\Big)\cos^2\varphi\sin\varphi
+\cos\varphi\sin^2\varphi+b_{02}^{(1)}\sin^3\varphi\Big],
\end{split}\\
\begin{split}
M_2(R, \varphi)
&=\frac{1}{2R}\Big[\frac{\cos^{6}\varphi\sin\varphi}{4(R^2+\cos\varphi)^2}
-\frac{25\cos^{4}\varphi\sin^3\varphi}{12(R^2+\cos\varphi)^2}
+\frac{625\cos^2\varphi\sin^{5}\varphi}{144(R^2+\cos\varphi)^2}\\
&\quad -\frac{\cos^{6}\varphi}{2(R^2+\cos\varphi)}
 -\Big(b_{02}^{(1)}+\frac{31}{12}\Big)
 \frac{\cos^{5}\varphi\sin\varphi}{2(R^2+\cos\varphi)}
+\frac{19\cos^{4}\varphi\sin^2\varphi}{12(R^2+\cos\varphi)} \\
&\quad +\Big(\frac{19}{12}b_{02}^{(1)}
+\frac{775}{144}\Big)\frac{\cos^3\varphi\sin^3\varphi}{R^2+\cos\varphi}
 +\frac{25\cos^2\varphi\sin^{4}\varphi}{12(R^2+\cos\varphi)} \\
&\quad +\frac{25}{12}b_{02}^{(1)}\frac{\cos\varphi\sin^{5}\varphi}
 {R^2+\cos\varphi}-b_{20}^{(2)}\frac{\cos^3\varphi\sin\varphi}{R^2+\cos\varphi}\\
&\quad -\frac{11\cos^2\varphi\sin^2\varphi}{12(R^2+\cos\varphi)}
+\Big(a_{11}^{(2)}-b_{02}^{(2)}\Big)\frac{\cos\varphi\sin^3\varphi}
 {R^2+\cos\varphi}+\Big(b_{11}^{(2)}-\frac{11}{12}\Big)\cos^3\varphi\\
&\quad +\Big(a_{11}^{(2)}+b_{20}^{(2)}\Big)\cos^2\varphi\sin\varphi
 +b_{11}^{(2)}\cos\varphi\sin^2\varphi+b_{02}^{(2)}\sin^3\varphi\Big].
\end{split}
\end{gather*}
It is not difficult to verify that for system \eqref{eq30}, the
first-order averaged function $M_1^0(R)$ is identically equal to
zero while the second-order averaged function $M_2^0(R)$ takes
the form
\begin{align*}
M_2^0(R)
=&-\frac{1}{8\pi R^3}\Big[-\frac{1}{2}
 \int_0^{2\pi}\frac{\cos^3\varphi\sin^2\varphi}{R^2+\cos\varphi} d\varphi
+\frac{25}{12}\int_0^{2\pi}\frac{\cos\varphi\sin^{4}\varphi}
 {R^2+\cos\varphi}d\varphi\\
&+\Big(b_{02}^{(1)}+\frac{31}{12}\Big)\int_0^{2\pi}\cos^2\varphi\sin^2\varphi
d\varphi +b_{02}^{(1)}\int_0^{2\pi}\sin^{4}\varphi d\varphi \\
&-R^2\int_0^{2\pi}\frac{\cos^3\varphi\sin^2\varphi}{(R^2+\cos\varphi)^2}
d\varphi
+\frac{25}{6}R^2\int_0^{2\pi}\frac{\cos\varphi\sin^{4}\varphi}
{(R^2+\cos\varphi)^2}d\varphi \\
&+\Big(\frac{31}{12}R^{4}-b_{02}^{(1)}-\frac{25}{12}\Big)
 \int_0^{2\pi}\cos^2\varphi d\varphi\\
&+\Big(-\frac{31}{12}R^{6}+\frac{25}{12}R^2\Big)\int_0^{2\pi}
 \cos\varphi \ln(R^2+\cos\varphi) d\varphi\Big]\\
&+\frac{1}{4\pi R}\Big[-\frac{1}{2}\int_0^{2\pi}
 \frac{\cos^{6}\varphi}{R^2+\cos\varphi}d\varphi
+\frac{19}{12}\int_0^{2\pi}\frac{\cos^{4}\varphi\sin^2\varphi}
 {R^2+\cos\varphi}d\varphi\\
&+\frac{25}{12}\int_0^{2\pi}\frac{\cos^2\varphi\sin^{4}\varphi}
 {R^2+\cos\varphi}d\varphi
-\frac{11}{12}\int_0^{2\pi}\frac{\cos^2\varphi\sin^2\varphi}
 {R^2+\cos\varphi}d\varphi\Big]\\
=&-\frac{1}{8R^3}\Big[-\frac{124}{3}R^{8}+27R^{4}
 +\Big(\frac{124}{3}R^{10}-\frac{143}{3}R^{6}
 +\frac{25}{3}R^2\Big)\frac{1}{\sqrt{R^{4}-1}}\Big]\\
&+\frac{1}{4R}\Big[\frac{10}{3}R^{6}-\frac{2}{3}R^2
 +\Big(-\frac{10}{3}R^{8}+\frac{7}{3}R^{4}\Big)\frac{1}{\sqrt{R^{4}-1}}\Big]\\
=&-\frac{1}{8R^3}\Big[-48R^{8}+\frac{85}{3}R^{4}
 +\Big(48R^{10}-\frac{157}{3}R^{6}+\frac{25}{3}R^2\Big)\frac{1}{\sqrt{R^{4}-1}}\Big]\\
=&-\frac{1}{8w(1+w^2)^{1/2}(1-w^2)^{5/2}}
 \Big[-48w(1+w^2)^3+\frac{85}{3}w(1+w^2)(1-w^2)^2\\
&+24(1+w^2)^{4}-\frac{157}{6}(1+w^2)^2(1-w^2)^2
 +\frac{25}{6}(1-w^2)^{4}\Big]\\
=&-\frac{(1-w)^{3/2}}{4w(1+w^2)^{1/2}(1+w)^{5/2}}
 \big(w-\frac{1}{2}\big)\big(w-\frac{1}{3}\big)
\Big(w-2\Big)\Big(w-3\Big),
\end{align*}
where $R$ and $w$ are defined as before. Apparently, $M_2^0(R)$
has exactly two positive zeros, denoted by
\begin{equation}
R_1^{(2)}=\frac{\sqrt{15}}{3},\quad R_2^{(2)}=\frac{\sqrt{5}}{2},
\end{equation}
corresponding to $w_1^{(2)}=1/2$ and $w_2^{(2)}=1/3$,
respectively, in $R\in(1, +\infty)$. Moreover, we have
\begin{equation}
\frac{dM_2^0\Big(R_1^{(2)}\Big)}{dR}=-\frac{5}{192}<0, \quad
\frac{dM_2^0\Big(R_2^{(2)}\Big)}{dR}=\frac{5}{27}>0.
\end{equation}
The proof is complete.
\end{proof}

On the basis of Lemma \ref{lem1}, Corollary \ref{coro1} and Proposition
\ref{prop3}, we have the following proposition

\begin{proposition}\label{prop4}
For $|\varepsilon|\neq 0$ sufficiently small, system \eqref{eq2} has
at most two limit cycles for the second order bifurcation arising
from the period annulus around the center of the unperturbed system
\eqref{eq2}$|_{\varepsilon=0}$, and this upper bound is sharp.
\end{proposition}


 Theorem \ref{thm1} follows immediately from
Propositions \ref{prop2} and \ref{prop4}.


\subsection*{Acknowledgements}

This research is supported by the National Science Foundation of
China under grants No. 11371046 and No. 11290141.


\begin{thebibliography}{00}

\bibitem{AI} V. Arnold, Y. Ilyashenko; 
\emph{Dynamical systems I: Ordinary differential equations}, 
Encyclopaedia Math. Sci, Vol. 1, Springer, Berlin, 1986.

\bibitem{AN} A. Atabaigi, N. Nyamoradi, H. Zangeneh; 
\emph{The number of limit cycles of a quintic polynomial system with a
center}, Nonlinear Anal., 71 (2009), 3008-3017.

\bibitem{BL1} R. Benterki, J. Llibre; 
\emph{Limit cycles of polynomial differential equations with quintic
homogeneous nonlinearities}, J. Math. Anal. Appl., 407 (2013), 16-22.

\bibitem{BL2} A. Buic\u{a}, J. Llibre; 
\emph{Averaging methods for finding periodic orbits via Brouwer degree}, 
Bull. Sci. Math., 128 (2004), 7-22.

\bibitem{BL3} A. Buic\u{a}, J. Llibre; 
\emph{Limit cycles of a perturbed cubic polynominal
differential center}, Chaos Solitons Fractals, 32 (2007), 1059-1069.

\bibitem{BP} T. Blows, L. Perko; 
\emph{Bifurcation of limit cycles from centers and separatrix cycles of 
planar analytic systems}, SIAM Rev., 36 (1994), 341-376.

\bibitem{CGP} B. Coll, A. Gasull, R. Prohens; 
\emph{Bifurcation of limit cycles from two families of centers},
 Dyn. Contin. Discrete Impuls. Syst., Ser. A (Math. Anal.), 12 (2005), 275-287.

\bibitem{CJ} C. Chicone, M. Jacobs; \emph{Bifurcation of limit cycles
from quadratic isochrones}, J. Differential Equations, 91 (1991),
268-326.

\bibitem{CLP} B. Coll, J. Llibre, R. Prohens; 
\emph{Limit cycles bifurcating from a perturbed quartic center}, 
Chaos Solitons Fractals, 44 (2011), 317-334.

\bibitem{GL} J. Gin\'e, J. Llibre; 
\emph{Limit cycles of cubic polynomial vector fields via the averaging theory},
 Nonlinear Anal., 66 (2007), 1707-1721.

\bibitem{GLV1} H. Giacomini, J. Llibre, M. Viano; 
\emph{On the nonexistence, existence and uniqueness of limit cycles},
Nonlinearity, 9 (1996), 501-516.

\bibitem{GLV2} H. Giacomini, J. Llibre, M. Viano; 
\emph{On the shape of limit cycles that bifurcate from Hamiltonian centers}, 
Nonlinear Anal., 41 (2000), 523-537.

\bibitem{GLV3} H. Giacomini, J. Llibre, M. Viano; 
\emph{On the shape of limit cycles that bifurcate from non-Hamiltonian centers}, 
Nonlinear Anal., 43 (2001), 837-859.

 \bibitem{L} J. Llibre; 
\emph{Averaging theory and limit cycles for  quadratic systems}, 
Radovi Matemati\v{c}ki, 11 (2002), 1-14.

 \bibitem{LL} C. Li, J. Llibre; 
\emph{Quadratic perturbations of a  quadratic reversible Lotka-Volterra system}, 
Qual. Theory Dyn.  Syst., 9 (2010), 235-249.

 \bibitem{LLZ} C. Li, J. Llibre, Z. Zhang; 
\emph{Weak focus, limit cycles and bifurcations for bounded quadrtic
 systems}, J. Differential Equations, 115 (1995), 193-223.

 \bibitem{LPR}J. Llibre, J. P\'erez del R\'io, J. Rodr\'iguez; 
\emph{Averaging analysis of a perturbed quadratic  center}, 
Nonlinear Anal., 46 (2001), 45-51.

 \bibitem{VLG} M. Viano, J. Llibre, H. Giacomini; 
\emph{Arbitrary order bifurcations for perturbed Hamiltonian planar systems via the
reciprocal of an integrating factor}, Nonlinear Anal., 48 (2002),
117-136.

\bibitem{XH} G. Xiang, M. Han; 
\emph{Global bifurcation of limit cycles in a family of polynomial systems}, 
J. Math. Anal. Appl., 295 (2004), 633-644.

\end{thebibliography}

\end{document}