Perturbed rank 2 Poisson systems and periodic orbits on Casimir invariant manifolds

A class of n-dimensional Poisson systems reducible to an unperturbed harmonic oscillator shall be considered. In such case, perturbations leaving invariant a given symplectic leaf shall be investigated. Our purpose will be to analyze the bifurcation phenomena of periodic orbits as a result of these perturbations in the period annulus associated to the unperturbed harmonic oscillator. This is accomplished via the averaging theory up to an arbitrary order in the perturbation parameter e. In that theory we shall also use both branching theory and singularity theory of smooth maps to analyze the bifurcation phenomena at points where the implicit function theorem is not applicable. When the perturbation is given by a polynomial family, the associated Melnikov functions are polynomial and tools of computational algebra based on Gr\"obner basis are employed in order to reduce the generators of some polynomial ideals needed to analyze the bifurcation problem. When the most general perturbation of the harmonic oscillator by a quadratic perturbation field is considered, the complete bifurcation diagram (except at a high codimension subset) in the parameter space is obtained. Examples are given.


Introduction
Finite-dimensional Poisson systems (see [15,16,19] and references therein for an overview) have a significant presence in most domains of physics and applied mathematics. The specific format of Poisson systems has allowed the development of many tools for their analysis (for instance, see [5]- [7], [16] and references therein for a sample). In addition, the relevance of Poisson systems arises from the fact that they constitute a generalization of classical Hamiltonian systems comprising nonconstant structure matrices as well as odd-dimensional vector fields. Additionally, Poisson system format is invariant under general diffeomorphic transformations, therefore not being restricted to the use of canonical transformations.
Consider a smooth vector field having a finite-dimensional Poisson structure dx dt = J (x) · ∇H(x) (1.1) of dimension n and rank r = 2 ≤ n constant in an open set Ω ⊆ R n . In (1.1) J (x) and H(x) are the structure matrix and Hamiltonian function, respectively. Then under these hypotheses for each point x 0 ∈ Ω there is (at least locally in a neighborhood Ω 0 ⊂ Ω of x 0 ) a complete set of functionally independent Casimir invariants {D 3 (x), . . . , D n (x)} in Ω 0 , as well as a transformation x → Φ D (x) = y where Φ D is a smooth diffeomorphism in Ω 0 bringing the system (1.1) into its Darboux canonical form. Thus, beyond the fact that Poisson systems are a formal generalization of classical Hamiltonian flows, Darboux Theorem provides the dynamical basis for such a generalization. The explicit construction of the Darboux coordinates may be a complicated task in general, specially in the case of their global determination, in which the transfer of results between the Poisson and the classical Hamiltonian formats is optimal for the applications. Similarly the explicit construction in closed form of the Casimir invariants D j is not an easy task, recall that these invariants are the solutions of the system of PDEs J (x) · ∇D j (x) = 0. In this article, a class of n-dimensional Poisson systems reducible to an unperturbed harmonic oscillator shall be considered. In such case, perturbations leaving invariant a given symplectic leaf shall be investigated. Our purpose will be to analyze the bifurcation phenomena of periodic orbits as a result of these perturbations in the period annulus associated to the unperturbed harmonic oscillator. This is accomplished via the averaging theory up to an arbitrary order in the perturbation parameter ε. In that theory we shall also use both branching theory and singularity theory of smooth maps to analyze the bifurcation phenomena at points where the implicit function theorem is not applicable. When the perturbation is given by a polynomial family, the associated Melnikov functions are polynomials. This sentence is a consequence of the fact that the Lagrange standard form (2.7) associated to the above perturbation problem can be written as dr/dθ = ∑ i≥1 G i (θ , r) ε i with G i (θ , r) polynomials in r for all i ∈ N, see Lemma 9 in [11] for a proof. There fore tools of computational algebra based on Gröbner basis are employed in order to reduce the generators of some polynomial ideals needed to analyze the bifurcation problem. When the most general perturbation of the harmonic oscillator by a quadratic perturbation field is considered, the complete bifurcation diagram (except at a high codimension subset) in the parameter space is obtained.

Reduction
Let us then assume, without loss of generality, a Hamiltonian quasi-harmonic for x 1 and x 2 . 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 with H * (y) = H • Φ −1 (y) =Ĥ 1 2 (y 2 1 + y 2 2 ), y 3 , . . . , y n , and J * (y) = η(y) · J D , where is the Darboux canonical form matrix for the rank-2 case, where O n−2 denotes the null square matrix of order n − 2. Finally, rescaling the time as t → τ with dτ = η dt we obtain the Darboux canonical form in Ω of the Poisson system (1.1). 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 isH 1 2 (y 2 1 + y 2 2 ) ≡Ĥ 1 2 (y 2 1 + y 2 2 ), c 3 , . . . , c n . Now denote asH (z) = dH(z)/dz. Then the reduction is completed by means of an additional time reparametrization τ → ρ with dρ = µ(y 1 , y 2 ) dτ, where µ(y 1 , y 2 ) =H 1 2 (y 2 1 + y 2 2 ) . The outcome is a one degree of freedom harmonic oscillator of Hamiltonian H (y 1 , y 2 ) = 1 2 (y 2 1 + y 2 2 ) and time ρ.

Perturbations leaving invariant a given simplectic leaf
We consider now the analytical perturbations of the initial Poisson system (1.1) where ε = 0 is a small perturbation real parameter and F is an analytic vector field in Ω depending analytically on the parameter ε and satisfying F(0; ε) = 0 and ∇ x F(0; ε) = 0. We will restrict our selves to those pertur= bation fields F(x; ε) that leave invariant for the flow of (2.3) a given simplectic leaf L c = ∩ n j=3 {D j (x) = c j } of the Poisson system (1.1) for certain c = (c 3 , . . . , c n ) ∈ R n−2 such that L c ∩ Ω = / 0. Performing the Darboux canonical form reduction of the previous subsection, we obtain that (2.3) becomes the analytic system dy dτ = J D · ∇H * (y) + εF * (y; ε) (2.4) defined in Ω * = Φ(Ω). Since L c becomes an invariant surface of the perturbed system (2.3), under these conditions, diffeomorphism Φ defined in (2.2) and the rescaling of time t → τ previously characterized transform (2.3) in Ω into a system in Ω * which can be restricted to Φ(L c ) leading to a two dimensional system because dim(L c ) = 2. More specifically (2.4) can be written as Finally, the restriction to Φ(L c ) combined with the time rescaling τ → ρ described in the previous subsection leads to where H (y 1 , y 2 ) = 1 2 (y 2 1 + y 2 2 ). The reduction to a perturbed harmonic oscillator is thus accomplished.

The Lagrange standard form of averaging theory
In polar coordinates, y 1 = r cos θ , y 2 = r sin θ , system (2.5) becomeṡ 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 by taking θ as the new independent variable. Recall that any 2π-periodic solution of (2.7) corresponds biunivocally with a periodic orbit of (2.3) on an arbitrarily large compact set included in L c ∩ Ω. Therefore, system (2.7) is 2π-periodic in variable θ and is in the Lagrange standard form of averaging theory.

Example: Maxwell-Bloch equations
The real-valued Maxwell-Bloch system (see [4] and references therein) is given by the following polynomial vector field in R 3 :ẋ Since rank(J ) = 2 everywhere it has one independent Casimir invariant which can be chosen as

3) if and only if there is a real analytic function
Note that the Maxwell-Bloch Hamiltonian is quasi-harmonic in terms of variables x 2 and x 3 , . Accordingly, we can perform the change of variables given by the diffeomorphism . This is the natural choice in order to arrive to a harmonic oscillator. Observe that under such transformation, the surface L c is transformed into the half-plane Π = {y ∈ Ω * ⊂ R 3 : y 3 = c} defined in Ω * = Φ(Ω) = {y ∈ R 3 : y 1 = 0, y 2 = 0, y 3 > y 2 }. The perturbed system (2.3) defined in Ω * adopts the forṁ and R(y) is an analytic function in Ω * . Now we restrict system (2.9) to its invariant plane Π and rescale the time t → τ with dτ = η dt to obtain the planar system which is defined on Π. Notice that in the particular case in which the perturbation (A(x), B(x),C(x)) is polynomial with A and B even and odd, respectively, in the variable x 1 , that is having the form is also a polynomial perturbation of the harmonic oscillator.

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 R 3 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 L c 2 ≡ {x ∈ R 3 : D(x) = c 2 }. The Hamiltonian is quasi-harmonic for every pair of variables. For instance, in terms of x 1 and According to the reduction procedure assumptions, we have ϕ i (x 1 , x 2 ) = 0 provided x 1 = 0 and x 2 = 0, and in addition we assume without loss of generality µ 3 > µ 1 and µ 3 > µ 2 . Let us also define the semispheres Consider now the most general analytic perturbation in R 3 \{x 3 = 0} of the Euler top, leaving invariant the semispheres L + c 2 and L − c 2 : 12) with P, Q and R analytic functions everywhere in R 3 . We then perform the following diffeomorphic change of variables: The perturbed system (2.12) restri= cted to L + c 2 adopts the forṁ y 1 = −κ 13 κ 23 y 3 − (y 1 /κ 13 ) 2 − (y 2 /κ 23 ) 2 ∂ H ∂ y 2 + εP(y 1 , y 2 , y 3 ) , 14) with H(y 1 , y 2 , y 3 ) = 1 2 (y 2 1 + y 2 2 ) + 1 2µ 3 y 3 . The perturbed system (2.12) restricted to the semispace x 3 < 0 is given by (2.14) changing the sign in the right-hand side ofẏ 1 andẏ 2 . Then, the restriction of system (2.14) to L + c 2 is given by the analytic systeṁ where H (y 1 , y 2 ) = 1 2 (y 2 1 + y 2 2 ). Finally, we introduce a time reparametrization t → τ of the form dτ = η dt, with η = −κ 13 κ 23 c 2 − (y 1 /κ 13 ) 2 − (y 2 /κ 23 ) 2 which completes the reduction to the form (2.5) of a perturbed harmonic oscillator.
The cyclicity of P under perturbations (P, Q) with |ε| 1 is the maximum number of limit cycles bifurcating from the circles that foliates P. A detailed analysis of the homogeneous case for which P and Q are homogeneous polynomials of the same arbitrary degree is given in [9] where its is shown that the cyclicity of P is zero.
Later in [13] the cyclicity of P under arbitrary polynomial perturbation fields (P(y 1 , y 2 ; ε), Q(y 1 , y 2 ; ε)) is analyzed but now allowing the coefficients to depend analytically on ε, that is P, Q ∈ R{ε}[y 1 , y 2 ]. In [13] it is derived the global upper bound [ (n − 1)/2] on the cyclicity of P where n = max{deg(P), deg(Q)} and is the order of the associated first Melnikov function which is not identically zero. Also in [13] some cases where the above upper bound is sharp are shown.
An interesting question arises if we assume that P, Q ∈ R m [ε][y 1 , y 2 ], that is, the coefficients of (P, Q) are polynomial functions of ε having some fixed maximum degree m: to find the bifurcation diagram of limit cycles in P in the parameter space. We consider here the simplest case with respect to the degrees, namely, (m, n) = (1, 2). Thus we consider the most general perturbation of a harmonic oscillator like (2.5) by a quadratic perturbation field (P, Q) whose coefficients are linear functions of the perturbation parameter ε. Moreover, the right hand side can be taken without loss of generality (after a rotation in the phase plane) in the called Bautin form (see [2]) with linear coefficients A i (ε) = a i0 + a i1 ε for i = 2, 3, 4, 5, 6. The resulting perturbation coefficients a i j are collected into the vector parameter λ ∈ R 10 .
Remark 3.1. After [2], it is well known that the origin is a center of family (3.1) for any ε ∈ R if and only if one of the following four conditions is fulfilled: Introducing polar coordinates y 1 = r cos θ , y 2 = r sin θ , and for |ε| sufficiently small, any systeṁ y 1 = y 2 + εP(y 1 , y 2 ; ε),ẏ 2 = −y 1 + εQ(y 1 , y 2 ; ε) and in particular system (3.1) is transformed into the analytic differential equation which is defined on the cylinder {(r, θ ) ∈ (R + ∪ {0}) × S 1 } with S 1 = R/2πZ and satisfies F (θ , r; λ , 0) ≡ 0. Therefore, equation (3.2) is written in the standard Lagrange form of the averaging theory with period 2π. The classical tool of averaging allows us to analyze the 2π-periodic solutions of (3.2), see for example the book [17] or, for recent advances, the papers [8] and [14]. The solution r(θ ; z, λ , ε) of (3.2) with initial condition r(0; z, λ , ε) = z ∈ R + admits the convergent power series expansion near ε = 0 like r(θ ; z, λ , ε) = z + ∑ j≥1 r j (θ , 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 : R + × R 12 × I → R + 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 0 ∈ R + 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 y 2 1 + y 2 2 = z 2 0 included in the period annulus P of the unperturbed harmonic oscillator.
In summary, the displacement map d is expressed as the following convergent series expansion and the coefficient functions f i (z; λ ) = r i (2π, z, λ ) can be computed by a recursive procedure, see [11] for the general structure. We call f i the i-th averaged function (also called i-th Melnikov function in the literature).
We say that a branch of limit cycles bifurcates from the circle y 2 1 + y 2 2 = z 2 0 with z 0 ∈ R + 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 0 and d(z * (λ , ε), λ , ε) ≡ 0. It is well known (see [17], for example) that in such a case z 0 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
, we can consider the polynomial ideal I generated by its coefficients ξ i j ∈ R[λ ] in the ring R[λ ]. We also can consider the ascending chain of ideals I 2 ⊆ I 3 ⊆ · · · ⊆ I k = I where I s = ξ i j : 2 ≤ i ≤ s . Since I is a Noetherian ring, the above chain stabilizes at, say, the moment k ∈ N. The former implies that if the parameters λ = λ * ∈ I k , then d(z, λ * , ε) ≡ 0 and the origin becomes a center of (3.1). Remark 3.3. We summarize here the classical averaging theory applied to the differential equation (3.2). Assume that z 0 ∈ R + 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 0 for |ε| 1. Then the following statements hold: 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 P 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 λ ∈ R 10 (here λ denote the vector parameter whose components are a i0 and a i1 for i = 2, 3, 4, 5, 6.): a 50 (a 30 − a 60 ), a 31 a 50 + a 30 a 51 − a 51 a 60 − a 50  Let N(λ ) be the number of limit cycles that bifurcate from its period annulus P = R 2 \{(0, 0)} as the perturbation parameter ε slightly varies from zero. Then the following holds: Proof. Straightforward computations produce the following averaged functions for system (3.1): where ξ i j ∈ R[λ ] are the polynomials in the parameters of family (3.1). In what follows we shall denote byξ i j the remainder of ξ i j upon division by a Gröbner basis of the ideal generated by all the ξ ks with k < i in the polynomial ring R[λ ]. This remainder can be computed, for instance, with the functions PolynomialReduce and GroebnerBasis of the computer algebra system Mathematica c . Another option is the use of reduce with the software Singular c . The nonidentically zero polynomialsξ i j ∈ R[λ ] are listed in the statement of the theorem. After such reduction, we will consider the polynomials: =46rom the expression of f 2 (z; λ ) andf 3 (z; λ ) we deduce statements (i) and (ii) respectively while from the expression off 4 (z; λ ) we obtain (iii) and (iv). Next (v) and (vi) are obtained from the expressions off 5 (z; λ ) andf 6 (z; λ ).
Notice that the complete limit cycle bifurcation diagram of P in the parameter space R 10 for family (3.1) when ξ 20 =ξ 30 =ξ 40 =ξ 42 =ξ 51 =ξ 61 =ξ 63 = 0 (that is for parameters λ = λ * lying in the real variety associated with I 6 ) is not presented. Unfortunately the massive computations to obtain f 7 (z; λ ), hencef 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 I 6 = I because there are parameters in I 6 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 ξ 20 , ξ 40 and ξ 42 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 0 of f (·; λ ) several methods can be employed. Among them branching theory and singularity theory applied to the reduced displacement map δ (z, λ , ε) = f (z; λ ) + ∑ i≥1 f +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 , 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 R 2 mapping the origin to (z 0 , 0) and preserving orientation whereas U(z, ε) > 0. A different approach dealing with the degenerate case for which z 0 is a multiple zero of f (·; λ ) and f k (z 0 ; λ ) = 0 for any k ∈ N can be found in [8].
In the next example, the analysis of multiple zeroes of f (·; λ ) is needed.