A&A 489, 1055-1063 (2008)
DOI: 10.1051/0004-6361:200810023

Periodic orbits in the logarithmic potential

G. Pucacco^1, - D. Boccaletti² - C. Belmonte³

1 - Physics Department, University of Rome ``Tor Vergata", Via della Ricerca Scientifica 1, 00133 Rome, Italy
2 - Mathematics Department, University of Rome ``La Sapienza", P.le A. Moro 2, 00185 Rome, Italy
3 - Physics Department, University of Rome ``La Sapienza", P.le A. Moro 2, 00185 Rome, Italy

Received 22 April 2008 / Accepted 16 July 2008

Abstract
Context. We investigate periodic orbits in galactic potentials by developing analytical methods.
Aims. We evaluate the quality of the approximation of periodic orbits in the logarithmic potential constructed using perturbation theory based on Hamiltonian normal forms.
Methods. The solutions of the equations of motion corresponding to periodic orbits are obtained as series expansions computed by inverting the normalizing canonical transformation. To improve the convergence of the series, a resummation based on a continued fraction may be performed. This method is analogous to the Prendergast method, which searches for approximate rational solutions.
Results. It is shown that with a normal form truncated at the lowest order incorporating the relevant resonance it is possible to construct accurate solutions both for normal modes and periodic orbits in general position.

Key words: galaxies: kinematics and dynamics - methods: analytical

1 Introduction

In his book on Dynamical Astronomy, Contopoulos (2004) encouraged investigation of the higher-order versions of the Prendergast (1982) method to solve non-linear differential equations. The original method was applied by Contopoulos & Seimenis (1990, hereafter CS90) to periodic orbits in the logarithmic potential and consists of approximating the exact solution with rational trigonometric functions. Even though the trigonometric series used in the rational approximation are truncated at the first non-trivial order, in CS90 it was shown that the quality of the fit to the exact result is quite good over a wide range of energy and ellipticity. On this basis, it is natural to presume that higher-order truncations would improve the quality of the prediction.

However, even the simplest version of the Prendergast (1982) method has two problematic aspects: 1) the choice of the dominant harmonic in the trigonometric series has to be made on the basis of knowledge about the orbit type under study; 2) the determination of the coefficients in the series, which depend on the parameters of the system and on initial conditions, originates from a non-linear algebraic system the solution of which must in general be performed numerically. This second aspect diminishes much of its simplicity, particularly if we attempt higher order truncations and consider the growth of the number of unknown coefficients.

In this paper, we would like to explore the link between the Prendergast-Contopoulos approach and the approximation of orbital solutions found with a resonant normal form. The motivation for this study stems from the idea of rooting a simplified version of the rational solution method into the frame of a modified normalization algorithm to devise a completely analytical approach. This step was proposed (Pucacco et al. 2008) to exploit a resummation technique based on continued fractions to speed up the convergence of series obtained in the framework of normal form perturbation theory. This technique extends the quality of predictions concerning the instability of normal modes and consequent bifurcations of families of boxlets (Belmonte et al. 2007).

In analogy with CS90, we apply this approach to investigate periodic orbits in the logarithmic potential (Binney & Tremaine 1987). We find analytical solutions to the equations of motion for the normal modes and the main low-order boxlets (``loops'' and ``bananas''). By inverting the normalizing transformation of coordinates, these solutions are either in the form of standard truncated power series or in a rational form constructed by a continued fraction truncated at the same order of the series. Knowing the ``normal form'' approximating the system under study, the procedure of creating those solutions is straightforward and does not require any further approximation or numerics.

We show that the analytic rational solutions obtained in this way offer a degree of reliability comparable, where data are available, to those of the semi-analytic treatment based on the Prendergast-Contopoulos approach. Both loops and bananas are quite well reconstructed in shape and dimension. We extend the analysis in CS90 to check also the energy conservation along the boxlets: we find that, whereas for normal modes energy is conserved within a few percent, for loops and bananas, at this level of approximation, it is not easy to go below 10%.

The plan of the paper is as follows: in Sect. 2, we briefly recall the method to construct normal forms for the logarithmic potential and in the Appendix we outline the explicit expressions of the 1:1 and 1:2 Hamiltonian and generating function. In Sect. 3, we analyze the approximation of the major-axis orbit and, in Sects. 4 and 5, we complete the same analysis for the loop and banana families respectively. In Sect. 6 we present our conclusions.

2 Normal forms for the logarithmic potential

We investigate the dynamics of the potential

$$V=\frac12 \log\left(R^{2} + x^2 + \frac{y^2}{q^2}\right)\cdot$$ (1)

For every finite values of the ``core radius'' R, the choice R=1 involves little loss of generality. However, the singular limit $R\rightarrow0$ associated with a central density cusp is also of relevance (Miralda-Escudé & Schwarzschild 1989). With the choice R=1, the energy E may take any non-negative value. Otherwise, the singular limit is ``scale-free'' and the dynamics are the same at every energy. The parameter q provides a measure of the ``ellipticity'' of the figure and has values ranging in the interval

$$0.6 \le q \le 1.$$ (2)

In principle, lower values of q can be considered but correspond to an unphysical density distribution. Values higher than unity are included in the treatment by reversing the role of the coordinate axes.

Normal forms for the Hamiltonian system corresponding to the potential described by Eq. (1) are constructed with standard methods (Boccaletti & Pucacco 1999; Giorgilli 2002) and were used to determine the main features of its orbit structure (Belmonte et al. 2006,2007). The starting point is the series expansion of Eq. (1) about the origin

$$V = \frac12 s - \frac{1 }{2 \cdot 2} s^{2} + \frac{1 }{2 \cdot 3} s^{3} - \frac{1 }{2 \cdot 4} s^{4} + \dots$$ (3)

where

$$s = x^2 + \frac{y^2}{q^2}\cdot$$ (4)

We briefly resume the procedure to fix notations. After a scaling transformation

$$p_y \longrightarrow {\sqrt{q}} ~ p_{y}, \quad y \longrightarrow y/\sqrt{q},$$ (5)

the original Hamiltonian

$$H(p_x,p_y,x,y)= \frac{1}{2}(p_x^2+p_y^2/q) + V(s(x,y))$$ (6)

undergoes a canonical transformation to new variables P_X,P_Y,X and Y, such that

$$K(P_{X},P_{Y},X,Y)=\sum_{n=0}^{N}K_n ,$$ (7)

with the prescription (K in ``normal form'')

$$\{K_0,K\}=0.$$ (8)

In these and subsequent formulas, we adopt the convention of labeling the first term in the expansion with the index zero: in general, the ``zero order'' terms are quadratic homogeneous polynomials and terms of order n are polynomials of degree n+2. The zero order (unperturbed) Hamiltonian,

$$K_{0} \equiv H_{0} = \frac12 (P_X^2 + X^2) + \frac{1}{2q} (P_Y^2 +Y^2),$$ (9)

with ``unperturbed'' frequencies $\omega_1=1,\omega_2 = 1/q$, plays, by means of the fundamental Eq. (8), the double role of determining the specific form of the transformation and assuming the status of the second integral of motion.

The generating function of the transformation is a series of the form

$$G=G_{1}+G_{2}+ \ldots$$ (10)

and, since the procedure is based on working at each order with quantities determined at lower orders, the normalization algorithm proceeds in steps up to the ``truncation'' order N. At each step n (with $1 \le n \le N$), the series are ``upgraded'' by expressing them in the new variables found with the normalizing transformation. In the Appendix A, we describe the expression of the normal forms and the generating function that we need to apply in the following.

It is customary to refer to the series constructed in this way as Birkhoff normal forms. The presence of terms with small denominators in the expansion forbids in general their convergence. It is therefore more effective to operate with resonant normal forms (Sanders et al. 2007; Gustavson 1966), which are still non-convergent, but have the advantage of avoiding the small divisors associated with a particular resonance. To determine the primary features of the orbital structure, we therefore approximate the frequencies with a rational number plus a small ``detuning'' (Contopoulos & Moutsoulas 1966; de Zeeuw & Merritt 1983)

$$\frac{\omega_1}{\omega_2} = \frac{m_1}{m_2} + \delta.$$ (11)

We refer to a detuned (m₁:m₂) resonance, where m₁+m₂ is the order of the resonance. Each resonance allows us to describe a set of possible periodic orbits appearing in the dynamics: we have the 1:1 ``loop'', the 1:2 ``banana'', the 2:3 ``fish'', and so forth (Miralda-Escudé & Schwarzschild 1989). Each orbit, if stable, is surrounded by a family of quasi-periodic orbits usually inheriting the same name.

A conservative strategy is of truncate at the lowest order $N_{\rm min}$ adequate to convey some non-trivial information about the system. In the resonant case, it can be shown (Tuwankotta & Verhulst 2000) that the lowest order in the normal form capable of including the primary characterisitics of the m₁:m₂ resonance with double reflection symmetries is

$$N_{\rm min} = 2\times(m_1+m_2-1).$$ (12)

Using ``action-angle''-like variables $\vec{J}, \vec{\theta}$ defined through the transformation

  

$$X = \sqrt{2 J_1} \cos \theta_1,\quad Y = \sqrt{2 J_2} \cos \theta_2,$$ (13)

$$P_X = \sqrt{2 J_1} \sin \theta_1,\quad P_Y = \sqrt{2 J_2} \sin \theta_2,$$ (14)

the typical structure of the doubly-symmetric resonant normal form truncated at $N_{\rm min}$ is (Sanders et al. 2007; Contopoulos 2004)

 

$$K=m_{1} J_1+m_{2} J_2+ \sum_{k=2}^{m_1+m_2} {\cal P}^{(k)}(J_1,J_... ...m_1 m_2} J_1^{m_{2}} J_2^{m_{1}} \cos~ [2(m_{2} \theta_{1}- m_{1} \theta_{2})],$$ (15)

where ${\cal P}^{(k)}$ are homogeneous polynomials of degree k whose coefficients may depend on $\delta$ and the constant a_{m1 m2}(q) is the only marker of the resonance. In these variables, the second integral is

$${\cal E}=m_{1} J_1+m_{2} J_2$$ (16)

and the angles appear only in the resonant combination

$$\psi=m_{2} \theta_{1}- m_{1} \theta_{2}.$$ (17)

For a given resonance, these two statements remain true for arbitrary $N>N_{\rm min}$. Introducing the variable conjugate to $\psi$,

$${\cal R}=m_{2} J_1-m_{1} J_2,$$ (18)

the new Hamiltonian can be expressed in the reduced form $K({\cal R}, \psi; {\cal E},q)$, that is a family of 1-d.o.f. systems parametrized by ${\cal E}$ and $\delta$.

We are interested in the solution of the equations of motion. For a non-resonant (Birkhoff) normal form, the problem is easily solved: the coefficient a_{m1 m2} vanishes and K no longer has a term containing angles. Therefore, the $\vec{J}$ are ``true'' conserved actions and the solutions are

$$X (t) = \sqrt{2 J_1} \cos \kappa_1 t ,\quad Y (t) = \sqrt{2 J_2} \cos_, (\kappa_2 t + \theta_0),$$ (19)

where

$$\vec{\kappa} = \nabla_{\vec{J}} K$$ (20)

is the frequency vector and $\theta_0$ is a suitable phase shift.

In the resonant case instead, it is not possible to write the solutions in closed form. It is true that the dynamics described by the 1-d.o.f. Hamiltonian $K({\cal R}, \psi; {\cal E},q)$ are always integrable, but, in general, the solutions cannot be written in terms of elementary functions. However, solutions can still be written down in the case of the main periodic orbits, for which $\vec{J}, \vec{\theta}$ are true action-angle variables. There are two types of periodic orbits that can be easily identified:

1.

The normal modes for which one of the $\vec{J}$ vanishes.

2.

The periodic orbits in general position characterized by a fixed relation between the two angles, $m_{2} \theta_{1}- m_{1} \theta_{2} \equiv \theta_0$.

In both cases, it is straightforward to check that the solutions retain a form analogous to Eq. (19) with known expressions of the actions and frequencies in terms of ${\cal E}$ and q such that $\kappa_1 / \kappa_2 = m_1 / m_2$.

By using the generating function Eq. (10), the solutions in terms of standard ``physical'' coordinates can be recovered (apart from possible scaling factors) inverting the canonical transformation defined by Eqs. (A.3) and (A.4). As discussed in the Appendix, the expansion Eq. (10) is composed of even-order terms only. Since in our applications we consider the 1:1 and 1:2 symmetric resonances, we have from Eq. (12) that at most $N_{\rm min} = 4$ so that the transformation back to the physical coordinates expressed as a series of the form

$$\vec{x} (t) = \vec{x}_{1} + \vec{x}_{2} + \vec{x}_{3} +\ldots$$ (21)

is given explicitly by

   

$$\vec{x}_{1} \vec{X},$$ = (22)

$$\vec{x}_{2}$$ = 0, (23)

$$\vec{x}_{3} L_{2}(\vec{X})=\{G_{2},\vec{X}\},$$ = (24)

$$\vec{x}_{4}$$ = 0, (25)

$$\vec{x}_{5} L_{4}(\vec{X}) + {\scriptstyle\frac12} L^{2}_{2}(\vec{X})= \{G_{4},\vec{X}\} + {\scriptstyle\frac12} \{G_{2},\{G_{2},\vec{X}\}\}.$$ = (26)

We again remark that the vanishing of terms of even degree is related to the double reflection symmetry embodied in the normal form. From a knowledge of the normalized solutions Eq. (19), we can therefore construct power series approximate solutions of the equations of motion of the original system

$$\frac{{\rm d}^{2} \vec{x}}{{\rm d}t^{2}} = -\nabla_{\vec{x}} V.$$ (27)

We are investigating a non-integrable system. This implies that any perturbation approach to cope with its dynamics is deemed to fail, since it produces series that do not in general converge. On the other hand, the normal form provides us with an efficient way of constructing series with an asymptotic character: this implies that at some point we should achieve an ``optimal'' value for the expansion order $N_{\rm opt}$ (hopefully $>N_{\rm min}$) that provides the best possible result (Efthymiopoulos et al. 2004). The optimal order depends on the size of the phase-space region in which we are interested. The larger the region, the lower are the value of $N_{\rm opt}$ and the accuracy of the approximation. In galactic dynamics (in contrast to celestial mechanics) it is in general preferable to obtain an overall picture of the dynamics at the expense of extreme accuracy and truncating at $N_{\rm min}$ appears a reasonable choice. Verifying if this conjecture is tenable is an aim of the present work.

3 Axial orbits

In systems of the form of Eq. (6), the orbits along the symmetry axes are simple periodic orbits. It can be verified readily that these orbits correspond to the two normal modes for which either J₁ or J₂ vanish. If the axial orbit is stable it parents a family of ``box'' orbits. A case that is both representative and useful in galactic applications is that of the stability of the x-axis periodic orbit, the ``major-axis orbit'', if q is in the range provide by Eq. (2)). Among its possible bifurcations, the most prominent is usually that due to the 1:2 resonance between the frequency of oscillation along the orbit and that of a normal perturbation, producing the ``banana'' and ``anti-banana'' orbits (Miralda-Escudé & Schwarzschild 1989). Therefore, to derive explicit solutions for both the major-axis orbit and stable bananas (the ``pendulum-like family'' in the denomination of CS90, see Sect. 5 below), we use the 1:2 symmetric normal form.

From the expression of K reported in the Appendix, we obtain on the normal mode J₂ = 0,

$$K_{A} =2q J_1 - \frac34 q J_1^2 + \frac12 q \left(\frac53 - \frac{17}{4} q(q-1) \right) J_1^3.$$ (28)

The value of the action can be computed by using the rescaling in Eq. (A.1), i.e. K_A =2qE and inverting the series. The original ``physical'' energy E can be expressed in terms of the amplitude A of the axial orbit

$$E = \frac12 \log~ (1+A^{2}).$$ (29)

The frequency is given by the usual differentiation

$$\kappa_{1} = \frac1{2q}\frac{\partial K_{A}}{\partial J_{1}},$$ (30)

where the rescaling of the energy is accounted for in order to be able to use t as the physical time. Therefore, in the normalization variables, we have a solution of the form (19) with Y=0 and (Belmonte et al. 2007)

   

$$J_{1} = E + \frac38 E^2 + \frac{25}{192} E^3,$$ (31)

$$\kappa_{1} = 1 - \frac{3}{4} E + \frac{11}{64} E^{2}.$$ (32)

Inserting this solution into the transformation formulas of Eqs. (22-26) and exploiting the terms of the generating functions of Eqs. (A.16) and (A.17), after some computer algebra we obtain

   

$$A \cos \kappa_1 t ,$$                          x₁ = (33)

$$\frac{A^{3}}{32} \left( \cos \kappa_1 t - \cos 3 \kappa_1 t \right),$$ x₃ = (34)

$$\frac{A^{5}}{64} \left( - \frac{59}{48} \cos \kappa_1 t + \cos 3 \kappa_1 t + \frac{11}{48} \cos 5 \kappa_1 t\right).$$ x₅ = (35)

This result coincides with that obtained by Scuflaire (1995) with an independent approach based on the Poincaré-Lindstedt method and provides the explicit time evolution of an oscillation starting at rest from x(0) = A.

To evaluate the quality of the approximation, a simple method is to follow the energy variation along the solution of the true potential of Eq. (1). We therefore compute

$${\tilde E}(t) = \frac12 \left(\frac{{\rm d}x}{{\rm d}t}\right)^{2} + \frac12 \log~ (1+x(t)^{2})$$ (36)

and compare this with the value of E determined by Eq. (29) for various amplitudes. To understand the question of the optimal order we can choose two different truncations of the prediction obtained with the normal form:

  

$$x^{(3)}_{\rm NF} = x_{1}+x_{3},$$ (37)

$$x^{(5)}_{\rm NF} = x_{1}+x_{3}+x_{5}$$ (38)

and compute the quantity

$$\frac{\Delta E}{E} = \frac{{\tilde E (t)} - E}{E}\cdot$$ (39)

	Figure 1: Relative energy error along the major-axis orbit for two different truncations of the normal form at E=0.1.
Open with DEXTER

	Figure 2: Relative energy error along the major-axis orbit for two different truncations of the normal form at E=0.5.
Open with DEXTER

In Fig. 1, we plot $\Delta E/E$ for E=0.1, corresponding to an amplitude A=0.47, over a half period: the curves replicate themselves in the subsequent half period. The solid line is computed using $x^{(3)}_{\rm NF}$ and the dashed line using $x^{(5)}_{\rm NF}$; the relative error in the energy conservation is almost three times smaller with the higher truncation and as low as $0.05 \%$. From Fig. 2 however, we see that with E=0.5 (A=1.31) the situation is upset: the lower order truncation, which corresponds to the first non-zero term in the normal form, has an error at least five times smaller than that of the higher truncation. We deduce that, between the two energy levels, the optimal order decreases by two and verify that, to obtain information about an orbit 3 times larger, we must accept a relative error of a few percent. In Table 1, we list the maximum absolute energy variation over a half period for various values of E and observe that the optimal order is equal to or higher than 4 up to E=0.2: for higher values of energy the optimal order is simply 2, that is for $x^{(3)}_{\rm NF}$ we achieve the most accurate result.

Table 1: Relative energy variations along the major-axis orbit with different analytic predictions.

After achieving the optimal order, it is disappointing to neglect terms evaluated by a costly high-order computation. There are however other rules for ``summing'' divergent series that adopt all terms (Bender & Orszag 1978), such as the construction of Padè approximant. A related approach is the construction of continued fractions: successive approximants obtained by truncating the fraction at various order can provide an improvement in the asymptotic convergence with respect to the original series (Khovanskii 1963). From the normal form series of Eqs. (37) and (38) we can compute the truncated fractions

  

$$x^{(3)}_{\rm CF} = \frac{x_{1}}{1-x_{3}/x_{1}},$$ (40)

$$x^{(5)}_{\rm CF} = \frac{x_{1}}{1-\frac{x_{3}/x_{1}}{1+\frac{x_{3}^{2} - x_{1}x_{5}}{x_{1}x_{3}}}}\cdot$$ (41)

These approximations produce rational solutions and it is natural to consider a relation with the Prendergast-Contopoulos approach of CS90. By using the explicit forms of Eqs. (33)-(35), we derive

 

$$x^{(3)}_{\rm CF} = \frac{ A \cos \kappa_1 t }{1+\frac{A^{2}}{16} (1 - \cos 2 \kappa_1 t) },$$ (42)

$$x^{(5)}_{\rm CF} = A \cos \kappa_1 t \ \frac{1+A^{2} (\frac{65}{9... ...os 2 \kappa_1 t) } {1+A^{2} (\frac{59}{96} + \frac{11}{48} \cos 2 \kappa_1 t) }$$ (43)

and the expression of $x^{(3)}_{\rm CF}$ is found to have the same structure as the trial rational approximation used in CS90:

$$\frac{ {\tilde A} \cos \kappa_1 t }{1+ B \cos 2 \kappa_1 t }\cdot$$ (44)

The main difference is that in $x^{(3)}_{\rm CF}$ the parameters A and $\kappa_1$ are known in analytic form by Eqs. (29) and (32), whereas in Eq. (44), ${\tilde A}, B$ and $\kappa_1$ must be computed numerically by solving a nonlinear algebraic system obtained by inserting the trial solution into the equations of motion.

	Figure 3: Relative energy error along the major-axis orbit for two different truncations of the continued fraction at E=0.5.
Open with DEXTER

In Fig. 3, we plot the same quantities as in Fig. 2 now derived via the continued fraction truncations: the solid line is computed with $x^{(3)}_{\rm CF}$ and the dashed line with $x^{(5)}_{\rm CF}$. At this energy level, the prediction using $x^{(5)}_{\rm CF}$ begins to outperform that using $x^{(3)}_{\rm NF}$. From Table 1 we observe that the performance of $x^{(5)}_{\rm CF}$ is the optimal one when going to energies higher than 0.5 and is at least as good as that of CS90 for the same range of energy.

4 Loop orbits

As a first example of a boxlet, we consider the ``loop'' orbits for which we can use the 1:1 symmetric normal form to derive explicit solutions. For moderate ellipticities (q>0.7), loops ensue as the lowest energy bifurcation due to the 1:1 resonance between the frequency of oscillation along the short (y-axis) periodic orbit and that of a normal perturbation (Miralda-Escudé & Schwarzschild 1989). From Eq. (12), for the 1:1 resonance we have $N_{\rm min} = 2$ so that a normal form truncated at K₂ is already able to produce loops. The bifurcation curve in the (q,E)-plane starts from the point (1,0) (Scuflaire 1995; Belmonte et al. 2007) and can be expressed as the series

$$E_{\rm c}(q) = 2 \left(1-q \right) + \left(1-q \right)^2 - \frac56 \left(1-q \right)^{3}\ldots$$ (45)

if the normal form is truncated at progressively higher orders. We limit ourselves to the case q=0.9 with transition energy $E_{\rm c}(0.9)=0.21$ and investigate the analytic prediction of the theory by fixing the energy level at E=1: with suitable rescaling, these are the same values of parameters used by CS90.

In the normalization variables, we have a solution in the form of Eq. (19) with

$$X (t) = \sqrt{2 J_1} \cos \kappa_{L} t ,\quad Y (t) = \sqrt{2 J_2} \sin \kappa_{L} t ,$$ (46)

with phase shift $\theta_0=\pi/2$. The actions and frequencies can be obtained from the following procedure: starting from the normal form (A.8)-(A.9), we determine the fixed points of the reduced Hamiltonian $K({\cal R}, \psi; {\cal E},q)$ with

   

$${\cal E}=J_1+ J_2,$$ (47)

$$\psi=\theta_{1}- \theta_{2},$$ (48)

$${\cal R}=J_1- J_2.$$ (49)

The fixed point corresponding to the loop is given by $\psi=\pi/2$ and

  

$$J_{1(L)} ({\cal E},q) =\left({\cal E} + {\cal R}_{L}({\cal E},q)\right)/2,$$ (50)

$$J_{2(L)} ({\cal E},q) =\left({\cal E} - {\cal R}_{L}({\cal E},q)\right)/2,$$ (51)

where ${\cal R}_{L}({\cal E},q)$ is the solution of the algebraic equation

$$\frac{\partial K}{\partial \cal R} \left({\cal R},\frac{\pi}{2}; {\cal E},q\right)=0.$$ (52)

The frequency is then given by

$$\kappa_{L} = \frac1{q} \frac{\partial K}{\partial J_{1(L)}} \left({\cal R}_{L},\frac{\pi}{2}; {\cal E},q\right).$$ (53)

Explicit expressions of actions and frequencies are (Belmonte et al. 2007)

   

$$J_{1(L)} ({\cal E},q) \frac{(3-q) {\cal E} + 4 q (q-1)}{3 q^{2}-2q+3},$$ = (54)

$$J_{2(L)} ({\cal E},q) \frac{q(3q-1) {\cal E} - 4 q (q-1)}{3 q^{2}-2q+3},$$ = (55)

$$\kappa_{L} \frac1{q} \left( 1 - \frac{3}{4 q} {\cal E} \right).$$ = (56)

For the solutions in the physical variables, we first detrmine the transformations (22)-(24) with the generating function represented by Eq. (A.11) obtaining

    

$$\sqrt{2 J_{1(L)}} \cos \kappa_L t ,$$                           x₁ = (57)

$$\frac{1}{16} \sqrt{2 J_{1(L)}} \cos \kappa_L t$$ x₃ =

$$\times\left(7(J_{2(L)} + q J_{1(L)}) -2(2 J_{2(L)} + q J_{1(L)}) \cos 2 \kappa_L t \right),$$ (58)

$$\sqrt{2 J_{2(L)}} \sin \kappa_L t ,$$ y₁ = (59)

$$\frac{1}{16 q} \sqrt{2 J_{2(L)}} \sin \kappa_L t$$ y₃ =

$$\times \left(7(J_{2(L)} + q J_{1(L)}) + 2 ( J_{2(L)} + 2q J_{1(L)} )\cos 2 \kappa_L t \right).$$ (60)

We observe that, in the first higher-order terms both actions appear, to imply that a strong coupling exists between the degrees of freedom. The inversion of the series

$$E=\frac1{q} K_{L}({\cal E},q)=\frac1{q} K\left({\cal R}_{L}({\cal E},q),\frac{\pi}{2}; {\cal E},q\right)$$ (61)

allows us to express actions and frequencies in terms of the physical energy. However, using the exact solutions of Eqs. (54)-(55) to evaluate the functions in Eqs. (57)-(60) would produce complicated expressions that would hinder the procedure of resummation with the continuous fraction. Therefore, in analogy with the series written for the axial orbit, a separation of terms of given orders is necessary and is achieved by linearly expanding the actions in the form

  

$$J_{1(L)}=a(E-E_{\rm c}),$$ (62)

$$J_{2(L)}=b+c(E-E_{\rm c}).$$ (63)

The first expansion provides an approximate value of J₁ above the bifurcation energy and must clearly be considered to be zero for $E<E_{\rm c}$. The second expression connects at $E=E_{\rm c}$ with the corresponding expression for the normal mode. Inserting these into the solutions above and expanding in powers of $E - E_{\rm c}$, we are able to group terms according to their order. Clearly, this grouping does not affect the series themselves (namely $x^{(3)}_{\rm NF}$ and $y^{(3)}_{\rm NF}$) but does influence the computation of the truncated fractions: we obtain

  

$$x^{(3)}_{\rm CF} = 5 \sqrt{2a(E-E_{\rm c})} \ \cos \kappa_ L t \ \frac{(16 + 7 b - 4 b \cos 2 \kappa_L t)^{2}}{A_{1} + A_{2} \cos 2 \kappa_L t },$$ (64)

$$y^{(3)}_{\rm CF} = b^{3/2} \sin \kappa_L t \ \frac{(72 + 35 b + 10 b \cos 2 \kappa_L t)^{2}}{B_{1} + B_{2} \cos 2 \kappa_L t },$$ (65)

where

$$4 \left(160 + 70 b -(63a+70c)(E-E_{\rm c})\right),$$                             A₁ = (66)

$$-8 \left(20 b -(9a+20c)(E-E_{\rm c})\right),$$ A₂ = (67)

$$18\sqrt{2} \left(2 b (72+35b)\!-\!3(21ab\!+\!35cb\!+\!34c)(E\!-\!E_{\rm c})\right),$$ B₁ = (68)

$$36 \sqrt{2} b \left(10 b -3(6a+5c)(E-E_{\rm c})\right).$$ B₂ = (69)

For moderate values of the bifurcation energy corresponding to large values of q in the range (2), a simple approximation is given by the linear term in Eq. (45), $E_{\rm c}=2(1-q)$. At this level of approximation, the constants appearing in the above solutions are

$$\frac32 -q ,$$              a = (70)

$$2\left(1-q \right) = E_{\rm c},$$ b = (71)

$$\frac{q}2$$ c = (72)

and can be used to plot the orbits and compare them with numerical computations. With the choice of the parameters mentioned above, in Fig. 4 we compare a numerical computation of the loop orbit (dots) with the analytic predictions given by $\vec{x}^{(3)}_{\rm NF}$ (dashed line) and $\vec{x}^{(3)}_{\rm CF}$ (continuous line). It appears to be clear how the rational solution originating in the continued fraction truncated at order 3 is accurate overall in locating both the shape and extrema of the orbit and achieves higher accuracy than the prediction of the standard truncated series. In Fig. 5 we compare the numerical computation of the same orbit (dots) with the analytic predictions given by $\vec{x}^{(5)}_{\rm NF}$ (dashed line) and $\vec{x}^{(5)}_{\rm CF}$ (continuous line). The prediction with the higher order truncation does not appear to be more accurate in terms of the shape of the loop. However, in Fig. 6, we compare the plots of $\Delta E/E$ versus time for the two truncations (order 3, continuous line; order 5, dashed line) over a half period: the relative error is some 25% at order 5, in contrast to more than 40% with $\vec{x}^{(3)}_{\rm CF}$.

	Figure 4: An orbit of the loop family at E=1.0 for q=0.9: dots correspond to the numerical solution; the continuous line corresponds to the prediction given by the continued fraction truncated at order 3; the dashed line that provided by the normal form truncated at order 3.
Open with DEXTER

	Figure 5: The same orbit of the previous figure (dots) compared with the predictions truncated at order 5 (continued fraction, continuous line; normal form, dashed line).
Open with DEXTER

	Figure 6: Relative energy error along the loop orbit with two different truncations of the continued fraction.
Open with DEXTER

5 Pendulum-like (banana) orbits

The bifurcation of the banana orbit from the major-axis occurs along a curve in the (q,E)-plane starting from the point (1/2,0) (Scuflaire 1995; Belmonte et al. 2007). It can be expressed in terms of the series

$$E_{\rm c}(q) = 8 \left(q - \scriptstyle\frac12 \right) - \fr... ...+ \frac{268}{9} \left(q - \scriptstyle\frac12 \right)^{3}\ldots$$ (73)

As before, we adopt the same values of the parameters (up to a suitable rescaling) used in CS90: the ellipticity is q=0.6 with transition energy $E_{\rm c}(0.6)=0.76$ and the energy level is fixed at E=1.15.

In the normalization variables we have a solution in the form of Eq. (19) with

$$X (t) = \sqrt{2 J_1} \cos \kappa_{B} t ,\quad Y (t) = \sqrt{2 J_2} \cos 2 \kappa_{B} t .$$ (74)

Starting from the normal form (A.13)-(A.15), to determine actions and frequencies, we locate the fixed points of the reduced Hamiltonian $K({\cal R}, \psi; {\cal E},q)$ with

   

$${\cal E}=J_1+ 2 J_2,$$ (75)

$$\psi=2 \theta_{1}- \theta_{2},$$ (76)

$${\cal R}=2 J_1- J_2.$$ (77)

The fixed point corresponding to the banana is given by $\psi=0$ and

  

$$J_{1(B)} ({\cal E},q) =\left({\cal E} + 2 {\cal R}_{B}({\cal E},q)\right)/5,$$ (78)

$$J_{2(B)} ({\cal E},q) =\left(2{\cal E} - {\cal R}_{B}({\cal E},q)\right)/5,$$ (79)

where ${\cal R}_{B}({\cal E},q)$ is the solution of the algebraic equation

$$\frac{\partial K}{\partial \cal R} ({\cal R},0; {\cal E},q)=0.$$ (80)

The frequency is then given by

$$\kappa_{B} = \frac1{2q} \frac{\partial K}{\partial J_{1(B)}} ({\cal R}_{B},0; {\cal E},q).$$ (81)

${\cal R}_{B}({\cal E},q)$ and, as a consequence J_1(B), J_2(B) and $\kappa_{B}$ are quite cumbersome algebraic expressions involving ${\cal E}$ and q. However, simple expressions to represent orbits in the initial physical coordinates can be obtained by replacing them with some suitable series expansions. We first compute the transformations (22)-(26) using the generating function of Eqs. (A.16)-(A.17) obtaining

    

$$\sqrt{2 J_{1(B)}} \cos \kappa_B t ,$$                            x₁ = (82)

$$\frac{1}{24} \sqrt{2 J_{1(B)}} \cos \kappa_B t \big(7(2 J_{2(B)} + 3 q J_{1(B)})$$ x₃ =

$$-(4 J_{2(B)} + 6 q J_{1(B)}) \cos 2 \kappa_B t + 2 J_{2(B)} \cos 4 \kappa_B t \big),$$ (83)

$$\sqrt{2 J_{2(B)}} \sin 2 \kappa_B t ,$$ y₁ = (84)

$$\frac{1}{48 q} \sqrt{2 J_{2(B)}} \big(24 q J_{1(B)} + 6(3 J_{2(B)} + 2 q J_{1(B)}) \cos 2 \kappa_B t$$ y₃ =

$$- 8 q J_{1(B)} \cos 4 \kappa_B t - 3 J_{2(B)} \cos 6 \kappa_B t \big)$$ (85)

and analogous expressions for x₅ and y₅. In analogy with the procedure followed for the loop orbit, a separation of terms of different low orders is useful and is achieved by linearly expanding the actions in the form

$$J_{1(B)}=a_{0}+a_{1}(E-E_{\rm c}),$$ (86)

$$J_{2(B)}=b_{1}(E-E_{\rm c}).$$ (87)

Inserting these into the solutions above and expanding in powers of $E - E_{\rm c}$, we are able to group terms according to their order. Here again, this grouping does not affect the series themselves (namely $x^{(k)}_{\rm NF}$ and $y^{(k)}_{\rm NF}$, with k=3,5) but rather influences the computation of the truncated fractions, namely the expressions in Eqs. (40-41) and the analogous expressions for the y coordinate.

For moderate values of the bifurcation energy (and of orbital energy), corresponding to small values of q in the range indicated by the expression in Eq. (2), a simple approximation is given by the linear term in Eq. (73), $E_{\rm c}=8(q-1/2)$, so that

  

$$J_{1(B)}=\frac12 E + \left(q - \scriptstyle\frac12 \right) \left(4 + \frac{35}{12} E \right),$$ (88)

$$J_{2(B)}=\frac14 E - \left(q - \scriptstyle\frac12 \right) \left(2 - \frac{37}{24} E \right).$$ (89)

In this way, the expansions can be written as series of the form

  

$$x^{(3)}_{\rm NF} = \sum_{j=1}^{4} A_{1j} \cos~(2j-1)\kappa_B t +(E-E_{\rm c})\cos\kappa_B t \sum_{j=0}^{4} A_{3j} \cos 2j \ \kappa_B t,$$ (90)

$$y^{(3)}_{\rm NF} = \sum_{j=0}^{5} A_{2j} \cos 2j \ \kappa_B t +(E-E_{\rm c}) \sum_{j=0}^{5} A_{4j} \cos 2j \ \kappa_B t,$$ (91)

so that, using Eq. (40), one can also construct $\vec{x}^{(3)}_{\rm CF}$. Analogously, we may proceed with the higher-order truncations $\vec{x}^{(5)}_{\rm NF}$ from which we obtain $\vec{x}^{(5)}_{\rm CF}$. A comparison of the structure of these predictions with the rational solutions based on the Prendergast-Contopoulos approach indicates that they have the same parity in the trigonometric parts: although in the expansions of Eqs. (90)-(91) many more harmonics appear, this is clearly not necessarily an indication of higher accuracy.

	Figure 7: An orbit of the pendulum-like (banana) family at E=1.15 for q=0.6: the dots correspond to the numerical solution; the continuous lines correspond to the predictions truncated at order 3.
Open with DEXTER

With the choice of parameters mentioned above, in Fig. 7, we compare a numerical computation of the banana orbit (dots) with the analytic predictions (continuous lines) given by $\vec{x}^{(3)}_{\rm NF}$ and $\vec{x}^{(3)}_{\rm CF}$. This rational solution that originates in the continued fraction truncated at order 3, is characterized by a pair of singularities in $y^{(3)}_{\rm CF} (t)$ due to the presence of poles. However, the prediction is accurate overall in locating both the shape and extrema of the orbit and exceeds the accuracy of the prediction with the standard truncated series. In Fig. 9, we plot the corresponding $\Delta E/E$ (continuous line) over a half period: the abrupt increase in the relative error is evident at the poles of the solution.

In Fig. 8, we compare the numerical computation of the same orbit (dots) with the analytic predictions given by $\vec{x}^{(5)}_{\rm NF}$ and $\vec{x}^{(5)}_{\rm CF}$. The two predictions now almost overlap but a superior performance is achieved by the continued fraction truncation at the extrema of the orbit. In Fig. 9, we plot the corresponding $\Delta E/E$ (dashed line) over a half period: the relative error is now lower than 20%, in contrast to the 30% for $\vec{x}^{(3)}_{\rm CF}$.

	Figure 8: The same orbit as in the previous figure (dots) compared with the predictions truncated at order 5 (continuous lines): for more clarity, the y-scale has been expanded.
Open with DEXTER

A comparison with the results of CS90 is possible only in terms of the reconstruction of the shape and location of the orbit (we have used the same values of the parameters, when properly rescaled). We are able to determine that the accuracy of our analytic predictions is at least as good as that in CS90. There is no statement in CS90 about the ability of their solution in conserving energy.

6 Conclusions

We have shown how to construct approximate solutions for the main periodic orbits in the cored logarithmic potential. The guiding line has been that of exploiting normal-form expansions truncated to the first order incorporating the resonance of the given family of periodic orbits. In this way, analytic approximate solutions can be developed by a complete algorithmic procedure. Although all series are truncated to the first non-trivial orders, the solutions have a simple form only in the case of the axial orbits (normal modes). For the low-order boxlets (loops and bananas), even truncations at the first non-trivial order are cumbersome and require the use of an algebraic manipulator. However, further simplifications can be achieved if the algebraic solutions for actions and frequencies are expanded about the energy and ellipticity corresponding to the bifurcation of the family. In this case, simple expressions of the expansions can be derived, both as standard series and continued fractions.

	Figure 9: Relative energy error along the banana orbit with two different truncations of the continued fraction.
Open with DEXTER

A comparison with the rational (Prendergast-Contopoulos) approach described by Contopoulos (2004), allows us to state the following conclusions: the two methods are almost equivalent in terms of the precision of the analytic prediction when evaluated to the same order (Contopoulos & Seimenis 1990). However, the normal form perturbation expansions, even if computationally heavy to compute, are completely algorithmic and analytic at every stage, whereas the evaluation of the coefficients in the rational expansions require the numerical solution of non-linear systems. We have shown that the analytic rational solutions obtained in this way offer a degree of reliability such that both loops and bananas are well reconstructed in both shape and dimension. We have extended the analysis of Contopoulos & Seimenis (1990) to check energy conservation along the boxlets: although for normal modes energy is conserved to within a few percent, for loops and bananas, at this level of approximation, it is not straightforward to achieve an accuracy better than 10%.

On the theoretical side, the usefulness of rational solutions can be explained in terms of the more rapid convergence performance of the truncated continued fractions. They are important as a resummation method of the series expansions produced in the usual way by the normalization approach. The generality of the approach enables us to confront with higher-order resonances and correspondingly higher commensurable boxlets.

In addition to the formal and algorithmic improvements, we remark on the relevance of this work to specific problems in galactic dynamics. The study of orbits in non-axisymmetric potentials is usually performed numerically; however, an exhaustive study with conventional integration methods is computer-intensive and difficult to interpret (Touma & Tremaine 1997). The availability of simple and accurate analytical recipes can be useful in several contexts in which periodic orbits and boxlets play an important role: we mention the study of the parameter space of non-axisymmetric discs (Zhao et al. 1999; Zhao 1999) and the orbit structure around massive black holes in galactic nuclei (Sridhar & Touma 1999). Even more promising is the possibility of deriving accurate solutions for periodic orbits in the triaxial case with and without rotation, for which the analysis is still at the level of the first-order averaging method applied to the 1:1:1 resonance by de Zeeuw (1985).

Acknowledgements

We thank G. Contopoulos for inspiring our original interest in this problem.

References

• Belmonte, C. Boccaletti, D., & Pucacco, G. 2006, CMDA, 95, 101 [NASA ADS] [CrossRef] (In the text)
• Belmonte, C. Boccaletti, D., & Pucacco, G. 2007, ApJ, 669, 202 [NASA ADS] [CrossRef] (In the text)
• Bender, C. M., & Orszag, S. A. 1978, Advanced Mathematical Methods for Scientists and Engineers (New York: McGraw-Hill) (In the text)
• Binney J., & Tremaine S. 1987, Galactic Dynamics (Princeton University Press) (In the text)
• Boccaletti D. and Pucacco G. 1999, Theory of Orbits: 2 [NASA ADS] (Berlin, Heidelberg: Springer-Verlag), 171 (In the text)
• Contopoulos, G. 2004, Order and Chaos in Dynamical Astronomy (Berlin, Heidelberg: Springer-Verlag), 92 (In the text)
• Contopoulos, G., & Moutsoulas, M. 1966, AJ, 71, 687 [NASA ADS] [CrossRef] (In the text)
• Contopoulos, G., & Seimenis, J. 1990, A&A, 227, 49 [NASA ADS] (In the text)
• de Zeeuw, T. 1985, MNRAS, 215, 731 [NASA ADS] (In the text)
• de Zeeuw, T., & Merritt, D. 1983, ApJ, 267, 571 [NASA ADS] [CrossRef] (In the text)
• Efthymiopoulos, C., Giorgilli, A., & Contopoulos, G. 2004, J. Phys. A: Math. Gen., 37, 10831 [NASA ADS] [CrossRef] (In the text)
• Gerhard, O., & Saha, P. 1991, MNRAS, 251, 449 [NASA ADS] (In the text)
• Giorgilli, A. 2002, Notes on Exponential Stability of Hamiltonian Systems, Centro di Ricerca Matematica E. De Giorgi, Pisa (In the text)
• Gustavson, F. 1966, AJ, 71, 670 [NASA ADS] [CrossRef] (In the text)
• Khovanskii, A. N. 1963, The Application of Continued Fractions and their Generalizations to Problems in Approximation Theory, Noordhoff, Groningen (In the text)
• Miralda-Escudé J., & Schwarzschild M. 1989, ApJ, 339, 752 [NASA ADS] [CrossRef] (In the text)
• Prendergast, K. 1982, in The Riemann Problem, eds. D. and G. Chudnowski (Berlin Heidelberg: Springer-Verlag), Lecture Notes Math., 925, 369 (In the text)
• Pucacco, G. Boccaletti, D., & Belmonte, C. 2008, CMDA, DOI 10.1007/s10569-008-9141-x. (In the text)
• Sanders, J. A. Verhulst, F., & Murdock, J. 2007, Averaging Methods in Nonlinear Dynamical Systems (Berlin, Heidelberg: Springer-Verlag) (In the text)
• Scuflaire, R. 1995, CMDA, 61, 261 [NASA ADS] [CrossRef] (In the text)
• Sridhar, S., & Touma, J. 1999, MNRAS, 303, 483 [NASA ADS] [CrossRef] (In the text)
• Touma, J., & Tremaine, S. 1997, MNRAS, 292, 905 [NASA ADS] (In the text)
• Tuwankotta, J. M., & Verhulst, F. 2000, SIAM J. Appl. Math. 61, 1369 (In the text)
• Yanguas, P. 2001, Nonlinearity, 14, 1 [NASA ADS] [CrossRef] (In the text)
• Zhao, H. S. 1999, MNRAS, 309, 636 [NASA ADS] [CrossRef] (In the text)
• Zhao, H. S., Carollo, C. M., & De Zeeuw, T. 1999, MNRAS, 304, 457 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

Appendix A: Resonant normal forms for the logarithmic potential

To implement the normalization algorithm, the original Hamiltonian (6) is rescaled according to

$${\cal H} := \frac{m_2 H}{\omega_2} = m_2 q H,$$ (A.1)

so that we redefine the Hamiltonian as the series

 

$${\cal H} (\vec{p},\vec{x}) = \sum_{k=0}^{\infty} {\cal H}_{k} \frac12 [m_1 (p_x^2+x^2) + m_2 (p_y^2+y^2)] + \frac12 m_2\delta (p_x^2+x^2) + \sum_{k=0}^{\infty} V_{k}(x^2,y^2) ,$$ (A.2)

where the detuning was introduced in Eq. (11) and the potential is the series expansion given by Eq. (3). The normalization is performed with the technique of the Lie transform (Gerhard & Saha 1991; Yanguas 2001). Considering a generating function G, new coordinates $\vec{P}$ and $\vec{X}$ are produced by the canonical transformation

$$(\vec{P},\vec{X}) = M_{G} (\vec{p},\vec{x}).$$ (A.3)

The Lie transform operator M_G is defined by (Boccaletti & Pucacco 1999)

$$M_{G} \equiv \sum_{k=0}^{\infty} M_k$$ (A.4)

where

$$M_0 = 1, \quad M_k = \sum_{j=1}^k \frac{j}{k} L_{G_j} M_{k-j}.$$ (A.5)

The linear differential operator L_G is defined in terms of the Poisson bracket, $L_{G}(\cdot)=\{G,\cdot\}$ and the functions G_j are terms in the expansion of the generating function. We fing that G₀=1 so that, in practice, the first term in its expansion can be neglected as in Eq. (10). The terms in the new Hamiltonian are determined through the recursive set of linear partial differential equations

$$K_n= {\cal H}_n +\sum_{j=0}^{n-1}M_{n-j}{\cal H}_j , \;\; n=1,2, \dots, N.$$ (A.6)

Solving the equation at the n-th step consists of a twofold task, to find K_n and G_n. We observe that, in view of the reflection symmetries of the Hamiltonian (A.2), the chain of Eq. (A.6) is composed only of members with even index and so the normal form and generating function are composed of even-index terms only. The unperturbed part of the Hamiltonian, ${\cal H}_0$, determines the specific form of the transformation. In fact, the new Hamiltonian K is said to be in normal form if, analogously to Eq. (8),

$$\{{\cal H}_0,K\}=0,$$ (A.7)

is satisfied. In the following formulas, we list the normal form and the generating function for the expansion of the logarithmic potential in the cases of the 1:1 and 1:2 resonances (Belmonte et al. 2007). The normal form is given in the more compact version given by using the action-angle-like variables $\vec{J}, \vec{\theta}$: the resulting expressions are in agreement with the general structure of (15). For the generating function, it is more useful to write the explicit version in standard $\vec{P},\vec{X}$ variables which, although cumbersome, is that used in the transformation back to the original $\vec{p},\vec{x}$ variables. For the 1:1 resonance, the terms of the normal form are

   

                                K₀ = J₁ + J₂, (A.8)

$$\delta J_{1} - \frac{3q}{8} J_1^2 \!+ \!\frac{3}{8q} J_2^2 \!+ \!\frac12 J_1 J_2 \!+\! \frac14 J_1 J_2 \cos(2 \theta_1 - 2 \theta_2) ,$$ K₂ = (A.9)

$$\frac{q}{4} \left(\frac53 - \frac{17q}{16} \right) J_1^3 + \frac{29}{192q^{2}} J_2^3$$ K₄ =

$$+ \frac18 \left(\frac{39}{8} - 3q \right) J_1^2 J_2 - \frac38 \left(\frac{3}{8} - \frac{1}{q} \right) J_1 J_2^2$$

$$+\frac18 \left[ \left(3\! -\! \frac{5q}{4} \right) J_1^2 J_2 - \l... ...! - \!\frac{11}{4q} \right) J_1 J_2^2 \right] \cos(2 \theta_1\! -\! 2 \theta_2)$$ (A.10)

and those of the generating function are

  

$$- \frac{3 q}{32} P_{X}^3 X - \frac{3}{32} P_{X} P_{Y}^2 X - \frac{5 q}{32} P_{X}\ X^3 - \frac{3}{32} P_{X}^2\ P_{Y} Y$$                           G₂ =

$$- \frac{3}{32 q} P_{Y}^3 Y - \frac{5}{32} P_{Y} X^2 Y - \frac{5}{32} P_{X} X Y^2 - \frac{5}{32 q} P_{Y} Y^3 ,$$ (A.11)

$$\frac{5 q}{96} P_{X}^5 X - \frac{13 q^2}{256} P_{X}^5 X + \frac{9}{128} P_{X}^3\ P_{Y}^2 X - \frac{13 q}{192} P_{X}^3\ P_{Y}^2 X$$ G₄ =

$$- \frac{13}{768} P_{X} P_{Y}^4 X \!+\! \frac{7}{384 q} P_{X} P_{Y}^4 X \!-\! \frac{19 }{768} P_{X} X Y^4$$

$$-\frac{19}{384} P_{Y} X^2 Y^3 + \frac{9}{128} P_{X} P_{Y}^2 X^3 - \frac{37 q}{384} P_{X} P_{Y}^2 X^3$$

$$+\frac{25}{192 q} P_{Y} X^2 Y^3- \frac{19 q^2}{256} P_{X} X^5 - \frac{9}{256} P_{X}^4 P_{Y} Y$$

$$- \frac{13 q}{384} P_{X}^4 P_{Y} Y- \frac{13}{384} P_{X}^2\ P_{Y}^3 Y+ \frac{7}{192 q} P_{X}^2\ P_{Y}^3 Y$$

$$+ \frac{25 }{384 q} P_{X} X Y^4 + \frac{9}{64} P_{X}^2 P_{Y} X^2 Y- \frac{5 q}{64} P_{X}^2 P_{Y} X^2 Y$$

$$- \frac{5}{384} P_{Y}^3 X^2 Y- \frac{5}{384 q} P_{Y}^3 X^2 Y+ \frac{23}{256} P_{Y} X^4 Y$$

$$- \frac{19 q}{384} P_{Y} X^4 Y+ \frac{9}{128} P_{X}^3 X Y^2- \frac{37 q}{384} P_{X}^3 X Y^2$$

$$-\frac{7}{64} P_{X} P_{Y}^2 X Y^2 + \frac{11}{64 q} P_{X} P_{Y}^2 X Y^2+ \frac{23}{128} P_{X} X^3 Y^2$$

$$-\frac{19 q}{192} P_{X} X^3 Y^2-\frac{5}{384} P_{X}^2 P_{Y} Y^3 - \frac{5}{384 q} P_{X}^2 P_{Y} Y^3$$

$$+ \frac{1}{288 q^2} P_{Y}^3 Y^3 +\frac{5 q}{36} P_{X}^3 X^3- \frac{13 q^2}{96} P_{X}^3 X^3+ \frac{11 q}{96} P_{X} X^5$$

$$+ \frac{31}{768 q^2} P_{Y} Y^5+ \frac{1}{768 q^2} P_{Y}^5 Y+\frac{\delta q}{32} \left(3 P_{X}^3 X + 5 P_{X} X^3 \right)$$

$$+\frac{\delta}{64} \left( 7 P_{X} P_{Y}^2 X - P_{X}^2 P_{Y} Y+ P_{Y} X^2 Y + 9 P_{X} X Y^2 \right).$$ (A.12)

For the 1:2 resonance, the terms of the normal form are

   

                       K₀ = J₁ + 2 J₂, (A.13)

$$2 \delta J_1 - \frac34 \left(q J_{1}^2 + \frac1{q} J_{2}^2 \right) - J_{1} J_{2} ,$$ K₂ = (A.14)

$$q \left(\frac56 - \frac{17}{16} q \right) J_{1}^3 + \frac{29}{96 q^{2}} J_{2}^3$$ K₄ =

$$+\left(\frac{13}{12} - \frac32 q \right) J_{1}^2 J_{2} - \left(\frac5{12} - \frac{3}{4q} \right) J_{1} J_{2}^2$$

$$+\frac18 q J_1^2 J_2 \cos(4 \theta_1 - 2 \theta_2)$$ (A.15)

and those of the generating function are

  

$$- \frac{3 q}{16} P_{X}^3 X - \frac13 P_{X} P_{Y}^2 X - \frac{5 q}{16} P_{X}\ X^3$$                                 G₂ =

$$+ \frac{1}{24} P_{X}^2\ P_{Y} Y -\frac{3}{32 q} P_{Y}^3 Y - \frac{7}{24} P_{Y} X^2 Y$$

$$- \frac16 P_{X} X Y^2 - \frac{5}{32 q} P_{Y} Y^3 ,$$ (A.16)

$$\frac{5 q}{48} P_{X}^5 X - \frac{13 q^2}{64} P_{X}^5 X + \frac{19}{144} P_{X}^3\ P_{Y}^2 X$$ G₄ =

$$- \frac{53 q}{192} P_{X}^3\ P_{Y}^2 X-\frac{4}{45} P_{X} P_{Y}^4 X + \frac{1}{32 q} P_{X} P_{Y}^4 X$$

$$- \frac{13 }{90} P_{X} X Y^4 - \frac{17}{2880} P_{Y} X^2 Y^3+\frac{25}{144} P_{X} P_{Y}^2 X^3$$

$$- \frac{67 q}{192} P_{X} P_{Y}^2 X^3 +\frac{65}{768 q} P_{Y} X^2 Y^3$$

$$- \frac{19 q^2}{64} P_{X} X^5-\frac{1}{144} P_{X}^4 P_{Y} Y - \frac{19 q}{384} P_{X}^4 P_{Y} Y$$

$$- \frac{67}{2880} P_{X}^2\ P_{Y}^3 Y + \frac{7}{256 q} P_{X}^2\ P_{Y}^3 Y+\frac{9 }{64 q} P_{X} X Y^4$$

$$+ \frac{1}{24} P_{X}^2 P_{Y} X^2 Y - \frac{5 q}{64} P_{X}^2 P_{Y} X^2 Y - \frac{83}{2880} P_{Y}^3 X^2 Y$$

$$-\frac{1}{256 q} P_{Y}^3 X^2 Y + \frac{17}{144} P_{Y} X^4 Y - \frac{q}{128} P_{Y} X^4 Y$$

$$+ \frac{5}{36} P_{X}^3 X Y^2 - \frac{73 q}{192} P_{X}^3 X Y^2 -\frac{13}{60} P_{X} P_{Y}^2 X Y^2$$

$$+ \frac{7}{64 q} P_{X} P_{Y}^2 X Y^2+ \frac{19}{72} P_{X} X^3 Y^2- \frac{83 q}{192} P_{X} X^3 Y^2$$

$$-\frac{73}{2880} P_{X}^2 P_{Y} Y^3 + \frac{25}{768 q} P_{X}^2 P_{Y} Y^3+ \frac{1}{288 q^2} P_{Y}^3 Y^3$$

$$+ \frac{5 q}{18} P_{X}^3 X^3 - \frac{13 q^2}{24} P_{X}^3 X^3 + \frac{11 q}{48} P_{X} X^5 + \frac{31}{768 q^2} P_{Y} Y^5$$

$$+ \frac{1}{768 q^2} P_{Y}^5 Y+ \frac{\delta q}{8} \left(3 P_{X}^3 X + 5 P_{X} X^3 \right)$$

$$+\frac{\delta}{9} \left( 2 P_{X} P_{Y}^2 X + 2 P_{X}^2 P_{Y} Y - 2 P_{Y} X^2 Y + 7 P_{X} X Y^2 \right).$$ (A.17)

Concerning these formulas, two remarks are in order:

1.

The monomials associated with the detuning (i.e. for which $\delta$ appears in the coefficients) are of degree 2 less than that of the specific term of a given series. For example, in G₄ (degree 6) they appear with degree 4. This is due to the choice to consider the detuning term in Eq. (A.2) of 2nd order in the perturbation: it can be shown (Pucacco et al. 2008) that this choice, in principle not unique, is the ``optimal'' one. Since they are present in G at order 4, they appear in K only at order 6.

2.

The normalizing variables $\vec{P},\vec{X}$ have to be considered as new canonical variables at each step of the normalization: so, for example, the $\vec{P},\vec{X}$ arguments of G₂, when truncating G at N=2, are different from the arguments of G₂ when truncating G at N=4. A notation able to represent these features could be introduced but it would be complicated and we prefer to stay with the standard practice of ignoring these subtleties. However, this observation justifies the extra Poisson brackets in the transformations of the form of Eq. (26).