High-resolution Extinction Map in the Direction of the Strongly Obscured Bulge Fossil Fragment Liller 1^*

The photometric data set used for this work consists in a combination of observations obtained at optical and NIR wavelengths. By using the Wide Field Channel (WFC) of the ACS on board the HST, we secured the first set of high-resolution optical images ever acquired for Liller 1 (GO 15231, PI: Ferraro). The observations have been performed in the filters F606W (V) and F814W (I), and they consist in a set of six deep images per filter (see Table 1). To study Liller 1 at NIR wavelengths we used a data set of high-resolution J and Ks (hereafter K) images acquired with the multiconjugate adaptive-optics assisted camera GeMS/GSAOI at the GEMINI South Telescope (GS-2-13-Q-23, PI: Geisler). These images were already used to obtain a first insight into the stellar population in Liller 1 by Saracino et al. (2015). All the images are dithered by few pixels, in order to avoid spurious effects due to bad pixels. In addition, the NIR data sets were secured following a dither pattern large enough to fill the inter-chip gaps, thus allowing a complete coverage of the innermost ∼100'' × 100'' region of the system.

The data reduction of the ACS-WFC images has been performed on the charge transfer efficiency corrected (flc) images after a correction for pixel-area-map. The photometric analysis was performed via the point-spread function (PSF) fitting method, by using DAOPHOT (Stetson 1987) and following the standard approach used in previous works (e.g., Pallanca et al. 2017; Dalessandro et al. 2018; Cadelano et al. 2019). Spatially variable PSF models were derived for each image by using some dozens of bright and nearly isolated stars. The PSF was modeled adopting an analytical MOFFAT function plus a spatially variable lookup table and it has been applied to all the identified sources with flux peaks at least 3σ above the local background. We then built a master list containing all the stellar-like sources detected in more than three images in at least one filter. In order to improve the level of completeness, we performed an accurate inspection of the residual images, in order to search for additional stellar sources missed during the first iteration of the analysis because of the serious crowding conditions. This procedure allowed us to recover a few dozens of stars, which were added to the master list.

We then forced a fit at the position of each star belonging to the master list, in each frame, by using the ALLFRAME package (Stetson 1987, 1994). For each star, multiple magnitude estimates from different exposures were homogenized by using DAOMATCH and DAOMASTER, and their weighted mean and standard deviation were finally adopted as the star magnitude and photometric error. The final optical HST instrumental catalog contains all the stars measured in at least three images in any of the two filters. For each star, it lists the instrumental coordinates, the mean magnitude in each filter and two quality parameters (chi and sharpness). ¹²

The instrumental optical magnitudes were calibrated onto the VEGAMAG photometric system by using the updated recipes and zero-points available in the HST websites. The instrumental coordinates were first corrected for geometric distortion and then reported to the absolute coordinate system (α, δ) defined by the World Coordinate System by using stars in common with the publicly available Gaia DR2 catalog (Gaia Collaboration et al. 2016a, 2016b). Although only disk field stars have been found in common with the Gaia catalog (due to the very high extinction toward this stellar system), the resulting 1σastrometric accuracy is ∼0.1''. The reduction of the ACS-WFC images has been limited to chip 1 because of the negligible overlap between chip 2 and the FOV covered by the GeMS observations (see Figure 1).

Zoom In Zoom Out Reset image size

At this point, the sample of stars detected in the optical and in the NIR catalogs were compared. From a first analysis it was evident that each catalog contained some objects missed in the other, mainly due to the faintness of the source and/or its vicinity to very bright objects with extreme (blue or red) colors. On the other hand, the analysis also demonstrated that, because of the large extinction in the direction of the system, the V-band exposures are significantly shallower that those secured through the I filter, and, similarly, the J images are much less deep than those obtained through the K filter. This clearly indicated that the optimal solution for the accurate study of the stellar populations in Liller 1 including optical and NIR filters consisted in combining the I and the K images. We thus proceeded with a second-level analysis starting from the independent reduction of the optical data set and the six best-quality (per each filter) NIR images. Two lists of sources have been obtained: one comprising all the stars detected in the I filter and one including the objects identified in the K band. The I-band list contained all the stars measured in at least three frames, over the six available. The K-band list contained all the stars measured in at least three over the six images analyzed. We thus combined the two data sets, building a single master list that contains all the detected stars (some having magnitude measured in both filters, others having only the I or the K measure). The sources detected only in the I images were then searched and analyzed also in the K images, and vice versa, thus maximizing the information provided by the two filters separately. The final magnitudes of this combined analysis were homogenized to the calibrated stand-alone HST and GeMS catalogs, in the Vegamag and Two Micron All Sky Survey systems, respectively. The magnitudes of the brightest stars (K < 12.5), which are saturated in the Gemini observations, have been recovered (up to K ∼ 10) from the PSF photometry of VVV Survey data performed by Mauro et al. (2013). The final catalog reports the I- and K-band magnitudes (together with the V and J ones, if measured) for 43,629 stars detected in both HST and GeMS images.

Figure 1 shows the FOV sampled by chip 1 of the HST/ACS-WFC (outlined in yellow) and by the GeMs observations (dashed blue line). As apparent, the overlapping region, from which we can construct the (I, I − K) CMD and determine the differential reddening map, corresponds to an area of about ∼90'' × 90'' (green area), almost centered onto the gravity center (i.e., determined as the photometric barycenter) of Liller 1 (Saracino et al. 2015). Hereafter, all the discussion and plots are made considering only the stars in this region. Their CMD distribution is shown in three different color and magnitude combinations in Figure 2: from left, to right, the optical (I, V − I) CMD, the hybrid (I, I − K) CMD, and the NIR (K, J − K) one. As expected, extinction along the line of sight limits the deepness of the sample especially in the optical diagram. The differential effect of the reddening is particularly visible in the deformed morphology of the red clump, which appears to be largely stretched along the direction of the reddening vector (from top left to bottom right), especially in the optical and hybrid CMDs. This clearly shows that any meaningful analysis of the stellar populations of Liller 1 requires an appropriate modeling of and correction for the differential reddening effects.

Zoom In Zoom Out Reset image size

For almost all the stars in the catalog, we also measured the relative proper motions by exploiting the time baseline of 6.3 yr separating the available observations. This allowed us to distinguish the stars belonging to Liller 1, from Galactic field interlopers (respectively, black and orange dots in Figure 2). The procedure used for the proper motion measurement and the adopted selection criteria for membership have been presented in F21.

As is well known, interstellar reddening is due to the absorption and scattering of the radiation caused by dust clouds along the light path to the observer. This phenomenon makes the flux of a given source systematically fainter and the color systematically redder than their true (emitted) values, with an efficiency that strongly depends on the wavelength (increasing at shorter wavelengths). The entity of reddening is usually parameterized by the color excess E(B − V), defined as the difference between the observed color (B − V) and the intrinsic color (B − V)₀ that, for stars, is mainly fixed by the temperature of the source.

The extinction law, i.e., the dependence of the absorption coefficient on wavelength (A_λ ) can be expressed as (see also the formulations by Cardelli et al. 1989; Fitzpatrick & Massa 1990; O'Donnell 1994; Fitzpatrick 1999):

$$\begin{eqnarray}&&{A}_{\lambda }={R}_{V}\times {c}_{\lambda ,{R}_{V}}\times E(B-V),\end{eqnarray} \tag{ 1 }$$

where c_λ , R_V = 1 is commonly adopted at the V-band wavelength (${c}_{V,{R}_{V}}=1$). The parameter R_V is usually set equal to 3.1, which is the standard value for diffuse interstellar medium (Schultz & Wiemer 1975; Sneden et al. 1978). This is indeed the value assumed in all previous works devoted to the estimate of reddening in the direction of Liller 1 (see the Introduction). Particular attention, however, has to be given to the extinction toward the inner Galaxy (Popowski 2000; Nataf et al. 2013; Casagrande & VandenBerg 2014; Alonso-García et al. 2017; Pallanca et al. 2021, and references therein), where the R_V value seems to vary along different directions. For example, Nataf et al. (2013) suggested that R_V = 2.5 is more appropriate to reproduce the bulge populations, while some other authors quote larger values (e.g., R_V = 3.2; Bica et al. 2016; Kerber et al. 2019).

Decreasing the value of R_V makes the parameter ${c}_{\lambda ,{R}_{V}}$ (thus, the extinction law A_λ /A_V ) more steeply dependent on wavelength. This is illustrated in the upper panel of Figure 3, where the two curves refer to the canonical R_V value (R_V = 3.1, blue line) and to R_V = 2.5 (red line). As shown in the central panel of the figure, the difference between c_λ,2.5 and c_λ,3.1 is quite small in the photometric bands considered in this work (marked by the dashed vertical lines). However, this difference increases by one order of magnitude when these coefficients are multiplied by their respective value of R_V : this is shown in the bottom panel of the figure, plotting ${R}_{\lambda ,{R}_{V}}\equiv {R}_{V}\times {c}_{\lambda ,{R}_{V}}$. In turn, recalling that $E(B-V)={A}_{\lambda }/{R}_{\lambda ,{R}_{V}}$, this implies that any given value of the absorption coefficient A_λ corresponds to a significantly different color excess E(B − V) if R_V = 2.5 is assumed in place of the standard value R_V = 3.1. Hence, the value of E(B − V) and the distance modulus (μ₀) necessary to superpose an isochrone onto an observed CMD change if R_V = 2.5 is adopted instead of R_V = 3.1. The bottom panel of Figure 3 also shows that the difference between R_{λ, 2.5} and R_{λ, 3.1} is not constant, but significantly increases (in absolute value) for decreasing wavelength in the range relevant for the present study: it varies between ∼0.1 at λ ≃ 21326.7 Å (the K band) and ∼0.6 at λ ≃ 5810.8 Å (the effective wavelength of the F606W filter; Rodrigo et al. 2012; Rodrigo & Solano 2020). ¹³ In turn, the increasing difference between R_{λ, 2.5} and R_{λ, 3.1} for decreasing wavelength necessarily implies that two isochrones, the one shifted to the observed plane by adopting R_V = 3.1, the other shifted under the assumption of R_V = 2.5, cannot simultaneously reproduce the NIR and the optical CMDs of a stellar system. This is illustrated in Figure 4, where the optical (I, V − I), hybrid (I, I − K), and NIR (K, J − K) CMDs of Liller 1 are shown together with an isochrone (Marigo et al. 2017) computed for an age t = 12 Gyr and a global metallicity [M/H] = −0.3. The black dashed line shows the location in the three CMDs of the isochrone obtained by assuming R_V = 3.1, E(B − V) = 3.30, and μ₀ = 14.55, which are the values needed to make it reproduce the observed NIR CMD (Saracino et al. 2015). However, no reasonable match of the observations can be simultaneously obtained in the optical and hybrid CMDs. Conversely, by adopting R_V = 2.5 (red line) all the evolutionary sequences of the cluster in any wavelength combination are well reproduced. This demonstrates that the canonical reddening law is not appropriate in that direction of the sky, and a much lower value of R_V is needed toward Liller 1, in agreement with that suggested in the pioneering work of Frogel et al. (1995) and that previously found for other bulge directions (e.g., Popowski 2000; Nataf et al. 2013; Casagrande & VandenBerg 2014; Alonso-García et al. 2017; Saha et al. 2019; Pallanca et al. 2021).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The manifest inadequacy of R_V = 3.1 in the direction of Liller 1 did not come to light in previous works because they were all based on NIR or optical data sets separately: it is indeed the simultaneous analysis of these CMDs that finally revealed the necessity of a much lower value of R_V . Of course, however, a proper comparison between theoretical isochrones and observations first requires to correct the CMDs for the effects of differential reddening, which is responsible for the distortions of the evolutionary sequences clearly visible in Figure 4.

Our result (i.e., R_V = 2.5) is perfectly in agreement with previous works aimed to estimate the reddening law in the bulge direction (Nataf et al. 2013; Saha et al. 2019). For example, R_V = 2.5 perfectly corresponds to the peak of the distribution shown in Figure 8 of Nataf et al. (2013): indeed, as explicitly discussed in the conclusions of that paper, the mean value of the distribution, A_I /E(V − I) ∼ 1.217, corresponds to R_V = 2.5. On the other hand, we found that the absorption coefficients derived by Saha et al. (2019) in the direction of the Baade's window and shown as black circles in their Figure 8 are well reproduced by the O'Donnell (1994) extinction law calculated for R_V ∼ 2.0, and the values obtained for R_V = 2.5 place in the middle between these points and those corresponding to the standard R_V = 3.1 (red circles in their figure). Note that a value as small as R_V = 2.0, even if surprising, is still in agreement with the left-hand tail of the distribution shown in Figure 8 of Nataf et al. (2013). Hence, our result is perfectly in agreement with previous works and it brings further evidence that the commonly adopted R_V = 3.1 is likely inappropriate toward the bulge, while a precise estimate of the extinction law in each specific direction is required (see also Pallanca et al. 2021).

The interstellar extinction along the Galactic plane and in the direction of the Galactic center can be highly spatially variable (differential reddening), on scales that can be as small as a few arcseconds. Since this spatial dimension is significantly smaller than the apparent size of most Galactic stellar systems, the CMD of star clusters affected by differential reddening shows a systematic elongation of all the evolutionary sequences along the direction of the reddening vector. In turn, this prevents the accurate characterization of their photometric properties, with a non negligible impact on the determination of fundamental parameters such as, for instance, the distance and age (e.g., Bonatto et al. 2013). To model and correct for the effect of differential reddening in GCs, two main approaches have been adopted in the literature. Schematically: (1) the "cell by cell" method, which consists in dividing the observed FOV in a number of cells, build the CMD of the stars included in each cell, and calculate the reddening value in each cell from the color and magnitude displacements with respect to the bluest CMD (Heitsch & Richtler 1999; Piotto et al. 1999; von Braun & Mateo 2001; McWilliam & Zoccali 2010; Nataf et al. 2010; Gonzalez et al. 2011, 2012; Massari et al. 2012; Bonatto et al. 2013; Saracino et al. 2015), and (2) the star-by-star approach, which consists in estimating the reddening of each star in the catalog from the relative displacement between its local CMD, i.e., the CMD built with the spatially closest objects, and the cluster mean ridge line (Milone et al. 2012; Bellini et al. 2013; Pallanca et al. 2019; Saracino et al. 2019; Cadelano et al. 2020). In the cited papers, different authors used stars in different evolutionary sequences as reference for the determination of the color and magnitude displacements: MS, horizontal branch, red giant branch (RGB) stars, and even variable RR Lyrae (Alonso-García et al. 2011, and references therein). In all cases, the spatial resolution is driven by the need for a statistically reliable sample of stars to build the cell or local CMDs. To determine the differential extinction map in the direction of Liller 1 we used the star-by-star method, adopting the RGB, subgiant branch (SGB) and bright MS stars as reference, as already tested in a previous work (Pallanca et al. 2019). In addition, because of the results discussed in the previous section, we adopted R_V = 2.5.

The first step consists in the definition of the reference mean ridge line (MRL) of the OP of Liller 1. This has been determined in the proper motion cleaned, hybrid (I, I − K) CMD, which is significantly deeper than any other available color-combination (even more than the (K, J − K) CMD, which is likewise deep but with a worst signal-to-noise ratio in the J filter at the bottom of the MS; see Figure 2), thus allowing us to also use stars along the SGB and MS phases. In turn, this implies a significant increase of the spatial resolution of the output reddening map, with respect to the case where only RGB stars could be used. In order to remove poorly measured objects that can bias the procedure, we first considered only stars with high-quality photometry, by adopting 3σ cuts in the chi and sharpness parameter distributions. The considered sample, comprising the RGB, SGB, and the upper MS of the OP, was divided in magnitude bins of variable size (ranging from 0.25–1 mag), in order to best sample the morphological changes of the sequence in the MS turnoff region. The mean color of each considered bin was then determined after a 3σ-clipping rejection, in order to remove any residual interlopers.

Once the MRL is defined, we considered all the stars in the catalog (i.e., not only the high-quality objects) and for each of them we selected, from the high-quality pool (and proper motion selected) sample, the N_* sources spatially lying in the vicinity of the star under investigation. The method consists in determining the value of δE(B − V) required to shift the MRL of the system along the direction of the reddening vector until it fits the local CMD defined by these N_* nearby sources. Of course, the larger the number of sources (N_*) is, the better defined the local CMD is. However, to keep the spatial resolution as high as possible, we constrained N_* below a maximum value and we assumed a limiting radial distance (r_lim) within which the nearby sources can be selected. We set N_*,max = 30 and r_lim = 3''. As soon as one of these two conditions is satisfied, the other one is no longer considered. Since the stellar density is higher in the center and decreases with radius, this approach naturally yields an adaptive grid producing reddening maps with higher (lower) spatial resolution in the central (external) regions.

As described in F21, the subsample with measured proper motion (PM) represents about 80% of the entire sample. In particular due to dithering patterns adopted in the GeMS/GSAOI observations (aimed at covering the inter-chip gaps), a sort of external corona and a central cross of the FOV was not properly sampled by all the considered NIR exposures to allow accurate PM measures. Thus, the procedure aimed at the construction of the reddening map has been applied twice. The first time, only stars classified as Liller 1 members from their measured proper motion were included in the high-quality pool. The second time, in order to cover the entire FOV, no proper motion selection was used. The differential reddening estimates obtained from these two approaches are fully compatible (as already found for NGC 6440 in Pallanca et al. 2019), demonstrating that the potential presence of a few field interlopers does not affect the result. Hence, the final differential reddening value assigned to each star has been obtained, whenever possible, from the (formally more rigorous) proper motion-based method, and complemented with the value derived neglecting the proper motion selection in the other cases.

Because of the large central densities, the number of stars within a circle of 3'' radius is typically much larger than 30 in the inner cluster regions; hence, the N_*,max closest sources are typically selected within a distance much smaller than r_lim from each star under investigation in this region. The opposite is true in the outskirts, where the driving parameter for the selection of nearby sources becomes r_lim. This is clearly visible in the left panel of Figure 5: the distance used for the selection of the nearby stars defining the local CMD can be as small as 1'' in the central regions, while it increases up to the maximum value (r_lim = 3'') close to the boundaries of the analyzed FOV, where only semicircles toward the center can be drawn. The few spots visible on the map correspond to local increases of the selection radius due to the presence of saturated stars. The discontinuities in the innermost spot and at the outermost edge of the FOV are due to the transition between reddening estimates that take into account/neglect the proper motion condition (which is unavailable in these regions: by keeping larger numbers of stars, this allows us to use a slightly smaller selection radius).

Zoom In Zoom Out Reset image size

The comparison between the MRL of the cluster OP and the local CMDs has been performed in the magnitude range 19.5 < I < 26, in order to discard potentially saturated and too faint (hence not well measured) objects. We also performed a 3σ rejection in order to exclude from the analysis all the stars with a color distance from the MRL significantly larger than that of the bulk of the selected objects. This is useful to exclude field interlopers not removed by the proper motion analysis and/or discard noncanonical stars with intrinsic colors different from those of the cluster main populations (Pallanca et al. 2010, 2013, 2014; Dalessandro et al. 2014; Cadelano et al. 2015; Ferraro et al. 2015, 2016a, 2018).

The comparison was performed by shifting the MRL along the reddening direction (which is fixed by the adopted extinction law) in steps of δE(B − V). For each considered step, we calculated the residual color Δ IK as

$$\begin{eqnarray}&&{\rm{\Delta }}{IK}=\sum _{i=1}^{{N}_{* }}\left(\left|{{IK}}_{\mathrm{obs},i}-{{IK}}_{\mathrm{MRL},i}\right|+{w}_{i}\cdot \left|{{IK}}_{\mathrm{obs},i}-{{IK}}_{\mathrm{MRL},i}\right|\right),\end{eqnarray} \tag{ 2 }$$

where IK_obs,i is the (I − K) color observed for each of the N_* selected objects constituting the local CMD, IK_MRL,i is that of the MRL at the same level of magnitude, and the weight w_i takes into account both the photometric error on the color (σ) and the spatial distance (d) of the ith source from the star under investigation, according to the following expression:

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{1}{{d}_{i}\cdot {\sigma }_{i}}{\left[\sum _{j=1}^{{N}_{* }}\left(\displaystyle \frac{1}{{d}_{j}\cdot {\sigma }_{j}}\right)\right]}^{-1}.\end{eqnarray} \tag{ 3 }$$

The second term is meant to give more weight to the closest stars, thus increasing the accuracy of the local reddening estimate, reducing the impact of using large selection distances where the level of dust absorption could change. Note that we cannot reduce Equation (2) to the second term only because a bright object (with small photometric errors) close to the central star would dominate the weight and make the reddening estimate essentially equal to the value needed to move the MRL exactly on that object. The final differential reddening value assigned to each star is the value of the step δE(B − V) that minimizes the normalized residual color (Δ IK/ N_*). The spatial distribution of differential extinction values thus derived for each star constitutes the differential reddening map shown in the right panel of Figure 5. We emphasize that in the adopted procedure, the value of δE(B − V) assigned to each star is not guided by its position in the CMD with respect to the cluster MRL, but by the position of its N_* neighboring stars (i.e., the N_* stars spatially adjacent on the plane of the sky to the object under analysis). This ensures that we are correcting for a phenomenon (differential reddening) affecting the region of sky where the star is located, not for possible intrinsic properties of the star itself (e.g., its chemical abundance).