THE HIGH-MASS STELLAR INITIAL MASS FUNCTION IN M31 CLUSTERS^*

Our best constraints on the high-mass IMF come from studies in the Milky Way (MW) and nearby galaxies. For the closest galaxies, it is possible to count individual stars and make a direct measurement of the high-mass IMF slope. While measuring the IMF from direct star counts is straightforward in principle, it is far more complicated in practice. Accurate measurements are challenging due to a variety of observational and physical effects including completeness, dynamical evolution, distance/extinction ambiguity, stellar multiplicity, the accuracy of stellar evolution models, degeneracies between star formation history (SFH) and the IMF slope, and a paucity of massive stars in the nearby universe, all of which complicate the translation of the present day mass function (MF) into the IMF (e.g., Lequeux 1979; Miller & Scalo 1979; Scalo 1986; Mateo et al. 1993; Massey et al. 1995a; de La Fuente Marcos 1997; Maíz Apellániz & Úbeda 2005; Elmegreen & Scalo 2006; Maíz Apellániz 2008; Ascenso et al. 2009; Kruijssen 2009; Bastian et al. 2010; De Marchi et al. 2010; Portegies Zwart et al. 2010; Kroupa et al. 2013; Schneider et al. 2015).

Observations of more distant, unresolved galaxies alleviate some of these challenges. For example, by studying an entire galaxy, effects of stellar dynamics (e.g., migration, mass segregation) and distance/extinction confusion are less important. Further, distant galaxies offer far more diversity than what is available in the very local universe and provide better opportunities to explore the high-mass IMF as a function of environment (metallicity, mass, redshift, etc.; e.g., Baldry & Glazebrook 2003; Fardal et al. 2007; Davé et al. 2011; Narayanan & Davé 2013).

However, the unresolved nature of these observations severely limits the precision and accuracy of any high-mass IMF measurement. Most integrated observations provide insufficient leverage to mitigate the SFH-IMF degeneracy, and we are forced to make simplifying assumptions about a galaxy's recent SFH (e.g., constant, single short burst) and dust (e.g., a single value for an entire galaxy) to provide any constraint on the high-mass IMF (e.g., Miller & Scalo 1979; Elmegreen & Scalo 2006). Although studies of unresolved clusters may provide a new avenue for IMF studies (e.g., Calzetti et al. 2010; Andrews et al. 2013; although see Weidner et al. 2014 for a contrasting perspective), high-mass IMF studies using integrated light typically provide coarse characterizations rather than precise measurements.

As a result of these challenges, our current knowledge of the high-mass IMF remains remarkably poor. The widespread adoption of a "universal" Salpeter (${\rm \Gamma }\;=\;+1.35$; Salpeter 1955) or Kroupa (${\rm \Gamma }\;=\;+1.3$; Kroupa 2001) IMF above 1 M_⊙, presumed not to vary with environment has given the impression that the high-mass IMF is a solved problem (e.g., Kennicutt 1998a). However, compilations of IMF studies in the literature from resolved star counts indicate an "average value" of ${\rm \Gamma }=1.7$ with a scatter $\gtrsim 0.5$ dex in IMF slope measurements that does not clearly support a "universal" IMF (e.g., Scalo 1986, 1998; Kroupa 2002; Scalo 2005).¹³ There also does not appear to be any unambiguous systematic dependence of IMF slope on environment and/or initial star-forming conditions (e.g., Bastian et al. 2010), which has significant implications for how stars form in a wide variety of physical conditions (e.g., Elmegreen 1999; McKee & Ostriker 2007; Hopkins 2013; Krumholz 2014).

Controversy over the form and universality of the IMF have persisted since the first studies in the 1950s. Following the seminal work of Salpeter (1955), a handful of studies reported similar IMFs in a number of MW clusters and field environments (Jaschek & Jaschek 1957; Sandage 1957; van den Bergh 1957), prompting early speculation about the universality of the "original MF" (e.g., Limber 1960; van den Bergh & Sher 1960). However, this early agreement with Salpeter gave way to large scatter in IMF measurements beginning in the 1970s, as several studies with more sophisticated completeness corrections and updated stellar physics highlighted the shortcomings in the first generation of IMF studies (e.g., Taff 1974; Lequeux 1979; Tarrab 1982). Furthermore, Miller & Scalo (1979) provided a comprehensive reassessment of the MW field star IMF, in light of the previously unaccounted for SFH-IMF degeneracy. They suggested a high-mass IMF of the with ${\rm \Gamma }=1.5$ for stars 1 ≲ M/M_⊙ ≲ 10 and ${\rm \Gamma }=2.3$ for stars M/M_⊙ ≳ 10. However, Kennicutt (1983) showed that such a steep IMF slope above ∼10 M_⊙ was incompatible with observations of $B-V$ colors and Hα equivalent widths in nearby disk galaxies, suggesting that extrapolating the value of ${\rm \Gamma }=1.5$ from 1–100 M_⊙provided a better match to the data. This conclusion was further solidified in Kennicutt et al. (1994).

The number of observational high-mass IMF studies peaked between the late 1980s and early 2000s. This golden age followed the comprehensive IMF review by Scalo (1986), which concluded that systematic uncertainties were the dominant source of confusion in our poor knowledge of the IMF. Shortly after, the proliferation of CCDs enabled significant progress in reducing these systematics, primarily by allowing for quantitative completeness corrections via artificial star tests (ASTs; e.g., Stetson 1987; Mateo et al. 1993). The 1990s saw a flurry of qualitatively new IMF studies including star counts in the LMC, spectroscopy of starburst galaxies, and IMF constraints from stellar abundance patterns (see Leitherer 1998; Massey 1998; Gilmore 2001 for extensive discussion of the literature in these areas). The sheer number and richness of IMF studies from this era precludes a detailed history in this paper, and we instead refer readers to a number of excellent papers that chronicle the history of the IMF (e.g., Scalo 1986, 2005; Massey 1998, 2003; Kennicutt 1998a; Elmegreen 1999, 2006; Eisenhauer 2001; Gilmore 2001; Zinnecker 2005; Bastian et al. 2010; Kroupa et al. 2013).

To briefly summarize, many studies from this era concluded that the high-mass IMF in the local universe was not drastically different from Salpeter over a wide range of environments (e.g., Leitherer 1998; Massey 1998, 2003; Gilmore 2001; Kroupa 2001). However, there was no clear explanation for the observed scatter in the measured IMF slopes, which contradicts a truly universal high-mass IMF. Some have argued that measurement uncertainty and systematics are the dominant source of scatter (e.g., Elmegreen 1999; Kroupa 2001, 2002; Massey 2003). However, without systematically revisiting each study, it is challenging to dismiss the scatter as solely a measurement uncertainty. As a clear counterexample, there are well-known studies of star clusters with some of the same physical properties (e.g., age, mass) that have been analyzed identically, but have very different high-mass IMF slopes (e.g., Phelps & Janes 1993), reinforcing that some of the scatter may be physical in nature. In a broader analysis of the literature, (Weisz et al. 2013) show that >75% of published resolved star IMF studies significantly underreport uncertainties relative to the theoretically possible precision, which underlines the difficulty in disentangling systematics from intrinsic IMF variation using IMF studies in the literature.

Since the early 2000s, studies of the high-mass IMF have diminished in number, as attention turned to the sub-solar regime (e.g., Chabrier 2003; Scalo 2005; Zinnecker 2005). However, focus on the high-mass IMF was recently rekindled following several reports of a systematic deficit of massive stars in lower-luminosity galaxies, based on observations of integrated broadband and/or Hα-to-UV flux ratios (e.g., Hoversten & Glazebrook 2008; Boselli et al. 2009; Lee et al. 2009; Meurer et al. 2009). Yet, whether these observations are indicative of a systematically varying IMF (e.g., Weidner & Kroupa 2005; Pflamm-Altenburg & Kroupa 2009; Kroupa et al. 2013) remains an open question, as other mechanisms such as stochastic sampling of the cluster MF and/or bursty SFHs can also explain these observed properties (e.g., Fumagalli et al. 2011; da Silva et al. 2012; Weisz et al. 2012). Other recent studies have challenged the universality of the high-mass IMF (e.g., Dib 2014), but broader conclusions are limited by the small sample sizes.

After decades of study, the most commonly used high-mass IMF is a universal Salpeter/Kroupa IMF that is perfectly known (i.e., uncertainties are typically not included). However, this de facto standard IMF is not unambiguously supported by a large number of literature IMF studies. This tension is the direct result of a continued reliance on a rich, but heterogenous set of literature IMF studies, and can only be resolved through a large, homogenous survey of individually resolved young, massive stars.

A primary goal of the Panchromatic Hubble Andromeda Treasury program (PHAT; Dalcanton et al. 2012) is a comprehensive exploration of the stellar IMF above ∼2–3 M_⊙. This 828-orbit Hubble Space Telescope (HST) multi-cycle program acquired near-UV through the near-IR imaging of ∼120 million resolved stars in ∼1/3 of M31's star-forming disk (Williams et al. 2014), providing a dataset of unprecedented size and homogeneity for measuring the high-mass IMF.

The focus of the present IMF study is on M31's young star cluster population. Star clusters provide a more direct path to measuring the IMF compared to field populations, which suffer from IMF-SFH degeneracies for stars more luminous than the oldest main sequence turnoff (e.g., Miller & Scalo 1979; Elmegreen & Scalo 2006); furthermore, the differential extinction found in the disk of M31 (e.g., Draine et al. 2013; J. J. Dalcanton 2015, in preparation), and similarly sized star-forming galaxies, adds a significant layer of complexity to modeling the field stellar populations (e.g., Lewis et al. 2015). In contrast, star clusters provide a relatively simple way to probe the IMF. Their coeval nature genuinely mitigates the effects of SFH, while their vulnerability to foreground differential extinction is minimal due to their small sizes.

In this paper, we focus on determining the MF slopes for PHAT clusters that are $\lesssim 25$ Myr old and ≳10³ M_⊙ in mass using resolved star color–magnitude diagrams (CMDs). For these clusters, the detectable portion of the main sequence is well-populated from ∼2 M_⊙ to ≳25 M_⊙¹⁴ and they are of order ∼0.1 relaxation times old, which mitigates the strongest effects of mass segregation, making them an ideal sample with which to study the high-mass IMF. We model the optical CMD of each cluster k to constrain its power-law slope ${{\rm \Gamma }}_{k}$ of the present day high-mass MF, $p({{\rm \Gamma }}_{k})$, marginalized over other parameters of that fit (e.g., cluster age, extinction, binarity).

Subsequently, we propose a simple model for the ensemble properties of the cluster MF, with a mean, Γ, and a scatter, ${\sigma }_{{\rm \Gamma }}$, and fit it to the MF probability distribution functions (PDFs) of the individual clusters. We also use this model to search for possible dependencies of Γ on cluster mass, age, and size. We verify this approach with an extensive set of artificial clusters that were inserted into PHAT images, processed, and analyzed identically to the real clusters.

This paper is organized as follows. In Section 2, we define our cluster sample and briefly describe the PHAT photometry, ASTs, and artificial cluster sample. We outline our method of modeling the cluster CMDs in Section 3 and present the results for the ensembles of real and artificial clusters in Section 4. Finally, in Section 5, we discuss our findings in the context of current and historical views on the high-mass IMF and illustrate the implications of our results.

Clusters in PHAT were identified as part of the Andromeda Project,¹⁵ a collaborative effort between PHAT and the Zooniverse¹⁶ citizen science platform. As described in Johnson et al. (2015), cluster candidates were assigned probabilities of actually being a cluster based on citizen scientist ranking. These probabilities were verified using synthetic cluster results and comparison with the preliminary catalog, which was assembled by professional astronomers (Johnson et al. 2012). Extensive details of the cluster identification methodology can be found in Johnson et al. (2015).

For the present study, we have limited our analysis to the 85 clusters with CMD-based ages $\lesssim 25$ Myr and masses $\gtrsim {10}^{3}$ M_⊙ (L. C. Beerman 2015, in preparation), and half-light radii ≳2 pc (Johnson et al. 2015). We focus on these young, more massive clusters in PHAT to avoid the very sparse sampling in the low-mass clusters; to limit the impact of stellar crowding in the densest clusters; to facilitate the extrapolation of the observable MF to an IMF; and to remain in the regime of high completeness (>90% for our sample; see Figure 10 in Johnson et al. 2015). We will explore each of these issues in detail in our comprehensive study of PHAT cluster MFs in a future paper.

Photometry of PHAT clusters was performed using DOLPHOT,¹⁷ a version of HSTPHOT (Dolphin 2000) with HST specific modules, following the description in Williams et al. (2014).¹⁸ For reasons discussed below, we made slight changes to the survey-wide photometry in that we re-reduced the entire survey using only the optical images (F475W, similar to SDSS g; and F814W, similar to I-band). The optical image and CMD of a typical cluster used in this study are shown in Figures 1 and 2.

Zoom In Zoom Out Reset image size

For the purposes of this paper, we only use optical photometry, which contains the most information about the IMF slope, particularly when modulated by observational and computational considerations. In principle, the full spectral energy distribution (SED) of each star can increase the degree of precision to which we know its mass, and, in turn provide improved constraints on the IMF. However, the UV and IR photometry in PHAT clusters typically contains fewer than ∼50% the number of stars detected in the optical, reducing their statistical utility for constraining the MF slope (see Weisz et al. 2013). PHAT UV imaging is significantly shallower (∼2 mag) than the optical, even along the main sequence, resulting in many fewer sources detected. The IR also yields far fewer main sequence stars as a result of significant crowding due to the lower angular resolution of the WFC3/IR camera.

There are also subtle systematics associated with non-optical bands. For example, the UV and IR bolometric corrections are not as certain as those in the optical, particularly for massive, metal-rich stars, which can introduce systematics into our analysis.

Finally, the modest gain (at best) in precision in the MF slope of each cluster must be balanced against the computational cost of processing >10 million additional ASTs. Thus, by only reducing and analyzing the optical data, we retain virtually all statistical leverage on measuring the MF slope, but save considerably on computational costs.

However, analysis of the full stellar SEDs is valuable for searching for the most massive stars in clusters. In addition to the slope, the maximum stellar mass plays a fundamental role in defining the high-mass IMF, but is unconstrained by the approach we use to measure the MF slope (e.g., Weisz et al. 2013). An investigation of the most massive stars in the PHAT area is the subject on an ongoing study and will be presented in a future paper.

For each cluster, we ran ∼50,000 ASTs by adding stars of known luminosities one at a time (to avoid artificially enhancing crowding effects) into the PHAT optical images of each cluster. The input ASTs were uniformly distributed over each cluster's CMD and spatially distributed approximately according to each cluster's light profile. Extensive testing showed that such a spatial scheme was necessary to accurately characterize the completeness and photometric uncertainty profiles of the clusters, as a simple uniform spatial scheme results in overly optimistic completeness and photometric uncertainty profiles. Similarly, our adoption of a uniform distribution of ASTs in CMD space (as opposed to only matching a given cluster's CMD) is necessary to evaluate all possible cluster models (i.e., all combinations of age, mass, metallicity, membership, IMF). In contrast, only using ASTs that roughly follow a cluster's observed CMD would result in a poor noise model for less obvious possible combinations of cluster properties (e.g., stochastic sampling of the IMF can be degenerate with cluster age), inhibiting a comprehensive search of likelihood space. Finally, we tested the effects of the number of ASTs (ranging from 10,000 to >100,000) used on the recovery of cluster parameters. We found that 50,000 ASTs provided an optimal balance between stability of the solutions and computational efficiency.

We use extensive sets of artificial clusters to verify the accuracy of each component of our MF determinations, including photometry, ASTs, and CMD modeling. The artificial clusters were designed to match the physical properties of our real cluster sample and span a range of physical parameters (age, mass, size, extinction) that encompass the real clusters. The artificial clusters were inserted into images at a range of galactocentric radii and environments to capture the full range of background densities in PHAT. Overall, the artificial cluster tests confirmed that we can accurately recover the true MF from clusters in PHAT imaging, along with all other cluster physical properties.

To measure the MF of a cluster we construct synthetic models of its optical CMD using MATCH, a CMD fitting program described in Dolphin (2002). Briefly, for a given stellar evolution library, binary fraction, IMF, age, metallicity, extinction, and distance, MATCH constructs a simple stellar population (SSP) in the desired filter combinations. This synthetic SSP is then convolved with completeness and photometric biases and uncertainties as determined by the ASTs. The resulting synthetic CMD is then linearly combined with a suitable CMD of foreground/background populations, with the relative contribution of each controlled by a scale parameter. For a given set of these parameters, the synthetic Hess diagram, the binned number density of stars in the CMD, is compared to the observed CMD using a Poisson likelihood function, which is necessary to account for the effects of sparse sampling. The procedure is repeated over a large grid of parameters in order to map all likelihood space. For this analysis, we adopted the Padova stellar evolution models (Girardi et al. 2002, 2010; Marigo et al. 2008) and we have listed the other model parameters and their ranges in Table 1.

Of these parameters, we limited the metallicity to a fairly small range near Z${}_{\odot }$. This was done because (a) the cluster CMDs are too sparsely populated to provide robust metallicity constraints and (b) the current metallicity of M31's disk is known to be approximately solar, with a very weak gradient based on H ii region and supergiant abundances (e.g., Venn et al. 2000; Trundle et al. 2002; Sanders et al. 2012; Zurita & Bresolin 2012).

Similarly, we adopted a fixed binary fraction of 0.35, where the mass ratio of the stars is drawn from a uniform distribution. Historically, stellar multiplicity has been problematic for high-mass IMF determinations from luminosity functions (e.g., Maíz Apellániz & Úbeda 2005; Maíz Apellániz 2008). However, for CMD analysis, it is much less of an issue for two reasons. First, the addition of color information provides clear separation between the single and binary star sequences, minimizing confusion between the two. In contrast, when only luminosity information is available, the binary fraction and IMF slope are essentially degenerate parameters. Second, stellar multiplicity is observationally less important for high-mass stars compared to sub-solar mass stars (e.g., Weidner et al. 2009). For stars above 1 M_⊙, the majority of companion stars are likely to be significantly less massive and less luminous than the primary (see Kobulnicky et al. 2014), resulting in minimal change to the CMDs. Recent work on binary mass transfer suggests that the precise details of interactions between binary (and multiple) companions may influence IMF determinations (e.g., de Mink et al. 2014; Schneider et al. 2015), but full stellar evolution models that include binary mass transfer are not yet available.

Within PHAT, we verified that our choice in binary fraction does not change the results. For example, in the case of AP 94 (Figures 2 and 3), we modeled the CMD with binary fractions of 0.35 and 0.70 and found differences in the resulting IMF slope constraints to be <0.02 dex. Tests including different binary fractions on artificial clusters yielded similarly small changes.

Zoom In Zoom Out Reset image size

However, we did find that extreme binary fraction values, i.e., near 0 or 1, resulted in significantly worse CMD fits. In these two cases, the model either entirely populated the single star sequence (binary fraction = 0) or only sparsely populated it (binary fraction = 1), which substantially reduced the resemblance of the model to the observed CMD. Thus, for the present work, we adopted a fixed binary fraction of 0.35, which appears to be a reasonable value from LMC cluster studies (e.g., Elson et al. 1998).

To model background populations, we construct an empirical CMD from an annulus around the cluster that is ten times its size. The large area is necessary to accurately constrain the background scaling parameter.

Given a set of cluster model parameters (age, extinction, metallicity, and MF slope) the resulting model CMD, and the empirically derived CMD for fore-/background sources are compared to the observed CMD in a likelihood sense (Dolphin 2002). The resulting posterior probabilities, $p(t,{A}_{V},{\rm \Gamma })$ are illustrated in Figure 3 for a typical cluster in the sample. The prior on metallicity tightly brackets solar metallicity and the cluster mass emerges as a normalization factor, rather than a fit parameter. The smoothness of the $p(t,{A}_{V},{\rm \Gamma })$ of each cluster allows us to interpolate in order to increase the resolution of our PDFs.

Zoom In Zoom Out Reset image size

Finally, we emphasize that our analysis does not account for all sources of systematic uncertainty. Perhaps the more important systematic is the assumed stellar evolution models. Different stellar models are known to produce CMDs that can look quite different (e.g., Gallart et al. 2005), particularly for evolved phases of evolution. In principle, measuring each cluster's MF should be done with multiple stellar models, to assess the robustness of the MF to stellar evolution effects. In practice, this is currently not possible. Aside from the Padova models, at this time, there are no other stellar libraries that include the range of ages and phases of stellar evolution necessary to analyze this diverse cluster population. If such models become available in the future, we will revisit our analysis to assess the impact of systematics. In addition to stellar models, other sources of systematics include the binary star model (discussed above), the (unlikely) possibility of significant differential extinction internal to each cluster, the assumption that the IMF is a power-law over this mass range (c.f. Scalo 1986 for discussion of various possible forms for the high-mass IMF), and complex dynamical effects which may result in the preferential ejection of stars of a particular mass (e.g., Oh et al. 2015a).

One challenge for past MF studies was how to combine measurements from individual clusters into a broader statement about the MF slope of "clusters in general." Typically, the variance weighted mean is used to compute the "average" MF from an ensemble (e.g., Scalo 1986). However, this can potentially be a biased estimator (e.g., assumptions of normally distributed uncertainties), and does not necessarily use all available information (e.g., the probability that an object is actually a cluster).

Here, we adopt a framework that allows us to infer population wide characteristics about the distribution of MF slopes, given a set of noisy measurements. This model follows the approach and implementation laid out in Hogg et al. (2010) and Foreman-Mackey et al. (2014). We briefly describe the model below and we refer the reader to those papers for more detail.

As a simple model for the distribution of cluster MFs, we assume that it can be described by a normal distribution, $\mathcal{N}({\rm \Gamma },{\sigma }_{{\rm \Gamma }})$, with a mean of Γ and an intrinsic dispersion of ${\sigma }_{{\rm \Gamma }}$. We initially assume that their parameters are independent of age, mass, and size of the cluster. If we assume the probability that each object is a cluster, ${\boldsymbol{Q}}_{k}$, we can then write down a Gaussian mixture model likelihood function for the MF of clusters, where we simply assume that the MFs of falsely presumed clusters have a different $\mathcal{N}({{\rm \Gamma }}_{\mathrm{false}},{\sigma }_{{\rm \Gamma },\mathrm{false}})$. This mixture model mitigates the impact of erroneously identified clusters, without forcing a binary decision on which clusters have been identified with sufficient confidence.

More explicitly, the mixture model can be written as

$$\begin{eqnarray}&&\begin{array}{c}p({\boldsymbol{\Gamma }}_{k}\;| \;\boldsymbol{\theta },{\boldsymbol{Q}}_{k})={\boldsymbol{Q}}_{k}\;\mathrm{exp}\left(-({\rm \Gamma }-{\boldsymbol{\Gamma }}_{k})/(2{\sigma }_{{\rm \Gamma }}^{2})\right)\\ \quad +(1-{\boldsymbol{Q}}_{k})\;\mathrm{exp}\left(-({{\rm \Gamma }}_{\mathrm{false}}-{\boldsymbol{\Gamma }}_{k})/(2{\sigma }_{\mathrm{false}}^{2})\right)\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{\Gamma }}_{k}$ is the marginalized PDF for the MF of the kth object and ${\boldsymbol{Q}}_{k}$ is the probability that the kth object is a cluster. Γ and ${\sigma }_{{\rm \Gamma }}^{2}$ are the mean and variance of the Gaussian of interest for the clusters and ${{\rm \Gamma }}_{\mathrm{false}}$ and ${\sigma }_{\mathrm{false}}^{2}$ are nuisance parameters for possible contaminating sources. Simply put, with this model, the more probable an object is a cluster the more it contributes to the MF parameters of interest. In this model context, a universal Kroupa MF would be represented by ${\rm \Gamma }=1.30$ and ${\sigma }_{{\rm \Gamma }}^{2}\to 0$, and universal Salpeter MF by ${\rm \Gamma }=1.35$ and ${\sigma }_{{\rm \Gamma }}^{2}\to 0$.

This same framework also readily allows us to explore whether the mean MF slope depends on any of the other cluster properties such as age, mass, and size. Specifically, we can generalize the model to

$$\begin{eqnarray}&&{\rm \Gamma }({t}_{k},{M}_{k},{r}_{k})=\bar{{\rm \Gamma }}+{a}_{m}\;\mathrm{log}{\displaystyle \frac{{M}_{k}}{{M}_{{\rm c}}}}\;+{a}_{t}\;\mathrm{log}{\displaystyle \frac{{t}_{k}}{{t}_{{\rm c}}}}+{a}_{r}\;\mathrm{log}{\displaystyle \frac{{r}_{k}}{{r}_{{\rm c}}}},\end{eqnarray} \tag{ 2 }$$

where $\{{t}_{k},{M}_{k},{r}_{k}\}$ are the most likely age, mass, and size for the kth object and $\{{t}_{{\rm c}},{M}_{{\rm c}},{r}_{{\rm c}}\}$ are the mean of our cluster sample in each case. We assume priors that are flat in $\bar{{\rm \Gamma }},\;\mathrm{log}{\sigma }_{{\rm \Gamma }},\;{a}_{m},\;{a}_{t},\;{a}_{r}$ over the ranges indicated in Table 1, and sample these parameters' PDFs with the affine invariant ensemble Markov chain Monte Carlo program emcee (Foreman-Mackey et al. 2013).¹⁹

It is important to place our results into context with other high-mass IMF studies in the LG. Among the strongest evidence for the universality of the high-mass IMF in the LG comes from a series of papers that combine spectroscopy and photometry of resolved, massive stars in clusters and OB associations in the LMC and MW (Massey et al. 1995a, 1995b; Massey 1998, 2003).

To summarize, these studies found (a) no drastic differences between the IMF slopes of very young, <5 Myr star-forming regions in the MW and LMC, despite large differences in metallicity; and (b) that the resulting IMF slopes are broadly consistent with Salpeter/Kroupa. Combined, these studies represented a major milestone in our understanding of the high-mass IMF and its (in-)sensitivity to environment.

Although re-analyzing these studies is beyond the scope of this paper, we can quantify what the ensemble high-mass IMF slopes are in the LMC and MW using data from Massey (2003). The purpose is to give a general sense of how different our M31 result is from other LG environments. To perform this analysis, we extracted the necessary data from Massey (2003) and have listed it in Table 3. We simply use the model from Equation (1) with ${\boldsymbol{Q}}_{k}=1$ for each object, i.e., we know they are all bone fide clusters, to infer the mean and variance for the MW and LMC clusters. We do not model possible age, mass, size dependencies.

Table 3.  Select Literature High-mass IMF Slopes

Cluster	Γ	Mupper
		(M${}_{\odot }$)
(1)	(2)	(3)
	LMC
LH 6	1.7 ± 0.4	85
LH 9	1.7 ± 0.4	55
LH 10	1.7 ± 0.4	90
LH 38	1.7 ± 0.4	85
LH 47/48	1.3 ± 0.2	50
LH 58	1.4 ± 0.2	50
LH 73	1.3 ± 0.4	65
LH 83	1.3 ± 0.5	50
LH 114	1.0 ± 0.1	50
LH 117/118	1.6 ± 0.2	>120
R136	1.3 ± 0.1	>120
	MW
NGC 6823	1.3 ± 0.4	40
NGC 6871	0.9 ± 0.4	40
NGC 6913	1.1 ± 0.6	40
Berkeley 86	1.7 ± 0.4	40
NGC 7280	1.7 ± 0.3	65
Cep OB5	2.1 ± 0.6	30
IC 1805	1.3 ± 0.2	100
NGC 1893	1.6 ± 0.3	65
NGC 2244	0.8 ± 0.3	70
NGC 6611	0.7 ± 0.2	75
Cyg OB2	0.9 ± 0.2	110
Tr 14/16	1.0 ± 0.2	>120
h & χ Per	1.3 ± 0.2	120

Note. High-mass IMF slopes from the massive star studies in the MW and LMC clusters as listed in Massey (1998) and Massey (2003). The dynamic mass range of the clusters extend from >1 M_⊙ to the mass listed in column (3).

Download table as:  ASCIITypeset image

From this exercise we find mean values of ${{\rm \Gamma }}_{\mathrm{MW}}={1.16}_{-0.10}^{+0.12}$ and ${{\rm \Gamma }}_{\mathrm{LMC}}={1.29}_{-0.11}^{+0.11}$. In both cases, ${\sigma }_{{\rm \Gamma }}\sim 0.3-0.4$. The large scatter in the posterior distributions is likely a reflection of underestimated uncertainties and/or systematics, as opposed to true physical variation.

Massey (2003) emphasizes that the reported uncertainties are only statistical in nature. However, even these are likely underestimated, as discussed in Weisz et al. (2013). Unfortunately, not all of these studies in the literature list both the number of stars and mass ranges needed to estimate degree of underestimation of the random uncertainties. There are also issues of systematic uncertainties in stellar models that can significantly change the masses of the most massive stars used in these studies (e.g., Massey 2011).

However, even from this cursory exercise, we see that the MW, LMC, and M31 do not appear to have drastically different IMF slopes. While this is quite remarkable given the diversity of environments, it should be tempered by the heterogenous nature of this comparison. If high-mass IMF variations are subtle, then homogenous and principled analyses are absolutely necessary to uncover them. At present, only our analysis in M31 is the result of a systematic study of the high-mass IMF over a large number of young clusters. Any comparison of the high-mass IMF slopes of LG galaxies must come from similarly homogenous data and principled analysis in order to minimize the systematics that have been persistent in IMF studies for the better part of six decades.

Given the analysis presented in this paper and comparison with LG literature, we do not find strong evidence that a universal Salpeter/Kroupa IMF is the best representation of the high-mass IMF slope in the LG. Although we only rule it out at a $\sim 2\sigma$ level, there appears to be little basis for this canonical value being correct in the first place.²⁰ Instead, for practical purposes the high-mass IMF can be represented by a single-sloped power-law with ${1.45}_{-0.06}^{+0.03}$, where the reported uncertainties are statistical only (see Section 3.1). Although this is only 0.1–0.15 dex steeper than Salpeter/Kroupa it does make a difference in the number of high-mass stars formed, and related quantities such as commonly used star formation rate (SFR) indicators. We discuss these implications further in Section 5.3.

For general usage, we recommend the adoption of the following forms of the the IMF over all stellar masses. With the M31 high-mass IMF, a modified Chabrier IMF Chabrier (2003) has the form

$$\begin{eqnarray}\xi (m){\rm \Delta }(m) & = & \;0.158\;(1/m)\;{\mathrm{exp}}^{\frac{-{(\mathrm{log}(m)-\mathrm{log}(0.08))}^{2}}{2\times {0.69}^{2}}}\\ \xi (m)=c\;{m}^{-({\rm \Gamma }+1)},\ {\rm \Gamma } & = & {1.45}_{-0.06}^{+0.03},{\rm for}\;{\rm }1.0\lt m\lt 100\ \;{M}_{\odot },\end{eqnarray} \tag{ 3 }$$

and a modified Kroupa (2001) IMF has the form

$$\begin{eqnarray}\xi (m)={\text{}}{{cm}}^{-({\rm \Gamma }+1)}\left\{\begin{array}{cc}{\rm \Gamma }=0.3 & {\rm for}\;{\rm }0.08\lt m\lt 0.5\;{M}_{\odot }\\ {\rm \Gamma }=1.3 & {\rm for}\;{\rm }0.5\lt m\lt 1.0\;{M}_{\odot }\\ {\rm \Gamma }={1.45}_{-0.06}^{+0.03}, & {\rm for}\;{\rm }1.0\lt m\lt 100\;{M}_{\odot }.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

Here we have made two extrapolations to match conventional IMF definitions in the literature. First, we have extrapolated the M31 IMF down to 1 M_⊙. For most practical purposes, this extension from 2 M_⊙ to 1 M_⊙ affects the resulting stellar mass spectrum at the level of a few percent.

Second, we have extrapolated the M31 IMF up to 100 M_⊙. While our data do not provide good constraints above ∼25 M_⊙ (see Section 1.2), there is little other information to go on for mass spectra of the most massive stars. Thus, as is commonplace in the literature, we suggest extrapolating the M31 IMF up to the highest masses, but recognize that systematic spectroscopic searches for the massive stars are absolutely critical for understanding the IMF of the most massive stars (e.g., Massey et al. 1995b; Massey 2003).

Finally, we have assumed that the sub-solar IMF in M31 is similar to that of the Galaxy. Although not confirmed by direct star counts, Conroy & van Dokkum (2012) have shown that spectral features sensitive to the low-mass IMF yield a result that is similar to the low-mass Galactic IMF.

In Figure 7 we illustrate the mass spectra of high-mass stars predicted by various commonly used IMF models. The plot shows the number of stars expected relative to a Kroupa IMF (${\rm \Gamma }=1.3$) as a function of mass. For fair comparison, the slopes have all been normalized at 2 M_⊙, the lower main sequence mass limit of the PHAT data, and extended to 100 M_⊙ on the upper end. We chose this normalization, as opposed to 1 M_⊙ to discuss our ability to distinguish different forms of the IMF in the literature over the dynamic mass range of the PHAT data, which is indicate by the vertical dashed lines. As previously discussed, the indicated upper limit (25 M_⊙) is likely a lower bound due to the inability of optical-only photometry to accurately characterize the most massive stars.

Zoom In Zoom Out Reset image size

The M31 high-mass IMF predicts the formation of fewer stars than Kroupa at all masses. The red solid line presents our median IMF model and shows that the fractional deficit varies from ∼0.9 at ∼10 M_⊙ to ∼0.7 at 100 M_⊙. The red shaded region reflects the 68% confidence interval in the number of predicted high-mass stars due to uncertainty on the mean high-mass IMF slope in M31. At most masses, the uncertainty deviates from the median by ± ∼5%, but increases to ±15% for stars >50 M_⊙.

The M31 high-mass IMF is similar to other high-mass IMF models in the literature. It is particularly close to the Kennicutt IMF (${\rm \Gamma }=1.5$; Kennicutt 1983). It is also the same as the IMF of Miller & Scalo (1979) over the 1–10 M_⊙ range. The steepness of the Miller & Scalo (1979) IMF above 10 M_⊙ is known to be too extreme (e.g., Kennicutt 1983; Kennicutt et al. 1994). Further, the high-mass IMF from Scalo (1986) has a similar slope (${\rm \Gamma }=1.6$) above ∼4 M_⊙, but a much steeper slope below it. Finally, the M31 high-mass IMF predicts fewer high-mass stars than a Salpter IMF for masses >2 M_⊙.

We further quantify differences in these IMF models in Table 4. Here we have computed the expected number of stars per 10⁶ M_⊙ formed at selected stellar masses assuming a Kroupa (2001) IMF below 2 M_⊙ and the listed IMFs above 2 M_⊙. This table solidifies many of the above points regarding differences in commonly used high-mass IMF models. It also highlights an important implication for core collapse supernovae, whose progenitors are believed to have masses ≳8 M_⊙(e.g., Smartt 2009). For example, the median M31 high-mass IMF model predicts ∼25% fewer stars with masses ≳8 M_⊙ compared to a Kroupa IMF.

Table 4.  Number and Fraction of Massive Stars for Different IMFs

IMF	${N}_{\star }\geqslant 2{M}_{\odot }$	${N}_{\star }\geqslant 8{M}_{\odot }$	${N}_{\star }\geqslant 20{M}_{\odot }$	${N}_{\star }\geqslant 50{M}_{\odot }$	${N}_{\star }\geqslant 90{M}_{\odot }$	Frac. $\geqslant 2$ ${M}_{\odot }$	Frac. $\geqslant 8$ ${M}_{\odot }$	Frac. $\geqslant 20$ ${M}_{\odot }$	Frac. $\geqslant 50$ ${M}_{\odot }$	Frac. $\geqslant 90$ ${M}_{\odot }$
	${\displaystyle \frac{{10}^{4}\mathrm{stars}}{{10}^{6}\;{M}_{\odot }}}$	${\displaystyle \frac{{10}^{4}\mathrm{stars}}{{10}^{6}\;{M}_{\odot }}}$	${\displaystyle \frac{{10}^{3}\mathrm{stars}}{{10}^{6}\;{M}_{\odot }}}$	${\displaystyle \frac{{10}^{2}\mathrm{stars}}{{10}^{6}\;{M}_{\odot }}}$	${\displaystyle \frac{{10}^{1}\mathrm{stars}}{{10}^{6}\;{M}_{\odot }}}$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
Kroupa	6.9	1.1	3.0	6.2	6.3	1	1	1	1	1
Salpeter	6.8	1.0	2.7	5.4	5.3	0.97	0.90	0.87	0.84	0.83
Weisz et al.	${6.7}_{-0.0}^{+0.1}$	${0.88}_{-0.04}^{+0.08}$	${2.2}_{-0.2}^{+0.3}$	${4.0}_{-0.4}^{+0.8}$	${3.8}_{-0.4}^{+0.9}$	${0.90}_{-0.01}^{+0.04}$	${0.74}_{-0.04}^{+0.09}$	${0.66}_{-0.06}^{+0.11}$	${0.60}_{-0.07}^{+0.13}$	${0.56}_{-0.06}^{+0.14}$
Kennicutt	6.7	0.82	1.9	3.5	3.2	0.87	0.67	0.57	0.50	0.46
Miller & Scalo	6.8	0.66	0.82	0.81	0.56	0.85	0.51	0.23	0.11	0.08
Scalo	5.6	0.20	0.43	0.72	0.65	0.54	0.14	0.10	0.08	0.07

Note. The number of expected stars above the specified mass limit per 10⁶ M_⊙. For these calculations, we have adopted a Kroupa (2001) IMF below 2 M_⊙, the specified IMF above 2 M_⊙, and assumed lower and upper mass bounds for the total IMF to be 0.08 M_⊙ and 100 M_⊙. Columns (7)–(11) give the fraction of stars formed above the listed mass, relative to the standard Kroupa IMF.

Download table as:  ASCIITypeset image

Finally, we consider the implications for commonly used SFR indicators (e.g., Kennicutt 1998b; Kennicutt & Evans 2012). Using the Flexible Stellar Population Synthesis code (Conroy et al. 2009, 2010), we compute the luminosity to SFR conversion coefficients for GALEX far- and near-UV luminosities, along with Hα, which we list in Table 5. To do this, we have assumed a constant SFH over the last 1 Gyr, the indicated IMF, mass limits from 0.08 to 100 M_⊙, and case B recombination to convert the number of ionizing photons to Hα luminosity. From Kennicutt & Evans (2012) the canonical conversion from luminosity to SFR can be written as

$$\begin{eqnarray}&&\mathrm{log}\;{\stackrel{\dot{}}{M}}_{\star }\;({M}_{\odot }\;{\mathrm{yr}}^{-1})=\mathrm{log}\;{L}_{x}\;-\;\mathrm{log}\;{C}_{x}\end{eqnarray} \tag{ 5 }$$

where ${\stackrel{\dot{}}{M}}_{\star }$ is the SFR, L_x is the luminosity over the bandpass x, and C_x is the conversion constant over the same wavelength range. We have provided updated values for C_x in Table 5.

Table 5.  Updated Broadband SFR Indicators

Band	Lx Units	Log(Cx)	Log(Cx)	${\stackrel{\dot{}}{M}}_{\star }(M31)/{\stackrel{\dot{}}{M}}_{\star }$(Kroupa)
		(Kroupa)	(M31)
(1)	(2)	(3)	(4)	(5)
FUV	erg s−1 ($\nu {L}_{\nu }$)	43.29	${43.18}_{-0.02}^{+0.05}$	${1.28}_{-0.12}^{+0.06}$
NUV	erg s−1 ($\nu {L}_{\nu }$)	43.14	${43.04}_{-0.02}^{+0.03}$	${1.24}_{-0.11}^{+0.06}$
Hα	erg s−1	40.91	${40.73}_{-0.07}^{+0.04}$	${1.52}_{-0.24}^{+0.13}$

Note. A comparison of commonly used SFR indicators for a Kroupa and M31 IMF values. The values of L_x, C_x, and ${\stackrel{\dot{}}{M}}_{\star }$ follow Kennicutt & Evans (2012) and Equation (5).

Download table as:  ASCIITypeset image

From column (5) of Table 5, we see modest differences in the SFR indicators due to an M31 high-mass IMF. The variations are such that for a fixed luminosity, using the M31 high-mass IMF results in SFRs that are a factor ∼1.3–1.5 larger.

Although the change in absolute calibration is modest, the IMF measurement in M31 provides a significant increase in the precision of the SFR. Including uncertainties the high-mass Kroupa IMF slope is ${\rm \Gamma }=1.3\pm 0.7$. Propagating this uncertainty through Equation (5) results in an order of magnitude uncertainty in Hα SFRs. In contrast, the precision of our IMF measurement in M31 reduces this uncertainty to be of the order unity, representing a substantial improvement in our knowledge of galaxy SFRs (under the assumption of a universal IMF).