Periodic orbits of multiplicity higher than one in an N -body barred galaxy potential

Aims. Periodic orbits (POs) have been exhaustively studied. On the contrary, to our knowledge, no complete and systematic study of higher-multiplicity ( M ) POs, that is, orbits that close after more than one revolution in phase space, exists. Here, we ﬁll 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 multi-aspect approach, using surfaces of section, stability diagrams, characteristic diagrams, studies of the shapes of individual orbits, and other properties of the POs. We modiﬁed and extended the standard tools used for M = 1, to M > 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 self-consistent, high-resolution numerical simulation. Results. We ﬁnd ﬁve main pairs of PO families with M = 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 four-fold symmetry, while all the M = 2 POs have only two-fold 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 four-fold symmetry. All results obtained here will be used in a following paper to study the e ﬀ ect of including M orbits in the disc. We also show that a given family can include orbits of more than one multiplicity.


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-dimensional1 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 M > 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 M > 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 M = 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 M > 1 POs.The second is to apply these methods so as to provide sufficient information on M > 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 M > 1 families that dominate in the bar region.Indeed, it is not possible to decide between the x1 and the M > 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.

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 ana-lytical 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. 1997Patsis et al. , 2009;;Kalapotharakos et al. 2010;Fragkoudi et al. 2015Fragkoudi et al. , 2016Fragkoudi et al. , 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(Athanassoula , 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 (2003Athanassoula ( , 2007)).

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 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 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.3and 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).

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 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.

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 3 More specifically, we used the Lars Hernquist SCF.F.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 M = 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 M = 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 M = 3 will contribute three invariant points, and so on.
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 orbits4 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 phasespace 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 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 Schwarzschild-type method (e.g.Schwarzschild 1979;Rix et al. 1997;Thomas et al. 2004;Wang et al. 2012Wang et al. , 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 −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.

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 A55, page 4 of 18 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.conditions directly from the simulation snapshot (which we refer to as improved initial conditions), so as to make comparisons.
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.

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 regu-lar 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 M = 3 orbit in the E J = [−1.9,−1.6] region and an M = 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 A55, page 5 of 18 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.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 M = 1 and retrograde.However, as we show in the next sections, there are also M = 2 or higher-M retrograde POs, but Fig. 2 argues that they are not encircled by any considerable number of regular orbits.

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.85there 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 A55, page 6 of 18 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 M = 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 M = 5, M = 8, or even M = 13 orbits.The former are located in the energy range −2.4 ≤ E J ≤ −2.2, similar to that of the M = 2 orbits discussed above.There are also similarities with respect to the location of these islands.Namely, the higher multiplicities are, similar to the M = 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.

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.

Periodic orbits with a multiplicity of one
The main M = 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 elon- gated.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.

Individual periodic orbits from the main M = 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.
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 A55, page 7 of 18 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.
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 .

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 M = 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.
The upper left panel gives five x1 orbits (i.e.M = 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 M = 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 M = 3 PO.
Similarly, the third (fifth) row shows four (three) examples of M = 4 (M = 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.

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 M = 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 M = 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 M = 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, M = 1 POs will intersect the y = 0 axis with v y > 0 only once, while M = 2 orbits will intersect it twice before closing, M = 3 orbits thrice, and so on.Thus, a family of M = 1 POs will contribute a single line to the CD, while an M = 2 PO family will contribute two such lines, and an M = 3 PO family three.It is important to note that these M > 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 M > 1 orbits, we decided it would be best to keep all M values and lines, and make it clear that they refer to one single family only.

Salient results from our M = 2 characteristic diagrams
In Fig. 6, we show the CD for these five pairs of M = 2 PO families.Contrary to what is customary for the CDs of M = 1 families, we plot every pair in a separate panel because the CDs of higher multiplicity POs are considerably more complex than those of M = 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 A55, page 9 of 18 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.A55, page 10 of 18 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 A55, page 11 of 18 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 M = 2 families whose orbits are symmetric with respect to both the x = 0 and the y = 0 axes.All M = 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 M = 1 POs, which are mainly elliptical-like, and M = 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 M = 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 M = 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.

Summary and discussion
In this paper, the first in a sequence of papers, we give information on POs with M ≥ 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 M ≥ 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 M = 1, as previous works, but with M = 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 M = 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 M = 5 or M = 8, or even M = 13 at some Jacobi energies.
(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 A55, page 12 of 18 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 M = 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 M = 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 M = 2 and a few other PO families with M 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.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 (M =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 M =1 (see lower right panel).Thus, as we move from lower to higher E J values, the multiplicity of the orbits changes from M =2 to M =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  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.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.
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.

Fig. 1 .
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.

Fig. 4 .
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.

Fig. 5 .
Fig. 5. Information on M = 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.
Fig. B.2 we compare characteristic lines, as calculated by HO and by CylSP,

Fig
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 M =2, while the one on the right has a single invariant point, i.e.M =1.This multiplicity shift occurs at E J =-2.0591.

B. 3 .
Comparison of m max = 4 and m max = 12 for CylSP potentials In Fig. B.

Fig. B. 2 .
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.

Fig. B. 3 .
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.