\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 111, pp. 1--27.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2015/111\hfil Limit cycles from a cubic reversible system] {Limit cycles from a cubic reversible system via the third-order averaging method} \author[L. Peng, Z. Feng \hfil EJDE-2015/111\hfilneg] {Linping Peng, Zhaosheng Feng} \address{Linping Peng (corresponding author) \newline School of Mathematics and System Sciences, Beihang University, LIMB of the Ministry of Education, Beijing 100191, China} \email{penglp@buaa.edu.cn Fax 86-(10)8231-7933} \address{Zhaosheng Feng \newline Department of Mathematics, University of Texas-Pan American\\ Edinburg, TX 78539, USA} \email{zsfeng@utpa.edu} \thanks{Submitted January 12, 2015. Published April 22, 2015.} \subjclass[2010]{34C07, 37G15, 34C05} \keywords{Bifurcation; limit cycles;homogeneous perturbation; \hfill\break\indent averaging method; cubic center; period annulus} \begin{abstract} This article concerns the bifurcation of limit cycles from a cubic integrable and non-Hamiltonian system. By using the averaging theory of the first and second orders, we show that under any small cubic homogeneous perturbation, at most two limit cycles bifurcate from the period annulus of the unperturbed system, and this upper bound is sharp. By using the averaging theory of the third order, we show that two is also the maximal number of limit cycles emerging from the period annulus of the unperturbed system. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction} \label{sect1} The Hilbert 16th problem proposes to find the maximal number of limit cycles of planar real polynomial differential equations $\dot{x}=f(x, y)$, $\dot{y}=g(x, y)$ in terms of the degree $n$ of the polynomials $f$ and $g$ \cite{H}. This is a longstanding problem, and many interesting and profound results have been established under various conditions. For example, the bifurcation of limit cycles from the periodic orbits around a center has been extensively studied in the literatures \cite{CLLZ, GI, GGI, I0, LL, LPR, L1} and the references therein. The weak Hilbert 16th problem in the quadratic Hamiltonian case was studied in \cite{CLLZ}. Simultaneously, quite a few 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, LL, LPR, L1} which is actually equivalent to the Abelian integrals in the plane. Although in the plane the method of Abelian integrals and the averaging theory are essentially equivalent, each has its own advantages. For example, when the associated Abelian integrals are complicated or we need to study the periodic orbits of the non-autonomous differential systems, the averaging method displays more flexibility. Roughly speaking, the averaging method gives a quantitative relation between the solutions of a non-autonomous periodic differential system and the solutions of its averaged differential equation, which is autonomous. Therefore, for some differential systems, the number of hyperbolic equilibrium points of their averaged differential equations can give a lower bound of the maximal number of limit cycles emerging from the periodic orbits around the center of the corresponding unperturbed system. As mentioned above, by using the averaging theory, the problem regarding the number of limit cycles of some differential systems can be reduced to the exploration of the number of hyperbolic equilibrium points of their averaged differential equations. Hence, the averaging theory has played a crucial role in the study of limit cycles of the differential systems and 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. However, it appears that these well-known results are focused on the first and second order bifurcations of limit cycles. As far as we know, analyzing the third order averaged function is generally very complicated, cumbersome and challenging, and even be out of the reach with the present stage of knowledge. Motivated by this fact and reference \cite{BL1}, in this paper we use the averaging theory of the first, second and third orders to study the bifurcation of limit cycles from the following cubic integrable and non-Hamiltonian system under any small cubic homogeneous perturbations \begin{equation}\label{eq1} \begin{gathered} \dot{x}=-y+x^{2}y,\\ \dot{y}=x+xy^{2}, \end{gathered} \end{equation} which has $$ H(x,y)=\frac{x^{2}+y^{2}}{1-x^{2}}=h $$ as its first integral with the integrating factor $2/(1-x^{2})^{2}$, and has the unique finite singularity $(0, 0)$ as its isochronous center. The period annulus, denoted by $$ \{(x, y)|H(x,y)=h, \quad h\in(0, +\infty)\} $$ starts at the center $(0, 0)$ and terminates at the unbounded separatrix formed by two invariant lines $x=\pm 1$ and the infinite degenerate singularities on the equator. The phase portrait of system \eqref{eq1} is shown in Figure \ref{fig1}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig1} \end{center} \caption{Phase portrait of system \eqref{eq1} in the Poincar\'{e} disk.} \label{fig1} \end{figure} We summarize our main results 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, 3; k=1, 2, 3)$, consider the following cubic homogeneous perturbation of system \eqref{eq1} \begin{equation}\label{eq2} \begin{gathered} \dot{x}=-y+x^{2}y+\sum_{k=1}^{3}\varepsilon^{k}\sum_{i+j=3}a_{ij}^{(k)}x^{i}y^{j},\\ \dot{y}=x+xy^{2}+\sum_{k=1}^{3}\varepsilon^{k}\sum_{i+j=3}b_{ij}^{(k)}x^{i}y^{j}. \end{gathered} \end{equation} Then the following two statements hold. \begin{enumerate} \item By using the averaging theory of first and second orders, system \eqref{eq2} has at most two limit cycles bifurcating from the period annulus around the center $(0, 0)$ of the unperturbed one, and in each case this upper bound is sharp. \item By using the averaging theory of third order, system \eqref{eq2} with $b_{03}^{(1)}$ being zero has at most two limit cycles bifurcating from the period annulus around the center $(0, 0)$ of the unperturbed one, and this upper bound is sharp. \end{enumerate} \end{theorem} The rest of this paper is organized as follows. In Section $2$, we give an introduction on the averaging theory of first, second and third orders, including some technical lemmas and methods employed in the averaging theory. Sections $3, 4$ and $5$ are dedicated to the study of the bifurcation of limit cycles by computing the first, second and third order averaged functions related to the equivalent system of system \eqref{eq2} and exploring the number of theirs simple zeros, respectively. 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, second and third orders, 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^{2}F_{2}(t, x)+\varepsilon^{3}F_{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 hypotheses $(i)$ and $(ii)$ 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}^{2}F_{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{align*} 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,\\ 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{align*} 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.\nonumber \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 (F_{1}^{0}+\varepsilon F_{2}^{0}+\varepsilon^{2} F_{3}^{0})(a)}{dx}\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{coro2.2} 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. If $F_{1}^{0}(x)$ and $F_{2}^{0}(x)$ are identically zero and $F_{3}^{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_{3}^{(0)}(x)$ for sufficiently small $|\varepsilon|$. In this case, conclusions in Lemma \ref{lem1} are true too. \end{corollary} \begin{remark} \rm To be convenient, we call the functions $F_{k}^{0}(x)~ (k=1, 2, 3)$, defined in Lemma \ref{lem1}, the first, second and third averaged functions associated with system \eqref{eq3}, respectively. \end{remark} Consider a 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 such continuous functions 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,\quad h_{c}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. In order 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(\rho(R, \varphi)\cos\varphi, \rho(R, \varphi)\sin\varphi)=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 \[ \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}, \] which is equivalent to \begin{align*} \frac{d R}{d\varphi} &= \Big\{\varepsilon\frac{\mu(x^{2}+y^{2})(Qp-Pq)}{2R(Qx-Py)}-\varepsilon^{2}\frac{\mu(x^{2}+y^{2})(Qp-Pq)(qx-py)}{2R(Qx-Py)^{2}}\\ &\quad +\varepsilon^{3}\frac{\mu(x^{2}+y^{2})(Qp-Pq)(qx-py)^{2}}{2R(Qx-Py)^{3}}\Big\} \Big|_{x=\rho(R, \varphi)\cos\varphi,\, y=\rho(R, \varphi)\sin\varphi}+O(\varepsilon^{4}),\nonumber \end{align*} where $P, Q, p$ and $q$ are defined as before. \end{lemma} \begin{remark} \rm It is notable that for the integrable and non-Hamiltonian systems, in general it is difficult to find the suitable transformations as described in Lemma \ref{lem2}. \end{remark} For $$ H(x, y)=\frac{x^{2}+y^{2}}{1-x^{2}}, $$ we choose the function $\rho=\rho(R, \varphi)$ as follows \begin{equation}\label{eq6} \rho(R, \varphi)=\frac{R}{\sqrt{1+R^{2}\cos^{2}\varphi}} \end{equation} such that $$ H(\rho(R, \varphi)\cos\varphi, \rho(R, \varphi)\sin\varphi)=R^{2}. $$ Applying Lemma \ref{lem2} to system \eqref{eq2}, we obtain the following result. \begin{lemma}\label{lem3} In 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+R^{2}\cos^{2}\varphi)^{2}}{R}\Big\{\varepsilon(Qp_{1}-Pq_{1})\\ &\quad +\varepsilon^{2}\Big[Qp_{2}-Pq_{2}-\frac{(Qp_{1}-Pq_{1})(xq_{1} -yp_{1})}{x^{2}+y^{2}}\Big] \\ &\quad +\varepsilon^{3}\Big[Qp_{3}-Pq_{3}-\frac{(Qp_{1}-Pq_{1})(xq_{2}-yp_{2}) +(Qp_{2}-Pq_{2})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\\ &\quad +\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})^{2}}{(x^{2}+y^{2})^{2}}\Big]\Big\}\Big|_{x=\rho(R, \varphi)\cos\varphi,\, y=\rho(R, \varphi)\sin\varphi}+O(\varepsilon^{4}), \end{split} \end{equation} where \begin{equation}\label{eq8} \begin{gathered} \begin{aligned} Qp_{k}-Pq_{k} &=a_{30}^{(k)}x^{4}+\Big(a_{21}^{(k)}+b_{30}^{(k)}\Big)x^{3}y +\Big(a_{12}^{(k)}+b_{21}^{(k)}\Big)x^{2}y^{2}\\ &\quad +\Big(a_{03}^{(k)}+b_{12}^{(k)}\Big)xy^{3}+b_{03}^{(k)}y^{4} -b_{30}^{(k)}x^{5}y+\Big(a_{30}^{(k)}-b_{21}^{(k)} \Big)x^{4}y^{2}\\ &\quad +\Big(a_{21}^{(k)}-b_{12}^{(k)}\Big)x^{3}y^{3} +\Big(a_{12}^{(k)}-b_{03}^{(k)}\Big)x^{2}y^{4}+a_{03}^{(k)}xy^{5}, \end{aligned}\\ \begin{aligned} xq_{k}-yp_{k} &=b_{30}^{(k)}x^{4}+\Big(b_{21}^{(k)}-a_{30}^{(k)}\Big)x^{3}y +\Big(b_{12}^{(k)}-a_{21}^{(k)}\Big)x^{2}y^{2}\\ &\quad +\Big(b_{03}^{(k)}-a_{12}^{(k)}\Big)xy^{3}-a_{03}^{(k)}y^{4}, \end{aligned} \end{gathered} \end{equation} and $a_{ij}^{(k)}$ and $b_{ij}^{(k)}$ $i, j=0, 1, 2, 3$; $k=1, 2, 3$ are real, and $\rho(R, \varphi)$ is given by \eqref{eq6}. \end{lemma} \section{Zeros of first order averaged functions} \label{sect3} \begin{proposition}\label{prop1} The first 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} The first order averaged equation corresponding to system \eqref{eq7} is \begin{equation}\label{eq9} \dot{R}=\varepsilon F_{1}^{0}(R), \end{equation} where \begin{equation} \begin{split} F_{1}^{0}(R)&= \frac{1}{2\pi}\int_{0}^{2\pi} \left[\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R}\Big(Qp_{1}-Pq_{1}\Big)\right] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi\\ &= \frac{1}{2\pi}\int_{0}^{2\pi}\Big\{R^{3} \left[a_{30}^{(1)}\cos^{4}\varphi+\Big(a_{12}^{(1)} +b_{21}^{(1)}\Big)\cos^{2}\varphi \sin^{2}\varphi+b_{03}^{(1)}\sin^{4}\varphi \right]\\ &\quad +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi} \Big[\Big(a_{30}^{(1)}-b_{21}^{(1)}\Big)\cos^{4}\varphi\sin^{2}\varphi\\ &\quad +\Big(a_{12}^{(1)}-b_{03}^{(1)}\Big)\cos^{2}\varphi\sin^{4}\varphi\Big] \Big\}d\varphi. \end{split} \label{eqq3.2} \end{equation} A straightforward calculation gives \begin{align*} &\int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{2}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi =\pi\big[\frac{1}{4R^{2}}-\frac{1}{R^{4}} -\frac{2}{R^{6}}+\big(\frac{2}{R^{4}}+\frac{2}{R^{6}}\big) \frac{1}{\sqrt{1+R^{2}}}\Big],\\ &\int_{0}^{2\pi}\frac{\cos^{2}\varphi\sin^{4}\varphi}{1+R^{2}\cos^{2} \varphi}d\varphi=\pi\big[\frac{3}{4R^{2}}+\frac{3}{R^{4}} +\frac{2}{R^{6}}-\big(\frac{2}{R^{2}}+\frac{4}{R^{4}}+\frac{2}{R^{6}}\big) \frac{1}{\sqrt{1+R^{2}}}\big]. \end{align*} Substituting the above formulas in \eqref{eqq3.2}, we find \begin{equation}\label{eq10} \begin{split} F_{1}^{0}(R) &= \frac{1}{2R}\Big\{\Big(a_{30}^{(1)}+a_{12}^{(1)}\Big)R^{4}+\Big(-a_{30}^{(1)}+3a_{12}^{(1)}+b_{21}^{(1)}-3b_{03}^{(1)}\Big)R^{2}\\ &\quad +2\Big(-a_{30}^{(1)}+a_{12}^{(1)}+b_{21}^{(1)}-b_{03}^{(1)}\Big) +\Big[-2\Big(a_{12}^{(1)}-b_{03}^{(1)}\Big)R^{4}\\ &\quad +2\Big(a_{30}^{(1)}-2a_{12}^{(1)}-b_{21}^{(1)}+2b_{03}^{(1)}\Big)R^{2}\\ &\quad +2\Big(a_{30}^{(1)}-a_{12}^{(1)}-b_{21}^{(1)}+b_{03}^{(1)}\Big)\Big] \frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{split} \end{equation} Recall that $\sqrt{1+R^{2}}>1$, and let $$ \sqrt{1+R^{2}}=\frac{1+w^{2}}{1-w^{2}} $$ for $00,\quad \frac{dG_{1}^{0}\big(R_{2}^{(1)}\big)}{dR}=-\frac{1}{832}<0. \] This completes the proof. \end{proof} \section{Zeros of second order averaged functions} \label{sect4} In this section, we consider the number of the zeros of second order averaged function associated with system \eqref{eq7}, in the case where the first order averaged function is $F_{1}^{0}(R)\equiv 0$. On the basis of formula \eqref{eq11}, we obtain the following result. \begin{lemma}\label{lem4} 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_{30}^{(1)}=-b_{03}^{(1)},\quad a_{12}^{(1)}=b_{03}^{(1)},\quad b_{21}^{(1)}=-b_{03}^{(1)}. \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{align*} F_{1}(R,\varphi) &=\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R}\Big[Qp_{1}-Pq_{1}\Big] \Big|_{x=\rho\cos\varphi, y=\rho\sin\varphi}\\ &=R^{3}\Big[-b_{03}^{(1)}\cos^{4}\varphi+\Big(a_{21}^{(1)}+b_{30}^{(1)} \Big)\cos^{3}\varphi\sin\varphi\\ &\quad +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{3}\varphi +b_{03}^{(1)}\sin^{4}\varphi\Big]\\ &\quad +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi} \Big[-b_{30}^{(1)}\cos^{5}\varphi\sin\varphi+\Big(a_{21}^{(1)}-b_{12}^{(1)} \Big)\cos^{3}\varphi\sin^{3}\varphi \\ &\quad +a_{03}^{(1)}\cos\varphi\sin^{5}\varphi\Big], \end{align*} \begin{align*} &F_{2}(R,\varphi)\\ &=\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R}\Big[Qp_{2}-Pq_{2} -\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\ &=\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[Qp_{2}-Pq_{2}\Big]\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\ &\quad +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi} \Big\{b_{30}^{(1)}b_{03}^{(1)}\cos^{8} \varphi -b_{30}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{7}\varphi\sin\varphi\\ &\quad +b_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \cos^{6}\varphi\sin^{2}\varphi -\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big) \Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\\ &\quad +b_{30}^{(1)}\Big(a_{03}^{(1)} +b_{12}^{(1)}\Big)\Big]\cos^{5}\varphi\sin^{3}\varphi -b_{03}^{(1)}\Big(a_{03}^{(1)}+b_{30}^{(1)}\Big)\cos^{4}\varphi\sin^{4}\varphi\\ &\quad +\left[a_{03}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big) -\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\right] \cos^{3}\varphi\sin^{5}\varphi \\ &\quad -b_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)} \Big)\cos^{2}\varphi\sin^{6}\varphi +a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \cos\varphi\sin^{7}\varphi+a_{03}^{(1)}b_{03}^{(1)}\sin^{8}\varphi\Big\}\\ &\quad +\frac{R^{7}}{(1+R^{2}\cos^{2}\varphi)^{2}}\Big\{\Big(b_{30}^{(1)} \Big)^{2}\cos^{9}\varphi\sin\varphi -2b_{30}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) \cos^{7}\varphi\sin^{3}\varphi\\ &\quad -\left(2a_{03}^{(1)}b_{30}^{(1)}-\Big(a_{21}^{(1)}-b_{12}^{(1)} \Big)^{2}\right)\cos^{5}\varphi\sin^{5}\varphi -2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \cos^{3}\varphi\sin^{7}\varphi\\ &\quad +\Big(a_{03}^{(1)}\Big)^{2}\cos\varphi\sin^{9}\varphi\Big\}, \end{align*} \begin{align*} y_{1}(R, \varphi) &=\int_{0}^{\varphi}F_{1}(R, \theta)d\theta\\ &=\int_{0}^{\varphi}R^{3}\Big[-b_{03}^{(1)}\cos^{4}\theta +\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{3}\theta\sin\theta\\ &\quad +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\theta\sin^{3}\theta +b_{03}^{(1)}\sin^{4}\theta\Big]d\theta\\ &\quad +R^{5}\Big[-b_{30}^{(1)}\int_{0}^{\varphi}\frac{\cos^{5}\theta\sin\theta}{1+R^{2}\cos^{2}\theta}d\theta+\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) \int_{0}^{\varphi}\frac{\cos^{3}\theta\sin^{3}\theta}{1+R^{2}\cos^{2}\theta}d\theta\\ &\quad +a_{03}^{(1)}\int_{0}^{\varphi}\frac{\cos\theta\sin^{5}\theta}{1+R^{2}\cos^{2} \theta}d\theta\Big], \end{align*} 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{lem5} The following integral equalities hold. \begin{gather*} \int_{0}^{\varphi}\frac{1}{1+R^{2}\cos^{2}\theta}d\cos^{2}\theta = \frac{1}{R^{2}}\ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \\ \int_{0}^{\varphi}\frac{\cos^{2}\theta}{1+R^{2}\cos^{2}\theta}d\cos^{2}\theta =-\frac{1}{R^{2}}+\frac{1}{R^{2}}\cos^{2}\varphi-\frac{1}{R^{4}} \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \\ \begin{aligned} \int_{0}^{\varphi}\frac{\cos^{4}\theta}{1+R^{2}\cos^{2}\theta}d\cos^{2}\theta &=-\frac{1}{2R^{2}}+\frac{1}{R^{4}}-\frac{1}{R^{4}}\cos^{2}\varphi +\frac{1}{2R^{2}}\cos^{4}\varphi\\ &\quad +\frac{1}{R^{6}} \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \end{aligned}\\ \begin{aligned} \int_{0}^{\varphi}\frac{\cos^{6}\theta}{1+R^{2}\cos^{2}\theta}d\cos^{2}\theta &=-\frac{1}{3R^{2}}+\frac{1}{2R^{4}}-\frac{1}{R^{6}} +\frac{1}{R^{6}}\cos^{2}\varphi-\frac{1}{2R^{4}}\cos^{4}\varphi\\ &\quad +\frac{1}{3R^{2}}\cos^{6}\varphi -\frac{1}{R^{8}}\ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big). \end{aligned} \end{gather*} \end{lemma} Based on Lemma \ref{lem5}, we obtain the following result. \begin{lemma}\label{lem6} The following integral equalities hold. \begin{gather*} \begin{aligned} \int_{0}^{\varphi}\frac{\cos^{5}\theta\sin\theta}{1+R^{2}\cos^{2}\theta}d\theta &= \frac{1}{4R^{2}}-\frac{1}{2R^{4}} +\frac{1}{2R^{4}}\cos^{2}\varphi-\frac{1}{4R^{2}}\cos^{4}\varphi\\ &\quad -\frac{1}{2R^{6}}\ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \end{aligned}\\ \begin{aligned} \int_{0}^{\varphi}\frac{\cos^{3}\theta\sin^{3}\theta}{1+R^{2}\cos^{2}\theta}d\theta &= \frac{1}{4R^{2}}+\frac{1}{2R^{4}} +\Big(-\frac{1}{2R^{2}}-\frac{1}{2R^{4}}\Big)\cos^{2}\varphi +\frac{1}{4R^{2}}\cos^{4}\varphi\\ &\quad +\Big(\frac{1}{2R^{4}}+\frac{1}{2R^{6}}\Big) \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \end{aligned}\\ \begin{aligned} \int_{0}^{\varphi}\frac{\cos\theta\sin^{5}\theta}{1+R^{2}\cos^{2}\theta}d\theta &= -\frac{3}{4R^{2}}-\frac{1}{2R^{4}} +\Big(\frac{1}{R^{2}}+\frac{1}{2R^{4}}\Big)\cos^{2}\varphi -\frac{1}{4R^{2}}\cos^{4}\varphi\\ &\quad +\Big(-\frac{1}{2R^{2}}-\frac{1}{R^{4}}-\frac{1}{2R^{6}}\Big) \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big). \end{aligned} \end{gather*} \end{lemma} Using Lemmas \ref{lem5} and \ref{lem6} and a straightforward computation, we have \begin{align*} &y_{1}(R,\varphi)\\ &= \frac{a_{21}^{(1)}-3a_{03}^{(1)}-b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2}R\\ &\quad +\Big[-b_{03}^{(1)}\sin\varphi\cos\varphi+\frac{a_{21}^{(1)} +b_{30}^{(1)}}{2}\sin^{2}\varphi +\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{4} \sin^{4}\varphi\Big]R^{3}\\ &\quad +\Big[\frac{-a_{21}^{(1)}+2a_{03}^{(1)}+b_{12}^{(1)}}{2}R^{3} +\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{2}R\Big] \cos^{2}\varphi\\ &\quad +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} \cos^{4}\varphi\\ &\quad +\Big[-\frac{a_{03}^{(1)}}{2}R^{3} +\frac{a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}}{2}R+\frac{a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}}{2R}\Big]\\ &\quad\times \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big). \end{align*} \begin{lemma}\label{lem7} The following integral equalities hold. \begin{gather*} \int_{0}^{2\pi}\frac{1}{1+R^{2}\cos^{2}\varphi}d\varphi =\frac{2\pi}{\sqrt{1+R^{2}}},\\ \int_{0}^{2\pi}\frac{\cos^{2}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi =\pi\Big[\frac{2}{R^{2}}-\frac{2}{R^{2}\sqrt{1+R^{2}}}\Big],\\ \int_{0}^{2\pi}\frac{\cos^{4}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi =\pi\Big[\frac{1}{R^{2}}-\frac{2}{R^{4}}+\frac{2}{R^{4}\sqrt{1+R^{2}}}\Big],\\ \int_{0}^{2\pi}\frac{\cos^{6}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi =\pi\Big[\frac{3}{4R^{2}}-\frac{1}{R^{4}}+\frac{2}{R^{6}} -\frac{2}{R^{6}\sqrt{1+R^{2}}}\Big],\\ \int_{0}^{2\pi}\frac{\cos^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi =\pi\Big[\frac{5}{8R^{2}}-\frac{3}{4R^{4}}+\frac{1}{R^{6}} -\frac{2}{R^{8}}+\frac{2}{R^{8}\sqrt{1+R^{2}}}\Big],\\ \int_{0}^{2\pi}(\cos^{4}\varphi-\sin^{4}\varphi) \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big)d\varphi =\pi\Big[\frac{4}{R^{2}}-\frac{4}{R^{2}}\sqrt{1+R^{2}}+2\Big]. \end{gather*} \end{lemma} \begin{proof} The first integral equalities can be obtained by a direct computation. Here we only show the derivation of the last integral. Let \begin{gather*} N_{1}(r) = \int_{0}^{2\pi}\cos^{4}\varphi \ln\Big(1-r^{2}\sin^{2}\varphi\Big)d\varphi,\\ N_{2}(r) = \int_{0}^{2\pi}\sin^{4}\varphi \ln\Big(1-r^{2}\sin^{2}\varphi\Big)d\varphi, \end{gather*} where $r^{2}=R^{2}/(1+R^{2})$. Since \begin{align*} N_{1}'(r)-N_{2}'(r) &= -2r\int_{0}^{2\pi}\frac{\sin^{2}\varphi}{1-r^{2}\sin^{2}\varphi}d\varphi +4r\int_{0}^{2\pi}\frac{\sin^{4}\varphi}{1-r^{2}\sin^{2}\varphi}d\varphi\\ &= \frac{4\pi r}{\sqrt{1-r^{2}}(1+\sqrt{1-r^{2}})^{2}}, \end{align*} and $N_{1}(0)=N_{2}(0)=0$, we get \[ N_{1}(r)-N_{2}(r)=\int_{0}^{r}\Big(N_{1}'(s)-N_{2}'(s)\Big)ds =4\pi\Big(\frac{1}{R^{2}}-\frac{1}{R^{2}}\sqrt{1+R^{2}}+\frac{1}{2}\Big). \] This enables us to arrive at the last integral equality. \end{proof} By Lemma \ref{lem7}, we obtain the following result. \begin{lemma}\label{lem8} The following integral equalities hold. \begin{gather*} \int_{0}^{2\pi}\frac{\cos^{6}\varphi\sin^{2}\varphi}{1+R^{2} \cos^{2}\varphi}d\varphi =\pi\Big[\frac{1}{8R^{2}}-\frac{1}{4R^{4}} +\frac{1}{R^{6}}+\frac{2}{R^{8}}+\Big(-\frac{2}{R^{6}} -\frac{2}{R^{8}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big], \\ \int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{4}\varphi}{1+R^{2} \cos^{2}\varphi}d\varphi=\pi\Big[\frac{1}{8R^{2}}-\frac{3}{4R^{4}} -\frac{3}{R^{6}}-\frac{2}{R^{8}} +\Big(\frac{2}{R^{4}}+\frac{4}{R^{6}}+\frac{2}{R^{8}}\Big) \frac{1}{\sqrt{1+R^{2}}}\Big], \\ \begin{aligned} &\int_{0}^{2\pi}\frac{\cos^{2}\varphi\sin^{6}\varphi}{1+R^{2}\cos^{2} \varphi}d\varphi\\ &=\pi\Big[\frac{5}{8R^{2}}+\frac{15}{4R^{4}} +\frac{5}{R^{6}}+\frac{2}{R^{8}} +\Big(-\frac{2}{R^{2}}-\frac{6}{R^{4}}-\frac{6}{R^{6}}-\frac{2}{R^{8}}\Big) \frac{1}{\sqrt{1+R^{2}}}\Big], \end{aligned}\\ \begin{aligned} &\int_{0}^{2\pi}\frac{\sin^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi\\ &=\pi\Big[-\frac{35}{8R^{2}}-\frac{35}{4R^{4}}-\frac{7}{R^{6}} -\frac{2}{R^{8}}+\Big(2+\frac{8}{R^{2}}+\frac{12}{R^{4}}+\frac{8}{R^{6}} +\frac{2}{R^{8}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big], \end{aligned}\\ \begin{aligned} &\int_{0}^{2\pi}\frac{\cos^{8}\varphi\sin^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}} d\varphi \\ &=\pi\Big[\frac{1}{8R^{4}}-\frac{1}{2R^{6}} +\frac{3}{R^{8}}+\frac{8}{R^{10}}+\Big(-\frac{7}{R^{8}} -\frac{8}{R^{10}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big], \end{aligned}\\ \begin{aligned} &\int_{0}^{2\pi}\frac{\cos^{6}\varphi\sin^{4}\varphi}{(1+R^{2} \cos^{2}\varphi)^{2}}d\varphi\\ &=\pi\Big[\frac{1}{8R^{4}}-\frac{3}{2R^{6}}-\frac{9}{R^{8}}- \frac{8}{R^{10}}+\Big(\frac{5}{R^{6}}+\frac{13}{R^{8}} +\frac{8}{R^{10}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big], \end{aligned}\\ \begin{aligned} &\int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{6}\varphi}{(1+R^{2} \cos^{2}\varphi)^{2}}d\varphi\\ &=\pi\Big[\frac{5}{8R^{4}}+\frac{15}{2R^{6}} +\frac{15}{R^{8}}+\frac{8}{R^{10}} +\Big(-\frac{3}{R^{4}}-\frac{14}{R^{6}}-\frac{19}{R^{8}} -\frac{8}{R^{10}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big], \end{aligned}\\ \int_{0}^{2\pi}\frac{\cos^{6}\varphi\sin^{2} \varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi =\pi\Big[\frac{1}{4R^{4}}-\frac{2}{R^{6}} -\frac{6}{R^{8}}+\Big(\frac{5}{R^{6}}+\frac{6}{R^{8}}\Big) \frac{1}{\sqrt{1+R^{2}}}\Big],\\ \int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{4}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}} d\varphi =\pi\Big[\frac{3}{4R^{4}}+\frac{6}{R^{6}} +\frac{6}{R^{8}} +\Big(-\frac{3}{R^{4}}-\frac{9}{R^{6}}-\frac{6}{R^{8}}\Big) \frac{1}{\sqrt{1+R^{2}}}\Big],\\ \begin{aligned} &\int_{0}^{2\pi}\frac{\cos^{2}\varphi\sin^{6} \varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\\ &=\pi\Big[-\frac{15}{4R^{4}}-\frac{10}{R^{6}} -\frac{6}{R^{8}}+\Big(\frac{1}{R^{2}}+\frac{8}{R^{4}} +\frac{13}{R^{6}}+\frac{6}{R^{8}}\Big)\frac{1}{\sqrt{1+R^{2}}}\Big]. \end{aligned} \end{gather*} \end{lemma} \begin{proposition}\label{prop2} 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).\nonumber \end{equation} \smallskip \noindent\textbf{Step 1:} Computation of the function $F_{21}^{0}(R)$. Let \begin{align*} A_{1} &= 3R^{2}\Big[-b_{03}^{(1)}\cos^{4}\varphi+\Big(a_{21}^{(1)} +b_{30}^{(1)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{3}\varphi\\ &\quad +b_{03}^{(1)}\sin^{4}\varphi\Big],\\ A_{2} &= \frac{5R^{4}+3R^{6}\cos^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big[-b_{30}^{(1)}\cos^{5}\varphi\sin\varphi +\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\cos^{3}\varphi\sin^{3}\varphi\\ &\quad +a_{03}^{(1)} \cos\varphi\sin^{5}\varphi\Big],\\ B_{1} &= R^{3}\Big[-b_{03}^{(1)}\sin\varphi\cos\varphi+\frac{a_{21}^{(1)} +b_{30}^{(1)}}{2}\sin^{2}\varphi +\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{4}\sin^{4}\varphi\Big], \\ B_{2} &= \frac{a_{21}^{(1)}-3a_{03}^{(1)}-b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2}R\\ &\quad +\Big[\frac{-a_{21}^{(1)}+2a_{03}^{(1)}+b_{12}^{(1)}}{2}R^{3} +\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{2}R\Big] \cos^{2}\varphi\\ &\quad +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} \cos^{4}\varphi\\ &\quad +\Big[-\frac{a_{03}^{(1)}}{2}R^{3}+\frac{a_{21}^{(1)}-2a_{03}^{(1)} -b_{12}^{(1)}}{2}R +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2R}\Big]\\ &\quad\times \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big). \end{align*} Then it gives $$ \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_{1}B_{1}+A_{1}B_{2} +A_{2}B_{1}+A_{2}B_{2})d\varphi. \end{equation} A direct calculation shows \begin{equation}\label{eq17} \int_{0}^{2\pi}A_{1}B_{1}d\varphi=0. \end{equation} Recalling that the function $A_{2}B_{2}$ is odd with respect to $\varphi$, we get \begin{equation}\label{eq18} \int_{0}^{2\pi}A_{2}B_{2}d\varphi=0. \end{equation} In addition, it is not difficult to verify that \begin{align} &\frac{1}{2\pi}\int_{0}^{2\pi}A_{1}B_{2}d\varphi \nonumber \\ &=\frac{1}{2\pi}\int_{0}^{2\pi} \Big\{\Big[\frac{3b_{03}^{(1)}\Big(a_{21}^{(1)}-3a_{03}^{(1)}-b_{30}^{(1)} -b_{12}^{(1)}\Big)}{4}R^{5} \nonumber\\ &\quad +\frac{3b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{2}R^{3}\Big] \Big(-\cos^{4}\varphi+\sin^{4}\varphi\Big) \nonumber\\ &\quad +\Big[\frac{3b_{03}^{(1)}\Big(-a_{21}^{(1)}+2a_{03}^{(1)} +b_{12}^{(1)}\Big)}{2}R^{5}+ \frac{3b_{03}^{(1)}\Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)} +b_{12}^{(1)}\Big)}{2}R^{3}\Big] \nonumber\\ &\quad \times \Big(-\cos^{6}\varphi+\cos^{2}\varphi\sin^{4}\varphi\Big) +\frac{3b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)} -b_{12}^{(1)}\Big)R^{5}}{4} \nonumber\\ &\quad\times \Big(-\cos^{8}\varphi +\cos^{4}\varphi\sin^{4}\varphi\Big) +\Big[\frac{3a_{03}^{(1)}b_{03}^{(1)}}{2}R^{5} -\frac{3b_{03}^{(1)}\Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big)}{2}R^{3} \nonumber\\ &\quad -\frac{3b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)} -b_{12}^{(1)}\Big)}{2}R\Big] \Big(\cos^{4}\varphi-\sin^{4}\varphi\Big) \nonumber\\ &\quad\times \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big)\Big\}d\varphi, \label{eq19} \end{align} and \begin{equation}\label{eq20} \begin{split} \frac{1}{2\pi}\int_{0}^{2\pi}A_{2}B_{1}d\varphi &=\frac{1}{2\pi}\Big\{5R^{7}\Big[b_{30}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi} \frac{\cos^{6}\varphi\sin^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\\ &\quad -b_{03}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\int_{0}^{2\pi} \frac{\cos^{4}\varphi\sin^{4}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\\ &\quad -a_{03}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi}\frac{\cos^{2} \varphi\sin^{6}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\Big]\\ &\quad +3R^{9}\Big[b_{30}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi} \frac{\cos^{8}\varphi\sin^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\\ &\quad -b_{03}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\int_{0}^{2\pi} \frac{\cos^{6}\varphi\sin^{4}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\\ &\quad -a_{03}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi}\frac{\cos^{4} \varphi\sin^{6}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}}d\varphi\Big]\Big\}. \end{split} \end{equation} Applying Lemmas \ref{lem7} and \ref{lem8} to \eqref{eq19} and \eqref{eq20} gives \begin{align} &\frac{1}{2\pi}\int_{0}^{2\pi}A_{1}B_{2}d\varphi \nonumber\\ &=\frac{1}{2}\Big\{\frac{3b_{03}^{(1)} \Big(a_{21}^{(1)}+5a_{03}^{(1)}-b_{30}^{(1)}-b_{12}^{(1)}\Big)}{8}R^{5} \nonumber\\ &\quad +\frac{3b_{03}^{(1)}\Big(-3a_{21}^{(1)}+15a_{03}^{(1)}+b_{30}^{(1)} +3b_{12}^{(1)}\Big)}{4}R^{3} \nonumber\\ &\quad +3b_{03}^{(1)}\Big(-3a_{21}^{(1)}+5a_{03}^{(1)}-b_{30}^{(1)}+3b_{12}^{(1)}\Big)R -\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R} \nonumber \\ &\quad +\Big[-6a_{03}^{(1)}b_{03}^{(1)}R^{3}+6b_{03}^{(1)} \Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big)R \nonumber\\ &\quad +\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big] \sqrt{1+R^{2}}\Big\}, \label{eq21} \end{align} and \begin{align} &\frac{1}{2\pi}\int_{0}^{2\pi}A_{2}B_{1}d\varphi \nonumber \\ &=\frac{1}{2}\Big\{\frac{3b_{03}^{(1)}\Big(-a_{21}^{(1)}-5a_{03}^{(1)} +b_{30}^{(1)}+b_{12}^{(1)}\Big)}{8}R^{5} \nonumber \\ &\quad +\frac{b_{03}^{(1)}\Big(3a_{21}^{(1)}-15a_{03}^{(1)}-b_{30}^{(1)} -3b_{12}^{(1)}\Big)}{4}R^{3} \nonumber\\ &\quad +b_{03}^{(1)}\Big(-3a_{21}^{(1)}+5a_{03}^{(1)}-b_{30}^{(1)} +3b_{12}^{(1)}\Big)R +\frac{6b_{03}^{(1)}\Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)} +b_{12}^{(1)}\Big)}{R} \nonumber\\ &\quad +\Big[4a_{03}^{(1)}b_{03}^{(1)}R^{5}+2a_{03}^{(1)}b_{03}^{(1)}R^{3} +2b_{03}^{(1)}\Big(3a_{21}^{(1)}-4a_{03}^{(1)}+2b_{30}^{(1)}-3b_{12}^{(1)}\Big)R \nonumber\\ &\quad +\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big]\frac{1}{\sqrt{1+R^{2}}}\Big\}. \label{eq22} \end{align} 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{1}{2}\Big\{\frac{b_{03}^{(1)}\Big(-3a_{21}^{(1)}+15a_{03}^{(1)} +b_{30}^{(1)}+3b_{12}^{(1)}\Big)}{2}R^{3}\\ &\quad +4b_{03}^{(1)}\Big(-3a_{21}^{(1)}+5a_{03}^{(1)}-b_{30}^{(1)}+3b_{12}^{(1)}\Big)R\\ &\quad +\frac{12b_{03}^{(1)}\Big(-a_{21}^{(1)} +a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}\Big)}{R}\\ &\quad +\Big[-6a_{03}^{(1)}b_{03}^{(1)}R^{3} +6b_{03}^{(1)}\Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big)R\\ &\quad +\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big]\sqrt{1+R^{2}}\\ &\quad +\Big[4a_{03}^{(1)}b_{03}^{(1)}R^{5}+2a_{03}^{(1)}b_{03}^{(1)}R^{3} +2b_{03}^{(1)}\Big(3a_{21}^{(1)}-4a_{03}^{(1)}+2b_{30}^{(1)}-3b_{12}^{(1)}\Big)R\\ &\quad +\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big]\frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{split} \end{equation} \smallskip \noindent\textbf{Step 2:} Computation of the Function $F_{22}^{0}(R)$. Similarly to Step 1, we have \begin{equation}\label{eq24} \begin{split} &\frac{1}{2\pi}\int_{0}^{2\pi}\Big[\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big(Qp_{2}-Pq_{2}\Big)\Big]\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi\\ &= \frac{1}{2R}\Big\{\Big(a_{30}^{(2)}+a_{12}^{(2)}\Big)R^{4}+\Big(-a_{30}^{(2)}+3a_{12}^{(2)}+b_{21}^{(2)}-3b_{03}^{(2)}\Big)R^{2}\\ &\quad +2\Big(-a_{30}^{(2)}+a_{12}^{(2)}+b_{21}^{(2)}-b_{03}^{(2)}\Big)+\Big[-2\Big(a_{12}^{(2)}-b_{03}^{(2)}\Big)R^{4}\\ &\quad +2\Big(a_{30}^{(2)}-2a_{12}^{(2)}-b_{21}^{(2)}+2b_{03}^{(2)}\Big)R^{2}\\ &\quad +2\Big(a_{30}^{(2)}-a_{12}^{(2)}-b_{21}^{(2)}+b_{03}^{(2)}\Big)\Big] \frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{split} \end{equation} Using Lemmas \ref{lem7} and \ref{lem8}, we deduce that \begin{align} &-\frac{1}{2\pi}\int_{0}^{2\pi}\Big[\frac{(1+R^{2} \cos^{2}\varphi)^{2}}{R}\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi \nonumber \\ &=\frac{R^{5}}{2\pi}\Big[b_{30}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi} \frac{\cos^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi +b_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \int_{0}^{2\pi}\frac{\cos^{6}\varphi\sin^{2}\varphi}{1+R^{2}\cos^{2} \varphi}d\varphi \nonumber\\ &\quad -b_{03}^{(1)}\left(a_{03}^{(1)}+b_{30}^{(1)}\right)\int_{0}^{2\pi} \frac{\cos^{4}\varphi\sin^{4}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi \nonumber\\ &\quad -b_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\int_{0}^{2\pi} \frac{\cos^{2}\varphi\sin^{6}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi \nonumber\\ &\quad +a_{03}^{(1)}b_{03}^{(1)}\int_{0}^{2\pi} \frac{\sin^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi\Big] \nonumber \\ &=\frac{1}{2}\Big\{\frac{b_{03}^{(1)}\Big(a_{21}^{(1)}-9a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{2}R^{3} +4b_{03}^{(1)}\Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big)R \nonumber \\ &\quad+\frac{4b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)} -b_{12}^{(1)}\Big)}{R} +\Big[2a_{03}^{(1)}b_{03}^{(1)}R^{5} \nonumber\\ &\quad -2b_{03}^{(1)} \Big(a_{21}^{(1)}-4a_{03}^{(1)}-b_{12}^{(1)}\Big)R^{3} -2b_{03}^{(1)}\Big(3a_{21}^{(1)}-5a_{03}^{(1)}+b_{30}^{(1)} -3b_{12}^{(1)}\Big)R \nonumber\\ &\quad -\frac{4b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big]\frac{1}{\sqrt{1+R^{2}}}\Big\}. \label{eq25} \end{align} It follows from \eqref{eq24} and \eqref{eq25} that \begin{align} F_{22}^{0}(R) &= \frac{1}{2}\Big\{\Big[\frac{b_{03}^{(1)}\Big(a_{21}^{(1)} -9a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{2}+a_{30}^{(2)}+a_{12}^{(2)} \Big]R^{3} \nonumber \\ &\quad +\Big[4b_{03}^{(1)}\Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big) -a_{30}^{(2)}+3a_{12}^{(2)}+b_{21}^{(2)}-3b_{03}^{(2)}\Big]R \nonumber \\ &\quad +\frac{4b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)} -b_{12}^{(1)}\Big) +2\Big(-a_{30}^{(2)}+a_{12}^{(2)}+b_{21}^{(2)}-b_{03}^{(2)}\Big)}{R} \nonumber\\ &\quad +\Big[2a_{03}^{(1)}b_{03}^{(1)}R^{5}-2\Big(b_{03}^{(1)} \Big(a_{21}^{(1)}-4a_{03}^{(1)}-b_{12}^{(1)}\Big) +a_{12}^{(2)}-b_{03}^{(2)}\Big)R^{3} \nonumber\\ &\quad +2\Big(b_{03}^{(1)}\Big(-3a_{21}^{(1)}+5a_{03}^{(1)}-b_{30}^{(1)} +3b_{12}^{(1)}\Big)+a_{30}^{(2)}-2a_{12}^{(2)}-b_{21}^{(2)}+2b_{03}^{(2)}\Big)R \nonumber \\ &\quad +\frac{4b_{03}^{(1)}\Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)} +b_{12}^{(1)}\Big)+2\Big(a_{30}^{(2)}-a_{12}^{(2)}-b_{21}^{(2)} +b_{03}^{(2)}\Big)}{R}\Big] \nonumber\\ &\quad\times \frac{1}{\sqrt{1+R^{2}}}\Big\}. \label{eq26} \end{align} Based on \eqref{eq23} and \eqref{eq26}, $F_{2}^{0}(R)$ can be re-expressed as \begin{align} F_{2}^{0}(R) &= \frac{1}{2}\Big\{\Big[b_{03}^{(1)}\Big(-a_{21}^{(1)}+3a_{03}^{(1)} +b_{30}^{(1)}+b_{12}^{(1)}\Big)+a_{30}^{(2)}+a_{12}^{(2)}\Big]R^{3} \nonumber \\ &\quad +\Big[4b_{03}^{(1)}\Big(-2a_{21}^{(1)}+3a_{03}^{(1)}-b_{30}^{(1)} +2b_{12}^{(1)}\Big)-a_{30}^{(2)}+3a_{12}^{(2)}+b_{21}^{(2)}-3b_{03}^{(2)}\Big]R \nonumber\\ &\quad +\frac{8b_{03}^{(1)}\Big(-a_{21}^{(1)}+a_{03}^{(1)} -b_{30}^{(1)}+b_{12}^{(1)}\Big)+2\Big(-a_{30}^{(2)}+a_{12}^{(2)}+b_{21}^{(2)} -b_{03}^{(2)}\Big)}{R} \nonumber\\ &\quad +\Big[-6a_{03}^{(1)}b_{03}^{(1)}R^{3}+6b_{03}^{(1)}\Big(a_{21}^{(1)} -2a_{03}^{(1)}-b_{12}^{(1)}\Big)R \nonumber\\ &\quad +\frac{6b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R}\Big] \sqrt{1+R^{2}} \nonumber\\ &\quad +\Big[6a_{03}^{(1)}b_{03}^{(1)}R^{5}+2\Big(b_{03}^{(1)} \Big(-a_{21}^{(1)}+5a_{03}^{(1)}+b_{12}^{(1)}\Big)-a_{12}^{(2)} +b_{03}^{(2)}\Big)R^{3} \nonumber\\ &\quad +2\Big(b_{03}^{(1)}\Big(a_{03}^{(1)}+b_{30}^{(1)}\Big) +a_{30}^{(2)}-2a_{12}^{(2)}-b_{21}^{(2)}+2b_{03}^{(2)}\Big)R \nonumber\\ &\quad +\frac{2b_{03}^{(1)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)+2\Big(a_{30}^{(2)}-a_{12}^{(2)} -b_{21}^{(2)}+b_{03}^{(2)}\Big)}{R}\Big] \nonumber\\ &\quad\times \frac{1}{\sqrt{1+R^{2}}}\Big\}. \label{eq27} \end{align} After making the transformations as before, \eqref{eq27} becomes \begin{equation}\label{eq28} \begin{split} F_{2}^{0}(R) &= \frac{w^{3}}{(1-w^{2})^{3}}\Big\{\Big[4b_{03}^{(1)} \Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}\Big) -a_{30}^{(2)}+a_{12}^{(2)} +b_{21}^{(2)} \\ &\quad -b_{03}^{(2)}\Big]w^{4} +\Big[8b_{03}^{(1)}\Big(a_{03}^{(1)}+b_{30}^{(1)}\Big) +2\Big(a_{30}^{(2)}+a_{12}^{(2)}-b_{21}^{(2)}-b_{03}^{(2)}\Big)\Big]w^{2}\\ &\quad +3a_{30}^{(2)}+a_{12}^{(2)}+b_{21}^{(2)}+3b_{03}^{(2)}\Big\}, \end{split} \end{equation} which shows that the second order averaged function $F_{2}^{0}(R)$ associated with system \eqref{eq7} has at most two zeros in $R\in(0, +\infty)$, by taking into account the multiplicity. Now we provide an example to demonstrate that this upper bound can be reached. Consider the system \begin{equation}\label{eq29} \begin{split} \dot{x} &= -y+x^{2}y+\varepsilon\Big[a_{21}^{(1)}x^{2}y+a_{03}^{(1)}y^{3}\Big] +\varepsilon^{2}\Big[-\frac{23}{4}x^{3}+a_{21}^{(2)}x^{2}y+\frac{19}{2}xy^{2} +a_{03}^{(2)}y^{3}\Big],\\ \dot{y} &= x+xy^{2}+\varepsilon\Big[b_{30}^{(1)}x^{3}+b_{12}^{(1)}xy^{2}\Big] +\varepsilon^{2}\Big[b_{30}^{(2)}x^{3}+\frac{35}{4}x^{2}y+b_{12}^{(2)}xy^{2}\Big], \end{split} \end{equation} where $a_{21}^{(k)}, a_{03}^{(k)}, b_{30}^{(k)}$ and $b_{12}^{(k)}~(k=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^{2}M_{2}(R, \varphi)+O(\varepsilon^{3}), \end{equation} where \begin{align*} M_{1}(R, \varphi) &= R^{3}\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{3}\varphi \sin\varphi+\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{3}\varphi\Big]\\ &\quad +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi}\Big[-b_{30}^{(1)} \cos^{5}\varphi\sin\varphi+\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) \cos^{3}\varphi\sin^{3}\varphi\\ &\quad +a_{03}^{(1)}\cos\varphi\sin^{5}\varphi\Big], \end{align*} \begin{align*} &M_{2}(R, \varphi)\\ &= R^{3}\Big[-\frac{23}{4}\cos^{4}\varphi+\Big(a_{21}^{(2)} +b_{30}^{(2)}\Big)\cos^{3}\varphi\sin\varphi +\frac{73}{4}\cos^{2}\varphi\sin^{2}\varphi \\ &\quad +\Big(a_{03}^{(2)}+b_{12}^{(2)}\Big)\cos\varphi\sin^{3}\varphi\Big] +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi}\Big\{-b_{30}^{(2)} \cos^{5}\varphi\sin\varphi \\ &\quad -\frac{29}{2}\cos^{4}\varphi\sin^{2}\varphi +\Big(a_{21}^{(2)}-b_{12}^{(2)}\Big)\cos^{3}\varphi\sin^{3}\varphi\\ &\quad +\frac{19}{2}\cos^{2}\varphi\sin^{4}\varphi+a_{03}^{(2)}\cos\varphi\sin^{5}\varphi -b_{30}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{7}\varphi\sin\varphi\\ &\quad -\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) + b_{30}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big]\cos^{5} \varphi\sin^{3}\varphi\\ &\quad +\Big[a_{03}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big) -\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big(b_{12}^{(1)} -a_{21}^{(1)}\Big)\Big]\cos^{3}\varphi\sin^{5}\varphi \\ &\quad +a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \cos\varphi\sin^{7}\varphi\Big\} +\frac{R^{7}}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big\{\Big(b_{30}^{(1)}\Big)^{2}\cos^{9}\varphi\sin\varphi\\ &\quad -2b_{30}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) \cos^{7}\varphi\sin^{3}\varphi -\Big[2a_{03}^{(1)}b_{30}^{(1)}-\Big(a_{21}^{(1)}-b_{12}^{(1)} \Big)^{2}\Big]\cos^{5}\varphi\sin^{5}\varphi\\ &\quad -2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \cos^{3}\varphi\sin^{7}\varphi+\Big(a_{03}^{(1)}\Big)^{2} \cos\varphi\sin^{9}\varphi\Big\}. \end{align*} It is not difficult to verify that for system \eqref{eq30}, the first order averaged function $M_{1}^{0}(R)\equiv 0$, while the second order averaged function $M_{2}^{0}(R)$ takes the form \begin{equation} \begin{split} M_{2}^{0}(R) &= \frac{1}{2\pi}\Big\{R^{3}\Big[-\frac{23}{4}\int_{0}^{2\pi}\cos^{4}\varphi d\varphi+\frac{73}{4}\int_{0}^{2\pi}\cos^{2}\varphi\sin^{2}\varphi d\varphi\Big]\\ &\quad +R^{5}\Big[-\frac{29}{2}\int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{2}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi +\frac{19}{2}\int_{0}^{2\pi}\frac{\cos^{2}\varphi\sin^{4}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi\Big]\Big\}\\ &= \frac{24w^{3}}{(1-w^{2})^{3}}\big(w-\frac{1}{4}\big) \big(w-\frac{1}{6}\big),\nonumber \end{split} \end{equation} 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{8}{15},\quad R_{2}^{(2)}=\frac{12}{35},\nonumber \end{equation} corresponding to $w_{1}^{(2)}=1/4$ and $w_{2}^{(2)}=1/6$, respectively, in $R\in(0, +\infty)$. Moreover, we have \[ \frac{dM_{2}^{0}\Big(R_{1}^{(2)}\Big)}{dR}=\frac{4}{255}>0, \quad \frac{dM_{2}^{0}\Big(R_{2}^{(2)}\Big)}{dR}=-\frac{6}{1295}<0. \] \end{proof} \section{Zeros of third order averaged function} \label{sect5} \begin{lemma}\label{lem9} For system \eqref{eq7}, the second order averaged function $F_{2}^{0}(R)\equiv 0$ (given by \eqref{eq28}) if and only if \begin{equation}\label{eq31} \begin{gathered} a_{30}^{(2)}=-b_{03}^{(2)}+b_{03}^{(1)}\Big(-a_{21}^{(1)} +a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}\Big),\\ a_{12}^{(2)}=b_{03}^{(2)}+2b_{03}^{(1)}\Big(a_{21}^{(1)} -2a_{03}^{(1)}-b_{12}^{(1)}\Big),\\ b_{21}^{(2)}=-b_{03}^{(2)}+b_{03}^{(1)}\Big(a_{21}^{(1)} +a_{03}^{(1)}+3b_{30}^{(1)}-b_{12}^{(1)}\Big). \end{gathered} \end{equation} \end{lemma} \begin{corollary}\label{coro1} Suppose $b_{03}^{(1)}=0$ in system \eqref{eq7}, then the second order averaged function $F_{2}^{0}(R)\equiv 0$ if and only if \begin{equation}\label{eq32} a_{30}^{(2)}=-b_{03}^{(2)},\quad a_{12}^{(2)}=b_{03}^{(2)},\quad b_{21}^{(2)}=-b_{03}^{(2)}. \end{equation} \end{corollary} Consider the third order averaged function $F_{3}^{0}(R)$ associated with system \eqref{eq7} in the case where conditions \eqref{eq14}, \eqref{eq32} and $b_{03}^{(1)}=0$ hold: \begin{equation}\label{eq33} F_{3}^{0}(R)=F_{31}^{0}(R)+F_{32}^{0}(R)+F_{33}^{0}(R)+F_{34}^{0}(R), \end{equation} where \begin{gather*} F_{31}^{0}(R)=\frac{1}{4\pi}\int_{0}^{2\pi}\frac{\partial^{2}F_{1}( R, \varphi)}{\partial R^{2}}y_{1}^{2}(R, \varphi)d\varphi,\\ F_{32}^{0}(R)=\frac{1}{4\pi}\int_{0}^{2\pi}\frac{\partial F_{1}(R, \varphi)}{\partial R}y_{2}(R, \varphi)d\varphi,\\ F_{33}^{0}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial F_{2}(R, \varphi)}{\partial R}y_{1}(R, \varphi)d\varphi,\\ F_{34}^{0}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}F_{3}(R, \varphi)d\varphi,\\ y_{2}(R, \varphi)=2\int_{0}^{\varphi}\Big[\frac{\partial F_{1}(R, \theta)}{\partial R}y_{1}(R, \theta)+F_{2}(R, \theta)\Big]d\theta, \end{gather*} \begin{align*} y_{1}(R, \varphi )&= \frac{a_{21}^{(1)}-3a_{03}^{(1)}-b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2}R\\ &\quad+\frac{a_{21}^{(1)}+b_{30}^{(1)}}{2}R^{3}\sin^{2}\varphi+\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{4}R^{3}\sin^{4}\varphi\\ &\quad+\Big[\frac{-a_{21}^{(1)}+2a_{03}^{(1)}+b_{12}^{(1)}}{2}R^{3} +\frac{-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}}{2}R\Big]\cos^{2}\varphi\\ &\quad+\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3}\cos^{4}\varphi\\ &\quad+\Big[-\frac{a_{03}^{(1)}}{2}R^{3} +\frac{a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}}{2}R+\frac{a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}}{2R}\Big]\\ &\quad\times \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big), \end{align*} \begin{align*} \frac{\partial F_{1}(R, \varphi)}{\partial R} &= 3R^{2}\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{3}\varphi\Big]\\ &\quad +\frac{5R^{4}+3R^{6}\cos^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big[-b_{30}^{(1)}\cos^{5}\varphi\sin\varphi\\ &\quad +\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\cos^{3}\varphi\sin^{3}\varphi +a_{03}^{(1)}\cos\varphi\sin^{5}\varphi\Big], \end{align*} \begin{align*} \frac{\partial^{2} F_{1}(R, \varphi)}{\partial R^{2}} &= 6R\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{3}\varphi\Big]\\ &\quad +\frac{2R^{3}(10+9R^{2}\cos^{2}\varphi+3R^{4}\cos^{4}\varphi)} {(1+R^{2}\cos^{2}\varphi)^{3}}\Big[-b_{30}^{(1)}\cos^{5}\varphi\sin\varphi\\ &\quad +\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\cos^{3}\varphi\sin^{3} \varphi+a_{03}^{(1)}\cos\varphi\sin^{5}\varphi\Big]. \end{align*} Here explicit expressions of $F_{2}(R, \varphi)$ and $F_{3}(R, \varphi)$ will be given in Steps 3 and 4 below, respectively. \begin{proposition}\label{prop3} Under conditions \eqref{eq14}, \eqref{eq32} and $b_{03}^{(1)}=0$, the third 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} To compute the third order averaged function $F_{3}^{0}(R)$, we divide our discussions into four steps. \smallskip \noindent\textbf{Step 1:} Computation of the function $F_{31}^{0}(R)$. Recall that $\partial^{2} F_{1}(R, \varphi)/\partial R^{2}$ is odd with respect to $\varphi$ while $y_{1}^{2}(R, \varphi)$ is even. Then we find \begin{equation}\label{eq34} F_{31}^{0}(R)=0. \end{equation} \smallskip \noindent\textbf{Step 2:} The Computation of the Function $F_{32}^{0}(R)$. According to Lemma \ref{lem1}, $y_{2}(R, \varphi)$ takes the form \begin{align*} y_{2}(R, \varphi) &=2\int_{0}^{\varphi}\Big(\frac{\partial F_{1}(R, \theta)}{\partial R}y_{1}(R, \theta)+F_{2}(R, \theta)\Big)d\theta\\ &=2\int_{0}^{\varphi}\frac{\partial F_{1}(R, \theta)}{\partial R} y_{1}(R, \theta)d\theta +2\int_{0}^{\varphi}\Big[\frac{(1+R^{2}\cos^{2}\theta)^{2}}{R} \Big(Qp_{2}-Pq_{2}\\ &\quad -\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big)\Big] \Big|_{x=\rho\cos\theta,\, y=\rho\sin\theta}d\theta. \end{align*} Since $\partial F_{1}(R, \theta)/{\partial R}$ is an odd function with respect to the variable $\theta$ and $y_{1}(R, \theta)$ is even, we obtain \[ \int_{0}^{2\pi}\frac{\partial F_{1}(R, \varphi)}{\partial \varphi}\Big(\int_{0}^{\varphi}\frac{\partial F_{1}(R, \theta)}{\partial \theta}y_{1}(R, \theta)d\theta\Big)d\varphi=0. \] Similar to the computation of $y_{1}(R, \varphi)$, we have \begin{align*} &\int_{0}^{\varphi}\frac{(1+R^{2}\cos^{2}\theta)^{2}}{R} \Big[Qp_{2}-Pq_{2}\Big]\Big|_{x=\rho\cos\theta,\, y=\rho\sin\theta}d\theta\\ &=\frac{a_{21}^{(2)}-3a_{03}^{(2)}-b_{30}^{(2)}-b_{12}^{(2)}}{4}R^{3} +\frac{a_{21}^{(2)}-a_{03}^{(2)}+b_{30}^{(2)}-b_{12}^{(2)}}{2}R\\ &\quad +\Big[-b_{03}^{(2)}\sin\varphi\cos\varphi +\frac{a_{21}^{(2)}+b_{30}^{(2)}}{2}\sin^{2}\varphi +\frac{-a_{21}^{(2)}+a_{03}^{(2)}-b_{30}^{(2)}+b_{12}^{(2)}}{4} \sin^{4}\varphi\Big]R^{3}\\ &\quad +\Big[\frac{-a_{21}^{(2)}+2a_{03}^{(2)}+b_{12}^{(2)}}{2}R^{3} +\frac{-a_{21}^{(2)}+a_{03}^{(2)}-b_{30}^{(2)}+b_{12}^{(2)}}{2}R\Big] \cos^{2}\varphi\\ &\quad +\frac{a_{21}^{(2)}-a_{03}^{(2)}+b_{30}^{(2)} -b_{12}^{(2)}}{4}R^{3}\cos^{4}\varphi\\ &\quad +\Big[-\frac{a_{03}^{(2)}}{2}R^{3} +\frac{a_{21}^{(2)}-2a_{03}^{(2)}-b_{12}^{(2)}}{2}R +\frac{a_{21}^{(2)}-a_{03}^{(2)} +b_{30}^{(2)}-b_{12}^{(2)}}{2R}\Big]\\ &\quad\times \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big). \end{align*} A straightforward computation yields \begin{equation}\label{eq35} \begin{split} &\int_{0}^{\varphi}-\frac{(1+R^{2}\cos^{2}\theta)^{2}}{R} \Big[\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\theta,\, y=\rho\sin\theta}d\theta\\ &=-R^{5}\Big\{b_{30}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)I_{7,1} +\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\Big(b_{12}^{(1)} -a_{21}^{(1)}\Big)\\ &\quad +b_{30}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big]I_{5,3} +\Big[-a_{03}^{(1)}\Big(a_{21}^{(1)} +b_{30}^{(1)}\Big)\\ &\quad +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\Big]I_{3,5} -a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)I_{1,7}\Big\}\\ &\quad -R^{7}\Big\{-\Big(b_{30}^{(1)}\Big)^{2}J_{9,1} +2b_{30}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)J_{7,3} +\Big[2a_{03}^{(1)}b_{30}^{(1)}\\ &\quad -\Big(a_{21}^{(1)} -b_{12}^{(1)}\Big)^{2}\Big]J_{5,5} +2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)J_{3,7} -\Big(a_{03}^{(1)}\Big)^{2}J_{1,9}\Big\}, \end{split} \end{equation} where \[ I_{k,l}=\int_{0}^{\varphi}\frac{\cos^{k}\theta\sin^{l} \theta}{1+R^{2}\cos^{2}\theta}d\theta, \] for $k, l=1, 3, 5, 7$ and \[ J_{k,l}=\int_{0}^{\varphi}\frac{\cos^{k}\theta\sin^{l}\theta} {(1+R^{2}\cos^{2}\theta)^{2}}d\theta, \] for $k, l=1, 3, 5, 7, 9$. Apparently, $I_{k,l}$ $(k, l=1, 3, 5, 7)$ and $J_{k,l}$ $(k, l=1, 3, 5,7, 9)$ are all even functions with respect to $\varphi$, so does the integral \eqref{eq35}. This fact together with the odd function $\partial F_{1}(R, \varphi)/\partial R$ leads to \begin{align*} &\int_{0}^{2\pi}\frac{\partial F_{1}(R, \varphi)}{\partial R}\Big(\int_{0}^{\varphi}-\frac{(1+R^{2}\cos^{2}\theta)^{2}}{R} \Big[\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{\substack{x=\rho\cos\theta,\\ y=\rho\sin\theta}} d\theta\Big)d\varphi\\ &=0. \end{align*} Hence, $F_{32}^{0}(R)$ can be simplified as \begin{align*} &F_{32}^{0}(R)\\ &= \frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial F_{1}(R, \varphi)}{\partial \varphi}\Big(\int_{0}^{\varphi} \frac{(1+R^{2}\cos^{2}\theta)^{2}}{R}\Big[Qp_{2} -Pq_{2}\Big]\Big|_{x=\rho\cos\theta, y=\rho\sin\theta}d\theta\Big)d\varphi\\ &= -\frac{b_{03}^{(2)}}{2\pi}R^{5}\Big\{\int_{0}^{2\pi}3 \Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{4}\varphi\sin^{2}\varphi +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos^{2}\varphi\sin^{4}\varphi\Big]d\varphi\\ &\quad +5R^{2}\Big[-b_{30}^{(1)}J_{6,2}+\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) J_{4,4}+a_{03}^{(1)}J_{2,6}\Big]\\ &\quad +3R^{4}\Big[-b_{30}^{(1)}J_{8,2}+\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) J_{6,4}+a_{03}^{(1)}J_{4,6}\Big]\Big\}, \end{align*} where \[ J_{k,l}=\int_{0}^{\varphi}\frac{\cos^{k}\theta\sin^{l} \theta}{(1+R^{2}\cos^{2}\theta)^{2}}d\theta, \] for $k=2, 4, 6, 8$; $l=2, 4, 6$. In view of Lemma \ref{lem8}, we have \begin{equation}\label{eq36} \begin{split} F_{32}^{0}(R) &= -\frac{b_{03}^{(2)}}{2}R\Big\{\frac{3\Big(a_{21}^{(1)}+3a_{03}^{(1)}\Big)}{4}R^{4} +\frac{-3a_{21}^{(1)}+15a_{03}^{(1)}+b_{30}^{(1)}+3b_{12}^{(1)}}{4}R^{2}\\ &\quad +3a_{21}^{(1)}-5a_{03}^{(1)}+b_{30}^{(1)}-3b_{12}^{(1)} +\frac{6\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R^{2}}\\ &\quad +\Big[-4a_{03}^{(1)}R^{4}-2a_{03}^{(1)}R^{2}+2\Big(-3a_{21}^{(1)} +4a_{03}^{(1)}-2b_{30}^{(1)}+3b_{12}^{(1)}\Big)\\ &\quad +\frac{6\Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)} +b_{12}^{(1)}\Big)}{R^{2}}\Big]\frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{split} \end{equation} \smallskip \noindent\textbf{Step 3:} Computation of the Function $F_{33}^{0}(R)$. Note that when conditions \eqref{eq14} and \eqref{eq32} hold, we derive that \begin{align*} &\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R}\Big[Qp_{2}-Pq_{2}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi} \\ &=R^{3}\Big[-b_{03}^{(2)}\cos^{4}\varphi +\Big(a_{21}^{(2)}+b_{30}^{(2)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(2)}+b_{12}^{(2)}\Big)\cos\varphi\sin^{3}\varphi\\ &\quad +b_{03}^{(2)}\sin^{4}\varphi\Big] +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi} \Big[-b_{30}^{(2)}\cos^{5}\varphi\sin\varphi\\ &\quad +\Big(a_{21}^{(2)}-b_{12}^{(2)}\Big)\cos^{3}\varphi\sin^{3}\varphi +a_{03}^{(2)}\cos\varphi\sin^{5}\varphi\Big]. \end{align*} A straightforward computation gives \begin{align*} &\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[-\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}} \Big]\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\ &=-\frac{R^{5}}{1+R^{2}\cos^{2}\varphi}\Big\{b_{30}^{(1)} \Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{7}\varphi\sin\varphi\\ &\quad +\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\Big(b_{12}^{(1)} -a_{21}^{(1)}\Big) +b_{30}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \Big]\cos^{5}\varphi\sin^{3}\varphi\\ &\quad +\Big[-a_{03}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big) +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\Big]\cos^{3}\varphi\sin^{5}\varphi\\ &\quad -a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \cos\varphi\sin^{7}\varphi\Big\} -\frac{R^{7}}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big\{-\Big(b_{30}^{(1)}\Big)^{2}\cos^{9}\varphi\sin\varphi\\ &\quad +2b_{30}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)\cos^{7}\varphi\sin^{3}\varphi\\ &\quad +\Big[2a_{03}^{(1)}b_{30}^{(1)}-\Big(a_{21}^{(1)} -b_{12}^{(1)}\Big)^{2}\Big]\cos^{5}\varphi\sin^{5}\varphi +2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\cos^{3}\varphi\sin^{7}\varphi\\ &\quad -\Big(a_{03}^{(1)}\Big)^{2}\cos\varphi\sin^{9}\varphi \Big\}. \end{align*} Hence, we have \begin{align*} & F_{2}(R,\varphi)\\ &= \frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R}\Big[Qp_{2}-Pq_{2} -\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})}{x^{2}+y^{2}} \Big]\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}\\ &= R^{3}\Big[-b_{03}^{(2)}\cos^{4}\varphi+\Big(a_{21}^{(2)} +b_{30}^{(2)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(2)}+b_{12}^{(2)}\Big)\cos\varphi\sin^{3}\varphi\\ &\quad +b_{03}^{(2)}\sin^{4}\varphi\Big] +\frac{R^{5}}{1+R^{2}\cos^{2}\varphi} \Big\{-b_{30}^{(2)}\cos^{5}\varphi\sin\varphi +\Big(a_{21}^{(2)}-b_{12}^{(2)}\Big)\cos^{3}\varphi\sin^{3}\varphi\\ &\quad +a_{03}^{(2)}\cos\varphi\sin^{5}\varphi -b_{30}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{7}\varphi\sin\varphi\\ &\quad -\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) +b_{30}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big]\cos^{5}\varphi\sin^{3}\varphi\\ &\quad +\Big[a_{03}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)-\Big(a_{03}^{(1)} +b_{12}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\Big] \cos^{3}\varphi\sin^{5}\varphi\\ &\quad +a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big) \cos\varphi\sin^{7}\varphi\Big\} +\frac{R^{7}}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big\{\Big(b_{30}^{(1)}\Big)^{2}\cos^{9}\varphi\sin\varphi\\ &\quad -2b_{30}^{(1)}\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big) \cos^{7}\varphi\sin^{3}\varphi -\Big[2a_{03}^{(1)}b_{30}^{(1)}-\Big(a_{21}^{(1)} -b_{12}^{(1)}\Big)^{2}\Big]\cos^{5}\varphi\sin^{5}\varphi\\ &\quad -2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \cos^{3}\varphi\sin^{7}\varphi +\Big(a_{03}^{(1)}\Big)^{2}\cos\varphi\sin^{9}\varphi\Big\}. \end{align*} Differentiating the function $F_{2}(R, \varphi)$ with respect to $R$ yields \begin{align*} &\frac{\partial F_{2}(R, \varphi)}{\partial R}\\ &= 3R^{2}\Big[-b_{03}^{(2)}\cos^{4}\varphi+\Big(a_{21}^{(2)} +b_{30}^{(2)}\Big)\cos^{3}\varphi\sin\varphi +\Big(a_{03}^{(2)}+b_{12}^{(2)}\Big)\cos\varphi\sin^{3}\varphi\\ &\quad +b_{03}^{(2)}\sin^{4}\varphi\Big] +\frac{5R^{4}+3R^{6}\cos^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{2}} \Big\{-b_{30}^{(2)}\cos^{5}\varphi\sin\varphi \\ &\quad +\Big(a_{21}^{(2)}-b_{12}^{(2)}\Big)\cos^{3}\varphi\sin^{3}\varphi +a_{03}^{(2)}\cos\varphi\sin^{5}\varphi -b_{30}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\cos^{7}\varphi\sin\varphi\\ &-\Big[\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) +b_{30}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big]\cos^{5}\varphi\sin^{3}\varphi\\ &\quad +\Big[a_{03}^{(1)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big) -\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\Big]\cos^{3}\varphi\sin^{5}\varphi\\ &\quad +a_{03}^{(1)}\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)\cos\varphi\sin^{7}\varphi\Big\}\\ &\quad +\frac{7R^{6}+3R^{8}\cos^{2}\varphi}{(1+R^{2}\cos^{2}\varphi)^{3}}\Big\{\Big(b_{30}^{(1)}\Big)^{2}\cos^{9}\varphi\sin\varphi -2b_{30}^{(1)}(a_{21}^{(1)}-b_{12}^{(1)})\cos^{7}\varphi\sin^{3}\varphi\\ &\quad -\Big[2a_{03}^{(1)}b_{30}^{(1)}-\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)^{2}\Big]\cos^{5}\varphi\sin^{5}\varphi -2a_{03}^{(1)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\cos^{3}\varphi\sin^{7}\varphi\\ &\quad +\Big(a_{03}^{(1)}\Big)^{2}\cos\varphi\sin^{9}\varphi \Big\}. \end{align*} Then %\begin{equation}\label{eq37} \begin{align} F_{33}^{0}(R) &=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial F_{2}(R, \varphi)}{\partial R}y_{1}(R, \varphi)d\varphi \nonumber \\ &=\frac{3R^{2}}{2\pi}\int_{0}^{2\pi}\Big\{b_{03}^{(2)} \Big[\frac{a_{21}^{(1)}-3a_{03}^{(1)}-b_{30}^{(1)}-b_{12}^{(1)}}{4}R^{3} +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2}R\Big] \nonumber\\ &\quad\times \Big(-\cos^{4}\varphi+\sin^{4}\varphi\Big) +\frac{b_{03}^{(2)}\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)R^{3}}{2} \Big(-\cos^{4}\varphi\sin^{2}\varphi+\sin^{6}\varphi\Big) \nonumber\\ &\quad +\frac{b_{03}^{(2)}\Big(-a_{21}^{(1)}+a_{03}^{(1)} -b_{30}^{(1)}+b_{12}^{(1)}\Big)R^{3}}{4} \Big(-\cos^{4}\varphi\sin^{4}\varphi+\sin^{8}\varphi\Big) \nonumber\\ &\quad +b_{03}^{(2)}\Big[\frac{-a_{21}^{(1)}+2a_{03}^{(1)} +b_{12}^{(1)}}{2}R^{3}+\frac{-a_{21}^{(1)}+a_{03}^{(1)} -b_{30}^{(1)}+b_{12}^{(1)}}{2}R\Big] \nonumber\\ &\quad\times \Big(-\cos^{6}\varphi+\cos^{2}\varphi\sin^{4}\varphi\Big) \nonumber\\ &\quad +\frac{b_{03}^{(2)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)R^{3}}{4}\Big(-\cos^{8}\varphi +\cos^{4}\varphi\sin^{4}\varphi\Big) \nonumber\\ &\quad +b_{03}^{(2)}\Big[-\frac{a_{03}^{(1)}}{2}R^{3} +\frac{a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}}{2}R +\frac{a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}}{2R}\Big] \nonumber\\ &\quad \times\Big(-\cos^{4}\varphi+\sin^{4}\varphi\Big) \ln\Big(1-\frac{R^{2}}{1+R^{2}}\sin^{2}\varphi\Big)\Big\}d\varphi. \label{eq37} \end{align} %\end{equation} By applying Lemma \ref{lem7}, the above equality becomes \begin{equation}\label{eq38} \begin{split} &F_{33}^{0}(R)\\ &= \frac{R}{2}\Big\{\frac{3b_{03}^{(2)}\Big(a_{21}^{(1)}+3a_{03}^{(1)}\Big)}{4}R^{4} +\frac{3b_{03}^{(2)}\Big(-3a_{21}^{(1)}+15a_{03}^{(1)} +b_{30}^{(1)}+3b_{12}^{(1)}\Big)}{4}R^{2}\\ &\quad +3b_{03}^{(2)}\Big(-3a_{21}^{(1)}+5a_{03}^{(1)} -b_{30}^{(1)}+3b_{12}^{(1)}\Big) -\frac{6b_{03}^{(2)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)} -b_{12}^{(1)}\Big)}{R^{2}}\\ &\quad +3b_{03}^{(2)}\Big[-2a_{03}^{(1)}R^{2}+2\Big(a_{21}^{(1)} -2a_{03}^{(1)}-b_{12}^{(1)}\Big)\\ &\quad +\frac{2\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R^{2}}\Big]\sqrt{1+R^{2}}\Big\}. \end{split} \end{equation} \smallskip \noindent\textbf{Step 4:} Computation of the function $F_{34}^{0}(R)$. Recall that \begin{align*} F_{3}(R, \varphi) &= \frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[Qp_{3}-Pq_{3}\\ &\quad -\frac{(Qp_{1}-Pq_{1})(xq_{2}-yp_{2}) +(Qp_{2}-Pq_{2})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\\ &\quad +\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})^{2}}{(x^{2}+y^{2})^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}. \end{align*} Similarly, we have \begin{align*} &\frac{1}{2\pi}\int_{0}^{2\pi}\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[Qp_{3}-Pq_{3}\Big]\Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi\\ &=\frac{1}{2R}\Big\{\Big(a_{30}^{(3)}+a_{12}^{(3)}\Big)R^{4} +\Big(-a_{30}^{(3)}+3a_{12}^{(3)}+b_{21}^{(3)}-3b_{03}^{(3)}\Big)R^{2}\\ &\quad +2\Big(-a_{30}^{(3)}+a_{12}^{(3)}+b_{21}^{(3)}-b_{03}^{(3)}\Big) +\Big[-2\Big(a_{12}^{(3)}-b_{03}^{(3)}\Big)R^{4}\\ &\quad +2\Big(a_{30}^{(3)}-2a_{12}^{(3)}-b_{21}^{(3)}+2b_{03}^{(3)}\Big)R^{2} +2\Big(a_{30}^{(3)}-a_{12}^{(3)}-b_{21}^{(3)}+b_{03}^{(3)}\Big)\Big] \frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{align*} Given the expressions \begin{gather*} \begin{aligned} Qp_{1}-Pq_{1}&=\Big(a_{21}^{(1)}+b_{30}^{(1)}\Big)x^{3}y +\Big(a_{03}^{(1)}+b_{12}^{(1)}\Big)xy^{3}-b_{30}^{(1)}x^{5}y\\ &\quad +\Big(a_{21}^{(1)}-b_{12}^{(1)}\Big)x^{3}y^{3}+a_{03}^{(1)}xy^{5}, \end{aligned}\\ xq_{k}-yp_{k}=b_{30}^{(k)}x^{4}+\Big(b_{12}^{(k)} -a_{21}^{(k)}\Big)x^{2}y^{2}-a_{03}^{(k)}y^{4},\quad k=1, 2,\nonumber \end{gather*} we assert that both $(Qp_{1}-Pq_{1})(xq_{2}-yp_{2})$ and $(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})^{2}$ are the odd functions with respect to $y$, which lead to \begin{gather*} \int_{0}^{2\pi}\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[\frac{(Qp_{1}-Pq_{1})(xq_{2}-yp_{2})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi=0,\\ \int_{0}^{2\pi}\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[\frac{(Qp_{1}-Pq_{1})(xq_{1}-yp_{1})^{2}}{(x^{2}+y^{2})^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi=0.\nonumber \end{gather*} Recall that \begin{gather*} \begin{aligned} Qp_{2}-Pq_{2} &=-b_{03}^{(2)}x^{4}+\Big(a_{21}^{(2)}+b_{30}^{(2)} \Big)x^{3}y+\Big(a_{03}^{(2)}+b_{12}^{(2)}\Big)xy^{3}+b_{03}^{(2)}y^{4}\\ &\quad -b_{30}^{(2)}x^{5}y+\Big(a_{21}^{(2)}-b_{12}^{(2)}\Big)x^{3}y^{3} +a_{03}^{(2)}xy^{5}, \end{aligned}\\ xq_{1}-yp_{1}=b_{30}^{(1)}x^{4}+\Big(b_{12}^{(1)} -a_{21}^{(1)}\Big)x^{2}y^{2}-a_{03}^{(1)}y^{4}. \end{gather*} Then \begin{align*} &-\frac{1}{2\pi}\int_{0}^{2\pi}\frac{(1+R^{2} \cos^{2}\varphi)^{2}}{R}\Big[\frac{(Qp_{2} -Pq_{2})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\varphi, y=\rho\sin\varphi}d\varphi\\ &=-\frac{R^{5}}{2\pi}\Big[-b_{30}^{(1)}b_{03}^{(2)} \int_{0}^{2\pi}\frac{\cos^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi -b_{03}^{(2)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big) \int_{0}^{2\pi}\frac{\cos^{6}\varphi\sin^{2}\varphi} {1+R^{2}\cos^{2}\varphi}d\varphi\\ &\quad +b_{03}^{(2)}\Big(a_{03}^{(1)}+b_{30}^{(1)}\Big) \int_{0}^{2\pi}\frac{\cos^{4}\varphi\sin^{4}\varphi} {1+R^{2}\cos^{2}\varphi}d\varphi\\ &\quad +b_{03}^{(2)}\Big(b_{12}^{(1)}-a_{21}^{(1)}\Big)\int_{0}^{2\pi} \frac{\cos^{2}\varphi\sin^{6}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi -a_{03}^{(1)}b_{03}^{(2)}\int_{0}^{2\pi} \frac{\sin^{8}\varphi}{1+R^{2}\cos^{2}\varphi}d\varphi\Big]. \end{align*} By using Lemmas \ref{lem7} and \ref{lem8}, the above function becomes \begin{align*} &-\frac{1}{2\pi}\int_{0}^{2\pi}\frac{(1+R^{2}\cos^{2}\varphi)^{2}}{R} \Big[\frac{(Qp_{2}-Pq_{2})(xq_{1}-yp_{1})}{x^{2}+y^{2}}\Big] \Big|_{x=\rho\cos\varphi,\, y=\rho\sin\varphi}d\varphi\\ &=-R\Big\{\frac{b_{03}^{(2)}\Big(-a_{21}^{(1)} +9a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}\Big)}{4}R^{2} +2b_{03}^{(2)}\Big(-a_{21}^{(1)}+2a_{03}^{(1)}+b_{12}^{(1)}\Big)\\ &\quad +\frac{2b_{03}^{(2)}\Big(-a_{21}^{(1)}+a_{03}^{(1)} -b_{30}^{(1)}+b_{12}^{(1)}\Big)}{R^{2}} \\ &\quad +\Big[-a_{03}^{(1)}b_{03}^{(2)}R^{4}+b_{03}^{(2)}\Big(a_{21}^{(1)} -4a_{03}^{(1)}-b_{12}^{(1)}\Big)R^{2}\\ &\quad +b_{03}^{(2)}\Big(3a_{21}^{(1)}-5a_{03}^{(1)}+b_{30}^{(1)}-3b_{12}^{(1)}\Big)\\ &\quad +\frac{2b_{03}^{(2)}\Big(a_{21}^{(1)}-a_{03}^{(1)} +b_{30}^{(1)}-b_{12}^{(1)}\Big)}{R^{2}}\Big]\frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{align*} Thus, we obtain \begin{equation}\label{eq39} \begin{split} &F_{34}^{(0)}(R)\\ &=\frac{1}{2\pi}\int_{0}^{2\pi}F_{3}(R, \varphi)d\varphi\\ &=\frac{1}{2R}\Big\{\Big[\frac{b_{03}^{(2)}\Big(a_{21}^{(1)}-9a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)}{2}+a_{30}^{(3)}+a_{12}^{(3)}\Big]R^{4}\\ &\quad +\Big[4b_{03}^{(2)}\Big(a_{21}^{(1)}-2a_{03}^{(1)}-b_{12}^{(1)}\Big)-a_{30}^{(3)}+3a_{12}^{(3)}+b_{21}^{(3)}-3b_{03}^{(3)}\Big]R^{2}\\ &\quad +\Big[4b_{03}^{(2)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big)+2\Big(-a_{30}^{(3)}+a_{12}^{(3)}+b_{21}^{(3)}-b_{03}^{(3)}\Big)\Big]\\ &\quad +\Big[2a_{03}^{(1)}b_{03}^{(2)}R^{6}+\Big(2b_{03}^{(2)}\Big(-a_{21}^{(1)}+4a_{03}^{(1)}+b_{12}^{(1)}\Big) +2\Big(-a_{12}^{(3)}+b_{03}^{(3)}\Big)\Big)R^{4}\\ &\quad +\Big(-2b_{03}^{(2)}\Big(3a_{21}^{(1)}-5a_{03}^{(1)}+b_{30}^{(1)}-3b_{12}^{(1)}\Big) +2\Big(a_{30}^{(3)}-2a_{12}^{(3)}-b_{21}^{(3)}+2b_{03}^{(3)}\Big)\Big)R^{2}\\ &\quad +\Big(-4b_{03}^{(2)}\Big(a_{21}^{(1)}-a_{03}^{(1)}+b_{30}^{(1)}-b_{12}^{(1)}\Big) +2\Big(a_{30}^{(3)}-a_{12}^{(3)}-b_{21}^{(3)}+b_{03}^{(3)}\Big)\Big)\Big] \frac{1}{\sqrt{1+R^{2}}}\Big\}. \end{split} \end{equation} Substituting \eqref{eq34}, \eqref{eq36}, \eqref{eq38} and \eqref{eq39} in \eqref{eq33}, and making the transformation $R=2w/(1-w^{2})$, we obtain \begin{align*} &F_{3}^{0}(R)\\ &= \frac{w^{3}}{(1-w^{2})^{3}}\Big\{\Big[4b_{03}^{(2)} \Big(-a_{21}^{(1)}+a_{03}^{(1)}-b_{30}^{(1)}+b_{12}^{(1)}\Big) -a_{30}^{(3)}+a_{12}^{(3)} +b_{21}^{(3)}-b_{03}^{(3)}\Big]w^{4}\\ &\quad +\Big[8b_{03}^{(2)}\Big(a_{03}^{(1)}+b_{30}^{(1)}\Big) +2\Big(a_{30}^{(3)}+a_{12}^{(3)}-b_{21}^{(3)}-b_{03}^{(3)}\Big)\Big]w^{2}\\ &\quad +3a_{30}^{(3)}+a_{12}^{(3)}+b_{21}^{(3)}+3b_{03}^{(3)}\Big\}, \end{align*} which has form similar to $F_{2}^{0}(R)$ given by \eqref{eq28}. Hence, $F_{3}^{0}(R)$ has at most two simple zeros in $R\in(0, +\infty)$, and this upper bound can be reached. For any sufficiently small $|\varepsilon|$, and any real constants $a_{21}^{(k)}, a_{03}^{(k)}, b_{30}^{(k)}$ and $b_{12}^{(k)}~(k=1, 2, 3)$, we take the following differential system as an example \begin{equation}\label{eq40} \begin{split} \dot{x} &= -y+x^{2}y+\varepsilon\Big[a_{21}^{(1)}x^{2}y+a_{03}^{(1)}y^{3}\Big] +\varepsilon^{2}\Big[a_{21}^{(2)}x^{2}y+a_{03}^{(2)}y^{3}\Big]\\ &\quad +\varepsilon^{3}\Big[-\frac{13}{2}x^{3}+a_{21}^{(3)}x^{2}y +\frac{21}{2}xy^{2}+a_{03}^{(3)}y^{3}\Big], \\ \dot{y}&= x+xy^{2}+\varepsilon\Big[b_{30}^{(1)}x^{3}+b_{12}^{(1)}xy^{2}\Big] +\varepsilon^{2}\Big[b_{30}^{(2)}x^{3}+b_{12}^{(2)}xy^{2}\Big]\\ &\quad +\varepsilon^{3}\Big[b_{30}^{(3)}x^{3}+10x^{2}y+b_{12}^{(3)}xy^{2}\Big]. \end{split} \end{equation} The third order averaged function corresponding to system \eqref{eq40} has exactly two simple zeros $R_{1}^{(3)}=3/4$ and $R_{2}^{(3)}=9/40$. \end{proof} Now we are in a position to prove Theorem \ref{thm1}. It follows immediately from Corollary \ref{coro2.2} and Propositions \ref{prop1}, \ref{prop2} and \ref{prop3} that system \eqref{eq7} has at most two periodic orbits, and there exist some systems which have exactly two periodic orbits shrinking to the corresponding hyperbolic equilibriums of their averaged equations, respectively. This implies that under the hypothesis of Theorem \ref{thm1}, system \eqref{eq2} has at most two limit cycles emerging from the period annulus around the center of the unperturbed system \eqref{eq2}$|_{\varepsilon=0}$, and the upper bound can be reached. \subsection*{Acknowledgments} This research is supported by the National Science Foundation of China under No. 11371046 and No. 11290141. \begin{thebibliography}{00} \bibitem{AI}V. I. Arnold, Y. S. Ilyashenko; \emph{Dynamical systems I: Ordinary differential equations, Encyclopaedia Math. Sci., Vol. 1}, Springer, Berlin, 1986. \bibitem{AN} A. Atabaigi, N. Nyamoradi, H. R. Z. 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. R. Blows, L. M. Perko; \emph{Bifurcation of limit cycles from centers and separatrix cycles of planar analytic systems}, SIAM Rev., 36 (1994), 341-376. \bibitem{CLLZ} F. D. Chen, C. Li, J. Llibre, Z. H. Zhang; \emph{A unified proof on the weak Hilbert 16th problem for $n=2$}, J. Differential Equations, 221 (2006), 309-342. \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{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{GI} L. Gavrilov, I. D. Iliev; \emph{Quadratic perturbations of quadratic codimension-four centers}, J. Math. Anal. Appl., 357 (2009), 69-76. \bibitem{GGI} S. Gautier, L. Gavrilov, I. D. Iliev; \emph{Perturbations of quadratic center of genus one}, Discrete Contin. Dyn. Syst., 25(2009), 511-535. \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{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{H} D. Hilbert; \emph{Mathematische probleme}, Arch. Math. Phys., 1(1901), 213-237. \bibitem{I0} I. D. Iliev; \emph{Perturbations of quadratic centers}, Bull. Sci. Math., 122 (1998), 107-161. \bibitem{LLZ} C. Li, J. Llibre, Z. Zhang; \emph{Weak focus, limit cycles and bifurcations for bounded quadartic systems}, J. Differential Equations, 115 (1995), 193-223. \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{LPR} J. Llibre, J. S. P\'erez del R\'{\i}o, J. A. Rodr\'{\i}guez; \emph{Averaging analysis of a perturbed quadratic center}, Nonlinear Anal., 46 (2001), 45-51. \bibitem{L1} J. Llibre; \emph{Averaging theory and limit cycles for quadratic systems}, Radovi Matemati$\breve{c}$ki, 11 (2002), 1-14. \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}