Measuring M31 globular cluster ages and metallicities using both photometry and spectroscopy

Globular clusters (GCs) are a powerful tool to study galaxy formation (see reviews by Brodie & Strader 2006 and Forbes et al. 2018) as the properties of GCs observed today reflect both the conditions of their formation but also the physics of GC survival. Found in virtually all galaxies with stellar masses above $10^{9}$ M_⊙ (and most galaxies more massive than $10^{7}$ M_⊙, Harris et al. 2013; Eadie et al. 2022), they are significantly brighter than individual stars, allowing them to be studied at much greater distances. The oldest GCs serve as fossils of the earliest stages of galaxy formation while the continuous formation of GCs until today (e.g. Martocchia et al., 2018; Bastian et al., 2019) allows them to trace the full history of galaxy formation. GC ages provide both an important test of theories of how GCs formed and evolved but also can provide important constraints on how galaxies form and assemble.

In this paper we use the term GC for any star cluster older than $\sim 1$ Gyr and young massive cluster (YMC) for any star cluster younger than this age. This division corresponds to the transition for a massive cluster from being dominated by formation physics to disruption physics and in stellar evolution terms to when the red giant branch starts to make a significant contribution to the luminosity of the cluster.

Two broad classes of models for GCs have been proposed. Motivated by the old ages of Milky Way (MW) GCs (e.g. Janes & Demarque, 1983; VandenBerg et al., 2013; Valcin et al., 2021) and the lack of easily observable star clusters forming today in the MW massive enough to survive a Hubble time, the first class of models (e.g. Peebles, 1984; Trenti et al., 2015) invoke special conditions in the early Universe. In these models GCs form during or even before the epoch of reionisation ($z\gtrsim 6$, lookback times $\gtrsim 13$ Gyr), often in their own dark matter halos. In the second class of models (e.g. Elmegreen & Efremov, 1997; Elmegreen, 2010; Kruijssen, 2015), motivated by observations of star clusters of GC mass or greater forming in starburst galaxies today (e.g. Whitmore & Schweizer, 1995; Holtzman et al., 1996; Adamo et al., 2020), GC formation is the natural outcome of intense star formation. In the first class of models, all GCs in all galaxies should have formed in the first Gyr of the Universe; in the second, the distribution of GC ages should be broader and should vary with galaxy assembly history. Quantitative models of GC formation (e.g. Li & Gnedin, 2014; Pfeffer et al., 2018; Valenzuela et al., 2021; Rodriguez et al., 2023) make predictions for the age distribution of GCs that can be compared to observations.

GC are also powerful probes of the formation and assembly history of galaxies. Since a large star formation rate density is required to form a massive, compact star cluster, the ages of GCs trace periods of intense star formation. If the age and metallicity of a GC are known, the galaxy mass-metallicity relation (see review by Maiolino & Mannucci, 2019) can be inverted to give the likely mass of the galaxy the GC formed in (Kruijssen et al., 2019). With multiple GCs, the mass assembly history of a galaxy can be constrained. The presence of GCs with a significant range of metallicities at the same age implies those GCs formed in different mass galaxies. The presence of at least two branches of the MW GC age-metallicity relationship was used to argue that the older branch was the in situ population and the younger branch was the accreted population (Forbes & Bridges, 2010; Leaman et al., 2013). This association was confirmed by the younger branch GCs having the same orbits as accreted satellite galaxies (e.g. Massari et al., 2019) By folding in constraints on the orbits of GCs, the merger tree of a galaxy can be reconstructed in detail as was done for the MW by Kruijssen et al. (2020) since earlier mergers and more massive mergers should deliver their GCs to smaller galactocentric radii (Pfeffer et al., 2020).

Beyond the Local Group, GCs cannot be resolved into their constituent stars and have to be studied using their integrated light. Despite their importance, measuring reliable GC ages from their integrated light remains a challenge. Optical photometry suffers from strong age-metallicity and age-extinction degeneracies (e.g. Worthey, 1994; Anders et al., 2004; de Meulenaer et al., 2014) in that old ages and higher metallicities both make a stellar population redder as can the effects of interstellar extinction. While UV or NIR photometry can be used to break this degeneracy, obtaining deep enough photometry for a significant sample of GCs can be observationally expensive. Spectroscopy can also provide ages but the quality of spectra required is likewise observationally expensive for significant samples. While spectroscopy does not suffer from the age-extinction degeneracy, it can suffer from a degeneracy between the morphology of the horizontal branch and age; see Cabrera-Ziri & Conroy (2022) for a detailed discussion.

Usher et al. (2019b) used optical photometry and spectroscopy of the near infrared calcium triplet (CaT) to derive ages and metallicities for GCs in three massive early-type galaxies from the SLUGGS survey (Brodie et al., 2014) and the MW. Usher et al. (2019a) showed that the strength of the CaT is a reliable and age-independent measure of metallicity for stellar populations older than a couple Gyr. Usher et al. (2019b) used these metallicities as priors when fitting the photometry with stellar population models to measure ages. In line with previous work (e.g. Usher et al., 2015; Powalka et al., 2016) that found that the GC colour-metallicity or colour-colour relationship varies between galaxies, Usher et al. (2019b) found these four galaxies show widely different distributions of GCs in age-metallicity space, suggesting different assembly histories and disfavouring models where the majority of GCs form in the early Universe.

As the nearest (780 kpc, Conn et al. 2012) massive (stellar mass $\sim 10^{11}$ M_⊙, Tamm et al. 2012) galaxy to the Milky Way, M31 provides an important test of techniques for studying GCs using their integrated light. M31 is close enough that its GCs can be observed in ways that are impractical for more distant galaxies - its GCs can be resolved into their constituent stars using space based observations and high signal-to-noise ratio, high resolution integrated spectroscopy can be obtained in a reasonable amount of time - yet it is distant enough that is relatively straight forward to obtain integrated photometry and spectroscopy of its GCs.

There is a long history of studying the metallicities and ages of GCs in M31 using their integrated light (e.g. van den Bergh, 1969; Burstein et al., 1984; Puzia et al., 2005; Beasley et al., 2005; Galleti et al., 2009; Caldwell et al., 2011; Sakari & Wallerstein, 2022). The superior spatial resolution of the Hubble Space Telescope (HST) allows the star clusters of M31 to be resolved into their individual stars. In most cases the resulting colour-magnitude diagrams only reach the upper red giant branch (e.g. Holland et al., 1997; Rich et al., 2005; Perina et al., 2009) but in the case of B379-G312, the photometry reaches below the main sequence turnoff, allowing the age to be measured directly (Brown et al., 2004; Ma et al., 2010). With shallower photometry, lower limits on the age can be placed and the presence of a blue horizontal branch can be used as evidence for an old stellar population (e.g. Perina et al., 2011). While most M31 GCs seem to be old ($\gtrsim 10$ Gyr) and M31 hosts a population of YMCs ($\lesssim 1$ Gyr e.g. Hodge 1979; Elson & Walterbos 1988; Caldwell et al. 2009; Johnson et al. 2016), the existence of a population of younger GCs (between 1 and 10 Gyr) has been debated (e.g. Burstein et al., 2004; Beasley et al., 2005; Caldwell et al., 2011; Perina et al., 2011).

In this paper we extend the analysis of Usher et al. (2019b) to M31 by using literature spectra described in Section 2 to measure the metallicity of 290 GCs using the strength of the CaT spectral feature (Section 3). In Section 4 we use these metallicities as a prior when deriving their ages from their optical photometry before summarising our work in 5.

We measure the metallicities of the M31 GCs using the strength of the calcium triplet (CaT) rather than relying on the Caldwell et al. (2011) Lick index based metallicities to both maintain commonality with Usher et al. (2019b) and since the strength of the CaT shows little dependence on age (e.g. Vazdekis et al., 2003; Usher et al., 2019a). We measure the strength of the CaT using the technique of Foster et al. (2010) and Usher et al. (2012, 2019a). We use the same Python code as in Usher et al. (2019a) to fit the observed spectra with a linear combination of stellar templates (we use the same templates as Foster et al. 2010; Usher et al. 2012, 2019a), normalise the continuum (using the same parameters as in Usher et al. 2019a) and measure the combined CaT strength of all three lines from the fitted templates using the same index definition (that of Armandroff & Zinn 1988) as in our previous work. We refer the interested reader to Usher et al. (2019a) for details of the measurement process and a detailed discussion of potential systematics. To avoid systematics due to poor sky subtraction or low signal-to-noise ratios (S/N), we only measure spectra with S/N greater than 20 Å^-1. For the Sakari et al. spectra we also excluded a handful of spectra with higher S/N that showed strong sky subtraction residuals.

Unfortunately, the measured strength of the CaT depends on the spectral resolution or velocity dispersion, with the CaT strength being systematically lower at lower spectral resolutions, especially at high metallicity (e.g. Cenarro et al., 2001; Usher et al., 2019a). While the resolution of the Sakari et al. spectra are high enough not to be significantly affected by this effect, the Caldwell et al. (2011) spectra are low enough resolution to be. We used high S/N spectra of the CaT region of 27 GCs from the WAGGS survey Usher et al. (2017, 2019a) which span the range of MW GC metallicities ($-2.4<$ [Fe/H] $<-0.1$) to derive an empirical correction for the effects of velocity dispersion on the measured strength of the CaT. We first smoothed each of the WAGGS spectra to a range of velocity dispersions between 10 and 100 km s^-1 and measured the strength of the CaT on each smoothed spectra. For each WAGGS GC we then fit a cubic spline between the measured velocity dispersion and the measured CaT strength. We used these splines to estimate the CaT strength at a velocity dispersion of 72.3 km s^-1, the median of the velocity dispersions fitted to the Caldwell et al. (2011) spectra, for each of the WAGGS GCs. We then fit a cubic relationship between the CaT strength measured on the spectra broadened to 72.3 km s^-1 and the CaT strength measured on the original, unbroadened spectra. This broadening is consistent with the line spread function of Hectospec ($\sigma=78\pm 4$ km s^-1) found by Fabricant et al. (2013) at 8600 Å once the line spread function ($\sigma=19$ km s^-1) of the DEIMOS template spectra is accounted for. We used this relationship to correct each of the CaT strengths measured from the Caldwell spectra. To estimate the systematic uncertainty introduced by this correction, we repeated the correction at 62.7 km s^-1 and 86.7 km s^-1, the 2.3 and 97.7 percentiles of the velocity dispersion distribution measured from the Caldwell et al. spectra. We used these respectively to correct the upper and lower uncertainties on the measured CaT strength. In terms of metallicity the strength of the correction varies from 0.09 dex at [Fe/H] $=-1$ to 0.42 dex at [Fe/H] $=0$ as seen in Figure 1.

We used the empirical correction of Usher et al. (2019a, their equation 1) to translate our CaT values into metallicities. The Usher et al. (2019a) relation is based on MW and MW satellite GCs. Since CaT is sensitive to the [$\alpha$/Fe] abundance ratio (Brodie et al., 2012; Sakari & Wallerstein, 2016; Usher et al., 2019a), if the [$\alpha$/Fe]-[Fe/H] relation of M31 GCs is different, our [Fe/H] measurements will be biased. In Figure 1, we show a comparison of our CaT based [Fe/H] with the [Fe/H] measured by Caldwell et al. (2011) using Lick Fe indices and an empirical relationship before and after performing the spectral resolution correction. The root mean squared (RMS) difference between the two metallicities is 0.25 dex, which is in agreement with the uncertainties. While the correction for spectral resolution improves the agreement at the highest metallicities, small ($<0.2$ dex) systematic differences remain, with the CaT metallicities being systematically higher near [Fe/H] $\sim-1.5$ and systematically lower at [Fe/H] $\sim-0.5$. We note that both the Caldwell et al. (2011) metallicities and our CaT metallicities are based on simple empirical relationships between the strength of spectral indices and the metallicities of MW GCs. Some level of systematic difference is unsurprising, given that the Fe Lick indices used by Caldwell et al. (2011) and the CaT used in this work show difference dependencies on [$\alpha$/Fe] and the simple empirical relations used by Caldwell et al. (2011, bilinear) and Usher et al. (2019a, linear).

In Figure 2 we show a comparison between the CaT [Fe/H] values measured from the Sakari & Wallerstein (2016) and Caldwell et al. (2011) spectra. The agreement between the two is excellent with a RMS difference of 0.11 dex in agreement with the statistical and expected systematic uncertainties. In Figure 3 we compare our CaT metallicities with those from Sakari & Wallerstein (2022). Again we see good agreement at most metallicities although the Sakari & Wallerstein (2022) metallicities are systematically higher at metallicities below [Fe/H] $=-2$. Our measurement for EXT8 ([Fe/H] $=-2.71\pm 0.07$) is closer to the value from high resolution spectroscopy (Larsen et al. 2020, [Fe/H] $=-2.91\pm 0.04$) than the CaT based value from Sakari & Wallerstein (2022) ([Fe/H] $=-2.27\pm 0.10$). In Figure 4 we compare our metallicities with those based on analysis of high resolution integrated light. We see good agreement with the optical studies of Colucci et al. (2014), Sakari et al. (2015) and Larsen et al. (2022) as well as the near-infrared study of Sakari et al. (2016) using either the Caldwell et al. or Sakari et al. spectra.

To determine ages we used the Monte Carlo Markov Chain (MCMC) code first presented in Usher et al. (2019b) which we now name mcmame (Monte Carlo Metallicity, Age, Mass and Extiction). In this analysis we sample age, metallicity, mass and reddening posteriors of a grid of stellar population models subject to the constraints provided by the photometry and optionally the CaT metallicity. As in Usher et al. (2019b) we use emcee (Foreman-Mackey et al., 2013) with a parallel-tempered ensemble sampler (Vousden et al., 2016) with 8 temperatures, 64 walkers, 200 burn-in steps and 2560 production steps. We use the same stellar populations models (FSPS Conroy et al. 2009; Conroy & Gunn 2010 calculated with the MIST isochrones Choi et al. 2016 and the MILES spectral library Sánchez-Blázquez et al. 2006), and mass priors (flat in log mass), but use slightly different age priors (flat between 1 Myr and 15 Gyr) and metallicity priors (flat between [Fe/H] $=-3.0$ and $+0.5$) than in Usher et al. (2019b). Unlike in Usher et al. (2019b) where the reddening vector from Schlafly & Finkbeiner (2011) was assumed for all ages and metallicities, we used FSPS to calculate the age and metallicity dependent reddening vector using the Cardelli et al. (1989) MW extinction curve. At low extinction values this assumption is valid but for larger values it breaks down since the effect of extinction depends on the shape of the spectral energy distribution within a filter passband, with a population with a bluer spectrum being more effected than a redder one. For reddening we need to adopt different priors than in Usher et al. (2019b) given that M31 has substantial internal extinction, unlike the massive early-type galaxies studied in Usher et al. (2019b).

As a first order approximation, we can think of most of the internal extinction in M31 coming from a thin disc of gas and dust. The majority of old and intermediate age GCs should be either in front of or behind this disc as star clusters would not survive long orbiting in the same plane as massive molecular clouds (Elmegreen, 2010; Kruijssen, 2015). Thus the prior on the extinction should be bimodal, with one peak corresponding to the foreground extinction and one peak corresponding to the extinction of the foreground plus the disc of M31. In principle a GC should have equal probability of being in front or behind of the disc. However, in a magnitude limited sample, GCs on the far side of the disc will be fainter and more likely to fall out of the selection.

The extinction can be modelled through various methods including from far infrared and sub-millimetre observations (e.g. Schlegel et al., 1998; Whitworth et al., 2019) and by modelling the colours of stars (e.g. Dalcanton et al., 2015). Unfortunately, the distribution of gas and dust is filamentary, with significant changes in the extinction value on comparable scales to the radius of a GC (e.g. Bonatto et al., 2013; Saracino et al., 2019; Bonatto & Chies-Santos, 2020). While for M31, extinction can be modelled on the scales of $\sim 30$ pc (e.g. Whitworth et al., 2019; Dalcanton et al., 2015), for more distant galaxies it can only be modelled on larger scales. To mimic the low spatial resolution available for more distant galaxies, we used the Schlegel et al. (1998) extinction maps which have a resolution 6.1 arcmin ($\sim 1.4$ kpc at the distance of M31). We modelled the contribution of the disc of M31 to the extinction as a Gaussian with a mean provided by the value of the Schlegel map at that location and a standard deviation equal to 50 % of the mean. We used the median extinction from the Schlegel map of GCs outside of the disc of M31 (0.08 mag) as the mean value of the foreground contribution and assumed a dispersion of 0.03 mag on the foreground. In our bimodal extinction prior we assigned equal weight to the foreground and M31 disc components. In all cases we enforced a further constraint that the extinction must be positive.

An example of the posterior distribution for B379-G312 is plotted in Figure 5. This presents the best case for sampling the posterior distribution as B379-G312 is bright, metal rich and outside of the projected disc of M31. A more typical case is show in Figure 6, a low metallicity ([Fe/H] $=-1.60$) GC with an apparent $g$-band magnitude equal to the median of our sample ($g=17.25$). In Table 1 we give our age, metallicity, mass and extinction posteriors, the full version of which is provided online as Supporting Information. Our median metallicity uncertainty is 0.10 dex and our median age uncertainty is 2.4 Gyr.

A number of GCs from the Caldwell sample show surprisingly young age posteriors, given that the Caldwell sample of spectra were selected to be older than a couple Gyr on the basis of their spectra between 3750 and 5400 Å (see section 5 of Caldwell et al. 2009). The lack of young GC in both of the Caldwell and Sakari samples is supported by the lack of GC spectra with significant Paschen absorption which is seen in stellar populations younger than $\sim 1$ Gyr (e.g. Usher et al., 2019a). In the case of some GCs such as M009, the age is poorly constrained, with the posterior distribution spanning a wide range of ages ($2.9_{-2.7}^{+8.8}$ Gyr in the case of M009 using the 16 to 84 % percentile confidence interval.). However, a number of GCs are inconsistent with being as old as 5 Gyr at a 95 % level. As can be seen in Figure 7, these GCs with apparently young ages are found at small galactocentric radii and have large priors on their extinctions. Thus these GCs are found in the crowded centre of M31, where the effects of the age-extinction degeneracy are strongest and both photometry and spectroscopy are most likely to suffer from systematic uncertainty.

Using the strength of the near infrared CaT spectral feature we have measured the metallicity of a large number of GCs in M31. In line with previous work (e.g. Sakari & Wallerstein, 2016; Usher et al., 2019a), we find that the CaT is a reliable measure of metallicity although its reliability declines at high metallicity and at the lower spectral resolution of the Caldwell et al. (2009) Hectospec spectra. We used these metallicities as priors when fitting optical photometry with stellar population models to measure the ages, metallicities and masses of the GCs. Due to the strong degeneracy between extinction and age, the ages of GCs projected on to the disc of M31 are less reliable than those in the halo of M31. We find good agreement between our age measurements with most but not all literature studies. Most of the GCs projected on to the disc and virtually all halo GCs are old (consistent with ages $>10.5$ Gyr, formation redshifts $z>2$). The distribution of GC ages and metallicities is similar but not identical to galaxies with a similar stellar mass. The old ages of the majority of M31’s GCs suggest that M31 formed much of its stellar mass early, inline with observations of the field star formation history. In the outer halo we find the most metal rich and youngest GCs are more likely to be associated with substructure rather than the smooth halo, in line with the predictions of Hughes et al. (2019).

The technique used in this paper of combining photometry with a prior from spectroscopy is likely not the best way of measuring the ages of M31 GCs where high S/N optical spectra are available (see Cabrera-Ziri et al. in prep.). Instead it serves as a test of a technique suitable for more distant galaxies. Upcoming highly multiplexed red and near infrared spectrographs, such as MOONS (Cirasuolo et al., 2014; Taylor et al., 2018) on the Very Large Telescope, will allow large numbers of extragalactic GCs to be observed efficiently. Like other techniques that rely on photometry, our method suffers from the age-extinction degeneracy and will perform worse when there is not a strong prior on extinction, such as in or near areas of active star formation. In addition, spectra, even when only available in the red or near infrared, do provide some age information. This technique should be refined by either fitting the spectra with stellar population models and using that age-metallicity posterior as a prior when fitting the photometry or by simultaneously fitting the photometry and spectroscopy with a stellar population model.

We thank the anonymous referee for their helpful comments which improved the manuscript. We also wish to thank Joel Pfeffer for useful discussions. CU acknowledges the support of the Swedish Research Council, Vetenskapsrådet. This study was supported by the Klaus Tschira Foundation. We wish to thank Charli Sakari for making her spectra available.

This work made use of numpy (van der Walt et al., 2011), scipy (Jones et al., 01), matplotlib (Hunter, 2007) and corner (Foreman-Mackey, 2016) as well as astropy, a community-developed core Python package for astronomy (Astropy Collaboration et al., 2013). The CaT fitting code of Usher et al. (2019a) relies on pPXF (Cappellari & Emsellem, 2004; Cappellari, 2017).

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics — Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The Caldwell et al. (2009) spectra are available from https://oirsa.cfa.harvard.edu/signature_program/. The Sakari & Wallerstein (2016), Sakari et al. (2021) and Sakari & Wallerstein (2022) spectra are available upon request from Charli Sakari. The Peacock et al. (2010) and Huxor et al. (2014) photometry was taken from those papers; the remaining SDSS DR17 (Abdurro’uf et al., 2022) photometry is available from http://skyserver.sdss.org/dr17.