The Effect of Age and Stellar Model Choice on Globular Cluster Color-to-Metallicity Conversions

The photometric colors of globular clusters (GCs) act as effective proxies for metallicity, since all normally used optical/IR color indices exhibit a nonlinear but monotonic relation between their integrated color and their metallicity. One color index, (g - z) or (F475W - F850LP), has been spectroscopically calibrated in several studies, providing leverage to define color-to-metallicity conversions for other indices. In this paper, building on the work of arXiv:2307.01863, we study the GC color-metallicity relation in more detail by testing the dependence of the relations on different suites of stellar models and different assumed GC ages. Though noticeable differences between models exist, we find that the net effect on the derived GCS metallicity distributions is small.


INTRODUCTION
Globular clusters (GCs)-old, massive, dense star clusters found in all but the least massive of dwarf galaxies (Harris 2010;Forbes & Remus 2018;Beasley 2020;Eadie et al. 2022)-are powerful tracers of early galaxy growth and chemical enrichment.Part of their appeal for observers stems from the scaling relations that they exhibit with galaxy properties, such as the M GCS -M h relation (e.g.Blakeslee et al. 1997;Forbes et al. 2016;Harris et al. 2017;Dornan & Harris 2023), and from the relations between fundamental GC properties and easyto-observe quantities.The relation between GC metallicity and color, for example, is monotonic (see Brodie & Strader 2006;Peng et al. 2006;Usher et al. 2015;Harris et al. 2017;Fahrion et al. 2020;Harris 2023, among others), making color an attractive proxy observable for metallicity.
Despite its usefulness, the specifics of the GC colormetallicity relation (CMR) have proven challenging to constrain.When spectroscopy is available, the relation is often modelled nonlinearly, as in Peng et al. (2006); Sinnott et al. (2010); Usher et al. (2012); Fahrion et al. (2020), andHarris (2023) (with the exception of Villaume et al. 2019).Furthermore, the relation has been spectroscopically studied primarily for the (g -z) color index (equivalent to the HST index (F475W -F850LP) in the Vegamag system); thus other color indices are useful metallicity indicators but must first be transformed to (g -z) with a color-color relation, as in Harris (2023) and Hartman et al. (2023).
Both Harris (2023) and Hartman et al. (2023) used the PARSEC stellar models (Marigo et al. 2017) to derive GC color-color relations.They adopted a fixed age of 12 Gyr for their simulated clusters and fitted simple nonlinear models to the resulting CMRs.Those studies focused on the GC metallicity distributions as they related to host galaxy properties; while the steps taken to construct color-color and color-metallicity relations were justifiable, the relations were a means to an end and were not themselves investigated in detail.Because GC integrated colors are dominated by relatively small numbers of red-giant and subgiant stars per cluster, stellar model specifics in these evolutionary stages could have a significant impact on the shape of GC color-color relations, with subsequent effects carried through the CMRs.
In this paper, we take further steps towards building the CMRs for old GCs by investigating the effect of stellar model choices on the shape and zeropoints of the relations.For the purposes of this study, we restrict the analysis to the particular color index (F475X-F110W) used in Hartman et al. (2023) for the GC systems in 15 giant early-type galaxies (ETGs), with our main goal to test the model-based transformations themselves.We will address extensions to other indices in later work.Section 2 details the properties of our simulated clusters, Section 3 outlines our procedures for deriving CMRs, and Section 4 looks at how the differences in CMRs affect analysis of observations.We discuss our findings and summarize our work in Section 5.
We derived GC CMRs using five suites of frequently used contemporary models: two versions of the PAR-SEC (PAdova and TRieste Stellar Evolution Code) models, two versions of the BaSTI (a Bag of Stellar Tracks and Isochrones) models, and the MIST (MESA Isochrones and Stellar Tracks) models.
First, we tested two recent versions of the PARSEC cluster simulator: CMD 3.6, the model used in Hartman et al. (2023), and CMD 3.7, the most up-to-date set available at the time of writing.We used the default settings from the Osservatorio Astronomico di Padova's online tool (CMD 3.6, CMD 3.7) for circumstellar dust, extinction, and long period variability, and specified the cluster age, metallicity, and mass.See Bressan et al. (2012); Chen et al. (2014Chen et al. ( , 2015)); Tang et al. (2014); Marigo et al. (2017); Pastorelli et al. (2019), andPastorelli et al. (2020) for details.
Next, we tested two abundance mixtures from BaSTI: solar-scaled, with [α/Fe] = 0.0, and alpha-enhanced, with [α/Fe] = +0.4,both available as drop-down options on the synthetic CMD page of BaSTI's online tool.We used the diffusionless grid for our solar-scaled clusters and the He = 0.247 grid for our alpha-enhanced clusters, along with default values for star formation rate (SFR), minimum mass, binary fraction, minimum binary mass ratio, treatment of variable stars, and photometric error.We also specified an SFR scale of 500,000 (this input is related to the mass of the resulting simulated cluster).BaSTI is the only model suite in this study with an alpha enhancement option as of the time of writing.See Hidalgo et al. (2018); Pietrinferni et al. (2021), andSalaris et al. (2022) for details.
Finally, we tested one set of MIST models with MIST's isochrone interpolation online tool.We used the default settings for MIST version, stellar rotation, extinction, and abundances.See Dotter (2016) and Choi et al. (2016), along with Paxton et al. (2011Paxton et al. ( , 2013Paxton et al. ( , 2015)), and Paxton et al. (2018) for details.
For each of the 5 suites of stellar models, we created a 5 × 23 grid of simulated clusters as they would appear in the HST WFC3 photometric system, with a metallicity range of -2.2 to 0.0 dex in 0.1-dex steps and an age range of 9 to 13 Gyr in 1-Gyr steps.Each simulated cluster had approximately the same mass (M ∼ 10 5 M ⊙ ), a Kroupa initial mass function, and no internal spread in metallicity [M/H].We used these simulated cluster sets to calculate predicted CMRs.

DERIVING CMRS
In order to recreate the color to metallicity conversion construction process from Hartman et al. (2023), we needed integrated magnitudes from simulated clus-ters for each of the filters in their color indices.The BaSTI online tool provides these automatically, producing two outputs for each model cluster: a large file containing magnitudes and stellar properties for each simulated star, and a small file with cluster properties and integrated magnitudes.Conversely, the MIST online tool produces isochrones rather than simulated clusters, and the PARSEC CMD 3.6 and 3.7 online tools' integrated magnitude option removes the cluster mass parameter, so we calculated integrated magnitudes for each PAR-SEC and MIST cluster manually.

Integrated cluster magnitudes
The output files for both PARSEC CMD 3.6 and 3.7 included magnitudes in all HST WFC3 filters for each individual star in the simulated cluster.We converted these magnitudes to luminosities, summed them to get an integrated cluster luminosity, and then converted that back into a magnitude for each relevant filter.
For MIST, an output file contains for each point on the isochrone a set of theoretical stellar parameters (e.g.mass, effective temperature, luminosity, etc.).For the MIST models, we calculated the number of stars at each point necessary for a 10 5 M ⊙ GC with a Kroupa IMF (Kroupa et al. 1993), and then followed the same luminosity summing procedures that we used with the PAR-SEC clusters.
Figure 1 shows (F475X -F110W) and (F475W -F850LP), the two color indices from Hartman et al. (2023), versus metallicity for individual simulated GCs in each of the five models.Because the MIST online tool does not produce simulated clusters like its PAR-SEC and BaSTI counterparts do, those clusters are not affected by stochasticity and do not exhibit any scatter.

Fitting
We fitted inverse exponential equations to each set of simulated GCs in both (F475X -F110W) and (F475W -F850LP), for a total of 10 fitted equations per stellar model with 5 per color index.We also fitted all simulated GCs regardless of age from each stellar model; the fitted parameters are defined in Equation 1 and listed in Table 1.An exponential model is preferable because unlike a quadratic model it remains monotonic even when used to extrapolate (although significant extrapolations should still be viewed with caution), and unlike a pair of joined linear models it assumes no sudden change in slope.Figure 2 shows model GCs of all ages with fitted exponential models overplotted, along with residuals for the CMRs for both color indices; Figure 3 places all fitted models in metallicity-color space; and Figure 4 shows differences when subtracting the color-metallicity and color-color relations of the other four models from the PARSEC CMD 3.6 color-metallicity and color-color relations that we use as a template.(That is, for a given value in (F475X -F110W), how much more metalpoor or metal rich, or bluer or redder in (F475W -F850LP), is the alternate model compared to PARSEC CMD 3.6?)The synthetic color-color relations show agreement within 0.15 magnitudes among the different sets of models, and the CMRs agree within 0.3 dex.
The full set of fitted parameters, broken out by GC age, for an equation of the form can be found in Appendix A.

CMR uncertainty
The uncertainties of metallicity values derived from the CMRs presented in this work can be estimated by multiplying a user-supplied color uncertainty ∆color by the slope of the CMR d[M/H] dcolor : where with B and C found in Table 1 as with Equation 1. Figure 5 shows the ratio of metallicity uncertainty to color uncertainty for each model in (F475X -F110W) and (F475W -F850LP); because of the change in slope of the CMRs, the red ends of the CMRs are more sensitive to metallicity than the blue ends, so the metallicity uncertainties at the red ends are smaller.

The GC age-metallicity relation
GC metallicity is expected on both theoretical and observational grounds to be correlated at least slightly with age, in the sense that more metal-rich clusters, on average, belong to a later stage of formation and are thus younger (Leaman et al. 2013;Choksi et al. 2018;Li & Gnedin 2020;Horta et al. 2021).To model the correlation between GC metallicity and age, we added an age-based correction to our color-metallicity procedure.
We calibrated the GC age-metallicity relation by combining data for Milky Way GCs from Dotter et al. Figure 6 shows the original data with Equation 4 overplotted in black.Though there is no guarantee that GCs in other large galaxies, such as the giant early-type galaxies studied in Harris (2023) and Hartman et al. (2023), will follow this same mean relation, its general property that metal-richer GCs should be younger on average is expected to hold (cf. the references cited above).The main purpose of the present study is, instead, to explore the sensitivity of the resulting GC metallicity distribution function to an age-metallicity relation.Equation 4 estimates a GC's age based on its metallicity.To determine the corresponding adjustment from the nominal CMRs listed in Table 1 for each stellar model, we divided the average color index range by our model age range (i.e. 4 Gyr) to determine an age-based color offset.After finding an individual GC's color offset based on its age and adding it to the measured GC color, we applied the nominal CMR again and obtained an age-corrected metallicity value.
Across all five stellar models, the age correction made very little difference to the CMR (see Figure 7): all corrections were less than 0.05 magnitudes.The agecorrected relations were slightly flatter than the original relations, with the extreme red and blue ends of the relation moving to less extreme metallicity values, but never leaving the color bounds set by the oldest and youngest model clusters.

APPLICATION TO OBSERVATIONS AND DISCUSSION
In order to compare the performance of each of our CMRs, we tested them on the data from Hartman et al. (2023).Their sample comprised 15 massive elliptical galaxies with densely populated GCSs, at distances ranging from approximately 60 to 110 Mpc.Because GCs appear as unresolved point sources to HST 's WFC3 camera at those distances, they extracted photometry using DOLPHOT (Dolphin 2000) and calculated integrated (F475X -F110W) color indices for each GC in We followed the same analysis procedure from that work, but substituted in our alternate CMRs before using GMM (Muratov & Gnedin 2010) to fit double Gaussian curves: where A mp and A mr are amplitudes, µ mp and µ mr are peak positions, and σ mp and σ mr are peak widths for the metal-poor and metal-rich Gaussian curves, respectively.Figure 8 shows metallicity histograms for each of the 15 galaxies in their sample, and Figure 9 shows double Gaussian fits to each metallicity distribution function (MDF) (standard for GCSs; see e.g.Kim et al. 2013;Cantiello et al. 2014;Brodie et al. 2014;Fensch et al. 2014;Escudero et al. 2015;Harris et al. 2017).Figures 8 and 9 contain the key results from this study: the shape of the MDF and the specifics of its Gaussian components are robust against the use of different stellar models and assumptions about GC age.Both the MDFs and the fitted double Gaussian curves are visually very similar to each other for each galaxy, producing the same distribution of metallicity values within the uncertainties.Figure 10 displays the parameter populations in box plot form (and the double Gaussian parameters themselves can be found in Appendix B).
Visual inspection (see Figure 11) indicates that, like stellar model choice, our built-in age correction makes no significant difference to subsequently fitted double Gaussian parameters; this was expected based on the sub-0.05-magnitudecolor changes that the correction produced.The underlying reason for this lack of sensitivity appears to be simply that in this high 9-to 13-Gyr age range, the positions of the stellar red giant tracks and thus the integrated GC colors are quite insensitive to age, and it is only the cluster metallicity that is most important.Said differently, at typical GC ages, most of the stars still present are low-luminosity and cool; beyond an age of 8 Gyr, the position of the red giant branch in CMDs changes so slowly that age effects are secondary to metallicity effects.
Even in metallicity ranges with the most disagreement between models, the differences in derived metallicities are no greater than ∼ 0.1 dex in magnitude.Based on expected color uncertainties and the corresponding metallicity uncertainties from Equation 2and Figure 5, for typical globular cluster metallicities, all five models agree with each other within uncertainty.At the GCS level, small differences in color-color transformations are effectively washed out, as demonstrated by the box plots in Figures 10 and 11.

CONCLUSIONS
In this study, we have compared five suites of stellar models, using each of them in the GC color-colormetallicity conversion process from Hartman et al. (2023) and quantifying their effects on the analysis of real GCS data.We have also probed the effect on the CMR of adding an age-metallicity correlation drawn from the Milky Way GCs.In all cases, regardless of stellar model choice or presence of the age correction, there was no significant change in the values used to characterize the MDFs derived from real data.
The data used in testing these models comes from 15 elliptical galaxies with very similar stellar masses and GCSs with of order 10 3 GCs; it would be prudent to investigate whether the results hold for smaller GCSs.In this work, we have studied the CMR primarily for the single WFC3 color index (F475X -F110W); in followup work, we will extend the transformations to other HSTbased colors including the ones used in Harris (2023).In an upcoming study with new HST multicolor imaging, the current heavy reliance on the stellar models will be reduced through the construction of strictly empirical transformations between color indices.Canada (NSERC) through a Discovery Grant to WEH.KH thanks Laura Greggio for her advice on working with MIST, and Veronika Dornan, Claude Cournoyer-Cloutier, and Jeremy Karam for helpful discussions.
B. FITTED DOUBLE GAUSSIAN PARAMETERS Double-Gaussian fits obtained through GMM, for the giant galaxies studied in Hartman et al. (2023), are listed in Table 3.

Figure 1 .
Figure 1.Model GCs from the PARSEC, BaSTI, and MIST SSPs.With all five models, we simulated GCs with metallicities ranging from -2.2 to 0.0 dex and ages ranging from 9 to 13 Gyr, and calculated integrated magnitudes in (F475X -F110W) and (F475W -F850LP).Top left: PARSEC CMD 3.6, the model used in Hartman et al. (2023).Top right: PARSEC CMD 3.7, the latest version of PARSEC available online as of the time of writing.Middle left: BaSTI's α-enhanced model.Middle right: BaSTI's solar-scaled model.Bottom left: MIST isochrone-derived model GCs.

Figure 3 .
Figure 3. Fitted exponential metallicity-color relations for the five models and the five GC age groups.

Figure 4 .
Figure 4. Differences between fitted exponential colormetallicity and color-color relations, comparing PARSEC CMD 3.6 to the other four models (i.e.subtracting each model from PARSEC CMD 3.6).As in previous figures, PARSEC CMD 3.7 is shown in blue, BaSTI alpha-enhanced in orange, BaSTI solar-scaled in red, and MIST in pink; similarly, relations based on the (F475X -F110W) color index are shown as solid lines and relations based on the (F475W -F850LP) index as dashed lines.

Figure 5 .
Figure 5. Ratio of metallicity uncertainty and color uncertainty by color index for each of the five stellar models and both color indices; see Equation2.Color is a less sensitive tracer of metallicity at the blue end of color-metallicity space than at the red end.

Figure 6 .
Figure 6.Age vs. metallicity for Milky Way GCs.Data are from VandenBerg et al. (2013), Forbes & Bridges (2010), and Dotter et al. (2010); the black line shows Equation 4, a fitted linear age-metallicity relation based on all three datasets.

Figure 7 .
Figure 7. GC age corrections for the fitted CMRs.The corrected relations are slightly flatter than the CMRs for specific GC ages in the PARSEC and MIST models, and slightly more curved in the BaSTI models.In all cases, the difference between fitted relations for individual ages and the age-corrected relation is very small.

Figure 8 .
Figure 8. MDFs for the sample from Hartman et al. (2023), transformed from color distribution functions using each of the five stellar models.

Figure 9 .
Figure 9. Double Gaussian fits for the MDFs from Figure 8.

Figure 10 .
Figure 10.A comparison of double Gaussian parameters from the fitted models in Figure 9.With the exception of the outlier σmp from the PARSEC CMD 3.7 model (see the lack of a local minimum in that model for NGC 4839 in Figure 9) the models produce very similar sets of double Gaussian parameters.

Figure 11 .
Figure 11.A comparison of double Gaussian parameters from the 5 stellar models with (lighter boxes) and without (darker boxes) age corrections.

Table 1 .
CMR parameters (as defined in Equation1) for each stellar model and both color indices, derived from all five simulated GC age groups.

Table 2 .
Parameters of the fitted color CMR (see Equation1) for each stellar model, color index, and GC age.These curves appear in Figure3.

Table 3 .
Hartman et al. (2023)n parameters from Equation 5 for theHartman et al. (2023)data transformed with our nominal color-metallicity relations.These curves appear in Figure9.