GALACTICNUCLEUS: A high-angular-resolution $JHK_s$ imaging survey of the Galactic centre. IV. Extinction maps and de-reddened photometry

The extreme extinction ($A_V\sim30$\,mag) and its variation on arc-second scales towards the Galactic centre hamper the study of its stars. Their analysis is restricted to the near infrared (NIR) regime, where the extinction curve can be approximated by a broken power law. Therefore, correcting for extinction is fundamental to analyse the structure and stellar population of the central regions of our Galaxy. We aim to, (1) discuss different strategies to de-redden the photometry and check the usefulness of extinction; (2) build extinction maps for the NIR bands $JHK_s$ and make them publicly available; (3) create a de-reddened catalogue of the GALACTICNUCLEUS (GNS) survey, identifying foreground stars; and (4) perform a preliminary analysis of the de-reddened $K_s$ luminosity functions (KLFs). We used photometry from the GNS survey to create extinction maps for the whole catalogue. We took red clump (RC) and red giant stars of similar brightnesses as a reference to build the maps and de-reddened the GNS photometry. We discussed the limitations of the process and analysed non-linear effects of the de-reddening. We obtained high resolution ($\sim3''$) extinction maps with low uncertainties ($\lesssim5$\,\%) and computed average extinctions for each of the regions covered by the GNS. We checked that our maps effectively correct the differential extinction reducing the spread of the RC features by a factor of $\sim2$. We assessed the validity of the broken power law approach computing two equivalent extinction maps $A_H$ using either $JH$ and $HK_s$ photometry for the same reference stars and obtained compatible average extinctions within the uncertainties. Finally, we analysed de-reddened KLFs for different lines of sight and found that the regions belonging to the NSD contain a homogeneous stellar population that is significantly different from that in the innermost bulge regions.


Introduction
The Galactic centre (GC) is an environment of fundamental importance for astrophysics beacuse it is the closest galactic nucleus, located at only 8 kpc from Earth (e.g. Gravity Collaboration et al. 2018;Do et al. 2019), and we can analyse its structure and stellar population, resolving individual stars down to milliparsec scales. It hosts a supermassive black hole at its centre, Sgr A*, which is embedded in a nuclear stellar cluster (e.g. Ghez et al. 1998;Schödel et al. 2002;Genzel et al. 2010;Schödel et al. 2014). At larger scales, the nuclear stellar disc (NSD) and the central molecular zone (CMZ) outline the limits of the GC (e.g. Morris & Serabyn 1996;Launhardt et al. 2002;Kruijssen et al. 2014;Nogueras-Lara et al. 2020a).
The extreme extinction towards the GC (A V 30 mag, A K s 2.5 mag, e.g. Nishiyama et al. 2008;Schödel et al. 2010;Nogueras-Lara et al. 2018a, 2020b impedes the study of its stellar populations in the optical. The extreme source crowding requires the use of high-angular-resolution observations to disentangle individual stars. Even previous studies using the VISTA The extinction maps and the catalogues described in this paper and specified in Sect. 2.1 and Table 3 are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.
Variables in the Via Lactea (VVV) survey (Minniti et al. 2010;Saito et al. 2012), which constitutes the most accurate up-to-date survey for the Galactic bulge, are not sufficient to properly determine the extinction variations in the innermost regions of our Galaxy (e.g. Gonzalez et al. 2012Gonzalez et al. , 2013Surot et al. 2020). To improve this situation, the GALACTICNUCLEUS (GNS) survey (Nogueras-Lara et al. 2018a, 2019a was optimised for the requirements of the GC (JHK s , high angular resolution, high dynamic range) and is currently the most complete near infrared (NIR) JHK s photometric catalogue of the GC. It provides highangular-resolution (∼ 0.2 ) photometry of more than 3 million stars distributed along the NSD, the inner bulge, and the transition region between them.
The extinction towards the GC varies significantly on arcsecond scales (Scoville et al. 2003;Schödel et al. 2010;Nogueras-Lara et al. 2018a) due to the presence of a significant amount of gas, dust, and molecular dark clouds associated with the CMZ (e.g. Dahmen et al. 1998;Pierce-Price et al. 2000). In order to quantify the magnitude of the extinction variations within the field of view and to trace the presence of the dark clouds, extinction maps can be extremely useful. They inform us about the properties and distribution of the interstellar medium (ISM), given that they are directly related to the dust content. Moreover, they can be used to determine the gas column density Article number, page 1 of 17 arXiv:2107.00021v2 [astro-ph.GA] 29 Jul 2021 A&A proofs: manuscript no. Nogueras-Lara distribution assuming a constant gas to dust ratio (e.g. Goodman et al. 2009;Rowles & Froebrich 2009).
To create extinction maps, it is key to understand the behaviour of the extinction curve. It has been generally assumed that the extinction curve in the NIR can be well approximated as a power law (e.g. Rieke & Lebofsky 1985;Nishiyama et al. 2006;Fritz et al. 2011). However, recent studies found a dependence of the extinction on the wavelength, which points towards a more complicated broken-power-law behaviour of the form A λ ∝ λ −α JH and A λ ∝ λ −α HKs , where α i is the extinction index and the i-subindinces indicate their range of validity between the photometric bands JH and HK s (Nogueras-Lara et al. 2018a;Hosek et al. 2018;Nogueras-Lara et al. 2019b, 2020b. Recent work by Nogueras-Lara et al. (2020b), carried out using the whole GNS survey, obtained α JH = 2.44 ± 0.05 and α HK s = 2.23 ± 0.05.
In this paper, we compare the de-reddening of NIR photometry using a star-by-star basis and an extinction map approach. Moreover, we also create the most detailed high-angularresolution extinction maps of the GC (covering a total region of ∼ 6000 pc 2 ) using the GNS survey. We validate the brokenpower-law approach and de-redden the photometry from GNS to obtain de-reddened colour-magnitude diagrams (CMDs) of the GC. We publish the de-reddened catalogue, analyse the limitations of the de-reddening, and study K s luminosity functions (KLFs) to carry out a preliminary study of the stellar population for the different regions covered.

Data
For the study presented in this paper, we used the GNS survey. This is a JHK s NIR high-angular-resolution (∼ 0.2 ) catalogue, obtained using the HAWK-I instrument (Kissler-Patig et al. 2008) at the ESO Very Large Telescope unit telescope 4. It covers an area of ∼ 6000 pc 2 distributed along seven regions that correspond to the nuclear stellar disc (Central, NSD East, and NSD West), two regions in the inner bulge (IB South and IB North), and two regions in the transition region between the NSD and the inner bulge (T East and T West). Figure 1 shows a scheme of the areas covered by the survey. The GNS catalogue was calibrated using the SIRIUS/Infrared Survey Facility telescope (IRSF) GC survey (e.g. Nagayama et al. 2003;Nishiyama et al. 2006) and has a systematic zero point (ZP) uncertainty of 0.04 mag in all three bands. It includes accurate photometry of ∼ 3.3 × 10 6 sources that reaches 5 σ detections for J ∼ 22, H ∼ 21, and K s ∼ 21 mag. The photometric uncertainties are below 0.05 mag at J 20, H 17, and K s 16 mag (Nogueras-Lara et al. 2019a).

Fixing duplicated stars
We noticed that some stars are duplicated in the final GNS catalogue. Slightly different positions for different detectors, pointings and/or bands, or the presence of very close neighbouring stars might explain this problem. The observing strategy to obtain the GNS survey includes 49 pointings. Each pointing is composed of the combination of the four HAWK-I detectors, which are separated by a cross-shaped gap between them. To cover this gap, we applied random jittering to the observations and computed the photometry on an individual detector basis (for further details, see Nogueras-Lara et al. 2019a). In this way, there are many overlapping regions between different detectors and pointings. The catalogues for each band were obtained combining all the detections from each independent detector and then  N). The white rectangles in the Central and I B S regions indicate a gap in the data due to bad quality (see text for details). The background image corresponds to a 3.6 µm Spitzer/IRAC image (Stolovy et al. 2006). from different pointings. We lastly combined all three bands into a final catalogue. Therefore, we believe that some stars escaped to be properly combined throughout the process.
We identified the duplicated stars and combined them into single sources. We checked whether they were detected in one or several bands and combined the coordinates and photometry accordingly, following the catalogue scheme shown in Table 2 of Nogueras-Lara et al. (2019a). The uncertainties were quadratically propagated from the individual detections. To identify these stars in the corrected catalogue, we modified the parameters i J , i H , and i K s , which correspond to the multiple detections in overlapping pointings for each band (see Sect. 4.4 for further details; Nogueras-Lara et al. 2019a), and set them to −1.
The number of duplicated stars is very low, so they do not affect any of the previous studies carried out using the GNS survey. Namely, for each of the catalogue regions, we found a rate of duplicated stars of: Central ∼ 0.002 %, NSD East ∼ 0.01 %, NSD West ∼ 0.02 %, T East ∼ 0.04 %, T West ∼ 0.01 %, IB South ∼ 0.02 %, and IB North ∼ 0.04 %. Figure 2 shows the colour-magnitude diagrams (CMDs) H versus J − H and K s versus H − K s for the different regions of the GNS survey. The CMDs are largely complete in the red clump (RC) region, which is key to produce extinction maps (e.g. Nogueras-Lara et al. 2018a, 2020b. Red clump stars are red giant stars in their helium burning phase whose intrinsic properties are well known (e.g. Girardi 2016). The over-densities around H ∼ 16 − 17 mag and K s ∼ 15 − 16 mag correspond to the RC features.

Colour-magnitude diagrams
To produce extinction maps, we removed foreground stars using the significantly different extinction between the stars from the Galactic disc and the GC (e.g. Nogueras-Lara et al. 2019a). We applied a colour cut, J − H ∼ 2.35 or H − K s ∼ 1.1 mag, depending on the CMD used (see Sect. 2 in Nogueras-Lara et al. 2021a).
We also defined another colour cut to remove stars belonging to the Galactic bulge in the GNS regions covering the NSD (Cen- Fig. 2. Colour-magnitude diagrams for each region of the GNS survey. The dashed lines represent the cuts applied in this work between the foreground population (in blue) and stars that we consider to lie inside the GC, behind two different layers of extinction (orange and red). The extinction layers were selected following the methodology explained in Sect. 3.2. The dashed rectangles indicate the selection of reference stars (mainly RC stars) to create the extinction maps. Only two layers are defined for the IBN and IBS regions. The low number of stars present for the NSD W region is due to the lower than average quality of the data for that region. Given the high number of stars, only a fraction of them is represented to not overcrowd the plots. tral, NSD E, NSD W, T E, and T W). This is possible because the stellar population from the GC lies behind a dense screen of extinction that may be related to the CMZ (e.g. Schultheis et al. 2009). We adopted H − K s ∼ 1.3 mag following previous results on the inner bulge (e.g. Nogueras-Lara et al. 2018b, 2020aSchultheis et al. 2021). Nevertheless, we treated each region in a different way due to the variation of the absolute extinction between regions. Given that we produced independent extinction maps using JH and HK s bands, we excluded the foreground population considering independent cuts in both CMDs H versus J − H and K s versus H − K s . Figure 2 shows the cuts applied to remove the foreground population in each case. We estimate that the possible contamination from stars belonging to the Galactic bulge/bar is below 25% using the model created by Schultheis et al. (2021), which applied similar colour cuts to analyse the NSD. In particular, they modelled the NSD, the bulge/bar, and the nuclear star cluster (NSC) using the results from Launhardt et al. (2002), Chatzopoulos et al. (2015), and Sormani et al. (2020).
We adopted a two-layers approach to create the extinction maps to overcome the extreme differential extinction in the GC (see Sect. 3.2). Figure 2 indicates the criteria to assign the stars to each of the presumed extinction layers.

Extinction maps
We created extinction maps for all the regions in the GNS catalogue following the approach described in detail in Ap-pendix A.2. We obtained extinction maps A J, JH and A H, JH , and A H, HK s and A K s , HK s , using JH and HK s photometry, as indicated by the subindices. We assumed the extinction curve A J /A H = 1.87 ± 0.03 and A H /A K s = 1.84 ± 0.03 derived in Nogueras-Lara et al. (2020b). We computed associated uncertainty maps using a Jackknife algorithm, estimating the extinction variation for a given pixel and systematically leaving out each of the stars used for its extinction calculation (for further details see Sect. 7 in Nogueras-Lara et al. 2018a). The systematic uncertainties were estimated by analysing the fluctuations of the extinction maps when varying the quantities involved in their calculation within their uncertainties. These are: A J /A H = 1.87 ± 0.03 and A H /A K s = 1.84 ± 0.03 (depending on the computed map), the ZP photometric uncertainty of the GNS data (0.04 mag for JHK s ) Nogueras-Lara et al. 2019a), and the intrinsic colours JH and HK s for RC stars and their associated uncertainties.

Intrinsic colour of red clump stars
The reference stars used to compute the extinction maps are mainly RC stars with a fraction of other red giant stars (red giant branch bump, RGBB, stars among others, e.g. Nogueras-Lara et al. 2018b), whose intrinsic colours are very similar. Therefore, we computed their intrinsic colours (J − H) 0 and (H − K s ) 0 using RC stars as a reference. We estimated these intrinsic colours and their uncertainties with three independent methods: (1) We selected stars in the RC features from the synthetic model created in Sect. A.1 for the NSD and estimated average intrinsic colours (J − H) 0 = 0.57 ± 0.04 mag and (H − K s ) 0 = 0.11 ± 0.01, where the uncertainties correspond to the standard deviation of the measurements.
(2) We used the values (J − H) 0 = 0.47 mag and (H − K s ) 0 = 0.09 mag computed following the method described in Sect. 3.2 of Nogueras-Lara et al. (2020b), using a Kurucz stellar model (Kurucz 1993) for a RC star (effective temperature of 4750 K, twice solar metallicity, and logg = +2.5) and the corresponding HAWK-I filters.
(3) Finally, we used the recent observational results (J − H) 0 = 0.52 ± 0.04 mag and (H − K s ) 0 = 0.10 ± 0.03 mag from Table 3 in Plevne et al. (2020), obtained for metal rich RC stars. The metallicity of the used metal rich RC stars is slightly slower than that expected for the GC (e.g. Feldmeier-Krause et al. 2017;Schultheis et al. 2019;Nogueras-Lara et al. 2020a), but is appropriate to obtain an estimation of the intrinsic colours by comparing the values with the previous ones.
We combined the previous values and estimated the uncertainties using the standard deviation of the distributions. The results are (J − H) 0 = 0.52 ± 0.04 mag and (H − K s ) 0 = 0.10 ± 0.01 mag.

Extinction layers
Building an extinction map involves the simplification of an actual tridimensional structure into a 2D map that corrects for extinction. In the GC, where the differential extinction is extreme and significantly varies along the line of sight (e.g. Nishiyama et al. 2006Nishiyama et al. , 2009bNogueras-Lara et al. 2018a), this is even more complex. Moreover, the steep rise of extinction towards short wavelengths in the NIR (>5 mag higher in J than in K s ) implies that the number of stars detected is significantly lower in J than in H and K s (∼ 20% of the stars were detected in J; ∼ 65%, in H; and ∼ 90%, in K s , Nogueras-Lara et al. 2019a). Therefore, the catalogue reaches stars in H and K s that are too extinguished to be detected in J.
In order to create extinction maps to de-redden the photometry in a consistent way, we decided to take a two-layers approach that allowed us to account for the variation of the extinction along a given line of sight and also to avoid over-corrections in the case of the J band. The two layers defined do not necessarily correspond to real extinction screens but are more related to the incompleteness of the J band and the robustness of the dereddening of the stellar population. In this way, and given that the stars used to create the extinction maps mainly correspond to the NSD (see Sect.2.2), both extinction layers are expected to lie at the same distance. The difference in extinction and obscuring material is then due to the extreme differential extinction and the 'granularity' characterising the NSD and the CMZ (e.g. Gosling et al. 2006).
The first layer of extinction is defined considering the stars detected in all three bands (JHK s ), whereas the second one is based on the stars that are not detected in the J band. The colour cut to distinguish the first and the second extinction layers in the CMD K s versus H − K s was chosen using the HK s counterparts of stars detected in the J band. In this way, we plotted those stars in the CMD K s versus H − K s and defined the cut for the second layer (see Fig. 2). Given the lower completeness of the J band and the intrinsic scatter of the data, not all the stars belonging to the first layer in the CMD K s versus H − K s exactly correspond to the stars in the same layer in CMD H versus J − H. Because of the steep extinction curve, completeness is a strong function of wavelength. This results in more stars belonging to the first layer detected towards faint magnitudes for the CMD K s versus H − K s . For the Central region in particular, around two times more stars are included in the first layer for the CMD K s versus H − K s than for the H versus J − H one.
For the first layer we computed four extinction maps. Using the ratio A J /A H and the JH photometry, we created a J-and an H-map. Using the A H /A K s and the HK s photometry, we created an H-and a K s -map. For the second layer we computed two extinction maps, corresponding to the H and K s bands, using the extinction ratio A H /A K s .
The reference stars to create each of the layers were selected in the CMDs K s versus H − K s (for each region), given that the vast majority of stars within the two layers are present there. For the first extinction layer, we considered all the stars observed in the J band that have a counterpart in H and K s . To select the limit between both layers, we considered the cut where the stars are observed in all three bands or only in H and K s . Figure 2 shows each of the established divisions that are depicted by the black dashed rectangles. The lower K s cut of the selected stars for the second layer is fainter given that the extinction is higher here, producing redder stars with fainter K s magnitudes.
The situation is simpler for the regions belonging to the inner bulge (I B S and I B N), and a single-layer approach is sufficient. This is due to the lower crowding and extinction, together with the less significant differential extinction (e.g. Nogueras-Lara et al. 2018b, 2020b). Figure 3 shows the A K s maps and their associated uncertainties obtained for the Central region of the GNS survey, as an example. The last panel of the figure shows the false colour (JHK s ) image of the analysed region, where it is possible to establish some correlation between the extinction maps and the dark regions in the image where the stellar density is lower. This correlation is much stronger for the second reddening layer, given that the first one correspond to stars whose extinction is lower. On the other hand, there is also some correlation between both extinction layers as it is not possible to completely distinguish stellar populations without any overlap (e.g. Nogueras-Lara et al. 2018a). The error maps show relatively constant uncertainties that are somewhat larger in regions where there are stars that might have different extinctions and are close in projection.

Results
We observed the strong variations caused by the differential extinction in both layers; they are more significant and vary on smaller scales for the second extinction layer. This is in agreement with the variations found for the central pointing of the GNS survey by Nogueras-Lara et al. (2018a) and corresponds to the high level of complex structure or 'granularity' observed by Gosling et al. (2006), associated to extinction variations with a characteristic scale of 5 − 15 .
All the maps obtained for the GNS regions are made publicly available at the Centre de Données astronomiques de Strasbourg (CDS). Table 1 indicates the main properties of the computed maps. We computed the uncertainties for each of the maps and ended up with statistical and systematic uncertainties 5 % and 6 %, respectively. We also computed the average extinction for the stars in each layer and its variation for a given band via its standard deviation (Table 1).
The maps for the Central region and the IBN include significant gaps (see Fig. 3 for the Central region), given that two of the pointings that form part of these regions were observed under very bad weather conditions and were therefore not included in the GALACTINUCLEUS catalogue.

Quality of the extinction maps and extinction curve assessment
Since we created two independent extinctions maps A H for the first extinction layer, we could compare them with each other to assess their reliability as well as the consistency of the extinction curve that we had assumed. Figure 4 shows both maps for the Central region. Qualitatively, they share the main extinction features without important disagreement within the uncertainties (dA H ∼ 0.20 mag, estimated by propagating the statistical and systematics uncertainties in Table 1). Moreover, we computed the average extinction obtained for the stars corresponding to each map ( Table 2) and concluded that the values agree within the uncertainties. The only exception is the NSD W region where the lower than average quality of the data for this field (see Table  A.3. in Nogueras-Lara et al. 2019a) makes the results unreliable.
In general, we observed some tendency towards larger values for the mean extinctions computed by using the HK s bands (so we used slightly different scales for each map in Fig. 2). We believe that a possible variation of the ZP, within the systematic uncertainties (∼ 0.04 mag in all three bands, Nogueras-Lara et al. 2019a), might explain this systematic variation. In particular, considering the mean A H values obtained for the Central region (

Previous 2D extinction maps
Previous extinction maps for the GC region were limited to the photometric surveys covering this region. In this way, they are normally affected by low photometric completeness that impedes a good coverage of the RC feature. In particular, Schultheis et al. (1999Schultheis et al. ( , 2009) used the Deep Near-Infrared Survey of the southern sky (DENIS) survey and Spitzer data to obtain extinction maps with resolutions ∼ 2 . In spite of the very different resolutions (∼ 3 for the extinction maps computed in this paper), we compared the extinction map obtained in Schultheis et al. (2009) with our results. We used all the stars in the Central region of the GNS survey and obtained their associated extinction using the combined extinction map from Schultheis et al. (2009) (see their Sect. 3.2 for further details). Averaging over the obtained values, we ended up with A K s = 2.74 ± 0.72 mag, where the uncertainty corresponds to the standard deviation of the distribution. We used the same stars and also computed the associated extinction using our two-layers approach. We obtained a mean value of A K s = 2.30 ± 0.46 mag, where the uncertainty was estimated via the standard deviation of the measurements. The somewhat larger value obtained when using the extinction map from Schultheis et al. (2009) can be explained due to the different extinction curve used to compute the extinction (Glass 1999). Converting the extinction value using the extinction curve by Nogueras-Lara et al. (2020b) gives a value of A K s = 2.29 ± 0.60 mag, which perfectly agrees with our result within the uncertainties. More recently, the use of the VVV survey (Minniti et al. 2010;Saito et al. 2012) improved the situation and allowed us to use RC stars to obtain high resolution extinction maps (∼ 2 ) for the whole Galactic bulge. Nevertheless, the high crowding for the GC region limits the detection of stars with the VVV survey (see Fig. 12 in Nogueras-Lara et al. 2019a), so the associated extinction map does not prove adequate to de-redden the GC stellar population. In spite of it, we attempted to use the extinction map computed by Gonzalez et al. (2012) using the VVV survey. We used the BEAM 1 (Bulge Extinction And Metallicity) calculator and a selection of ∼ 1 % of the stars in the Central region of the GNS survey (∼ 20, 000 stars). We could not use more stars due to the limit on the number of stars to compute the extinction using the BEAM calculator. Given the limitations of the VVV survey for the very central region of the Galaxy, the extinction was computed for only ∼ 60 % of the total number of stars. We used the extinction curves by Nishiyama et al. (2009a) and Cardelli et al. (1989) and obtained A K s = 2.18 ± 0.24 and A K s = 2.85 ± 0.31 mag, respectively. The uncertainties were estimated as the standard deviation of the distributions. The significantly different values obtained illustrate how different extinction curves influence the determination of the absolute extinction (e.g. Nogueras-Lara et al. 2019b). To better compare the results with the work presented here, we translated these values into A K s = 1.69 ± 0.18 mag using the extinction ratio A J /A K s = 3.44 obtained in Nogueras-Lara et al. (2020b). This is in good agreement with the value A K s = 1.82 ± 0.14 mag obtained in Table 1 for the first extinction layer. This lower extinction, which only reaches the first extinction layer as defined in this paper, is not surprising given the limiting magnitude when using aperture photometry from the VVV survey and RC stars to compute the extinction in the NSD (see Fig. 12 in Nogueras-Lara et al. 2019a). White regions indicate that there were not enough stars to assign a value. The lower panel is a false colour image using JHK s bands. The dashed white lines indicate dark clouds that can be identified in the extinction map corresponding to the second layer.
Dedicated high-angular-resolution observations for particular GC regions make the situation different for small interesting regions like the nuclear star cluster. In this way, using NACO (NAOS-CONICA) data, Schödel et al. (2010) obtained a high resolution extinction map, ∼ 1 − 2 , for a region with a radius ∼ 20 centred on Sagittarius A*. The mean extinction towards the central parsec around Sgr A* was estimated to be A K s = 2.54 ± 0.12 mag (Schödel et al. 2010). To compare this value with our results, we considered stars in a radius of ∼ 13 and estimated their extinction combining the first and the second extinction layers defined in this paper. We obtained a mean value of A K s = 2.40 ± 0.28 mag, in good agreement with the previous value.

Comparison with 3D extinction maps
We aimed to compare our extinction maps with more general 3D extinction maps, in spite of their different resolutions and depths. In this way, this comparison is limited by their typical angular resolution of several arcminutes and by the sensitivities of the surveys used (see Sect. 1 in Schultheis et al. 2009). This implies that some 3D extinction maps, like the one derived by Marshall et al. (2006), are limited by relatively low extinctions (A v ∼ 25 mag) and by the underlying models used to build them (e.g. Drimmel et al. 2003). Therefore, they cannot be directly compared with our results.
On the other hand, we used the more recent 3D extinction map obtained by Chen et al. (2013), whose extinction values for the inner bulge and the central region appear to be in agreement with deeper 2D studies (see Sect. 5.3 in Chen et al. 2013). The extinction map is presented at a resolution of 15 × 15 for six different combinations of colours in distance bins of 1 kpc. We chose the lines of sight towards the Central region of the GNS survey and computed average values of the colour excesses E(J − K) = 3.99 ± 0.23 mag and E(J − K) = 1.34 ± 0.12 mag, at the GC distance of ∼ 8 kpc. The uncertainties correspond to the standard deviation of the values obtained for the different lines of sight. Using the extinction curve obtained by Nogueras-Lara et al. (2020b), we transformed the colour excesses and obtained A K s JK s = 1.64±0.09 mag, and A K s HK s = 1.59±0.14 mag. The subindices indicate the used colour excess. These values are in agreement with A K s = 1.82 ± 0.14 obtained in Table 1 for Notes. The super-indices 1 and 2 indicate the first and the second extinction layer, respectively. The mean extinction for the stars corresponding to each extinction map is indicated byĀ λ . The associated standard deviation, σĀ λ , shows the variation of the extinction within a layer for a given band. The symbols ∆ statĀλ and ∆ systĀλ indicate the average statistical uncertainty computed using a Jackknife algorithm for a given band and layer, and the systematic uncertainty of the corresponding extinction map, respectively. In all the cases the first H 1 corresponds to the extinction layer obtained using the JH bands, whereas the second one was computed using HK s . The results obtained for NSD W are less reliable due to bad data quality caused by bad weather conditions. the first layer of the Central GNS region, where the uncertainty corresponds to the standard deviation of the extinction values. We adopted the first layer for comparison due to the limitations imposed by the Two Micron All Sky Survey (2MASS) and the Galactic Legacy Infrared Midplane Survey Extraordinaire II (GLIMPSE-II) used to compute the 3D extinction maps.

Previous work using the GNS survey
We also compared our results with previous work carried out using the GNS survey, which calculated average extinction values for different regions: (2) Nogueras-Lara et al. (2018b) computed A K s = 1.14 ± 0.09 mag and A K s = 1.39 ± 0.09 mag for pointings #B1 and #B2 from the inner bulge regions I B N and I B S, respectively, where the uncertainties correspond to the systematics.
We found that the average values agree well within the uncertainties with the extinction values presented in Table 1 for the Central region (Ā K s = 1.82 ± 0.12 mag), the IBN (Ā K s = 1.18 ± 0.09 mag), and the IBS (Ā K s = 1.49 ± 0.10 mag). The uncertainties were obtained propagating the statistical and systematic uncertainties of the average extinction values.

De-reddening
We de-reddened the GNS photometry using the two-layers approach described previously, and obtained a catalogue of dereddened sources.

Results
We excluded the foreground population and assigned each star to its corresponding extinction layer following the criterion explained in Sect. 3.2, as indicated in Fig. 2. For each band and layer, we obtained de-reddened photometry independently applying the obtained extinction maps. For stars in the first extinction layer, we applied the extinction maps A J, JH and A H, JH , (obtained using JH bands, as the subindices indicate), and A H, HK s and A K s , HK s , (obtained using HK s bands). The result for each star is an associated de-reddened value for J and K s photometry, and two de-reddened values for the H band, corresponding to the extinction maps computed using either JH or HK s photometry.
For stars in the second layer, we applied the second layer extinction maps A H and A K s , and computed the de-reddened photometry for the H and K s bands.
The results obtained for the Central region are shown in Fig. 5. The spread of the data points and the width of the standard distribution of the RC sources is ∼ 2 times lower in the corrected diagrams. Therefore, the differential extinction is significantly corrected in all the cases.
We over-plotted a Parsec isochrone (release v1.2S + COL-IBRI S_35 + PR16, Bressan et al. 2012;Chen et al. 2014Chen et al. , 2015Tang et al. 2014;Marigo et al. 2017;Pastorelli et al. 2019) of an old (∼ 10 Gyr) metal rich (twice solar) stellar population, which is the dominant stellar population in the NSD (Nogueras-Lara et al. 2020a), to check that the de-reddening was appropriate. We observed a deviation of the de-reddened data from the model at the faint end (K s ∼ 15 − 16 mag) due to the incompleteness of the data (Nogueras-Lara et al. 2020b). The de-reddened CMDs corresponding to the remaining regions are shown in Appendix B (Figs. B.1 and B.2).

Catalogue
We created a new version of the GNS catalogue by adding the results obtained for extinction. We present an independent catalogue for each of the regions of the survey in the same way as we did in Nogueras-Lara et al. (2019a). We followed the scheme shown in Table 2 of Nogueras-Lara et al. (2019a) and used the new catalogue corrected for duplications (see Sect. 2.1). For simplicity, we only included the absolute coordinates of each source (RA, ∆RA, DEC, and ∆DEC) and omitted the coordinates corresponding to the detection of each source in any given band (they are available in the original catalogue corrected for duplications). We also included the photometry and the corresponding uncertainties for all three bands. A value of 99 for the photometry indicates non-detection for a given source in a given band.
We added two columns indicating whether a star is detected as foreground population using the methodology explained in Sect. 3.2. Since the identification was carried out using the CMDs H versus J − H and K s versus H − K s , we specified the detection as foreground population for a given source in these two diagrams using F JH and F HK s , respectively. A value of 1 indicates that a source is detected as foreground population, whereas a value of -1 indicates that the source corresponds to the target population (GC stars for Central, NSD E, NSD W, T E, and T W, and inner bulge population for I B S and I B N). Due to the scatter associated to the photometric uncertainties and the degeneracies between different populations because of differential extinction, it is unavoidable that the results for some sources disagree between F JH and F HK s . This means that some sources might be catalogued as foreground using F JH , but are not using F HK s and vice versa.
We also added six columns corresponding to the extinction values for each band and layer (A i band diagram ), and six columns with the associated uncertainties (A i band diagram ). The super-index i refers to the extinction layer, whereas the subindex band diagram gives information about the band of the extinction maps and the photometric bands used to compute the extinction maps from the reference stars. As it happens to the foreground population, it might occur that some stars are detected as belonging to the first layer using a given reference diagram and to the second one using the other. For instance, a star might belong  to the first layer using JH as a reference, but be identified as from the second one when using HK s . This is due to the intrinsic scatter of the GC population and cannot be avoided. A value of -1 indicates that a source does not belong to a given extinction layer. A value of 0 indicates that a source is too close to the edge of the corresponding extinction map and thus the extinction value was not computed. A value of 'nan' indicates that there is no associated value for the extinction map given the low number of reference stars.
The two-layers approach was not necessary in the case of the inner bulge regions (I B S and I B N), as was explained previously (see Sect. 3.2 for further details). Therefore, only one layer is specified (four columns for extinction and four for its associated uncertainties). Table 3 shows part of the catalogue obtained for the Central region as an example. All the catalogues presented in this paper are made publicly available at the CDS.

Limitations of the de-reddening
We studied the limitations of the de-reddening process, giving the adopted strategy, and its application to different stellar types analysing possible non-linear effects.

Extinction layers
Extinction is a 3D phenomenon that needs to be simplified to be taken into account in an affordable way. The reference stars used to create the extinction maps are distributed along the line of sight. Although the distance does not affect the measured stellar colours directly, the amount of reddening material towards the stars is a function of their distance. We used the two-layers approach, dividing the stars by means of a colour cut in the CMDs. In this way, we avoided mixing reference stars whose absolute extinctions are too different. This approach can be generalised to be applied to the whole survey, as explained in previous sections. It gives reasonably good average values for the whole stellar population (see Table 1). Nevertheless, it also produces over-or under-corrections for some of the stars, as shown by the scatter of the de-reddened CMDs. To correct specific regions where the differential extinction is too large, therefore, a more convenient approach would be to impose constraints on the colour range to accept reference stars for a given extinction-map pixel in order to avoid mixing stars with too different extinctions (see Methods section in Nogueras-Lara et al. 2020a).
The extinction layers that we present here do not necessarily correspond to independent physical layers. They are the result of the incompleteness of the J band data and serve as a boundary to avoid computing the extinction maps using stars that are close in projection but suffer from very different extinctions. We recommend therefore that de-reddened stars from independent layers should not be mixed in order to avoid possible systematics related to the cuts introduced when defining the layers. In those cases, the best approach is to create dedicated extinction maps following the approach described in Nogueras-Lara et al. (2020a), imposing constraints on the colour differences for the neighbouring stars to avoid too different extinctions.
The GC is a complex and clumpy region characterised by the presence of dark clouds that produce strong extinction variations on arc-second scales (e.g. Launhardt et al. 2002;Nogueras-Lara et al. 2018a;Battersby et al. 2020). To account for these variations, we combined a number of reference stars per pixel, restricting the distance from a given pixel to the reference stars and introducing the IDW method to weight the stars depending on their distance to a given pixel. Moreover, the two-layers approach helps us to avoid the combination of stars with very different extinctions but close in projection.

Non-linear effects
There are some non-linear effects, also known as bandwidth effects, which might affect the extinction correction that we present here (Jones & Hyland 1980;Straižys & Lazauskaitė 2008;Maíz Apellániz et al. 2020).
One of the most relevant effects is that the extinction is not linear with the stellar type. We built extinction maps using RC stars and applied them to the target stars to correct their photometry. This means that other types of stars, different from the red

Notes.
Example of the final catalogue. Rows 120-142 correspond to the Central region of the GNS survey. The position and photometry of the sources were taken from the original catalogue (see Table 2 of Nogueras-Lara et al. 2019a). The right ascension and the declination of the sources are RA, ∆RA, DEC, and ∆DEC. The coordinates were rounded to two decimals in this table. The photometry and the associated uncertainties are indicated by J, dJ, H, dH, K s , and dK s . A 99 indicates that there is no detection for a given source. F i indicates whether a source is detected as foreground population (=1) using the bands specified as a reference (i). A value of -1 refers to a non-foreground population. A i j and dA i j correspond to the value of the extinction and its associated uncertainty. The subindices i and j are the bands used to compute the extinction maps and the extinction layer, respectively. A value of -1 indicates that a star does not correspond to a given layer. A value of 0 indicates a star too close to the borders of the corresponding extinction map. A 'nan' value indicates that there is no associated value. giants used as a reference, might be affected by the extinction in a slightly different way that is not considered in our correction. To analyse this effect, we computed the extinction corresponding to the JHK s bands for stars with different temperatures, using the intrinsic and reddened H − K s colours. We used colours instead of magnitudes because they do not depend on distance and stellar radius. To compute them, we assumed the corresponding Kurucz models (Kurucz 1993), the HAWK-I filter curves (Kissler-Patig et al. 2008), the extinction curve derived in Nogueras-Lara et al. (2020b), and a monochromatic amount of extinction of A 1.61 = 3.40 mag for all of them (Nogueras-Lara et al. 2020b). In the latter, the subindex 1.61 indicates the monochromatic wavelength 1.61 µm and the selected value corresponds to the average extinction to the central region of the GNS survey following Table 2 in Nogueras-Lara et al. (2020b). We used Eq. A.1 to calculate the A K s values, substituting the intrinsic and the de-reddened colours previously computed. Table 4 shows the results. We estimated whether the differences between stellar types are significant considering an uncertainty of ∼ 10 % for the intrinsic colours (as for RC stars, see Sect. 3.1). We obtained that there is no significant difference within the uncertainties for any of the stellar types considered here in comparison with RC stars. Moreover, we translated the A K s values into A H and A J extinctions using the ratios A H /A K s = 1.84 ± 0.03 and A J /A K s = 3.44 ± 0.08 (Table 4). The difference between values is larger for shorter wavelengths but the uncertainties also grow accordingly (∆A H ∼ 0.14 mag and ∆A J ∼ 0.28 mag) when considering the corresponding extinction ratios. Thus, the effect is not significant within the uncertainties.
A further small side effect is the influence of the effective wavelength when using wide-band filters to produce the extinction maps. The ratios between A H /A K s and A J /A K s that we used in this paper were computed using RC stars as a reference (see Nogueras-Lara et al. 2019b, 2020b. In reality, these quantities would vary for each star depending on its stellar type. To estimate the variation, we computed the effective wavelength for each of the stellar types in Table 4 using the method describe in Nogueras-Lara et al. (2018a). We then used the equation A λ 1 /A λ 2 = (λ 1 /λ 2 ) −α λ 1 λ 2 (Nogueras-Lara et al. 2020b) to estimate the influence of the effective wavelength variation on the extinction ratios. Averaging over the values obtained for the used stellar types, we ended up with A J /A H = 1.88 ± 0.01 and A H /A K s = 1.84 ± 0.01, where the uncertainties refer to the standard deviation of the distributions. The results fully agree with the extinction ratios derived in Nogueras-Lara et al. (2020b). Therefore, we concluded that the effect of the effective wavelength on the de-reddening is not significant within the uncertainties.

De-reddened K s luminosity functions
To check the usefulness of the extinction maps and the dereddened catalogues, we performed a preliminary analysis by building KLFs for each of the regions in the GNS catalogue. We restricted the analysis to stars from the first extinction layer (A 1 K s , HK s > 0, as indicated by the super-index), where the completeness is higher than in the secondary layer, and the foreground population is excluded (see Fig. 2). We obtained de-reddened KLFs with a completeness larger than 80 % at K s0 ∼ 14.5 mag (according to the completeness ∼ 80 % at K s ∼ 16.3 mag, computed for the Central region in Nogueras-Lara et al. 2020b). Only the KLF from the NSD W shows lower completeness and a significantly lower number of detected stars due to lower than average data quality for this region (see Table  A .3. in Nogueras-Lara et al. 2019a). Figure 6 shows the corresponding KLFs. The regions used to produce the KLFs are larger than those used in previous studies by a factor of ∼ 3 for the NSD regions (Nogueras-Lara et al. 2020a) and ∼ 2 for the innermost bulge (Nogueras-Lara et al. 2018b). The transition regions had not been previously studied using the KLFs.
The main features appearing in the KLFs correspond to the asymptotic giant branch bump (AGBB), the RC, and the red giant branch bump (RGBB) (see Fig. 2 in supplementary material of Nogueras-Lara et al. 2020a). Analysing the relative fraction of stars between them, it is possible to infer the star formation history (SFH) of the regions (e.g. Nogueras-Lara et al. 2020a).
A preliminary analysis indicates significant differences between the KLFs corresponding to the different regions studied. We observed that the shape of the KLFs varies with latitude, being similar for regions within the same structure. Namely, the KLFs corresponding to regions in the NSD, the inner bulge, and the transition region between the NSD and the bulge present similar features. This points towards physically different components with different SFHs, in agreement with previous work (e.g. Launhardt et al. 2002;Nishiyama et al. 2013;Nogueras-Lara et al. 2018b;Gallego-Cano et al. 2020;Nogueras-Lara et al. 2020a). In particular, the RC feature appears to be thicker for the NSD regions (excluding the NSD W, whose lower completeness avoids a good coverage for faint K s magnitudes), which is probably the result of the presence of a secondary RC feature due to a star formation event ∼ 1 Gyr ago that was responsible for the formation of ∼ 5 % of the stellar mass (Nogueras-Lara et al. 2020a).
The double RC feature does not appear in the case of the inner bulge regions, where the RGBB feature is more prominent (and the AGBB less important), indicating a different SFH dominated by mainly old stars according to previous work on the inner bulge (Nogueras-Lara et al. 2018b). Moreover, the AGBB is more prominent in the NSD regions, indicating a larger contribution of young stars (several hundred million years) to the KLF. On the other hand, the KLFs corresponding to the transition regions share characteristics from both, the NSD and the inner bulge fields, though are more influenced by the latter. Thus, this region appears as a border between these two different components of the innermost region of the Galaxy.  (Harris et al. 2020). The RC, AGBB, and the RGBB are indicated in the panels, where the features are detected. The coloured error bars on the top of the panels indicate the systematic uncertainty due to the de-reddening process and the ZP systematics. The green KLF corresponding to the NSD W suffers from low completeness due to lower than average data quality. The zoom-in panel indicates the double RC feature associated with the Central region.

Conclusions
In this paper we present high-angular-resolution extinction maps (∼ 3 ) for the whole GNS catalogue and de-redden the GNS photometry by discriminating the foreground population. We make the maps and the de-reddened catalogues publicly available to the community.
We identified duplicated sources in the original GNS catalogue and corrected them to produce a clean second data release. We created CMDs with the new data and identified the foreground stellar population via colour cuts in the CMDs H versus J−H and K s versus H−K s . Given that the J band is more affected by extinction than the H and K s bands, we defined a two-layer approach to build extinction maps to de-redden the data in a consistent way. We used RC and red giant stars of similar brightness as a reference for the extinction maps and computed mean intrinsic colours (J − H) 0 = 0.52 ± 0.04 and (H − K s ) 0 = 0.10 ± 0.01, averaging over previous values and computing new ones using a Parsec synthetic stellar population. We obtained mean values for the extinction (A J , A H , and A K s ) for all the regions covered by the GNS, de-reddening the stars belonging to each of the two extinction layers considered (Table 1).
To check the quality of the extinction maps, we compared our results with previous values for the analysed regions and found a good agreement between them. Moreover, we also assessed the differential extinction correction and computed de-reddened CMDs. We found that the scatter of the RC features is significantly reduced (a factor ∼ 2) after applying the extinction maps. Finally, we estimated the statistical and systematic uncertainties for the extinction maps and ended up with values 5 % of the mean extinction in both cases.
The presented extinction maps supersede previous ones for the innermost regions of the Galaxy because they used the GNS survey, which constitutes the highest angular resolution photometric survey in the NIR for the innermost regions of the Milky Way (Nogueras-Lara et al. 2019a). They also comprise the whole GNS survey covering a region (∼ 6000 pc 2 ), which is approximately four times larger than previous extinction maps using the GNS survey.
Given that we computed two A H extinction maps for the first extinction layer (using JH and HK s photometry), we assessed the extinction curve used (Nogueras-Lara et al. 2019b, 2020b and compared the mean values obtained for the extinction A H for each of the extinction maps. They agree within the uncertainties, being consistent with the extinction curve given by the ratios A J /A H = 1.87 ± 0.03 and A H /A K s = 1.84 ± 0.03.
We discussed the limitations of the de-reddening process using the computed extinction maps and, paying special attention to the complexity of translating a 3D process into two 2D extinction layers, and we analysed the possible over-or under-corrections for individual stars. We also discussed non-linear effects affecting the extinction maps, such as the de-reddening of different stellar types and the influence of the effective wavelength when using broad-band filters. We find that these effects are negligible given the uncertainties of the de-reddening process.
To check the effectiveness of the extinction correction, we compared the extinction-map approach with de-reddening on a star-by-star basis (Appendix A). We did this by simulating a stellar population according to the one expected for the NSD (Nogueras-Lara et al. 2020a), also assuming the presence of variable stars. We concluded that the extinction-map approach is more appropriate to de-redden the stars. In particular, the performance of the extinction maps is significantly better when considering the influence of variable stars. Moreover, the extinction maps using RC stars suffer less from systematics effects associated to the a-priori unknown intrinsic colours of each of the de-reddened stars. On the other hand, the extinction maps also allow us to improve the correction for 3D effects because they used more than a single star close in projection.
Finally, we used the de-reddened stars from the first extinction layer to build KLFs and compared the stellar populations and SFHs from the different regions covered by the GNS survey. We concluded that the regions corresponding to the NSD present a significantly different KLF in comparison with the inner bulge regions, pointing towards different stellar populations and SFHs, in agreement with previous work (Launhardt et al. 2002;Nishiyama et al. 2013;Nogueras-Lara et al. 2018b;Gallego-Cano et al. 2020;Nogueras-Lara et al. 2020a).
We accounted for variability considering two simple toy models corresponding to two different scenarios. The first one assumes that ∼ 15 % of the stars are variable following the results of Dong et al. (2017) for a 2.3 × 2.3 region centred on Sagittarius A*. The second one is a more extreme case, ∼ 50 % of variable stars, based on the results obtained by Gautam et al. (2019) for the innermost region of the NSC (∼0.4 pc). We assumed a mean amplitude of H = 0.25 mag for the variable stars (see Fig. 9 in Dong et al. 2017). To account for the variability in K s , we transformed the mean amplitude using Eq. B5 in Feast et al. (2008). This equation was derived for RR Lyraes but given that we considered that the observations for different bands were not taken simultaneously, we did not assume any time correlation between the amplitudes in both bands for each varying star. Thus, it is a way of considering that a variable star in H band is also variable in K s , with an amplitude varying within the limits of the transformation used. In this way, we used two independent Gaussian distributions with standard deviations equal to the mean amplitude for each band to simulate the variability. The described approach is independent of the type of variable stars and does not consider any dependence on the variability period but on the amplitude.
To simulate the extinction and a realistic distribution of the stars following the crowding conditions, we used the extinction map derived in Nogueras-Lara et al. (2020a) and the GNS catalogue. We selected a region of 1.3 × 1.3 and placed the synthetic stars following the positions of real stars from the survey and the extinctions corresponding to the map. Given that too high extinction values are mainly associated with the faint end of the RC features used to create the original extinction map, and thus might be associated with data incompleteness, we removed them to keep a quasi-Gaussian distribution of the extinctions. The uncertainties for the H and K s photometry were simulated considering a Gaussian distribution with a standard deviation of 0.05 mag in agreement with the GNS uncertainties (Nogueras-Lara et al. 2019a).
To simulate the distribution of the stars in the NSD, we assumed a Gaussian distribution centred on 8 kpc (e.g. Gravity Collaboration et al. 2018;Do et al. 2019), with a standard distribution of 150 pc, corresponding to the effective radius of the NSD (e.g. Nishiyama et al. 2013;Gallego-Cano et al. 2020 We kept all the stars with K s < 18.3 mag, corresponding to the 50 % completeness limit of the GNS survey (Nogueras-Lara et al. 2020b). Figure A.1 shows the obtained CMDs (K s vs. H − K s ) for the considered scenarios of variability.

Appendix A.2: Comparison between methods
We used two methods to de-redden the synthetic photometry: -Star by star approach. We used the NICE method (e.g. Lada et al. 1994;Lombardi & Alves 2001) to de-redden each star, following the equation where i and j refer to two NIR photometric bands, A i is the extinction for the i band, and the subindex 0 indicates intrinsic colour. We used j = H and i = K s bands and the ratio A H /A K s = 1.84 (Nogueras-Lara et al. 2020b). The intrinsic colour (H − K s ) 0 is very similar for different stellar types and, in particular, for the red giant stars (see Fig. 29 in Nogueras-Lara et al. 2018a). Moreover, the majority of stars that can be observed in the GC are red giants, given the SFH of the NSD and the extreme source crowding that impedes the observation of main sequence stars (e.g. Nogueras-Lara et al. 2020a, 2021b. Therefore, we assumed a constant intrinsic colour for all the stars (H − K s ) 0 = 0.1 mag (see Fig. 29 in Nogueras-Lara et al. 2018a). A more refined assignation of intrinsic colour for early type star would require spectroscopic observations that are normally restricted to small key regions in the GC given the extreme number of sources.
-Extinction maps. We applied the approach described in Nogueras-Lara et al. (2018a,b, 2020a. We used RC stars as a reference to compute the extinction, given that they are well distributed across the studied fields. To increase the number of reference stars, we also added red giant stars, whose intrinsic colours are similar to RC stars (see Figs. 33 and 34 of Nogueras-Lara et al. 2018a). The extinction for each reference star was computed using Eq. A.1. We determined the extinction at each point by computing an inverse distance The upper and lower panels correspond to the extinction map and the star-by-star approach, respectively. The legend indicates the different temperature ranges considered: T1 = 2500-6500 K, T2 = 6500-10,500 K, and T3 = 10,500-14,500 K. The mean and standard deviations of the distributions are included for each panel.
Article number, page 15 of 17 A&A proofs: manuscript no. Nogueras- Lara   Fig. B.1. Colour-magnitude diagrams before and after the application of the extinction maps for the NSD E, NSD W, T E, and TW. The red and black dots are the original and the de-reddened CMDs, respectively. Only a fraction of stars are plotted given the high number of sources. The labels above each panel indicate the corresponding extinction layer. The green line depicts a Parsec isochrone of ∼ 10 Gyr with twice solar metallicity. The blue and green error bars indicate the systematic uncertainties of the de-reddening and the ZP, respectively.