Issue 
A&A
Volume 685, May 2024



Article Number  A93  
Number of page(s)  19  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202349000  
Published online  08 May 2024 
Orbital analysis of stars in the nuclear stellar disc of the Milky Way
^{1}
Université Côte d’Azur, Observatoire de la Côte d’Azur, Laboratoire Lagrange, CNRS, Blvd de l’Observatoire, 06304 Nice, France
^{2}
Division of Astrophysics, Department of Physics, Lund University, Box 43 22100 Lund, Sweden
email: niels.nieuwmunster@oca.euemail: lsmith@ast.cam.ac.ukemail: jason.sanders@ucl.ac.uk
^{3}
Department of Physics, University of Surrey, Guildford GU2 7XH, UK
^{4}
Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK
^{5}
MaxPlanck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
^{6}
Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain
^{7}
School of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK
^{8}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{9}
Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
Received:
18
December
2023
Accepted:
21
February
2024
Context. While orbital analysis studies were so far mainly focused on the Galactic halo, it is possible now to do these studies in the heavily obscured region close to the Galactic Centre.
Aims. We aim to do a detailed orbital analysis of stars located in the nuclear stellar disc (NSD) of the Milky Way allowing us to trace the dynamical history of this structure.
Methods. We integrated orbits of the observed stars in a nonaxisymmetric potential. We used a Fourier transform to estimate the orbital frequencies. We compared two orbital classifications, one made by eye and the other with an algorithm, in order to identify the main orbital families. We also compared the Lyapunov and the frequency drift techniques to estimate the chaoticity of the orbits.
Results. We identified several orbital families as chaotic, ztube, xtube, banana, fish, saucer, pretzel, 5:4, and 5:6 orbits. As expected for stars located in a NSD, the large majority of orbits are identified as ztubes (or as a subfamily of ztubes). Since the latter are parented by x_{2} orbits, this result supports the contribution of the bar (in which x_{2} orbits are dominant in the inner region) in the formation of the NSD. Moreover, most of the chaotic orbits are found to be contaminants from the bar or bulge which would confirm the predicted contamination from the most recent NSD models.
Conclusions. Based on a detailed orbital analysis, we were able to classify orbits into various families, most of which are parented by x_{2}type orbits, which are dominant in the inner part of the bar.
Key words: stars: kinematics and dynamics / Galaxy: bulge / Galaxy: center / Galaxy: nucleus / Galaxy: structure
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The nuclear stellar disc (NSD) is a dense stellar structure in the centre of the Milky Way and surrounds the massive nuclear star cluster (NSC) with its central massive black hole (Launhardt et al. 2002). The NSD is a flattened disc with a radius of ∼200 pc and a scale height of ∼50 pc (Launhardt et al. 2002; Nishiyama et al. 2013, GallegoCano et al. 2020). Increasing evidence is reported that the NSD is a distinct structure from the NSC and the nuclear bulge: NoguerasLara et al. (2020) determined the SFH in the NSD using the GALACTICNUCLEUS data (NoguerasLara et al. 2018) and analysing the luminosity function together with stellar evolutionary models. They found that ∼80% of the stars formed more than 8 Gyr ago, followed by a quenching phase and then by a recent star formation activity, about 1 Gyr ago, in which about 5% of the NSD mass was formed. While most of studies agree that the NSD has a relatively early formation time, the detailed SFH is still under discussion (see e.g. NoguerasLara et al. 2023, Sanders et al. 2024), and much more work is clearly needed.
By using a large sample of KMOS observations in the NSD (Fritz et al. 2021), Schultheis et al. (2021) found a difference in the chemistry, that is, in the metallicity distribution function, between the NSC, NSD, and the nuclear bulge that reinforces a different formation scenario of the NSD. Furthermore, they found some evidence that metalrich stars may have formed in the central molecular zone, while metalpoor stars show more similarities to the surrounding Galactic bulge.
Kinematic studies relying on radial velocity measurements or proper motion studies show evidence that the NSD is rotating (see e.g. Lindqvist et al. 1992, Schönrich et al. 2015, Fritz et al. 2021, Shahzamanian et al. 2022). Linking the rotation to the chemistry, Schultheis et al. (2021) found that metalrich stars rotate faster than metalpoor stars, with some hints of counterrotation for the most metalpoor stars.
Extragalactic studies showed that many barred galaxies host nuclear discs or rings (Gadotti et al. 2019, 2020). So far, the most likely formation scenario of a nuclear disc, called insideout formation, is linked to the galactic bar (Bittner et al. 2020). According to this scenario, a nuclear disc is a built up from a series of gaseous rings (i.e. nuclear rings) that grow in radius over time. The growth is caused by the gas that is moved towards the galactic centre by the bar.
Based on the 3D velocities, Sormani et al. (2022) constructed axisymmetric selfconsistent equilibrium dynamical models of the NSD providing the full 6D distribution function (position and velocity) of the NSD. These models provide the best description of the rotation curve in the innermost few hundred parsecs of the Milky Way, and they are implemented in the AGAMA (Vasiliev 2019) software package.
A powerful method for obtaining a complete picture of the properties of the individual orbits is the socalled frequency analysis (Laskar 1993, Valluri & Merritt 1998, Vasiliev 2013) where the three fundamental frequencies of the orbit oscillation can be extracted accurately. This frequency analysis can be used to distinguish between regular and chaotic orbits and to classify the main orbital families. So far, this technique has mainly been used for studies in the Galactic halo and disc. Amarante et al. (2020) used frequency maps to show that the proposed wedges in the R_{apo} − z_{max} plane identified by Haywood et al. (2018) as possible signs of accretion can be explained by the existence of different orbital families. Koppelman et al. (2021) revealed the prominent presence of resonances. According to this, ∼30% of the halo stars are associated with resonant families.
In this paper, we calculate the orbital parameters of a representative sample of stars belonging to the NSD and we use frequency analysis to classify the different types of orbits present in the NSD.
2. Observations
2.1. Data sample
We used the NSD data obtained with the KMOS (Sharples et al. 2013) spectrometer at the ESO VLT. The detailed survey strategy and data reduction procedures are described by Fritz et al. (2021). As well as their radial velocities, we used the proper motions derived from the preliminary VIRAC2 (Smith et al. 2018) photometric and astrometric reduction of the VVV data (Minniti et al. 2010). We refer to Sormani et al. (2022) for a more thorough description of the data. Our sample includes only stars that are primary sources of the survey leaving us a total of 2501 stars.
2.2. Catalogue selection
Sormani et al. (2022) pointed out that stars belonging to the Galactic bar contribute significantly in the outermost fields of the NSD sample of Fritz et al. (2021). In their Table 10, they quantified the contamination of the bar: It ranges from 20–30% in the inner fields to up to 75% in the outermost fields. For this reason, we decided to use only the innermost fields, that is, l< 1.5^{o} and b< 0.3^{o}, where the probability of NSD membership is higher than 70% (see Table 2 of Sormani et al. 2022). In addition, we used the same colour cut (H − Ks) > max(1.3, −0.0233 K + 1.63) as in Schultheis et al. (2021) to remove foreground stars.
In order to obtain highly reliable orbital parameters, we constructed a golden sample and chose to keep stars with small proper motion uncertainties: μ_{l,err} < 0.6 mas yr^{−1} and μ_{b,err} < 0.6 mas yr^{−1} (corresponding to the 98th percentile; see the vertical dashed line in Fig. 1), leading to a total sample of 1130 stars. As pointed out by Sormani et al. (2022), this proper motion cut also removes very bright stars (K < 10), for which the proper motion errors become large due to saturation effects. Figure 2 shows the corresponding HK colourmagnitude diagram after removing the bright stars. The final sample contains stars in the interval 6.6575 < K − 1.37 (H − K) < 9.1575, corresponding to a dereddened magnitude of 7.0 < K_{0} < 9.5. The main purpose of our selection criteria was to obtain an unbiased sample in the metallicity distribution as discussed in Schultheis et al. (2021). We caution that our sample is far from complete, which affects the fraction of the different orbital families (see Sect. 5). To overcome this, we plan to address this issue in a forthcoming paper by using Nbody simulations of the NSD.
Fig. 1. Upper panel: histograms of the radial velocities and proper motions in μ_{l} and μ_{b}, respectively. Lower panel: same, but for the uncertainties. The grey sample is the full sample, and our final sample, used for the analysis, is shown in magenta. 
Fig. 2. K vs. H–K colour magnitude diagram of the total Fritz et al. (2021) sample in grey. Our sample after application of the proper motion cuts is shown in magenta. Colour cut 1: (H − K) > max(1.3, −0.0233K + 1.63). Colour cuts 2 and 3: 6.6575 < K − 1.37(H − K) < 9.1575. The typical error for each axis is indicated in the lower left corner of the figure. 
3. Analysis
3.1. Orbit integration
We used the software package AGAMA (Vasiliev 2019) to determine the orbital parameters. It is both fast in terms of computation time and provides methods for the easy handling of the computation of different potentials. This gives us the flexibility to test different potentials (e.g. axisymmetric, nonaxisymmetric, or different bar pattern speeds).
One main source of uncertainties when calculating the orbital parameters is the distance uncertainty. While previous studies of the GC used a constant distance (e.g. 8.25 kpc, GRAVITY Collaboration 2020), we implemented a distance distribution in our analysis. As the NSD is an extended stellar feature with a scale length of ∼100 pc (e.g. GallegoCano et al. 2020, Sormani et al. 2022, NoguerasLara 2022), we assumed this typical scale length and integrated the orbits 100 times by using different distances chosen from a Gaussian distribution (μ = 8.2 kpc, σ = 50 pc) between 8.1 kpc and 8.3 kpc (Launhardt et al. 2002, NoguerasLara 2022).
In addition, we carried out two other tests in which we assigned the distances of the stars in two different ways: (i) We used the colour cuts made by NoguerasLara et al. (2023) to identify, the stars in our sample that belong to the close edge (i.e. bluer in H–K) and those that belong to the inner region of the NSD (i.e. redder in H–K). Typical distances of 8.05 kpc and 8.2 kpc with a standard deviation of 50 pc were used for the closest edge and the inner population, respectively. (ii) We computed the most likely distances for each star, by using their probability to belonging to the NSD. The probabilities were derived from the distribution function provided by the selfconsistent model from Sormani et al. (2022). The results obtained with these different distance values are discussed in Sects. 4 and 5.
For the purpose of this paper, we constructed a nonaxisymmetric potential by combining the potentials that correspond to the main components of the Galaxy that affect the dynamics of stars in the NSD: the inner bulge/bar, NSD, and NSC. The total density is ρ_{tot} = ρ_{bar} + ρ_{NSD} + ρ_{NSC} where (i) the bar/bulge density is taken from Launhardt et al. (2002); (ii) We adopted the NSD density from Sormani et al. (2020) (see their Eq. (27)); (iii) We used the NSC density from Chatzopoulos et al. (2015) (see their Eq. (17)). We show a comparison with other potentials in Appendix B.
We took into account a typical (clockwise) rotation pattern speed of Ω_{b} = 40 km s^{−1} kpc^{−1} (Portail et al. 2017) for the bar component of the potential which is at an angle of α = 25° from the line of sight towards the Galactic Centre (GC).
For the purpose of this study, it is not necessary to set a long integration time because stars located near the GC rotate quickly around it. As explained in Valluri et al. (2016), an accurate determination of the fundamental orbital frequencies requires that orbits are integrated for at least 20 orbital periods. During 500 Myr, most of the stars of our sample have therefore made hundreds to thousands of periods, which is enough to estimate orbital frequencies (see Sect. 3.2). We chose a short timestep of 4000 yr for a good sampling of the orbits.
3.2. Orbital frequency determination
Bounded regular orbits in a triaxial potential have three fundamental frequencies, Ω, that determine the periodic behaviour of motion (Binney & Tremaine 2008). This motion can be decomposed as a Fourier sum, where the Fourier frequencies are linear combinations of the fundamental frequencies. These fundamental frequencies can be recovered, as shown by Laskar (1993), using the socalled numerical approximation of fundamental frequencies (NAFF), which has been applied in planetary dynamics. PriceWhelan (2015) developed an opensource code called SuperFreq, which is a Python implementation similar to the NAFF code. This code finds the fundamental frequencies (starting with the highest amplitude) by computing the Fourier spectra for the phasespace coordinates used to describe the orbit. A plot of these fundamental frequencies of a set of orbits gives then a frequency map that allows identifying resonances. Resonant orbits are those for which the fundamental frequencies Ω = Ω_{x}, Ω_{y}, Ω_{z} (or Ω_{R}, Ω_{Φ}, Ω_{z} depending on the selected coordinate system) can be measured and where n ⋅ Ω = 0 for n = (n_{x}, n_{y}, n_{z}), where the vector n only contains integer numbers and at most one zero (Koppelman et al. 2021). For example, in Cartesian coordinates, n_{x}Ω_{x} + n_{y}Ω_{y} + n_{z}Ω_{z} = 0, and the resonance is then called n_{x} : n_{y} : n_{z}. Moreover, boxlets are special cases of resonant orbits in which one of the integers n_{x}, n_{y}, n_{z} is zero.
Frequency maps are a powerful tool for obtaining an automatic classification of the different orbital families with a much clearer separation when plotting the fundamental frequencies in Cartesian coordinates (Valluri et al. 2016).
We can classify orbits into two main categories: tubes (shortaxis or longaxis)^{1}, and box orbits (chaotic orbits correspond to orbits that cannot be identified with one of these categories, see Sect. 3.3). Based on the frequencies in Cartesian coordinates (Ω_{x}, Ω_{y}, Ω_{z}), we can identify most of the orbital families corresponding to the different resonances. However, since the chosen coordinate system only allows us to trace symmetries within the system, shortaxis and longaxis tube orbits in the Cartesian case are clustered to a trivial resonance (e.g. 1:−1:0 resonance for shortaxis tube orbits). It is then more effective to study orbital frequencies in several coordinate systems even though the Cartesian frequencies give us a general picture of the orbital families. To study tube orbits more precisely and identify their true resonances, it can be better to examine the frequencies in cylindrical coordinates (Ω_{R}, Ω_{Φ}, Ω_{z}). However, Valluri et al. (2016) showed that Cartesian orbital frequencies are more reliable for bar orbit classification.
3.3. Chaoticity
To study chaoticity and particularly identify chaotic orbits, we used two different methods that we introduce in this section.
3.3.1. Lyapunov exponent
The Lyapunov exponent (Lyapunov 1992) is a fundamental concept in the field of nonlinear dynamics and chaos theory. It quantifies the sensitivity of a dynamical system to initial conditions. In a global context, the Lyapunov exponent characterizes how the trajectories in a system diverge or converge as time progresses, providing insights into the longterm behaviour of complex systems.
When considering a dynamical system, even tiny differences in initial conditions can lead to drastically different trajectories over time. The Lyapunov exponent captures this phenomenon by measuring the rate of exponential separation between initially close trajectories. A high Lyapunov exponent indicates chaotic behaviour, where trajectories diverge exponentially over time, making longterm predictions inherently uncertain. Conversely, a low Lyapunov exponent indicates convergence towards a stable equilibrium or periodic behaviour.
We used the highest Lyapunov exponent estimation method provided by AGAMA (Vasiliev 2019). For more details about the method, see Sect. 4.3 of Vasiliev (2013).
3.3.2. Frequency drift
Another way to study chaoticity is to directly use the orbital frequencies. Valluri et al. (2010) (see their Sect. 3.1) showed that it is possible to measure the stochasticity of an orbit based on the change in the fundamental frequencies over two consecutive time intervals. For each frequency component f_{i}, they computed what they called the frequency drift:
where i defines the frequency component in Cartesian coordinates (i.e. log(Δf_{x}), log(Δf_{y}) and log(Δf_{z})). The highest value of the three frequency drift parameters log(Δf_{i}) is then associated with the frequency drift parameter log(Δf). The higher the value of log(Δf), the more chaotic the orbit. However, as shown by Valluri et al. (2010), the accuracy of the frequency analysis requires at least 20 oscillation periods in order to avoid a misclassification of the orbits as chaotic. Figure 3 shows the distribution of the frequency drift parameter log(Δf). In order to define a threshold value of log(Δf) at which orbits are classified as chaotic, we followed a similar approach as in Valluri et al. (2010). 1.5 % of the orbits have log(Δf) > − 0.2 which we consider as the threshold of being chaotic orbits. In total, we obtained 18 chaotic orbits in our sample. Figure 4 shows the comparison between the frequency drift parameter and the highest Lyapunov exponent. These two measurements of the chaoticity are generally correlated but the dispersion is high. Our visual classification of chaotic orbits (see Sect. 4.2.2, red points in Fig. 4) nicely shows that chaotic orbits can indeed be identified based on their large Lyapunov exponent and high frequency drift parameter.
Fig. 3. Histogram of the frequency drift. The dashed red area denotes the 98.5 percentile limit above which all orbits are classified as chaotic orbits. 
Fig. 4. Comparison between highest Lyapunov exponent and frequency drift. The red points show the chaotic orbits identified with the visual classification (see Sect. 4.2.2) while the red shaded area indicates the zone containing only chaotic orbits according to the frequency drift method (see Sect. 3.3.2). Only orbits with a highest Lyapunov exponent greater than 0 are shown in this diagram. 
4. Results
4.1. R_{max} vs. z_{max} diagram
Figure 5 shows the apocentric radius (R_{max}) versus the maximum height (z_{max}) of our NSD sample colourcoded by the eccentricity.
Fig. 5. R_{max} vs. z_{max} diagram colourcoded with eccentricity. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 
This diagram shows several features: (i) One striking feature is that the stars are not homogeneously distributed, but congregate in distinct diagonal wedges in which z_{max} increases with R_{max}. (ii) Highly eccentric orbits are confined to one wedge. (iii) Some of these highly eccentric orbits extend beyond the typical radius of the NSD and are likely stars related to the Galactic bar. Figure 6 shows the impact for a relaxed assumption of a constant distance and when the distances of the stars are instead allowed to vary within a spread of 100 pc for the NSD, using MCMC simulations (as detailed in Sect. 3). The key characteristics mentioned earlier clearly persist. This indicates that the uncertainty in the distances of NSD stars can safely be disregarded.
Fig. 6. Average R_{max} vs. z_{max} diagram (colourcoded with eccentricity) by using 100 MCMC distances. The standard deviation for each star is indicated with the error bars. 
Figure 7 shows the same feature, but as a function of rotational velocity V_{ϕ}. The NSD clearly rotates with typical velocities of 80 km s^{−1}, which agrees with the works of Lindqvist et al. (1992), Schönrich et al. (2015), Schultheis et al. (2021), Shahzamanian et al. (2022), and Sormani et al. (2022). The figure also shows stars with slower rotation and even counterrotating stars, which have been identified previously by Schultheis et al. (2021).
Fig. 7. R_{max} vs. z_{max} diagram colourcoded with the rotational velocity. V_{ϕ} > 0 means clockwise (i.e. the same direction as the bar). The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 
4.2. Orbit classification
We discuss and compare the two methods we used to classify orbits in various orbital families.
4.2.1. Automatic method
We call this method “automatic” because no visual inspection of the orbits is used to classify them (compared to the other classification method presented in Sect. 4.2.2).
Firstly, to differentiate tube orbits from chaotic/box orbits, we determined whether the orbit circulates around a specific axis (e.g. a circulation around the xaxis results in an xtube) or not (box orbit or chaotic orbit). To do this, we verified whether there was a change in the sign of the angular momentum about an axis. We find 34% box orbits and 66% tube orbits among our 1130 stars.
Secondly, to identify the different orbital families, we studied the resonances on the frequency map as explained in Sect. 3.2. In the Cartesian version of the frequency map, Fig. 8, we observe two main structures: (i) a diagonal line called the 1:−1:0 resonance, where shortaxis tubes (called here after “ztubes”) lie, and (ii) a cloud located above the latter resonance and below the 0:1:−1 resonance that is thought to contain mainly box and chaotic orbits. In addition to these two large structures that clearly dominate our frequency map, a few orbits lie close to the 0:1:−1 resonance where longaxis tubes, called here after xtubes, are located. An example of the morphology of an xtube is given in Fig. A.2.
Fig. 8. Frequency map in Cartesian coordinates vs. highest Lyapunov exponent. For reasons of legibility, the banana, fish, pretzel, 5:4 and 5:6 resonances are only shown here for the (x, z) plane case. The typical error on both frequency ratios arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 
Based on the three main types of orbits (box/chaotic, ztube, xtube), we proceeded with the classification in order to distinguish the resonant orbits that populate the two large structures of the map.
Classical tube orbits are simply orbits with a resonance 1:1 between two of the three frequency ratios. In the same way, we obtained banana orbits with a 2:1 resonance (see Fig. A.3), fish orbits with 3:2 (see Fig. A.4), pretzel orbits with 4:3 (see Fig. A.6), and so on. By plotting lines that correspond to resonances on the frequency map, we gained a first idea of the orbital families that populate our sample. For the same global resonance, for instance, 4:3, this kind of orbit can be observed in different planes (x, y), (x, z), or (y, z), that match the resonances 4:−3:0, 4:0:−3, and 0:4:−3, respectively. In addition to the previously introduced families, we also searched for other pretzel varieties: 5:4 orbits (see Fig. A.7) and 5:6 orbits (see Fig. A.8).
As explained Sect. 3.2, we also considered the frequencies in cylindrical coordinates (see Fig. 9), which were only computed for ztubes (because our sample contains very few xtubes, we did not consider it necessary to study the frequency map in their corresponding cylindrical coordinates). There is no bisymmetry around Ω_{ϕ}/Ω_{R} = 0, which is another way of detecting the rotation of the NSD. More stars orbit in the positive than in the negative sense. The resonances also appear as straight lines here (horizontal lines for a resonance between the vertical oscillation frequency Ω_{z} and the radial oscillation frequency Ω_{R} and vertical lines between the azimuthal oscillation frequency Ω_{ϕ} and the radial oscillation frequency Ω_{R}). We discern then some resonances, Ω_{z} : Ω_{R} = 1 : 1, 3 : 4, 1 : 2, and Ω_{ϕ} : Ω_{R} = 1 : 2. Saucer orbits (Sambhus & Sridhar 2000;Vasiliev 2014) (see an example Fig. A.5) were found via the resonance Ω_{z} : Ω_{R} = 1 : 1 (Yavetz et al. 2023), and we therefore have an estimate of their number.
Fig. 9. Frequency map in cylindrical coordinates for ztubes alone (orbits for which a circulation around the zaxis has been detected), i.e. orbits from the resonance (1:−1:0) in Cartesian coordinates (see Fig. 8). The horizontal and vertical lines correspond to resonances between Ω_{z} and Ω_{R} and between Ω_{ϕ} and Ω_{R} respectively. The typical error on both frequency ratios arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 
We call a chaotic orbit a weakly chaotic or sticky orbit when it lies near a resonance (Valluri et al. 2016). It might appear regular for some time, but will finally exhibit chaotic behaviour. With this method, we selected (see Figs. 8 and 9) all orbits close enough to the resonances in the frequency maps. More precisely, we made the selection at ±0.01 (in frequency ratio), which is enough for contamination by sticky chaotic orbits. However, we limited this contamination by removing orbits from the selection of tubes, that were identified as box/chaotic orbits (without any circulation around an axis). We list in Table 1 the percentage for each family with and without these latter contaminants.
Family membership results for the automatic and visual methods
4.2.2. Visual method
To evaluate the validity of the automatic orbit classification, we carried out a visual classification by studying each of our 1130 orbits. We were able to identify orbits belonging to the families presented in Sect. 4.2.1 that are chaotic, xtube, ztube, saucer, banana, fish, pretzel, 5 : 4 and 5 : 6 orbits. In addition to this, unlike the automatic method, the visual identification allowed us to differentiate centrophobic and centrophilic orbits of the same family (i.e. of the same resonance). Therefore, we have an estimate of the quantity of antibanana,fish, and pretzel orbits. We are also able to identify orbits precisely whose own wellknown resonance is lacking in the Cartesian/cylindrical frequency maps as for example saucer orbits. Because it was too difficult to visually identify box from chaotic orbits, we decided to class them in the same category for the rest of the orbital analysis as in Sect. 4.2.1.
As explained in Sect. 4.2.1, resonant orbits have specific ratios of frequencies that correspond to their resonance. This information is clearly visible in the orbital shape. Therefore, we can visually determine the resonance of an orbit by counting how many times the orbit crossed each of the axes and obtain ratio from this. We identified the different orbital families by using this method.
The results of this classification show in the frequency maps of the Figs. A.1–A.8 that the stars belonging to the identified orbital families lie near the corresponding resonances with a larger scattering than the maximum distance from the resonance we took for the automatic method in Sect. 4.2.1. However, most of the stars for each family are still located very close to the resonance. The identified chaotic/box orbits populate the frequency map in a less regular manner, as expected.
The location of the stars belonging to the identified orbital families in the R_{max} vs. z_{max} diagram shows that they are arranged in preferential zones:
(i) Chaotic orbits that are scattered all over the frequency map, are only located along the observed filament in the R_{max} vs. z_{max} diagram (see Fig. A.1).
(ii) xtubes only populate the uppermost part of the R_{max} vs. z_{max} diagram, that is, they have the highest z_{max} values for a wide range of R_{max} (see Fig. A.2).
(iii) Banana (and antibanana) orbits are found to occupy a few wedges, but seem to be mainly present in the upper part of the R_{max} vs. z_{max} diagram (see Fig. A.3).
(iv) Fish (and antifish) orbits are scattered and lack a strong tendency for any wedge (see Fig. A.4).
(v) Saucer orbits populate the middle part of the diagram and only appear to correspond to these R_{max}, z_{max} values (see Fig. A.5).
(vi) Pretzel (and antipretzel) orbits and their derivatives (5:4 and 5:6 orbits) are mainly located in the upper wedges (see Figs. A.6–A.8 respectively).
5. Discussion
We made the first orbital analysis of stars located in the NSD and also using a nonaxisymmetric potential. The purpose of this study therefore is to obtain a general picture of the dynamical signatures in the NSD by studying the orbital resonances.
Wedges in the R_{max} vs. z_{max} diagram have been detected in the Galactic halo by Haywood et al. (2018), who attributed these substructures to some heating process related to the early phase of the Galactic disc. Koppelman et al. (2021) investigated these features in the halo in detail using orbital frequencies and axisymmetric potentials. They demonstrated that these structures are due to resonant families and that the depletion around these resonances is related to nonintegrable potentials with some indication of chaotic orbits (PriceWhelan et al. 2016).
As explained Sect. 3.1, we integrated orbits using four different distance estimations: a case with a fixed distance value, another case with distances chosen from a Gaussian distribution, one case with two different distances depending on the estimated position of stars along the line of sight, that is, in front or behind, and a final case where the distance of each star was derived from its position probability along the line of sight. For all the cases, we found the similar R_{max} vs. z_{max} diagrams and frequency maps, showing a very similar orbit distribution. These wedges therefore do not depend on the assumed distance inputs. This confirms that the observed substructures are indeed real and that our frequency maps are reliable. Therefore, we clearly detect wedges in the NSD, as shown in Fig. 5, where the resonances are visible as thin straight lines. In addition, we conducted many tests with different Galactic potentials, for instance, assuming a classical Milky Way potential with and without a rotating bar (nonaxisymmetric and axisymmetric; see Fig. B.1) with which the resonances still occur. This confirms that the conclusions of Koppelman et al. (2021) are relevant in this case as well.
The R_{max} vs. z_{max} diagrams Figs. 10 and A.1 show a “filament” at high z_{max} that dives into the structure at low R_{max}. It is mainly composed of chaotic orbits, regardless of the method used to identify them (highest Lyapunov exponent Sect. 3.3.1, frequency drift Sect. 3.3.2 and visual classification Sect. 4.2.2). Because this structure is only visible when the bar potential is used in the combination of potentials and because it becomes denser with increasing integration time, we can conclude that the bar disturbs orbits to the point of making them chaotic. With an axisymmetric potential (i.e. without a bar), this filament disappears. This supports the argument of a clear bar signature of this filament. In addition, a projection of these stars in the (l,b) plane shows that these stars are mostly located in the outer parts of the NSD (l> 1.0°), where the contamination of the bulge/bar stars in the NSD increases significantly (see Fig. 10 of Sormani et al. 2022). Our observed filament consists of ∼20% stars of our sample, which agrees with the predicted contamination rate from the NSD models of ∼25% (Sormani et al. 2022).
Fig. 10. Upper panel: R_{max} vs. z_{max} vs. Lyapunov. Lower panel: R_{max} vs. z_{max} vs. frequency drift. Red indicates the most chaotic orbits in both methods. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner of each plot. 
The distribution of orbits in the Cartesian frequency map Fig. 8, shows a large majority of ztubes that is expected for stars forming a disc structure such as the NSD and a minority of xtubes. In addition, we were able to identify a variety of families that mostly belong to ztubes. We introduced in Sect. 4 two orbital classification methods used in this study: an automatic (see Sect. 4.2.1) and visual method (see Sect. 4.2.2). Based on previous works about orbit classification (Valluri et al. 2010, 2012, 2016), we were able to search several orbital families. Therefore, the automatic and visual classifications contain the same families and we compare their number in Table 1. The results of the two methods agree very well for the different identified families. However, the automatic way seems to slightly overestimate chaotic, saucer, 5 : 4, and 5 : 6 orbits and underestimate ztube, xtube, banana, fish and pretzel orbits. Even though the visual method can lead to false positives, it still remains better than the automatic one because of the visual check that allowed us to identify every orbital family without any dependence on resonance knowledge. In the era of machine learning (ML) techniques, the visual classification might in future be performed using different algorithms. de los Rios et al. (2021) tested different ML techniques such as a random forest technique, which is a supervised decision tree algorithm, a support vector machine, or Knearest neighbours to classify galaxies in, and around clusters according to their projected phasespace position. In order to obtain a precise classification, a large training set is necessary on which Nbody simulations can be used.
Figure 11 shows that fish and saucer orbits dominate the middle and low part of the R_{max} vs. z_{max} diagram, unlike banana and pretzel (and 5:4, 5:6) orbits, which cover the high part. In addition, the few xtubes are located at the highest z_{max}. Therefore, the results given by the visual classification confirm the link between the observed wedges in the R_{max} vs. z_{max} diagram and orbital resonances, and more precisely, orbital families.
Fig. 11. R_{max} vs. z_{max} diagram with the different identified orbital families. Separate diagrams for each family are given in Appendix A. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 
One direct application of the orbital parameter is the determination of the radial and vertical extent of the NSD. Figure 12 shows the median apocentric radius R_{max} (left panel) and z_{max} for stars with a different minimum eccentricity threshold (ecc_{min}). For ecc_{min} < 0.25, the truncation radius of the NSD (indicated by the grey area) is 138 ± 5 pc and the vertical extent is 53 ± 1 pc. Schönrich et al. (2015) constructed a toy model based on APOGEE kinematic data assuming a disk in which the best fit suggest a truncation radius of 150 pc (see their Table 1). This is very consistent with the values found in this work but smaller than the radial scale of ∼230 pc found by Launhardt et al. (2002). This could be due to the heavy differential reddening as they used photometric data. The vertical extent is about 50 pc in Schönrich et al. (2015), but they did not confine this parameter from their toy model.
Fig. 12. Median apocentric radius R_{max} and maximum height z_{max} for stars with a different minimum eccentricity threshold. For instance, the median R_{max} for a minimum eccentricity of 0.5 means that we only computed the median of the R_{max} values of stars with an eccentricity above 0.5. 
We stress that our sample size is limited due to our selection cuts (see Sect. 2.2), leading to statistically small sample sizes for each derived family of orbits. A first rough estimate of the completeness of our sample by using luminosity functions (see e.g. NoguerasLara et al. 2023) shows that our working sample is indeed highly incomplete (∼15% completeness at K_{S} ∼ 11.5 mag).
Complementary highresolution Nbody simulations of the NSD are necessary to (i) extend our study to a much larger sample, (ii) compare simulations and our observed data set in detail, and (iii) study the effect of different Galactic potentials on the derived orbital families.
In the vicinity of the Galactic bar, orbits perpendicular to the bar are called x_{2}type, and orbits parallel to the bar are called x_{1}type orbits. x_{2}type orbits are parents of the ztube family, occupying the 1:−1:0 main resonance line in the frequency map (see Fig. 8), where about twothirds of our sample lies. The formation of nuclear stellar discs is strongly connected to the properties of the galactic bar. Hydrodynamical simulations suggest that nuclear stellar discs form close to the inner Lindblad resonance of the main bar, where the x_{2}type orbital family dominates (Athanassoula 1992a,b, Li et al. 2015, Sormani et al. 2018, 2024). After the Galactic bar was formed, gas was funnelled along the bar towards the GC, where it settled down to form a nuclear disc. The gas then formed stars, which keep their resulting x_{2} orbits, and the resulting stellar population resembles a disc. Additional evidence comes from the spatial and kinematical overlap between the central molecular zone and the NSD (Schönrich et al. 2015, Schultheis et al. 2021) in the Milky Way, which supports this scenario that the stars in the NSD should be more metalrich and dynamically cooler than the surrounded bar/bulge populations. This feature can also be seen in external galaxies with NSDs (Bittner et al. 2020, Gadotti et al. 2019).
Planar periodic orbits x_{1}, x_{2}, and x_{4} share the 1 : 2 resonance between the tangential oscillation frequency (Ω_{ϕ}) and the radial frequency (Ω_{R}). Unfortunately, they are not anticipated to be abundant in Nbody simulations (Valluri et al. 2016) and it is a delicate task to identify these orbits visually. This prevented us from distinguishing them properly. However, we were able to detect bifurcations from the x_{1}tree (Athanassoula 2005) with banana orbits that were called “x_{1}v_{1}” by Skokos et al. (2002). The presence of the latter would suggest the existence of an inner bar embedded in the NSD, as observed in others galaxies (MéndezAbreu et al. 2019), but a more indepth study is needed.
6. Conclusions
We presented a detailed orbital analysis of stars in the NSD by using orbital frequencies and a visual classification of the orbits obtained from the orbital integration. A comparison between these two methods shows very similar results in the classification of the different orbital families, such as chaotic/box, ztube, xtube, banana, fish, saucer, pretzel, 5:4, and 5:6 resonances. The large majority of sources are ztubes, which is expected for stars located in an NSD. We used two different methods: the Lyapunov exponent, and the frequency drift. We estimated the chaoticity of the orbits with these methods and showed by using in addition a visual classification that chaotic orbits can be identified by their high Lyapunov exponent as well as by their high frequency drift parameter. They occupy in the R_{max} vs. z_{max} the filament at high z_{max} where the bar most likely causes the highly chaotic orbits. About 20% of our stars are contaminated by the bar/bulge population which is very close to the predicted 25% contamination from the most recent NSD models of Sormani et al. (2022). We emphasize that we performed here a statistical approach and that the orbits of the individual stars can be affected by the presence of molecular clouds as well as by the disruption of their trajectory by tidal forces (e.g. Portegies Zwart et al. 2002; Kruijssen et al. 2014).
We detected clear substructures in the R_{max} vs. z_{max} diagram in the NSD that were identified as wedges that are related to different resonances and therefore different orbital families. 68.2% of our sample show ztube orbits that are parented by the x2type orbits. This is indeed expected if the formation of the NSD is coupled with the formation of the Galactic bar where x2type orbits are the dominant population. As a followup work, a comparison with selfconsistent models of the NSD (e.g. Sormani et al. 2022) is necessary for a detailed comparison between the observations and the predictions from the models (e.g. fraction of different orbital families) and to improve our orbital classification.
We used the R_{max} vs. z_{max} diagram to constrain the radial and vertical extent of the NSD with 138 ± 5 pc and 53 ± 1 pc, respectively, which is consistent with the values found in Schönrich et al. (2015). Due to our limited sample size, no obvious trends between the orbital families and the chemistry (e.g. metallicities) is detected. Future large surveys of the NSD, such as the upcoming MOONS survey (Gonzalez et al. 2020) will clearly benfit from the large number of stars to establish the strong connection between dynamics and chemistry.
Acknowledgments
We want to thank M. Valluri, A. PriceWhelan, E. Vasiliev for fruitful discussions. R. Schödel acknowledges financial support from the Severo Ochoa grant CEX2021001131S funded by MCIN/AEI/ 10.13039/501100011033, from grant EUR2022134031 funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR, and from grant PID2022136640NBC21 funded by MCIU/AEI 10.13039/501100011033 and by the European Union. MCS acknowledges financial support of the Royal Society (URF\R1\221118). N.N. thanks Elisa Denis and Gabriele Contursi for useful comments. Softwares: AGAMA (Vasiliev 2019), SuperFreq (PriceWhelan 2015), Gala (PriceWhelan et al. 2022).
References
 Amarante, J. A. S., Smith, M. C., & Boeche, C. 2020, MNRAS, 492, 3816 [Google Scholar]
 Athanassoula, E. 1992a, MNRAS, 259, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E. 1992b, MNRAS, 259, 345 [Google Scholar]
 Athanassoula, E. 2005, Ann. New York Acad. Sci., 1045, 168 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: SecondEdition, (Princeton, NJ USA: Princeton University Press) [Google Scholar]
 Bittner, A., SánchezBlázquez, P., Gadotti, D. A., et al. 2020, A&A, 643, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948 [Google Scholar]
 de los Rios, M., Martínez, H. J., Coenda, V., et al. 2021, MNRAS, 500, 1784 [Google Scholar]
 Fritz, T. K., Patrick, L., FeldemierKrause, A., et al. 2021, A&A, 649, A83 [Google Scholar]
 Gadotti, D. A., SánchezBlázquez, P., FalcónBarroso, J., et al. 2019, MNRAS, 482, 506 [Google Scholar]
 Gadotti, D. A., Bittner, A., FalcónBarroso, J., et al. 2020, A&A, 643, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GallegoCano, E., Schödel, R., NoguerasLara, F., et al. 2020, A&A, 634, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gonzalez, O. A., Mucciarelli, A., Origlia, L., et al. 2020, The Messenger, 180, 18 [NASA ADS] [Google Scholar]
 GRAVITY Collaboration (Abuter, R., et al.) 2020, A&A, 636, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113 [Google Scholar]
 Koppelman, H. H., Hagen, J. H. J., & Helmi, A. 2021, A&A, 647, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G., et al. 2014, MNRAS, 440, 3370 [Google Scholar]
 Laskar, J. 1993, Celestial Mech. Dyn. Astron., 56, 191 [Google Scholar]
 Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, Z., Shen, J., & Kim, W.T. 2015, ApJ, 806, 150 [Google Scholar]
 Lindqvist, M., Habing, H. J., & Winnberg, A. 1992, A&A, 259, 118 [NASA ADS] [Google Scholar]
 Lyapunov, A. M. 1992, Int. J. Control, 55, 531 [Google Scholar]
 MéndezAbreu, J., de LorenzoCáceres, A., Gadotti, D. A., et al. 2019, MNRAS, 482, L118 [Google Scholar]
 Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433 [Google Scholar]
 Nishiyama, S., Yasui, K., Nagata, T., et al. 2013, ApJ, 769, L28 [Google Scholar]
 NoguerasLara, F. 2022, A&A, 668, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 NoguerasLara, F., GallegoCalvente, A. T., Dong, H., et al. 2018, A&A, 610, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 NoguerasLara, F., Schödel, R., GallegoCalvente, A. T., et al. 2020, Nat. Astron., 4, 377 [Google Scholar]
 NoguerasLara, F., Schultheis, M., Najarro, F., et al. 2023, A&A, 671, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 2002, ApJ, 565, 265 [NASA ADS] [CrossRef] [Google Scholar]
 PriceWhelan, A. M. 2015, https://doi.org/10.5281/zenodo.18787 [Google Scholar]
 PriceWhelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079 [NASA ADS] [CrossRef] [Google Scholar]
 PriceWhelan, A., Sipőcz, B., Starkman, N., et al. 2022, https://doi.org/10.5281/zenodo.6325733 [Google Scholar]
 Sambhus, N., & Sridhar, S. 2000, ApJ, 542, 143 [Google Scholar]
 Sanders, J. L., Smith, L., Evans, N. W., & Lucas, P. 2019, MNRAS, 487, 5188 [CrossRef] [Google Scholar]
 Sanders, J. L., Kawata, D., Matsunaga, N., et al. 2024, MNRAS, 530, 2972 [Google Scholar]
 Schönrich, R., Aumer, M., & Sale, S. E. 2015, ApJ, 812, L21 [Google Scholar]
 Schultheis, M., Fritz, T. K., Nandakumar, G., et al. 2021, A&A, 650, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shahzamanian, B., Schödel, R., NoguerasLara, F., et al. 2022, A&A, 662, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sharples, R., Bender, R., Agudo Berbel, A., et al. 2013, The Messenger, 151, 21 [NASA ADS] [Google Scholar]
 Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847 [Google Scholar]
 Smith, L. C., Lucas, P. W., Kurtev, R., et al. 2018, MNRAS, 474, 1826 [Google Scholar]
 Sormani, M. C., Treß, R. G., Ridley, M., et al. 2018, MNRAS, 475, 2383 [NASA ADS] [CrossRef] [Google Scholar]
 Sormani, M. C., Magorrian, J., NoguerasLara, F., et al. 2020, MNRAS, 499, 7 [Google Scholar]
 Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022, MNRAS, 512, 1857 [CrossRef] [Google Scholar]
 Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, MNRAS, 528, 5762 [Google Scholar]
 Valluri, M., & Merritt, D. 1998, ApJ, 506, 686 [Google Scholar]
 Valluri, M., Debattista, V. P., Quinn, T., & Moore, B. 2010, MNRAS, 403, 525 [Google Scholar]
 Valluri, M., Debattista, V. P., Quinn, T. R., Roškar, R., & Wadsley, J. 2012, MNRAS, 419, 1951 [Google Scholar]
 Valluri, M., Shen, J., Abbott, C., & Debattista, V. P. 2016, ApJ, 818, 141 [Google Scholar]
 Vasiliev, E. 2013, MNRAS, 434, 3174 [Google Scholar]
 Vasiliev, E. 2014, CQG, 31, 244002 [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [Google Scholar]
 Yavetz, T. D., Johnston, K. V., Pearson, S., PriceWhelan, A. M., & Hamilton, C. 2023, ApJ, 954, 215 [Google Scholar]
Appendix A: Orbital families
We present here typical examples of orbits in the x, y, z plane as well as the location in the R_{max} vs. z_{max} diagram and in the frequency map in Cartesian coordinates. We also show in Fig.A.9 the location of certain families in the cylindrical frequency map.
Fig. A.1. Example of a chaotic orbit. Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the chaotic orbits identified with the visual method (see 4.2.2). 
Fig. A.2. Example of an xtube, i.e. 1 : 1 resonance between Ω_{y} and Ω_{z}, also called (0 : 1 : −1). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the xtube orbits identified with the visual method (see 4.2.2). 
Fig. A.3. Example of a banana orbit (here, a (x, z) banana), i.e. 2 : 1 resonance between two orbital frequencies (here, Ω_{z} and Ω_{x}, also called (2 : 0 : −1)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and antibanana orbits identified with the visual method (see 4.2.2). 
Fig. A.4. Example of a fish orbit (here, a (y, z) fish), i.e. 3 : 2 resonance between two orbital frequencies (here, Ω_{z} and Ω_{y}, also called (0 : 3 : −2)). Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the fish and antifish orbits identified with the visual method (see 4.2.2). 
Fig. A.5. Example of a saucer orbit, i.e. 1 : 1 resonance between Ω_{z} and Ω_{R}. Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the saucer orbits identified with the visual method (see 4.2.2). 
Fig. A.6. Example of a pretzel orbit, i.e. 4 : 3 resonance between two orbital frequencies (here, Ω_{z} and Ω_{y}, also called (0 : 4 : −3)). Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the pretzel and antipretzel orbits identified with the visual method (see 4.2.2). 
Fig. A.7. Example of a 5:4 orbit, i.e. 5 : 4 resonance between two orbital frequencies (here Ω_{z}, and Ω_{y}, also called (0 : 5 : −4)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:4 orbits identified with the visual method (see 4.2.2). 
Fig. A.8. Example of a 5:6 orbit, i.e. 5 : 6 resonance between two orbital frequencies (here, Ω_{y} and Ω_{z}, also called (0 : 6 : −5)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:6 orbits identified with the visual method (see 4.2.2). 
Fig. A.9. Frequency map in cylindrical coordinates for different families (identified with the visual method; see Section.4.2.2) with computable cylindrical frequencies. 
Appendix B: Comparison of potentials
To assess the sensitivity of the study outcomes to the selection of the mean potential, we integrated the orbits using an axisymmetric potential or slightly varied combinations of (nonaxisymmetric) potentials. As detailed in section 3.1, our approach involves considering three components: the bulge/bar, the NSD, and the NSC. In the following figures (Fig. B.1, Fig. B.2, B.3 and B.4), we compare the R_{max} vs. z_{max} diagram presened in this paper (Fig. 5) with those obtained with an axisymmetric potential (MWPotential14 from Bovy 2015) or by changing the NSD potential (Sormani et al. 2022 and Launhardt et al. 2002), the bar potential (Portail et al. 2017), or the bar angle.
Fig. B.1. R_{max} vs. z_{max} diagram comparison for the axisymmetric and nonaxisymmetric cases. Left: Reference diagram, with the nonaxisymmetric combination of potentials presented in Section 3.1. Right: Using the axisymmetric potential MWPotential14 from Bovy (2015). 
Fig. B.2. R_{max} vs. z_{max} diagram comparison with different NSD potentials. Upper panel: Reference diagram, with the NSD potential from Sormani et al. (2020) (model 3) used in our study. Middle panel: Using the NSD potential from Sormani et al. (2022). Lower panel: Using the NSD potential from Launhardt et al. (2002). In the three cases, the bar and the NSC potentials are not changed. 
Fig. B.3. R_{max} vs. z_{max} diagram comparison with different bar potentials. Left: Reference diagram, with the bar potential from Launhardt et al. (2002). Right: Using the bar potential from Portail et al. (2017). In both cases, the NSD and NSC potentials are not changed. 
Fig. B.4. R_{max} vs. z_{max} diagram comparison with a different bar angle with respect to the line of sight towards the GC. Left: Reference diagram as presented in Section 3.1. Right: Case with a 10% variation of the angle. 
All Tables
All Figures
Fig. 1. Upper panel: histograms of the radial velocities and proper motions in μ_{l} and μ_{b}, respectively. Lower panel: same, but for the uncertainties. The grey sample is the full sample, and our final sample, used for the analysis, is shown in magenta. 

In the text 
Fig. 2. K vs. H–K colour magnitude diagram of the total Fritz et al. (2021) sample in grey. Our sample after application of the proper motion cuts is shown in magenta. Colour cut 1: (H − K) > max(1.3, −0.0233K + 1.63). Colour cuts 2 and 3: 6.6575 < K − 1.37(H − K) < 9.1575. The typical error for each axis is indicated in the lower left corner of the figure. 

In the text 
Fig. 3. Histogram of the frequency drift. The dashed red area denotes the 98.5 percentile limit above which all orbits are classified as chaotic orbits. 

In the text 
Fig. 4. Comparison between highest Lyapunov exponent and frequency drift. The red points show the chaotic orbits identified with the visual classification (see Sect. 4.2.2) while the red shaded area indicates the zone containing only chaotic orbits according to the frequency drift method (see Sect. 3.3.2). Only orbits with a highest Lyapunov exponent greater than 0 are shown in this diagram. 

In the text 
Fig. 5. R_{max} vs. z_{max} diagram colourcoded with eccentricity. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 

In the text 
Fig. 6. Average R_{max} vs. z_{max} diagram (colourcoded with eccentricity) by using 100 MCMC distances. The standard deviation for each star is indicated with the error bars. 

In the text 
Fig. 7. R_{max} vs. z_{max} diagram colourcoded with the rotational velocity. V_{ϕ} > 0 means clockwise (i.e. the same direction as the bar). The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 

In the text 
Fig. 8. Frequency map in Cartesian coordinates vs. highest Lyapunov exponent. For reasons of legibility, the banana, fish, pretzel, 5:4 and 5:6 resonances are only shown here for the (x, z) plane case. The typical error on both frequency ratios arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 

In the text 
Fig. 9. Frequency map in cylindrical coordinates for ztubes alone (orbits for which a circulation around the zaxis has been detected), i.e. orbits from the resonance (1:−1:0) in Cartesian coordinates (see Fig. 8). The horizontal and vertical lines correspond to resonances between Ω_{z} and Ω_{R} and between Ω_{ϕ} and Ω_{R} respectively. The typical error on both frequency ratios arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 

In the text 
Fig. 10. Upper panel: R_{max} vs. z_{max} vs. Lyapunov. Lower panel: R_{max} vs. z_{max} vs. frequency drift. Red indicates the most chaotic orbits in both methods. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner of each plot. 

In the text 
Fig. 11. R_{max} vs. z_{max} diagram with the different identified orbital families. Separate diagrams for each family are given in Appendix A. The typical error on both parameters arising from the propagation of the observational uncertainties and the distance uncertainty is indicated in the lower right corner. 

In the text 
Fig. 12. Median apocentric radius R_{max} and maximum height z_{max} for stars with a different minimum eccentricity threshold. For instance, the median R_{max} for a minimum eccentricity of 0.5 means that we only computed the median of the R_{max} values of stars with an eccentricity above 0.5. 

In the text 
Fig. A.1. Example of a chaotic orbit. Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the chaotic orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.2. Example of an xtube, i.e. 1 : 1 resonance between Ω_{y} and Ω_{z}, also called (0 : 1 : −1). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the xtube orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.3. Example of a banana orbit (here, a (x, z) banana), i.e. 2 : 1 resonance between two orbital frequencies (here, Ω_{z} and Ω_{x}, also called (2 : 0 : −1)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and antibanana orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.4. Example of a fish orbit (here, a (y, z) fish), i.e. 3 : 2 resonance between two orbital frequencies (here, Ω_{z} and Ω_{y}, also called (0 : 3 : −2)). Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the fish and antifish orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.5. Example of a saucer orbit, i.e. 1 : 1 resonance between Ω_{z} and Ω_{R}. Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the saucer orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.6. Example of a pretzel orbit, i.e. 4 : 3 resonance between two orbital frequencies (here, Ω_{z} and Ω_{y}, also called (0 : 4 : −3)). Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the pretzel and antipretzel orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.7. Example of a 5:4 orbit, i.e. 5 : 4 resonance between two orbital frequencies (here Ω_{z}, and Ω_{y}, also called (0 : 5 : −4)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:4 orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.8. Example of a 5:6 orbit, i.e. 5 : 6 resonance between two orbital frequencies (here, Ω_{y} and Ω_{z}, also called (0 : 6 : −5)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel:R_{max} vs. z_{max} diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:6 orbits identified with the visual method (see 4.2.2). 

In the text 
Fig. A.9. Frequency map in cylindrical coordinates for different families (identified with the visual method; see Section.4.2.2) with computable cylindrical frequencies. 

In the text 
Fig. B.1. R_{max} vs. z_{max} diagram comparison for the axisymmetric and nonaxisymmetric cases. Left: Reference diagram, with the nonaxisymmetric combination of potentials presented in Section 3.1. Right: Using the axisymmetric potential MWPotential14 from Bovy (2015). 

In the text 
Fig. B.2. R_{max} vs. z_{max} diagram comparison with different NSD potentials. Upper panel: Reference diagram, with the NSD potential from Sormani et al. (2020) (model 3) used in our study. Middle panel: Using the NSD potential from Sormani et al. (2022). Lower panel: Using the NSD potential from Launhardt et al. (2002). In the three cases, the bar and the NSC potentials are not changed. 

In the text 
Fig. B.3. R_{max} vs. z_{max} diagram comparison with different bar potentials. Left: Reference diagram, with the bar potential from Launhardt et al. (2002). Right: Using the bar potential from Portail et al. (2017). In both cases, the NSD and NSC potentials are not changed. 

In the text 
Fig. B.4. R_{max} vs. z_{max} diagram comparison with a different bar angle with respect to the line of sight towards the GC. Left: Reference diagram as presented in Section 3.1. Right: Case with a 10% variation of the angle. 

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