Digging for Relics of the Past: The Ancient and Obscured Bulge Globular Cluster NGC 6256

The differential reddening-corrected CMD discussed above has been used to determine the cluster properties (distance, absolute color excess, and stellar population age) via isochrone fitting. To this aim, it is helpful to obtain independent estimates of both the cluster distance modulus and color excess to be used as a starting point in the isochrone fitting procedure. As discussed in Section 1, the values for these parameters available in the literature are still extremely uncertain.

Unfortunately, no RRLyrae stars (that would be suitable to determine the distance modulus) are known in this cluster. Moreover, the publicly available parallaxes in the Gaia Data Release 2 (DR2) catalog (Gaia Collaboration et al. 2018) are extremely uncertain for the stars in the direction of this cluster and they cannot be used to constrain the cluster distance. The same holds for the absolute color excess: the spatial resolution of the publicly available all-sky maps of interstellar extinction (e.g., Schlafly & Finkbeiner 2011) are too poor to get a reliable average value of this parameter, given the extreme variations on very small scales shown in Figure 3.

We thus estimated the cluster distance modulus and average color excess by comparing the observed CMD with a catalog of stars of the GC NGC 6752, obtained from images acquired under GO 11904 (PI: Kalirai) with the same instrument and combination of filters as for NGC 6256. The data reduction and calibration of these images are exactly the same as those described in Section 2.

NGC 6752 is a ∼13 Gyr old system (Gratton et al. 2003; Dotter et al. 2010), with a metallicity [Fe/H] = −1.56 ± 0.04 (Carretta et al. 2009) comparable to that of NGC 6256, and very well constrained distance modulus and color excess: (m − M)₀ = 13.14 ± 0.06 and E(B − V) = 0.04 ± 0.01 (Ferraro et al. 1999; Baumgardt et al. 2019). For two systems with approximately the same metallicity and age, it is expected that any magnitude difference between their MRLs is mainly due to the relative difference in distance and color excess. The left panel of Figure 4 shows that, assuming the metallicity, distance modulus, and color excess quoted above for NGC 6752, a 13 Gyr isochrone extracted from the Dartmouth Stellar Evolutionary Database (DSED; Dotter et al. 2008) provides a good match of the cluster CMD. We then applied color and magnitude shifts to the MRL of NGC 6256 obtained in Section 3 (cyan dashed line in the right panel of Figure 4) to superimpose it onto the MRL of NGC 6752 (created using the same procedure described in Section 3). To do this, we performed a fit of the two curves in the magnitude range $15\lt {m}_{{\rm{F}}814{\rm{W}}}\lt 18$ and chose the best-fit shifts on the basis of the χ² statistics, finding ${\rm{\Delta }}({m}_{{\rm{F}}555{\rm{W}}}-{m}_{{\rm{F}}814{\rm{W}}})=-1.56$ and ${\rm{\Delta }}{m}_{{\rm{F}}814{\rm{W}}}=-3.3$. Adopting the same extinction coefficients used in Section 3, these shifts correspond to a color excess E(B − V) = 1.18 ± 0.05 and a distance modulus (m − M)₀ = 14.25 ± 0.1 for NGC 6256. We stress that these are only first-guess estimates, just used as starting points for the accurate isochrone fitting of the cluster CMD discussed in the next section. In fact, these relative measurements could be biased by the possible differences in the cluster's relative age and in the light-element abundances of the cluster's subpopulations. Moreover, these results could also be biased by the effects of a very different E(B − V) on the observed sequences (see Section 4.2 and Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To derive the absolute age of the cluster via isochrone fitting, we followed a Bayesian procedure similar to that used by Saracino et al. (2019, see also Correnti et al. 2016; Kerber et al. 2018, 2019; Cadelano et al. 2019). This approach allows us to estimate the stellar system age through a one-to-one comparison between the observed CMD and a set of theoretical models, simultaneously exploring reasonable grids of values for the relevant cluster parameters (not only the age, but also the distance modulus and color excess).

The observed CMD was compared to three different sets of α-enhanced isochrones: the DSED models (Dotter et al. 2008), the Victoria–Regina Isochrone Database (VR, VandenBerg et al. 2014), and the BaSTI stellar evolution models (Pietrinferni et al. 2004, 2006). For each family of models we assumed [α/Fe] = +0.4, which is the typical value for bulge GCs, and a standard He abundance Y = 0.25. In the case of the DSED and VR databases, isochrones can be generated at different metallicities, while for the BaSTI database only a selection of [Fe/H] values is available. Therefore, in order to perform a proper comparison between the results obtained from these different models, we assumed in all cases [Fe/H] = −1.62, a value that is extremely close to that derived through spectroscopy (Vásquez et al. 2018) and that is available for all three sets of models. The adopted extinction coefficients in the F555W and F814W bands take into account the dependence of the reddening on the effective temperature of the stars, following the extinction laws of Cardelli et al. (1989), the equations in Girardi et al. (2002) and assuming the alpha-enhanced spectral library used in Pietrinferni et al. (2006). This is indeed necessary for performing reliable isochrone fitting in the case of highly reddened systems, because neglecting such a temperature dependence makes the red giant branch and the main sequence shallower, thus increasing the overall color range spanned by the isochrone in the CMD. This is shown in Figure 5, where we plot, as an example, the same DSED isochrone (with an age of 12.5 Gyr, [Fe/H] = −1.62 and the first-guess values of distance modulus and E(B − V) estimated for NGC 6256 in the previous section) before and after such a correction: the black line corresponds to the isochrone obtained by using constant extinction coefficients (R_F555W = 3.207 and R_F814W = 1.842), while the red line is derived by correcting the extinction coefficients for the dependence on the surface temperature of the stars. As can be seen, the significant color difference between the two curves demonstrates that reliable age estimates for heavily reddened systems, such as NGC 6256, cannot be performed by assuming constant extinction coefficients, especially when dealing with optical filters (in fact, this effect is larger for decreasing wavelengths).

To compare the observed CMD of NGC 6256 with the three families of stellar models we applied a Markov Chain Monte Carlo (MCMC) sampling technique. We assumed that the probability of an nth star to belong to an isochrone (corrected for a given reddening and distance) can be expressed in terms of a Gaussian distribution:

$$\begin{eqnarray}&&{p}_{n}=\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{n}}\exp \left[-\displaystyle \frac{1}{2}\displaystyle \frac{{d}_{n}^{2}}{{\sigma }_{n}^{2}}\right]\propto \displaystyle \frac{\exp \left(-\tfrac{{d}_{n}^{2}}{2{\sigma }_{n}^{2}}\right)}{{\sigma }_{n}},\end{eqnarray} \tag{ 2 }$$

where d_n is the minimum distance of the star from the isochrone and σ_n the uncertainty of such a distance expressed as $\sigma =\sqrt{{\sigma }_{\mathrm{col}}^{2}+{\sigma }_{\mathrm{mag}}^{2}+{\sigma }_{\mathrm{DR}}^{2}}$, where σ_col, σ_mag are the photometric uncertainties on the color and F814W magnitude, respectively, while σ_DR are the uncertainties on the differential correction derived as discussed in Section 3. Therefore, the logarithmic likelihood function of a given isochrone can be expressed as the logarithm of the sum of p_n over all the nth stars:

$$\begin{eqnarray}&&\mathrm{ln}{\mathscr{L}}=\mathrm{ln}\sum _{n=1}^{{N}_{\mathrm{stars}}}{p}_{n}\propto -\displaystyle \frac{1}{2}\sum _{n=1}^{{N}_{\mathrm{stars}}}\left[\displaystyle \frac{{d}_{n}^{2}}{{\sigma }_{n}^{2}}+\mathrm{ln}({\sigma }_{n}^{2})\right].\end{eqnarray} \tag{ 3 }$$

To sample the posterior probability distribution in the n-dimensional parameter space, we used the emcee code (Foreman-Mackey et al. 2019). Since no age estimates are available in the literature for NGC 6256, we explored a wide range of old ages, uniformly distributed from 10 to 15 Gyr, in steps of 0.25 Gyr. To put the isochrones into the observational plane, we used values of the color excess and distance modulus following Gaussian prior distributions peaked at E(B − V) = 1.2 ± 0.1 and (m − M)₀ = 14.25 ± 0.1, respectively. In order to minimize the contamination by field interlopers, we extended these calculations only to stars within the cluster half-light radius (∼50''; Harris 2010). Moreover, we restricted the analysis only to stars in the magnitude range $18.5\lt {m}_{{\rm{F}}814{\rm{W}}}\lt 20.5$, a CMD region surrounding the turn-off and therefore sensitive to stellar age variations.

The results obtained in terms of age, distance modulus, and color excess are shown in Figures 6 and 7 for the three adopted sets of theoretical models. In each figure, the left-hand panel shows the $({m}_{{\rm{F}}814{\rm{W}}},{m}_{{\rm{F}}555{\rm{W}}}-{m}_{{\rm{F}}814{\rm{W}}})$ CMD and the best-fit isochrones. In all the cases, the best-fit model very well reproduces the cluster evolutionary sequences. The one- and two-dimensional posterior probabilities for all of the parameter combinations are presented in the right-hand panels as corner plots. The best-fit values and their uncertainties (based on the 16th, 50th, and 84th percentiles) obtained for the age, E(B − V) and (m − M)₀ from the three sets of models are also listed in Table 1. Since we are using different models (adopting slightly different solar abundances, opacities, reaction rates, efficiency of atomic diffusion, etc.), the resulting best-fit values are somewhat expected to be not exactly the same. However, they are all mutually consistent within the errors. NGC 6256 turns out to be a very old cluster, with an age of 13.0 Gyr (obtained as the average of the three derived values). The average color excess, E(B − V) = 1.19, is in good agreement with that quoted by Valenti et al. (2007), and it implies that, within the field of view covered by our observations, this parameter varies from a minimum of ∼1.0 up to ∼1.5. The obtained distance modulus corresponds to a distance of 6.8 ± 0.1 kpc from the Sun and turns out to be significantly different from some of the values quoted in the literature. Indeed, the cluster results in a value 3.5 kpc closer than that reported in Harris (2010) and 2.3 kpc closer than the value quoted in Valenti et al. (2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a consistency check, we verified that the results remain basically unchanged if a uniform prior spanning a large interval of values for both the color excess and the distance modulus is assumed. As an additional test, we repeated the isochrone fitting procedure including the cluster metallicity as a fit parameter. This was, however, feasible only for the DSED and the VR models. We allowed the metallicity of the cluster to vary from −1.1 to −1.8 dex, following a normal distribution with peak value [Fe/H] = −1.61 and standard deviation equal to 0.2, as derived by Vásquez et al. (2018). The results in terms of age, color excess, and distance modulus are basically unchanged with respect to those quoted in Table 1, while the best-fit metallicity values are [Fe/H] = −1.6 ± 0.1 for both the DSED and VR models. They are in excellent agreement with the values derived through spectroscopy.

Throughout the whole analysis, we assumed that the cluster stellar population has a standard He content Y = 0.25. The presence of a He-enhanced subpopulation of stars could introduce a broadening of the sequences, which should be visible in optical CMDs. We compared the observed width of the evolutionary sequences in the CMD with that derived from a synthetic stellar population having a standard He content. This was created generating stars with F555W and F814W magnitudes randomly extracted from the best-fit BaSTI model convoluted with the observed error distribution in each filter: $\sigma =\sqrt{{\sigma }_{\mathrm{mag}}^{2}+{\sigma }_{\mathrm{DR}}^{2}}$. We found that the width of the synthetic CMD is comparable to the observed one along all the evolutionary sequences, thus confirming that the possible presence of a He-enhanced subpopulation cannot be assessed by using the current data set, as the sequence broadening is dominated by the residuals of the differential reddening corrections. While for a cluster with a mass around 1 × 10⁵M_⊙, such as NGC 6256 (Baumgardt & Hilker 2018), the presence of subpopulations is expected with a maximum He spread δY_max = 0.01–0.02 (Milone et al. 2018), we repeated the MCMC analysis using BaSTI models for a stellar population exclusively composed of stars with Y = 0.3 (δY = 0.05). We found that the derived age is consistent within the uncertainties with those derived previously (Table 1). Therefore we conclude that, with the available photometry, variations of He abundances expected for this clusters do not have significant impact on the results presented here.

The gravitational center of the cluster was determined following an iterative procedure based on the position of resolved stars and described in Lanzoni et al. (2010, see also Montegriffo et al. 1995; Cadelano et al. 2017; Lanzoni et al. 2019). The procedure starts by determining the distance from a first-guess center of all the stars included within a circle of radius r centered on it, and then adopts as new center the average of the star coordinates. The procedure is then iteratively repeated: during each iteration the starting center is that derived in the previous iteration and the procedure stops when convergence is reached (i.e., when the newly determined center differs by less than $0\buildrel{\prime\prime}\over{.} 01$ from the 10 previous determinations). We performed this procedure several times, adopting different values of r and selecting stars in different magnitude ranges, chosen as a compromise between having high enough statistics and avoiding spurious effects due to incompleteness and saturation. In particular, the radius r was chosen in the range from 20'' to 30'', spaced by 5''. For each radius r, we explored different magnitude ranges, from ${m}_{{\rm{F}}814{\rm{W}}}\gt 16.0$ (in order to exclude stars close to the saturation limit), down to m_F814W = 20–21, in steps of 0.5 mag. We also excluded from the procedure stars with color $2\lt ({m}_{{\rm{F}}555{\rm{W}}}-{m}_{{\rm{F}}814{\rm{W}}})\lt 3$, in order to minimize the contamination due to field interlopers. As a first-guess center, we used the value quoted by Harris (2010). The final coordinates adopted for the cluster gravitational center are the mean of the different values obtained in the different runs, and the uncertainty is their standard deviation of the measures: R.A. = 16^h59^m32668 and decl. = −37°7'15139, with an uncertainty of about 04. Our newly determined center differs by ∼2'' from that quoted in Harris (2010).

To compute the bulk absolute motion of NGC 6256 on the plane of the sky, we made use of the absolute PMs of individual stars available in the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018). First, we selected all the stars in common between our HST catalog and the Gaia DR2 data set having an absolute PM measure. Then, we rejected all the objects with poorly measured PM according to the prescriptions given in Arenou et al. (2018). The CMD and the vector-point diagram (VPD) of all the stars in common between our observations and the Gaia catalog are shown in Figure 8. We further refined the selection by considering only those stars included within a circle centered on the overdensity of points clearly visible in the VPD, whose center and radius were derived as the mean and standard deviations of the best-fit Gaussians of the histograms shown in the right panel of Figure 8 (${\mu }_{\alpha }\cos \delta =3.6\pm 1.0$ mas yr⁻¹, μ_δ ≈ −1.6 ± 1.0 mas yr⁻¹). These stars are plotted as red dots in the CMD of Figure 8). The absolute PM of NGC 6526 was finally measured from this fiducial sample of 189 stars, by using the Gaussian maximum-likelihood method described in Walker et al. (2006, see also Equation (3) of Pryor & Meylan 1993) stars. We found ${\mu }_{\alpha }\cos (\delta )=-3.7\pm 0.2$ mas yr⁻¹ and μ_δ = −1.6 ± 0.2 mas yr⁻¹, in the J2000.0 system. These values are consistent with those determined by Vasiliev (2019) and Baumgardt et al. (2019), which are both based on Gaia DR2 PMs as well.

Zoom In Zoom Out Reset image size

The cluster absolute PM thus derived, combined with the radial systemic velocity from Vásquez et al. (2018, v_r = −103.4 ± 0.5 km s⁻¹), was used to back-integrate the cluster orbit within the Milky Way potential well. To this purpose, we used GravPot16⁶ (Fernández-Trincado et al. 2017), a package that generates stellar orbits in a semianalytic, steady-state, and three-dimensional Galaxy model, based on the gravitational potential derived from the Besançon Galactic Model (Robin et al. 2003) and including also a prolate bar (Robin et al. 2012). This tool has been recently used to compute the orbits of several Galactic star clusters (Gaia Collaboration et al. 2018; Libralato et al. 2018; Bellazzini et al. 2019). The Galaxy parameters (such as the bar properties, Sun distance from the Galactic center, etc.) were set to the same standard values adopted by Gaia Collaboration et al. (2018). Figure 9 shows the resulting orbit of the cluster during the last gigayear. We verified that nothing significantly changes if the cluster initial conditions are randomly varied within the uncertainty ranges of its position and velocity. It can be seen that NGC 6256 is in a low-eccentricity orbit (e ∼ 0.11), strongly confined within the Galactic bulge. In fact, during its ∼50 Myr orbital period, it reaches a minimum and maximum distance from the Galactic center of about 0.8 kpc and 2.9 kpc, respectively. This confirms that, despite its relatively low metallicity, this system is not a halo intruder crossing the bulge during a fraction of its orbit (as suggested for other bulge GCs, like NGC 6681, Terzan 10 and Djorgovski 1, Massari et al. 2013; Ortolani et al. 2019), but it is a GC genuinely belonging to the bulge population. As a consistency check, we compared our results with those obtained using different orbit integrators (see, e.g., Cadelano et al. 2017) that adopt different shapes of the Galactic potential. We found that all the results are consistent within the uncertainties. These results are also qualitatively in agreement with those reported by Baumgardt et al. (2019) and Perez-Villegas et al. (2020). It is worth stressing that the simulated cluster orbit cannot be used to firmly asses if the cluster birthplace was indeed the Galactic bulge, since a back-integration of the orbit for a time as long as the cluster age should take into account the (basically unknown) variation of the Galactic potential as a function of time.

Zoom In Zoom Out Reset image size

Massari et al. (2019) calculated the integrals of motion for a large sample of Galactic GCs to infer their birthplaces (i.e., to infer which clusters belong to the main Galactic population and which ones have been likely accreted during merger events that the Galaxy experienced in the past). To this aim they used the cluster current positions from Harris (2010) and kinematics from the Gaia DR2 catalog. They found that NGC 6256 is unlikely to belong to the Galaxy in situ population and it is most likely associated with an accreted low-energy component (see Figure 3 in Massari et al. 2019). However, their result is based on a distance value significantly different from the one computed here, while the cluster kinematic is in agreement with the results derived here and by Vásquez et al. (2018). We therefore repeated the same procedure followed by Massari et al. (2019) using the updated cluster position and found that NGC 6256 has an orbital energy of E ∼ −2.3 × 10⁵ km² s⁻² and an angular momentum of L_z ∼ −300 km s⁻¹ kpc. These updated values place NGC 6256 in a region of the integrals of motion space occupied by systems that genuinely belong to the population of bulge clusters formed in situ (Figure 10).

Zoom In Zoom Out Reset image size

These dynamical arguments strongly suggest that NGC 6256 is an old and low-metallicity GC that was formed in the bulge during the very early epoch of the Galaxy assembly.