Issue 
A&A
Volume 668, December 2022



Article Number  A55  
Number of page(s)  18  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202243699  
Published online  01 December 2022 
Periodic orbits of multiplicity higher than one in an Nbody barred galaxy potential
^{1}
National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Beijing 100101, PR China
email: wangyg@bao.ac.cn
^{2}
AixMarseille Univ., CNRS, CNES, LAM, UMR 7326, 13388 Marseille, France
^{3}
Research Center for Astronomy, Academy of Athens, Soranou Efessiou 4, 115 27 Athens, Greece
^{4}
Department of Astronomy, Tsinghua University, Beijing 100084, PR China
Received:
2
April
2022
Accepted:
28
September
2022
Aims. Periodic orbits (POs) have been exhaustively studied. On the contrary, to our knowledge, no complete and systematic study of highermultiplicity (ℳ) POs, that is, orbits that close after more than one revolution in phase space, exists. Here, we fill this gap and also extend the standard tools used for studies of the x1 POs to studies of higher multiplicity POs.
Methods. We adopted a multiaspect approach, using surfaces of section, stability diagrams, characteristic diagrams, studies of the shapes of individual orbits, and other properties of the POs. We modified and extended the standard tools used for ℳ = 1, to ℳ > 1 cases, allowing them to use the snapshot information more fully. Our potential is more realistic than those of most previous studies, as it is obtained directly from a snapshot of a fully selfconsistent, highresolution numerical simulation.
Results. We find five main pairs of PO families with ℳ = 2. Two of these bifurcate from the x1 family and are direct rotators, and one bifurcates from the x4 family and rotates retrograde. We suggest that the remaining families do not bifurcate, but form parts of bubbles. The POs of the x1 family have fourfold symmetry, while all the ℳ = 2 POs have only twofold symmetry, with respect to either the x or the y axis. Furthermore, two orbits of the same Jacobi constant and from families of the same PO pair are mirror images of each other. Thus, by considering them together, it is possible to achieve fourfold symmetry. All results obtained here will be used in a following paper to study the effect of including ℳ orbits in the disc. We also show that a given family can include orbits of more than one multiplicity.
Key words: galaxies: structure / galaxies: general / galaxies: evolution / Galaxy: structure / Galaxy: general / Galaxy: evolution
© Y. Wang et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Periodic orbits (POs) are orbits that close in phase space. The first PO orbits to be studied in barred galaxy potentials were those in the equatorial plane of the bar. The main family of such orbits is generally known as the x1 family (e.g. Contopoulos & Papayannopoulos 1980; Contopoulos & Grosbol 1989; Athanassoula 1992a; Contopoulos 2002; Skokos et al. 2002a,b), although in some studies it is referred to as family B (from Bar), (e.g. Athanassoula et al. 1983; Pfenniger 1984; McGough et al. 2020).
Its orbits are elongated along the bar major axis. The main retrograde family in Contopoulos’ notation is called x4 (e.g. Contopoulos 1983a). There is also the pair of x2 and x3 orbits, which are both elongated perpendicular to the bar major axis and are confined to a central region much less extended than that of the x1 family. Finally, there are families at the 3:1 radial resonance, which have a triangular shape since they have three radial oscillations for one rotation, and are at the centre of the galaxy. Studies of orbital structures in threedimensional^{1} potentials have revealed a large number of 3d orbits that bifurcate from the x1 family and have considerable vertical extensions. These are often referred to as the x1 tree and are used to explain the formation, features, and properties of boxy and peanut bulges (e.g. Skokos et al. 2002a,b; Patsis et al. 2002; Athanassoula 2016).
Periodic orbits constitute the backbone of all orbital structures (Athanassoula et al. 1983). If stable, they are surrounded by regular quasiperiodic orbits, while, if unstable, they trigger chaos in their immediate surroundings in phase space. Thus, the morphology of a bar is strongly dependent on its POs, and the next question to be addressed is which of the PO families is the main contributor to the bar and can be considered as its backbone. The first study to tackle this problem was by Athanassoula et al. (1983), who found that it had to be the x1 family, or, in their notation, family B. This was generally accepted for more than three decades, until about six years ago when three studies came more or less at the same time (Gajda et al. 2016; Valluri et al. 2016; Wang et al. 2016; Abbott et al. 2017), with a rather surprising result. By Fourier analysing the temporal evolution of the x and y Cartesian coordinates of a large number of simulation orbits, they showed that there were relatively few, if any, of the x1, x2, x3, x4, or 3:1 orbits that we presented at the beginning of this section. On the contrary, they found a very large number of orbits of a quite different and much more complex morphology. These orbits have a 3:2 ratio for their main frequencies. Thus, according to these studies, it should be these orbits, and not the x1 orbits that should be the backbone of the bar.
Higher multiplicity orbits have been witnessed in a number of previous studies (e.g. Katsanikas et al. 2011; Patsis & Katsanikas 2014a; Aumer & Schönrich 2015; Tsigaridi & Patsis 2015), but were not the focus of the study and in some cases were even totally ignored. Patsis & Athanassoula (2019, hereafter PA19) were the first exception, as they gave an extensive ‘dictionary’ of a number of families of higher multiplicity orbits, giving a fiducial orbit, its multiplicity, and some references. However, so far there have been no studies of characteristic diagrams (CDs), and only a very limited number of surfaces of section (hereafter SoSs, or SoS for singular) have been examined for their ℳ > 1 contributions (see e.g. Patsis et al. 1997; Patsis & Athanassoula 2019); in other words, there are no systematic studies of SoSs for higher multiplicities. Furthermore, none of the previous studies coupled an extensive orbital structure work with frequency analysis, or with quantitative morphology of the individual orbits. Thus, there are important gaps in our knowledge on PO families with ℳ > 1 in bar potentials, and our understanding is quite incomplete.
In this paper, we aim to fill this gap by making an indepth study of POs with ℳ = 2 in a potential obtained directly from a simulation snapshot. Since we derived the potential is obtained directly from a simulation snapshot, our model is selfconsistent, and we also have a distribution function from the simulation particles. Our first aim here is to see how tools such as the CDs and the SoSs can be extended and modified to enable us to use them for higher multiplicity orbits in general. This will also be of use for studies of higher multiplicity orbits in other problems and topics. We also introduce a more realistic potential, which is fully selfconsistent because it is obtained from an Nbody simulation.
To summarise, our main aim is twofold. The first is to extend the methods used so far, so that they can be applied to ℳ > 1 POs. The second is to apply these methods so as to provide sufficient information on ℳ > 1 orbits, and thus allow comparisons between different multiplicity orbits as backbones of various discgalaxy components, for example, bars. We use these tools on a complete sample of orbits from a more realistic barred galaxy potential. These results, together with results from further analyses, will be used in a future paper to decide whether it is the x1 family and tree, or the ℳ > 1 families that dominate in the bar region. Indeed, it is not possible to decide between the x1 and the ℳ > 1 orbits as bar building blocks if we have not studied the latter to a similar degree as the former.
The structure of this paper is as follows. In Sect. 2, we introduce the barred galaxy simulation that we use, and describe how the potential and accelerations are calculated. In Sect. 3 we present how SoSs can be used for Nbody snapshots and we calculate them for a number of energy values. These values are chosen so as to cover the range useful for the bar region. In Sect. 4 we calculate and present the CDs for the POs with multiplicity 2. We summarise and discuss our results in Sect. 5. All results found here will be used in a following paper to find the backbone(s) of bars.
2. Model
2.1. Choosing a model
The first studies of orbital structures in barred galaxies used simple analytical potentials of the form A_{2}(R)cos(2θ), or, if a more complex modelling was necessary, A_{2}(R)cos(2θ)+A_{4}(R)cos(4θ), where R is the cylindrical radius and θ the azimuthal angle. These had the big advantage of enabling analytical work if simple expressions, such as polynomials, were adopted for the A(R). Their use also made it possible to reach a number of interesting conclusions and to obtain useful insights into bar dynamics (e.g. Contopoulos & Papayannopoulos 1980; Contopoulos 1980; Athanassoula 1990). However, they also had the disadvantage of leading, in some cases, to negative surface densities.
The next step was to use Ferrer’s ellipsoids for the bar and simple analytic models for the disc and the spheroid, such as the MiyamotoNagai disc (Miyamoto & Nagai 1975) and the Plummer sphere (Plummer 1911), respectively. Although these are a big improvement with respect to previous works, they still have a number of serious drawbacks and cannot claim to be realistic. For example, the disc and spheroid stay axisymmetric even when the bar is very strong, the bar isodensities are ellipsoidal, and there is no boxy, peanut, or Xshaped bulge. For further such drawbacks, we direct the reader to Athanassoula et al. (1990), the discussion in Sect. 7.6 of Athanassoula et al. (2015), or Athanassoula (2016) and references therein.
Potentials can also be calculated from appropriate images in the infrared or nearinfrared, and have thus been used for orbital structure work or other dynamical calculations (e.g. Quillen et al. 1994; Patsis et al. 1997, 2009; Kalapotharakos et al. 2010; Fragkoudi et al. 2015, 2016, 2017; Querejeta et al. 2015). This approach gives, by definition, realistic potentials, but has the important drawback that observations give only 2d projections and, although going from 2d to 3d does not present any major problems for nonbarred galaxies, it does for bars, whose third dimension has an interesting but quite complex structure (Athanassoula 2005, 2016), introducing a number of free parameters.
Last but not least, it is possible to obtain a barred galaxy potential from a highresolution Nbody simulation. The most straightforward way of doing this is by directly using the particle orbits, as calculated from the simulation codes. This, however, introduces some difficulties since orbital structure work is much more demanding in terms of accuracy than the running of the simulation itself. For example, in the case of grid codes, there may be discontinuities of higherorder derivatives as a particle moves from one cell to another. This may not influence the simulation, but could well influence the amount of chaotic orbits.
We found that a better approach is to freeze the simulation at a given time and obtain the potential and accelerations directly from the positions of the particles at that time. The potential and accelerations can be calculated using basisfunction expansion techniques (e.g. RojasNiño et al. 2016), as described in Wang et al. (2020). The description of the simulation can be seen in our Appendix, as well as in the works of Athanassoula & Misiriotis (2002) and Athanassoula (2003, 2007).
2.2. Potential and accelerations
Several methods can be used to calculate the potential and accelerations for a frozen simulation snapshot. In Wang et al. (2020), we compared two such expansion methods for galactic potentials obtained from a density distribution given by N particles of equal mass. The first method, named HO, uses spherical coordinates and spherical harmonics, and was initially proposed by Hernquist & Ostriker (1992)^{2}. The second one, named CylSP, uses cylindrical coordinates, Bessel functions, and a spline approximation for the centremost region. We found that CylSP fares much better than HO in relevant cases. In this paper we extend this comparison to SoSs, CDs, and to calculations of individual orbits.
An analytic bar in two dimensions has a fourfold symmetric structure, which our Nbody bar does not have. To obtain it, we therefore add new particles that are just reflections of the existing particles with respect to the main axes of symmetry.
More specifically, to calculate the potential and the accelerations, we used the publicly available CylSP code from AGAMA (Vasiliev 2019). AGAMA is a publicly available software library that can be widely used in stellar dynamics. In the cylindrical coordinates (R, z, ϕ), the gravitational potential and accelerations are first computed in some fixed grids (N_{R} × N_{z}) within a fixed region R ∈ [R_{min}, R_{max}], z ∈ [z_{min}, z_{max}]. Then, the potential and accelerations in any position are obtained by the spline interpolation. The accuracy of the accelerations depends on the truncated number of azimuthal Fourier harmonic terms m_{max} (Vasiliev 2019; Wang et al. 2020).
In all calculations shown here and unless otherwise noted, we used the following values for the constants: N_{R} = 50, N_{z} = 50, R_{min} = 0.1, z_{min} = 0.1, R_{max} = 392.3, z_{max} = 387.3 and m_{max} = 12. These values can give high accuracy of the force.
Wherever useful, we made comparisons of the CylSP results to those obtained with the HO method (Hernquist & Ostriker 1992)^{3}. We adopted as maximum values for the radial expansion terms n_{max} = 16 and for the angular terms l_{max} = 12, which give a good compromise between high speed and high accuracy. More detailed discussions of the two methods and of their comparison can be found in Wang et al. (2020).
2.3. Rotating frame of reference
All calculations in this and the remaining papers of this sequence are made in a frame of reference that rotates with the bar pattern speed (i.e. where the bar is at rest). Since we are using a single snapshot from this simulation, both the potential and the bar pattern speed are constant with time, so that the system is autonomous. In this frame of reference, the Jacobi energy (or Jacobi constant) is a conserved quantity. We often refer to it as ‘the energy in the rotating system’ or even shorten it to ‘the energy’. In Cartesian coordinates, it is given by
Here x and y are the coordinates in the Cartesian frame of reference corotating with the bar with angular velocity Ω_{p},ẋ and ẏ are their time derivatives, Φ(x, y) is the potential at the position (x, y), and E_{J} is the Jacobi constant. The related equations of motion are derived from the Hamiltonian given above and can be found, for example, in Binney & Tremaine (2008), or Skokos et al. (2002a). For the orbit integration, we used a 7/8 RungeKutta algorithm and we focused on orbits in the equatorial plane of the model. For the z extension, we direct the reader to the studies of Pfenniger (1984), Portail et al. (2015), PA19, and references therein.
3. Surfaces of section
3.1. General
An orbit in a 2d configuration space can be considered as a time sequence of the four phasespace coordinates (x, y, v_{x}, v_{y}). As discussed in Sect. 2.3, in our model the potential in the frame of reference rotating with the bar pattern speed is independent of time. Thus, the energy in this frame is constant and we can eliminate one of the four phasespace coordinates, say v_{y}, and the orbit can be considered as a sequence of the coordinates (x, y, v_{x}). By considering only the intersections of the orbit with the plane y = 0, and more specifically those with v_{y} > 0, we obtain a sequence of points on the 2d (x,v_{x}) plane, which is called an SoS. Any pair of quantities from phase space can be used for an SoS. Here, we consider the intersections of the orbit with the plane y = 0 and having v_{y} > 0, that is, leading to a (x, v_{x}) SoS.
Given that an ℳ = 1 PO will close after one single turn around the centre of the galaxy, it will contribute only one point to the SoS, its invariant point. For example, in Fig. 1, the x1 orbit contributes the blue filled circle on each SoS. Similarly, POs of multiplicity ℳ = 2 close after two turns around the galactic centre, and will thus cross y = 0 twice and contribute two invariant points on the SoS. Orbits with ℳ = 3 will contribute three invariant points, and so on.
Fig. 1.
Surfaces of section in the potential calculated with CylSP and for different Jacobi energies, whose values are given in the upper right corner of each panel. The blue filled circle in each panel denotes the position of the x1 invariant point. The SoS in each panel was obtained by selecting the initial conditions of each orbit on a grid covering the available part of the SoS, i.e. using the standard method. The number of orbits included in a given SoS is shown in the bottom right corner of each panel. The red and blue boxes in each panel indicate the corresponding phase regions in Figs. 2 and 3, respectively. 
Orbits encircling a stable PO (often called quasiperiodic, or regular orbits) are represented on an SoS by a sequence of black points lying on a curve, called the invariant curve. On the contrary, truly chaotic orbits^{4} will eventually fill, more or less homogeneously, the space on the SoS that is available to them. Thus, SoSs can, at a glance, give information on whether a given orbit is periodic, regular, or truly chaotic.
Surfaces of section were first introduced by Poincare (1899) and are often used in orbital structure studies because they are a powerful tool that can provide crucial information on the orbital structure in phase space (e.g. Hénon & Heiles 1964). They are mainly used for 2d problems, but have been often extended to 3d cases by considering projections from 3d to 2d subspaces, or by using techniques that visualise the distribution of the points in the 4d space (Patsis & Zachilas 1994).
It is important to note that there is a crucial difference between the SoSs as they have been used so far, and how we propose to use them here. In previous studies, the initial conditions for the orbits were obtained by setting a grid, usually Cartesian, on the SoS and starting the orbits from all the nodes of the grid that are allowed (i.e. are within the curve of zero velocity). Depending on the resolution one needs, it is possible to take a more tight or a less tight grid, or to favour a given phasespace region by using subgrids and nonequidistant nodes. We hereafter refer to this method as the standard method.
However, in our study it is possible to improve this method because, from an Nbody simulation, we can directly obtain, for any given time t, not only the potential but also the positions and velocities of each individual particle in the region of interest (here the bar and the region surrounding it, i.e. the inner part of the disc). It is all these, and only these orbits that we use for the initial conditions of the SoS orbits, instead of the nodes of the grid. This is particularly advantageous for the analysis we have set out to perform. Namely, it allows us to see the phasespace structure as it actually is, not as it would have been had all the orbits compatible with this potential been present. Unfortunately, it is possible to have such an SoS only if the potential is issued from a simulation (e.g. Sparke & Sellwood 1987; Harsoula & Kalapotharakos 2009) or from a Schwarzschildtype method (e.g. Schwarzschild 1979; Rix et al. 1997; Thomas et al. 2004; Wang et al. 2012, 2013; Magorrian 2019). We refer to this method as the Nbody method.
Strictly speaking, an SoS, corresponds to a single E_{J} value. If we applied this to our case, however, there would be only one orbit per SoS, or none (i.e. too little information). To avoid this, we grouped orbits into very narrow bins in E_{J}. Thus, every SoS corresponds to all orbits with E_{J} − dE_{J} < E_{J, orb} < E_{J} + dE_{J}, where E_{J, orb} is the Jacobi energy of the orbit and E_{J} that of the SoS. The value of dE_{J} has to be sufficiently large so that the SoS is reasonably populated, and sufficiently small so that we do not include in a given energy bin the orbits with widely different E_{J} values. Thus, dE_{J} depends on the total number of orbits available from the simulation and on the complexity of the phase space of the problem under consideration. We found the value dE_{J} = 10^{−5} to be a good compromise for our case and we adopt it throughout this work.
To the orbits obtained from the initial conditions, we added a posteriori on each SoS a filled circle (in blue in all SoS figures given here), centred on the corresponding location of the x1 orbit invariant point on the SoS. We did this in order to be able to visualise, on each SoS, the position of the higher multiplicity orbits with respect to the x1 orbit. The reader, however, should bear in mind that the x1 orbit, and indeed any PO, is only a mathematical object, as it takes infinite precision to obtain it, and thus has zero probability of being on any SoS. So what we actually mean above is that the blue point gives us the location of the x1 orbit to within the accuracy of our calculations. This, however, is several orders of magnitude better than that of the plotting, so, after having made these remarks, we can now freely talk about positions of POs on the SoS.
3.2. Salient results from the surfaces of section
We start by examining the orbital structures from a series of 15 SoSs, which span the Jacobi energy range from −2.95 to −1.5. In Fig. 1 we use the standard initial conditions and in Fig. 2 the initial conditions directly from the simulation snapshot (which we refer to as improved initial conditions), so as to make comparisons.
Fig. 2.
Same as for Fig. 1, but now using as initial conditions only the simulation orbits having a Jacobi energy between E_{J} − dE and E_{J} + dE, where dE = 10^{−5}. See the text for more information. The orbit number to obtain the SoS with similar E_{J} is labelled in the bottom right corner in each panel. 
As expected, we find on any given SoS both islands of stability and/or chaotic sea(s). Each of these occupies a fraction of the SoS surface, depending on the value of E_{J}. We now examine the structures of the SoS and of the corresponding orbits.
3.2.1. Retrograde orbits
The first thing to note is that with the standard method (Fig. 1), we find that a considerable fraction of the allowable SoS surface is occupied by retrograde orbits. This fraction is small at the innermost parts of the galaxy, that is, deep in the potential well and at low energies, but increases very fast with increasing energy, till roughly E_{J} = −2.0, where the retrograde orbits cover a very considerable part of the SoS (see from left to right and from top to bottom in Fig. 1). Furthermore, as the retrograde orbits are mainly regular, the ratio of the area occupied by regular retrograde orbits to that of the area covered by regular direct orbits is large. Furthermore, the regular retrograde area has a considerable structure, as it involves not only the wellknown x4 orbit, but also higher multiplicity retrograde orbits, such as the ℳ = 3 orbit in the E_{J} = [ − 1.9, −1.6] region and an ℳ = 5 orbit for E_{J} = −1.5.
This relatively large area covered by retrograde orbits in the SoS was also found for potentials with a Ferrer’s bar (e.g. Athanassoula et al. 1983), but it does not imply that bars are predominantly made of stars on retrograde orbits. Indeed, if we compare the SoSs obtained with the standard initial conditions (Fig. 1) to the SoSs at the same E_{J} values, but now using the improved initial conditions, we get a very different picture (Fig. 2). In the latter there are very few, if any, retrograde orbits on the lowenergy side and then the number of retrograde orbits remains of the order of 10% of the total number in the E_{J} region [−2.85, −2.4]. At yet higher energies, there are no signs of x4 orbits, or of orbits trapped around it, and also no sign of retrograde higher multiplicity orbits. As our Nbody model is fully selfconsistent, we can assert that, despite the size of their invariant curves, there are only a few corresponding retrograde orbits, that is, little retrograde trapping in realistic bar models, and thus, presumably, in the bar region of real barred galaxies.
All the invariant curves in the retrograde part of the SoS in Fig. 2 seem to be linked to x4 POs, which are ℳ = 1 and retrograde. However, as we show in the next sections, there are also ℳ = 2 or higherℳ retrograde POs, but Fig. 2 argues that they are not encircled by any considerable number of regular orbits.
3.2.2. x1 POs and other orbital types
An approximate visual estimate of the position of the x1 invariant point can be found in most of our plots, and one would have been able to localise it even if the blue filled circles were not included in the plots.
In the SoSs with standard initial conditions (Fig. 1) and within the E_{J} range [−2.75, −2.5], the x1 invariant point is surrounded by a large number of invariant curves, which, if the density is not far from constant, argues for considerable trapping around the x1 orbit. It is, however, not easy to argue for such constancy.
For the same energy region, but with SoSs constructed with Nbody initial conditions (Fig. 2), it is not clear that this statement could be made, and if it could, it would be definitely less clear cut. Furthermore, in many of the SoSs shown in Fig. 2, x1 is surrounded by an empty space. Good examples of cases with such gaps can be seen for E_{J} < −2.75. But it should be noted that such SoSs contain relatively few orbits in total; for example, for E_{J} = −2.85 there are only eight quasiperiodic orbits in total, that is, the lowest value among those shown in Fig. 2. Such gaps cannot be seen in cases where we have better statistics, such as E_{J} = −1.6, with 38 orbits in total. On the contrary, in these cases there are a number of invariant curves near and around x1 orbits, that is, a number of orbits tightly encircling the x1 orbits.
As we increase, on the SoS, the distances from x1 orbits, we find that the quasiperiodic orbits encircling the x1 orbits cease where the higher multiplicity islands become sizeable and that this, for our case, occurs roughly at E_{J} = −2.4. At larger distances from the x1, we witness the two islands characteristic of the ℳ = 2 orbital families. They have been found in several previous works with models including a Ferrer’s bar, and their corresponding families are referred to as rm21 and rm22 in PA19, and are referred to as 2A and 2A_{s}, respectively, in this paper (Sect. 4.2).
At larger distances from the x1, we encounter a number of yet higher multiplicity orbits, such as ℳ = 5, ℳ = 8, or even ℳ = 13 orbits. The former are located in the energy range −2.4 ≤ E_{J} ≤ −2.2, similar to that of the ℳ = 2 orbits discussed above. There are also similarities with respect to the location of these islands. Namely, the higher multiplicities are, similar to the ℳ = 2 orbits, confined on the SoS in a narrow asymmetric ring around the x1 orbits. At the highest couple of energies – E_{J} = −1.5 for Fig. 1, and E_{J} = −1.6 and −1.5 for Fig. 2 – we see some dissolution of the invariant curves.
3.3. Comparing CylSP and HO surfaces of section
Comparing Figs. 2 and 3, it is clearly seen that the most important difference between the two potentials is the amount of chaos. In the potential calculated in the HO way, there is considerable chaos both at the smallest energy (E_{J} = −2.95) and the second largest (E_{J} = −1.6) energy, which is absent when we use CylSP. Moreover, with HO a considerable fraction of the invariant curves are in the dissolution phase, while this is much less pronounced for the CylSP potential. This can be understood by the difference between the CylSP and the HO in the centremost region, where the latter is poorly described and leads to much larger forces and, at some level, discontinuities, while CylSP gives a much better representation. Such problems for HO can introduce chaos for orbits crossing this centremost region (centrophylic orbits). This explains why we find more chaos and more invariant curve dissolution for the HO potential. We thus confirm what was noted in Wang et al. (2020), where this difference between the two potentials has been discussed exhaustively. Looking more closely, we can see a number of further differences between Figs. 2 and 3. They are, however, much less conspicuous and quantitative than the differences discussed above concerning chaos. For example, the retrograde orbits are, in the lower energies, more pronounced for the CylSP than for the HO potential.
Fig. 3.
Same as for Fig. 2, but in the potential calculated with the HO method. It should be noted that there are similarities, but also considerable differences between these two figures. One of the strongest differences is that, for small E_{J} values, the CylSP method shows mainly regular orbits, while the HO method shows a lot of chaos and dissolving invariant curves. See the text for further discussion. 
3.4. Periodic orbits with different multiplicities
3.4.1. Periodic orbits with a multiplicity of one
The main ℳ = 1 family of POs in barred galaxy potentials is the well known x1 family, from which bifurcate a large number of other POs. Figure 4 shows a sequence of such orbits for a number of E_{J} values. These are the same values that we used for the SoSs in Sect. 3, increasing from left to right and top to bottom. Near the galactic centre, the size of the orbits is very small and increases with an increasing Jacobi constant. At the apocentres, namely at the locations where the y is maximum, most of the orbits have cusps, a sign that the bar is quite strong (Athanassoula 1992a). At the larger energy values, the cusps disappear while the orbits become considerably less elongated. These features are in good agreement with those of the orbital shapes of strong bars in Athanassoula (1992a). Since in Athanassoula (1992a) the bar potential was modelled by a simple Ferrer’s ellipsoid, we can argue that these features are generic of the 2d orbital structure in a bar.
Fig. 4.
Shape of the x1 POs in our model potential, as a function of the Jacobi constant. The Jacobi constant increases from left to right and from top to bottom. It is given in the upper right corner of each panel. It should be noted that the major axis of the bar is along the y axis. 
3.4.2. Individual periodic orbits from the main ℳ = 2 families
We find ten different families that come in pairs; in other words, for each of the families, we find there is a second one whose orbits are symmetric to those of the first with respect to the x or the y axis. This is due to the imposed symmetry of the potential. We named these five pairs 2A/2A_{s}, 2B/2B_{s}, 2C/2C_{s}, 2D/2D_{s}, and 2E/2E_{s}, where, for example, family 2A_{s} has orbits that are symmetric to those of 2A with respect to the x axis, and 2B_{s} has orbits that are symmetric to those of 2B with respect to the y axis, and so on.
In Fig. 5, the two first rows from the top give characteristic examples of the orbits of the five families of POs and of their symmetric counterparts. From left to right, we give the pairs of families 2A/2A_{s}, 2B/2B_{s}, 2C/2C_{s} (for E_{J} = −2.4), 2D/2D_{s} (for E_{J} = −1.7), and 2E/2E_{s} (for E_{J} = −2.2). The small filled circles show the positions where the orbit intersects the x axis and the arrows show the directions of the corresponding velocity vectors. Thus, they give the positions at which the orbit will cross the (x,v_{x}) SoS and are marked in the appropriate colour in the third row of panels.
Fig. 5.
Information on ℳ = 2 POs. Uppermost row (in red) shows the (x, y) view of families 2A, 2B, 2C, 2D, and 2E, as given in the upper right corner of the corresponding panel. The red filled circles show the location where the orbit intersects the x axis and the two red arrows show the directions of the corresponding velocity vectors. Middle row shows the same information but now for their symmetric orbits, i.e. from families 2A_{s}, 2B_{s}, 2C_{s}, 2D_{s}, and 2E_{s}, respectively (in green). Bottom row: we show the corresponding SoSs and the invariant points of the orbits from the first two rows. In all cases we use red for the first family and green for its symmetric counter part. See the text for more information. The blue filled circles indicate the position of the x1 orbit invariant point. 
As already mentioned, each member of the 2A or 2A_{s} families is symmetric with respect to the x axis and each member of 2A_{s} is symmetric to the member of the 2A family with the same E_{J} value, now with respect to the y axis. This is now clearly seen by comparing the orbits of the first and second rows. From the velocity arrows in these same panels, we see that 2A orbit will contribute two invariant points to the SoS, both with v_{x} = 0. One has a very small x value, while the other is located at a considerably larger x value. The x values of the two positions at which the 2A_{s} orbit crosses the x axis coincide, while their v_{x} components have equal amplitudes and opposite signs. Thus, four points in total are contributed to the SoS (the lowest panel), two from the 2A orbit and two from the 2As orbit. However, since the two crossings of the A_{s} orbit will occur at the same x value, this family will contribute only one line to the CD, as is indeed witnessed in the upper left panel of Fig. 6, and in total there would be three characteristic curves for the 2A/2A_{s} pair.
Fig. 6.
Characteristic diagram calculated using the CylSP potential. Each panel gives information on a pair of symmetric families, identified at the top of the panel. Information on the form of the orbits is given in small inserts, as described in Sect. 4.2. We note the name of the various families in the left of the panel. We also distinguish between stable (S on the plot) and unstable orbits (U on the plot), using a full line for the former and a dashed line for the latter. We also plot the zerovelocity curve (thick dashed black line), and the characteristics of the x1 and x4 orbits (thin solid and dotted lines, respectively). The cyan line in the bottom right panel is for the 1E family. More detailed discussions for family 1E can be found in Appendix B.1. 
The individual orbits of the 2D/2D_{s} families have topological characteristics similar to those of the 2A/2A_{s} families. Each individual orbit of these families is symmetric with respect to the x axis. Thus, as expected, these families also contribute four points to the SoS, two from the 2D family and the other two from the 2D_{s} family. And, again since the two x values of the 2D_{s} family are the same, they will contribute only one line to the CD, as indeed found in the middle left panel of Fig. 6. It should also be noted that the intersections of the 2D/2D_{s} families on the SoS are considerably further from the x1 orbit than the corresponding intersections of the 2A/2A_{s} families, and that the size of the individual orbits are considerably larger.
Orbital family 2B, however, follows the above only partly. Contrary to the 2A/2A_{s} and the 2D/2D_{s} cases, each 2B orbit is symmetric with respect to the y axis and although each member of B_{s} is symmetric to the member of 2B that has the same E_{J}, the symmetry now is with respect to the x axis. Thus, although these orbits have a different symmetry from that of the 2A/2A_{s} pair, by superposing two corresponding 2B and 2B_{s} orbits, we again obtain a shape symmetric with respect to both the major and the minor axes of the bar. A further difference is that of the matching of the x values of the invariant points of these orbits to the SoS. Contrary to the two previously described families, now the x value of one of the two invariant points of the 2B family is equal to one of the two invariant points of the 2B_{s} family, and the remaining 2B invariant point will match with the remaining 2B_{s} one. Thus, the four contributions of the 2B and 2B_{s} together will define only two distinct x values, which means that there should be only two 2B lines on the CD, as can indeed be witnessed in the upper row of the middle panel of Fig. 6. Family pairs 2C/2C_{s} and 2E/2E_{s} have the same symmetries and behaviour as the 2B/2B_{s} pair, and, as was the case for the 2A/2A_{s} and the 2D/2D_{s} family pairs, their contributions are much further from the x1 orbit than those of the 2B/2B_{s} pair. In fact the sequence of increasing distances is 2B/2B_{s}, 2C/2C_{s} and 2E/2E_{s}. On the contrary, the size of the individual orbits of the three families is very similar, but that could well be due to the fact that the corresponding E_{J} values are not the same in the three cases.
The SOS points for each E_{J} shown here are produced by the initial conditions of a simulation particle. From the positions of these five families in the SOS, we can see that the 2A_{s} and 2D families are relatively popular in the bar potential, while the 2B/2B_{s}, 2C/2C_{c}, and 2E/2E_{s} families are not popular in these selected E_{J}.
3.4.3. Orbits with higher multiplicities
We viewed a number of POs in our model with multiplicities ranging from one to five, and a few POs with even higher multiplicities, and noted that, in general, higher multiplicity POs have more loops and a more complex morphology than the lower multiplicity ones. The difference gets particularly strong when we compare them to ℳ = 1 POs. Typical examples are shown in Fig. 7. For each PO we show its (x, y) view (upper subpanels) and the SoS (lower subpanels) at the E_{J} of the orbit(s). The value of the latter is given in the upper right corner of each orbit panel and the multiplicity is given in the upper left. To distinguish between POs of the same E_{J} and multiplicity, we added after the multiplicity a lower case letter, for example ‘a’, or ‘b’. In each case, the location where the orbit intersects the x axis with v_{y} > 0 is given with a red filled circle and the direction of the velocity is shown by an arrow of the same colour. The location of the intersection of the orbit with the SoS is given in the lower subpanel, also with a red filled circle.
Fig. 7.
Comparison of orbits of various multiplicities. See text. 
The upper left panel gives five x1 orbits (i.e. ℳ = 1), for E_{J} = −2.2, −1.8, −1.7, −1.6, and −1.5. As there is more than one E_{J} value involved, we cannot plot an SoS. Instead, we plot in the lower subpanel the location of the x1 orbit on the (x, v_{x}) plane for each of the E_{J} values considered. We note that these orbits are well elongated and have, at least for our model, cusps. This can set a constraint on the bar strength as the strongest of bars have loops, the intermediate ones have cusps, and the weak bars have neither of these features (Athanassoula 1992a,b).
Still in the upper row, but now in the second and third columns from the left, we give two ℳ = 2 examples from the 2A and 2B families. In this and all subsequent cases, the contribution of the orbit shown in the upper subpanel to the (x, v_{x}) SoS in the lower subpanel is shown with red filled circles. In the upper row and fourth column, we show an ℳ = 3 PO.
Similarly, the third (fifth) row shows four (three) examples of ℳ = 4 (ℳ = 5) orbits. As in previous cases, red filled circles on the corresponding SoS in the lower subpanels give the contribution of the orbit.
Even a fast perusal of Fig. 7 shows us that with increasing multiplicity, the orbits acquire more loops and their morphology becomes more complex. We note that this is quite clear also in the retrograde orbits.
4. Characteristic diagrams of the periodic orbits with a multiplicity of two
4.1. General
Characteristic diagrams provide crucial information on POs and are very often used in studies of orbital structures in bars. To start simple, we assume we work with 2d orbits that are symmetric with respect to the x axis and have ℳ = 1. Then, at the intersection of the orbit with the x axis, we have y = v_{x} = = 0 and every PO can be represented by a single point, with its energy and two coordinates (E_{J}, x, v_{y}). Thus, for a given E_{J} value, we can have for every PO its x = f(E_{J}). A continuous sequence of POs forms a family and is represented by a line, usually called a characteristic curve on the (E_{J}, x) plane, in turn usually called a characteristic diagram. This is of course the most simplified case, and we could have cases in which v_{x} is nonzero. In that case, the number of dimensions of the CD increases accordingly, and the visualisation of the characteristic curves becomes increasingly difficult to plot and to comprehend. In such cases, it is usual to project the characteristic curves on planes and in our case we have chosen the (E_{J}, x) projection, which we loosely call the characteristic diagram.
Following previous studies and for the same reasons, in this section, we also restrict our CD studies to 2d potentials. We cannot, however, restrict them to POs with v_{x} = 0, because of the nature of the higher multiplicity orbits. Indeed, the ℳ = 1 orbits, have only one intersection with the x axis, and are thus either symmetric or asymmetric, that is, they either have v_{x} = 0 or v_{x} ≠ 0. This, however, is not true for higher multiplicity orbits, which will have a number of intersections with the x axis equal to their multiplicity. So a ℳ = 2 orbit can intersect the x axis once with v_{x} = 0 and once with v_{x} ≠ 0. To define such a 2d nonsymmetric PO uniquely, we need three values, for example, E_{J}, x, and v_{x}. Thus, CDs including such POs are constituted not by lines on a plane, but by lines in a 3d volume. In most cases, however, only the projection on the (E_{J}, x) plane is shown. Here we show (E_{J}, x) projections in Fig. 6.
The information obtained from the SoSs and the information from CDs is complementary. Indeed, an SoS gives information on a single E_{J}, but on all orbits, periodic or not, regular or chaotic, and for all multiplicities and families. On the contrary, a CD includes the whole relevant E_{J} range, but includes only POs, and often focuses on a single multiplicity and/or on given families.
All POs of a given family have some common morphological characteristics that must be found and analysed before one can assert whether and how a given family can contribute to a given component or morphological feature. Such analyses were performed, for example, by Athanassoula (1992a) and Patsis et al. (2003), for the faceon projection of bars, and by Skokos et al. (2002a) and Patsis et al. (2002) for their boxy and peanut central parts. We thus need to perform such analyses for orbit families with multiplicities of two or higher, before being able to assess their role in the formation and evolution of bars. To achieve this we must first extend the CDs to higher multiplicity orbits, as there are no previous studies on this subject. We thus have to start by introducing how CDs should be defined in such cases and what information can be obtained.
As already mentioned, ℳ = 1 POs will intersect the y = 0 axis with v_{y} > 0 only once, while ℳ = 2 orbits will intersect it twice before closing, ℳ = 3 orbits thrice, and so on. Thus, a family of ℳ = 1 POs will contribute a single line to the CD, while an ℳ = 2 PO family will contribute two such lines, and an ℳ = 3 PO family three. It is important to note that these ℳ > 1 lines do not correspond to several, different families, but to one single family only. Thus, for a given E_{J} there are two or more values of x, all of which correspond to a single orbit. Showing only one line on the CD plot is not a reasonable option, since there is no primary line and no specific reason to choose one rather than another. Furthermore, it could lead to important misunderstandings if two independent studies chose to plot different lines for the same family. For this reason, and as there are no previous CD studies for ℳ > 1 orbits, we decided it would be best to keep all ℳ values and lines, and make it clear that they refer to one single family only.
4.2. Salient results from our ℳ = 2 characteristic diagrams
In Fig. 6, we show the CD for these five pairs of ℳ = 2 PO families. Contrary to what is customary for the CDs of ℳ = 1 families, we plot every pair in a separate panel because the CDs of higher multiplicity POs are considerably more complex than those of ℳ = 1 and including them all in the same panel makes it very difficult to discern any structure. The names of the families in each panel are given at the top of the corresponding panel. We also include in each panel small inserts showing individual orbits. Each insert is linked by a thin, vertical dashed line to the corresponding characteristic, so that it is clear to which E_{J} each orbit corresponds. All these subpanels have the same extent and scale, so as to allow comparisons between orbits of different families, as well as the evolution of the orbital size of a given family with E_{J} at a glance. Their scales for the x and for y axes are also the same, to allow an immediate visual estimate of the elongation of each orbit.
Still in Fig. 6, we distinguish between segments of the characteristic curves where the POs are 2d stable, which are plotted with a full line, and segments where the orbits are unstable, where we use a dashed line. We also plot the zerovelocity curve (thick black dashed line), and the characteristics of the x1 and x4 POs (thin black lines, which are solid or dotted, respectively), in order to visualise whether the family described in each panel emanates from one of these two families, or not.
As described in Sect. 2, we purposely symmetrised the density distribution with respect to the symmetry axes of the bar, so as to be able to rely on previously introduced techniques and to compare with previous results, whose potentials were, in the vast majority of cases, symmetric. We have, however, found no ℳ = 2 families whose orbits are symmetric with respect to both the x = 0 and the y = 0 axes. All ℳ = 2 POs are symmetric either with respect to the x = 0 planes, or with respect to the y = 0 one. This constitutes a major difference between ℳ = 1 POs, which are mainly ellipticallike, and ℳ = 2 POs, whose orbits have much more complex shapes. Furthermore, it should be noted that for every PO that is symmetric with respect to the x = 0 plane, there is another PO that is symmetric to the first one with respect to the y = 0 plane, so that adding the two together results in a full fourfold symmetry. Therefore, the shape obtained from the superposition of the two will be symmetric with respect to both axes, as is required to build a symmetric bar. More specifically, the orbits whose characteristic curves are shown in the left panels of Fig. 6 (i.e. orbits 2A/2A_{s} and 2D/2D_{s}) are symmetric with respect to the y = 0 line. Then the two symmetric families of each pair are symmetric with respect to the x = 0 plane. On the contrary, the orbits whose characteristics are plotted in the right panels (i.e. orbits 2B/2B_{s}, 2C/2C_{s}, and 2E/2E_{s}) are symmetric with respect to the x = 0 plane, while the two orbits on the two symmetric families are symmetric with respect to y = 0. Thus, all ℳ = 2 families come in pairs, as already hinted in Sect. 3.
The characteristic curves in Fig. 6 argue that families 2A/2A_{s} and 2B/2B_{s} bifurcate from the x1 family, in good agreement with what was already found by PA19. On the contrary, family 2E/2E_{s} bifurcates from family x4, that is, the retrograde ℳ = 1 family (see lower right panel of Fig. 6), in good agreement with the fact that this is a retrograde family. We have found no family from which the 2C/2C_{s} and 2D/2D_{s} families could bifurcate from.
In Fig. 6, we see that the 2A family bifurcates from the x1 family, being unstable, but shortly after that it becomes stable. On the contrary, the 2B family bifurcates as stable and shortly after that it becomes unstable. Thus, the 2A family is mainly stable and the 2B family is mainly unstable. The 2D family is generally unstable, while the 2C family is partly stable and partly unstable. The 2E family is stable, and the 1E family is unstable if E_{J} is larger than −1.841. For a more detailed discussion on stability, we refer the interested readers to Appendix C.
5. Summary and discussion
In this paper, the first in a sequence of papers, we give information on POs with ℳ ≥ 2, a subject on which little information has been so far published. Indeed, it would have been unfair to compare the x1 family, on which a number of papers has been written, to families with ℳ ≥ 2, to which only a few sentences have been devoted, but not a single full paper.
We studied the orbital structure in a barred galaxy potential, as did a number of previous works. Ours, however, is not ‘yet another’ such study, as it has three considerable differences from and improvements over previous such studies.
Firstly, we explored orbits not with ℳ = 1, as previous works, but with ℳ = 2 or higher, and their structures and their properties. Our study includes CDs, stability diagrams, SoSs, and studies of individual orbits and their properties. This is, to our knowledge, the first attempt at a thorough and complete study of higher multiplicity orbits and their properties. Such allround studies have so far aimed only at ℳ = 1 orbits, which leads to a partial coverage of the orbits, and may thus miss essential information.
Secondly, we did not model the bar using the by now standard Ferrer’s ellipsoid. We obtained the potential directly from an Nbody snapshot, which is fully selfconsistent and much more realistic than the Ferrer’s ellipsoid.
Thirdly, using an Nbody simulation has the further advantage that it provides us also with a distribution function, that is, full information on how the particles are distributed in phase space. This was useful for the SoS (compare Figs. 1 and 2), but will also be useful when we Fourier analyse the orbits to obtain their main frequencies.
We combined the SoSs and CDs with the stability diagrams to study the shapes of individual orbits and other properties of the POs with multiplicity larger than one. We also extended the results of Wang et al. (2020), who compared the accuracy of potential and force calculations of two expansion methods, to find which one is preferable for orbit integration in barred galaxy simulations. The first such method is that of Hernquist & Ostriker (1992, HO for short), which uses a spherical coordinate system and was built specifically for the Hernquist model (Hernquist 1990). The second one is described by Vasiliev (2019), is known as CylSP, and uses a cylindrical coordinate system and a central spline. We find that the CylSP provides a better orbit integration method. The main results of this paper can be summarised as follows:
(1) From the SOS generated by the initial conditions of the real simulated particles only, it is found that the x1 orbit is encircled by lots of quasiperiodic orbits, which indicates that x1 is still an important family in the bar structure. If we increase the distance from the x1 orbit, we find that higher multiplicity islands become sizeable if E_{J} ≥ −2.4. At larger distances from the x1 orbit, there are a number of yet higher multiplicity orbits, such as ℳ = 5 or ℳ = 8, or even ℳ = 13 at some Jacobi energies.
(2) We find five main pairs of PO families with ℳ = 2, which we named 2A/2A_{s}, 2B/2B_{s}, 2C/2C_{s}, 2D/2D_{s}, and 2E/2E_{s}. The 2A/2A_{s} and 2D/2D_{s} families are symmetric to the x axis, while the 2B/2B_{s}, 2C/2C_{s}, and 2E/2E_{s} families are symmetric to the y axis.
(3) The 2A/2A_{s} and 2B/2B_{s} families bifurcate from the x1 and their existence was already known from previous works. The next three families have not been previously discussed. The 2E/2E_{s} family bifurcates from the retrograde x4 family, and its orbits are retrograde. For the 2C/2C_{s} and 2D/2D_{s} families, we suggest that they do not bifurcate from any family, but form a bubble, as the x2 and x3 families do (for the latter, see Contopoulos 1983b,c).
(4) The orbits of the 2A family are generally stable, except for a few relatively narrow ranges of E_{J}, the two widest being [−2.5296, −2.4708] and [−1.9947,−1.8950]. On the contrary, we found that the 2B and 2D families are generally unstable, while the 2C family is partly stable and partly unstable.
(5) We find a number of POs in our model, with multiplicities ranging from one to five, and a few POs of even higher multiplicities. Higher multiplicity POs have more loops and a more complex morphology than the lower multiplicity POs.
(6) We find that the SOSs calculated by the HO and CylSP methods have the strongest differences for small E_{J} values; the CylSP method shows mainly regular orbits, while the HO method shows a lot of chaos and dissolving invariant curves. The CDs calculated by these two methods have relatively small quantitative differences, except at the lowest and highest EJ values, where the difference becomes sizable.
For the POs in the bar structure, two candidates have been put forward: the x1 POs (Athanassoula et al. 1983) and the ℳ = 2 orbits (Gajda et al. 2016; Valluri et al. 2016; Wang et al. 2016; Abbott et al. 2017). A vast literature is available for the x1 family, providing all the necessary information for our task. However, little is known for the ℳ = 2 orbits. Therefore, our first step was to fill this gap and provide at least the basic information necessary for our task. This information, as well as the techniques we have developed and used here, should be useful also for other galactic dynamics studies.
Although we have found five main pairs of PO families with ℳ = 2 and a few other PO families with ℳ larger than two, it is difficult give the precise fraction of each orbit family. The first reason is that there are only a few real POs in the simulation bar model and most orbits are quasiperiodic. How to describe the difference and transition between periodic and quasiperiodic orbits is a challenging topic. The second reason is that we need a reliable orbit classification method that is valid for the general orbits. This is precisely what we try to accomplish in this work. In future papers of this series, building on the work here, we will study the contributions of different orbit families more closely.
It should be noted that this was found to be quite adequate for modelling dark matter haloes in MilkyWaylike galaxies (Sanders et al. 2020).
We use the term ‘truly chaotic’ to distinguish it from the ‘sticky’ or ‘confined chaos’, which are discussed in Appendix C.
Acknowledgments
We thank the referee for comments and suggestions that improved the paper. We acknowledge the support by the National Key R&D Program (No. 2017YFA0402603, 2018YFA0404501), the National Science Foundation of China (Grant No.11821303, 11773034, 11761141012, 11761131004 and 11633004), the Chinese Academy of Sciences (CAS) Strategic Priority Research Program XDA15020200 and and the CAS Interdisciplinary Innovation Team (JCTD 201905). EA thanks the CNES for financial support. This work was granted access to the HPC resources of CINES under the allocations 2017A0020407665 and 2018A0040407665 attributed by GENCI (Grand Equipement National de Calcul Intensif). Centre de Calcul Intensif AixMarseille is acknowledged for granting access to its high performance computing resources. These two CPU time allocations allowed EA to run and analyse the simulation used here. A considerable fraction of this work was done while EA was a “Distinguished Guest” of Tsinghua University and she thanks the Department of Astronomy for their warm hospitality and for financial support that made this visit possible. Similarly, PP and YW thank the LAM for the financial support and hospitality that made their visit possible and fruitfull. Most of the computing was performed on the high performance computing cluster at the Information and Computing Center at National Astronomical Observatories, Chinese Academy of Sciences. This work is also supported by Astronomical Big Data Joint Research Center, cofounded by National Astronomical Observatories, Chinese Academy of Sciences and Alibaba Cloud.
References
 Abbott, C. G., Valluri, M., Shen, J., & Debattista, V. P. 2017, MNRAS, 470, 1526 [Google Scholar]
 Athanassoula, E. 1990, Ann. New York Acad. Sci., 596, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E. 1992a, MNRAS, 259, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E. 1992b, MNRAS, 259, 345 [Google Scholar]
 Athanassoula, E. 2003, MNRAS, 341, 1179 [Google Scholar]
 Athanassoula, E. 2005, MNRAS, 358, 1477 [Google Scholar]
 Athanassoula, E. 2007, MNRAS, 377, 1569 [Google Scholar]
 Athanassoula, E. 2016, in Boxy/Peanut/X Bulges, Barlenses and the Thick Part of Galactic Bars: What Are They and How Did They Form?, eds. E. Laurikainen, R. Peletier, & D. Gadotti, 418, 391 [Google Scholar]
 Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35 [Google Scholar]
 Athanassoula, E., Bienayme, O., Martinet, L., & Pfenniger, D. 1983, A&A, 127, 349 [NASA ADS] [Google Scholar]
 Athanassoula, E., Morin, S., Wozniak, H., et al. 1990, MNRAS, 245, 130 [NASA ADS] [Google Scholar]
 Athanassoula, E., Lambert, J. C., & Dehnen, W. 2005, MNRAS, 363, 496 [Google Scholar]
 Athanassoula, E., RomeroGómez, M., Bosma, A., & Masdemont, J. J. 2009, MNRAS, 400, 1706 [Google Scholar]
 Athanassoula, E., RomeroGómez, M., Bosma, A., & Masdemont, J. J. 2010, MNRAS, 407, 1433 [Google Scholar]
 Athanassoula, E., Laurikainen, E., Salo, H., & Bosma, A. 2015, MNRAS, 454, 3843 [Google Scholar]
 Aumer, M., & Schönrich, R. 2015, MNRAS, 454, 3166 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 Broucke, R. 1969, NASA Techn. Rep., 32, 1360 [Google Scholar]
 Contopoulos, G. 1980, A&A, 81, 198 [NASA ADS] [Google Scholar]
 Contopoulos, G. 1983a, ApJ, 275, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G. 1983b, Phys. D Nonlinear Phenom., 8, 142 [Google Scholar]
 Contopoulos, G. 1983c, Celestial Mech., 31, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G. 2002, Order and Chaos in Dynamical Astronomy (Berlin: Springer) [Google Scholar]
 Contopoulos, G., & Grosbol, P. 1989, A&ARv, 1, 261 [Google Scholar]
 Contopoulos, G., & Harsoula, M. 2008, Int. J. Bifurcation Chaos, 18, 2929 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G., & Magnenat, P. 1985, Celestial Mech., 37, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G., & Papayannopoulos, T. 1980, A&A, 92, 33 [NASA ADS] [Google Scholar]
 Dehnen, W. 2000, AJ, 119, 800 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 2002, J. Comput. Phys., 179, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Elmegreen, B. G., & Elmegreen, D. M. 1985, ApJ, 288, 438 [Google Scholar]
 Fragkoudi, F., Athanassoula, E., Bosma, A., & Iannuzzi, F. 2015, MNRAS, 450, 229 [Google Scholar]
 Fragkoudi, F., Athanassoula, E., & Bosma, A. 2016, MNRAS, 462, L41 [Google Scholar]
 Fragkoudi, F., Athanassoula, E., & Bosma, A. 2017, MNRAS, 466, 474 [Google Scholar]
 Gajda, G., Łokas, E. L., & Athanassoula, E. 2016, ApJ, 830, 108 [Google Scholar]
 Harsoula, M., & Kalapotharakos, C. 2009, MNRAS, 394, 1605 [Google Scholar]
 Hénon, M. 1965, Ann. Astrophys., 28, 992 [Google Scholar]
 Hénon, M., & Heiles, C. 1964, AJ, 69, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
 Hernquist, L. 1993, ApJS, 86, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375 [Google Scholar]
 Kalapotharakos, C., Patsis, P. A., & Grosbøl, P. 2010, MNRAS, 403, 83 [Google Scholar]
 Katsanikas, M., Patsis, P. A., & Pinotsis, A. D. 2011, Int. J. Bifurcation Chaos, 21, 2331 [NASA ADS] [CrossRef] [Google Scholar]
 Magorrian, J. 2019, MNRAS, 484, 1166 [Google Scholar]
 McGough, D. P., Evans, N. W., & Sanders, J. L. 2020, MNRAS, 493, 2676 [Google Scholar]
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
 Patsis, P. A., & Athanassoula, E. 2019, MNRAS, 490, 2740 [Google Scholar]
 Patsis, P. A., & Katsanikas, M. 2014a, MNRAS, 445, 3525 [Google Scholar]
 Patsis, P. A., & Katsanikas, M. 2014b, MNRAS, 445, 3546 [Google Scholar]
 Patsis, P. A., & Zachilas, L. 1994, Int. J. Bifurcation Chaos, 6, 1399 [NASA ADS] [CrossRef] [Google Scholar]
 Patsis, P. A., Athanassoula, E., & Quillen, A. C. 1997, ApJ, 483, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578 [Google Scholar]
 Patsis, P. A., Skokos, C., & Athanassoula, E. 2003, MNRAS, 342, 69 [Google Scholar]
 Patsis, P. A., Kaufmann, D. E., Gottesman, S. T., & Boonyasait, V. 2009, MNRAS, 394, 142 [Google Scholar]
 Pfenniger, D. 1984, A&A, 134, 373 [NASA ADS] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
 Poincare, H. 1899, Les méthodes nouvelles de la mécanique céleste [Google Scholar]
 Portail, M., Wegg, C., & Gerhard, O. 2015, MNRAS, 450, L66 [Google Scholar]
 Querejeta, M., Meidt, S. E., Schinnerer, E., et al. 2015, ApJS, 219, 5 [Google Scholar]
 Quillen, A. C., Frogel, J. A., & Gonzalez, R. A. 1994, ApJ, 437, 162 [Google Scholar]
 Rix, H.W., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & Carollo, C. M. 1997, ApJ, 488, 702 [NASA ADS] [CrossRef] [Google Scholar]
 RojasNiño, A., Read, J. I., Aguilar, L., & Delorme, M. 2016, MNRAS, 459, 3349 [Google Scholar]
 Sanders, J. L., Lilley, E. J., Vasiliev, E., Evans, N. W., & Erkal, D. 2020, MNRAS, 499, 4793 [Google Scholar]
 Schwarzschild, M. 1979, ApJ, 232, 236 [Google Scholar]
 Skokos, C., Patsis, P. A., & Athanassoula, E. 2002a, MNRAS, 333, 847 [Google Scholar]
 Skokos, C., Patsis, P. A., & Athanassoula, E. 2002b, MNRAS, 333, 861 [Google Scholar]
 Sparke, L. S., & Sellwood, J. A. 1987, MNRAS, 225, 653 [NASA ADS] [Google Scholar]
 Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391 [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
 Tsigaridi, L., & Patsis, P. A. 2015, MNRAS, 448, 3081 [Google Scholar]
 Valluri, M., Shen, J., Abbott, C., & Debattista, V. P. 2016, ApJ, 818, 141 [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [Google Scholar]
 Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429 [Google Scholar]
 Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437 [Google Scholar]
 Wang, Y., Athanassoula, E., & Mao, S. 2016, MNRAS, 463, 3499 [Google Scholar]
 Wang, Y., Athanassoula, E., & Mao, S. 2020, A&A, 639, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Simulation
This simulation was run by one of us (EA) with initial conditions including two components, a halo and a disc, and were generated as described in Athanassoula (2007). Both components were live, thus allowing the evolution of the bar via angular momentum redistribution within the galaxy (for a review see Athanassoula 2016). The initial volume density distribution of the disc is given by
where R is the cylindrical radius, h is the disc radial scale length, z_{0} is the disc vertical scale thickness, and M_{d} is the disc mass. For the halo, the initial volume density is
where r is the spherical radius, M_{h} is the halo mass, r_{h} and r_{c} are the halo core and cutoff radii, respectively, and the constant α is given by
where q = r_{h}/r_{c} (Hernquist 1993). We used here the standard units used in stellar dynamics (i.e. in simulations that do not include gas). They have G = M_{d} = h = 1, where M_{d} is the disc mass, h is the initial disc scale length, and G is the gravitational constant. We also took M_{h} = 5., r_{h} = 0.5, r_{c} = 10, z_{0} = 0.2 and for the initial Toomre Q parameter (Toomre 1964) Q = 1.2. To convert these to astronomical units, one needs to adopt relevant values for M_{d} and R_{d}; see for example Athanassoula & Misiriotis 2002. For standard, MilkyWaylike galaxies, the disc mass is M_{d} ≈ 5 × 10^{10}, and the scale length is h ≈ 3.5 kpc. As a result, the corresponding velocity unit is 248 km s^{−1} and the unit of time is 1.4 × 10^{7} yr. The simulation was run using gyrfalcON (Dehnen 2000); Dehnen.02. A simulation with the same initial condition radial profiles and with the same parameter values was run by Athanassoula (2003)  see model Mγ3. There is, nevertheless, a difference since for orbital structure studies, a much smoother potential is required. We thus increased very considerably the number of particles for this run.
The mass resolution of the present bar simulation is quite high, and the mass of each particle is 6.25 × 10^{−7}. The total number of particles in the simulation is 9,304,278, of which 1,600,000 particles are disc particles. The remaining particles are for the dark matter halo.
As we are studying an autonomous system, we need only one time and its corresponding snapshot. For that we used a time t = 450 (see Athanassoula 2003). Fig. A.1 displays the faceon view of the disc component at this time. The bar is rotated so that the major axis of the bar coincides with the y axis. It is seen that there is a clear bar structure in the disc centre at our chosen time. The pattern speed of the bar is calculated by the evolution of the position angle of the major axis of the bar Ω_{p} = dθ/dt.
Fig. A.1.
Isodensities of the disc particles at snapshot at our chosen time. Viewed faceon. 
In order to obtain a measure of the bar strength, we follow Athanassoula & Misiriotis (2002) to project the disc particles on the (x, y), each Fourier component is given by
and
We divide the surface into annuli of equal width ΔR = 0.1 and calculate the A_{m} and B_{m} for each annulus. The ratio of the amplitude is defined as .
The left panel of Fig. A.2 shows the projected density profile along the bar major and minor axes. As expected, the projected density decreases more quickly along the bar minor axis than that along the major axis. A comparison with observations (e.g. Elmegreen & Elmegreen 1985 ) or simulations (e.g. Athanassoula & Misiriotis 2002) shows clearly that our model is realistic. The right panel of Fig. A.2 shows the ratios of the amplitude of Fourier components of the density for the even terms. It is seen that the m = 2 component has the largest ratio of the amplitude, which is consistent with the typical characteristic, as shown in Athanassoula & Misiriotis (2002).
Fig. A.2.
Basic information on our simulated model. Left: Projected density profiles along the bar major (dashed line) and minor (solid line) axes. Right: Ratio of the amplitude of for the even terms of the mass of the disc particles. 
We project all particles in the x − y plane. When we calculate the Jacobi energy, we ignore the vertical motion, the Jacobi energy is calculated by . The equipotentials of the effective potential are shown in Fig. A.3.
Fig. A.3.
Equipotentials of the effective potential of the bar at snap t = 450. Five red stars denote the positions of five Lagrange points. 
Appendix B: More information on the characteristic diagrams
B.1. A change in multiplicity along a characteristic curve
The 2E family presents a very interesting feature: for values of E_{J} less than −2.0591, its POs have ℳ = 2, as shown in the left panels of Fig. B.1. On the contrary, for larger E_{J} values, the POs have a multiplicity of ℳ = 1. This shift is abrupt and easily understood by comparing the two panels of the upper row. Indeed, for the lower of the two energies (E_{J}=2.2), the PO intersects the y = 0 line four times, two with a positive v_{y} and two with negative (see lower left panel). It will, thus, contribute two points to the SoS and will be of a multiplicity of two (ℳ = 2). As the energy increases, the orbit also gradually changes. There is a gradual increase in size and, in particular, its quasihorizontal section will gradually shift to lower y values, so that it will intersect the y = 0 line only two times, one with a positive and the other with a negative v_{y}. Thus, this PO will have ℳ = 1 (see lower right panel).
Fig. B.1.
Change in multiplicity along the characteristic curve of the 2E family. The top row of panels shows the (x,y) view of two POs whose E_{J} values are given in the bottom right corner of the corresponding panel. The thin red horizontal lines are just to guide the eye and to show the location of y = 0. The bottom row of panels shows the location of the corresponding invariant points on the SoS. The blue filled circle in each panel corresponds to the x1 invariant point, while the red ones correspond to those of the 2E POs. It should be noted that the orbit to the left has two invariant points, i.e. an ℳ = 2, while the one on the right has a single invariant point, i.e. ℳ = 1. This multiplicity shift occurs at E_{J}=2.0591. 
Thus, as we move from lower to higher E_{J} values, the multiplicity of the orbits changes from ℳ = 2 to ℳ = 1. This shows that the multiplicity is not a property linked to a given family, but that of the individual POs that constitute it. Furthermore, it is also linked to the adopted SOS.
B.2. Comparison of HO and CylSP characteristic diagrams
As already mentioned, Wang et al. (2020) compared the adequacy of two expansion methods to calculate galactic potentials obtained from a Nbody snapshot including several components such as a halo, a disc, and a bar. They found that the CylSP method is more adequate for the potential and accelerations calculations.
Here we extend this comparison to the CDs. In Fig. B.2 we compare characteristic lines, as calculated by HO and by CylSP, for our five pairs families. We note that the difference for most cases is quite small. There are, however, cases where it is noticeable, namely for both the 2A and 2B families, two regions, one with the smallest and the other with the largest E_{J} values. Thus, E_{J} extends to somewhat larger values when calculated with the CylSP method compared to the HO method.
Fig. B.2.
Comparison of the CDs obtained with the CylSP method with those obtained with the HO method of calculating the potential and forces. See the text for further information. 
B.3. Comparison of m_{max} = 4 and m_{max} = 12 for CylSP potentials
In Fig. B.3 we make a further comparison, now between the CD as calculated using m_{max} = 4 and as calculated with m_{max} = 12, both for CylSP potentials. We find that for the x1, 2C, 2D, and 2E families, the differences range between negligible and very small. For the 2A and 2B families, however, we find a larger span in differences. For most cases there is good agreement, but for a few cases, the difference is quite large, particularly at the highest E_{J} values. Thus, for the x1 family, the difference is very small, or negligible. On the contrary, for the 2A family, the maximum E_{J} value is 1.7612 for m_{max} = 12, and shrinks to 1.8617 for m_{max} = 4. Similar differences are found for the 2B family. The m_{max} = 4 presumably gives insufficient azimuthal resolution, which is particularly necessary in the higher multiplicity orbits as their structure is quite complex.
Fig. B.3.
Effect of the value of m_{max} on the CDs obtained using CylSP. The thin lines give the results calculated with m_{max} = 12, while the thick lines show those with m_{max} = 4. See the text for further information. 
It is easy to understand why the outermost regions need to have higher values for m_{max}, as we find here. Indeed, as can be seen from Fig. 7 of Athanassoula & Misiriotis (2002), or from the left panel of Fig. 4 or 5 of Athanassoula et al. (2005), the cylindrical radius at which the amplitude of the mth Fourier component reaches a maximum increases with increasing m value. We repeated the above, but now with the potential used here, and find qualitatively the same result. Thus, the largest m values are necessary for the highest distances from the centre.
Appendix C: Stability and stability curves
The stability of a given PO plays a considerable role in the dynamics of other orbits in its immediate surroundings in phase space. Thus, stable POs are surrounded by regular orbits, which reinforce to a certain degree their structure. On the contrary, unstable orbits are linked to chaos.
The above description is a rather simplified description. It is indeed possible for unstable POs to have other orbits in their immediate neighbourhood, but only over limited times. They could thus temporarily reinforce a structure despite their instability, but ultimately generate chaos. Such chaos was initially found in barred galaxy dynamics by, for example, Athanassoula et al. (1983). It was later studied in depth in several works such as Contopoulos & Harsoula (2008), Patsis & Katsanikas (2014a,b), who named it sticky chaos, and Athanassoula et al. (2009, 2010), who named it confined chaos.
Hénon (1965) described a simple and straightforward method by which one can distinguish stable from unstable POs. It simply follows and compares orbits starting in the very close neighbourhood of the POs in question. Here we follow this idea, but using the formalism of Broucke (1969) and Contopoulos & Magnenat (1985), in order to facilitate comparison with previous works on orbital structures in a barred galaxy. The nature of a PO can then be directly obtained from two stability indices: b_{1} and b_{2}, one of which is associated with radial perturbations and the other with vertical perturbations. A PO is stable if both of its indexes are within the range [−2, 2]. For a more detailed description see Skokos et al. (2002a).
Fig. C.1 shows the stability indicators b_{1} and b_{2} as a function of the E_{J} for the x1 family. It shows that, over most of the E_{J} range, this family is stable with, nevertheless, a number of tangencies. In particular, there is a tangency with b=+2 at E_{J}=2.5296, whose location we have marked with a vertical dashed line. For our study this is the main tangency, since both the 2A and 2B families bifurcate from the x1 family at this point. Further information linked to stability is given in Fig. 6, with the help of which it is possible to link the stability with location of the orbit in space.
Fig. C.1.
Stability diagram for x1 orbits. We plot b_{1} and b_{2} as a function of E_{J} with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the E_{J} at which the A and B families bifurcate from the x1 family. 
It is impotant to note that this is not the only tangency, just the first one. It is followed by others, which get closer to one another as the E_{J} increases. The second tangency of that curve occurs at E_{J}=−1.7612 and should be the bifurcation point of one or two new ℳ = 2 families.
Figure 6 gives information on the stability of the five main pairs of families with ℳ = 2. In order to discern the locations of the various stable and unstable sections, we plot them with different line styles, namely a solid line for the stable sections and a dashed one for the unstable ones. This allows us to see where orbits are stable or unstable, both as a function of E_{J} and as a function of position in the bar. We will use this information in the next paper of this series (Athanassoula et al., in prep.). We see that the 2A family bifurcates from the x1 family, being unstable, but shortly after that it becomes stable. On the contrary, the 2B family bifurcates as stable and shortly after that it becomes unstable. Thus, the 2A family is mainly stable and the 2B family is mainly unstable. The 2D family is generally unstable, while the 2C family is partly stable and partly unstable.
Figure C.2 shows the stability indicators b_{1} and b_{2} as a function of the E_{J} for the x4 family. It shows that, over most of the E_{J} range, this family is stable with, nevertheless, instability at high energies and a tangency at E_{J} = −2.3798. This is the position where the ℳ = 2 2E family bifurcates from the x4 family, and is shown by a thin vertical black line. Further information linked to this is given in Fig. 6, with the help of which it is possible to link the stability with the value of the Jacobi constant and with the location of the orbit in space.
Fig. C.2.
Stability diagram for x4 orbits. We plot b1 and b2 as a function of E_{J} with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the E_{J} at which the 2E family bifurcates from the x4 family. 
All Figures
Fig. 1.
Surfaces of section in the potential calculated with CylSP and for different Jacobi energies, whose values are given in the upper right corner of each panel. The blue filled circle in each panel denotes the position of the x1 invariant point. The SoS in each panel was obtained by selecting the initial conditions of each orbit on a grid covering the available part of the SoS, i.e. using the standard method. The number of orbits included in a given SoS is shown in the bottom right corner of each panel. The red and blue boxes in each panel indicate the corresponding phase regions in Figs. 2 and 3, respectively. 

In the text 
Fig. 2.
Same as for Fig. 1, but now using as initial conditions only the simulation orbits having a Jacobi energy between E_{J} − dE and E_{J} + dE, where dE = 10^{−5}. See the text for more information. The orbit number to obtain the SoS with similar E_{J} is labelled in the bottom right corner in each panel. 

In the text 
Fig. 3.
Same as for Fig. 2, but in the potential calculated with the HO method. It should be noted that there are similarities, but also considerable differences between these two figures. One of the strongest differences is that, for small E_{J} values, the CylSP method shows mainly regular orbits, while the HO method shows a lot of chaos and dissolving invariant curves. See the text for further discussion. 

In the text 
Fig. 4.
Shape of the x1 POs in our model potential, as a function of the Jacobi constant. The Jacobi constant increases from left to right and from top to bottom. It is given in the upper right corner of each panel. It should be noted that the major axis of the bar is along the y axis. 

In the text 
Fig. 5.
Information on ℳ = 2 POs. Uppermost row (in red) shows the (x, y) view of families 2A, 2B, 2C, 2D, and 2E, as given in the upper right corner of the corresponding panel. The red filled circles show the location where the orbit intersects the x axis and the two red arrows show the directions of the corresponding velocity vectors. Middle row shows the same information but now for their symmetric orbits, i.e. from families 2A_{s}, 2B_{s}, 2C_{s}, 2D_{s}, and 2E_{s}, respectively (in green). Bottom row: we show the corresponding SoSs and the invariant points of the orbits from the first two rows. In all cases we use red for the first family and green for its symmetric counter part. See the text for more information. The blue filled circles indicate the position of the x1 orbit invariant point. 

In the text 
Fig. 6.
Characteristic diagram calculated using the CylSP potential. Each panel gives information on a pair of symmetric families, identified at the top of the panel. Information on the form of the orbits is given in small inserts, as described in Sect. 4.2. We note the name of the various families in the left of the panel. We also distinguish between stable (S on the plot) and unstable orbits (U on the plot), using a full line for the former and a dashed line for the latter. We also plot the zerovelocity curve (thick dashed black line), and the characteristics of the x1 and x4 orbits (thin solid and dotted lines, respectively). The cyan line in the bottom right panel is for the 1E family. More detailed discussions for family 1E can be found in Appendix B.1. 

In the text 
Fig. 7.
Comparison of orbits of various multiplicities. See text. 

In the text 
Fig. A.1.
Isodensities of the disc particles at snapshot at our chosen time. Viewed faceon. 

In the text 
Fig. A.2.
Basic information on our simulated model. Left: Projected density profiles along the bar major (dashed line) and minor (solid line) axes. Right: Ratio of the amplitude of for the even terms of the mass of the disc particles. 

In the text 
Fig. A.3.
Equipotentials of the effective potential of the bar at snap t = 450. Five red stars denote the positions of five Lagrange points. 

In the text 
Fig. B.1.
Change in multiplicity along the characteristic curve of the 2E family. The top row of panels shows the (x,y) view of two POs whose E_{J} values are given in the bottom right corner of the corresponding panel. The thin red horizontal lines are just to guide the eye and to show the location of y = 0. The bottom row of panels shows the location of the corresponding invariant points on the SoS. The blue filled circle in each panel corresponds to the x1 invariant point, while the red ones correspond to those of the 2E POs. It should be noted that the orbit to the left has two invariant points, i.e. an ℳ = 2, while the one on the right has a single invariant point, i.e. ℳ = 1. This multiplicity shift occurs at E_{J}=2.0591. 

In the text 
Fig. B.2.
Comparison of the CDs obtained with the CylSP method with those obtained with the HO method of calculating the potential and forces. See the text for further information. 

In the text 
Fig. B.3.
Effect of the value of m_{max} on the CDs obtained using CylSP. The thin lines give the results calculated with m_{max} = 12, while the thick lines show those with m_{max} = 4. See the text for further information. 

In the text 
Fig. C.1.
Stability diagram for x1 orbits. We plot b_{1} and b_{2} as a function of E_{J} with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the E_{J} at which the A and B families bifurcate from the x1 family. 

In the text 
Fig. C.2.
Stability diagram for x4 orbits. We plot b1 and b2 as a function of E_{J} with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the E_{J} at which the 2E family bifurcates from the x4 family. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.