A ﬁrst glimpse at the line-of-sight structure of the Milky Way’s nuclear stellar disc

,


Introduction
The nuclear stellar disc (NSD) is a dense stellar structure at the centre of the Galaxy with a radius of ∼150 pc and a scale height of ∼40 pc (e.g., Launhardt et al. 2002;Gallego-Cano et al. 2020;Sormani et al. 2020Sormani et al. , 2022)).It partially overlaps with the central molecular zone (CMZ), a dense accumulation of gas in orbits around the Galactic centre (e.g., Henshaw et al. 2022), and constitutes a distinct structure from the larger Galactic bulge/bar and from the Milky Way's nuclear star cluster located at its centre (e.g., Schultheis et al. 2021;Nogueras-Lara et al. 2021c;Nogueras-Lara 2022).The NSD hosts a total stellar mass of ∼10 9 M (e.g., Launhardt et al. 2002;Nogueras-Lara et al. 2020a;Sormani et al. 2022), and is characterised by a predominantly old stellar population (∼90% of the stellar mass is older than 8 Gyr), but also presents a significant mass of stars (∼5% of the total stellar mass) with an age of ∼1 Gyr (Nogueras-Lara et al. 2020a).Moreover, there is evidence of current star formation that makes the NSD the most active starforming region in the Galaxy when averaged over volume (e.g., Mezger et al. 1996;Matsunaga et al. 2011;Nogueras-Lara et al. 2020a, 2022).
The NSD is located at only 8 kpc from Earth, which makes it a unique template to understand stellar nuclei and their role in the context of galaxy evolution.Its formation seems to be closely related to the Galactic bulge/bar, which funnelled gas towards the Galactic centre originating the formation of the NSD (e.g., Nogueras-Lara et al. 2020a).Simulations suggest that the NSD formation is indeed related to the bar, and can be even used to date the Galactic bar (e.g., Baba & Kawata 2020).Therefore, the NSD is related to the Galaxy as a whole and a detailed analysis of its structure and stellar population is crucial to understanding the origin and properties of the Galactic bulge/bar.However the high extinction and the extreme source crowding towards the Galactic centre challenge the analysis of the NSD stellar population (e.g., Nishiyama et al. 2008;Nogueras-Lara et al. 2021b).This means that the NSD structure is not clear yet given the complicated determination of the line-of-sight distance towards individual stars.
In this work we use red clump (RC) stars, red giant stars in their helium core burning sequence, whose intrinsic properties are well known (e.g., Girardi 2016) and kinematic data (Libralato et al. 2021) to compute the line-of-sight distance and extinction towards different regions in the NSD.We obtained the first direct estimate of the NSD depth along the line of sight and conclude that our results are compatible with an axisymmetric morphology of the NSD, with relatively low extinction within it, although based on the present results we cannot rule out that the NSD is barred.

Proper motions
We used the publicly available NSD proper motion catalogue obtained with the Wide-Field Camera 3 (WFC3/IR, filter F153M) by Libralato et al. (2021).The catalogue contains absolute proper motions for more than 800 000 stars calibrated with the Gaia Data Release 2 (Gaia Collaboration 2016, 2018).The proper motions were calculated using point spread function photometry on images from two different epochs (2012 and 2015).Figure 1 shows the region covered by the catalogue.For our analysis we removed stars with high uncertainties, and only used proper motions with ∆µ RA , ∆µ Dec < 1.5 mas yr −1 .We also applied a rotation to transform the equatorial coordinates into Galactic coordinates to make the analysis of the NSD kinematics easier.

Photometry
We cross-correlated the H and K s data from the GALAC-TICNUCLEUS survey (Nogueras-Lara et al. 2018, 2019) with the proper motion catalogue to build a colour-magnitude diagram (CMD) and identify RC stars for the subsequent analysis.The GALACTICNUCLEUS survey was specially designed to observe the Galactic centre, and contains accurate point spread function photometry for more than three million stars in the NSD and the innermost Galactic bulge/bar.The photometric uncertainties are below 0.05 mag at H ∼ 19 mag and K s ∼ 18 mag.The zero point systematic uncertainty is below 0.04 mag.
Figure 2 shows the resulting CMD K s vs. H − K s .To select the RC stars from the NSD, we firstly removed foreground stars, which mainly belong to the Galactic disc and the Galactic bulge/bar, by applying a colour cut H − K s = 1.3 mag.This accounts for the significantly different extinction between these Galactic components and the NSD (e.g., Nogueras-Lara et al. 2021a).We also assumed a colour cut H − K s = 2.5 mag to avoid regions whose completeness is too low.The red parallelogram in Fig. 2 shows the resulting selection of RC stars.

Proper motion distribution
Recent work (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022) computed the proper motion distribution of the NSD distinguishing between its parallel and perpendicular components to the Galactic plane (µ l and µ b , respectively).They found that the µ l distribution of the NSD presents three stellar components: two of them due to the rotation of the NSD, and the third probably due to the contamination from the Galactic bulge.On the other hand, the µ b distribution shows two components centred around 0 mas yr −1 that are probably due to the combined contributions from the NSD and the Galactic bulge/bar.
Here we use a deeper data set (Libralato et al. 2021), which allows us to reach the RC feature and use these stars to trace the stellar kinematics, instead of using a mixture of brighter stars with different stellar types (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022).We obtained the distribution of proper motions and applied the SCIKIT-LEARN Python function GaussianMixture (GMM, Pedregosa et al. 2011) to compute the probability density function that originates the underlying data distribution based on Gaussian models.We tried up to three Gaussians for the µ l distribution, and up to two for the µ b , in agreement with previous work (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022).We computed the Bayesian information criterion (BIC, Schwarz 1978) and the Akaike information criterion (AIC, Akaike 1974) to choose the model that best represents the observed distribution.We obtained that the µ l distribution is best represented by a three-Gaussian model, whereas the µ b distribution is best represented by a two-Gaussian model.Figure 3 and Table 1 show the obtained results.We estimated the uncertainties of the GMM parameters resorting to Monte Carlo simulations and generating 1000 samples of synthetic data assuming that each data point can vary following a Gaussian distribution with a standard deviation equal to its uncertainty.The three-Gaussian distribution found for µ l agrees with previous work (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022) and indicates the presence of two stellar populations moving eastwards and westwards, which we interpret as NSD rotation.Our results are computed using absolute proper motions, which explains why the central peak is not close to zero.We measured a difference of ∼3 mas yr −1 between the mean values of the µ l components moving eastwards and westwards, which is somewhat smaller that the difference obtained in previous work (∼4 mas yr −1 , in Shahzamanian et al. 2022).Moreover, the Gaussian model corresponding to the stellar population moving westwards is smaller (∼50% lower) than the model corresponding to the stellar population moving eastwards, whereas they should be similar in size according to previous work (Shahzamanian et al. 2022).We believe that this is because the stellar population moving westwards corresponds to the farthest edge of the NSD (Shahzamanian et al. 2022), and then the extinction towards these stars is higher and their completeness (also affected by the extreme source crowding), is lower than for the eastward-moving stars.The same effect can also cause the difference between the central peak and the extreme ones to be smaller than measured in previous work.On the other hand, the differences between our results and the work by Shahzamanian et al. (2022) can also be explained because their data set contains the nuclear star cluster, which is not present in our data, whose proper motion distribution is different from the NSD (Nogueras-Lara 2022).Moreover, the discrepancies can also be due to the different techniques used to analyse the data.The work by Shahzamanian et al. (2022) applied a direct fit to the data instead of the GMM modelling that we use in this work.Our approach has the advantage of being independent of the data binning, and informs about the underlying Gaussian model that originates the observed data distribution.
On the other hand, the two-Gaussian model found for the µ b distribution agrees well with previous work (e.g., Shahzamanian et al. 2022).We found that the wider Gaussian, probably due to stars from the Galactic bulge/bar, accounts for ∼15% of the stars in the distribution.This is smaller than the ∼50% found by Shahzamanian et al. (2022).We believe that this is due to the smaller region closer to the Galactic plane that we consider in this work, which reduces the relative fraction of contaminant stars from the Galactic bulge/bar.Furthermore, our results agree with the expected contamination (∼20%) for the innermost fields of the NSD from Galactic bulge/bar stars estimated by Sormani et al. (2022) for stars with a colour cut of H − K s > 1.3 mag (see their Table 2).The fact that we obtain a contribution of ∼5% when considering the central Gaussian of the µ l distribution is probably due to the more challenging identification of Galactic bulge/bar stars in this proper motion space in comparison to the µ b distribution.Shahzamanian et al. (2022) also measured a lower contribution of Galactic bulge/bar stars in the µ l distribution (∼38%) in comparison with the µ b distribution (∼55%).

Testing the asymmetry of the stellar kinematics
It is well known that the gas distribution in the innermost 300 pc of the Galaxy is strongly asymmetric (e.g., Bally et al. 1988;Sormani et al. 2018).We tested whether the stellar kinematics is also different when considering different regions of the NSD in the following way.We divided our proper motion data set into two different regions (1+2 and 3+4 in Fig. 1), corresponding to different longitudes of the NSD. Figure 4 shows the proper motion distribution for each of the regions.We conclude that the µ l and µ b distributions are similar for both regions of the NSD, so there is no remarkable difference between the stellar kinematics depending on the observed regions.

Stellar sample
Using the proper motion distribution derived in Sect.3.1, we obtained a clean sample of RC stars belonging to the eastwardand westward-moving stars from the NSD; they trace the stellar population from the closest edge (eastward-moving stars) and the farthest one (westward-moving stars) of the NSD.For this, we firstly removed a significant fraction of the contaminant stars from the Galactic bulge/bar by choosing stars with L8, page 3 of 7  |µ b | < 5 mas yr −1 (green shaded region in Fig. 3), which reduces their contamination below 10%.We then computed the probability of a given star to be part of each of the three-Gaussian models that originate the µ l distribution.We defined two regions in the µ l distribution that are dominated (i.e., ∼80 % of the stars) by stars moving eastwards (grey shaded area in Fig. 3), and by stars moving westwards (red shaded area in Fig. 3).These regions are not symmetric and have different sizes because the corresponding underlying Gaussian distributions are also not symmetric due to the lower completeness of the stellar sample moving westwards (see Sect. 3.1).To analyse the behaviour of the NSD population for different lines of sight, we divided the stellar samples into four different regions, each with a similar number of stars, as shown in Fig. 1.

Extinction
We computed extinction to de-redden the stars in each region using a star-by-star approach applying the equation (e.g., Nogueras-Lara et al. 2021c,d) where A K s is the K s extinction, H and K s indicate the photometric data, (H − K s ) 0 = 0.10 ± 0.01 is the intrinsic colour (Nogueras-Lara et al. 2021b), and A H /A K s = 1.84 ± 0.03 is the extinction curve (Nogueras-Lara et al. 2020b).Table 2 shows the mean extinction values obtained by applying a 3σ clipping algorithm to remove outliers.The uncertainties were obtained quadratically propagating the uncertainties of the quantities involved in the calculation.Moreover, we also computed the relative difference in extinction between the stellar populations moving eastwards and westwards.In this case the uncertainties refer to the standard error of the distribution and to the uncertainty associated with the used extinction curve.Other systematic uncertainties from the quantities involved in Eq. ( 1) do not affect the difference between the extinctions towards each of the stellar populations.We obtained that the extinction is larger for the westwardmoving stellar population, which is in agreement with this population belonging to the farthest edge of the NSD.Comparing the extinctions between the eastward-and the westward-moving stars, we found that <10% of the total extinction occurs within the NSD.This suggests that the majority of gas and dust causing the extinction in the NSD is not inside it, which may indicate that most of the gas in the CMZ surrounds the NSD (e.g., Henshaw et al. 2022).This can be related with the inside-out formation of nuclear stellar discs suggested by Bittner et al. (2020) based on observations of external galaxies.In this picture, the CMZ would be the star-forming outer edge of the NSD.
An alternative explanation for the relatively low extinction within the NSD might be an extremely clumpy distribution of the dust and gas inside the NSD that would leave the majority of stars from the farthest edge of the NSD relatively unobscured.In this way, some regions of the NSD could contain a large quantity of dust and gas.To check this scenario, we produced extinction maps for the closest and the farthest NSD edges using the analysed RC stars.We followed the process described in Nogueras-Lara et al. (2021b), and chose a pixel size of 10 for the extinction maps.We computed an extinction value for each pixel using the five closest RC stars in a radius of 15 (if fewer than five stars are present within the reference radius for a given pixel, we did not compute an associated value).Figure 5 shows the obtained extinction maps.We also computed the difference between the two extinction maps (Fig. 5c).The main extinction structures appearing in both extinction layers are not present in the difference map, and thus affect the two stellar components in the same way, probably being due to the presence of dust and gas along the line of sight towards the Galactic centre, and likely belonging to the closest edge of the CMZ.We computed the mean value of the difference map and its standard deviation, 0.18 ± 0.16 mag, and obtained that less than 0.5% of the pixels have an extinction >3σ away from the mean (corresponding to regions with a potential clumpy distribution of dust and gas), and thus the extinction map is quite homogeneous.Even for the regions 3σ away from the mean value, the extinction within the NSD is only 30% of the total extinction along the line of sight.
L8, page 4 of 7 Notes.Region indicates the corresponding field indicated in Fig. 1.A i indicates the absolute extinction for each of the i stellar populations.∆A Ks is the relative extinction between the analysed stellar populations.∆K s corresponds to the de-reddened magnitude difference between the stellar populations.∆d indicates the relative distance between the stellar populations.

Relative distance
We used the previously computed extinctions to build dereddened K s luminosity functions for the RC stars of the eastward-and westward-moving stellar populations for each of the four NSD regions analysed.Given the stellar population present in the NSD, we expected to find two peaks in the RC bump.The brightest peak, around K s ∼ 13 mag, is caused by a predominantly old stellar population (more than 80% of the stellar mass is older than 8 Gyr).The secondary peak is around ∼0.5 mag fainter than the main peak, and it is produced by RC stars with an age of around 1 Gyr (Nogueras-Lara et al. 2021b).
We fitted the underlying distributions for each region and stellar population with different kinematics using a Gaussian model.
To account for the possible influence of the double peak in the luminosity functions, we computed the maximum (max) of each distribution and restricted the analysis to stars with K s ∈ [max − 0.55, max + 0.45] mag.In this way we guarantee that the fit is consistent and that the possible influence of the secondary peak is similar for all the regions and stellar groups with different kinematics.Figure 6 shows the results.We obtained that the mean magnitude of the RC bump is always brighter for the eastward-moving stellar population in comparison with the westward-moving one.Our results agree with a different distance between eastward-and westward-moving stars, as suggested by Shahzamanian et al. (2022).
We computed the relative distance between the eastward-and the westward-moving stellar populations for each of the regions, by applying the equation where r d is the ratio of the distances towards the two stellar populations, ∆µ is the difference in their distance modulus, and ∆K s is their de-reddened magnitude difference.To compute the distance difference between the two stellar populations, we assumed a Galactic centre distance of 8.1 kpc (GRAVITY Collaboration 2018), and considered that the closest edge of the NSD is 150 pc closer than that (e.g., Launhardt et al. 2002;Gallego-Cano et al. 2020;Sormani et al. 2020Sormani et al. , 2022)).In this way we placed the eastward-moving stellar population at a distance of 7.95 kpc and used the ratio r d to estimate the relative distance between the components.Table 2 shows the obtained results.The associated uncertainties were computed considering the uncertainty of the relative distance modulus (∆µ), which includes the uncertainty of the Gaussian fits, and also the relative uncertainty due to the de-reddening process.Given that we are interested in the relative distance between the stellar populations, the systematic uncertainties related to the photometric zero point, the absolute distance used to scale the values, the extinction curve, and the intrinsic colour do not affect the computed uncertainty.
We obtained that the distance between the two stellar populations is similar and compatible within the uncertainties for all the analysed regions.Given that the 2D projected radius of the NSD was estimated to be ∼150 pc (e.g., Launhardt et al. 2002;Gallego-Cano et al. 2020;Sormani et al. 2020Sormani et al. , 2022)), our results suggest a similar structure along and across the line of sight.This, in combination with the symmetric proper motion distribution in the analysed region (see Fig. 4) and the rotation pattern of the NSD obtained from line-of-sight velocities using APOGEE data (Schönrich et al. 2015), is compatible with an axisymmetric geometry of the NSD, although a nuclear bar as found in external galaxies (Bittner et al. 2021) cannot be excluded.We computed a mean relative distance value between the two stellar populations averaging over the results for each of the four analysed regions.We obtained d = 330 ± 50 pc, where the uncertainty corresponds to the standard deviation of the distribution.This value is compatible with the previously estimated NSD radius of about 150 pc (e.g., Launhardt et al. 2002;Schödel et al. 2014;Gallego-Cano et al. 2020;Sormani et al. 2020Sormani et al. , 2022)).Nevertheless, the obtained value is instead a lower limit of the NSD diameter because it corresponds with the relative distance between the stellar populations considered in this work, and might be lower than the actual diameter because our sample of RC stars is probably not complete at its faint end.In any case, the depth along the line of sight is comparable to the extension across the line of sight, consistent with an axisymmetric structure of the NSD.On the other hand, this analysis does not imply that the NSD is necessarily a ring because we did not consider the transition region where stars from the eastward-and westward-moving populations are mixed (white space between the red and grey shaded regions in Fig. 3), and might correspond to intermediate distances, building a disc-like geometry or even a nuclear bar.

Conclusion
In this Letter we studied the kinematics of RC stars from the NSD.We analysed their proper motion components parallel (µ l ) and perpendicular (µ b ) to the Galactic plane by applying a GMM modelling.We obtained that the µ l component is best described by a three-Gaussian model, which we interpret as a result of the rotation of the NSD in combination with some residual contamination from Galactic bulge/bar stars.The µ b distribution is described by a two-Gaussian model that accounts for the NSD stars and also the contaminant stars from the Galactic bulge/bar.Our results are consistent with previous work (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022), with the exception of a somewhat smaller contribution from Galactic bulge/bar stars.We believe that this is due to the region used in this work, smaller and closer to the Galactic plane, which reduces the contamination from the Galactic bulge/bar.Moreover, we divided the data set into two regions with different Galactic longitudes and checked that the proper motion distributions are similar, pointing towards a symmetric NSD structure.
We also created a clean sample of RC stars belonging to the eastward-and westward-moving stellar populations, and computed their relative distance and extinction towards four subregions covering different Galactic longitudes.We obtained that the relative distance between the two stellar populations is constant within the uncertainties.These results, in combination with previous measurements of the NSD extension across the line of sight (e.g., Launhardt et al. 2002;Gallego-Cano et al. 2020) and the symmetric distribution of the proper motions, are consistent with the NSD being axisymmetric, although the presence of a nuclear bar cannot be ruled out with the current measurements.Moreover, we find that the extinction within the NSD is relatively small in comparison with the overall extinction of the NSD stars.Our findings suggest that the most of the extinction in the Galactic centre is due to the presence of the CMZ that surrounds the NSD, being the amount of gas and dust inside the NSD relatively small.

Fig. 1 .
Fig. 1.GALACTICNUCLEUS false colour image of the central region of the NSD.The white rectangles indicate the regions covered by the proper motion catalogue.The white dashed circle shows the effective radius of the nuclear star cluster (∼5 pc, e.g., Gallego-Cano et al. 2020) with Sagittarius A* at its centre.The white dashed lines indicate the four regions used for the distance and extinction analysis in Sect. 4.

Fig. 2 .
Fig. 2. CMD K s vs. H −K s for common stars between the proper motion catalogue and the GALACTINUCLEUS survey.The red parallelogram shows the selection of RC stars.

Fig. 3 .
Fig.3.GMM analysis of the proper motion distributions for all the RC stars in the analysed region.The left and right panels show the distribution of the proper motion components parallel (µ l ) and perpendicular (µ b ) to the Galactic plane, respectively.The blue histograms correspond to the data, whereas the orange lines are the total models (solid lines), and each of their Gaussian components (dashed lines).The grey and red shaded areas in the left panel, and the green shaded area in the right panel correspond to the regions analysed in Sect. 4.

Fig. 4 .
Fig. 4. Proper motion distributions of east and west NSD regions.The orange dots and histograms correspond to the proper motions distribution of the region 1+2 in Fig. 1, whereas the blue dots and histograms correspond to region 3+4 in Fig. 1.

Fig. 5 .
Fig. 5. Extinction maps towards the analysed regions (panels a and b, built using eastward-and westward-moving RC stars, respectively), and their difference (panel c).The black star indicates the position of Sagittarius A*.The blue arrows in panels a and b correspond to regions with similar extinction structures in both maps.The white pixels indicate that there is not any extinction value associated, due to the low number of reference stars.

Fig. 6 .
Fig. 6.De-reddened RC K s luminosity function for different NSD regions.The histograms, Gaussian fits, and numbers in each panel refer to the eastward-moving stellar population (in grey) and the westward-moving stars (in red).The numbers in each panel indicate the associated region in Fig. 1.The K s e and K s w values indicate the mean value of the Gaussian fits and its associated uncertainty for the eastward-and westward-moving stars.

Table 1 .
Results from the GMM analysis of the proper motion distribution.

Table 2 .
Results from the analysis of the de-reddened RC luminosity functions.