© Y. Wang et al. 2022

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 three-dimensional¹ 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 quasi-periodic 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 in-depth 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 self-consistent, 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 self-consistent because it is obtained from an N-body 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 disc-galaxy 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 N-body 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₂(R)cos(2θ), or, if a more complex modelling was necessary, A₂(R)cos(2θ)+A₄(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 Miyamoto-Nagai 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 X-shaped 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 near-infrared, 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 non-barred 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 high-resolution N-body 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 higher-order 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 basis-function expansion techniques (e.g. Rojas-Niñ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)². The second one, named CylSP, uses cylindrical coordinates, Bessel functions, and a spline approximation for the centre-most 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 four-fold symmetric structure, which our N-body 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)³. 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

(1)

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 Runge-Kutta 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 phase-space 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 phase-space 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 quasi-periodic, 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⁴ 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 sub-spaces, 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 phase-space region by using sub-grids and non-equidistant 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 N-body 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 phase-space 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 Schwarzschild-type 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 N-body 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⁻⁵ 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 EJ − dE and EJ + dE, where dE = 10−5. See the text for more information. The orbit number to obtain the SoS with similar EJ 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 well-known 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 low-energy 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 N-body model is fully self-consistent, 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 N-body 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 quasi-periodic 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 quasi-periodic 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 centre-most 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 centre-most 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 EJ 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 2As, 2Bs, 2Cs, 2Ds, and 2Es, 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 zero-velocity 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 sub-panels) and the SoS (lower sub-panels) 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 sub-panel, 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 sub-panel 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 sub-panel to the (x, v_x) SoS in the lower sub-panel 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 sub-panels 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 non-zero. 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 non-symmetric 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 face-on 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 sub-panels 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 zero-velocity 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 elliptical-like, 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 four-fold 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 all-round 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 N-body snapshot, which is fully self-consistent and much more realistic than the Ferrer’s ellipsoid.

Thirdly, using an N-body 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 quasi-periodic 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 quasi-periodic. How to describe the difference and transition between periodic and quasi-periodic 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.

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- 2019-05). EA thanks the CNES for financial support. This work was granted access to the HPC resources of CINES under the allocations 2017-A0020407665 and 2018-A0040407665 attributed by GENCI (Grand Equipement National de Calcul Intensif). Centre de Calcul Intensif Aix-Marseille 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, co-founded by National Astronomical Observatories, Chinese Academy of Sciences and Alibaba Cloud.

References

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

(A.1)

where R is the cylindrical radius, h is the disc radial scale length, z₀ is the disc vertical scale thickness, and M_d is the disc mass. For the halo, the initial volume density is

(A.2)

where r is the spherical radius, M_h is the halo mass, r_h and r_c are the halo core and cut-off radii, respectively, and the constant α is given by

(A.3)

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.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, Milky-Way-like galaxies, the disc mass is M_d ≈ 5 × 10¹⁰, and the scale length is h ≈ 3.5 kpc. As a result, the corresponding velocity unit is 248 km s⁻¹ and the unit of time is 1.4 × 10⁷ 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⁻⁷. 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 face-on 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 face-on.

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

(A.4)

and

(A.5)

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 quasi-horizontal 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 EJ 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 EJ=-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 N-body 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 mmax on the CDs obtained using CylSP. The thin lines give the results calculated with mmax = 12, while the thick lines show those with mmax = 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₁ and b₂, 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₁ and b₂ 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 b1 and b2 as a function of EJ with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the EJ 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₁ and b₂ 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 EJ with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the EJ 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 EJ − dE and EJ + dE, where dE = 10−5. See the text for more information. The orbit number to obtain the SoS with similar EJ 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 EJ 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 2As, 2Bs, 2Cs, 2Ds, and 2Es, 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 zero-velocity 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 face-on.
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 EJ 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 EJ=-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 mmax on the CDs obtained using CylSP. The thin lines give the results calculated with mmax = 12, while the thick lines show those with mmax = 4. See the text for further information.
In the text

	Fig. C.1. Stability diagram for x1 orbits. We plot b1 and b2 as a function of EJ with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the EJ 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 EJ with a solid line and indicate the limits of the stability range by two horizontal dashed lines. The vertical dashed line shows the EJ at which the 2E family bifurcates from the x4 family.
In the text