{'fraction_non_alphanumeric': 0.027753894830966054, 'fraction_numerical': 0.014444406231729545, 'mean_word_length': 4.623283697543996, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of axial orbits in galactic potentials

We investigate the dynamics in a galactic potential with two reflection symmetries. The phase-space structure of the real system is approximated with a resonant detuned normal form constructed with the method based on the Lie transform. Attention is focused on the stability properties of the axial periodic orbits that play an important role in galactic models. Using energy and ellipticity as parameters, we find analytical expressions of bifurcations and compare them with numerical results available in the literature.


Introduction
To determine salient features of the orbital structure of non-integrable potentials is an important topic in dynamical astronomy. Techniques based on the various versions of perturbation theory have been applied to several examples and with various degrees of approximation (for a review, see, e.g., Contopoulos, 2002). Of particular interest is to understand motion in potentials which seem to be suitable to describe elliptical galaxies. Among other features, the knowledge of the stability properties of the main periodic orbits is of paramount importance, since the bulk of density distribution is shaped by the stars in regular phasespace regions around stable periodic orbits (Binney and Tremaine, 1987). In particular, for triaxial ellipsoids, periodic orbits along symmetry axes play a special role. An enormous effort has therefore been devoted to investigate families and bifurcations of periodic orbits, starting with the study of models based on perturbed oscillators (again, for a review, Contopoulos, 2002) and gradually exploring more realistic galactic potentials with numerical (Miralda-Escudé and Schwarzschild, 1989;Fridman and Merritt, 1997) and semi-analytical (de Zeeuw and Merritt, 1983;Scuflaire, 1995) approaches.
One of the most powerful analytic tools is the normal form approximation of a non integrable system. Although the normal form approach is quite widespread in galactic dynamics, its use in studying stability of periodic orbits has not been as systematic as the theory could allow (Sanders and Verhulst, 1985). Aim of the present paper is to apply the Lie transform normalization method (Dragt and Finn, 1976;Finn, 1984;Koseleff, 1994) to approximate the dynamics of a Binney logarithmic potential (Binney and Tremaine, 1987). We compare the findings to that of Miralda-Escudé and Schwarzschild (1989), who employ purely numerical techniques to implement the Floquet method and to that of Scuflaire (1995), who studies the stability of axial orbits by solving the Hill-like perturbation equation with the Lindstedt-Poincaré approach. Another example that is briefly treated is provided by the galactic Schwarzschild (1979) potential with a comparison to the results of de Zeeuw and Merritt (1983). These authors based their approach on the averaging procedure of normalization: it is therefore interesting a comparison with that method also. We remark that, in careful numerical computations, the accuracy of predictions is usually much higher than in approximate analytical approaches. However, a reliable analytic tool is of invaluable help to gather a global overview of the behavior of the system.
To study the linear orbital stability of the main periodic orbits with a truncated normal form one can proceed in essentially two ways: the most general and exhaustive is that of determining the explicit form of the normal modes and solve the equation of their perturbation. A less general but easier approach is that of determining the nature of the fixed points on a surface of section. This is constructed with the aid of the approximate integral of motion provided by the normalization. The first method is in general quite cumbersome and can be applied when the procedure of reduction to a single degree of freedom Hamiltonian system and the use of action-angle variables lead to a reasonably simple system of equations. The second one is clearly less general but relies on simple geometric arguments related to the Hessian of a polynomial in its critical points and is, at least in principle, quite easy to implement. In this work we are going to apply both methods to perform the comparison mentioned in the paragraph above.
In galactic dynamics, the periodic orbits along the axes of symmetry (axial orbits) play a particularly important role; moreover they are easily identified both as normal modes of the reduced system and as "central" fixed points on the surfaces of section. Therefore, we will limit the detailed evaluation of the stability characteristics in the parameter space to these axial orbits. However, both procedures we have followed are quite general and can be directly applied to all periodic orbits of sufficiently low commensurability.
From the results obtained, we can state that the predictive power of the normal form ranges well outside the neighborhood in which the expansion of the original Hamiltonian is performed. It is rather related to the extent of the asymptotic convergence radius of the approximate integrals of motion. However, in concrete applications, the validity of the prediction has to be corroborated with an independent evaluation of the best suited resonant normal form of the problem at hand. A criterion for the choice is illustrated in the last section devoted to the applications and is connected to the ratio between the frequency of the periodic orbit and that of a normal perturbation to it.
The plan of the paper is as follows: in Section 2 we recall the procedure of normalization as applied to reflection symmetry potentials; in Sections 3 and 4 we study the 1:1 and 1:2 resonances respectively; in Section 5 we reconstruct approximate integrals of motion; in Section 6 we compare our analytical results with those available in the literature.

Suppose the original system is given by a Hamiltonian
with V a smooth potential with an absolute minimum and reflection symmetry with respect to both axes. In our case we will use the Schwarzschild and logarithmic potentials described below. We expand the potential up to some given degree so that and look for a new Hamiltonian given by where P, Q result from the canonical transformation (P, Q) = M g (p, q).
By expanding (3) in power series of ε and equating the coefficients of the same order, one has . . .
(5) The linear differential operator M g is defined by where the functions g n are the coefficients in the expansion of the generating function of the canonical transformation and the linear differential operator L S is defined through the Poisson bracket The exponentials in the definition of M g are intended as the formal sum of a power series so that it gives rise to a near identity coordinate transformation known as Lie series. The unperturbed part of the Hamiltonian, H 0 , determines the form of the transformation. In fact, the new Hamiltonian K is said to be in normal form if {H 0 , K} = 0.
This condition is used at each step of the procedure to determine each function g n in order to eliminate as much as possible terms in the new Hamiltonian. The only terms of which K is made of are those staying in the kernel of the operator L H 0 associated to H 0 through the definition above. The procedure is stopped at some "optimal" order and therefore in all ensuing discussion we refer to a "truncated" normal form. H 0 must be considered a function of the new coordinates at each step in the process: it is therefore an integral of the motion for the new Hamiltonian K. The function can be therefore used as a second integral of the motion conveying approximate informations on the dynamics of the original system. For practical applications (for example to compare results with numerical computations) it is useful to express approximating functions in the original physical coordinates. Inverting the coordinate transformation, the new integral of motion can be expressed in terms of the original variables. Denoting it as the power series its terms can be recovered by means of We remark that in all subsequent applications involving series expansions, the role of the perturbation parameter can also be played by the size of the neighbourhood of the origin where the Hamiltonian is considered. Therefore the powers of the parameter ε are left in all expansion formulas just to indicate their order and are treated as unity in the computations.

Galactic potentials
The model potentials we will consider are the Binney logarithmic potential (Binney and Tremaine, 1987) and the Schwarzschild (1979) potential. In both cases we will actually need an expansion of the form (2) and we will assume that each term can be written as a homogeneous polynomial of degree k + 2 of the form The logarithmic potential is and plays a very important role in galactic dynamics because, despite its very simple form, it has realistic features like a density distribution compatible at large radii with flat rotation curves. The form written here is simplified by the choice of fixing the length scale (the "core radius" R c ) equal to one, but this is not a limitation due to the invariance in both the length scale and the energy scale. With these units, the energy E may take any non negative value The parameter giving the "ellipticity" of the figure ranges in the interval 0.6 ≤ q ≤ 1.
Lower values of q can in principle be considered but correspond to a non physical density distribution. Values greater than unity are included in the treatment by reversing the role of the coordinate axes. The series expansion of the logarithmic potential is so that the lowest order coefficients are a 60 = 1, a 42 = 3/q 2 , a 24 = 3/q 4 , a 06 = 1/q 6 .
The Schwarzschild (1979) potential is to be considered more for its historical role rather than for its practical usefulness. However, it has been deeply investigated and is therefore a good benchmark for comparison. It can be written as where are two functions of r = x 2 + y 2 such that 0 < u, w < 1 and the c's are fixed constants. With the choice of de Zeeuw and Merritt (1983) 1 the lowest order coefficients are The energy range in the Schwarzschild potential is

Detuning the normal form
The natural setting in which one can perform a low order normalization is therefore that of a perturbed quadratic Hamiltonian with a potential starting with a harmonic term. In the general case in which the frequencies are rationally independent, the kernel of the operator associated to is trivial, consisting only of functions of the partial energies in the harmonic potential: it is customary to refer to the normal form constructed in this case as a "Birkhoff" normal form (Birkhoff, 1927). The presence of terms with small denominators in the expansion, forbids in general its convergence. It is therefore more effective to work since the start with a "resonant" normal form, with a "richer" kernel that allows to reconstruct the main natural resonances shaping the phasespace of the system. To catch the main features of the orbital structure, we therefore approximate the frequencies with an integer ratio plus a small "detuning" that we assume O(ε 2 ) and we speak of a detuned (m:n) resonance, with m + n the order of the resonance. We have to put the system in a form suitable to apply the normalization procedure: we rescale variables in order to put the Hamiltonian in the form where we have used the same notation for the rescaled variables and The procedure is now that of an ordinary resonant "Birkhoff-Gustavson" normalization (Gustavson, 1966;Moser, 1968) with two variants: the coordinate transformations are performed through the Lie series and the detuning quadratic term is treated as a term of higher order and put in the perturbation.

Choice of the resonance
Given an arbitrary pair of unperturbed frequencies, it could seem better to approximate their ratio as close as possible with a suitable pair of integers. However, beside possible computational problems, there are arguments on which a more effective choice can be based. Actually, the resonance should be of the lowest possible commensurability giving rise to a frequency ratio compatible with the dynamics of the actual system. The reason for this is that, the lower the order of the resonance, the richer the family of terms compatible with it that are available to construct the normal form. Moreover, another argument in favor of low order resonances comes from their role in the stability properties of periodic orbits. A typical situation is that in which a family of periodic orbits becomes unstable when a low order resonance occurs between its fundamental frequency and that of a normal perturbation: the simplest case is given by an axial orbit that, depending on the specific form of the potential, can be unstable through bifurcation of loop orbits (1:1 resonance), "banana" orbits (1:2 resonance), "fish" orbits (2:3 resonance), etcetera. Therefore a detuned low-order resonant normal form can be quite accurate in describing the corresponding bifurcations.
Finally, it must be emphasized that the structure of a resonant normal form is also affected by the symmetries of the original system. The normal form must preserve these symmetries and this in general also leads to a criterion for truncation. In the present instance of a double reflection symmetry, given a resonance ratio m/n, the normal form must contain at least terms of degree 2(m + n) (see, e.g. Tuwankotta and Verhulst, 2000). Therefore, the criterion we have adopted in this paper has been that of working with the lowest order truncated normal form incorporating the symmetries of a typical galactic potential: the 1:1 symmetric resonance which allows to truncate the normal form to degree 4 and the 1:2 symmetric resonance which requires to truncate the normal form to degree 6. A systematic investigation of the optimal order of truncation has recently been performed by Contopoulos et al. (2003) and Efthymiopoulos et al. (2004). Their results confirm the rapid decrease of the optimal order with the radius of the phase-space domain in which expansions are computed: we may conjecture that if we are interested in the global dynamics and accept a moderate level of accuracy, with this very conservative approach we can get reliable information up to the breakdown of the regular dynamics.

1:1 symmetric resonance and first order normalization
A Lie transform normalization truncated to the second order gives the following expression of the first-order normal form (cfr. Belmonte et al. 2006) where k 1 , k 2 are rational numbers dependent on m and n, Eq. (29) has been used, that in the present instance reads and the canonical variables P, Q are as in Eq. (3). We see that in this case we have the same situation as in the first order averaging approach (Verhulst, 1996): the 1:1 resonance or all other resonances.
To show more clearly how the symmetries influence the structure of the normal form, we have that it can be written as (33) where an overall rescaling by a factor ω 2 has been performed and the action-angle variables are introduced according to In fact, inverting these expressions and putting them into (33) we get The structure of (33) displays the effect of the symmetries on the resonant part: angles appear only through the combination 2θ 1 − 2θ 2 and this shows why the symmetric 1:1 resonance can also be dubbed a "2:2" resonance. We can use (33) to identify the main periodic orbits. The procedure is the following (Sanders and Verhulst, 1985, sect.7.4). We perform the following canonical transformation to "adapted resonance coordinates" In this way we get where the new action E is the additional integral of motion and has therefore been subtracted to get the effective Hamiltonian The coefficients appearing in K, give it the simplest form.
Considering the dynamics at a fixed value of E, we have that K defines a one-degree of freedom (ψ, R) system. We get the following equations of motioṅ Let us determine the fixed points of this system: these in turn give the periodic orbits of the original system. The right hand of (49) vanishes either for R = ±E or for ψ = 0, ±π. In the first case, the right hand of (48) vanishes when δ + 2 [B ± 2 (A − C(2 + cos ψ))] E = 0 (50) celmecIV.tex; 30/03/2022; 19:56; p.10 and the two periodic orbits ensue. In the second case, the right hand of (48) vanishes either when The fixed point in (53) determines the "inclined" orbit Note that 0 ≤ J 1 ≤ E and this range determines the condition for existence of the orbit of Type II. The fixed point in (54) determines the elliptic orbit The range (56) still determines the condition for existence of the orbit of Type III. Let us now consider the question of the stability of the periodic orbits. In particular, we are interested in what happens in the case of axial orbits of Type I: unfortunately, action-angle variables have singularities on these orbits and these affect also the adapted resonance coordinates. However, the remedy is quite straightforward: to use a mixed combination of action-angle variables on the normal mode and Cartesian variables for the other degree of freedom. The ensuing procedure is then first to determine the condition for the normal mode to be a critical curve of the Hamiltonian in these coordinates. Second, to assess its nature (Kummer, 1977;Contopoulos, 1978;Sanders and Verhulst, 1985, sect.7.4.4): the condition is found by considering the function where µ has to be considered as a Lagrange multiplier to take into account that there is the constraint H 0 = E. The Lagrange multiplier is found by imposing d K (µ) = 0 (59) on the normal mode. Its nature is assessed by computing the matrix of the second derivatives of K (µ) : if the Hessian determinant of the second variation is positive definite the mode is elliptic stable; if it is negative definite the mode is hyperbolic unstable.
In the case of the y-axis orbit of Eq. (52), good coordinates are so that the periodic orbit is given by The terms in the normal form are then and It is straightforward to check that, in this case, imposing Eq. (59) on the periodic orbit defined by Eq. (64), we get the equation which allows to find the required value of the Lagrange multiplier. With this result, the matrix of the second derivatives of K (µ) on the normal bundle to the periodic orbit is Recalling the rescaling in (33), the physical energy is given by If, as in the application cases that will be examined later, the first coefficient in Eq. (69) is positive, the range of instability of the y-axis orbit is Proceeding in the same way in the case of the x-axis orbit of Eq. (51), the analogous expressions 8 3 can be written for the other conditions.

1:2 symmetric resonance and second order normalization
In the case of the 1:2 resonance in presence of reflection symmetries about both axes, we know that the normal form must be pushed at least to degree 2 × (1 + 2) = 6. We therefore have to perform a further step of normalization and include K 4 in the normal form. Its expression is quite involved (cfr. Belmonte et al. 2006), but we can exploit the change of variables to action-angle coordinates to see that the normal form has the structure K = J 1 +2J 2 +ε 2 (2δJ 1 +P 2 (J 1 , J 2 ))+ε 4 (P 3 (J 1 , J 2 )+kJ 2 1 J 2 cos(4θ 1 −2θ 2 )), (74) where Eq. (29) has been used, that in the present instance reads the polynomials P 2 and P 3 are homogeneous of degree 2 and 3 respectively The study of existence and stability of periodic orbits proceeds in the same way as in the previous section: we limit the presentation to the axial orbits. In the present case they are given by The procedure to determine the condition of stability of axial orbits lead to analyze critical curves of the function (58) where K is given by (74). Considering the x-axis (type Ia, Eq. (77)) orbit, good coordinates are so that the periodic orbit is given by and The condition Eq. (59) allows to find the Lagrange multiplier The equation detK (µ) (E) = 0, obtained by computing the matrix of the second derivatives of K (µ) on the normal bundle to the periodic orbit (83), is where We see that the inequalities to be satisfied in this case are quadratic in the energy, contrary to the linear case of (72) and (73). Below we discuss the specific example of the x-axis orbit in the logarithmic potential.

Approximate integral in the original variables
The general procedure illustrated so far can be applied to study the dynamics associated to every reduced normal form. Specifically, it can be extended to forms truncated to an arbitrary high order and to all normal modes and periodic orbits associated to the given resonance. However, if the 1:1 and 1:2 cases treated here are rather simple, computations become quite cumbersome for higher orders. It can be useful to exploit an alternative, less general, but easier procedure. Recalling the generic expression of the second invariant (11), if we truncate at first order, we have This is the best approximate integral of motion of a symmetric perturbed 1:1 oscillator to order ε 4 in the perturbation parameter: in fact, due to the symmetries of the problem, odd-degree terms are absent both in the normal form and in the approximate integral. To assess the stability of, e.g., the y-axis orbit, we may proceed in this way: to use I 2 and the conserved energy to construct an x − p x Poincaré section by means of the intersection of the function I 2 (x, y, p x ; E) with the y = 0 hyperplane and the level curves of the function to study the nature of the origin as a critical point determining if it is either an extremum (elliptic fixed point = stable periodic orbit) or a saddle (hyperbolic fixed point = unstable periodic orbit) by using the Hessian determinant Clearly, for resonant periodic orbits even the location of the critical points can be already quite difficult and this limits the generality of the approach. However, in the case of axial orbits, the approach is straightforward: in the case of the y-axis orbit, we have that the second derivatives in the origin are The mixed second derivative F xpx is identically zero in the origin. The range of instability (namely, the range where the two second derivatives have different sign) is Comparing with (72) we have an order δ 2 discrepancy whose origin is connected with the construction of the function (92) where the term with the detuning is not considered as a perturbation. Analogous inequalities can be written for the x-axis orbit by constructing a y − p y Poincaré section and studying the level curves of the function obtained by means of the intersection of the function I 2 (x, y, p y ; E) with the x = 0 hyperplane. On the other side, we can develop this parallel approach of determining the nature of fixed points on the surface of section for higher order approximate integrals too. Using (11), if we truncate at second order, we have where g 2 is the first order generating function determined in the first step of the normalization. This is the best approximate integral of motion of a symmetric perturbed 1:2 oscillator to order ε 6 in the perturbation parameter. Let us consider the x-axis orbit: the construction of the y − p y Poincaré section and the study of the critical points of the function F = I 4 (0, y, p y ; E), obtained by means of the intersection of the function I 4 (x, y, p y ; E) with the x = 0 hyperplane, proceeds as above. Concerning the nature of the fixed point in the origin (linked to the stability of the x-axis orbit), we get quadratic inequalities in the energy analogous to those arising from (86): the discrepancy between the two methods is still of order δ 2 and is discussed in the subsequent section about a specific example.

Comparison of analytical and numerical results
We have applied the theory discussed in the previous sections to the logarithmic and the Schwarzschild potentials. We can compare these analytical predictions concerning the stability energy ranges of axial orbits with other results, mainly numerical, from the literature. In particular, we have chosen the work by Miralda-Escudé and Schwarzschild (1989, MES in what follows) who made an accurate numerical exploration of the phase-space structure of the logarithmic potential: they give the existence parameter ranges of the main families of periodic orbits and the bifurcation ensuing from instability thresholds in two models (corresponding to ellipticity values q = 0.7 and q = 0.9) determined by solving the perturbation equations with the Floquet method. The authors use the core radius R c to parametrize the sequence of periodic orbits because they were interested in comparing the results with the case of the singular scale-free case R c = 0: in fact they fix the energy E = 0 in all their computations and vary R c and q. For our purposes, it is more natural to use the energy as a parameter and therefore we have made the conversion to compare the results. Another approach to the study of the stability of axial orbits in the logarithmic potential has been followed by Scuflaire (1995, SC in what follows) who solved the Hill-like equation of the perturbation by means of a Poincaré-Lindstedt series expansion up to order 20. SC uses a, the amplitude of the axial orbit, as a parameter: the conversion to the energy is given by In tables I. and II. we display the results in the case of the logarithmic potential: in the first column are the two values of the parameter q for which reliable numerical estimates are available; in the second are the analytical predictions according to the procedures of Section 3 for table  I and 4 for table II; in the third are the analytical predictions according to the procedures of Section 5; in the fourth the numerical prediction of MES; in the last the semi-analytical prediction of SC. The values shown in table I. represent the energy thresholds to instability of the short y-axis orbit: this orbit becomes unstable when its frequency falls to that of a normal perturbation and as soon as this happens a loop orbit bifurcates. The values of the second and third columns are therefore computed using the results of the 1:1 resonance treatment of Sections 3 and 5, namely the first inequality in (72) and the first in (96). The computations are performed by using the coefficients given in (18) recalling definition (31). Clearly the accuracy of both predictions is better for q = 0.9 corresponding to a smaller detuning of the 1:1 resonance. Overall, we see that the prediction using the invariant curves performs significantly better than that obtained directly working with the normal form.
The long x-axis orbit can suffer instability, but not through a 1:1 resonance. This orbit becomes unstable when its frequency falls to 1/2 of that of a normal perturbation and, as soon as this happens, a banana orbit bifurcates. Therefore, we turn the attention to table II. where the values of the second and third columns are computed using the results of the 1:2 resonance treatment of Sections 4 and 5: namely, the condition in order to the quartic in (86) is negative, using the coefficients computed from (18) and (19) and the nature of the fixed point in the origin as determined through function (99). An important prediction of the numerical approaches is that, if the ellipticity is less than a given value, there is an instability interval in energy: actually in the figure of SC it is not possible to determine the upper extremum of it, but this interval is very well determined in the computations of MES. We see that the analytical estimates of the instability interval are not very far from them, with an acceptable prediction of the upper boundary of the interval too. Our treatment also provides a critical value of q c = 0.78 of the ellipticity, beyond which the interval disappears. It would be interesting to numerically test this prediction. For higher values of q the analytical estimates of the threshold of instability are not so good and the reason is that the detuning of the 1:2 resonance becomes excessively large. Finally, in table III. we report the comparisons concerning the short y-axis orbit in the Schwarzschild potential, with the coefficients listed in (24) and (25): in the third column we report the prediction given by de Zeeuw and Merritt (1983, DZM) using the first order averaging of the equations of motion and in the fourth the corresponding prediction of the same authors using a numerical integration. We see that, as usual, the analytical prediction based on the invariant is quite accurate.
We conclude observing that applications like those presented here confirm that this very conservative strategy provides sufficient qualitative and quantitative agreement with other more accurate but less general approaches. An important issue to clarify is if, in order to improve the degree of accuracy, it is enough to simply truncate the normal form to higher orders or if it is rather necessary to incorporate in a more careful way the detuning in higher order terms of the perturbation.