2.1. Darboux canonical form and harmonic oscillator form

For the Poisson system (1.1) and under the assumptions previously stated, we shall say that the Hamiltonian function H(x) is quasi-harmonic for variables x_i and x_j in a domain Ω ⊂ ℝⁿ if, by definition, it can be written in the form H(x) ≡ H(φ(x_i, x_j), D₃(x),...,D_n(x)) where φ(x_i, x_j) admits, in the region of interest Ω, at least one decomposition of the kind φ(x_i, x_j) = φ₁(x_i, x_j) + φ₂(x_i, x_j) with φ_k(x_i, x_j) > 0 for k = 1, 2, and with application (x_i, x_j) ↦ (φ₁(x_i, x_j), φ₂(x_i, x_j)) being invertible. Let us then assume, without loss of generality, a Hamiltonian quasi-harmonic for x₁ and x₂. We finally assume that the following non-degeneracy Jacobian condition holds

Then the following change of variables is to be performed:

By definition, Ω is the open set such that the Poisson system (1.1) is defined and has rank r = 2, and in addition Φ|_Ω is a diffeomorphism by (2.1). Then the transformed system can be written as

Moreover, we can proceed further and reduce the system completely to a classical harmonic oscillator. For this, we first rectrict ourselves to one symplectic leaf y_i = y_i(0) = c_i, for i = 3,...,n. We are thus left with a planar classical Hamiltonian system for which the structure matrix is the 2 × 2 symplectic matrix, and the Hamiltonian is H˜(12(y12+y22))≡H^(12(y12+y22),c3,…,cn). Now denote as H˜′(z)=dH˜(z)/dz. Then the reduction is completed by means of an additional time reparametrization τ ↦ ρ with dρ = μ(y₁, y₂) dτ, where μ(y1,y2)=H˜′(12(y12+y22)). The outcome is a one degree of freedom harmonic oscillator of Hamiltonian ℋ(y1,y2)=12(y12+y22) and time ρ.

2.2. Perturbations leaving invariant a given simplectic leaf

We consider now the analytical perturbations of the initial Poisson system (1.1)

2.3. The Lagrange standard form of averaging theory

In polar coordinates, y₁ = r cosθ, y₂ = r sinθ, system (2.5) becomes

Notice that this system is only well defined for r > 0. Moreover, in this region, since for sufficiently small ε we have θ˙<0 in an arbitrarily large ball centered at the origin, we can rewrite the differential system (2.6) in such ball into the form

2.4. Example: Maxwell-Bloch equations

The real-valued Maxwell-Bloch system (see [4] and references therein) is given by the following polynomial vector field in ℝ³:

Equations (2.8) can be written as a Poisson system (1.1) with Hamiltonian H(x)=12(x22+x32) and structure matrix

Since rank(𝒥) = 2 everywhere it has one independent Casimir invariant which can be chosen as D(x)=x3+12x12.

Let F(x; ε) = (A(x), B(x), C(x)) be the perturbation vector field in (2.3). We make the following statement: the perturbed field (2.3) has the invariant surface ℒ_c = {x ∈ ℝ³ : D(x) = c} for some arbitrary real constant c ∈ ℝ if and only if D(x) − c divides the analytic function xA(x) +C(x). The proof is as follows: S_c is an invariant surface of (2.3) if and only if there is a real analytic function K in ℝ³ such that 𝒴 (D(x)−c) = K(x)(D(x)−c) where 𝒴 = A(x)∂_x1 +B(x)∂_x2 +C(x)∂_x3 is the vector field (linear differential operator) associated to F. Then, direct computations give xA(x) +C(x) = K(x)(D(x) − c) thus proving the claim.

Note that the Maxwell-Bloch Hamiltonian is quasi-harmonic in terms of variables x₂ and x₃, namely H(x)φ₁(x₂, x₃)+φ₂(x₂, x₃), with ϕ1(x2,x3)=12x22 and ϕ2(x2,x3)=12x32. Accordingly, we can perform the change of variables given by the diffeomorphism x = (x₁, x₂, x₃) ↦ y = (y₁, y₂, y₃) = Φ(x) = (x₂, x₃, D(x)) defined in the region Ω = {x ∈ ℝ³ : x_i ≠ 0, i = 1, 2, 3}. This is the natural choice in order to arrive to a harmonic oscillator. Observe that under such transformation, the surface ℒ_c is transformed into the half-plane Π = {y ∈ Ω^* ⊂ ℝ³ : y₃ = c} defined in Ω^* = Φ(Ω) = {y ∈ ℝ³ : y₁ ≠ 0, y₂ ≠ 0, y₃ > y₂}. The perturbed system (2.3) defined in Ω^* adopts the form

2.5. Example: Euler top

As a second instance of the reduction procedure consider the Euler equations, which describe the rotation of a rigid body:

In system (2.11) each variable x_i denotes the ith component of angular momentum, and constants μ_i are the moments of inertia about the coordinate axes, both for i = 1, 2, 3. Energy is conserved for this vector field, and actually this is a Poisson system [1, 16] in terms of the following structure matrix:

Obviously the rank of the structure matrix is 2 everywhere in ℝ³ except at the origin. The Hamiltonian, which is the total (kinetic) energy, can be written as:

Euler top has received a significant attention in the Poisson system framework, for instance see [5] and references therein. From the point of view of the study of periodic solution bifurcations after perturbations of the Euler top, see [3]. Excluding the origin, there is one independent Casimir invariant, which can be taken as the square of the angular momentum norm:

Accordingly, we shall denote the symplectic leaves as ℒ_c² ≡ {x ∈ ℝ³ : D(x) = c²}. The Hamiltonian is quasi-harmonic for every pair of variables. For instance, in terms of x₁ and x₂ we have H(x)=ϕ1(x1,x2)+ϕ2(x1,x2)+12μ3D(x), where ϕi(x1,x2)=12κi32xi2, for i = 1, 2, and

According to the reduction procedure assumptions, we have φ_i(x₁, x₂) ≠ 0 provided x₁ ≠ 0 and x₂ ≠ 0, and in addition we assume without loss of generality μ₃ > μ₁ and μ₃ > μ₂. Let us also define the semispheres

Consider now the most general analytic perturbation in ℝ³\{x₃ = 0} of the Euler top, leaving invariant the semispheres ℒc2+ and ℒc2−:

Remark 3.1.

After [2], it is well known that the origin is a center of family (3.1) for any ε ∈ ℝ if and only if one of the following four conditions is fulfilled:

1. (a)

   A₄(ε) = A₅(ε) ≡ 0;

2. (b)

   A₃(ε) − A₆(ε) ≡ 0;

3. (c)

   A5(ε)=A4(ε)+5(A3(ε)−A6(ε))=A3(ε)A6(ε)−2A62(ε)−A22(ε)≡0;

4. (d)

   A₂(ε) = A₅(ε) ≡ 0.

Introducing polar coordinates y₁ = r cosθ, y₂ = r sinθ, and for |ε| sufficiently small, any system ẏ₁ = y₂ + εP(y₁, y₂; ε), ẏ₂ = −y₁ + εQ(y₁, y₂; ε) and in particular system (3.1) is transformed into the analytic differential equation

The solution r(θ; z, λ, ε) of (3.2) with initial condition r(0; z, λ, ε) = z ∈ ℝ⁺ admits the convergent power series expansion near ε = 0 like r(θ;z,λ,ε)=z+∑j≥1rj(θ,z,λ)εj where the coefficient functions r_j are real analytic. The function r(·; z, λ, ε) is defined on the interval [0, 2π] provided that ε is close enough to 0, hence we can define the displacement map d : ℝ⁺ × ℝ¹² × I → ℝ⁺ with I some real interval containing the origin as d(z, λ, ε) = r(2π; z, λ, ε) − z. From this definition we see that the isolated positive zeros z₀ ∈ ℝ⁺ of d(·, λ, ε) are just the initial conditions for the 2π-periodic solutions of (3.2), which clearly are in one-to-one correspondence with the limit cycles of system (3.1) bifurcating from the circle y12+y22=z02 included in the period annulus 𝒫 of the unperturbed harmonic oscillator.

In summary, the displacement map d is expressed as the following convergent series expansion

We say that a branch of limit cycles bifurcates from the circle y12+y22=z02 with z₀ ∈ ℝ⁺ if there is a function z^*(λ, ε) (which may be defined only for values of ε on a half-neighborhood of zero) such that z^*(λ, 0) = z₀ and d(z^*(λ, ε), λ, ε) ≡ 0. It is well known (see [17], for example) that in such a case z₀ must be a zero of the function f_ℓ(·; λ) where ℓ is the first subindex such that f_ℓ(z; λ) ≢ 0, that is the first non-identically zero averaged function is the ℓ-th.

Remark 3.2.

Since the averaged functions fi(z;λ)=zmj∑j=0njξij(λ)zj∈ℝ[λ][z], we can consider the polynomial ideal ℐ generated by its coefficients ξ_ij ∈ ℝ[λ] in the ring ℝ[λ]. We also can consider the ascending chain of ideals

Remark 3.3.

We summarize here the classical averaging theory applied to the differential equation (3.2). Assume that z₀ ∈ ℝ⁺ is a zero of f_ℓ(·; λ^*), the first non identically zero averaged function and let N be the number of isolated branches of 2π-periodic solutions of (3.2) with parameters λ = λ^* bifurcating from z₀ for |ε| ≪ 1. Then the following statements hold:

1. (i)

   If z₀ is simple then N = 1.

2. (ii)

   If z₀ is multiple of multiplicity k¯, then N≤k¯.

Notice that (i) is a simple consequence of the Implicit Function Theorem while for (ii) it is required the Weierstrass Preparation Theorem.

The following result is new in the literature and is a first approach to obtain the complete limit cycle bifurcation diagram of 𝒫 in the parameter space of the quadratic family with linear polynomials of ε in their coefficients. The analysis includes the orders ℓ ≤ 6 of the associated first Melnikov functions.

Theorem 3.1.

Let us consider the perturbed harmonic oscillator given by family (3.1) and the following set of polynomials in their parameters λ ∈ ℝ¹⁰ (here λ denote the vector parameter whose components are a_i0 and a_i1 for i = 2, 3, 4, 5, 6.):

Let N(λ) be the number of limit cycles that bifurcate from its period annulus P = ℝ²\{(0, 0)} as the perturbation parameter ε slightly varies from zero. Then the following holds:

1. (i)

   If ξ₂₀ ≠ 0 then N = 0;

2. (ii)

   If ξ₂₀ = 0 and ξ^30≠0 then N = 0;

3. (iii)

   If ξ20=ξ^30=ξ^42 then N = 0;

4. (iv)

   If ξ20=ξ^30 but ξ^42≠0 then, defining s1=ξ^40/ξ^42, we have that N = 1 or N = 0 according to wether s₁ < 0 or s₁ ≥ 0, respectively;

5. (v)

   If ξ20=ξ^30=ξ^40=ξ^42=0 and ξ^51≠0 then N = 0;

6. (vi)

   If ξ20=ξ^30=ξ^40=ξ^42=ξ^51=0 but ξ^63≠ then, defining s2=ξ^61/ξ^63, we have that N = 1 or N = 0 according to whether s₂ < 0 or s₂ ≥ 0, respectively.

Proof.

Straightforward computations produce the following averaged functions for system (3.1):

Notice that the complete limit cycle bifurcation diagram of 𝒫 in the parameter space ℝ¹⁰ for family (3.1) when ξ20=ξ^30=ξ^40=ξ^42=ξ^51=ξ^61=ξ^63=0 (that is for parameters λ = λ^* lying in the real variety associated with ℐ₆) is not presented. Unfortunately the massive computations to obtain f₇(z; λ), hence f^7(z;λ), in the proof of Theorem 3.1 do not seem to be possible in our computer. In other words, for family (3.1) we are unable to get the ideal stabilization explained in Remark 3.2. The reason is that we can check that ℐ₆ ≠ ℐ because there are parameters in ℐ₆ for which the origin is not a center of (3.1) as can be easily seen by using Remark 3.1. Anyway, the bifurcation diagram can be made complete with a further case-by-case explicit analysis of the 10 subcases that arise after the vanishing of the factors in the expressions of ξ₂₀, ξ₄₀ and ξ₄₂ which are the simpler ones.

We remark on the other hand that in all the cases exposed in Theorem 3.1 we have obtained simple zeroes of the corresponding averaged function f_ℓ(·; λ). In order to compute the actual value (not only its upper bound as in part (ii) of Remark 3.3) of the number of branches bifurcating from a multiple zero z₀ of f_ℓ(·;λ) several methods can be employed. Among them branching theory and singularity theory applied to the reduced displacement map δ(z,λ,ε)=fℓ(z;λ)+∑i≥1fℓ+i(z;λ)εi are worth mentioning. Branching theory uses the Newton’s diagram of δ (see [18]) to analyze the local structure of the zeroes of δ near (z, ε) = (z₀, 0). The approach of singularity theory of smooth functions (see for example [12]) is completely different: the goal is to find when λ = λ^* a normal form δ^(z,ε) of δ(z, ε) such that U(z,ε)δ(Z(z,ε),Λ(ε))=δ^(z,ε) where (z, ε) ↦ (Z(z, ε), Λ(ε)) is a local diffeomorphism of ℝ² mapping the origin to (z₀, 0) and preserving orientation whereas U(z, ε) > 0. A different approach dealing with the degenerate case for which z₀ is a multiple zero of f_ℓ(·;λ) and f_k(z₀; λ) = 0 for any k ∈ ℕ can be found in [8].

In the next example, the analysis of multiple zeroes of f_ℓ(·;λ) is needed.

Proposition 3.1.

Let us consider the perturbed harmonic oscillator given by system ẏ₁ = y₂ + εP₃(y₁, y₂; ε), ẏ₂ = −y₁ + εQ₃(y₁, y₂; ε) with the cubic perturbation

Proof.

Straightforward computations produce the following averaged functions for system (3.1):

Therefore, the reduced displacement map δ(z, ε) = d(z, ε)/ε³ has the form δ(z, ε) = f₃(z) + f₄(z)ε + 𝒪(ε²) where z0=2/2∈ℝ+ is a multiple zero of f₃ of multiplicity k¯=2. We know then that at most 2 limit cycles can bifurcate from the circle x2+y2=z02. The following analysis will show that actually this bound is sharp. Indeed, since f₄(z₀) ≠ 0, using singularity theory of smooth maps (see [12]), we deduce that δ is strongly equivalent to the normal form δ˜(z,ε)=δ1z2+δ2ε where δ_j are ±1 according to the signs

We recall here that δ˜(z,ε) and δ(z, ε) are strongly equivalent if they are related by U(z,ε)δ(Z(z,ε),ε)=δ˜(z,ε) where z ↦ Z(z, ε) is a local diffeomorphism of ℝ mapping the origin to z₀ and preserving orientation, and U(z, ε) is a positive function. Notice that if N_δ (ε) denotes the number of local zeros of δ(·, ε) near z₀ and Nδ˜(ε) the number of local zeros of δ˜(⋅,ε) near 0 then we arrive at the important consequence for our purpose that Nδ=Nδ˜(ε).

In our case δ₁ = 1 and δ₂ = −1 so that δ˜(⋅,ε) has exactly two zeros z±*(ε)=±ε which only appear when ε > 0 so that z±*∈ℝ. Therefore, going back we conclude that exactly two limit cycles bifurcate from the circle x2+y2=z02 when ε > 0 and no limit cycle bifurcation occurs with the contrary sign of ε.