text:
# Theory of Current and Shot Noise Spectroscopy in Single-Molecular Quantum Dots with Phonon Mode

## Abstract

Using the Keldysh nonequilibrium Green function technique, we study the current and shot noise spectroscopy of a single molecular quantum dot coupled to a local phonon mode. It is found that in the presence of electron-phonon coupling, in addition to the resonant peak associated with the single level of the dot, satellite peaks with the separation set by the frequency of phonon mode appear in the differential conductance. In the “single level” resonant tunneling region, the differential shot noise power exhibit two split peaks. However, only single peaks show up in the “phonon assisted” resonant-tunneling region. An experimental setup to test these predictions is also proposed.

Recent progress in engineering and fabrication of nanoscale electronic devices has made it possible to study the transport properties of molecular devices Aviram98 ; Langlais99 ; Park00 , the linear dimension of which is at least an order smaller than semiconductor quantum dots. In addition to their potential industrial application, these devices provide an ideal test ground for the study of basic physics including the quantum size and many-body effects. The mechanism for electron conduction on such a small scale is not well understood yet. However, the evidence for quantum nature of transport properties has been observed in differential conductance for molecular wires by several experimental groups Andres96 ; Reed97 ; Kergueris99 ; Reichert02 ; Hong00 ; Rosink00 ; Chen00 ; Porath00 ; Smit02 . Theoretically, a lot of effort has been focused on the study of current-voltage characteristics of molecular wires based on either semi-empirical Tian98 ; Magoga99 ; Hall00 ; Paulsson01 ; Emberly00  or ab initio Ventra00 ; Taylor01 ; Palacios01 ; Damle02  methods. In contrast to semiconductor quantum dots, which is quite rigid in space, molecules involved in the electron tunneling process naturally possess the vibrational degrees of freedom which will inevitably react to the electron transfer through the molecules. So far, the influence of inelastic scattering process on the current-voltage characteristics of molecular wires has been considered Emberly00  while the influence on the current and its fluctuation (i.e., shot noise) has not been addressed for molecular dots. In this Letter, we use the Keldysh nonequilibrium Green function technique to calculate the current and shot noise through a single molecular quantum dot, for the first time. We focus on the effect of inelastic scattering process. A simplest Holstein-type model with a single phonon mode is employed to address the vibrational degrees of freedom in the molecular dot. All other complexity of real molecular devices, apart from interaction with phonon mode, is ignored. We find that in addition to the resonant peak associated with the single level of the dot, satellite peaks with the separation set by the frequency of phonon mode appear in the differential conductance. In the “single level” resonant tunneling region, the differential shot noise power exhibits two split peaks. However, only single peaks show up in the “phonon assisted” resonant-tunneling region. Therefore, due to the quantum nature of tunneling, the Fano factor is dramatically different from the Poisson limit even in the presence of inelastic process. These results may also be applied to the STM-based inelastic tunneling spectroscopy around a local vibrational mode on surfaces Stipe99 .

The model system under consideration is illustrated in Fig. 1. It consists of a quantum dot connected with two normal conducting leads. The electrons on the dot is also coupled to a single phonon mode. The single level of the dot is tuned by a gate voltage. The system Hamiltonian is written as:

$$H=H_L+H_R+H_X+H_D+H_T.$$
(1)
The first two terms are respectively the Hamiltonian for electrons in the left and right non-interacting metallic leads:

$$H_L+H_R=\underset{kL,R;\sigma }{}ϵ_kc_{k\sigma }^{}c_{k\sigma },$$
(2)
where we have denoted the electron creation (annihilation) operators in the leads by $`c_{k\sigma }^{}`$ ($`c_{k\sigma }`$). The quantity $`k`$ is the momentum and $`\sigma `$ is the spin index and $`ϵ_k`$ is the single particle energy of conduction electrons. The third term describes the phonon mode:

$$H_X=\omega _0a^{}a,$$
(3)
where $`\omega _0`$ is the frequency of the single phonon mode, $`a^{}`$ ($`a`$) is phonon creation (annihilation) operator. The fourth term describes the electron in the quantum dot:

$$H_D=\underset{\sigma }{}[ϵ_0+\lambda (a+a^{})]d_\alpha ^{}d_\alpha ,$$
(4)
where $`d_\alpha ^{}`$ ($`d_\alpha `$) are the creation (annihilation) operators of dot electrons, $`ϵ_0`$ is the single energy level of the dot, $`\lambda `$ is the coupling constant between the dot electrons and the phonon mode. The last term represents the coupling of the dot to the leads:

$$H_T=\underset{kL,R;\sigma ,\alpha }{}[V_{k\sigma ,\alpha }c_{k\sigma }^{}d_\alpha +H.c.],$$
(5)
where the tunneling matrix elements $`V_{k\sigma ,\alpha }`$ transfer electrons through an insulating barrier out of the dot.

The current operators from the leads to the dot are, respectively, defined as:

$$\widehat{I}_{L(R)}=ie\underset{kL(R);\sigma ,\alpha }{}(V_{k\sigma ,\alpha }c_{k\sigma }^{}d_\alpha H.c.).$$
(6)
The terminal current is given by $`I=(\widehat{I}_L\widehat{I}_R)/2`$, where $`\mathrm{}`$ denotes the statistical average of physical observables. Using the Keldysh nonequilibrium Green function formalism Caroli71 ; Keldysh65 , the terminal current can be obtained as Meir92 ; Jauho94 :

$`I`$ $`=`$ $`{\displaystyle \frac{ie}{4\pi }}{\displaystyle }dϵ\{\text{Tr}[(f_L𝚪^Lf_R𝚪^R)(𝐆^r𝐆^a)]`$ (7)
$`+\text{Tr}[(𝚪^L𝚪^R)𝐆^<]\},`$
where $`f_{L(R)}`$ are the Fermi distribution function of the left and right leads, which has different chemical potential upon a voltage bias $`\mu _L\mu _R=eV`$, the coupling of the dot to the leads is characterized by the parameter:

$$\mathrm{\Gamma }_{\alpha \alpha ^{}}^{L(R)}=2\pi \underset{\sigma }{}\rho _{L(R),\sigma }(ϵ)V_{\sigma ,\alpha }^{}(ϵ)V_{\sigma ,\alpha ^{}}(ϵ)$$
(8)
with $`\rho _{L(R),\sigma }`$ the spin-$`\sigma `$ band density of states in the two leads, $`G_{\alpha \alpha ^{}}^{r(a)}`$ and $`G_{\alpha \alpha ^{}}^<`$, are the Fourier transform of the dot electron retarded (advanced), $`G_{\alpha \alpha ^{}}^{r(a)}(t,t^{})=i\theta (\pm tt^{})\{d_\alpha (t),d_\alpha ^{}^{}(t^{})\}`$, and lesser Green function, $`G_{\alpha \alpha ^{}}^<(t,t^{})=id_\alpha ^{}^{}(t^{})d_\alpha (t)`$. Note that the boldface notation indicates that the level-width function $`𝚪`$ and the dot electron Green functions are matrices in the spin space of dot electrons.

The electrical current through a device is fluctuating in time even under a dc bias Blanter00 . Its fluctuation is characterized by the spectral density, which is given by the Fourier transform $`S(\omega )`$ of the current correlation function Note1 :

$$S(t,t^{})=\frac{1}{2}\{\mathrm{\Delta }\widehat{I}_L(t),\mathrm{\Delta }\widehat{I}_L(t^{})\},$$
(9)
where $`\mathrm{\Delta }\widehat{I}_L(t)=\widehat{I}_L(t)\widehat{I}_L(t)`$. The quantum statistical (non-equilibrium) average involved in the current correlation can be evaluated conveniently in the Keldysh Green function formalism. Following a tedious but standard procedure, we obtain the spectral density of shot noise in the zero-frequency limit:

$`S(\omega 0)={\displaystyle \frac{e^2}{2\pi }}{\displaystyle }dϵ\{f_L(1f_L)(\text{Tr}[(𝚪^L𝐆^r)^2]`$
$`+\text{Tr}[(𝚪^L𝐆^a)^2])+if_L\text{Tr}[𝚪^L𝐆^>]i(1f_L)\text{Tr}[𝚪^L𝐆^<]`$
$`+f_L\text{Tr}[𝚪^L𝐆^>𝚪^L(𝐆^r𝐆^a)]+\text{Tr}[𝚪^L𝐆^>𝚪^L𝐆^<]`$
$`(1f_L)\text{Tr}[𝚪^L(𝐆^r𝐆^a)𝚪^L𝐆^<]\},`$ (10)
where $`G_{\alpha \alpha ^{}}^>`$ is the Fourier transform of the greater Green function, $`G_{\alpha \alpha ^{}}^>(t,t^{})=id_\alpha (t)d_\alpha ^{}^{}(t^{})`$. Equation (10) is a basic formula for the low-frequency shot noise through a quantum dot, which can take into account the many-body effects conveniently. It expresses the fluctuations of current through the quantum dot, an interacting region, in terms of the distribution functions in the leads and local properties of the quantum dot, such as the occupation and density of states. This formalism can be viewed as a generalized version of the two-terminal shot-noise formula for the non-interacting case Chen91 ; Buttiker92 , which will become clear in the following discussion. It can also be used to study the current fluctuation in the case of spin-dependent transport, which is of much interest in the study of spintronics and single spin detection.

Once the dot electron Green functions are known, the electrical current and shot noise can be calculated using Eqs. (7) and (10). In the following, we calculate these Green functions, which should be carried out in the presence of the leads. By performing a canonical transformation, $`\overline{H}=e^𝒮He^𝒮`$ with $`𝒮=(\lambda /\omega _0)_\alpha d_\alpha ^{}d_\alpha (a^{}a)`$ Mahan00 , one can obtain $`\overline{H}=\overline{H}_{el}+\overline{H}_{ph}`$. Here $`\overline{H}_{el}=\overline{H}_L+\overline{H}_R+\overline{H}_D+\overline{H}_T`$ and $`\overline{H}_{ph}=\overline{H}_X`$, where $`\overline{H}_{L(R)}=H_{L(R)}`$, $`\overline{H}_D=_\alpha (ϵ_0\lambda ^2/\omega _0)d_\alpha ^{}d_\alpha `$, $`\overline{H}_TH_T`$, and $`\overline{H}_X=H_X`$. Consequently, the original dot-electron Green function can be decoupled as:

$`G_{\alpha \alpha ^{}}^{r(a)}(t)`$ $`=`$ $`i\theta (\pm t)\{\stackrel{~}{d}_\alpha (t),\stackrel{~}{d}_\alpha ^{}^{}(0)\}_{el}X(t)X^{}(0)_{ph}`$ (11)
$`=`$ $`\stackrel{~}{G}_{\alpha \alpha ^{}}^{r(a)}(t)X(t)X^{}(0)_{ph},`$
where $`\stackrel{~}{d}_\alpha (t)=e^{i\overline{H}_{el}t}d_\alpha e^{i\overline{H}_{el}t}`$, and $`X(t)=e^{i\overline{H}_{ph}t}Xe^{i\overline{H}_{ph}t}`$ with $`X=\mathrm{exp}[(\lambda /\omega _0)(a^{}a)]`$. The renormalization factor due to the electron-phonon interaction is evaluated to be Mahan00 : $`X(t)X^{}(0)_{ph}=e^{\mathrm{\Phi }(t)}`$, where $`\mathrm{\Phi }(t)=(\lambda /\omega _0)^2[N_{ph}(1e^{i\omega _0t})+(N_{ph}+1)(1e^{i\omega _0t})]`$ with $`N_{ph}=[\mathrm{exp}(\beta \omega _0)1]^1`$. We then apply the equation-of-motion approach to compute the Green function $`\stackrel{~}{G}_{\alpha \alpha ^{}}^{r(a)}`$. It consists of differentiating the Green function with respect to time, thereby generating higher-order Green functions which eventually have to close (by truncation in the presence of electron-electron interaction) the equation for the Green function $`\stackrel{~}{G}^{r(a)}`$. A little algebra gives rise to the Fourier transform:

$$\stackrel{~}{G}_{\alpha \alpha ^{}}^{r(a)}(\omega )=\frac{\delta _{\alpha \alpha ^{}}}{\omega (ϵ_0\mathrm{\Delta })\mathrm{\Sigma }_{\alpha \alpha }^{r(a)}},$$
(12)
where the energy shift $`\mathrm{\Delta }=\lambda ^2/\omega _0`$ and the retarded (advanced) self-energy due to the tunneling into the electrical leads are given by:

$$\mathrm{\Sigma }_{\alpha \alpha ^{}}^{r(a)}(\omega )=\underset{kL,R}{}\frac{|V_{k\alpha ,\alpha }|^2\delta _{\alpha \alpha ^{}}}{\omega ϵ_k\pm i0^+}.$$
(13)
The full width of the resonance is then just the sum of elastic couplings to the two leads which is given by Eq. (8). The Fourier transform of the full Green function given by Eq. (11) can be obtained as:

$`G_{\alpha \alpha ^{}}^{r(a)}`$ $`=`$ $`e^{g(2N_{ph}+1)}{\displaystyle \underset{l=\mathrm{}}{\overset{\mathrm{}}{}}}J_l(2g\sqrt{N_{ph}(N_{ph}+1)})`$ (14)
$`\times {\displaystyle \frac{\delta _{\alpha \alpha ^{}}e^{l\omega _0\beta /2}}{\omega (ϵ_0\mathrm{\Delta })l\omega _0\mathrm{\Sigma }_{\alpha \alpha }^{r(a)}}},`$
where $`J_l(x)`$ are the Bessel functions of complex argument, the parameter $`g=(\lambda /\omega _0)^2`$. The Green function $`G_{\alpha \alpha ^{}}^{<(>)}`$ cannot be obtained by directly using the above equation-of-motion approach. By instead applying the operational rules as given by Langreth Langreth76  to the Dyson equation for the contour-ordered Green function, one can show the following Keldysh equation for the lesser and greater functions Jauho94 : $`𝐆^{<(>)}(\omega )=𝐆^r(\omega )𝚺_T^{<(>)}(\omega )𝐆^a(\omega ).`$ In the weak electron-phonon coupling limit as we are considering here, the contribution to the self-energy $`𝚺_T^{<(>)}`$ from the electron-phonon interaction is negligible such that $`𝚺_T^{<(>)}`$ can be approximated by the lesser (greater) self-energy due to the tunneling into the two leads: $`𝚺^<=i[f_L𝚪^L+f_R𝚪^R],`$ and $`𝚺^>=i[(1f_L)𝚪^L+(1f_R)𝚪^R].`$ This relation enables us to rewrite the current and shot noise formulae as:

$$I=\frac{e}{2\pi }𝑑ϵ[f_Lf_R]\text{Tr}[𝐓],$$
(15)
and

$`S`$ $`=`$ $`{\displaystyle \frac{e^2}{2\pi }}{\displaystyle }dϵ\{[f_L(1f_L)+f_R(1f_R)]\text{Tr}[𝐓]`$ (16)
$`+(f_Lf_R)^2\text{Tr}[(1𝐓)𝐓]\},`$
where we have defined the transmission coefficient matrix $`𝐓=𝐆^a𝚪^L𝐆^r𝚪^R`$. We remark that for the noninteracting case (no electron-phonon and electron-electron interactions), the above two expressions are exact. These are just the Landauer-Büttiker formalisms developed for the noninteracting electron transport based on the scattering matrix theory. The connection between the two formalisms for the current was first established by Meir and Wingreen Meir92 .

For simplicity, we consider the coupling of the dot to the two leads to be symmetric $`𝚪^𝐋=𝚪^R`$, and assume that the leads have broad and flat density of states (i.e., the wideband limit). For this case, the elastic couplings, from which the self-energy (13) can be determined, are independent of energy. In Figs. 2 and 3, we plot the zero-temperature differential conductance and shot noise power as a function of the Fermi energy (that is, the equilibrium chemical potential in the leads) measured relative to the single level $`ϵ_0`$ in the presence (with typical value $`\lambda =0.4\omega _0`$, red-solid lines in the figures) and absence ($`\lambda =0`$, blue-dashed lines) of electron-phonon coupling. Here all energy quantities are measured in units of the phonon mode frequency $`\omega _0`$, and the elastic couplings $`\mathrm{\Gamma }^L=\mathrm{\Gamma }^R`$ are fixed at $`0.2\omega _0`$. When the electron is not coupled to the phonon mode, only one resonant conductance peak shows up as the Fermi energy matches the single level in the quantum dot. In the presence of electron-phonon coupling, the overall spectrum is shifted by $`\mathrm{\Delta }=\lambda ^2/\omega _0`$. In addition to the main peak related to the single level, new satellite resonant peaks appear at the positive energy side. The separation between the conductance peaks is set by the frequency of the phonon mode $`\omega _0`$. Since at zero temperature no phonon modes are excited on the quantum dot, the electrons tunneling onto the dot can only excite phonon modes, which explains why the satellite peaks are located at the positive energy region. It can be expected that there will appear satellite peaks at negative energy region at finite temperatures, where thermally excited phonons can be absorbed by the electrons tunneling onto the dot. Moreover, the intensity of the satellite peaks is much smaller than the main resonant peak because they are evolving from the emission of phonon modes, which is controlled by the electron-phonon coupling. Correspondingly, the differential shot noise power exhibits two peaks located symmetrically around the position where the conductance peak associated with the single level of the dot is located. The origin of this behavior is that no noise is generated when the transmission probability $`T=0`$ or $`1`$, while it is generated maximally in between. In this “single level” resonant-tunneling region, the Fano factor defined as the ratio of shot noise to current is significantly enhanced. Since the transmission probability corresponding to those satellite conductance peaks (from the phonon emission process) is small, the differential noise power is approximately proportional to the conductance and only show a single peak rather than two peaks. Therefore, in the “phonon assisted” resonant-tunneling region, the Fano factor is smaller than one, in contrast to the case of “single level” resonant tunneling. Therefore, due to the quantum tunneling nature, the Fano factor for the present system differs the Poisson limit (equal to 2), where electrons diffuse in an uncorrelated way.

To experimentally investigate the current and shot noise spectroscopy of molecular quantum dots, we suggest a setup where a gate electrode is attached to the dot capacitively, as shown in Fig. 1. A gate voltage is applied solely to tune the level position on the dot while a bias voltage is applied in the linear response regime. The merit of the proposed setup is that the internal electronic and ionic charge buildup on the dot, which is generated by a nonlinear response, could be avoided safely. Experimental effort along this direction is in progress Park02 .

To summarize, using the Keldysh nonequilibrium Green function technique, we have studied in this work, to the best of our knowledge, for the first time the current and shot noise spectroscopy of a single molecular quantum dot coupled to a local phonon mode. We show that, in the presence of electron-phonon coupling, in addition to the resonant peak associated with the single level of the dot, satellite peaks with the separation set by the frequency of phonon mode appear in the differential conductance. In the “single level” resonant tunneling region, the differential shot noise power exhibits two split peaks. However, only single peaks show up in the “phonon assisted” resonant-tunneling region. The Fano factor is different from the Poisson limit even in the presence of inelastic process. The derived current and noise formalism can also be applied to more complicated systems.

Acknowledgments: We thank A. Abanov and R. Lu for useful discussions. This work was supported by the Department of Energy.


format: text

text:
# Number and amplitude of limit cycles emerging from topologically equivalent perturbed centers

## 1 Introduction

Dissipation is present in almost every physical phenomenon. The only possibility to create self-sustained oscillations in a system is by re-injecting the energy dissipated into it. For instance, the pumping in a switched on laser must provide the necessary energy to compensate the lost energy ought to the interaction between the single frequency light beam and its surrounding (cavity, mirrors, etc.). In other context, the solar nuclear bomb must radiate to the space enough energy to maintain all the living cycles on the earth. In this case, many different frequencies are associated to the thousands of different species of animals and plants . Thus, every oscillation in nature is generated by a balance between input and output energies, in definitive between amplification and dissipation. This dynamical state (characterized by a preferred period, wave form and amplitude) is represented in phase space by a isolated closed curve called limit cycle . This isolated periodic motion can emerge from the perturbation of a conservative system where a continuum of periodic solutions exists. Only those closed orbits where the energetic balance after perturbation vanishes remain as periodic motions. Take, for instance, the harmonic oscillator $`\ddot{x}+x=0`$. A pattern of circles filling the whole plane $`(x,\dot{x})`$ are all the solutions of this system. A weak nonlinear perturbation proportional to the velocity, such as that in the Liénard equation $`\ddot{x}+ϵf(x)\dot{x}+x=0`$, destroys the geometrical aspect of the phase space. The system can balance the pumping and damping caused by the nonlinear term, $`ϵf(x)\dot{x}`$, only on some special periodic orbits. In the van der Pol case, where $`f(x)=(1x^2)`$, it has been proved the existence, uniqueness and non-algebraicity of a limit cycle for every value of the control parameter $`ϵ`$ . Depending on the conditions imposed over the strength of the nonlinearity $`ϵ`$ and over the properties of $`f(x)`$: degree, parity, etc., many results on the number and amplitude of the limit cycles for general Liénard systems are scattered in the literature (see, for instance, Ref.  and references therein). This behavior is the mathematical evidence of how a small perturbation in a physical system can destroy the coexistence of infinitely many periodic motions and leaves only a finite number of them verifying strict energetic balance conditions.

Following this line of reasoning, the aim of this work is to exploit the topological similarities of different systems presenting an infinite number of periodic oscillations which give rise to the emergence of isolated closed orbits under small perturbations. This is an alternative perspective of the problem of identifying the limit cycles growing from centers respect to other methods existing in the literature . In Section 2 we sketch the theory behind this type of systems and we develop the method to find those periodic solutions. In section 3 we apply this method to three different problems and calculate the number and amplitude of the limit cycles emerging from their perturbations. The comparison with previous investigations on these systems is also performed. Finally, we present our conclusions.

## 2 Geometrical and Topological Equivalence

The plane is the natural space for the representation of the integral curves of a second order one-dimensional or first order two-dimensional dynamical system. A periodic motion of this system draws a closed orbit on the plane. If the relation $`H(x,y)=h`$, with $`h`$ a fixed real number, is verified by every point $`(x,y)`$ on this type of orbit, a continuum set of periodic orbits is obtained when $`h`$ runs on an interval of real values, $`h[h_1,h_2]`$. This orbit structure receives the name of period annulus. Suppose, without lose of generality, that the boundary orbit given by $`H=h_1`$ is a degenerate orbit, say the origin $`(0,0)`$. This equilibrium point surrounded in its immediate neighborhood by closed paths is called a center . Obviously, every orbit of the period annulus verifies the differential relation

$$dH=0.$$
(1)
Depending on the time parameterization of the curves $`H(x,y)=h`$, different dynamical systems having the same pattern of integral paths are obtained. At a first sight, the difficulty to find the integrating factor in each case prevent us from establishing the phase plane geometrical equivalence between them: we say that two systems are phase plane geometrically equivalent when their integral curves on the plane are the same. Think, for instance, on the set of circles $`H_c(x,y)=\frac{1}{2}(x^2+y^2)=h`$. They obey the differential relation $`dH_c=2xdx+2ydy=0`$. Different time parameterizations of this expression produce systems with different time behaviors (e.g.: (1) $`\ddot{x}+x=0`$ ; (2) $`\dot{x}=ky`$, $`\dot{y}=kx`$, $`k=cte.`$, and (3) $`\dot{x}=y^2`$, $`\dot{y}=xy`$), although with the circles as their integral curves. If the interest does not reside on the time evolution of the system, the last concept of geometrical equivalence will take importance in the study of the system.

One step further. It is also possible that two different planar systems (having each one a period annulus, given by $`H(x,y)=h`$ and $`\overline{H}(x,y)=\overline{h}`$, with $`h[h_1,h_2]`$ and $`\overline{h}[\overline{h}_1,\overline{h}_2]`$, respectively) which are not phase plane geometrical equivalent, are yet topologically equivalent. By this we mean that there exists a bijective and continuous transformation of coordinates (homomorphism), $`\mathrm{\Gamma }:𝐑^2𝐑^2`$, transforming one period annulus into the other one. That is, if we represent the integral curves of the first system as $`C_h=\{(x,y)H(x,y)=h\}`$ and of the second one as $`\overline{C}_{\overline{h}}=\{(x,y)\overline{H}(x,y)=\overline{h}\}`$ then $`\mathrm{\Gamma }(\overline{C}_{\overline{h}})=C_h`$ in a continuous and monotone way when $`h`$ and $`\overline{h}`$ run in their intervals of existence. In particular, the boundary conditions over $`\mathrm{\Gamma }`$ are: $`\mathrm{\Gamma }(\overline{C}_{\overline{h}_1})=C_{h_1}`$ and $`\mathrm{\Gamma }(\overline{C}_{\overline{h}_2})=C_{h_2}`$. Take for example the system $`\dot{x}=ay`$, $`\dot{y}=bx`$ with $`a,b>0`$. It presents a period annulus formed by ellipses, $`H(x,y)=\frac{x^2}{a}+\frac{y^2}{b}=h`$. Clearly, this system is topologically equivalent to another one having a circular period annulus: $`\overline{H}_c(x,y)=\frac{1}{2}(x^2+y^2)=\overline{h}`$. It is straightforward to verify that the function $`\mathrm{\Gamma }(x,y)=(\frac{x}{\sqrt{a}},\frac{y}{\sqrt{b}})`$ establishes the homomorphism between both systems. Obviously, if two systems are geometrical equivalent they are topologically equivalent because in this case $`H=\overline{H}`$ and then we can choose $`\mathrm{\Gamma }`$ as the identity mapping.

The relevance of this property is the following. Imagine that we are able to find the limit cycles that emerge from the perturbation of the first system,

$$dH+\text{perturbation}=0.$$
(2)
Then, the application of $`\mathrm{\Gamma }`$ to the second perturbed system, $`d\overline{H}+\text{perturbation’}=0`$, bring this last expression to the form (2), even if we do not know $`\overline{H}`$. Now the limit cycles of the second system can be found by solving equation (2). Undoing the change of variables with $`\mathrm{\Gamma }^1`$, we will obtain the periodic motions of the perturbed system in the original coordinates.

## 3 Limit Cycles emerging from some Centers

### 3.1 Perturbation of a circular period annulus

As an application of the method above explained, we study here the limit cycles emerging from the perturbation of what we could consider a paradigmatic case in the center’s typology: a continuum set of circles, $`H_c(x,y)=\frac{1}{2}(x^2+y^2)=h`$. The solution of many different dynamical systems is this kind of period annulus. As a simple representation of them we take the linear center:

$$\begin{array}{c}\dot{x}=y,\hfill \\ \dot{y}=x,\hfill \end{array}$$
(3)
whose integral curves clearly verify $`dH_c=x\dot{x}+y\dot{y}=0`$. This path diagram is destroyed when we slightly perturb this conservative system with two general nonlinear terms, $`f(x,y)`$ and $`g(x,y)`$, controlled by the parameter $`ϵ1`$,

$$\begin{array}{c}\dot{x}=y+ϵf(x,y),\hfill \\ \dot{y}=x+ϵg(x,y).\hfill \end{array}$$
(4)
It is a time parameterized version of the differential relation: $`dH_c+ϵ(fdygdx)=0`$, in which the time variable has been removed. Following the procedure proposed in the previous section, we need to study this differential relation in order to find the limit cycles of system (4). If its integral curves are represented by $`y(x)`$ and we define $`y^{}(x)=dy/dx`$, we obtain:

$$(yy^{}+x)+ϵ[f(x,y)y^{}g(x,y)]=0.$$
(5)
If we suppose the origin $`(0,0)`$ is the only fixed point of equation (4), a limit cycle $`C_l(x,y_\pm (x))`$ of it, with a positive branch $`y_+(x)>0`$ and a negative branch $`y_{}(x)<0`$, cut the $`x`$-axis in two points $`(a_{},0)`$ and $`(a_+,0)`$ with $`a_{},a_+>0`$. Every limit cycle $`C_l`$ solution of equation (4) encloses the origin and the oscillation $`x`$ runs in the interval $`a_{}<x<a_+`$.

The amplitudes of oscillation $`a_{},a_+`$ identify the limit cycle. The result for $`ϵ0`$ is a nested set of closed curves that defines the qualitative distribution of the integral curves in the plane $`(x,y)`$. The stability of the limit cycles is alternated. For a given stable limit cycle, the two neighboring limit cycles, the closest one in its interior and the closest one in its exterior, are unstable, and conversely (see Figure 1). If we determine the stability of the origin, then the stability of every limit cycle remains fixed by the alternation property. The stability of the origin is determined by the sign of the real part of the eigenvalues of the Jacobian at this point. Hence, its calculation shows that if $`ϵ(f_x(0,0)+g_y(0,0))<0`$ then the origin is stable and if $`ϵ(f_x(0,0)+g_y(0,0))>0`$ then it is unstable.

Another property of a limit cycle can be derived from the fact that the mechanical energy $`E=H_c=\frac{1}{2}(x^2+y^2)/2`$ is conserved in a whole oscillation:

$$_{C_l}𝑑H_c=_{C_l}\frac{dE}{dx}𝑑x=0.$$

Thus, if equation (5) is integrated along a limit cycle, between the amplitudes of oscillation, we obtain:

$$\begin{array}{cc}_a_{}^{a_+}[g(x,y_+(x))f(x,y_+(x))y_+^{}(x)]𝑑x+& \\ _{a_+}^a_{}[g(x,y_{}(x))f(x,y_{}(x))y_{}^{}(x)]𝑑x=0.& \end{array}$$
(6)
Every couple of solutions $`\{y_+(x),y_{}(x)\}`$ of equation (6), vanishing in the extremes, $`y_+(a_{})=y_{}(a_{})=y_+(a_+)=y_{}(a_+)=0`$, and verifying equation (5), constitute the finite set of limit cycles of equation (4). As a consequence of the continuity in the parameter $`ϵ`$, these conditions are also valid for $`ϵ=0`$. For this parameter value, the limit cycles are circles: $`y_+(x)=\sqrt{a^2x^2}`$ and $`y_{}(x)=\sqrt{a^2x^2}`$. Here $`a_+=a_{}=a>0`$ represents the amplitudes of the limit cycles. At this order of approximation, the condition (5) is verified, and condition (6) reads:

$$\begin{array}{c}_a^a\left[g(x,\sqrt{a^2x^2})+\frac{xf(x,\sqrt{a^2x^2})}{\sqrt{a^2x^2}}\right]𝑑x+\\ _a^a\left[g(x,\sqrt{a^2x^2})\frac{xf(x,\sqrt{a^2x^2})}{\sqrt{a^2x^2}}\right]𝑑x=0.\end{array}$$
(7)
that is,

$$\beta (a)_a^a\left[\stackrel{~}{g}(x,\sqrt{a^2x^2})+x\frac{\stackrel{~}{f}(x,\sqrt{a^2x^2})}{\sqrt{a^2x^2}}\right]𝑑x=0,$$
(8)
where

$$\begin{array}{c}\stackrel{~}{f}(x,y)\frac{1}{2}[f(x,y)+f(x,y)f(x,y)f(x,y)],\\ \stackrel{~}{g}(x,y)\frac{1}{2}[g(x,y)g(x,y)+g(x,y)g(x,y)].\end{array}$$
(9)
If $`\beta (a)`$ does not vanish identically, each solution $`a>0`$ of the equation $`\beta (a)=0`$ is the amplitude of a limit cycle of the system (4) in the weak nonlinear regime. And conversely, at order zero in $`ϵ`$, the amplitudes of all the limit cycles emerging from the period annulus are solutions of equation (8). These results are exact for $`ϵ=0`$. Therefore, equations (7) or (8) determine the amplitudes of the period motions surviving to a slight perturbation of a period annulus formed by circles. A typical portrait of limit cycles in this regime is given in Figure 1.

We remark that $`\beta (a)`$ is the first Melnikov function of the system (4). Here it has been obtained by following a different approach to the usual one (based on the calculation of the displacement function of the first return mapping on the period annulus). Observe also that $`\stackrel{~}{f}(x,y)`$ is an odd function of $`x`$ and an even function of $`y`$, whereas $`\stackrel{~}{g}(x,y)`$ is an even function of $`x`$ and an odd function of $`y`$. Therefore, if $`f(x,y)`$ and $`g(x,y)`$ are polynomials in $`x`$ and $`y`$, only the odd terms in $`x`$ and even in $`y`$ of $`f(x,y)`$ and the even terms in $`x`$ and odd in $`y`$ of $`g(x,y)`$ survive in (9) and contribute to $`\beta (a)`$ in (8).

As an example, we integrate equation (8) when $`f(x,y)`$ is a polynomial of degree $`2n_1`$ or $`2n_1+1`$ in $`x`$ and degree $`2m_1`$ or $`2m_1+1`$ in $`y`$ and $`g(x,y)`$ is a polynomial of degree $`2n_2`$ or $`2n_2+1`$ in $`x`$ and degree $`2m_2`$ or $`2m_2+1`$ in $`y`$. Then,

$$\begin{array}{c}\stackrel{~}{f}(x,y)=_{j=0}^{n_1}_{k=0}^{m_1}a_{j,k}x^{2j+1}y^{2k},\\ \stackrel{~}{g}(x,y)=_{j=0}^{n_2}_{k=0}^{m_2}b_{j,k}x^{2j}y^{2k+1}.\end{array}$$

where $`a_{j,k}`$ and $`b_{j,k}`$ are real coefficients. The result is

$$\beta (a)=\frac{a^2}{2}\underset{j=0}{\overset{n}{}}\underset{k=0}{\overset{m}{}}\frac{\mathrm{\Gamma }(j+1/2)\mathrm{\Gamma }(k+1/2)}{(k+j+1)!}\left[\left(j+\frac{1}{2}\right)a_{j,k}+\left(k+\frac{1}{2}\right)b_{j,k}\right]a^{2(k+j)},$$

where $`n`$Max$`\{n_1,n_2\}`$ and $`m`$Max$`\{m_1,m_2\}`$. Here, $`a_{j,k}=0`$ for $`j>n_1`$ or $`k>m_1`$ and $`b_{j,k}=0`$ for $`j>n_2`$ or $`k>m_2`$. The root $`a=0`$ of $`\beta (a)`$ corresponds to the fixed point $`(0,0)`$ and the factor $`a^2`$ can be eliminated. Thus, the possible amplitudes $`a`$ are the zeros of $`\beta (a)/a^2`$, a polynomial in $`a^2`$ of degree $`Max\{n_1+m_1,n_2+m_2\}`$. There are no more than $`Max\{n_1+m_1,n_2+m_2\}`$ different solutions $`a>0`$ and therefore, the maximum number of limit cycles in this case is Max$`\{n_1+m_1,n_2+m_2\}`$ as it was proved by Iliev . This result can also be written as $`\left[\frac{Max\{degf,degg\}1}{2}\right]`$ where $`[]`$ means the integer part and $`degf`$ (resp. $`degg`$) denotes the degree of $`f(x,y)`$ (resp. of $`g(x,y)`$). If we restrict ourselves to the case of Liénard equations where $`f(x,y)=0`$ and $`g(x,y)`$ is linear in $`y`$, we recover the known result that the number of limit cycles is less or equal than $`n_2`$ in the weak nonlinear regime. The extension of this relation for the whole range of $`ϵ`$ is the still not proved Lins-Melo-Pugh conjecture on the number of limit cycles for Liénard systems .

### 3.2 Reduction of some examples to a perturbed circular period annulus

Example 1: (a particular circular period annulus)

$$\begin{array}{c}\dot{x}=yϵ(a_1x+a_2x^2+a_3x^3),\hfill \\ \dot{y}=xϵb_3x^3,\hfill \end{array}$$
(10)
Yuquan and Zhujun  have investigated the number of limit cycles of a similar and more general system. The conditions of their theorems can be translated for this particular case of weak perturbations ($`ϵ1`$) as following: (i) If $`ϵb_3>0`$, $`ϵa_1<0`$ and $`ϵa_3>0`$ then system (10) has a unique stable limit cycle (Theorem 1 in Ref. ); (ii) If $`ϵb_3>0`$, $`ϵa_1>0`$, and $`ϵa_3<0`$ then system (10) has a unique unstable limit cycle (Theorem 2 in Ref. ).

These results are confirmed by the method introduced in the previous section. Moreover, as it is explained in that section, the amplitudes $`a`$ of the limit cycles in the weakly nonlinear regime are the nontrivial positive solutions of the equation $`\beta (a)=0`$. In this case, $`f(x,y)=a_1xa_2x^2a_3x^3`$ and $`g(x,y)=b_3x^3`$. We obtain:

$$\beta (a)=\frac{\pi }{2}a^2(a_1+\frac{3}{4}a_3a^2),$$

which shows that the system has at most one limit cycle of amplitude

$$a=\sqrt{\frac{4a_1}{3a_3}}.$$

Therefore, the system has not limit cycles if $`a_1a_3>0`$ and just one limit cycle if $`a_1a_3<0`$. This periodic motion is stable if $`ϵa_1<0`$ and unstable if $`ϵa_1>0`$. Numerical simulations are in strong agreement with this analytical result. Observe that there are certain conditions of Yuquan and Zhujun theorems that can be removed in this regime in order to have just one limit cycle.

Example 2: (a system topologically equivalent to a circular period annulus)

$$\begin{array}{c}\dot{x}=y^{2l1}+ϵP(x,y),\hfill \\ \dot{y}=x^{2k1}+ϵQ(x,y),\hfill \end{array}$$
(11)
where $`k`$ and $`l`$ are positive integers and $`P(x,y)`$ and $`Q(x,y)`$ are polynomials in $`x`$ and $`y`$. The number of limit cycles of this system for $`ϵ1`$ has been studied in Ref. . The constant of motion for the unperturbed problem is $`\overline{H}(x,y)=\frac{1}{2k}x^{2k}+\frac{1}{2l}y^{2l}`$. When $`k,l>1`$ it corresponds to a degenerate center topologically equivalent to the circular period annulus obtained when $`k=l=1`$ (then $`H_c=\frac{1}{2}(x^2+y^2)`$). If we remove the time variable, we obtain

$$y^{2l1}\frac{dy}{dx}+x^{2k1}+ϵ\left[Q(x,y)P(x,y)\frac{dy}{dx}\right]=0.$$
(12)
In order to apply the method introduced in the previous section, we perform the change of variables $`\mathrm{\Gamma }:(x,y)(X,Y)`$:

$$\begin{array}{c}X=sign(x)\frac{x^k}{\sqrt{k}},\\ Y=sign(y)\frac{y^l}{\sqrt{l}}.\end{array}$$

Then, equation (12) reads

$$Y\frac{dY}{dX}+X+ϵ\left[\frac{X^{(1/k)1}}{k^{1(1/2k)}}Q(x,y)\frac{Y^{(1/l)1}}{l^{1(1/2l)}}P(x,y)\frac{dY}{dX}\right]=0,$$
(13)
with

$$\begin{array}{c}x=sign(X)(\sqrt{k}X)^{1/k}\\ y=sign(Y)(\sqrt{l}Y)^{1/l}.\end{array}$$
(14)
Equation (13) has the form (5) with $`(x,y)`$ replaced by $`(X,Y)`$ and

$$\begin{array}{c}f(X,Y)=\frac{Y^{(1/l)1}}{1l^{(1/2l)}}P(sign(X)(\sqrt{k}X)^{1/k},sign(Y)(\sqrt{l}Y)^{1/l})\\ g(X,Y)=\frac{X^{(1/k)1}}{1k^{(1/2k)}}Q(sign(X)(\sqrt{k}X)^{1/k},sign(Y)(\sqrt{l}Y)^{1/l}).\end{array}$$
(15)
Therefore, in these coordinates, if $`A`$ is the amplitude of a limit cycle $`(X,Y(X))`$ in the weakly nonlinear regime, it verifies the equation $`\overline{\beta }(A)=0`$, where

$$\overline{\beta }(A)_A^A\left\{\stackrel{~}{g}(X,Y_A(X))+\frac{X}{Y_A(X)}\stackrel{~}{f}(X,Y_A(X))\right\}𝑑X,$$

$`\stackrel{~}{f}(X,Y)`$ and $`\stackrel{~}{g}(X,Y)`$ are defined in (9) with $`f(X,Y)`$ and $`g(X,Y)`$ given in (15) and $`Y_A(X)\sqrt{A^2X^2}`$. If we undo the change of variable $`(X,Y)(x,y)`$ in the above integral and write $`a^{2k}kA^2`$, then we can write $`\overline{\beta }(A)\beta (a)`$ in terms of the amplitudes $`a`$ of the original variables $`(x,y)`$:

$$\beta (a)_a^a\left\{\stackrel{~}{Q}(x,y_a(x))\stackrel{~}{P}(x,y_a(x))\frac{dy_a(x)}{dx}\right\}𝑑x,$$

where $`y_a(x)\left[\frac{l}{k}\left(a^{2k}x^{2k}\right)\right]^{1/(2l)}`$ and $`\stackrel{~}{P}(x,y)`$ and $`\stackrel{~}{Q}(x,y)`$ are defined in (9) for $`f(x,y)=P(x,y)`$ and $`g(x,y)=Q(x,y)`$. Hence, the amplitudes, $`a`$, of the limit cycles in the weakly nonlinear regime are the nontrivial solutions of the equation $`\beta (a)=0`$.

If $`k=l`$, then $`\beta (a)`$ reads

$$\beta (a)=a_1^1\left\{\stackrel{~}{Q}(ax,a\stackrel{~}{y}(x))+\frac{x^{2k1}\stackrel{~}{y}(x)\stackrel{~}{P}(ax,a\stackrel{~}{y}(x))}{1x^{2k}}\right\}𝑑x,$$

where $`\stackrel{~}{y}(x)\left(1x^{2k}\right)^{1/(2k)}`$. If $`P(x,y)`$ and $`Q(x,y)`$ are polynomials of degree at most $`n`$, $`a^2\beta (a)`$ is a polynomial of degree $`\frac{n1}{2}`$ in $`a^2`$ and we conclude that the maximum number of limit cycles in this regime is $`\frac{n1}{2}`$, as it has been shown by Coll et al. .

Take, for instance, the particular case:

$$\begin{array}{c}\dot{x}=y^3+ϵbxy^2,\\ \dot{y}=x^3+ϵcx^2y^3.\end{array}$$
(16)
with $`k=l=2`$ and $`n=5`$. The equation $`\beta (a)=0`$ has at most one non trivial solution given by

$$a=\frac{(\pi /2)^{3/4}}{\mathrm{\Gamma }(3/4)}\sqrt{\frac{b}{c}},$$

In this case, the system has at most one limit cycle of amplitude $`a`$ emerging from the slightly perturbed period annulus. This analytical result is in agreement with the direct verification of the dynamics by integrating the equations (16).

Example 3: (a system topologically equivalent to a circular period annulus)

$$\begin{array}{c}\dot{x}=y+yx^2+ϵF(x,y),\hfill \\ \dot{y}=x+xy^2+ϵG(x,y),\hfill \end{array}$$
(17)
where $`F(x,y)=a_0x^3+a_1x^2y+a_2xy^2+a_3y^3`$ and $`G(x,y)=b_0x^3+b_1x^2y+b_2xy^2+b_3y^3`$. It has been shown that this system has, at most, two limit cycles (theorem 2.6 in Ref. ).
The constant of motion of equation (17) is $`H(x,y)=\frac{1+y^2}{1x^2}=h`$. A period annulus topologically equivalent to $`\overline{H}_c=\frac{1}{2}(x^2+y^2)=\overline{h}`$ is defined in the region $`1<x<1`$ when $`h`$ runs from $`1`$ to $`\mathrm{}`$. The origin $`(0,0)`$ is found for $`h=1`$ and the boundaries of this region are the lines $`x=\pm 1`$ which are approached when $`h\mathrm{}`$. The interior is formed by closed curves topologically equivalent to circles. The integral paths exterior to the period annulus are open curves obtained when $`h`$ increases from $`\mathrm{}`$ to zero.

Following the same steps as in the preceding example, we proceed to remove the time variable of system (17). Writing $`x^{}\frac{dx}{dy}`$:

$$(1+y^2)xx^{}+(1x^2)y+ϵ[G(x,y)x^{}F(x,y)]=0.$$
(18)
The homomorphism with the circular period annulus is established by means of the change of variables $`\mathrm{\Gamma }:(x,y)(X,Y)`$:

$$\begin{array}{cc}X=Sign(x)\sqrt{\mathrm{log}(1x^2)},& \\ Y=Sign(y)\sqrt{\mathrm{log}(1+y^2)}.& \end{array}$$

Then, the equation (18) reads

$$\begin{array}{cc}X\frac{dX}{dY}+Y+ϵ[X\frac{dX}{dY}\frac{G(Sign(X)\sqrt{1e^{X^2}},Sign(Y)\sqrt{e^{Y^2}1})}{e^{Y^2}\sqrt{1e^{X^2}}}& \\ \frac{YF(Sign(X)\sqrt{1e^{X^2}},Sign(Y)\sqrt{e^{Y^2}1})}{e^{X^2}\sqrt{e^{Y^2}1}}]=0.& \end{array}$$

Now this equation has the form (5) with the variables $`(x,y)`$ replaced by the variables $`(Y,X)`$ and

$$\begin{array}{c}f(Y,X)=X\frac{G(Sign(X)\sqrt{1e^{X^2}},Sign(Y)\sqrt{e^{Y^2}1})}{e^{Y^2}\sqrt{1e^{X^2}}},\\ g(Y,X)=Y\frac{F(Sign(X)\sqrt{1e^{X^2}},Sign(Y)\sqrt{e^{Y^2}1})}{e^{X^2}\sqrt{e^{Y^2}1}}.\end{array}$$
(19)
Therefore, we can apply the method of the previous section to this equation. Then, the amplitude $`A`$ of a limit cycle $`(Y,X(Y))`$ of this system in the weakly nonlinear regime, $`X(Y)\sqrt{A^2Y^2}`$, satisfies $`\overline{\beta }(A)=0`$, with

$$\begin{array}{c}\overline{\beta }(A)_A^A\left\{\frac{Y\stackrel{~}{F}(\sqrt{1e^{Y^2A^2}},\sqrt{e^{Y^2}1})}{e^{Y^2A^2}\sqrt{e^{Y^2}1}}+\frac{Y\stackrel{~}{G}(\sqrt{1e^{Y^2A^2}},\sqrt{e^{Y^2}1})}{e^{Y^2}\sqrt{1e^{Y^2A^2}}}\right\}𝑑Y,\end{array}$$

$`\stackrel{~}{F}(x,y)a_0x^3+a_2xy^2`$ and $`\stackrel{~}{G}(x,y)b_1x^2y+b_3y^3`$. If we undo the change of variable $`(Y,X)(x,y)`$ in the above integral and write $`A^2=\mathrm{log}(1+a^2)`$, then we can write $`\overline{\beta }(A)\beta (a)`$ in terms of the amplitudes $`a`$ of the original variables $`(x,y)`$:

$$\beta (a)_a^a\left\{\frac{y\sqrt{1+a^2}\stackrel{~}{G}(\sqrt{\frac{a^2y^2}{1+a^2}},y)}{\sqrt{a^2y^2}}+(1+a^2)\stackrel{~}{F}(\sqrt{\frac{a^2y^2}{1+a^2}},y)\right\}\frac{dy}{(1+y^2)^2}.$$

After straightforward computations we obtain that the real roots of $`\beta (a)`$ are the real solutions of the equation $`\stackrel{~}{\beta }(a)=0`$, where

$$\begin{array}{c}\stackrel{~}{\beta }(a)[b_1+a_2(1+a^2)][2+a^22\sqrt{1+a^2}]+\\ \\ b_3[2(1+a^2)^{3/2}3a^22]+a_0[2\sqrt{1+a^2}2a^2+a^4].\end{array}$$

We replace the variable $`a`$ in this equation by the variable $`B\sqrt{1+a^2}`$. We look for solutions $`a>0`$, that is, $`B>1`$. The two only possible solutions are

$$B^\pm \frac{(a_0+b_3)\pm \sqrt{(a_0+b_3)^2(a_0+a_2)(b_1+b_3)}}{a_0+a_2}.$$

Writing

$$\begin{array}{c}u\frac{a_0+b_3}{a_0+a_2},v\frac{b_1+b_3}{a_0+a_2},\end{array}$$
(20)
we have that $`B^\pm =u\pm \sqrt{u^2v}`$. Then, the number of limit cycles in the weakly nonlinear regime is $`0`$, $`1`$ or $`2`$ depending on $`a_0,a_2,b_1,b_3`$ (see figure 2):

(i) If $`u>1`$ and $`v+2u+1>0`$ or $`v>u^2`$, neither $`B^+`$ neither $`B^{}`$ are greater than $`1`$. Therefore, the system has not limit cycles.
(ii) If $`v+2u+1<0`$, then $`B^+>1`$, $`B^{}<1`$. Therefore, the system has a unique limit cycle with amplitude $`a=\sqrt{(B^+)^21}`$.
(iii) If $`v+2u+1>0`$ and $`v<u^2`$, then $`B^+,B^{}>1`$. Therefore, the system has two limit cycles with amplitudes $`a^+=\sqrt{(B^+)^21}`$ and $`a^{}=\sqrt{(B^{})^21}`$.

Computer integration of equation (17) is in agreement with these analytical results. Theorem 2.6 of Ref.  states that the maximum number of limit cycles of this system is two. By applying the method proposed in Section 2 we obtain a more detailed information: the exact number and amplitudes of the limit cycles of this system as a function of the parameters $`\{a_k,b_kk=1,2,3\}`$.

## 4 Conclusions

Nowadays it is well known that the inclusion of a slight dissipation in a conservative system destroys almost completely the geometrical aspect of the phase space of the unperturbed motion. But the analytical methods available for the study of such perturbed systems suffer a lack of predictive power when the nonlinearities start to dominate the dynamics.

In this work, we have carried out the calculation of the precise number and amplitude of the periodic orbits surviving to a weak perturbation of different topologically equivalent period annulus. First, the equation $`\beta (a)=0`$ containing that information about the limit cycles emerging from a general slight perturbation of a continuum set of circles is stated. Essentially this equation corresponds to the first Melnikov function associated to the system, although it has been obtained by an alternative line of reasoning. Second, the dynamical equations of other topologically equivalent perturbed centers have been reduced to equivalent equations of the former case. Then, the calculation of the number and amplitude of their limit cycles has also been carried out in the same way. The numerical integration of these different systems supports the analytical calculations and confirms the theoretical predictions.

We must stress that, on the one hand, our method has given the correct maximum number of limit cycles in all the examples analyzed. On the other hand, it has allowed us to predict the exact number and amplitude of the limit cycles emerging from different topologically equivalent perturbed period annulus. It seems that in order to find the number and amplitude of limit cycles, it is not important to maintain the condition of a planar system to be a polynomial system. Perhaps it would be more interesting to try to identify the class of all its geometrical and topologically equivalent systems and to perform the calculations in the simplest system of this class.

Acknowledgements: J.L. López acknowledges DGCYT (BFM2000-0803) and Gobierno de Navarra (Res. 92/2002) for financial support.

Figure Captions

Figure 1: A typical phase portrait of equation (5). The limit cycles of amplitudes $`(a_1,a_2)`$, $`(b_1,b_2)`$, $`\mathrm{}`$, enclose the origin. Stable and unstable limit cycles alternate. The order of this alternation depends on the stability of the origin.

Figure 2: The complete bifurcation diagram of system (17). The system has no periodic solutions in the white region, one limit cycle in the grey region and two limit cycles in the black region. See example 3 for details of the calculation of $`(u,v)`$ from the original parameters of equation (17).