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

N. Nieuwmunster, M. Schultheis, M. Sormani, F. Fragkoudi, F.

Nogueras-Lara, R. Schödel, P. Mcmillan, L. C. Smith, J. L. Sanders

► To cite this version:

N. Nieuwmunster, M. Schultheis, M. Sormani, F. Fragkoudi, F. Nogueras-Lara, et al. Orbital analysis of stars in the nuclear stellar disc of the Milky Way. Astronomy and Astrophysics - A&A, 2024, 685, pp.A93. 10.1051/0004-6361/202349000. hal-04571845

HAL Id: hal-04571845

https://hal.science/hal-04571845

Submitted on 9 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Orbital analysis of stars in the nuclear stellar disc of the Milky Way

N. Nieuwmunster, M. Schultheis, M. Sormani, F. Fragkoudi, F.

Nogueras-Lara, R. Schödel, P. Mcmillan, L. C. Smith, J. L. Sanders

► To cite this version:

N. Nieuwmunster, M. Schultheis, M. Sormani, F. Fragkoudi, F. Nogueras-Lara, et al. Orbital analysis of stars in the nuclear stellar disc of the Milky Way. Astronomy & Astrophysics - A&A, 2024, 685, pp.A93. 10.1051/0004-6361/202349000. hal-04571845

HAL Id: hal-04571845

https://hal.science/hal-04571845v1

Submitted on 9 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Orbital analysis of stars in the nuclear stellar disc of the Milky Way

N. Nieuwmunster<fnref n="1" />,<fnref n="2" />, M. Schultheis<fnref n="1" />, M. Sormani<fnref n="3" />, F. Fragkoudi<fnref n="4" />, F. Nogueras-Lara<fnref n="5" />, R. Schödel<fnref n="6" />, P. McMillan<fnref n="2" />,<fnref n="7" />, L. C. Smith<fnref n="8" />, and J. L. Sanders<fnref n="9" />

Received 18 December 2023 / Accepted 21 February 2024

## ABSTRACT

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

**Key words.** stars: kinematics and dynamics – Galaxy: bulge – Galaxy: center – Galaxy: nucleus – Galaxy: structure

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

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, Laboratoire Lagrange, CNRS, Blvd de l'Observatoire, 06304 Nice, France e-mail: niels.nieuwmunster@oca.eu</footnote>

<footnote n="2">² Division of Astrophysics, Department of Physics, Lund University, Box 43, 22100 Lund, Sweden</footnote>

<footnote n="3">³ Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote n="4">⁴ Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5">⁵ Max-Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain</footnote>

<footnote n="7">⁷ School of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote n="8">⁸ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University College London, London WC1E 6BT, UK</footnote>
# Orbital analysis of stars in the nuclear stellar disc of the Milky Way

N. Nieuwmunster<fnref n="1" />,<fnref n="2" />, M. Schultheis<fnref n="1" />, M. Sormani<fnref n="3" />, F. Fragkoudi<fnref n="4" />, F. Nogueras-Lara<fnref n="5" />, R. Schödel<fnref n="6" />, P. McMillan<fnref n="2" />,<fnref n="7" />, L. C. Smith<fnref n="8" />, and J. L. Sanders<fnref n="9" />

Received 18 December 2023 / Accepted 21 February 2024

## ABSTRACT

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

**Key words.** stars: kinematics and dynamics – Galaxy: bulge – Galaxy: center – Galaxy: nucleus – Galaxy: structure

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

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, Laboratoire Lagrange, CNRS, Blvd de l'Observatoire, 06304 Nice, France e-mail: niels.nieuwmunster@oca.eu</footnote>

<footnote n="2">² Division of Astrophysics, Department of Physics, Lund University, Box 43, 22100 Lund, Sweden</footnote>

<footnote n="3">³ Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote n="4">⁴ Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5">⁵ Max-Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain</footnote>

<footnote n="7">⁷ School of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote n="8">⁸ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University College London, London WC1E 6BT, UK</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Upper panel: histograms of the radial velocities and proper motions in $\mu_l$ and $\mu_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.</figcaption></figure>

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^{\circ}$ and $|b| < 0.3^{\circ}$, 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 - K_s$) > 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: $\mu_{l,\text{err}} < 0.6$ mas yr$^{-1}$ and $\mu_{b,\text{err}} < 0.6$ mas yr$^{-1}$ (corresponding to the 98th percentile; see the vertical dashed line in Fig. 1), leading to a total sample of 1130 stars. As pointed out by Sormani et al. (2022), this proper motion cut also removes very bright stars ($K < 10$), for which the proper motion errors become large due to saturation effects. Figure 2 shows the corresponding *HK* 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 Sect. 5). To overcome this, we plan to address this issue in a forthcoming paper by using *N*-body simulations of the NSD.

## 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 dis-
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Upper panel: histograms of the radial velocities and proper motions in $\mu_l$ and $\mu_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.</figcaption></figure>

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^\circ$ and $|b| < 0.3^\circ$, 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 - K_s$) > 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: $\mu_{l,\text{err}} < 0.6$ mas yr$^{-1}$ and $\mu_{b,\text{err}} < 0.6$ mas yr$^{-1}$ (corresponding to the 98th percentile; see the vertical dashed line in Fig. 1), leading to a total sample of 1130 stars. As pointed out by Sormani et al. (2022), this proper motion cut also removes very bright stars ($K < 10$), for which the proper motion errors become large due to saturation effects. Figure 2 shows the corresponding *HK* 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 Sect. 5). To overcome this, we plan to address this issue in a forthcoming paper by using *N*-body simulations of the NSD.

## 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
<figure><img src="image_3.png" /><figcaption>**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.0233*K* + 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.</figcaption></figure>

tribution 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 distribu-
tion (μ = 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 clos-
est 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 cor-
respond 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 ρ<sub>tot</sub> = ρ<sub>bar</sub> + ρ<sub>NSD</sub> + ρ<sub>NSC</sub>
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 com-
parison with other potentials in Appendix B.

We took into account a typical (clock-wise) rotation pattern
speed of Ω<sub>b</sub> = 40 km s<sup>-1</sup> kpc<sup>-1</sup> (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 funda-
mental frequencies, Ω, that determine the periodic behaviour of
motion (Binney & Tremaine 2008). This motion can be decom-
posed as a Fourier sum, where the Fourier frequencies are lin-
ear combinations of the fundamental frequencies. These funda-
mental frequencies can be recovered, as shown by Laskar
(1993), using the so-called numerical approximation of funda-
mental frequencies (NAFF), which has been applied in plane-
tary dynamics. Price-Whelan (2015) developed an open-source
code called SuperFreq, which is a Python implementation sim-
ilar to the NAFF code. This code finds the fundamental fre-
quencies (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 res-
onances. Resonant orbits are those for which the fundamental
frequencies Ω = Ωₓ, Ωᵧ, Ωz (or ΩR, ΩΦ, Ωz depending on the
selected coordinate system) can be measured and where n·Ω = 0
for n = (nₓ, nᵧ, nz), where the vector n only contains integer num-
bers and at most one zero (Koppelman et al. 2021). For example,
in Cartesian coordinates, nₓΩₓ + nᵧΩᵧ + nzΩz = 0, and the reso-
nance is then called nₓ:nᵧ:nz. Moreover, boxlets are special cases
of resonant orbits in which one of the integers nₓ, nᵧ, nz is zero.

Frequency maps are a powerful tool for obtaining an auto-
matic 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)<fnref n="1" />, and box orbits (chaotic orbits correspond to orbits that cannot be identified with one of these categories, see Sect. 3.3). Based on the frequencies in Cartesian coordinates (Ωₓ, Ωᵧ, Ω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.

<footnote n="1">¹ Tube orbits circulate around a certain axis: *z* for short-axis tubes and *x* for long-axis tubes.</footnote>
<figure><img src="image_3.png" /><figcaption>**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.0233*K* + 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.</figcaption></figure>

tribution 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 distribu-
tion (μ = 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 clos-
est 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 cor-
respond 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 ρ<sub>tot</sub> = ρ<sub>bar</sub> + ρ<sub>NSD</sub> + ρ<sub>NSC</sub>
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 com-
parison with other potentials in Appendix B.

We took into account a typical (clock-wise) rotation pattern
speed of Ω<sub>b</sub> = 40 km s<sup>-1</sup> kpc<sup>-1</sup> (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 funda-
mental frequencies, Ω, that determine the periodic behaviour of
motion (Binney & Tremaine 2008). This motion can be decom-
posed as a Fourier sum, where the Fourier frequencies are lin-
ear combinations of the fundamental frequencies. These funda-
mental frequencies can be recovered, as shown by Laskar
(1993), using the so-called numerical approximation of funda-
mental frequencies (NAFF), which has been applied in plane-
tary dynamics. Price-Whelan (2015) developed an open-source
code called SuperFreq, which is a Python implementation sim-
ilar to the NAFF code. This code finds the fundamental fre-
quencies (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 res-
onances. Resonant orbits are those for which the fundamental
frequencies Ω = Ωₓ, Ωᵧ, Ωz (or ΩR, ΩΦ, Ωz depending on the
selected coordinate system) can be measured and where n·Ω = 0
for n = (nₓ, nᵧ, nz), where the vector n only contains integer num-
bers and at most one zero (Koppelman et al. 2021). For example,
in Cartesian coordinates, nₓΩₓ + nᵧΩᵧ + nzΩz = 0, and the reso-
nance is then called nₓ:nᵧ:nz. Moreover, boxlets are special cases
of resonant orbits in which one of the integers nₓ, nᵧ, nz is zero.

Frequency maps are a powerful tool for obtaining an auto-
matic 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)<fnref n="1" />, and box orbits (chaotic orbits correspond to orbits that cannot be identified with one of these categories, see Sect. 3.3). Based on the frequencies in Cartesian coordinates (Ωₓ, Ωᵧ, Ω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.

<footnote n="1">¹ Tube orbits circulate around a certain axis: *z* for short-axis tubes and *x* for long-axis tubes.</footnote>
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

### 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,

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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:

$$ \log(\Delta f_i) = \log \left| \frac{\Omega_i(t_1) - \Omega_i(t_2)}{\Omega_i(t_1)} \right| \quad (1) $$

where $i$ defines the frequency component in Cartesian coordinates (i.e. $\log(\Delta f_x)$, $\log(\Delta f_y)$ and $\log(\Delta f_z)$). The highest value of the three frequency drift parameters $\log(\Delta f_i)$ is then associated with the frequency drift parameter $\log(\Delta f)$. The higher the value of $\log(\Delta 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(\Delta f)$. In order to define a threshold value of $\log(\Delta 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(\Delta 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.

## 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.
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

### 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,

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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:

$$ \log(\Delta f_i) = \log \left| \frac{\Omega_i(t_1) - \Omega_i(t_2)}{\Omega_i(t_1)} \right| \quad (1) $$

where $i$ defines the frequency component in Cartesian coordinates (i.e. $\log(\Delta f_x)$, $\log(\Delta f_y)$ and $\log(\Delta f_z)$). The highest value of the three frequency drift parameters $\log(\Delta f_i)$ is then associated with the frequency drift parameter $\log(\Delta f)$. The higher the value of $\log(\Delta 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(\Delta f)$. In order to define a threshold value of $\log(\Delta 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(\Delta 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.

## 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.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Average $R_{max}$ vs. $z_{max}$ diagram (colour-coded with eccentricity) by using 100 MCMC distances. The standard deviation for each star is indicated with the error bars.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** $R_{max}$ vs. $z_{max}$ diagram colour-coded with the rotational velocity. $V_{\phi} > 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.</figcaption></figure>

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.

Figure 7 shows the same feature, but as a function of rotational velocity $V_{\phi}$. 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).

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

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 $\Omega_{\phi}/\Omega_{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 $\Omega_{z}$ and the radial oscillation frequency $\Omega_{R}$ and vertical lines between the azimuthal oscillation frequency $\Omega_{\phi}$ and the radial oscillation frequency $\Omega_{R}$). We discern then some resonances, $\Omega_{z} : \Omega_{R} = 1:1, 3:4, 1:2$, and $\Omega_{\phi}:\Omega_{R} = 1:2$. Saucer orbits (Sambhus & Sridhar 2000; Vasiliev 2014) (see an example Fig. A.5) were found via the resonance $\Omega_{z}:\Omega_{R} = 1:1$ (Yavetz et al. 2023), and we therefore have an estimate of their number.

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 $\pm 0.01$ (in frequency ratio), which is enough for contamination by sticky chaotic orbits. However, we
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Average $R_{max}$ vs. $z_{max}$ diagram (colour-coded with eccentricity) by using 100 MCMC distances. The standard deviation for each star is indicated with the error bars.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** $R_{max}$ vs. $z_{max}$ diagram colour-coded with the rotational velocity. $V_{\phi} > 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.</figcaption></figure>

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.

Figure 7 shows the same feature, but as a function of rotational velocity $V_{\phi}$. 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).

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

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 $\Omega_{\phi}/\Omega_{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 $\Omega_{z}$ and the radial oscillation frequency $\Omega_{R}$ and vertical lines between the azimuthal oscillation frequency $\Omega_{\phi}$ and the radial oscillation frequency $\Omega_{R}$). We discern then some resonances, $\Omega_{z} : \Omega_{R} = 1:1, 3:4, 1:2$, and $\Omega_{\phi}:\Omega_{R} = 1:2$. Saucer orbits (Sambhus & Sridhar 2000; Vasiliev 2014) (see an example Fig. A.5) were found via the resonance $\Omega_{z}:\Omega_{R} = 1:1$ (Yavetz et al. 2023), and we therefore have an estimate of their number.

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 $\pm 0.01$ (in frequency ratio), which is enough for contamination by sticky chaotic orbits. However, we
<figure><img src="image_9.png" /></figure>

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

<figure><img src="image_10.png" /></figure>

**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 Ω<sub>z</sub> and Ω<sub>R</sub> and between Ω<sub>φ</sub> and Ω<sub>R</sub> 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.
<figure><img src="image_9.png" /></figure>

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

<figure><img src="image_10.png" /></figure>

**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 Ω<sub>z</sub> and Ω<sub>R</sub> and between Ω<sub>φ</sub> and Ω<sub>R</sub> 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.
<table><caption>**Table 1.** Family membership results for the automatic and visual methods</caption><thead><tr><th></th><th></th><th>Chaotic/Box</th><th>x-tube</th><th>z-tube <sup>(1)</sup></th><th>Banana</th><th>Fish</th><th>Saucer</th><th>Pretzel</th><th>5:4</th><th>5:6</th></tr></thead><tbody><tr><td rowspan="2">Automatic</td><td>Total</td><td>34.0%</td><td>1.1%</td><td>64.9%</td><td>4.7%</td><td>9.0%</td><td>13.6%</td><td>3.5%</td><td>1.5%</td><td>1.6%</td></tr><tr><td>With identified chaotics</td><td>—</td><td>2.1%</td><td>—</td><td>—</td><td>11.2%</td><td>—</td><td>5.0%</td><td>1.9%</td><td>1.8%</td></tr><tr><td rowspan="3">Visual</td><td>Total</td><td>24.7%</td><td>1.5%</td><td>68.2%</td><td>5.7%</td><td>11.2%</td><td>8.7%</td><td>4.5%</td><td>0.6%</td><td>1.1%</td></tr><tr><td>Centrophobic</td><td>—</td><td>—</td><td>—</td><td>4.0%</td><td>9.6%</td><td>—</td><td>1.8%</td><td>—</td><td>—</td></tr><tr><td>Centrophilic</td><td>—</td><td>—</td><td>—</td><td>1.7%</td><td>1.6%</td><td>—</td><td>2.6%</td><td>—</td><td>—</td></tr></tbody></table>

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

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
<table><caption>**Table 1.** Family membership results for the automatic and visual methods</caption><thead><tr><th></th><th></th><th>Chaotic/Box</th><th>x-tube</th><th>z-tube <sup>(1)</sup></th><th>Banana</th><th>Fish</th><th>Saucer</th><th>Pretzel</th><th>5:4</th><th>5:6</th></tr></thead><tbody><tr><td rowspan="2">Automatic</td><td>Total</td><td>34.0%</td><td>1.1%</td><td>64.9%</td><td>4.7%</td><td>9.0%</td><td>13.6%</td><td>3.5%</td><td>1.5%</td><td>1.6%</td></tr><tr><td>With identified chaotics</td><td>—</td><td>2.1%</td><td>—</td><td>—</td><td>11.2%</td><td>—</td><td>5.0%</td><td>1.9%</td><td>1.8%</td></tr><tr><td rowspan="3">Visual</td><td>Total</td><td>24.7%</td><td>1.5%</td><td>68.2%</td><td>5.7%</td><td>11.2%</td><td>8.7%</td><td>4.5%</td><td>0.6%</td><td>1.1%</td></tr><tr><td>Centrophobic</td><td>—</td><td>—</td><td>—</td><td>4.0%</td><td>9.6%</td><td>—</td><td>1.8%</td><td>—</td><td>—</td></tr><tr><td>Centrophilic</td><td>—</td><td>—</td><td>—</td><td>1.7%</td><td>1.6%</td><td>—</td><td>2.6%</td><td>—</td><td>—</td></tr></tbody></table>

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

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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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, 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 learn-

ing (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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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, 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 learn-

ing (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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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.

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_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, 2024). After the Galactic bar was formed, gas was funnelled along the bar towards the GC, where it settled down to form a nuclear disc. The gas then formed stars, which keep their resulting $x_2$ orbits, and the resulting stellar population resembles a disc. Additional evidence comes from the spatial and kinematical overlap between the central molecular zone and the NSD (Schönrich et al. 2015, Schultheis et al. 2021) in the Milky Way, which supports this scenario that the stars in the NSD should be more 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 ($\Omega_\phi$) and the radial frequency ($\Omega_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_1v_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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

<caption>**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.</caption>

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.

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_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 for-

mation of nuclear stellar discs is strongly connected to the properties of the galactic bar. Hydrodynamical simulations suggest that nuclear stellar discs form close to the inner Lindblad resonance of the main bar, where the $x_2$-type orbital family dominates (Athanassoula 1992a,b, Li et al. 2015, Sormani et al. 2018, 2024). After the Galactic bar was formed, gas was funnelled along the bar towards the GC, where it settled down to form a nuclear disc. The gas then formed stars, which keep their resulting $x_2$ orbits, and the resulting stellar population resembles a disc. Additional evidence comes from the spatial and kinematical overlap between the central molecular zone and the NSD (Schönrich et al. 2015, Schultheis et al. 2021) in the Milky Way, which supports this scenario that the stars in the NSD should be more 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 ($\Omega_\phi$) and the radial frequency ($\Omega_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_1v_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.

## 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 benefit from the large number of stars to establish the strong connection between dynamics and chemistry.

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

Amarante, J. A. S., Smith, M. C., & Boeche, C. 2020, MNRAS, 492, 3816

Athanassoula, E. 1992a, MNRAS, 259, 328

Athanassoula, E. 1992b, MNRAS, 259, 345

Athanassoula, E. 2005, Ann. New York Acad. Sci., 1045, 168

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition, (Princeton, NJ USA: Princeton University Press)

Bittner, A., Sánchez-Blázquez, P., Gadotti, D. A., et al. 2020, A&A, 643, A65

Bovy, J. 2015, ApJS, 216, 29

Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948

de los Rios, M., Martínez, H. J., Coenda, V., et al. 2021, MNRAS, 500, 1784

Fritz, T. K., Patrick, L., Feldemier-Krause, A., et al. 2021, A&A, 649, A83

Gadotti, D. A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2019, MNRAS, 482, 506

Gadotti, D. A., Bittner, A., Falcón-Barroso, J., et al. 2020, A&A, 643, A14

Gallego-Cano, E., Schödel, R., Nogueras-Lara, F., et al. 2020, A&A, 634, A71

Gonzalez, O. A., Mucciarelli, A., Origlia, L., et al. 2020, The Messenger, 180, 18

GRAVITY Collaboration (Abuter, R., et al.) 2020, A&A, 636, L5

Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113

Koppelman, H. H., Hagen, J. H. J., & Helmi, A. 2021, A&A, 647, A37

Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G., et al. 2014, MNRAS, 440, 3370

Laskar, J. 1993, Celestial Mech. Dyn. Astron., 56, 191

Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112

Li, Z., Shen, J., & Kim, W.-T. 2015, ApJ, 806, 150

Lindqvist, M., Habing, H. J., & Winnberg, A. 1992, A&A, 259, 118

Lyapunov, A. M. 1992, Int. J. Control, 55, 531

Méndez-Abreu, J., de Lorenzo-Cáceres, A., Gadotti, D. A., et al. 2019, MNRAS, 482, L118

Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433

Nishiyama, S., Yasui, K., Nagata, T., et al. 2013, ApJ, 769, L28

Nogueras-Lara, F. 2022, A&A, 668, L8

Nogueras-Lara, F., Gallego-Calvente, A. T., Dong, H., et al. 2018, A&A, 610, A83

Nogueras-Lara, F., Schödel, R., Gallego-Calvente, A. T., et al. 2020, Nat. Astron., 4, 377

Nogueras-Lara, F., Schultheis, M., Najarro, F., et al. 2023, A&A, 671, L10

Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233

Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 2002, ApJ, 565, 265

Price-Whelan, A. M. 2015, https://doi.org/10.5281/zenodo.18787

Price-Whelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079

Price-Whelan, A., Sipőcz, B., Starkman, N., et al. 2022, https://doi.org/10.5281/zenodo.6325733

Sambhus, N., & Sridhar, S. 2000, ApJ, 542, 143

Sanders, J. L., Smith, L., Evans, N. W., & Lucas, P. 2019, MNRAS, 487, 5188

Sanders, J.L., Kawata, D., Matsunaga, N., et al. 2024, MNRAS, 530, 2972

Schönrich, R., Aumer, M., & Sale, S. E. 2015, ApJ, 812, L21

Schultheis, M., Fritz, T. K., Nandakumar, G., et al. 2021, A&A, 650, A191

Shahzamanian, B., Schödel, R., Nogueras-Lara, F., et al. 2022, A&A, 662, A11

Sharples, R., Bender, R., Agudo Berbel, A., et al. 2013, The Messenger, 151, 21

Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847

Smith, L. C., Lucas, P. W., Kurtev, R., et al. 2018, MNRAS, 474, 1826

Sormani, M. C., Treß, R. G., Ridley, M., et al. 2018, MNRAS, 475, 2383

Sormani, M. C., Magorrian, J., Nogueras-Lara, F., et al. 2020, MNRAS, 499, 7

Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022, MNRAS, 512, 1857

Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, MNRAS, 528, 5762

Valluri, M., & Merritt, D. 1998, ApJ, 506, 686

Valluri, M., Debattista, V. P., Quinn, T., & Moore, B. 2010, MNRAS, 403, 525

Valluri, M., Debattista, V. P., Quinn, T. R., Roškar, R., & Wadsley, J. 2012, MNRAS, 419, 1951

Valluri, M., Shen, J., Abbott, C., & Debattista, V. P. 2016, ApJ, 818, 141

Vasiliev, E. 2013, MNRAS, 434, 3174

Vasiliev, E. 2014, CQG, 31, 244002

Vasiliev, E. 2019, MNRAS, 482, 1525

Yavetz, T. D., Johnston, K. V., Pearson, S., Price-Whelan, A. M., & Hamilton, C. 2023, ApJ, 954, 215
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 benefit from the large number of stars to establish the strong connection between dynamics and chemistry.

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

Amarante, J. A. S., Smith, M. C., & Boeche, C. 2020, MNRAS, 492, 3816

Athanassoula, E. 1992a, MNRAS, 259, 328

Athanassoula, E. 1992b, MNRAS, 259, 345

Athanassoula, E. 2005, Ann. New York Acad. Sci., 1045, 168

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition, (Princeton, NJ USA: Princeton University Press)

Bittner, A., Sánchez-Blázquez, P., Gadotti, D. A., et al. 2020, A&A, 643, A65

Bovy, J. 2015, ApJS, 216, 29

Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948

de los Rios, M., Martínez, H. J., Coenda, V., et al. 2021, MNRAS, 500, 1784

Fritz, T. K., Patrick, L., Feldemier-Krause, A., et al. 2021, A&A, 649, A83

Gadotti, D. A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2019, MNRAS, 482, 506

Gadotti, D. A., Bittner, A., Falcón-Barroso, J., et al. 2020, A&A, 643, A14

Gallego-Cano, E., Schödel, R., Nogueras-Lara, F., et al. 2020, A&A, 634, A71

Gonzalez, O. A., Mucciarelli, A., Origlia, L., et al. 2020, The Messenger, 180, 18

GRAVITY Collaboration (Abuter, R., et al.) 2020, A&A, 636, L5

Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113

Koppelman, H. H., Hagen, J. H. J., & Helmi, A. 2021, A&A, 647, A37

Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G., et al. 2014, MNRAS, 440, 3370

Laskar, J. 1993, Celestial Mech. Dyn. Astron., 56, 191

Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112

Li, Z., Shen, J., & Kim, W.-T. 2015, ApJ, 806, 150

Lindqvist, M., Habing, H. J., & Winnberg, A. 1992, A&A, 259, 118

Lyapunov, A. M. 1992, Int. J. Control, 55, 531

Méndez-Abreu, J., de Lorenzo-Cáceres, A., Gadotti, D. A., et al. 2019, MNRAS, 482, L118

Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433

Nishiyama, S., Yasui, K., Nagata, T., et al. 2013, ApJ, 769, L28

Nogueras-Lara, F. 2022, A&A, 668, L8

Nogueras-Lara, F., Gallego-Calvente, A. T., Dong, H., et al. 2018, A&A, 610, A83

Nogueras-Lara, F., Schödel, R., Gallego-Calvente, A. T., et al. 2020, Nat. Astron., 4, 377

Nogueras-Lara, F., Schultheis, M., Najarro, F., et al. 2023, A&A, 671, L10

Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233

Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 2002, ApJ, 565, 265

Price-Whelan, A. M. 2015, https://doi.org/10.5281/zenodo.18787

Price-Whelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079

Price-Whelan, A., Sipőcz, B., Starkman, N., et al. 2022, https://doi.org/10.5281/zenodo.6325733

Sambhus, N., & Sridhar, S. 2000, ApJ, 542, 143

Sanders, J. L., Smith, L., Evans, N. W., & Lucas, P. 2019, MNRAS, 487, 5188

Sanders, J.L., Kawata, D., Matsunaga, N., et al. 2024, MNRAS, 530, 2972

Schönrich, R., Aumer, M., & Sale, S. E. 2015, ApJ, 812, L21

Schultheis, M., Fritz, T. K., Nandakumar, G., et al. 2021, A&A, 650, A191

Shahzamanian, B., Schödel, R., Nogueras-Lara, F., et al. 2022, A&A, 662, A11

Sharples, R., Bender, R., Agudo Berbel, A., et al. 2013, The Messenger, 151, 21

Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847

Smith, L. C., Lucas, P. W., Kurtev, R., et al. 2018, MNRAS, 474, 1826

Sormani, M. C., Treß, R. G., Ridley, M., et al. 2018, MNRAS, 475, 2383

Sormani, M. C., Magorrian, J., Nogueras-Lara, F., et al. 2020, MNRAS, 499, 7

Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022, MNRAS, 512, 1857

Sormani, M. C., Sobacchi, E., & Sanders, J. L. 2024, MNRAS, 528, 5762

Valluri, M., & Merritt, D. 1998, ApJ, 506, 686

Valluri, M., Debattista, V. P., Quinn, T., & Moore, B. 2010, MNRAS, 403, 525

Valluri, M., Debattista, V. P., Quinn, T. R., Roškar, R., & Wadsley, J. 2012, MNRAS, 419, 1951

Valluri, M., Shen, J., Abbott, C., & Debattista, V. P. 2016, ApJ, 818, 141

Vasiliev, E. 2013, MNRAS, 434, 3174

Vasiliev, E. 2014, CQG, 31, 244002

Vasiliev, E. 2019, MNRAS, 482, 1525

Yavetz, T. D., Johnston, K. V., Pearson, S., Price-Whelan, A. M., & Hamilton, C. 2023, ApJ, 954, 215
**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*<sub>max</sub> vs. *z*<sub>max</sub> diagram and in the fre-

quency map in Cartesian coordinates. We also show in Fig.A.9 the location of certain families in the cylindrical frequency map.

<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Example of a chaotic orbit. *Upper panel:* Orbit plotted in the (x, y), (x, z) and (y, z) planes. *Lower panel:* *R*<sub>max</sub> vs. *z*<sub>max</sub> diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the chaotic orbits identified with the visual method (see 4.2.2).</figcaption></figure>
**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*<sub>max</sub> vs. *z*<sub>max</sub> diagram and in the fre-

quency map in Cartesian coordinates. We also show in Fig.A.9 the location of certain families in the cylindrical frequency map.

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Example of a chaotic orbit. *Upper panel:* Orbit plotted in the (*x*, *y*), (*x*, *z*) and (*y*, *z*) planes. *Lower panel:* *R*<sub>max</sub> vs. *z*<sub>max</sub> diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the chaotic orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. A.2.** Example of an x-tube, i.e. 1 : 1 resonance between $\Omega_y$ and $\Omega_z$, also called (0 : 1 : -1). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the x-tube orbits identified with the visual method (see 4.2.2).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. A.3.** Example of a banana orbit (here, a (x, z) banana), i.e. 2 : 1 resonance between two orbital frequencies (here, $\Omega_z$ and $\Omega_x$, also called (2 : 0 : -1)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and anti-banana orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. A.2.** Example of an x-tube, i.e. 1 : 1 resonance between $\Omega_y$ and $\Omega_z$, also called (0 : 1 : -1). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the x-tube orbits identified with the visual method (see 4.2.2).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.3.** Example of a banana orbit (here, a (x, z) banana), i.e. 2 : 1 resonance between two orbital frequencies (here, $\Omega_z$ and $\Omega_x$, also called (2 : 0 : -1)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the banana and anti-banana orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. A.4.** Example of a fish orbit (here, a *(y, z)* fish), i.e. 3 : 2 resonance between two orbital frequencies (here, *Ω*<sub>*z*</sub> and *Ω*<sub>*y*</sub>, also called (0 : 3 : −2)). *Upper panel:* Orbit plotted in the *(x, y)*, *(x, z)* and *(y, z)* planes. *Lower panel:* *R*<sub>max</sub> vs. *z*<sub>max</sub> 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).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.5.** Example of a saucer orbit, i.e. 1 : 1 resonance between $\Omega_z$ and $\Omega_R$. *Upper panel:* Orbit plotted in the $(x, y)$, $(x, z)$, and $(y, z)$ planes. *Lower panel:* $R_{max}$ vs. $z_{max}$ diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the saucer orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. A.4.** Example of a fish orbit (here, a (y, z) fish), i.e. 3 : 2 resonance between two orbital frequencies (here, $\Omega_x$ and $\Omega_y$, also called (0 : 3 : -2)). Upper panel: Orbit plotted in the (x, y), (x, z) and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the fish and anti-fish orbits identified with the visual method (see 4.2.2).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. A.5.** Example of a saucer orbit, i.e. 1 : 1 resonance between $\Omega_x$ and $\Omega_y$. Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the saucer orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. A.6.** Example of a pretzel orbit, i.e. 4 : 3 resonance between two orbital frequencies (here, $\Omega_z$ and $\Omega_y$, also called (0 : 4 : -3)). **Upper panel:** Orbit plotted in the (x, y), (x, z) and (y, z) planes. **Lower panel:** $R_{max}$ vs. $z_{max}$ diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the pretzel and anti-pretzel orbits identified with the visual method (see 4.2.2).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. A.7.** Example of a 5:4 orbit, i.e. 5 : 4 resonance between two orbital frequencies (here $\Omega_z$, and $\Omega_y$, also called (0 : 5 : -4)). **Upper panel:** Orbit plotted in the (x, y), (x, z), and (y, z) planes. **Lower panel:** $R_{max}$ vs. $z_{max}$ diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the 5:4 orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. A.6.** Example of a pretzel orbit, i.e. 4 : 3 resonance between two orbital frequencies (here, $\Omega_z$ and $\Omega_y$, also called (0 : 4 : -3)). **Upper panel:** Orbit plotted in the (x, y), (x, z) and (y, z) planes. **Lower panel:** $R_{max}$ vs. $z_{max}$ diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the pretzel and anti-pretzel orbits identified with the visual method (see 4.2.2).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. A.7.** Example of a 5:4 orbit, i.e. 5 : 4 resonance between two orbital frequencies (here $\Omega_z$, and $\Omega_y$, also called (0 : 5 : -4)). **Upper panel:** Orbit plotted in the (x, y), (x, z), and (y, z) planes. **Lower panel:** $R_{max}$ vs. $z_{max}$ diagram (**left**) and Cartesian frequency map (**right**). The coloured markers correspond to the 5:4 orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. A.8.** Example of a 5:6 orbit, i.e. 5 : 6 resonance between two orbital frequencies (here, $\Omega_y$ and $\Omega_z$, also called (0 : 6 : -5)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:6 orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. A.8.** Example of a 5:6 orbit, i.e. 5 : 6 resonance between two orbital frequencies (here, $\Omega_y$ and $\Omega_z$, also called (0 : 6 : -5)). Upper panel: Orbit plotted in the (x, y), (x, z), and (y, z) planes. Lower panel: $R_{max}$ vs. $z_{max}$ diagram (left) and Cartesian frequency map (right). The coloured markers correspond to the 5:6 orbits identified with the visual method (see 4.2.2).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>
**Appendix B: Comparison of potentials**

To assess the sensitivity of the study outcomes to the selec-
tion 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 Rmax vs. zmax diagram pres-
ened in this paper (Fig. 5) with those obtained with an axisym-
metric potential (MWPotential14 from Bovy 2015) or by chang-
ing the NSD potential (Sormani et al. 2022 and Launhardt et al.
2002), the bar potential (Portail et al. 2017), or the bar angle.

<figure><img src="image_23.png" /><figcaption>**Fig. B.1.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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).</figcaption></figure>
**Appendix B: Comparison of potentials**

To assess the sensitivity of the study outcomes to the selec-
tion 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 Rmax vs. zmax diagram pres-
ened in this paper (Fig. 5) with those obtained with an axisym-
metric potential (MWPotential14 from Bovy 2015) or by chang-
ing the NSD potential (Sormani et al. 2022 and Launhardt et al.
2002), the bar potential (Portail et al. 2017), or the bar angle.

<figure><img src="image_22.png" /><figcaption>**Fig. B.1.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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).</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. B.2.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. B.2.** *R<sub>max</sub> vs. z<sub>max</sub> 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.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. B.3.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. B.4.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. B.3.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. B.4.** *R*<sub>max</sub> vs. *z*<sub>max</sub> 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.</figcaption></figure>