The ACS Survey of Galactic Globular Clusters XIV: Bayesian Single-Population Analysis of 69 Globular Clusters

TL;DR

This study employs Bayesian methods to analyze HST imaging data, deriving ages, distances, and chemical compositions for 69 globular clusters, providing insights into their formation and evolution.

Contribution

It introduces a consistent Bayesian framework for single-population analysis of globular clusters, improving parameter estimation accuracy and reproducibility.

Findings

Identified correlations between helium abundance and metallicity.

Revised age-metallicity relation for globular clusters.

Found links between helium and cluster properties like binary fraction.

Abstract

We use Hubble Space Telescope (HST) imaging from the ACS Treasury Survey to determine fits for single population isochrones of 69 Galactic globular clusters. Using robust Bayesian analysis techniques, we simultaneously determine ages, distances, absorptions, and helium values for each cluster under the scenario of a "single" stellar population on model grids with solar ratio heavy element abundances. The set of cluster parameters is determined in a consistent and reproducible manner for all clusters using the Bayesian analysis suite BASE-9. Our results are used to re-visit the age-metallicity relation. We find correlations with helium and several other parameters such as metallicity, binary fraction, and proxies for cluster mass. The helium abundances of the clusters are also considered in the context of CNO abundances and the multiple population scenario.

Tables

Table 1. Table 1: BASE-9 Model

Notation Symbol	Description
$M_i$	Primary initial stellar mass
$R_i$	Ratio of secondary to primary initial stellar masses
$Z_i$	Cluster membership indicator
$\Theta$	Cluster-level parameters
	(${\theta }_{\mathrm{age}}$, ${\theta }_{\left[Fe/H\right]}$, ${\theta }_{m-M_V}$, ${\theta }_{A_V}$)
${\Phi }_Y$	Helium abundance of the cluster
Xi	Photometry for $N$ stars in $n$ filters
${𝚺}_i$	Photometric uncertainties for $N$ stars in $n$ filters

Table 2. Table 5: Primordial Helium Values from Literature

Primordial Helium	Reference
0.232${}_{-0.047}{}^{}{}_{}{}^{+0.044}$ to 0.293${}_{--0.048}{}^{}{}_{}{}^{+0.046}$	Planck Collaboration et al. (2014)
0.254$\pm 0.003$	Izotov, Stasińska & Guseva (2013)
0.278${}_{-0.032}{}^{}{}_{}{}^{+0.034}$	Hinshaw et al. (2013)
0.2485$\pm$0.0002	Aver et al. (2013)
0.2477$\pm$0.0029	Peimbert, Luridiana & Peimbert (2007)
0.250$\pm$0.006	Salaris et al. (2004)

Equations

$$\mathbf{\mbox{\boldmath{$\mu$}}_{i}}=-2.5\log_{10}(10^{-\mbox{\boldmath{$G$}}(M_{i},\mbox{\boldmath{$\Theta$}},\phi_{Y})/2.5}$$

$$L(\mathbf{M},\mathbf{R},\mathbf{Z},\mathbf{\Theta},\phi_{Y}\lvert\mathbf{X},\mathbf{\Sigma})=\prod_{i=1}^{N}\left[Z_{i}\times\frac{1}{\sqrt{(2\pi)^{n}\left\lvert\mathbf{\Sigma}_{i}\right\lvert}}exp\left(-\frac{1}{2}(\mathbf{X}_{i}-\mbox{\boldmath{$\mu$}}_{i})^{\top}\times\mathbf{\Sigma}_{i}^{-1}(\mathbf{X}_{i}-\mbox{\boldmath{$\mu$}}_{i})\right)+(1-Z_{i})\times P(\mathbf{X}_{i}\lvert Z_{i}=0)\right].$$

$$L(\mathbf{M},\mathbf{R},\mathbf{Z},\mathbf{\Theta },{\varphi }_Y|\mathbf{X},\mathbf{\Sigma })=\prod _{i=1}^N\left[Z_i\times P({\mathbf{X}}_i|{\mathbf{\Sigma }}_i,M_i,R_i,\mathbf{\Theta },{\varphi }_Y,Z_i=1)+(1-Zi)\times P({\mathbf{X}}_i|Z_i=0)\right],$$

$$M_{RGBB}=0.161(\pm 0.030)[Fe/H{]}^2+$$

$$1.240(\pm 0.088)[Fe/H]+1.917(\pm 0.064)$$

$$M_{RGBB}=-0.080(\pm 0.014)(Y\times 100)+3.120(\pm 0.452)$$

$$\Delta Log(P)=-0.072(\pm 0.024)[Fe/H{]}^2$$

$$-0.308(\pm 0.080)[Fe/H]-0.286(\pm 0.065)$$

$$Resid=0.085(\pm 0.143)Y-0.027(\pm 0.046)$$

Full text

Keywords: (Galaxy:) globular clusters: general, (stars:) Hertzsprung-Russell and colour-magnitude diagrams, Galaxy: formation

††pagerange: The ACS Survey of Galactic Globular Clusters XIV: Bayesian Single-Population Analysis of 69 Globular Clusters–The ACS Survey of Galactic Globular Clusters XIV: Bayesian Single-Population Analysis of 69 Globular Clusters††pubyear: 2014

1 Introduction

Although it is clear that most, if not all, globular clusters harbor two or more populations of stars, these multiple populations do not obviously affect the color-magnitude diagram (CMD) at the visual wavelengths. This is true for many clusters even with CMDs of high photometric quality, such as those obtained from HST imaging. After decades, the visual CMD is well-studied and well-modeled theoretically. Multiple population scenarios agree on a mutual distance and absorption for all populations of stars within the cluster, and most suggest little to no spread in the iron abundances for most clusters. Internal differences in age are also expected to be small for the majority of clusters. Therefore, value remains in using the visual filters F606W and F814W to determine single population CMD-based estimates of various cluster parameters, especially to obtain relative ages, chemical composition, and distances to a high degree of accuracy and precision. This will continue to be true until theoretical models at ultraviolet wavelengths are able to model clusters with similar quality as the visual wavelengths.

For decades, CMDs have been fit with isochrones “by eye”. Occasionally, various numerical techniques (including Bayesian approaches) have been implemented to attempt a determination of the best isochrone (Andreuzzi et al. 2011; Jørgensen & Lindegren 2005; Naylor & Jeffries 2006; Hernandez & Valls-Gabaud 2008; Janes et al. 2013; Valls-Gabaud 2014). However, uncertainties in distances, reddenings, and metallicities continue to complicate isochrone fitting attempts, making absolute age measurements of Galactic globular clusters (GGCs) difficult. Many previous studies have focused on comparing the relative ages, which are still capable of increasing our understanding of the formation and evolution of the Galaxy. Studies often focus on measuring differences in color or magnitude from key points in the CMD, such as the main sequence turnoff, red giant branch, and zero age horizontal branch (Gratton 1985; Sarajedini & King 1989; VandenBerg 2000). An alternative method by Marín-Franch et al. (2009) used main sequence fitting of deep photometry from Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) to determine relative ages and distances.

In this paper, we focus on Bayesian analysis, which incorporates many of the above insights as well as a broad range of what we know about globular clusters. We employ principled probability-based methods that embed physics-based models into an overall statistical model that accounts for the complexity of the astronomical sources, thereby leading to greater reliability in practice (van Dyk et al. 2009; De Gennaro et al. 2009; Stein et al. 2013). With the computational power available in the modern era, implementing a Bayesian technique is preferable. Besides being more objective than the chi-by-eye approach, Bayesian methods allow simultaneous fitting of multiple characteristics of a cluster along with principled measures of uncertainty (van Dyk et al. 2009; Jeffery et al. 2016). Bayesian methods also can produce the joint posterior distribution for two or more parameters, allowing us to fully characterize complex non-linear correlations among the fitted parameters. We use a large grid of theoretical models and simultaneous fitting of multiple filters and multiple parameters to obtain greater precision compared to traditional methods. Additionally, unlike standard point estimates and error bars, Bayesian methods allow us to recover the full posterior distribution for each parameter. Together, this makes Bayesian analysis a powerful tool for examining globular clusters and for accurately and precisely determining their characteristics.

With well-determined ages and metallicities derived using Bayesian methods, we can re-visit the age-metallicity relation. The grouping of globular clusters into two groups in the age-metallicity relation is often cited as evidence for the Milky Way halo’s two-part formation. A group of old clusters appears to have no dependence on metallicity, while the remaining group of clusters are old and metal-poor, then become younger at higher metallicities (Sarajedini & King 1989; Chaboyer, Sarajedini & Demarque 1992; Richer et al. 1996; Rosenberg et al. 1999; Dotter et al. 2010; Dotter, Sarajedini & Anderson 2011; Leaman, VandenBerg & Mendel 2013). The older clusters likely formed during a rapid initial assembling of the inner halo, while the other clusters were possibly accreted from the cannibalization of dwarf galaxies over time to the outer halo.

Our Bayesian analysis provides estimates of the overall helium abundances of the clusters, and we use this to explore possible relationships among helium enrichment and other global cluster parameters. Helium is tightly tied to the multiple population scenario in globular clusters (e.g.: Bedin et al. 2004; Milone et al. 2009; Gratton, Carretta & Bragaglia 2012; Milone et al. 2012; Piotto et al. 2015) and we investigate the possible connections herein.

The paper is organized as follows: Section 2 discusses the data and Section 3 provides an overview of the Bayesian software suite. Section 4 presents the results of the Bayesian fitting, compares our results to recent studies, and explores the correlations we find for the clusters. In Section 5 we discuss the possible implications of our results. We present our conclusions in Section 6.

2 Data

The data used in this work are primarily from the ACS Globular Cluster Treasury program (GO Cycle 14 Proposal 10775; Sarajedini et al. 2007). This program observed 65 Galactic globular clusters over 132 orbits and produced a consistent set of deep photometry in the HST F606W and F814W filters. Additional data of 6 clusters were obtained from GO-11586 (PI: Dotter) obtained during Cycle 17 (see Dotter et al. 2010). We employ a Bayesian analysis to characterize 69 of the 71 clusters. Palomar 2 is left out of our analysis due to high differential reddening and E3 is removed from the sample due to lack of red giant branch (RGB) stars.

The ACS Globular Cluster Treasury program has provided observations of several hundred thousand stars. To help BASE-9 perform effectively, we make modest quality cuts on the photometry for the clusters before randomly selecting a subsample of stars. To rid the data of any stars with poorly determined photometry, we remove stars for which both filters fall into the outer 5% tail of the photometric error distributions. Additionally, we remove stars in the outer 2.5% tails of the distributions of X and Y pixel location errors from frame to frame, as high pixel location errors may indicate non-cluster members. The exceptions to this are the clusters NGC 5986, NGC 6397, and NGC 6779. For these clusters, we ignore the location error constraint as it removes the majority of bright stars above the main sequence turn-off that were observed in the short exposure images. However, we still remove stars in the outer 5% tail of the photometric error distributions for both filters.

With the cleaned photometry, we randomly select a subsample of $\lesssim$ 3000 stars, with half above the main sequence turnoff point (MSTOP) of the cluster and half below the MSTOP. If there are fewer than 1500 stars above the MSTOP, we match the number of stars above and below the MSTOP. This procedure is adopted to ensure a reasonable sample of stars on the sub-giant and red-giant branches, while balancing their contribution with MS stars, without the computational cost of running 10,000 or more stars per cluster.

3 Methods

3.1 Bayesian Framework

Our software suite, known as BASE-9, is a tool for fitting and characterizing observations of open clusters (von Hippel et al. 2006; De Gennaro et al. 2009; Jeffery et al. 2011; Stein et al. 2013; Hills et al. 2015), globular clusters (Stenning et al. 2016; Wagner-Kaiser et al. 2016a, b), and even individual stars (O’Malley, von Hippel & van Dyk 2013). BASE-9 provides reproducible and precise fits to stellar clusters, determines the cluster-wide parameters of age, distance, metallicity, absorption, and helium fraction, as well as individual stellar parameters of membership, binarity, and mass.

Whether using a numerical method or chi-by-eye, fitting isochrones to stellar clusters often relies on previous studies to assume fixed values of one or several parameters (e.g.: a study interested in age may use [Fe/H] or distance determinations from previous work, or may assume helium abundances). BASE-9 allows us to fix one or more parameters in this way, to use these externally-derived parameters when specifying prior distributions, or, if the data are rich enough, to fit the cluster’s parameters simultaneously. Our results presented here also provide a baseline measurement to compare to multiple population studies using similar Bayesian methods (Stenning et al. 2016; Wagner-Kaiser et al. 2016a).

BASE-9 uses adaptive Metropolis (AM) Markov chain Monte Carlo (MCMC) to estimate model parameters and map their full posterior distributions (Stenning et al. 2016). It does so by sampling from the joint posterior distribution of distance, metallicity, absorption, age, and helium of a cluster as well as individual parameters of membership, binarity, and mass. In our analyses we fix two of these parameters: binarity and metallicity (as discussed below) and marginalize over the individual stellar parameters.. The adaptive MCMC sampler is designed to explore the joint posterior distribution and use the empirical variance-covariance matrix of past MCMC draws to optimize the proposal distribution. One advantage of the AM algorithm is that the user need only provide starting values for the chain and adaptation is done automatically, such that additional tuning is generally unnecessary. After an initial burn-in period (typically 1000 iterations), we run the AM MCMC chain for at least 10,000 additional iterations. We use one chain per cluster due to computational limitations. Additional details regarding our statistical computation methods and implementation of the adaptive Metropolis algorithm are provided in Stenning et al. (2016).

We assume a fixed value of 0 for $R_{i}$, thus not including binaries in our model due to computational limitations. Globular clusters tend to have low binary fractions, thus we do not expect that treating stars as singletons should have a significant effect on the final results (see Section 4.3 for further discussion). Future development of the BASE-9 software will enable the inclusion of binaries in large stellar clusters at reasonable computational cost.

Metallicity values are assumed for our clusters from the spectroscopic measurements in the Harris (2010) catalog. Because the theoretical models have a strong correlation between metallicity and helium, we cannot gain enough leverage to constrain both simultaneously with BASE-9 using only F606W and F814W. However, fixing metallicity at a reasonable value obtained from the literature allows us to gain precise determinations of helium abundances for a large sample of clusters.

We define a statistical model based on a hierarchy of properties belonging to individual stars or the entire cluster (De Gennaro et al. 2009; van Dyk et al. 2009; Stein et al. 2013). The individual parameters are the zero-age main sequence masses of the stars ($M_{i}$; additional individual parameters are defined shortly). The cluster level parameters are age (base-10 logarithm, $\theta_{\rm age}$), metallicity ($\theta_{[Fe/H]}$), distance modulus ($\theta_{m-M_{V}}$), absorption ($\theta_{A_{V}}$), and helium fraction ($\phi_{Y}$). We represent the observed photometric data as $X_{ij}$, for the ith star in the jth filter (for $N$ stars and $n$ filters), as described by De Gennaro et al. (2009), van Dyk et al. (2009), Stein et al. (2013), and Stenning et al. (2016). The known Gaussian measurement errors are in the variance-covariance matrix $\mathbf{\Sigma}_{i}$ for each star $i$.

For a given set of proposed individual and cluster level parameters, we predict photometry for each of $N$ stars. The function $G$ is the stellar evolution model that uses $M_{i}$, $\Theta$, and $\phi_{Y}$ to predict an $n$-length vector of magnitudes for 1, $\ldots$, $N$ stars (where $\Theta$ = ($\theta_{\rm age}$, $\theta_{[Fe/H]}$, $\theta_{m-M_{V}}$, $\theta_{A_{V}}$), represented as \mbox{\boldmath{\mu}}_{i} = $G$($M_{i}$, $\Theta$, $\phi_{Y}$). Although for our purposes helium is a cluster-level parameter, it is treated separately because it is allowed to vary within the cluster in other contexts (Stenning et al. 2016; Wagner-Kaiser et al. 2016a). As in van Dyk et al. (2009), binaries can be included in the model by treating each star as a binary system and fitting the component masses. The ratio of the companion mass relative to the primary mass is referred to as $R_{i}$.

We represent the predicted magnitude for star $i$ as

[TABLE]

where \mbox{\boldmath{\mu}}_{i} is the vector of predicted magnitudes. We also account for the possibility that each star could belong to the Milky Way field population rather than to the cluster itself with a two-component finite mixture distribution (van Dyk et al. 2009). As such, we define indicator variables, $Z_{i}$ for $i$=1, …, $N$, that take the value of 0 for a field star and 1 for a cluster star. We provide a brief description of the main components of the statistical model in Table 1 for reference.

The likelihood function for a single population of stars is

[TABLE]

Equation 2 may be rewritten as

[TABLE]

where M = ($M_{1}$, …, $M_{N}$), R = ($R_{1}$, …, $R_{N}$), Z = ($\textbf{Z}_{1}$, …, $\textbf{Z}_{N}$), X = ($\textbf{X}_{1}$, …, $\textbf{X}_{N}$), and $\Sigma$ = ($\mathbf{\Sigma}_{1}$, …, $\mathbf{\Sigma}_{N}$). The conditional distribution of observed photometric magnitudes for a star belonging to the cluster is represented by P(Xi $\lvert$ $\mathbf{\Sigma}_{i}$, $M_{i}$, $R_{i}$, $\Theta$, $\phi_{Y}$, $Z_{i}$=1). Essentially, the first term inside the bracket in Equation 3 models the star as if it were a cluster star ($Z_{i}$=1), while the second term models it as belonging to the field star population.

In this way, photometry is predicted at each iteration for a proposed set of stellar and cluster parameters. Over many iterations, we use BASE-9 to explore the parameter space of age, distance, absorption, and helium. The estimates we report are posterior medians of these parameters, from which we can generate a single isochrone for each cluster. Parameters that are not of interest to this study are either marginalized over or fixed, as discussed above. The user inputs to BASE-9 include the photometry and photometric uncertainty for each star in each filter, prior distributions, and starting values for each cluster parameter. The Dartmouth Stellar Evolution Database theoretical isochrones (Dotter et al. 2008*http://stellar.dartmouth.edu/) presents a multi-dimensional space for the MCMC sampler to step through. We have generated a grid of isochrones to cover a range of ages from 1 to 15 Gyr, a range of [Fe/H] metallicity from –2.5 to 0.5 dex, and a range of helium values from $\sim$0.23 to 0.4. This DSED grid has been generated at a [$\alpha$/Fe] value of 0.0 (solar) and we assume an RV=3.1 reddening law (Cardelli, Clayton & Mathis 1989) during the analysis.

3.2 Prior Distributions

Because we perform a Bayesian analysis, we must specify prior distributions that summarize our knowledge regarding the model parameters before considering the current data. Previous studies are used to generate Gaussian prior distributions for distance and absorption for each cluster (Harris 2010). Generally speaking these priors are centered on published estimates but we use conservative prior dispersions (e.g., relatively large prior standard deviations) to reduce the influence of the priors on our final estimates. For the distance prior, we choose $\sigma$(m$-$M) = 0.05 mag. For absorption, we conservatively use $\sigma$($A_{\rm V}$) = (1/3) $A_{\rm V}$, using a Gaussian truncated at zero. Because there is no useful prior information on helium or absolute age, a uniform prior on helium of 0.15 to 0.45 is assumed and a uniform prior on age of 1 to 13.5 Gyr. We impose an upper limit of 13.5 Gyr, because although the models extend to 15 Gyr, we do not expect GGCs to be older than the epoch of galaxy formation (Pilipenko 2013; Oesch et al. 2016).

For helium, Y=0.24 is adopted as a starting point for all clusters. Ages from Dotter et al. (2010) and Roediger et al. (2014) are used as starting values for the age. For the clusters not in either of these studies, we adopt a starting point of 12 Gyr for the age. A Miller & Scalo (1979) initial mass function defines the prior distribution for initial stellar mass, over the range of 0.1 to 8 solar masses (von Hippel et al. 2006, De Gennaro et al. 2009, Stenning et al. 2016).

The cluster membership represents whether a star belongs to the cluster, as opposed to the field star population. As the cleaned photometry discussed above is expected to be largely free of field star contamination, we adopt a incidence rate of field stars of 5%, as we expect different values of the membership probability in the corresponding range of 0.9 to 1.0 produce consistent results (Stenning et al. 2016). Blue straggler stars are not included in our model. As previously mentioned, because the primary interest is in helium, which is strongly correlated with [Fe/H], the spectroscopic metallicities from Harris (2010) are used to set the metallicity to previously measured values.

The clusters are listed in Table 2 with their adopted prior distributions and starting values. We provide an example of the posterior draws from BASE-9 as well as a 2-D correlation plot in Figures 1 and 2 for NGC 6752.

4 Results

In Table 3, we provide the results from our Bayesian analyses. The values are the median of the posterior distributions for each parameter and the errors reflect 90% central Bayesian credible intervals from the BASE-9 posterior distribution. These errors represent statistical uncertainty, assuming the validity of the models (discussed further in Section 5). We use the medians and credible intervals in subsequent analyses and tables, with the quoted errors representing these statistical uncertainties.

Figure 3 presents CMDs for each cluster. All ACS photometry are shown as gray points, and the subsamples of stars used in the BASE-9 analysis are are shown in black. The cyan isochrones are generated using the median of the posterior distributions for the parameters, as presented in Table 3.

Many clusters in our sample have light element abundances at non-solar ratios. Dotter et al. (2010) classifies clusters with abundances of [$\alpha$/Fe] = 0.0, 0.2, and 0.4. However, DSED models are only available with varying helium for the two [alpha/Fe] endpoints, at [$\alpha$/Fe] = 0.0 and 0.4. For the clusters identified as [$\alpha$/Fe] = 0.4 in Dotter et al. (2010), we run them with a grid of [$\alpha$/Fe] = 0.4 instead of [$\alpha$/Fe] = 0.0. We present these results in Table 4, and show the resulting isochrone in red in Figure 3. Aside from NGC 288, the [$\alpha$/Fe] = 0.4 model grid does not appear to result in visually satisfying fits.

4.1 Comparisons to Previous Work

In this section, we present a brief comparison of our results to previously published results, specifically those of Dotter et al. (2010) and Harris (2010). Because we allow the helium fraction to vary as opposed to previous studies which fix or assume its value, other parameters change to maintain a best fit of the models to the data, so differences with published studies are not unexpected.

In Figure 4, we show a comparison of the median of our posterior distribution for distance and absorption results to the estimates for each cluster from Harris (2010) in the left and middle panels. We calculate the differences as the Harris (2010) estimates subtracted from our own. In these panels, the errorbars represent our BASE-9 90% central Bayesian credible intervals and we indicate the error estimates from Harris (2010) as dotted lines. Specifically, an error of $\pm$0.1 for distance and $\pm$0.1*Av for absorption. Our distance estimates largely fall within the range of error from of distances from the Harris (2010) catalog. As one might expect, it appears that for larger values of absorption, the scatter increases. In many cases, we estimate moderately larger absorption than the values from Harris (2010).

In the rightmost panel of Figure 4 we compare the BASE-9 and Dotter et al. (2010) estimates of age†††We do not plot $\Delta$(BASE-9 - Dotter) on the vertical axis of the rightmost panel of Figure 4 because errors on $\Delta$(BASE-9 - Dotter) are not readily computable. This is because the Dotter et al and the BASE-9 estimates are compiled using the same ACS photometry from Sarajedini et al. (2007), leading to correlation between the estimates so that the error on $\Delta$(BASE-9 - Dotter) is not simply the sum (in quadrature) of the errors of the two estimates.. While the two studies largely agree within the (correlated) errors, it is not unexpected that we should find slightly different results from Dotter et al. (2010) as we explore variations in helium and not in metallicity. We provide relative ages where available for the clusters with respect to Dotter et al. (2010) in the final column of Table 3.

4.2 Age$-$Metallicity Relation

We revisit the age-metallicity relationship with our results, as seen in Figure 5. The age-metallicity relation (AMR) has long been held as support for a two-phase formation scenario of the Milky Way (Chaboyer, Demarque & Sarajedini 1996, Marín-Franch et al. 2009, VandenBerg et al. 2013). This view of the formation of the Galaxy suggests that the bulk of the inner halo formed early-on along with the metal-poor clusters of the outer halo. Cluster formation then continued mainly in the outer halo, forming more metal-rich and younger clusters.

Consistent with previous studies, we find a group of old clusters with a relatively narrow age range almost exclusively found in the inner 8 kpc of the Galaxy, with a error-weighted mean age of 13.433$\pm$1.005 Gyr (log(age) of 10.124$\pm$0.033). The group of clusters beyond 8 kpc follow a different trend, such that as these clusters become more metal-rich, they also become younger with more significant scatter. For these clusters, we find an error-weighted mean age of 12.320$\pm$1.755 Gyr (log(age) of 10.067$\pm$0.069). We see a clear distinction between the two sequences and consider this strong evidence for the two-part formation history of the Milky Way. The old, inner halo clusters likely formed quickly early on in the history of the Galaxy, while the younger trend of clusters slowly formed or joined the outer halo over time via chaotic accretion. This is also thought to be tied to the horizontal branch morphology–age relationship and the second parameter problem (Searle & Zinn 1978; van den Bergh 1965; Sarajedini & King 1989; Dotter et al. 2010). However, our results do not show clear trends with helium in the context of the age-metallicity relation.

Many cluster ages are pinned against the age of the Universe, an external prior built into our analysis. This indicates that the isochrones prefer an older fit. Nonetheless, despite the limitations in the absolute ages of isochrone fitting, we are focusing on relative ages, specifically on distinguishing the median age and age spread of the inner clusters versus the outer clusters.

There are several inner clusters within 8 kpc of the Galactic center that we find along the “wrong” sequence. One of these is $\omega$ Cen (NGC 5139), which is both the most massive globular cluster and which likely contains the largest number of stellar populations (Bedin et al. 2004; Bekki & Norris 2006; Bellini et al. 2010; Johnson et al. 2008). NGC 6584 has a large rotational velocity, making it probable that it does not truly belong to the “inner” galactic clusters. Instead, it is likely an outer halo cluster that is observed to be passing through the inner halo. It’s unclear why the other clusters, NGC 5904, NGC 6254, and NGC 6752 follow the sequence of outer halo clusters rather than that of the inner clusters. We indicate these clusters with a red outline in Figure 5 to show their locations.

4.3 Helium

In this section, we present a variety of correlations we find among helium and other cluster parameters, and discuss the possible implications of these findings.

The distribution of fitted helium fractions in the clusters is approximately symmetric, as seen in Figure 6, with a total range of about 0.17, from Y = 0.235 (NGC 6652) to Y = 0.409 ($\omega$ Cen). The peak in the distribution is at Y = 0.321$\pm$0.035. Because we are probing the overall (“average”) helium abundance of the cluster, we do not expect helium abundances lower than the primordial value. Our lowest estimated helium abundance is 0.235, similar to the lower limit on the primordial value (we list a few recent values in Table 5 for comparison). As almost, if not all, these clusters harbor multiple populations, the second (or third, fourth, etc.) generation stars in each cluster may be more enriched in helium. The helium abundances measured should be population-weighted values (e.g.: a higher measurement of helium could be caused by a higher proportion of later generations of stars, a higher level of helium enrichment, or both). Our helium fraction estimates are also sensitive to our assumption of metallicity abundances, with more metal-rich clusters being more sensitive to the assumed metallicity.

We find a statistically significant correlation between the helium abundance of a cluster and its metallicity, as seen in Figure 7, with errors given by the 90% central Bayesian credible intervals from Table 3. We fit a broken line to the data (solid line), with Y = 0.348 for [Fe/H]$\textless$–1.5 and Y = –0.122($\pm$0.022) [Fe/H] + 0.165($\pm 0.024$) for [Fe/H]$\textgreater$-1.5. The relationship has an error-weighted Pearson correlation coefficient of –0.545$\pm$0.001.

Overall, we find that the metal-poor clusters have similar helium abundances. As clusters become more metal rich, they appear to become less helium-rich. However, the clusters NGC 6304, NGC 6366, NGC 6388, NGC 6441, NGC 6624, and Terzan 7 do not follow this trend, instead being well-explained by the model assumption from Dotter et al. (2008) of $\Delta$Y/$\Delta$Z = 1.54.

The spread in helium for metal-rich clusters could be due to a high percentage of chemically unenriched stars or due to less helium enrichment compared to other clusters. In our study, we are unable to distinguish between the two effects.

More massive clusters are predicted to attain higher enrichments of helium (Milone et al. 2014, Milone 2015). As shown in Figure 8, we find that the absolute magnitude of the clusters from Harris (2010) (a proxy for the mass of the cluster) is likely correlated to the helium of the cluster, such that brighter clusters tend to have marginally higher values of helium. The weighted Pearson correlation coefficient between helium and the absolute magnitude is moderate at –0.261$\pm$0.032, but the trend suggests support for previous claims.

We also see a relationship between the overall helium content and the binary fraction in the core of the cluster, as seen in Figure 9, using the observed binary fractions from Milone et al. (2012). As the helium abundances of the clusters decrease, the incidence of binaries increases. This trend continues until the binary fraction reaches about 40%, after which the correlation appears to be weaker. However, there are few clusters beyond a 40% binary fraction, and the errors tend to be larger.

The binary fractions were determined by Milone et al. (2012), who noted that the more massive clusters (using absolute magnitude as a proxy for mass) had lower binary fractions. From Figure 8, we find that the less massive clusters have lower helium content, and Figure 9 shows that these same clusters have higher binary fractions, consistent with Milone et al. (2012). While our observed trend between helium and binary fraction may be driven primarily by the anti-correlated mass-binary fraction, it may be important to consider what the binary fraction means with respect to helium, as it may inform the multiple population scenario. Because the binary fraction determinations and helium estimates derive from the same ACS photometry (Sarajedini et al. 2007), they are expected to be correlated. Although our tests suggest that binaries are largely classified as field stars in our BASE-9 analysis and do not contribute significantly to the posterior distribution, it is possible that remaining binaries classified as cluster members could weight the analysis towards slightly lower helium abundances. Thus, the relationship between binary fraction and helium abundance requires further investigation.

Both Figures 8 and 9 indirectly suggest a relationship between the helium of a cluster and its mass. Thus, it appears that there is a preference for clusters of higher mass to also have a greater overall enhancement of helium. This relationship has been previously seen for a smaller sample of clusters (Milone et al. 2014, Milone 2015), but we find that it holds true for the majority of Galactic globular clusters. A recent study by Lucatello et al. (2015) showed that second generation stars, which are presumably more helium enhanced, have a lower incidence of binaries. We also find that the more enhanced clusters, which likely harbor a greater proportion of second generation stars, have fewer binary systems.

4.4 Light Elements

Additionally, we also find a possible correlation between helium and light element abundances. Understanding the connection between helium and the levels of carbon, oxygen, and nitrogen is a crucial step to disentangling the progenitors of the multiple population phenomena. Using the results from Roediger et al. (2014) obtained from compiling spectral data from the literature, we compare our helium estimates to carbon, oxygen, and nitrogen abundances in Figure 10 for a small sample of clusters with spectroscopic measurements.

Despite the relatively small sample of clusters with CNO spectral observations and the large measurement errors, clusters with higher helium tend to have lower carbon abundances and greater nitrogen abundances. There is a strong correlation with a Pearson coefficient of –0.742 $\pm$ 0.004 for carbon and a moderate correlation (0.436 $\pm$ 0.240) for nitrogen. For oxygen, there is no apparent relationship with helium, with a correlation coefficient of essentially zero (0.060 $\pm$ 0.812).

For many proposed multiple population scenarios, the first generation of evolved stars provides processed material to the formation of the later generations of stars (Renzini 2008, Gratton, Carretta & Bragaglia 2012, Bastian, Cabrera-Ziri & Salaris 2015). The CNO cycle taking place in these progenitors would be expected to deplete carbon and oxygen abundances and boost nitrogen abundances, while the same processes produce helium. Our observations of these effects for carbon and nitrogen supports this scenario of multiple population formation, except we find little to no change in the oxygen abundances. While in Figure 10 we are comparing the average abundances for the cluster against the average helium value of the cluster, future work should aspire to compare the light element abundances and the helium enrichment of each population. Although there are only a handful of clusters with measured abundances for carbon, oxygen, and nitrogen, we do not see a relationship between helium and total CNO.

4.5 Red Giant Branch Bump

The red giant branch bump (RGBB) occurs due to the hydrogen-burning shell in the star moving outward during the first ascent of the RGB. Stars undergoing this event become temporarily brighter, until the hydrogen-burning shell reaches the chemical discontinuity near the convective envelope in the star. After this, the stars grow fainter again, causing a “pile-up” in the CMD along the RGB. Awareness of the RGBB, clearly visible in many CMDs along the RGB, dates back to its prediction by Thomas (1967) and its confirmation by King, Da Costa & Demarque (1985). The RGBB of a cluster has long been predicted to be sensitive to helium, age, and metallicity (Cassisi & Salaris 1997, Bono et al. 2001, Bjork & Chaboyer 2006, Salaris et al. 2006, Di Cecco et al. 2010, Cassisi et al. 2011, Nataf et al. 2013). Such studies purport that an increase in initial helium abundance will increase the brightness of the RGBB.

Using measurements from Nataf et al. (2013), we compare the absolute V-magnitude of the RGBB to the metallicities and helium abundances in Figure 11. There is a clear trend in both, supporting predictions and previous work. The correlation between [Fe/H] and the RGBB magnitude is very strong, and appears to be non-linear. We fit an X and Y error-weighted quadratic to the data:

[TABLE]

We also fit a linear model, again accounting for errors in both directions, to helium and absolute RGBB magnitude:

[TABLE]

These fits are shown as solid lines in the panels of Figure 11 ([Fe/H] on the left and helium on the right).

The correlation between helium and RGBB absolute magnitude is evident, with a weighted Pearson coefficient of –0.364$\pm$0.009. However, metallicity is clearly more important, and may drive a significant portion of the relationship between Y and RGBB magnitude. In the left panel of Figure 12, we plot the residuals from the fit between metallicity and RGBB magnitude against the helium of each cluster. The residuals have a median of –0.019 and a standard deviation of 0.068. An X and Y error-weighted linear model to the residuals results in a slope of –0.429 $\pm$ 0.271 and an intercept of 0.137 $\pm$ 0.088, suggesting that some effect due to helium remains. With the effect on the RGBB magnitude due to metallicity removed, our fit predicts that a change in 0.01 in Y results in a brighter RGBB by 0.008 magnitudes. This is similar to the result of Cassisi & Salaris (1997), who suggests that for the same increase in helium, the magnitude of $\Delta$(Vbump-VTOP) becomes brighter by 0.011 magnitudes. Salaris et al. (2006) also suggests a minor shift in the RGBB magnitude for an increase in helium.

Using our linear fit (shown in the left panel in Figure 12), we remove the remaining dependence of the RGBB magnitude on Y to examine any remaining residual after removing the dependencies on metallicity and helium. The remaining scatter is plot against the age of the clusters in the right panel of Figure 12. If we attempt to fit an error-weighted linear relation, the intercept is consistent consistent with zero. The remaining scatter could be from uncertainties in the distances or absorptions. The median residual is –0.004 with a standard deviation of 0.015; the standard deviation of the residuals is reduced by a factor of 4.5 after the removal of the remaining trend with helium.

Interestingly, we do not see any clear trends between the brightness of the RGBB and the light element abundances.

4.6 RR Lyrae Periods

We also examine our results in the context of Lee, Demarque & Zinn (1994), who explored the theoretical dependence of the horizontal branch morphology on parameters such as age, metallicity, and helium. Lee, Demarque & Zinn (1994) predicted that increased helium would lead to a greater $\Delta$log(P), which they defined as the difference between the fundamentalized average periods of clusters’ RR Lyraes compared to M3. The fundamentalized period for each cluster was defined by (log(Pf) = log(Pc)+0.13), with Pc being the average period of RR Lyrae c-type stars in each cluster (Castellani & Quarta 1987). Consequently, the study concluded that the fundamentalized periods of clusters were not affected by helium because the observational data did not follow the theoretical predictions, which suggested that an increase in helium would be reflected in an increased in $\Delta$log(P).

The fundamentalized periods of RR Lyrae type c observations from Castellani & Quarta (1987) and the Harris (2010) metallicities are used to recreate Figure 10 from Lee, Demarque & Zinn (1994), as seen in the left panel of Figure 13. In this panel, the fundamentalized RR Lyrae periods, with respect to the fundamentalized period of M3, are plot against the metallicities. A quadratic is fit to the relation

[TABLE]

We exclude the two outliers, NGC 2419 (lower) and NGC 6333 (upper) in the fit. Leaving out these same two clusters, we find that residuals of Equation 6 are related to helium by the following relation:

[TABLE]

Qualitatively, our results match the general trend predicted by the theoretical calculations from Lee, Demarque & Zinn (1994) that an enriched Y value leads to a greater $\Delta$log(P). However, the trend is not statistically significant, suggesting that this approach cannot provide any hard conclusions on the possible effects of helium on the horizontal branch morphology.

4.7 Anomalous Clusters

While most globular clusters harbor distinct populations that vary in helium and light element abundances, a smaller subset of clusters has been found to also show variations in [Fe/H] abundances and s-process element abundances (Bedin et al. 2004; Gratton, Carretta & Bragaglia 2012; Marino et al. 2015). Although only a handful of clusters have been categorized as anomalous clusters, there are likely more that have not yet been the target of sufficient high-resolution spectroscopy to be identified as such.

Several anomalous clusters in the Milky Way overlap our cluster sample: $\omega$ Centauri (NGC 5139), NGC 1851, M22 (NGC 6656), M2 (NGC 7089), M54 (NGC 6715), and NGC 5286 (Milone et al. 2009, 2013; Marino et al. 2015). Except for NGC 1851, the most metal-rich and the youngest of the anomalous clusters in our sample, the clusters have above average helium abundances. Otherwise, we do not notice anything particularly unusual about the characteristics of the anomalous clusters in our sample, probably due to only a rough sensitivity to the complex multiple population characteristics from visual photometry.

5 Discussion

We stress that our results are dependent on two primary assumptions. First, we assume the clusters analyzed have a single iron abundance and a single age. For some clusters with multiple populations this may not be the case. In a few clusters, ages of different populations of stars could vary and the internal abundance spread may be significant. However, for the handful of clusters where this may be the case, we expect that the internal variation in age is within the fitting errors.

Secondly, our methodology is dependent on the validity of the theoretical models, which are assumed to be accurate (von Hippel et al. 2006; De Gennaro et al. 2009; Stein et al. 2013; Hills et al. 2015). Our results, while precise, can therefore only be as accurate as the models. We note that the quoted uncertainties on the parameters derived herein represent the statistical uncertainty. This does not incorporate uncertainties in an astronomical sense that may arise from model misspecification, unincorporated effects, incorrect assumptions about physical processes, or observational systematics. Once models improve, the approach presented here will allow us to focus on absolute quantities, such as absolute helium abundances and absolute ages. In the meantime, the high precision results we obtain from BASE-9 fits with HST data allow us to perform numerous relative analyses.

The fitted helium values can be interpreted as “average” helium values, which depend not only on the extent of the helium enhancement of stars in the cluster, but also the fraction of stars that are enhanced in helium. In this study, these two are degenerate. Further work is needed to disentangle these effects and study their individual influences on the correlations we see.

The relationships we find among the helium fraction and the light element abundances of carbon and nitrogen, and the binary fraction, are important avenues for further investigation with respect to the multiple population scenario. Ultraviolet observations of these clusters are now a common avenue for analyzing multiple populations. The ultraviolet wavelengths encode the relative abundances of the light elements carbon, nitrogen, and oxygen, rather than directly probing helium of a cluster (Piotto et al. 2015). Thus, it is very important to understand any relationships between abundances of various light elements and the helium content. This is especially true as helium and $\alpha$-enhancement tend to have off-setting effects on the shape and location of theoretical isochrones. These effects are more dramatic in the ultraviolet than in the visual wavelengths, but in either case, are generally more exaggerated for the more metal-rich clusters.

The correlation between helium abundance and binary fraction is also of interest to multiple population research. In multiple population formation models for clusters that invoke multiple star formation episodes, it is often suggested that the first generation of stars was significantly ($/geq$ 10$\times$) more massive than the present-day cluster masses (Bastian, Cabrera-Ziri & Salaris 2015; Cottrell & Da Costa 1981; Gratton, Sneden & Carretta 2004; D’Ercole et al. 2008; Decressin et al. 2007). The first generation of stars is then largely lost to the field star population. Previous work has suggested that first generation stars in clusters have higher binary fractions (up to 15%) compared to later generations, which have only a few percent (Milone et al. 2008, Lucatello et al. 2015). Our results show the more massive clusters with a lower binary fraction tend to have more helium, as we may expect from stars belonging to the chemically enriched population. We also observe that the less enriched clusters appear to be less massive and have a higher incidence of binaries, as is predicted for first generation dominated clusters. Further investigation of the binary fraction with respect to helium may be able to shed more light on the different binary fractions in globular clusters (and their sub-populations) and field stars (Milone et al. 2008, D’Orazi et al. 2010, Gratton, Carretta & Bragaglia 2012), thereby providing information on the source of the field star population.

6 Conclusions

We analyzed 69 Galactic globular clusters to simultaneously determine their their precise relative ages, distances, extinctions, and helium fractions with statistical rigor. Based on the results, we reach the following conclusions:

1. We find a symmetric distribution of helium fractions for the clusters, ranging from approximately 0.235 to 0.409 with a peak at Y $\sim$ 0.321.

2. In the age-metallicity relation, there is a clear demarcation between two unique age-metallicity sequences, with older clusters being primarily located in the inner halo, and younger clusters in the outer halo. This suggests a multi-phase formation and evolution of the Milky Way, consistent with previous studies.

3. Correlations between the helium fractions of the clusters and both their absolute magnitudes and their binary fractions lend additional evidence to the suggestions that more massive clusters harbor higher helium contents.

4. Carbon and nitrogen abundances are found to be correlated with the overall helium fractions of several clusters. Although spectral measurements are still relatively sparse, it appears that higher helium is associated with a lower [C/Fe] ratio and a greater [N/Fe] ratio. We find no statistically significant relationship between helium and oxygen abundance.

5. The absolute magnitude of the RGBB and the helium fraction are related, as predicted from previous studies; a change of 0.01 in helium results in a brighter RGBB by 0.008 magnitudes.

Acknowledgments

We are grateful to the referee, Raffaele Gratton, whose comments greatly improved the paper. Support for this work (proposal number GO-10775) was provided by NASA through a grant from the Space Telescope Science Institute which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This material is based upon work supported by the National Aeronautics and Space Administration under Grant NNX11AF34G issued through the Office of SpaceScience and through the University of Central Florida’s NASA Florida Space Grant Consortium. DS work was supported by NSF grants DMS 1208791. DvD was partially supported by a Wolfson Research Merit Award (WM110023) provided by the British Royal Society and by Marie-Curie Career Integration (FP7-PEOPLE-2012-CIG-321865) and Marie-Skodowska-Curie RISE (H2020-MSCA-RISE-2015-691164) Grants both provided by the European Commission.