THE SUPERNOVA PROGENITOR MASS DISTRIBUTIONS OF M31 AND M33: FURTHER EVIDENCE FOR AN UPPER MASS LIMIT

First, we select SNRs from catalogs which overlap with HST fields in M31 and M33. We limited our search to fields observed with the ACS and WFPC-2 instruments, as the optical CMDs tend to be deeper and offer the best SFH constraints. For M31, we combine three SNR catalogs to select targets from: Braun & Walterbos (1993); Magnier et al. (1995); Williams et al. (1995). All three catalogs make use of [S ii]-to-Hα ratios to identify SNR. Our M33 SNR catalog comes from Long et al. (2010), who selected SNRs using both [S ii]-to-Hα ratios and X-ray observations. Long et al. (2010) incorporated all previous SNR catalogs from the literature in their analysis. As mentioned in Section 1, all new M31 data analyzed in this work come from newly available PHAT data sets and therefore have uniform filters and exposure times, while the M33 data comes from various archival data sets coincident with SNR. The locations of SNRs analyzed in this paper and J12 are plotted in Figure 1. The full list of SNRs analyzed, including the associated HST data sets used, are in Table 1.

Zoom In Zoom Out Reset image size

We used the photometry pipeline from the ACS Nearby Galaxy Treasury Program (Dalcanton et al. 2009) and the PHAT Program (Dalcanton et al. 2012) to perform resolved stellar photometry on all stars in the selected fields. The full details of how this photometry is performed are provided in Dalcanton et al. (2009, 2012). Briefly, the pipeline uses DOLPHOT (Dolphin 2000) to fit the well-characterized ACS point-spread function to all of the point sources in the image. Photometric measurements are then converted to Vega magnitudes using zero-points from the ACS handbook. We use fake star tests to assess photometric uncertainties and completeness; 10⁵ fake stars of known color and magnitude are inserted into the full fields and blindly recovered using the identical software.

We assume distance moduli of 24.47 to M31 (McConnachie et al. 2005) and of 24.69 to M33 (Barker et al. 2011). Typical uncertainties of from these surveys are ∼0.05 mag for M31 and ∼0.1 mag for M33. We fix the distance moduli to these values for the remainder of the paper.

We found in J12 that it was difficult to constrain older SF events in shallow fields. We defined the "depth" of the field as the point at which photometric completeness dropped below 50% and adopted the same depth cut as J12 for this work. The values used are F475W = 24.5 for M31 data, shifted to F475W = 24.7 for the M33 fields. In practice, nearly all the fields examined in this work have depths significantly below these depth limits, making the precise placement of the cut largely unimportant.

After performing full-field photometry, we select all stars in an annulus of ∼50 pc around the location of the SNR of interest, which we use to then measure the SFH of the region. We use the software package MATCH (Dolphin 2002) to fit the observed CMD with synthetic ones based on the models of Marigo et al. (2008); Girardi et al. (2010). For purposes of populating the models, we assume a Salpeter IMF (dN/dM∝M^−2.35) and a binary fraction of 0.35. J12 and Gogarten et al. (2009a) found no appreciable effect on measured SFH values from reasonable variations of these parameters. MATCH also makes use of fake stars to evaluate completeness and photometric uncertainties, which we extract from an annulus ∼2.5 times the size to ensure a sufficient number of fake stars. The distance modulus is fixed to the respective value for each galaxy from above. The CMD is binned with bin sizes 0.3 in magnitude and 0.15 in color. Metallicity is constrained to have a spread of ∼0.15 dex, and to increase or stay constant with time.

Our age bins for SFH determination are in log-space from log (Age) = 6.60 (4 Myr) to log (Age) = 10.10 (12.5 Gyr) in steps of 0.05 dex. For a solar metallicity system, these ages would correspond to supernova progenitors of roughly 52 M_☉ and 7.3 M_☉ respectively. Note that any SF MATCH finds for ages younger than 4 Myr is included in the 4 Myr to 4.4 Myr bin. Thus SF found in this youngest bin actually acts as an upper limit on the age, constraining the SF to be younger than 4.4 Myr. The final result from our CMD fitting analysis is an SFH from the present back to 12.5 Gyr.

Our treatment of reddening in our CMD fitting was the same as in J12, and merits further discussion. Our prescription for reddening is to assume that A_V follows a top-hat distribution, the width of which is specified by a dA_V parameter defined by the user. MATCH then fits for the minimum value of this A_V distribution, defining the top-hat in the course of the fitting procedure. To determine the value of the parameter most appropriate for each SNR location, we increased the width of the dA_V parameter until the fit returned a reddening value consistent with the Schlegel et al. (1998) foreground reddening value. In other words, we treat all reddening in addition to foreground MW extinction as broadening of the CMD due to differential reddening.

For purposes of age/mass conversion, we use the same models as those fitted to the CMDs. For each age bin, we assume the star with the largest zero-age main sequence mass (M_ZAMS) remaining in the isochrone will be the next one to become a SN. This mass is then taken as the progenitor mass for that given age bin. We perform a simple linear interpolation across age bins to get an estimated progenitor mass for an arbitrary age. We adopt the Z = 0.019 isochrones for M31 and the Z = 0.008 isochrones for M33, which are consistent with gas-phase metallicity measurements for the two galaxies (Blair et al. 1982; Barker et al. 2011). However, as demonstrated in J12, the choice of metallicity is largely unimportant in converting age to mass as the M_ZAMS masses are very similar across reasonable ranges of metallicity. For a given age, the difference in mass between the two isochrones is typically a few percent.

Numerous sources of theoretical and observational evidence point toward a minimum mass necessary for a star to undergo core collapse. We wish to identify the age range where a population may still produce an SN. Coeval populations that are too old will have no stars remaining massive enough to become core-collapse SNe. In J12, we found our data most consistent with a minimum mass of 7.3 M_☉, corresponding to an age of 50 Myr. We verified that our new data are also consistent with this minimum mass, and therefore only consider ages ⩽50 Myr in our age distributions. However, we were not able to improve this constraint using our current data set.

We do not assume any maximum mass at this stage. However, for ages younger than 4.4 Myr, we may only quote a limit on the progenitor mass of >52 M_☉ for the Z = 0.019 isochrone (i.e., M31) and >55 M_☉ for the Z = 0.008 isochrone (i.e., M33) because the optical CMDs are degenerate for ages younger than 4.4 Myr.

To estimate the age of the progenitor, we identify the age at which, over the past 50 Myr, 50% of the stars have formed. This median age is then converted to a median mass using the interpolation defined above. We also assign uncertainties to the age from two sources using the methods discussed in the next section. Furthermore, we calculate full probability distributions for each progenitor. We assume that the progenitor originates from any SF that has occurred in the region over the past 50 Myr, where the fraction of stellar mass formed in a given age bin corresponds to the probability that the progenitor star was of that age. Thus, all of our probability distributions sum to one because they include all SF that has occurred over the past 50 Myr.

Uncertainties on the age estimate may arise either from uncertainties in the CMD fitting process, or from the age spread present in the stellar population.

To characterize uncertainties from fitting, we use a Monte Carlo (MC) technique in which we resample the CMD to account for Poisson noise using the MATCH software (see Dolphin 2002). We also apply random shifts to the isochrones in temperature, with σ = 0.02, and bolometric luminosity, with σ = 0.17, during this procedure. These random shifts are intended to mimic potential uncertainties in the stellar evolution models themselves. The magnitude of the luminosity shift is also larger than the uncertainty in the distance, thereby incorporating those uncertainties. Each resampling is then refit using our identical CMD fitting procedure, but with fixed reddening. These resulting SFHs provide uncertainties on the SF in a given bin. With these uncertainties, we calculate a probability distribution for the value of the median age. We adopt the 16% and 84% range of this distribution as our uncertainty in the median age.

Another source of uncertainty is the intrinsic spread of the SFH across multiple age bins. There are many SNR populations in which there are two (or more) distinct SF events, and certainly some where a single SF event may spread across two or more adjacent bins. To account for this, we find the locations where 16% and 84% of the stellar mass was formed and adopt these as our uncertainties due to the spread of SF. These age-spread uncertainties therefore include an age range that contains 68% of the young stellar mass present in the region.

We add these two sources of uncertainty in quadrature to find our final (potentially asymmetric) estimates of uncertainty on the median age. Finally, we argue that we cannot truly constrain our progenitors to any better precision that that afforded by the age bins in which we have performed our CMD fitting. In other words, we always round our uncertainties to the nearest age bin boundary.

In practice, the second source of uncertainty due to the spread of SF across multiple age bins tends to dominate that provided by the random uncertainties in the fitting process. This is because, while the fitting process may affect the amount of SF found in a given bin, the change in SF must be very large to modify the estimate of the median age. The random uncertainty must be enough to significantly change the relative prominence of a preferred burst of SF. For most of the MC realizations, it is not, and as a result, the spread of SF tends to be the largest factor.

The final source of uncertainty, which we neglect in determining our final answers, is the conversion from age to mass, which depends on the models applied. We do not include this uncertainty in our mass tables, as our uncertainties assume the Marigo et al. (2008) and Girardi et al. (2010) models to be consistent with our SFH fitting as they are used in both processes. In J12, we estimated that the effects of applying different models for the age/mass conversion would potentially be ∼0.5–1.0 M_☉ (younger ages have larger uncertainty) by comparing the predicted masses of Pietrinferni et al. (2004) with those of Marigo et al. (2008) and Girardi et al. (2010). However, as we have used the Marigo et al. (2008) and Girardi et al. (2010) models to derive these ages in the first place, it is not clear that comparing the Pietrinferni et al. (2004) masses to those of our chosen isochrones is reflective of the actual data.

We make multiple assumptions in our technique based upon justifications detailed in J12. In this section, we briefly review these assumptions and their potential impacts on our conclusions.

First, naturally, in using the surrounding stellar population to infer information about the progenitor star, we are assuming that the two are evolutionarily linked. Given that the vast majority of stars form in clusters (Lada & Lada 2003) with coeval populations, and that that these populations remain spatially associated on 50 pc scales for 50 Myr, we expect this to be a reasonable assumption (see also Bastian & Goodwin 2006; Gogarten et al. 2009b; Eldridge et al. 2011).

Related to this point is that our methodology is also contingent on the accuracy, completeness, and selection effects of the SNR catalogs used. If the SNR catalogs are biased toward one specific type of progenitor or environment, this could bias the overall inferred distribution. Braun & Walterbos (1993) and Magnier et al. (1995) both consider extinction to be a minor problem. Both studies note that detection biases may exist as a function of the age of the remnant, but do not note any biases toward any particular age of the surrounding stellar population. Related to this, we have also assumed that SNRs behave similarly in different environments. If SNRs in differently aged populations had significantly different lifetimes or luminosities, we would naturally be subject to selection effects from this.

We assume for our study that the SNR catalogs used are not biased toward any particular progenitor age, and we have included as many SNR locations as possible in our analysis. We note the consistency of the M33 and M31 results and the consistency between different subsets of the M31 sample, which depend on different SNR catalogs. This suggests that there are no significant selection biases imposed by the detection method or survey location. However, we have no real way to test the second caveat, that SNRs from differently aged populations have different lifetimes or detectability thresholds. We discuss this more in Section 4.3.

Next, our method provides no information about the type of SN which created the SNR. For CCSNe, this ambiguity does not affect the mass estimation process, since all CCSNe are linked to the deaths of massive stars and therefore young stellar populations (although the mass distribution of progenitors as a function of SN type is certainly an interesting question, we are unable to investigate this with our methodology). We note that there could be contaminants from thermonuclear Type Ia SNe. However, Type Ia SNe will generally be associated with older stellar populations. The absence of any recent SF from the CMD analysis will allow one to remove some fraction of these contaminants. We identified six such SNRs with zero recent SF in J12. In this work, we only find one such SNR in our new M31 sample (K413) and zero such SNRs in our M33 sample. The new M31 data analyzed in this work essentially all comes from the star-forming ring of M31 due to the distribution of the PHAT survey, so it is not entirely surprising that nearly all fields analyzed return significant recent SF and that nearly all of the SNe are CCSNe. The result of all analyzed M33 SNR populations having young SF is also likely due to the fact that M33 has a higher SF intensity than M31.

Typical fractions of Type Ia SNe compared to all SNe are ∼25%, although such a figure is dependent on the galaxy type in question (Li et al. 2011). We would expect a higher fraction for M31 and a lower fraction for M33 based on morphological distinction. In addition, there are many reasons to assume that this Ia fraction may actually be an upper limit, with a bias resulting from the faintness of type CCSNe compared to Type Ia SNe (see discussions in Thompson et al. 2009; Horiuchi et al. 2009; Horiuchi et al. 2011). Finally, both Braun & Walterbos (1993) and Magnier et al. (1995) explicitly note that they will be significantly more biased toward identifying CCSNe remnants over Type Ia remnants due to the spatial distribution of their surveys. We do not include the M31 bulge in our sample, where most of the Ia events would be expected to occur. Thus we actually expect Ia SNe in our sample to be a significantly smaller fraction than would be expected based on galaxy-wide SN rate studies.

Furthermore, in J12, we demonstrated that our results are robust to inclusion of additional misidentified Type Ia SNe by removing additional potential contaminants even up to 25%. We identified those SNR locations with minimal GALEX FUV flux and removed those objects from our sample. The resulting progenitor mass distribution was essentially unchanged (the 95% confidence interval on the power-law exponent changed from 2.7 < α < 4.4 to 2.6 < α < 4.3). We conclude that, while there is uncertainty associated with an individual progenitor mass measurement due to the possibility of it being a Type Ia SNe, the overall distributions of progenitor masses are largely robust to the inclusion of a few progenitors with random ages since the vast majority of our SNR will be the result of CCSNe. In other words, while there are almost certainly some Type Ia contaminants masquerading as CCSNe in our sample, they will be a small number and will not significantly modify the overall conclusions.

Finally, our method is contingent on the accuracy of the stellar evolution models used to populate our model CMDs. Our CMD-derived SFHs include estimates of these model uncertainties (see Section 2.4), although we neglect them during our mass conversion procedure. It is worth noting that studies of directly identified progenitor stars suffer from similar uncertainties in that they must fit the spectral energy distributions to derive a luminosity and temperature and use stellar evolution models to infer the mass of late-stage massive stars (e.g., Smartt et al. 2009). These are thought to be the most uncertain stages of the model predictions. Our approach makes use of the entire CMD, making our results less sensitive to the details of the late-stage stellar evolution models.

To properly include uncertainties in our mass distribution fitting, we adopted a probabilistic framework developed by Weisz et al. (2013) to fit a power-law function to our mass distributions. A full probability distribution for the entire sample is created from the individual mass measurements and uncertainties. We then fit a power-law distribution using the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013) to sample the posterior distributions. We then extracted the functional parameters and robust uncertainties from these posterior distributions. Quoted uncertainties correspond to the 16th and 84th percentiles of the probability distributions. When performing our fitting, we set a prior that M_max  > 20 M_☉, and fixed the minimum mass to 7.3 M_☉. As an illustration of the distribution of power-law indices, we plot 100 randomly selected MCMC steps in Figure 4.

We found that we were unable to constrain the maximum mass beyond our prior if both the mass function slope and the maximum mass were left as free parameters (all values >20 M_☉ were found to be consistent with the data if α  ∼  4). Therefore, we used two different approaches to fit our distributions. For one, we assumed the distribution followed a power-law function of the form dN/dM∝M^−α, with α left as a free parameter and with M_max fixed to 90 M_☉. We also performed an MCMC run in which we fixed the functional form of the mass function to be a Salpeter IMF, dN/dM∝M^−2.35, as one may expect to recover if all massive stars produce SNe. In these fits, M_max was constrained by the data. Weisz et al. (2013) explicitly notes that a meaningful limit on the maximum mass cannot be constrained by the data beyond simply returning the most massive star observed, so we emphasize that fit values for M_max in particular should be considered illustrative. The precise values are not necessarily meaningful, whereas the fits for a free slope are more reliable.

The results from these fits are listed in Table 4. We do not attempt to evaluate which model description is a better fit to the data. There is evidence that the regions of progenitor masses over which SNe may explode may actually be quite complex (e.g., Sukhbold & Woosley 2014), making either parameterization of the mass distribution a likely oversimplification. However, the result of a paucity of massive progenitors is robust regardless of model selection. The use of the MCMC technique properly incorporates our measured uncertainties and rules out a Salpeter mass function for our progenitor sample.

Table 4. MCMC Inferred Mass Distribution Parameters

Sample	Best Fit αa	Maximum Mass
		(M☉)a
Index α as Free Parameter, Mmax = 90 M☉
M31 Sample (J12 + This Work)	$\alpha =4.4^{+0.4}_{-0.4}$	Mmax = 90
M33 Sample (This Work)	$\alpha =3.8^{+0.5}_{-0.4}$	Mmax = 90
Full Sample (J12 + This Work)	$\alpha =4.2^{+0.3}_{-0.3}$	Mmax = 90
Index Fixed to Salpeter (α = 2.35), Mmax as Free Parameter
Full M31 Sample (J12 + This Work)	2.35	$M_{\rm max} = 38^{+15}_{-6}$
Full M33 Sample (This Work)	2.35	$M_{\rm max} = 41^{+45}_{-12}$
Full Distribution (J12 + This Work)	2.35	$M_{\rm max} = 35^{+5}_{-4}$

Note. ^aQuoted uncertainties come from 16th and 84th percentiles of the probability distributions.

Download table as:  ASCIITypeset image

We find that, independent of the sample that we fit, all distributions analyzed favor a model in which some fraction of the massive star population do not explode as SNe, whether it is parameterized as a very steep mass function or a Salpeter mass function with a cutoff. Furthermore, the M31 and M33 fit values indicate that both the slope and the maximum mass cutoff are consistent to within the uncertainties, suggesting that any metallicity or star formation intensity effects do not change the index by more than ±  ∼1. While this range in power-law indices is still large, the distributions are all still steeper than Salpeter regardless. A more precise comparison of the M31 and M33 distributions would naturally be interesting, but we are unable to make a more constraining statement than this given the current size of the data set and the uncertainties involved.

Both the individual M31 and M33 progenitor mass distributions, as well as the full progenitor mass distribution, are steeper than a Salpeter IMF with at least 95% confidence, regardless of the method used to fit the distribution. Thus, both the full sample and the individual M31 and M33 distributions display a paucity of massive stars compared to a simple model IMF. If we assume that a Salpeter distribution is a reasonable description of the massive star population, then this paucity suggests that some fraction of massive stars are not in our SNR progenitor sample. Recent theoretical and observational work has suggested that in certain mass ranges, massive stars may not undergo SN explosions (e.g., Smartt et al. 2009; Horiuchi et al. 2011; Dessart et al. 2011; Sukhbold & Woosley 2014). Our observation is consistent with the interpretation that there are certain mass ranges in the massive star population in which stars fail to end their lives as CCSN.

One important caveat is that the interpretation of an observed upper-mass CCSNe cutoff relies on sampling an unbiased population of massive stars, i.e., that our SNR catalogs are unbiased to any one type of progenitor star environment. Since we would expect higher-mass progenitors to be preferentially found in more-extincted regions, this may cause SNR from higher-mass progenitors to be systematically under-sampled in SNR catalogs. We have now expanded the sample to include SNR from M33, which is less inclined to our line of sight and lower metallicity (Bresolin 2011; Zurita & Bresolin 2012). Both of these qualities should reduce the amount of extinction, somewhat reducing these biases. However, we also note that a commonly used method for identifying SNRs is by their high [S ii]/Hα flux ratios (e.g., Gordon et al. 1998). The youngest regions will have more photoionized H ii, increasing the Hα flux and making SNR identification more difficult. We made some attempt to reduce this bias by incorporating SNR catalogs that use different methodology. In particular, the Long et al. (2010) M33 SNR catalog incorporates observations from radio to X-ray to catalog SNRs. In any case, given that SNRs do not stand out if they are located within a photo-ionized interstellar medium, we are currently unable to definitively distinguish between the possibility of SNRs from more-massive progenitors being more difficult to identify, and very high-mass stars not producing SNe.

In addition to difficulties in detection due to obscuration, any inherent differences in luminosity or lifetime of SNRs as a function of progenitor mass could also introduce some bias. These properties are dependent upon both the energy of the SN and the density of the medium into which the SNR expands. If one would appeal to this as an explanation for the lack of high mass progenitors, this would imply that SNR from higher mass stars have inherently lower luminosities and/or inherently shorter lifetimes. At present, we have no real ability to test either of these cases with our methodology. If there are significant differences in SNR frequencies or properties as a function of progenitor mass, then naturally we are no longer dealing with an unbiased sample of all exploded massive stars, and our conclusions about the SNe progenitor mass distribution are weakened.