Evidence of an age gradient along the line of sight in the nuclear stellar disc of the Milky Way

The nuclear stellar disc (NSD) is a flat dense stellar structure at the heart of the Milky Way. Recent work has shown that analogous structures are common in the nuclei of external spiral galaxies, where there is evidence of an age gradient that indicates that they form inside-out. However, the characterisation of the age of the NSD stellar population along the line of sight is still missing due to its extreme source crowding and the high interstellar extinction towards the Galactic centre. We aim to characterise the age of the stellar population at different average Galactocentric NSD radii to investigate for the first time the presence of an age gradient along the line of sight. We selected two groups of stars at different NSD radii via their different extinction and proper motion distribution. We analysed their stellar population by fitting their de-reddened $K_s$ luminosity functions with a linear combination of theoretical models. We find significant differences in the stellar population at different NSD radii, indicating the presence of an age gradient along the line of sight. Our sample from the closest edge of the NSD contains a significant fraction ($\sim40 \%$ of its total stellar mass) of intermediate-age stars (2-7 Gyr) that is not present in the sample from stars deeper inside the NSD, in which $\sim90 \%$ of the stellar mass is older than 7 Gyr. Our results suggest that the NSD age distribution is similar to the one found in external galaxies and they imply that bar-driven processes observed in external galaxies are similarly at play in the Milky Way.


Introduction
The nuclear stellar disc (NSD) is a flat dense stellar structure that roughly outlines the Galactic centre 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)).Its study is hampered by the extreme source crowding and the high extinction towards the innermost regions of the Galaxy (e.g., Nishiyama et al. 2006;Nogueras-Lara et al. 2021b).The analysis of its structure and stellar population is therefore mainly limited to near-infrared high-resolution photometry, while spectroscopy is restricted to very small regions of high interest (e.g., Lohr et al. 2018;Clark et al. 2018).
A photometric analysis of the innermost ∼1600 pc 2 of the NSD using the GALACTICNUCLEUS catalogue (a highangular, ∼0.2 , resolution photometric survey in the nearinfrared, Nogueras-Lara et al. 2018a, 2019, 2020a), determined that the NSD is mainly dominated by old stars (more than 80% of its total stellar mass was formed more than 8 Gyr ago).A similar analysis of Sgr B1, a region of intense HII emission located ∼100 pc away from the centre of the NSD (e.g., Simpson et al. 2021), showed evidence of an intermediate-age stellar population (∼40% of the stellar mass was formed ∼2−7 Gyr ago) in addition to the old stars (∼40% of the stellar mass is older than 7 Gyr, Nogueras-Lara et al. 2022).This is consistent with the inside-out formation scenario proposed for the NSD (based on recent studies on nearby spiral galaxies, Gadotti et al. 2020;Bittner et al. 2020) that would originate a radial age gradient.In this picture, the formation of the NSD is related to the Galactic bar, which funnels gas from the Galactic disc towards the innermost regions of the Galaxy (e.g., Sormani & Barnes 2019), producing an accumulation of gas and dust that grows the NSD from inside-out.
In this Letter, we follow up on previous work that studied stellar populations at different lines of sight (Nogueras-Lara et al. 2022) and investigate the age distribution of stellar populations located at different NSD radii along the line of sight.We aim to assess the presence of an age gradient that would support the inside-out formation scenario proposed for the NSD.

Photometry
We used H and K s photometry from the GALACTICNU-CLEUS survey (Nogueras-Lara et al. 2018a, 2019).This is a high-angular (∼0.2 ) resolution near-infrared catalogue specially designed to observe the NSD.It contains accurate point spread function photometry for more than three million sources.The zero point systematic uncertainty is below 0.04 mag in all bands, whereas the statistical uncertainties are below 0.05 mag at H ∼ 19 mag and K s ∼ 18 mag.
We corrected potentially saturated sources in K s (K s < 11.5 mag, e.g., Nogueras-Lara et al. 2019) using K s data from the SIRIUS IRSF survey (Nishiyama et al. 2008).We also included saturated sources which were not detected in the GALACTICNUCLEUS catalogue.
Figure 1 indicates the target region, which was chosen to overlap with a high-precision proper motion catalogue of the Galactic centre (Libralato et al. 2021).We excluded regions close to the nuclear star cluster because its different star formation history (e.g., Schödel et al. 2020;Nogueras-Lara et al. 2021c) could affect our analysis.

Proper motions
We used the Galactic centre proper motion catalogue obtained by Libralato et al. (2021) using the Wide-Field Camera 3 (WFC3/IR, filter F153M) at the HST.This catalogue contains absolute high-precision proper motions for more than 800 000 stars calibrated with Gaia DR2 (Gaia Collaboration 2016, 2018).

Colour-magnitude diagram and targets selection
Figure 2 shows the colour magnitude diagram (CMD) K s versus H − K s of the target region.To remove foreground stars from the Galactic disc and bulge/bar, we applied a colour cut H − K s 1.3 mag (e.g., Nogueras-Lara et al. 2021a,b).
Recent results on the NSD structure indicate that the extinction within the NSD correlates with the position of the stars along the line of sight (Nogueras-Lara 2022a).Hence, stars belonging to the closest edge of the NSD suffer a lower reddening than stars from its farthest edge.It is then possible to statistically distinguish between stars located at different NSD radii along the line of sight by applying a colour cut in the CMD K s versus H − K s .We selected two target groups corresponding to stars from the closest edge of the NSD and stars deeper inside the NSD (hereafter NSD outer and inner regions, respectively) by applying the colour cuts indicated in Fig. 2.
To check our target selection and assess whether the target stellar groups are dominated by stars located at different NSD radii, we cross-correlated the photometry from the GALAC-TICNUCLEUS survey (Nogueras-Lara et al. 2018a, 2019) with the proper motion catalogue by Libralato et al. (2021).Figure 3 shows the proper motion distribution of stars from each of the groups.We clearly distinguish a different distribution of the proper motion component parallel to the Galactic plane (µ l ) for each of the groups.We computed a mean value of µ l = −4.66 ± 0.02 mas yr −1 for the group of stars with the bluest L10, page 2 of 6  colour (i.e. the lowest extinction), where the uncertainty corresponds to the standard error of the distribution.This value is in agreement with the rotation of stars from the closest edge of the NSD, µ l = −4.70 ± 0.03 mas yr −1 , obtained in Nogueras-Lara (2022a).On the other hand, we obtained a mean value of µ l = −5.77± 0.03 mas yr −1 for the group of stars with the reddest colour (i.e. the highest extinction).This value is lower than the one computed in previous work for stars from the farthest edge of the NSD (µ l = −7.65 ± 0.03 mas yr −1 , see Table 1 in Nogueras-Lara 2022a), but significantly larger than the result for stars from the closest edge.This indicates that the stars in this group mainly belong to the inner regions of the NSD, being placed at more internal radii than stars from the farthest edge of the NSD and the outer group previously defined.
We found that the proper motion distribution of the component perpendicular to the Galactic plane is µ b ∼ 0 mas yr −1 , which is in agreement with previous work (Shahzamanian et al. 2022;Martínez-Arranz et al. 2022;Nogueras-Lara 2022a).Small differences in the shape of the µ b distribution of the target groups might be caused by a different contribution of contaminant stars from the Galactic bulge/bar which can depend on the applied colour cut for the sample selection.

Stellar population analysis
To analyse the stellar population in each of the target groups, we created K s luminosity functions and fitted them with a linear combination of theoretical models following the method in Nogueras-Lara et al. (2020a), Schödel et al. (2020), andNogueras-Lara et al. (2022).

Reddening correction
We de-reddened the stars in each of the target groups by creating dedicated extinction maps following the approach described in Nogueras-Lara et al. (2021b) and using the extinction curve A H /A K s = 1.84 ± 0.03 (Nogueras-Lara et al. 2020b).To build the extinction maps, we used red clump stars (giant stars in their helium core-burning sequence, e.g., Girardi 2016) and other red giant stars with a similar brightness (dashed boxes in Fig. 2), as reference stars, and assumed that they have an intrinsic colour of (H−K s ) 0 = 0.10±0.01(Nogueras-Lara et al. 2021b).We defined a pixel size of ∼10 and used the five closest reference stars within a maximum radius of 15 from the centre of each pixel to compute the extinction value applying an inverse distance weight method (for details see Appendix A.2 of Nogueras-Lara et al. 2021b).We only computed an extinction value for a given pixel if at least five reference stars were found within the maximum radius.
Figure 4 shows the obtained extinction maps.We computed mean extinction values of A K s ∼ 1.7 mag and A K s ∼ 2.4 mag for the extinction maps corresponding to stars from the NSD outer and inner regions, respectively.We applied the extinction maps to de-redden the photometry of each of the target groups.

Luminosity function
We created a K s luminosity function using the de-reddened photometry for each of the target groups.We only used stars within the defined colour cut in Fig. 2 for the group of stars from the NSD outer region, whereas we also included stars that were detected in K s but not in H for the group of stars from the inner region of the NSD for the sake of completeness.Not detecting these stars in H indicates that they are affected by a larger extinction and then belong to the target group of stars.
To avoid over-de-reddened stars, we excluded stars with a de-reddened H − K s colour more than 2σ bluer than the mean value of the de-reddened distribution of the red clump features (e.g., Nogueras-Lara et al. 2020a, 2022).To build the luminosity functions, we chose the bin width that maximises the Freedman-Diaconis (Freedman & Diaconis 1981) and Sturges (Sturges 1926) estimators, using the 'auto' option in the python function numpy.histogram(Harris et al. 2020).We assumed Poisson uncertainties (i.e. the square root of the number of stars) for each luminosity function bin.L10, page 3 of 6 The faint end of the luminosity functions is significantly affected by incompleteness due to the extreme source crowding in the NSD (e.g., Nogueras-Lara et al. 2020a, 2022).In this way, we computed a crowding completeness solution determining the critical distance at which a star can be detected around a brighter star, following the method described in Eisenhauer et al. (1998).We divided the analysed field into small subregions of 2 × 1.4 and computed a completeness solution for each of them.We then combined the results and computed their mean and standard deviation to obtain the final completeness and its uncertainty, as explained in Nogueras-Lara et al. (2020b).We used our completeness solution to correct the luminosity functions, setting a lower limit of 70% of data completeness.On the other hand, due to possible saturation and incompleteness of the bright end of the luminosity function, even after using the SIRIUS IRSF catalogue to correct saturated stars (see Sect. 2.1), we restricted the bright end of the luminosity function to K s = 8.5 mag for the subsequent analysis (e.g., Nogueras-Lara et al. 2022).

Stellar populations
Luminosity functions contain fundamental information about the star formation history.In this way, the presence of particular features such as the asymptotic giant branch bump, the red clump, or the red giant branch bump, and their relative weights, allows us to characterise their stellar population (e.g., Nogueras-Lara et al. 2018b, 2020a;Schödel et al. 2020).
We fitted the obtained K s luminosity functions with a linear combination of theoretical models to determine the stellar populations present in each target group, as explained in Nogueras-Lara et al. (2020a, 2022).We used 14 Parsec models 1 (Bressan et al. 2012;Chen et al. 2014Chen et al. , 2015;;Tang et al. 2014;Marigo et al. 2017;Pastorelli et al. 2019Pastorelli et al. , 2020)), choosing ages that homogeneously sample the possible stellar populations in the analysed luminosity functions (14, 11, 8, 6, 3, 1.5, 0.6, 0.4, 0.2, 0.1, 0.04, 0.02, 0.01, and 0.005 Gyr).The age-spacing between the models was selected to account for the differences that appear in the luminosity function for stellar populations with different ages.In this way, we increased the number of models towards younger ages, where the luminosity function is more affected by age variations (see Sect.Methods in Nogueras-Lara et al. 2022).
We assumed a Kroupa initial mass function corrected for unresolved binaries (Kroupa et al. 2013), and twice solar metallicity (Z = 0.03) for our models, in agreement with previous results for the NSD (e.g., Schultheis et al. 2019Schultheis et al. , 2021;;Nogueras-Lara et al. 2020a, 2022).Our fitting algorithm included a parameter to consider the distance towards the target stellar populations (including potential offsets introduced in the de-reddening process), and also a Gaussian smoothing factor accounting for possible distance differences between stars and also some residual uncorrected differential extinction.Figure 5 shows the best fits obtained by minimising a χ 2 for each of the luminosity functions, where we combined models with similar ages to minimise the degeneracy.
We resorted to Monte Carlo simulations to determine the contribution of each model to the luminosity function and estimate the uncertainties, as explained in Nogueras-Lara et al. (2022).We created 1000 Monte Carlo samples by computing the number of stars per bin in each luminosity function by randomly varying the original value assuming a Gaussian distribution with a standard deviation equal to each bin uncertainty.We then fit-1 Generated by CMD 3.6 (http://stev.oapd.inaf.it/cmd).ted each Monte Carlo sample following our method and average over the results to compute a mean value and its standard deviation.
We repeated the process using also MIST models (Paxton et al. 2013;Dotter 2016;Choi et al. 2016) with similar ages and a Salpeter initial mass function, to account for possible systematic effects caused by the choice of different stellar evolutionary tracks (e.g., Nogueras-Lara et al. 2022).We obtained small differences between the analyses carried out with Parsec and MIST models, as shown in Fig. 6.The final results were computed by combining the values obtained when using both models (Fig. 6).The uncertainties were calculated quadratically propagating the ones obtained for each age bin with Parsec and MIST models.
From Fig. 6, we conclude that the stellar populations in each of the target regions are significantly different.While most stars across the NSD have ages above 7 Gyr, we find that in the outer region of the NSD there is a factor 3 more stars with ages between 2 and 7 Gyr than in its inner region, strongly suggesting the presence of an age gradient.

Systematic uncertainties
To assess our results, we studied potential sources of systematic uncertainties: 1. Luminosity function bin width.We repeated the described process for each of the target populations assuming half and double the previously computed bin width.We did not observe any significant difference within the uncertainties.
2. Completeness solution and faint end of the luminosity function.We created two additional completeness solutions considering the completeness values per magnitude bin, subtracting and adding the corresponding uncertainty.We corrected the luminosity functions applying these two completeness solutions and repeated the fitting process.We slightly changed the faint end of the luminosity function to avoid problems related to incompleteness due to sensitivity.We did not observe any significant variation within the obtained uncertainties.We also repeated the process assuming a completeness limit of 75% instead of 70% without any significant change.
3. Bright end of the luminosity function.We set a bright limit of K s = 8.5 mag (de-reddened) for our study.Nevertheless, we also repeated our analysis including stars up to K s = 8 mag and conclude that it does not affect our results in any significant way.Only a somewhat larger contribution ( 2% variations) is observed for some young age bins (<2 Gyr) in the sample from the NSD outer region.
4. Metallicity of the stellar population.We assumed twice solar metallicity for our analysis in agreement with the super solar metallicities found for NSD stars (e.g., Schultheis et al. 2021;Nogueras-Lara 2022b;Feldmeier-Krause 2022).Nevertheless, we also explored solar metallicity and 1.5 solar metallicity using Parsec models.In both cases we determined that the old contribution is shifted towards the age bin 2-7 Gyr for the NSD outer region.On the other hand, there is not any variation in the results obtained for the NSD inner region.This result still agrees with the presence of different stellar populations for different NSD radii and is also compatible with an inside-out formation of the NSD.However, we determined that the reduced χ 2 is larger for both target populations when considering solar or 1.5 solar metallicity, indicating that assuming twice solar metallicity gives a better result.
5. Influence of the extinction maps.We repeated the process for each target group assuming different parameters to build the extinction map (seven reference stars, and a maximum radius of 20 to select reference stars).We did not observe any significant change within the uncertainties.
6. Width of the colour cut for target selection.We selected stars in the NSD outer region by applying a colour cut width ∆(H − K s ) = 0.3 mag, whereas the width of the colour cut was larger for the NSD inner region (Fig. 2).Our choice was motivated by the lower completeness when considering the stars from the inner regions of the NSD.To assess our results, we also tested the effect of assuming a colour cut with a similar width (H − K s = 0.3 mag) for stars in the NSD inner region.This implies that the maximum variation of the differential extinction is the same for both stellar samples.We determined that the contribution of stars older than 7 Gyr is ∼80%, whereas the contribution from the age bin 2−7 Gyr is ∼10%.This is probably because this stellar sample contains fewer stars from the innermost regions of the NSD than the original one.Therefore, our results are compatible with an age gradient in which stars in the colour range H − K s ∼ 1.8−2.1 mag correspond to an intermediate case in between the previously analysed ones (Fig. 6.) 7. Initial mass function.Recent results on the known young clusters in the NSD point towards a top-heavy initial mass function (e.g., Hosek et al. 2019).Nevertheless, our stellar samples are dominated by giant stars with a mass of 1 M and therefore assuming a top-heavy initial mass function does not have any observable effect on our results (e.g., Schödel et al. 2020).
8. Contamination from the Galactic bulge/bar in the target groups.We used a colour cut H − K s 1.3 mag to remove foreground stars from the NSD sample.However, residual contamination from the innermost regions of the Galactic bulge/bar could remain in our target groups ( 20% of the stars in the sample, as estimated by Sormani et al. 2022).Given that the stellar population from the innermost bulge/bar is mainly old and metal rich (e.g., Clarkson et al. 2011;Nogueras-Lara et al. 2018b;Renzini et al. 2018), its presence would not affect the detected radial age gradient, and could only increase the number of stars in the oldest age bin, without affecting our conclusions.9. Contamination between the NSD outer and inner regions.We chose our target stellar groups statistically distinguishing between stars located at different NSD radii using their different extinctions.Nevertheless, it is not possible to completely separate stars from the outer and the inner regions of the NSD given their proximity (the radius of the NSD is ∼150 pc, Nogueras-Lara 2022a) and their small kinematic difference (∼3−4 mas yr −1 between the NSD edges, e.g., Shahzamanian et al. 2022;Nogueras-Lara 2022a).To estimate the contamination between the target groups, we chose stars with well-defined proper motions (i.e.uncertainty below 0.5 mas yr −1 ).For the NSD outer region, we assumed that stars with µ l < −5.77 + 3 × 0.02 are contaminant stars from the inner region, where the reference value corresponds to the mean µ l for the inner region plus 3σ (see Sect. 3).Analogously, for the NSD inner region, we estimated the fraction of stars with µ l > −4.66−3 × 0.03 and considered that they likely belong to the outer region.We ended up with a contamination of 30% for each of the target groups from the other one.Nevertheless, this is probably an overestimation given the low number of stars with proper motion uncertainties below 0.5 mas yr −1 (∼15% of the stars with proper motions).This is particularly clear for the inner NSD, in which having 30% of stars from the outer NSD would imply a much larger stellar population in the 2-7 Gyr age bin than the detected one.Actually, an alternative way to estimate the contamination from the outer NSD in the inner NSD sample is to assume that the detected stellar contribution from stars in the age bin 2-7 Gyr is completely due to a stellar population as the one measured for the outer NSD.In this way we conclude that ∼5% of contamination from the outer NSD stellar population can account for the estimated ∼1% uncertainty in the age bin 2-7 Gyr (middle panel Fig. 6).L10, page 5 of 6 In any case, in spite of the potential contamination between the target populations, we determined that their stellar populations are significantly different.

Discussion and conclusion
A previous study of the central regions of the NSD (the innermost ∼20 pc × 90 pc excluding the nuclear star cluster) showed that it is dominated by old stars (more than 80% of the total stellar mass was found to be older than 8 Gyr, Nogueras-Lara et al. 2020a).In this Letter, we have analysed a similarly concentrated region towards the centre of the NSD, but considered stars located at different NSD radii along the line of sight.We selected our target samples of stars by using their different extinction, and found that they are placed on average at different Galactocentric NSD radii by analysing their proper motion distribution.Hence, the current study was only possible after Nogueras-Lara (2022a) found that the extinction within the NSD correlates with the position of the stars along the line of sight.
Our results indicate that the stellar population at different NSD radii along the line of sight is significantly different.We detected an age gradient in the NSD along the line of sight, so that the stellar population from the innermost regions is significantly older than stars from the outer regions of the NSD.In this way, we determined that there is a significant fraction of stars with ages between 2-7 Gyr in the NSD outer region that is not present in our sample deeper inside the NSD, where ∼90% of the originally formed stellar mass is older than 7 Gyr.We also observed that the contribution of even younger stars (age bins 0.5-2 Gyr and 0-0.06 Gyr) is also larger in the NSD outer region in comparison to the inner one.
To compare our results with the predominantly old stellar population detected in the central region of the NSD (Nogueras-Lara et al. 2020a), we also studied the K s luminosity function derived from all the stars in the analysed region.Figure 6 (right panel) shows the result.We found that the stellar mass is dominated by old stars (>7 Gyr) in agreement with previous studies (Nogueras-Lara et al. 2020a, 2022).We detected a somewhat lower contribution of the age bins between 0.06-2 Gyr in comparison with Nogueras-Lara et al. (2020a), which can be explained due to the different regions analysed, the different models used, and the different age bins.The contribution from the intermediate age stellar population (2-7 Gyr) is not significant when analysing all the stars in the region because the majority of stars in our sample belong to the innermost regions of the NSD.Therefore, the stars from the closest NSD edge do not contribute to the total luminosity function much.To further assess this, we computed the stellar mass from the NSD outer region in the age bin 2−7 Gyr, using Parsec models (see Sect.Methods in Nogueras-Lara et al. 2022).We ended up with a total mass of (5.3 ± 1.8) × 10 6 M in the age bin 2−7 Gyr.Comparing this value with the total stellar mass obtained when using all the stars in the target region, (1.3 ± 0.1) × 10 8 M , we conclude that the mass of stars in the age bin 2-7 Gyr only accounts for ∼4% of the total stellar mass, which is compatible with the results in Fig. 6 (right panel), and also with previous work, where the old stellar component dominates (Nogueras-Lara et al. 2020a, 2022).This is in agreement with an exponential radial scale length of the NSD, in which the stellar mass density increases towards its innermost regions (e.g., Sormani et al. 2022).
Our results, in combination with the detection of a significant stellar mass in the age bin of 2-7 Gyr in the Sgr B1 region (belonging to the NSD and located at ∼100 pc in projection from its centre), indicate the presence of an age gradient in the NSD.
Therefore, the age distribution of the NSD stellar population is similar to the one found for other NSDs in external spiral galaxies (Bittner et al. 2020), and it supports an inside-out formation of the NSD (Gadotti et al. 2020;Bittner et al. 2020).

Fig. 1 .
Fig. 1.Spitzer false colour image using 3.6, 4.5, and 5.8 µm, as red, green, and blue, respectively (Stolovy et al. 2006).The white shaded boxes indicate the target region.The position of Sgr B1 as well as the nuclear star cluster (NSC) are indicated.The circular region indicates the effective radius of the nuclear star cluster (∼5 pc, e.g., Gallego-Cano et al. 2020).

Fig. 2 .
Fig. 2. CMD K s versus H − K s .The blue and salmon coloured regions indicate the two target groups of stars with different reddening dominated on average by stars from the closest edge of the NSD (NSD outer region) and stars deeper inside the NSD (NSD inner region), respectively.The dashed boxes show the reference stars used to build the extinction maps for each of the stellar groups analysed (see Sect. 4.1).The black arrow indicates the reddening vector.

Fig. 3 .
Fig. 3. Proper motion distribution of stars from each of the target groups.The blue and salmon dots indicate stars from the NSD outer and inner regions, respectively.Only a fraction of the stars is shown to not overcrowd the plot.

Fig. 4 .
Fig. 4. Extinction maps obtained for the stars from the NSD outer (upper panel) and inner regions (lower panel).The black star indicates the position of the supermassive black hole at the centre of the Galaxy.The white pixels indicate that there is no associated extinction value due to the lack of reference stars.

Fig. 5 .
Fig.5.Analysis of the de-reddened luminosity functions corresponding to the NSD outer (left panel) and inner (right panel) regions.The reduced χ 2 of the fit is shown in each panel.The 14 theoretical models used are grouped into five broader age bins to decrease the degeneracy between models with similar ages.

Fig. 6 .
Fig. 6.Stellar populations present in the NSD outer region (left panel), inner region (central panel), and all the NSD stars in the sample (right panel).The numbers indicate the percentage of mass due to a given age bin and its standard deviation.