THE PANCHROMATIC HUBBLE ANDROMEDA TREASURY. V. AGES AND MASSES OF THE YEAR 1 STELLAR CLUSTERS^*

Synthetic spectral energy distributions (SEDs) of clusters are constructed with the population synthesis code Pégase.2n (M. Fouesneau et al., in preparation), which is derived from Pégase (Fioc & Rocca-Volmerange 1997). As in the original population synthesis code, the underlying stellar evolution tracks are those of the Padova group (Bressan et al. 1993; Bertelli et al. 1994), with a simple extension through the thermally pulsating AGB based on the prescriptions of Groenewegen & de Jong (1993). The input stellar spectra are taken from the library of Lejeune et al. (1997). The stellar initial mass function (IMF) is taken from Kroupa et al. (1993) and extends from 0.1 to 120 M_☉. Nebular emission (lines and continuum) is computed as in Fioc & Rocca-Volmerange (1997) and is included in the calculated spectra and broadband fluxes under the assumption that no ionizing photons escape. The photometry for the synthetic clusters is computed using the response curves of the PHAT HST/ACS and HST/WFC3 filters. A reference spectrum of Vega provides zero-magnitude fluxes (Bohlin 2007).

The code uses MC methods to populate the stellar mass function (SMF) with a finite number of stars (the simulations draw an explicit number of stars instead of a target mass, which avoids the possible biases linked to the latter case; (Kroupa et al. 2013)). As a result, models explicitly account for stochastic variations in the stellar content of clusters. We have extended the simulated cluster set to lower masses than available in Fouesneau et al. (2012), and the coverage now ranges from ∼50 to 5 × 10⁵ M_☉, and ages from 1 Myr to 20 Gyr. With a few × 10⁸ individual models, the collection of synthetic clusters is large enough to include all reasonably likely cluster properties.

We also include a transition to continuous models above 2 × 10⁵ M_☉ (i.e., models assuming a continuously populated SMF), motivated by the presence of about ∼30 massive globular clusters in the PHAT catalog (>10⁵ M_☉; Caldwell et al. 2011). For computational reasons discussed in Fouesneau et al. (2012), the mass distribution in the collection of models follows a power law of index −1. The ages of the synthetic clusters are drawn from a power-law distribution with index −1 (equal numbers per logarithmic bin), rounded to integer multiples of 10⁶ yr. The extinction, A_V, is allowed to vary uniformly from 0 to 3 mag with a fixed R_V of 3.1 assuming the standard extinction law of Cardelli et al. (1989).

Based on H ii region abundances (e.g., Zurita & Bresolin 2012), we expect young clusters in M31 to have approximately solar metallicity. Thus, we fix the metallicity, Z, to 0.02 (solar) for the discrete part of the collection. In contrast, globular clusters are known to have lower metallicities (e.g., Caldwell et al. 2011; Cezario et al. 2013), and thus the metallicity of the continuous models for massive clusters is allowed to vary (equiprobably) between the fixed values of [0.004, 0.008, 0.02, 0.05], corresponding approximately to Small and Large Magellanic Cloud, solar, and super-solar metallicities, respectively.

Individual cluster estimates are collected in Table 2 in Appendix A.

The probabilistic method developed in Fouesneau & Lançon (2010) calculates posterior probability distributions in the age–mass–extinction space (marginalized over the metallicity dimension), using multi-wavelength photometric observations and the large collection of MC simulations of clusters of finite stellar masses, described above. As in all Bayesian inference approaches, the results are stated in probabilistic terms, and they depend on a priori probability distributions of some model parameters (priors).

In our context, the probability for one cluster to have a specific age, mass, extinction, and metallicity, given a set of photometric observations, is

$$\begin{equation} \mathcal {P}(\, \overrightarrow{\theta }\,\mid \,\overrightarrow{d}, \overrightarrow{\sigma } \,) \propto \mathcal {P}(\, \overrightarrow{d}\,\mid \,\overrightarrow{\theta }, \overrightarrow{\sigma } \,) \times \mathcal {P}(\, \overrightarrow{\theta } \,), \end{equation} \tag{ 1 }$$

where the $\overrightarrow{d}$ and $\overrightarrow{\sigma }$ are the six-band photometric measurements and associated uncertainties, $\overrightarrow{\theta }$ is the ensemble of parameter values (i.e., age–mass–extinction–metallicity), $\mathcal {P}(\, \overrightarrow{\theta } \,)$ are the priors on the values of the parameters $\overrightarrow{\theta }$, and $\mathcal {P}(\, \overrightarrow{d}\,\mid \,\overrightarrow{\theta }, \overrightarrow{\sigma } \,)$ is the likelihood of measuring an energy distribution for a given set of $\overrightarrow{\theta }$.

Likelihood. We adopt a normal based-likelihood function, assuming that the photometric errors are Gaussian and independent flux measurements:

$$\begin{eqnarray} &&\mathcal {P}(\, \overrightarrow{d}\,\mid \,\overrightarrow{\theta } \,) = \prod _{k=1}^6 \frac{1}{\sqrt{2\pi }\sigma _k} \exp \left[-\frac{\left(d_k - \hat{d_k}(\overrightarrow{\theta })\right)^2}{2\,\sigma _k^2}\right],\quad\quad \end{eqnarray} \tag{ 2 }$$

in which k indicates properties for the kth filter and $\hat{d_k}(\overrightarrow{\theta })$ the predicted flux in the filter k for the given set of $\overrightarrow{\theta }$.

Priors. The posterior probability distribution, $\mathcal {P}(\, \overrightarrow{\theta }\,\mid \,\overrightarrow{d} \,)$, depends on the prior constraints on the values of the parameters $\overrightarrow{\theta }$. This prior information translates the age and mass distributions of the synthetic clusters, together with the values allowed for extinction. As described in Section 3.1, our priors can be explicitly expressed as an independent combination of age, mass, and extinction functions:

$$\begin{eqnarray} &&\mathcal {P}(\, \overrightarrow{\theta } \,) \propto (M/M_{\odot}) ^ {-1} \times (A/\mathrm{yr}) ^ {-1}, \end{eqnarray} \tag{ 3 }$$

corresponding to independent power laws in age and mass and implicitly uniform in A_V, respectively. The completeness limits of the observations are not used in any manner during the determination of individual ages and masses. The power-law distributions we have adopted mimic two major qualitative trends seen in star-forming galaxies: low-mass clusters are more numerous than high-mass clusters, and, because of a variety of efficient disruption mechanisms, young clusters are more numerous than old ones. As demonstrated in Fouesneau et al. (2012), the priors do not globally dominate the resulting behavior of $\mathcal {P}(\, \overrightarrow{\theta }\,\mid \,\overrightarrow{d}, \overrightarrow{\sigma } \,)$ (Equation (1)). A study based on the full survey will further optimize the prior distributions, in particular, age–mass, to derive cluster disruption efficiencies across the covered area of M31.

An illustrative example of the method for a typical cluster is given in Figure 2. A composite image of the cluster is shown in panel (a). We show the cluster color–magnitude diagram (CMD) in panel (b), in which the gray density map in the background represents the possible contamination by the field stars. On this CMD, we overlaid a blue isochrone corresponding to the best-fit value from the integrated photometry analysis accounting for age–mass–extinction correlations. We represent the correlation between age and mass by the joined age–mass probability distribution function (pdf) in panel (c) of Figure 2, marginalized over the extinction parameter, A_V.

Zoom In Zoom Out Reset image size

In panel (d) of Figure 2, we demonstrate the possible effect of A_V by considering the joined age–mass pdf at two distinct slices of extinction, before integrating over all A_V as in panel (c). We first consider a slice a zero extinction, which is close to the best-fit extinction of 0.3. We represent this pdf in green in panel (d) and the corresponding green isochrone in panel (b). We also show a slice at high extinction (A_V = 2.0) represented in red with a corresponding red isochrone in panel (b). The latter slice contains a very small fraction of the complete pdf, but has been renormalized to allow it to be visible in panel (d).

Figure 2 also demonstrates the subtleties in interpreting pdfs. First, the best value does not always correspond to the peak in a marginalized pdf. The dotted lines indicate the position of the best fit in the age–mass-A_V space, i.e., of the triplet that maximizes the posterior distribution in the full parameter space. Although very close to the peak of the age–mass distribution in panel (c), the best fit suggests a slightly younger and more massive cluster. Second, distributions are complex. The strong correlation between all the different parameters leads to distributions far from Gaussian and sometimes multi-modal.

From this analysis of the integrated light of this particular cluster in the PHAT sample, we obtain a qualitatively good fit, according to the CMD locus of the blue isochrone on panel (b). This panel suggests that the fit indeed accounts for the presence of two bright red stars, a significant number of bright main-sequence stars, and a reddened main sequence. We also expect that the extended distribution toward younger ages (see Figure 2(d)) captures the presence of the two brightest stars on the CMD.

Appendix A gives a table of individual estimates for each cluster of the sample, for the discrete models (including continuous extension at the higher-mass end). The quoted values in this table are the age–mass–extinction triplets that maximize the posterior probability over the full parameter space. Quoted uncertainties are based on the percentiles of the posterior distribution. As we mentioned above, only the full posterior distribution keeps the complexity of the information. While pdfs from this preliminary study are available upon request, full pdfs will be made available with the cluster catalog when the survey will be completed.

We characterized the potential artifacts and biases in the determination of cluster ages, masses, and extinctions using samples from the synthetic cluster collection. Briefly, for a sample of synthetic clusters spanning the full range of cluster ages, masses, and extinctions, we generated six-filter photometry for one sample of synthetic clusters and perturbed the "measured magnitudes" according to uncertainties distributed as in the actual PHAT cluster data. We then re-derived the properties for this synthetic set of clusters and compared the recovered values to the input values.

Overall, we do not find significant age or mass biases from the analysis of synthetic data. The dispersions we obtain are ∼0.14 dex in mass and ∼0.18 dex in age. Those dispersions are consistent with the scatter expected on the basis of the derived pdfs. Based on these tests, we will not venture to interpret features smaller than 0.2 dex in either age or mass (conservatively, since the test data set corresponds to perfectly modeled data). To reflect these limitations, we bin the posterior probability distributions to 0.2 dex in all the subsequent figures.

We find that large errors in age occur for a few percent of the different realizations of the synthetic clusters, mostly at ages of one or a few gigayears. These errors are due to the age–extinction degeneracy, coupled with the addition of the metallicity as a new parameter for the high-mass regime. However, these failures can easily be detected through a visual inspection of the CMDs. Therefore, we used panel (b) in Figure 2 as our baseline to visually inspect the CMDs of each observed cluster and added a caution flag in Table 2 when we estimate a potential failure of the fit. These cases include ∼60 clusters, mainly in the bulge. In further studies, we will include independent determinations of the extinction and/or metallicity when available, providing a better set of initial priors.

In the present study, we estimate ages and masses by analyzing the integrated photometry measurements of the clusters. We could argue that such a method is uncertain and may not give accurate estimates, especially when the energy distributions of the clusters are deeply affected by stochasticity and stars in rare evolutionary phases with unusual colors. Although limited to nearby galaxies where clusters can be resolved into stars, the analysis of CMDs offers more robust estimates.

Taking advantage of the high-resolution HST imaging, PHAT provides individual star photometry for a significant number of clusters. To verify the ages and masses obtained from integrated light, we conducted a CMD analysis of 100 clusters. This CMD analysis will be presented in detail in L. Beerman et al. (in preparation); we briefly summarize the method here. The analysis method is based on Dolphin (2002), which optimizes the likelihood of the observed CMDs with theoretical ones created from stellar evolution tracks, or isochrones, for a variety of ages, mass functions, binary fractions, metallicity, etc., including all observed errors and contamination from a background field population. We selected the sample of 100 clusters for the CMD analysis from various regions of the survey, selecting only clusters with a visible main sequence to ensure the accuracy of the recovered ages. Because of these choices, we restricted the CMD models to have ages less than 10⁹ yr.

Figure 3 presents a one-to-one comparison of age estimates derived from integrated photometry (i.e., the present study) with estimates derived from the CMD analysis. Ellipses are centered on the best-fit values, and their sizes encode the uncertainties from the CMD and integrated light analyses, respectively. While both methods have their own caveats and failure regimes, in general, the CMD fits have smaller uncertainties. Our choice to limit fits to ages younger than 10⁹ yr is likely to be responsible for the few massive and old outliers. Estimated ages are in agreement within the uncertainties. The dispersion around the 1:1 line is symmetric, suggesting that integrated light estimates do not present any significant bias toward younger or older estimates relative to the CMD-derived estimates.

Zoom In Zoom Out Reset image size

In Figure 4, we compare the loci of the observations with a set of unreddened discrete models, in four projections of color–magnitude space. The data are shown together with half of the unreddened models described in Section 3.1, which will be used to assign age, mass, and extinction estimates to each individual cluster. These panels illustrate the dispersions in color and flux that result from the stochasticity inherent to the discrete nature of the IMF. Models cover broad regions of the diagrams, and complex overlap between ages exists and may not be easily visible at first sight. This complexity will only increase with the inclusion of reddened models. Moreover, even though the models plotted in Figure 4 do not include reddening, the majority of the observed clusters lie well within the regions covered by the synthetic clusters. In contrast, continuous population synthesis models (shown by solid lines) are unable to reproduce some of the observations, even when possible reddening is considered. For example, the top left panel of Figure 4 shows a significant fraction of robust measurements (black dots) lying on the left side of the continuous age sequence, i.e., the solid line. Allowing for extinction will not help the continuous models to predict such colors. These colors correspond to those expected for relatively low-mass clusters (< a few × 10⁴ M_☉) that lack any post-main-sequence stars.

Zoom In Zoom Out Reset image size

However, there are also observations even bluer than any models in the F475W − F814W color, lying outside of the ranges of predicted colors and fluxes of both models. We have colored all points with large photometric uncertainties (>0.4 mag) with gray, which shows that these outliers are likely to be due to photometric errors. The UV in particular, F275W and F336W, may still be adversely affected by cosmic-ray artifacts as cautioned in Johnson et al. (2012).

The locations of the clusters in color and magnitude space draw the unsurprising picture of a cluster population spanning a wide range of age and mass. The entire length of the age sequence seems to be populated, and the mass range suggests the presence of a significant number of clusters with masses well below 10³ M_☉. The two bottom panels reveal a general trend that the most luminous clusters in this sample have old ages; these very old clusters fall above the 10⁵ M_☉ curve and are very likely to be old globular clusters. The typical completeness limit of −3 in the F475W band is manifested as the lower limit to the data in the bottom left panel.

Figure 5 shows the age–mass–extinction distributions resulting from the analysis described in Section 3. Each panel on the left-hand side represents the joint probability distribution of two of the variables, after marginalizing over the third. The right panels indicate the best-fit values and their relative accuracy; the ellipses are arbitrarily normalized to reduce the clutter in the figures. The derived ages are distributed between a few Myr and about 10 Gyr, with the few massive candidates appearing only for ages older than 5 Gyr, as expected for the globular cluster population. The derived masses range from a little above the lower limit of our model catalog (50 M_☉) to a bit more than 10⁶ M_☉.

Zoom In Zoom Out Reset image size

The significant number of clusters found to have masses below 10⁴ M_☉ is a strong argument against using traditional continuous models that ignore the discrete nature of the IMF. Multiple studies have shown the impact of stochasticity on the determination of ages and masses of clusters (e.g., Barbaro & Bertelli 1977; Bruzual & Charlot 2003; Cerviño & Luridiana 2004; Maíz Apellániz 2009; Piskunov et al. 2009; Fouesneau & Lançon 2010; Popescu & Hanson 2010). In particular, they identify artifacts that translate into unphysical overdensities at particular ages for young and old clusters, while underestimating the number of moderate age clusters. These artifacts of continuous models are not present in our analysis (right panels of Figure 5), for instance, we do not predict an accumulation of clusters at young ages and instead a relatively smooth age distribution.

Figure 5 shows almost no old and low-mass clusters, nor highly reddened massive old clusters. At old ages, we do not expect many highly reddened clusters to be present, since they should have drifted far from their birth sites. However, even unreddened old clusters are unlikely to be detected unless they were very massive, given that less massive clusters would have faded below the detection limit of the sample. Cluster disruption processes further lower the number of old clusters. Both arguments explain the lack of old, low-mass clusters in the top panel of Figure 5.

There are poor constraints on the number of clusters at the high-mass end, mainly due to their low birthrates. If we suppose the canonical power law with a −2 index, a cluster of 10⁵ M_☉ is 100 times less likely to form than a cluster of 10⁴ M_☉. Although a large galaxy like M31 may have formed a few dozen clusters above 10⁵ M_☉, they are intrinsically rare. Given that we expect to observe all the young massive clusters in the galaxy, the lack of clusters younger than 50 Myr compared to older ages (>1 Gyr) in the age–mass distribution suggests that M31 may have been more efficient at producing massive clusters in the past. We discuss the decrease in the rate of massive cluster production in more detail in Section 4.2.

The bottom panel of Figure 5 shows that very young clusters come with a large range of extinction values, as has been seen in many star-forming galaxies (e.g., Whitmore & Zhang 2002; Kim et al. 2012). However, at ages older than 10⁸ yr, clusters with more than one magnitude of extinction become rare. Between 10⁷ and 10⁸ yr, the absence of reddened clusters most likely reflects a real lack of highly reddened objects, rather than selection effects, given that the detection limits of the data would allow us to detect 5 × 10⁷ yr old clusters with up to more than 2 mag of extinction, for masses above log (M/M_☉) ≈ 4.

Figure 6 shows the inferred dereddened optical colors of the clusters in the sample. These colors are generated by dereddening the clusters according to their best extinction estimate at their best age (and mass) estimates. This two-panel figure distinguishes clusters with typical photometric errors (top panel) from the ones with abnormally large uncertainties, due primarily to uncertainties in background determinations. In both plots we indicate the color–age space covered by the models in blue. The comparison between observations and best-fit values from our models confirms the self-consistency of the analysis: clusters find logical estimates when uncertainties are typical (top panel), while large uncertainties may result in unsatisfactory estimates (bottom panel).

Zoom In Zoom Out Reset image size

The age probability distribution represents the apparent age distribution of the clusters, resulting from the combined underlying cluster formation history, the cluster destruction rate, and observational selection effects. The distinction between the observed age distribution and the history of the cluster formation rate is analogous to the difference between the present day and the IMF.

We can derive the age distribution of the ensemble cluster population. In principle, we derive this distribution by optimizing our prior age distribution, leading to a more complex framework that we will further develop in the context of cluster disruption. In the context of this present work, we approximate this distribution with the co-add of the age pdfs of each individual cluster. This procedure tends to increase the presence of tails in the distributions; however, it produces a more accurate global age distribution than assigning each cluster to a single "best-fit" age, especially in the light of complex and often multi-modal pdfs. The combination of pdfs also preserves the relative quality of the fits between clusters. Moreover, probability distributions are more robust to binning effects.

The composite distribution of ages for the ensemble of clusters is shown in Figure 7 (black solid line). This figure shows two representations of the age distribution of the ensemble of clusters using both pdf representations as a function of logarithmic age. These probability distributions show the relative number of clusters (not mass in clusters) as a function of age. The top panel shows the probability distribution of logarithmic ages, dN/dlog(A), in contrast with the bottom panel showing the distribution of age, dN/dA; $\mathcal {P}(\log (A))$ is defined such that the average number of stars in a given logarithmic interval [log (A), log (A) + dlog (A)] is ${ dN} = N\, \mathcal {P}(\log (A))\, d\log (A)$, while $\mathcal {\mathcal {P}(A)}$ is defined such that the average number of stars found in the linear age interval [A, A + dA] is ${ dN} = N\, \mathcal {P}(A)\, dA$. Hence, an equal number of clusters at every age would appear as a constant horizontal line in the bottom panel and as an exponential function in the top one. In both representations, we also include dotted lines showing the composite distributions when restricted to masses above 10³ M_☉, the estimated completeness limit of our sample over at least a few × 100 Myr.

Zoom In Zoom Out Reset image size

The uncertainties in the age distribution (represented by the shaded region in Figure 7) are dominated by the random sampling of a finite number of a relatively small (∼600) sample of clusters. At old ages (e.g., 5 Gyr) where there are few clusters, including or removing a single cluster will lead to large variations. In contrast, such an alteration of the sample at 100 Myr will have less influence. To estimate these sampling uncertainties, we characterize the variations of the posterior distributions by bootstrapping (Efron 1987; Rubin 1981). Specifically, we make 1000 realizations of the cluster sample, randomly drawing 601 clusters for each but allowing duplications, and re-derive the ensemble age distribution for each realization. The shaded region indicates the range containing 95% of the realizations of the age posterior distribution. Each individual realization is also used to further assess uncertainties when comparing different cluster age distribution models.

The distribution in Figure 7 shows that the present-day ages of clusters span from a few  Myr up to 10 Gyr. The young cluster age distribution (<100 Myr) is relatively flat. Although we note a small decay, the age distribution is statistically consistent with a uniform distribution over the last 100 Myr, as indicated by the thick orange line. If the birthrate of clusters has been relatively steady over this same interval, then the nearly constant age distribution suggests that the cluster disruption processes are likely to be inefficient over ∼100 Myr timescales. We revisit this point fully in Section 5.2. If we apply a mass cut at 10³ M_☉, at the expected mass completeness limit, the resulting cluster age distribution (dotted line in Figure 7) shows a better agreement with a constant rate at younger ages and little change at older ages.

The distribution drops off at ages older than 10⁸ yr, as expected from cluster disruption and observational selection limits (as described in Section 4.1). We fit a power law with spectral index β to the observed age distribution for ages between 10⁸ and 10⁹ yr, using a χ² likelihood statistics:

$$\begin{equation} \mathcal {P}(A/{\rm yr})\propto A^{\beta }. \end{equation} \tag{ 4 }$$

We find that the observed present-day age distribution in this interval can be approximated with a power law of index β = −1.15 ± 0.1. However, note that this power law does not map directly into the cluster formation rate since no cuts have been applied to ensure that the same range of cluster masses is detectable at all ages. The distribution of ages in the older regime is consistent with observations of other galaxies in the literature (e.g., Hunter et al. 2003; de Grijs & Anders 2006; Chandar et al. 2010; Bastian et al. 2012; Fouesneau et al. 2012), which also find power-law present-day age distributions with spectral indices close to −1.

The marginal distribution of cluster masses (i.e., the distribution summed over all ages and extinctions) is shown in Figure 8. We derived the composite distribution of the whole cluster sample from their individual probability distributions as it was done for the age distribution in Section 4.2. The resulting mass distribution is uncertain in the low-mass end (∼10³ M_☉), where the sample is less than 50% complete. We adopt this limit as a lower-mass limit for this distribution. There are also significant uncertainties in the high-mass regime, where the rare presence of one single massive cluster can induce large variations in the distribution. To minimize stochastic sampling of the cluster mass function, we conservatively define a mass upper limit as the mass where the uncertainties calculated from bootstrap resampling (see Section 4.2) are more than 50% of the median value. The mass function is most reliable between these two regimes, which are delimited by the thick vertical lines in Figure 8.

Zoom In Zoom Out Reset image size

Over the interval of 10³–10^5.5 M_☉, the observed present-day distribution of masses for the entire sample can be approximated by a power law,

$$\begin{equation} \mathcal {P}(M/M_{\odot})\propto M^{\alpha }, \end{equation} \tag{ 5 }$$

with index α = −1.73 ± 0.11. Alternatively, if we restrict the sample to only clusters with ages between 10⁷ and 10⁹ yr in order to avoid possible incompleteness, the distribution becomes steeper with an index of −1.89 ± 0.12, over the same mass range, but remains statistically compatible with the fit of the full distribution. Overall, we find that the present-day cluster mass function is well described by a power law, which agrees well with previous analyses of other galaxies (e.g., in Figure 10 of the review from Portegies Zwart et al. 2010). In particular, we find that the younger (<1 Gyr) cluster mass distribution is a close match to a power-law distribution with spectral index −2, in agreement with several other cluster mass function determinations that find indices close to −2 (e.g., Zhang & Fall 1999; McCrady & Graham 2007; Bik et al. 2003; Chandar et al. 2010; Popescu et al. 2012). Variation from the canonical −2 index of the present-day mass function may be the result of stochastic formation of the cluster population, in which a variation of one massive cluster could induce such variation. Variations may also be the result of systematics that we explore in Section 5.

The mass distribution has an apparent peak near 10³ M_☉, which falls off toward lower masses. Based on the MC collection, we know that log (M/M_☉) = 3.1 is the threshold mass under which more than half of the clusters have fluxes below the 50% completeness limit (defined in Section 2) in the F475W photometric passband. This peak is therefore more likely to be due to the incompleteness of the cluster sample at low fluxes than a real feature in the cluster mass distribution. The exact characterization of this peak is uncertain, because (1) being faint, the low-mass objects have larger observational errors than massive ones; (2) clusters that have low masses while remaining above the detection limits must be young, and therefore are in the regime most sensitive to the stochastically sampled cluster models; and (3) the detection of such low-mass clusters becomes inefficient due to their low contrast against the field star background.

The present-day age and mass distributions discussed in Sections 4.2 and 4.3 were derived for the sample as a whole. They therefore include clusters in the bulge (B01), which has a distinct cluster population observed with a different detection efficiency. In this section, we explore potential variations in distributions of cluster properties at different positions within the galaxy to make a first assessment of environmental dependencies in the cluster population.

Figure 9 shows the 2D joint parameter distributions for each of the bricks (i.e., equivalent to Figure 5, but subdivided by regions). Figure 9 highlights the variable completeness of our sample across the galaxy; indeed, it is very difficult to detect low-mass and/or highly extinguished clusters in the bulge, because of the high-luminosity background. In contrast, it becomes straightforward to find them in the outer disk. From the four top panels of Figure 9, one can see that the effective 95% completeness moves ∼1 dex in both age and mass, such that older and lower-mass objects are more easily detected at larger radii. This effect is reflected in the translation to the left of the diagonal limit in the bottom right corner of the upper row of the plots.

Zoom In Zoom Out Reset image size

Figure 9 shows that the difference between the bulge and the other regions is very strong. Most of the massive old clusters are in the bulge, making up the bulk of the globular clusters from our sample. In contrast, the lowest-mass clusters are mainly in the outer regions of M31 (B15 and B21), which have the lowest stellar background density, and thus the best contrast for detecting low-luminosity clusters. Moreover, all of the disk fields are quite similar, beyond the small variations in sensitivity. As a result, for the rest of this paper we will consider only the disk fields (B09, B15, and B21).

Figure 10 compares the marginalized age and mass distributions for each of the three disk regions, following the same conventions as Figures 7 and 8. The black line reproduces the distribution of the whole sample, for reference. There are roughly twice as many clusters in B15 relative to the two other regions, leading to larger uncertainties in the inner (B09) and outer (B21) regions.

Zoom In Zoom Out Reset image size

Age distributions. The top panel of Figure 10 shows that there are no significant radial changes in the local age distributions, when averaged over the scale of the PHAT bricks (1.5 × 3 kpc at the distance of M31). We can see this similarity more clearly in Figure 11, where we overlay the present-day age distribution for all three disk regions.

Zoom In Zoom Out Reset image size

At young ages (<100 Myr), all three age distributions are consistent with a uniform distribution over the last 100 Myr, suggesting that any variations in the cluster formation rate were coherent over the galaxy.

At older ages (>100 Myr), the age distributions in B09, B15, and B21 are all compatible with a power law. The power-law spectral indices of the age distributions from 10⁸ to 10⁹ yr in Figure 10 are −1.39 ± 0.1, −1.11 ± 0.1, and −1.21 ± 0.1 for bricks 9, 15, and 21, respectively, which can be compared to the power law of index −1.15 ± 0.1 derived for the whole sample (Section 4.2). The broad similarity among the age distributions (illustrated in Figure 11 by the orange shaded region) suggests a common cluster formation history within the three different forming regions. In other words, the net results of changes in the cluster formation rate and destruction rate were coherent across the galaxy, in spite of the fact that the disk cluster sample spans a broad range of environments. In particular, we would have expected that variations in the local gas density should lead to different cluster disruption efficiencies (Boutloukos & Lamers 2003; Lamers et al. 2005). However, we do not observe statistically significant differences in the age distributions of clusters from the ring (B15) and the outermost region (B21), where the gas density is the highest and the lowest, respectively. Note that this spatially resolved analysis suggests that the slight deviation from consistency with Figure 7 was primarily due to the inclusion of clusters in the bulge field.

Unfortunately, the current data are not sufficient to distinguish statistically significant variations from brick to brick. While there is a hint of a radial variation in the position of the rollover of the distribution, a reliable interpretation of the rollover is difficult without larger samples and a better characterization of the effects of the completeness of the sample. We therefore postpone this analysis to the full PHAT data set.

Mass distributions. The bottom panels of Figure 10 show the mass distributions of all clusters in the three disk fields. The shaded regions follow the same conventions as in Figure 8. The lower-mass regime is defined to be the region where the sample is less than 50% complete (i.e., M < 10^3.3 M_☉). By excluding this mass region, we remove sensitivity and ensure that we have the same high completeness in all three disk regions. We also set an upper mass limit to be where the uncertainties are more than 50% of the median value (as explained in Section 4.3). The upper mass limit is very different in the 10 kpc star-forming ring (B15) due to the presence of more massive clusters than in the two other regions (B09 and B21). This difference is primarily due to the factor of 2–3 larger number of clusters in brick 15. For stochastic sampling of the low-frequency tail of a power-law distribution, the number of massive clusters in each region is expected to vary much more than this factor of 2–3, and instead should vary by up to a factor of 10 (see Table 1). In other words, while B15 has a greater number of massive clusters, it still has fewer than we would expect to follow the same mass function. We further discuss the lack of massive clusters in Section 5.3.

Qualitatively, B09 and B21 share a similar present-day mass distribution within the mass range where the estimated distributions are reliable. At first glance, however, the mass distribution in B15 appears to be different, although this difference may partially reflect the much larger mass range that can be probed reliably in brick 15.

To quantify the differences among the bricks, we have fit each distribution with a power-law distribution. Based on the bottom row of Figure 10, a fit to the mass distribution of brick 15 may be more consistent with the two other regions if we limit the mass interval to a common range. Therefore, we restrict the fits to both the reliable mass ranges (see Section 4.3) and a common mass range (3.2 < log (M/M_☉) < 4.2), while accounting for their associated uncertainties. We repeat this exercise using bootstrap resampling of the cluster sample in each brick. The top panel of Figure 12 shows the resulting mass spectral index probability distributions for each individual brick. We obtain broad distributions for α in B09 and B21, as expected from the limited mass range used during the fit and the smaller number of clusters. The dashed green line shows the fit of brick 15 mass distribution when including the full reliable mass interval, while the solid version shows the resulting distribution when restricting B15 to the same mass range as B09 and B21, of 3.2 < log (M/M_☉) < 4.2. We find a relatively narrow distribution for B15 given the large number of clusters, and as expected, limiting the mass interval favors steeper mass functions by ∼0.2 dex. When fitting over the common mass range, we obtain the following values and standard deviations: −2.2 ± 0.17 for B09, −1.8 ± 0.1 for B15, and −2.1 ± 0.14 for B21 (−2.24 ± 0.1 for the full reliable interval). Although all of the three distributions of α are compatible (within 2σ) with a single power law of index of −1.8 (found in Section 4.3), the differences from this overall description are close to 3σ.

Zoom In Zoom Out Reset image size

Given that both B09 and B21 share a similar environment and that they also seem to follow a similar age-mass distribution, we combined those samples into a super-region probing spiral arms. This combination increases the number of clusters in the statistics, allowing us to compare two distinct types of star-forming regions in M31: spiral arms and the ring. The bottom panel of Figure 12 compares the distributions for B15 from the top panel with the distribution obtained for the combination of B09 and B21 cluster samples. As expected, the combination of the two regions narrows down the spectral index dispersion. However, the resulting distribution does not reconcile with what we find in B15 (best fit of −1.8 ± 0.12), and the resulting discrepancy is statistically increased. As the PHAT survey will eventually cover inter-regions and include many more clusters, we will be able to fully characterize the mass variations with environmental conditions.

We have compared the cluster populations in three star-forming regions in M31, and we may have identified the first clear evidence of variation of the present-day cluster mass function within one galaxy. On the other hand, we find very similar age distributions in the same three regions. We now explore the possible observational artifacts that could be affecting our analysis.

Systematics from the analysis method. We have made multiple assumptions when deriving ages and masses for the cluster samples. Some assumptions could potentially lead to systematic errors in the age and mass determinations.

First, we assumed in the models that no ionizing photons escape from the cluster. If some photons do escape, then we will overestimate the flux in the nebular emission (continuum + lines), and thus, the model colors will vary. These color variations will shift clusters to different apparent ages, with little change in the inferred masses. Moreover, these effects are in the opposite sense of what is needed to explain the data where we find very similar age distributions but differences in the observed mass functions.

Second, there may be some level of inconsistency between the stellar evolution models and the actual clusters, which could lead to biases in the derived ages and masses. However, such an effect would be apparent in all regions and would not produce radial variations.

Third, the lower-mass limit of the current collection is sufficiently low that the choice of the stellar IMF or its sampling method can affect the derived SEDs. If real clusters have an SMF different from that assumed, then the cluster ages and masses would be biased. Again this would not produce a radial variation in the mass distribution of the clusters, unless the stellar IMF were also environmentally dependent.

Finally, our choice of priors during the analysis of cluster colors may not be optimal. We assumed uniform expectation in age and mass on logarithmic scales, but our tests have shown that varying the prior assumptions does not significantly affect the resulting distributions. Therefore, our choice of prior is unlikely to produce variations in the cluster mass distribution across the galaxy.

Apart from an environmental variation of the initial SMF, none of the above possibilities appear likely to produce the observed radial variations in the cluster mass function while keeping similar age distributions. All seem likely to affect all of the regions in a similar manner (outside of the bulge).

Variations in the fraction of bound clusters. Our analysis is based on the cluster sample from the Johnson et al. (2012) catalog. Like any catalog, the resulting sample has biases that reflect how the clusters were selected. Johnson et al. (2012) adopted a definition of a "star cluster" to be a group of stars assumed to form a coeval population (i.e., single age, metallicity, etc.). This definition includes any clustered stars regardless of whether they are gravitationally bound or not. It is therefore possible that the observed variations are due to radial changes in the relative numbers of gravitationally bound clusters and unbound associations. Associations are most likely young because of their intrinsic instability (Gieles & Portegies Zwart 2011). When their size is comparable to stellar clusters, they are also likely to be relatively low mass. Therefore, increasing the proportion of associations will lead to a higher fraction of clusters with young ages and low masses.

However, we do not observe a significant change at young ages between the three different regions. Including many low-mass objects could result in a steeper mass function. If there is a systematic radial change in the fraction of associations, we expect this fraction to be the highest in B21, which is the ideal environment for finding faint objects because of its low background and relatively low density of sources and extinction. However, as shown in Section 4.4, the present-day mass function of B21 is similar to that of B09, which is the innermost region used in this comparison. When combined with the fact that we see no obvious variations in the observed age distributions at recent times, it seems highly unlikely that variations in the fraction of bound clusters are systematically affecting the observed age and mass distributions.

Completeness variations and missing low-mass clusters. The completeness of the PHAT cluster sample is a complex function of six-filter photometry affected by measurement errors, which thus does not translate into a single mass cut in the age–mass plane. Moreover, clusters in this sample are partially resolved into stars, which increases the complexity of estimating the completeness accurately.

From the observational limits given in Johnson et al. (2012), we derived an approximate mass completeness of ∼10³ M_☉ using our collection of synthetic models.

Instead of deriving a complex completeness function, we instead applied a conservative mass cut of 10^3.2 M_☉, based on combining the observational limits given in Johnson et al. (2012) with our collection of synthetic models. Above this mass, we expect the sample to be essentially 100% complete, except perhaps at the very highest extinctions. This expectation is born out by the distributions in Figure 9, which show no obvious signs of incompleteness at the 10^3.2 M_☉ level.

Brick 15 is the most gas rich of the three star-forming regions. It is possible that B15 contains more dust, and this results in higher incompleteness due to dust extinction. We can estimate the number of potentially "missing" clusters in brick 15 by deriving samples to make an intrinsic mass power law with index α = −2 appear to have a distribution of α = −1.8 (as the median value in brick 15). We find that brick 15 would need to have a factor of ∼5 more clusters between 10³ and 10⁵ M_☉ to recover a spectral index of −2. This corresponds to having 20 missing clusters above 10⁴ M_☉, which were undetected because of more than 2 mag of extinction. We find this possibility to be unlikely, given that the A_V distribution of detected clusters falls off steadily toward high extinctions. Moreover, we have visually inspected all F160W images for embedded clusters at the locations of dense molecular cloud from high-resolution CARMA maps of B15 (A. Schruba et al., in preparation), and we find no evidence for highly embedded massive clusters.

In this study, we find that the present-day cluster age distribution, dN/dt (bottom panel of Figure 7), is globally flat over the last 100 Myr, particularly where confined to the disk fields (Figure 11). There are other cluster samples from the literature that show the same constant distribution at young ages (e.g., Lamers et al. 2005, limited to 600 pc from the solar neighborhood; Hodge 1987; Chiosi et al. 2006, flat over 1 Gyr in the SMC). After this initial 100 Myr period, we find a power-law drop-off at older ages consistent with index −1.15. This index may change when we eventually include the complete sample from the PHAT survey, which will provide greater weight to regions that lack recent star formation than the Year 1 Johnson et al. (2012) catalog.

The interpretation of the observed dN/dt distributions is not straightforward, because these distributions are heavily dependent on the observational completeness of the sample, the birthrate of the clusters as a function of time, and the efficiency of the cluster disruption processes as a function of time, mass, and environment. We are therefore deferring a full analysis of the age distribution until a subsequent paper, when we will have a larger number of clusters in the sample and a better characterization of its completeness.

As an intermediate step, however, we can draw an initial comparison between the observed present-day age distribution and that expected under a few assumptions. If one assumes that the cluster population formed at a constant rate and with a power-law initial cluster mass function (ICMF), and no dissolution, Gieles & Bastian (2008) demonstrate that the age distribution can be analytically estimated. Based on continuous population synthesis flux predictions for cluster fading, they show that the observed age distribution should follow a power-law distribution with an index of ∼ − 0.7, if the sample is limited by one optical band detection, and dN/dt should be constant, if the sample is mass limited.

We find a uniform distribution for the first 100 Myr. If we assume that the young clusters from Johnson et al. (2012) are complete down to our 10^3.2 M_☉ mass limit, then we should obtain a constant distribution until the age that the fading or disruption starts to remove clusters from the sample.

We find that the drop-off occurs at ∼100 Myr, independent of the region of the galaxy, suggesting that cluster disruption has little to no effect prior to this timescale. In addition, the 100 Myr extent of the flat dN/dt distribution appears to be the same in all regions, suggesting that the environmental dependence of the cluster disruption must be weak.

At older ages (>100 Myr), we find a power-law decrease in the observed dN/dt distribution, with an index of β = −1.15. This index is steeper than the predicted value by the cluster fading model of Gieles & Bastian (2008) in the presence of a magnitude limit. This may suggest that there are some cluster disruption effects at work above the 10^3.2 M_☉ mass limit of our analysis. On the other hand, the role of selection effects has not yet been fully qualified and can potentially lead to the steepening that is observed. We will work on this issue more fully in an upcoming paper.

The reliable mass range over which we can fit the present-day mass distribution with a power-law function varies from one region to another. We adopted different limits motivated by the lack of massive clusters in B09 and B21 in contrast with B15. This difference could simply reflect the smaller number of clusters overall in B09 and B21, or it could reflect that the cluster mass function in these bricks is "truncated," given that other studies have found evidence for a truncation of this power law at the high-mass end (e.g., Bastian & Gieles 2008; Larsen 2009; Vansevičius et al. 2009, in M31 for the latter). To test for a possible truncation, we can estimate the most massive cluster we would expect from an untruncated power-law mass function and compare to our observations.

To do this test, we assume that the cluster formation rate is constant and that the cluster mass function follows a power-law mass distribution over a mass range [M₁, M₂],

$$\begin{eqnarray} \mathcal {P}(M/M_{\odot}) &\propto & M^{\alpha }, \end{eqnarray} \tag{ 6 }$$

where α is the spectral index of the cluster mass function. With these assumptions, we can define the probability of obtaining a cluster with a mass m above a given mass M as

$$\begin{eqnarray} \mathcal {P}(m \ge M \mid \alpha) & = & \frac{\int _{M}^{M_2}{ x^\alpha \,dx}}{\int _{M_1}^{M_2}{ x^\alpha \,dx}}. \end{eqnarray} \tag{ 7 }$$

We need then to introduce the total number of clusters, N_cl, formed from the mass distribution. If we consider independent draws from this mass distribution, then the expected maximum mass M_max satisfies the condition that we have one and only one cluster for m = M_max; therefore,

$$\begin{eqnarray} 1 &=& N_{cl} \times \mathcal {P}(m > M_{{\rm max}} \mid \alpha) \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} M_{{\rm max}} &=& \left\lbrace M_2^{\alpha +1} - \frac{1}{N_{cl}} \times \left(M_2^{\alpha +1} - M_1^{\alpha +1}\right) \right\rbrace ^ {-(\alpha +1)}.\quad\quad \end{eqnarray} \tag{ 9 }$$

If we consider the upper mass limit M₂ to be infinite (i.e., no truncation), and α < −1 for convergence, then

$$\begin{eqnarray} &&M_{{\rm max}} = {M_1} \times N_{cl} ^ {\alpha +1}. \end{eqnarray} \tag{ 10 }$$

Figure 13 shows the different expected distributions of M_max as a function of the power-law index α, for each field brick in our study. The different curves are constructed from Equation (8), assuming that the mass function is defined within 10³–10⁷ M_☉, to match our typical completeness limits. The number of clusters N_cl above 10³ M_☉ is taken from Table 1 for each region. We compare the predicted maximum cluster mass with the observed values derived from the data in Section 4.4. Because the masses were derived from pdfs rather than single best-fit values, our estimated M_max is taken to be the point, where the uncertainties on the mass distribution become more than 50% of the estimated value.

Zoom In Zoom Out Reset image size

In Figure 13, the values for each brick are indicated by the black crosses, for which the spectral indices correspond to the best values obtained in Section 4.4.

The comparison between the predicted and observed values of M_max shows that the lack of massive clusters in B09 and B21 is consistent with the expectation for a stochastic sampling of a power-law mass distribution. In other words, there is no need to invoke mass truncation to explain the observed upper mass limits of the present-day cluster mass function. The mass cut we observe in B15 (M_max = 1.7 × 10⁵ M_☉) appears to be a factor of 2 smaller than the theoretical sampling prediction of the maximum mass (M_max = 5.2 × 10⁵ M_☉). There are only three to four clusters between the adopted M_max and the theoretical prediction value, which leads to large uncertainties at this regime. We therefore have adopted our more conservative value of M_max.