© The Authors 2024

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, Gallego-Cano et al. 2020). Increasing evidence is reported that the NSD is a distinct structure from the NSC and the nuclear bulge: Nogueras-Lara et al. (2020) determined the SFH in the NSD using the GALACTICNUCLEUS data (Nogueras-Lara 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. Nogueras-Lara 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 metal-rich stars may have formed in the central molecular zone, while metal-poor 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 metal-rich stars rotate faster than metal-poor stars, with some hints of counter-rotation for the most metal-poor 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 inside-out 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 self-consistent 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 so-called 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⁻¹ and μ_b,err < 0.6 mas yr⁻¹ (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 colour-magnitude 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₀ < 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 N-body 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, non-axisymmetric, 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. Gallego-Cano et al. 2020, Sormani et al. 2022, Nogueras-Lara 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, Nogueras-Lara 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 Nogueras-Lara 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 self-consistent 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 non-axisymmetric 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 (clock-wise) rotation pattern speed of Ω_b = 40 km s⁻¹ kpc⁻¹ (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 so-called numerical approximation of fundamental frequencies (NAFF), which has been applied in planetary dynamics. Price-Whelan (2015) developed an open-source 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 phase-space 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 (short-axis or long-axis)¹, 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, short-axis and long-axis tube orbits in the Cartesian case are clustered to a trivial resonance (e.g. 1:−1:0 resonance for short-axis 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 non-linear 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 long-term 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 long-term 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:

$$\begin{array}{c}\hfill log(\mathrm{\Delta }{f}_i)=log|\frac{{\mathrm{\Omega }}_i(t_1)-{\mathrm{\Omega }}_i(t_2)}{{\mathrm{\Omega }}_i(t_1)}|\end{array}$$(1)

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 colour-coded by the eccentricity.

Fig. 5. Rmax vs. zmax diagram colour-coded 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 Rmax vs. zmax diagram (colour-coded 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⁻¹, 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 counter-rotating stars, which have been identified previously by Schultheis et al. (2021).

Fig. 7. Rmax vs. zmax diagram colour-coded 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 x-axis results in an x-tube) 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 short-axis tubes (called here after “z-tubes”) 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 long-axis tubes, called here after x-tubes, are located. An example of the morphology of an x-tube 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, z-tube, x-tube), 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 z-tubes (because our sample contains very few x-tubes, 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 z-tubes alone (orbits for which a circulation around the z-axis 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.

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, x-tube, z-tube, 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 anti-banana,-fish, and -pretzel orbits. We are also able to identify orbits precisely whose own well-known 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) x-tubes 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 anti-banana) 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 anti-fish) 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 anti-pretzel) 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 non-axisymmetric 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 non-integrable potentials with some indication of chaotic orbits (Price-Whelan 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 (non-axisymmetric 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: Rmax vs. zmax vs. Lyapunov. Lower panel: Rmax vs. zmax 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 z-tubes that is expected for stars forming a disc structure such as the NSD and a minority of x-tubes. In addition, we were able to identify a variety of families that mostly belong to z-tubes. 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 z-tube, x-tube, 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 K-nearest neighbours to classify galaxies in, and around clusters according to their projected phase-space position. In order to obtain a precise classification, a large training set is necessary on which N-body 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 x-tubes 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. Rmax vs. zmax 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 Rmax and maximum height zmax for stars with a different minimum eccentricity threshold. For instance, the median Rmax for a minimum eccentricity of 0.5 means that we only computed the median of the Rmax 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. Nogueras-Lara et al. 2023) shows that our working sample is indeed highly incomplete (∼15% completeness at K_S ∼ 11.5 mag).

Complementary high-resolution N-body 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₂-type, and orbits parallel to the bar are called x₁-type orbits. x₂-type orbits are parents of the z-tube family, occupying the 1:−1:0 main resonance line in the frequency map (see Fig. 8), where about two-thirds 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₂-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₂ 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 metal-rich 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₁, x₂, and x₄ share the 1 : 2 resonance between the tangential oscillation frequency (Ω_ϕ) and the radial frequency (Ω_R). Unfortunately, they are not anticipated to be abundant in N-body 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₁-tree (Athanassoula 2005) with banana orbits that were called “x₁v₁” 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éndez-Abreu et al. 2019), but a more in-depth 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, z-tube, x-tube, banana, fish, saucer, pretzel, 5:4, and 5:6 resonances. The large majority of sources are z-tubes, 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 z-tube orbits that are parented by the x2-type orbits. This is indeed expected if the formation of the NSD is coupled with the formation of the Galactic bar where x2-type orbits are the dominant population. As a follow-up work, a comparison with self-consistent 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. Price-Whelan, E. Vasiliev for fruitful discussions. R. Schödel acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033, from grant EUR2022-134031 funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR, and from grant PID2022-136640NB-C21 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 (Price-Whelan 2015), Gala (Price-Whelan et al. 2022).

References

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:Rmax vs. zmax 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 x-tube, 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the x-tube 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and anti-banana 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the fish and anti-fish 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:Rmax vs. zmax 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the pretzel and anti-pretzel 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:Rmax vs. zmax 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:Rmax vs. zmax 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 (non-axisymmetric) 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. Rmax vs. zmax diagram comparison for the axisymmetric and non-axisymmetric cases. Left: Reference diagram, with the non-axisymmetric combination of potentials presented in Section 3.1. Right: Using the axisymmetric potential MWPotential14 from Bovy (2015).

Fig. B.2. Rmax vs. zmax 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. Rmax vs. zmax 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. Rmax vs. zmax 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. Rmax vs. zmax diagram colour-coded 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 Rmax vs. zmax diagram (colour-coded with eccentricity) by using 100 MCMC distances. The standard deviation for each star is indicated with the error bars.
In the text

	Fig. 7. Rmax vs. zmax diagram colour-coded 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 z-tubes alone (orbits for which a circulation around the z-axis 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: Rmax vs. zmax vs. Lyapunov. Lower panel: Rmax vs. zmax 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. Rmax vs. zmax 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 Rmax and maximum height zmax for stars with a different minimum eccentricity threshold. For instance, the median Rmax for a minimum eccentricity of 0.5 means that we only computed the median of the Rmax 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:Rmax vs. zmax 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 x-tube, 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the x-tube 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and anti-banana 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the fish and anti-fish 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:Rmax vs. zmax 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:Rmax vs. zmax diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the pretzel and anti-pretzel 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:Rmax vs. zmax 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:Rmax vs. zmax 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. Rmax vs. zmax diagram comparison for the axisymmetric and non-axisymmetric cases. Left: Reference diagram, with the non-axisymmetric combination of potentials presented in Section 3.1. Right: Using the axisymmetric potential MWPotential14 from Bovy (2015).
In the text

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