Orbital analysis of stars in the nuclear stellar disc of the Milky Way

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 non-axisymmetric 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, z -tube, x -tube, 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 z -tubes (or as a sub-family of z -tubes). 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.


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. 2023), 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 counterrotation for the most metal-poor stars.
Extragalactic studies showed that many barred galaxies host nuclear discs or rings (Gadotti et al. 2019(Gadotti et al. , 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 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.

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 VIRAC2 (Smith et al. in prep) photometric and astrometric reduction of the VVV data (Minniti et al. 2010).These proper motions were rescaled to the Gaia absolute reference frame (Sanders et al. 2019).We also refer here to Sormani et al. (2022).Our sample includes only stars that are primary sources of the survey leaving us a total of 2501 stars.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 Tab.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 Tab. 2 of Sormani et al. 2022).In addition, we used the same colour cut (H − Ks) > max(1.3,−0.0233K + 1.63) as in Schultheis et al. (2021) to remove foreground stars.

Sormani
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 98 th 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 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 Section 5).To overcome this, we plan to address this issue in a forthcoming paper by using N-body simulations of the NSD.

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 et al. 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 Sections 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 Equation 27); (iii) We used the NSC density from Chatzopoulos et al. (2015) (see their Equation 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 −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 Section 3.2).We chose a short timestep of 4000 yr for a good sampling of the orbits.

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 , Ω 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 long-axis)1 , and box orbits (chaotic orbits correspond to orbits that cannot be identified with one of these categories, see Section 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.

Chaoticity
To study chaoticity and particularly identify chaotic orbits, we used two different methods that we introduce in this section.Highest Lyapunov exponent

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 Section 4.3 of Vasiliev (2013).

Frequency drift
Another way to study chaoticity is to directly use the orbital frequencies.Valluri et al. (2010) (see their Section 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.

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.
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 Section 3).The key characteristics mentioned earlier clearly persist.This indicates that the uncertainty in the distances of NSD stars can safely be disregarded.
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), andSormani 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).

Orbit classification
We discuss and compare the two methods we used to classify orbits in various orbital families.

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 Section 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 Section 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. 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.As explained Section 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)   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 Fig. 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.

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 Section 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 Section.4.2.1.As explained in Section 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 figures A.1 to 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 Section 4.   Notes.1: the "z-tube" family contains all orbits with a circulation around the z-axis and includes therefore fish, saucer, pretzel, 5 : 4, 5 : 6 orbits.
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.

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 As explained Section 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 Other z-tubes Chaotic Saucer Fish/anti-fish x-tube Banana/anti-banana Pretzel/anti-pretzel 5:4 5:6 (l,b) plane shows that these stars are mostly located in the outer parts of the NSD (|l| > 1.0 o ), 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).
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(Valluri et al. , 2012(Valluri et al. , 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.
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 Tab.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.
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 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.
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 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 2 -type orbital family dominates (Athanassoula 1992a,b, Li et al. 2015, Sormani et al. 2018, Sormani et al. 2023).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 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 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 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 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éndez-Abreu et al. 2019), but a more in-depth study is needed.

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

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

Fig. 6 .
Fig. 5. R max vs. z max 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.

Fig. 7 .
Fig. 7. R max vs. z max 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.

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

Fig. 9 .
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.
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 (seeFig.A.6,A.7,A.8 respectively).

Fig. 10 .
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.
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 Fig.10 and Fig.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

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

Table 1 .
Family membership results for the automatic and visual methods