The Masses of Supernova Remnant Progenitors in M33

Using resolved optical stellar photometry from the Panchromatic Hubble Andromeda Treasury Triangulum Extended Region (PHATTER) survey, we measured the star formation history (SFH) near the position of 85 supernova remnants (SNRs) in M33. We constrained the progenitor masses for 60 of these SNRs, finding the remaining 25 remnants had no local SF in the last 56 Myr consistent with core-collapse SNe (CCSNe), making them potential Type Ia candidates. We then infer a progenitor mass distribution from the age distribution, assuming single star evolution. We find that the progenitor mass distribution is consistent with being drawn from a power-law with an index of $-2.9^{+1.2}_{-1.0}$. Additionally, we infer a minimum progenitor mass of $7.1^{+0.1}_{-0.2}\ M_{\odot}$ from this sample, consistent with several previous studies, providing further evidence that stars with ages older than the lifetimes of single 8 $M_{\odot}$ stars are producing supernovae.


INTRODUCTION
Stellar evolution theory predicts that stars with a mass above ∼8 M will end their lives as core-collapse supernovae (CCSNe; Woosley et al. 2002). The precise value of this lower limit has been the subject of multiple studies, some of which have found that red supergiants with masses as low as 7 M can be progenitors for Type IIP SNe (Smartt et al. 2009;Fraser et al. 2011;Jennings et al. 2012).
In addition to the uncertainty surrounding the lower mass limit of supernova (SN) progenitors, there has been growing evidence that not all stars with masses >8 M experience a canonical CCSNe (e.g., Smartt et al. 2009;Van Dyk 2017). For example, stellar evolution theory suggests that there is an upper mass cutoff of ∼30 M for Type II SNe (Massey et al. 2017). Observations, however, point to a new smaller upper limit of 17 − 19 M (Smartt 2015). This discrepancy, along with the lack of high-mass progenitors in observations, was dubbed "the red supergiant problem" by Smartt et al. (2009), who first argued that the discrepancy is statistically significant. More recently, Davies & Beasor (2020) suggested that some of the discrepancy can be attributed to the steepness of the luminosity distribution of red supergiants as well as the small sample size. They suggest the problem may not be as significant as previously believed, but also find a similar upper limit of 19 +4 −2 M . Above this mass limit, some stars may instead collapse straight into a black hole (Pejcha & Thompson 2015). These "direct collapse" systems would not produce a visible SN or leave behind a visible remnant. The first observed black hole formation candidate occurred in the star-forming galaxy NGC 6946 (Gerke et al. 2015;Adams et al. 2017). Murphy et al. (2018) and Koplitz et al. (2021) used the local stellar populations of the vanishing star to determine that the progenitor was likely ∼10.6 Myr old, which a single star progenitor points to an initial mass of ∼17 M .
Further progress on understanding the fates of massive stars requires increasing the number of CCSNe pro-genitors with mass constraints and expanding the measured distribution of progenitor masses to wider ranges of galaxy properties. The traditional method for determining the mass of SNe progenitors is by directly imaging the progenitor stars (e.g. Smartt et al. 2003Smartt et al. , 2004Van Dyk et al. 2003;Li et al. 2006;Gal-Yam et al. 2007;Kilpatrick et al. 2021). This technique requires high resolution (better than ∼0. 1) images of the SN site both before and after the event, which involves a large amount of serendipity. The difficult requirement of having spatially resolved photometry of the location before the explosion has resulted in only 34 SNe having their progenitor masses determined by this method, along with 40 upper limits constrained (Van Dyk 2017; Kilpatrick & Foley 2018;Van Dyk et al. 2018;O'Neill et al. 2019;Kilpatrick et al. 2021;Tinyanont et al. 2022;Vazquez et al. 2022). While the number of cataloged SNe has increased in recent years (e.g. Guillochon et al. 2017;Holoien et al. 2019), few of these SNe have had their progenitor constrained due to insufficient precursor imaging.
An alternative method, which does not require preexplosion images, uses an age-dating technique of the stellar populations surrounding an SN event (Gogarten et al. 2009;Murphy et al. 2011). This technique leverages the stellar populations surrounding an SN to measure the local star formation history (SFH) by finding the model age distribution that best fits the colormagnitude diagram (CMD) of the resolved local stars. By assuming the progenitor star belongs to the median population near the event, we are able to place statistical constraints on the age of the SN progenitor. We can then infer the most likely mass of the progenitor by assuming that it was the most massive star that survives to that age according to the models.
This age-dating technique was shown to be a reliable way to infer progenitor ages for distances out to ∼8 Mpc (Murphy et al. 2011). Assuming only stars with masses 7 M become CCSNe requires photometry that is sensitive to populations as old as 56 Myr (Girardi et al. 2002). Because the technique does not require precursor imaging, it can be applied to any location where an SN has occurred in the recent past, including any known SN remnants (SNRs). As a result, several previous works have shown that most young stars within 50 pc of an SN event are associated with the progenitor (Bastian & Goodwin 2006;Badenes et al. 2009;Gogarten et al. 2009;Jennings et al. 2012;Williams et al. 2014). For example, this technique was used to constrain the masses of SNR progenitors in the local star-forming galaxies M31 (Jennings et al. 2012), NGC 6946 (Koplitz et al. 2021), as well as the Magellanic Clouds (Badenes et al. 2009;Auchettl et al. 2019). This technique has also been used to constrain the mass of observed CC-SNe (Williams et al. 2014Díaz-Rodríguez et al. 2021;Koplitz et al. 2021). Progenitor masses in M83 have also been constrained, including one with a most likely mass of 59 M whose errors exclude ages older than 8 Myr, the highest mass progenitor inferred from the technique to date (Williams et al. 2019).
M33, or the Triangulum Galaxy, is an excellent target for our technique. It is nearby, relatively face on (i = 56 • ; Zaritsky et al. 1989), and is known to host over 200 SNRs (Long et al. 2010;Lee & Lee 2014a). Jennings et al. (2014), hereafter J14, have already applied this technique to 33 SNRs in M33, finding that the distribution was well fit by power-law distributions with indices that were significantly steeper than a standard Salpeter initial mass function power-law index of −2.35 (Salpeter 1955). However, their analysis in M33 was limited by the heterogeneous set of archival Hubble Space Telescope (HST) images available. This heterogeneous coverage resulted in inconsistent filter coverage and photometric depths between observations containing SNRs. Furthermore, they did not fit a separate field star component to their ages, which could have resulted in age biases. Here, we follow up on their work using the deep, uniform coverage provided by the Panchromatic Hubble Andromeda Treasury Triangulum Extended Region (PHATTER) survey ) as well as updated fitting techniques.
The analysis we present here takes advantage of the work by Díaz-Rodríguez et al. (2018), hereafter DR18, who developed a Bayesian hierarchical analysis capable of constraining the progenitor mass distribution index with an improved method for accounting for background effects as well as the minimum and maximum mass at which a star is able to undergo a CCSNe event from a set of SFHs. They reanalyzed the SFHs from J14 as well as those from Lewis et al. (2015) which correspond to likely SNRs from Lee & Lee (2014b), finding a progenitor mass index closer to, but not consistent with, a Salpeter index (−2.96 +0.45 −0.25 ). This combined M31 and M33 distribution pointed to a minimum mass of ∼7.3 M and a maximum mass of >59 M . However, they found the SFHs from M31 led to a Salpeter progenitor mass distribution index (−2.35 +0.36 −0.48 ) with a minimum mass of 6.5 M and a maximum mass of >46 M .
In this paper, we take an updated look at the ages of SNR progenitors in M33 using resolved stellar photometry from the PHATTER survey. Our larger sample and more homogeneous photometry catalog allow us to compare different fitting methods and quantify the impact these changes have on the age and mass results. Ad-ditionally, we compare our custom SNR-centered SFHs to those measured by Lazzarini et al. (2022) in grids, allowing us to determine whether grid SFHs are sufficient for inferring a progenitor age and mass. The rest of the paper is outlined as follows: Section 2 details our SNR source catalog, as well as how our SFHs were measured. Section 3 presents our progenitor age and mass estimates. In Section 4, we discuss our constraint on the lower mass limit for CCSNe as well as the results of Kolmogorov − Smirnov (KS) tests on our observed distribution, then compare our results to similar studies in the literature. Finally, Section 5 provides a short summary of our results. Throughout this paper, we assume a distance to M33 of 859 kpc (de Grijs et al. 2017).

DATA AND ANALYSIS
Our technique has two main data requirements. First, we need to know the locations of past SN activity. Second, we require resolved stellar photometry of the current populations within 50 pc of the SNe, as stars tend to remain spatially correlated within about 100 pc of their siblings for about 100 Myr, even if the cluster is not gravitationally bound (Bastian & Goodwin 2006). Using these, we can measure the star formation rate as a function of lookback time, known as the SFH, at each SNR location. The SFH provides the age distribution of the stars near each SN. We then apply this age distribution to constrain the age and mass of the progenitor star. We detail each of these steps below.

SNR Locations
For the locations of past SN activity in M33, we take the locations of SNRs from the catalogs of Long et al. (2010) and Lee & Lee (2014a), hereafter L10 and LL14, respectively. L10 identified candidates based on their X-ray spectrum as well as having [S II]:Hα ratios ≥0.4. The candidates in LL14 were identified based on their lack of blue stars, remnant morphology, and [S II]:Hα ratios ≥0.4. Of the 137 SNR candidates in L10, 120 are included in LL14's catalog of 199 candidates. The remaining 17 locations were classified as likely superbubbles or H II regions, leading LL14 to exclude them from their final catalog. Of these 17 potential SNRs, 4 (L10-043, L10-050, L10-079, L10-098) are within the PHATTER survey footprint. We include these 4 locations in our catalog since they may be SNRs located within larger star forming complexes. Of the 199 candidates from LL14, 81 reside in the PHATTER footprint, leading to our catalog of 85 SNR candidates.
In addition to the SNR locations, we produced 2 control catalogs of locations not associated with SNRs. The first sample is 85 locations randomly distributed within the PHATTER footprint. The second sample is 2500 random draws of the grid SFHs from Lazzarini et al. (2022), which do not contain an SNR. Differences between the random and SNR samples provide additional evidence that the stellar populations near the SNRs are likely related to the progenitors, and not chance spatial coincidences of young stars (see Section 4.3 for details).

Photometry
Once we had determined the historical SN locations, the second requirement was resolved stellar photometry at those locations. This photometry was obtained from the PHATTER survey . The survey measured resolved stellar photometry for 22 million stars within M33 in optical (Advanced Camera for Surveys F 475W and F 814W ), near-ultraviolet (Wide Field Camera 3 (WFC3) F 275W and F 336W ), and nearinfrared (WFC3 F 110W and F 160W ) bands. Our photometry is derived from the optical images (F 475W and F 814W ) of the PHATTER survey, rather than measuring photometry in all 6 bands simultaneously. We took samples from this photometry catalog for each SNR location and each control location, with the samples consisting of all of the stars within 50 pc (12 ) from the SNR or random position. We also collected samples of the widespread young populations surrounding each SNR from 50 to 1000 pc (12 to 4 ). These "background" samples allow us to identify young populations unique to the region containing the SNR.
To fit stellar evolution models to the photometric data, we require artificial star tests (ASTs) to correctly model the photometric completeness and uncertainty as a function of color and magnitude. We used the ASTs from these data that were created by Lazzarini et al. (2022), who used them to measure grid SFHs in M33, as discussed in Section 2.5. These tests are obtained by adding stars of known flux to an image and blindly rerunning the photometry routine to measure the photometric bias, uncertainty, and completeness as a function of color and magnitude when fitting models to the data. This is done at least 50,000 times within a region of interest. Williams et al. (2017) and Koplitz et al. (2021) found that one set of artificial stars could be used for all locations of similar stellar density, rather than creating a set for each location. This greatly reduces the computation time required. Lazzarini et al. (2022) used this technique to optimize the number of ASTs that needed to be created. Since we are using the same photometry catalog as Lazzarini et al. (2022), we are able to use the same ASTs when analyzing the SNRs in our catalog. These tests and the optical photometry catalog are described in further detail in Lazzarini et al. (2022).

CMD Fitting
Once we had the photometry and ASTs necessary to study each SNR location, we used the CMD fitting program MATCH (Dolphin 2002(Dolphin , 2012(Dolphin , 2013 to measure SFHs near the SNRs in our catalog. MATCH has been used to constrain the age of SN progenitors (e.g., Jennings et al. 2012;J14;Williams et al. 2018Williams et al. , 2019Koplitz et al. 2021) and SFH for nearby galaxies (e.g. Williams et al. 2009;Weisz et al. 2014;Skillman et al. 2017). MATCH fits the observed CMD using the PARSEC stellar evolution models (Bressan et al. 2012). For each model age and metallicity, it creates a model CMD by assuming a Kroupa initial mass function (Kroupa 2001). It then finds the highest likelihood linear combination of those models that provides the best fit to the observed CMD using a maximum likelihood estimator and taking into account the bias, uncertainty, and completeness of the photometry as determined by ASTs. This combination of models yields the distribution of ages and metallicities for the stars in the observed CMD, which we refer to as the SFH of the region.
Below, we provide a brief description of our technique for running MATCH. A more detailed account of the process can be found in Koplitz et al. (2021), which is identical to how we ran MATCH here. In short, for each SNR location, we fit the CMD of the resolved stellar photometry with a grid of model CMDs generated from the PARSEC stellar evolution models (Bressan et al. 2012). Our model grid had time bins of size 0.05 dex from log 10 (t/yr) = 6.6 − 8.0 while bins of size 0.1 dex were used from log 10 (t/yr) = 8.0 − 10.2. Since M33 is known to have a subsolar metallicity (e.g., Barker et al. 2011), we limited the metallicities MATCH applies to the model grid to be −0.5 ≤ [Fe/H] ≤ 0.1 using the zinc flag. Multiple massive stars can reside in the same location on CMDs even though they have different metallicities. As a result, using the zinc flag forces MATCH to use models for the young stars that are within the known metallicity range of M33.
As in Koplitz et al. (2021), our model also includes a "background" or "contamination" CMD of the stars in an annulus between 50 and 1000 pc (12 − 4 ) from the SNR. The contamination CMD is scaled to the size of our regions before fitting. This allows us to identify young populations that are sparse in the surrounding field and more heavily weight the populations concentrated within the regions being fit.
Furthermore, the fitting routine accounts for the effects of dust on the photometry. We allowed MATCH to find the combination of reddening parameters along with the combination of ages and metallicities, which provided the best fit to the observed CMD. Since young populations are often found in dusty regions, MATCH applies three types of extinction to the model CMDs when fitting the stellar populations. The first, A V , is the total foreground extinction over the full region. The second, dA V , is the extinction spread due to the stars along the line of sight. The third, dA V Y , is additional differential extinction added to populations younger than 100 Myr old. The default dA V Y value of 0.5 was used. To determine A V and dA V for an SNR, we fit a range of possible values at the location. We allowed A V to be between 0.05 and 1.00 in steps of 0.05 while dA V could be between 0.0 and 2.0 in steps of 0.2.
On average, our locations returned an A V value of 0.30, higher than the Schlafly & Finkbeiner (2011) value of 0.114. This higher A V is not surprising given that MATCH takes into account the Milky Way and M33 reddening, whereas Schlafly & Finkbeiner (2011) only account for the Milky Way. The vast majority of dA V values in our sample were 0.0, meaning the default differential reddening for the youngest stars (dA V Y = 0.5) was sufficient to account for differential reddening in most cases.

Uncertainty Estimation
Random and systematic uncertainties are inherent to fitting stellar models to CMDs. Most of the random uncertainties in our fits arise from photometric errors as well as the number of stars used to determine the most likely progenitor age. The systematic uncertainties are from any deficiencies present in the stellar evolution models used during the SFH fits. Lazzarini et al. (2022) have shown that there is good agreement between model sets for fits to young ages, and that the random uncertainties dominate the error budget in these fits. Thus, we use the random uncertainties determination to estimate the uncertainties in our SFHs.
To estimate our random uncertainties, we used the hybridMC tool within MATCH (Dolphin 2013). This task uses a hybrid Monte Carlo algorithm to accept or reject potential SFHs around the best fit SFH based on likelihood. We report the narrowest 68% of the distribution of accepted SFHs that decreases with look-back time in columns (4) and (5) of Table 1. A detailed description of how our uncertainties are estimated can be found in Section 2.6.

SFHs from Previous Work
Recently, Lazzarini et al. (2022) published recent SFH maps of the PHATTER region of M33. They used the same PHATTER optical photometry to measure the SFH of M33's inner disk in a grid of 100 × 100 pc (24 × 24 ) cells, which they have released to the community. Lazzarini et al. (2022) largely used the same MATCH fitting technique as we have, but there were a few differences. In their analysis, time bins of size 0.1 dex were used for all bins (log 10 (t/yr) = 6.6 − 10.2). Since they were measuring the total amount of star formation in each location, and not attempting to isolate very localized populations, a contamination CMD was not included during their fits.
Being able to constrain progenitor masses using such a grid of spatially resolved age distributions would be very powerful, since it would avoid having to access the original photometry and ASTs and run custom fitting for each SNR location. Thus, we also attempted to age date the SNR locations using this grid of published SFHs by assigning an SFH from Lazzarini et al. (2022) that corresponded to the location of each SNR in our sample. We then compare the ages and masses of custom fits to those taken from a less optimized, but more easily accessible, source.

Constraining Progenitor Mass
The next step in constraining the masses of SNR progenitors is to convert the recent SFH from MATCH into a probability distribution for the age of the progenitor. This calculation is done by determining the fraction of the total stellar mass present in each age bin. We take this fraction to be equal to the probability that the progenitor is associated with that age. We also take the error on that fraction as the error on the probability. We provide an example of such a probability distribution in Table 1.
The age probability distribution presented in Table 1 is for the SNR LL14-060. Similar tables for each SNR with SF in the last 56 Myr are combined into one and made available in the online supplemental material.
While the age probability distribution derived from the SFH is the most complete constraint on the progenitor age, we also provide a single progenitor mass estimate with uncertainties. This age simplifies the mass inference, as well as comparisons with other measurements and mass distribution analysis. To derive the most likely progenitor age, we use the SFHs and uncertainties produced by MATCH to calculate the median age of the stellar populations younger than 56 Myr surrounding each SNR. We then take that age as the most likely progenitor age. We determine the uncertainties on the median age as follows. We recalculate the median age a million times by accounting for the uncertainties and resampling the SF rates in each time bin, then determine the narrowest 68th percentile of this distribution of ages that contain the best fit. We use a 56 Myr cutoff for our SNR-centered SFHs, rather than the 50 Myr used by other works, because of the results of our Bayesian analysis presented in Section 4.1. A 56 Myr (log(t/yr) = 7.75) cutoff is not possible for the grid SFHs since Lazzarini et al. (2022) ran MATCH with time bins of 0.1 dex. As a result, we must decide whether to use a 50 or 63 Myr cutoff (log(t/yr) = 7.7 or 7.8) for the grid SFH samples. We adopt a 50 Myr cutoff as this limits the number of contamination populations being included in our analysis. To infer the progenitor mass for each age bin, we assume that the SNR progenitor is the highest surviving mass on the PARSEC stellar isochrone (Bressan et al. 2012).
We present an example of our progenitor age fitting results for SNR LL14-160 in Figure 1. Similar summary plots are available for all of the SNR locations in the online supplemental material associated with the paper.
Past studies have shown that progenitor masses estimated from the SFHs produced by MATCH are consistent with estimates from other techniques (e.g., Jennings et al. 2012;Williams et al. 2019;Koplitz et al. 2021). DR18 found that their combined M31 and M33 distribution pointed to a minimum mass for CCSN progenitors of ∼7 M , which corresponds to an age of ∼50 Myr assuming single star evolution. Populations older than this are more likely to be unrelated to the SNR since they have had more time to distance themselves from their parent cluster. Older stars in binaries have been shown to be possible SN progenitors (e.g., Xiao et al. 2019); however, our current inference from age to mass requires that we assume single star evolution. Fortunately, this assumption should not impact our age constraints, which come from the surrounding population, but it could significantly impact our conversions between age and progenitor mass if the progenitor system was a mass-exchanging binary.

RESULTS
We present and provide the progenitor mass results from our own custom SFHs. We then compare to results that we obtain from previously published SFHs, as well as control samples and results from SNR studies of other nearby galaxies. These comparisons suggest that custom SFHs with a contamination CMD included in the fit to account for the more widespread populations are required to isolate the ages of the stars most likely associated with each SNR.

Comparing Grid SFHs to SNR-Centered SFHs
We present our progenitor mass constraints for the SNRs in our catalog in Table 2 and compare the resulting age distributions in Figure 2, which reveals that the masses from Lazzarini et al. (2022) are systematically lower than our custom measurements. For 42 of the 85 locations in our catalog (∼49%) the best fit masses were not consistent with each other. KS tests between these samples returned a p−value of 0.11, suggesting we cannot rule out that they are from the same parent distribution.
To isolate the cause of the observed difference, we reran our custom fits without including a contamination component, which returned a distribution similar to the one from the Lazzarini et al. (2022) grid SFHs. Performing KS tests between the SNR-centered distributions returned a p−value of 0.09 while 0.27 was returned when comparing the grid distribution to the SNR-centered without a contamination CMD sample. Figure 3 is a histogram comparing the distribution of progenitor masses that resulted from using the grid SFHs as well as the SNR-centered SFHs with and without a contamination CMD. Each distribution is normalized such that they integrate to one. The overall distribution from the grid SFHs is similar to that of our centered SFHs without a contamination CMD, which is expected given that both SFHs were fit without a contamination CMD and the sample populations overlap.
Even though none of the distributions contain a progenitor that excludes masses <20 M , the grid SFHs produced 10 locations consistent with being more massive than 20 M while the SNR-centered SFHs returned 15 with a contamination CMD and 9 without one. A similar fraction of locations with masses between 7 − 15 and 15−25 M were found in the grid and SNR-centered without a contamination component distributions (86%, 13% and 86%, 11% respectively). These show that the inclusion of the contamination CMD impacts the resulting distribution the most, though the high-precision custom location does play a large role.

Type Ia Candidates
Of the 85 locations in our catalog, we classify 25 as Type Ia candidates. Zapartas et al. (2017) showed that binaries with ages down to 200 Myr can produce delayed CCSNe; however, these systems have had enough time to move a significant distance away from their parent cluster, making the SF we measure older than ∼56 Myr likely contaminated by nearby populations that are not associated with the SN event. Thus, any location without SF in the last ∼56 Myr we classify as Type Ia candidates since our technique cannot reliably determine the progenitor age beyond this.
Including contamination CMDs in our SFH fits forces MATCH to only fit for SF above any background young stellar populations. While this requirement can be helpful in isolating populations more likely to be associated with an SNR, it can also cause some SNRs to be classified as Type Ia candidates when they are actually Type II or Type Ibc in origin, because their associated young population may be too similar to that of the larger surroundings. To check how many of our Type Ia candidates could actually be CCSNe, we can use our results from the Lazzarini et al. (2022) SFHs, which measured the total star formation in each location. The SNRs with mass estimates in column (9) of Table 2 but without a constraint in column (7) are less likely to be Type Ia in origin, as there are relatively high-mass populations nearby, just not above the background level. Of the 25 Type Ia candidates from our SNR-centered SFHs, only LL14-103's grid SFH contained no SF within the last 50 Myr, making it our best Type Ia candidate. The other 24 Type Ia candidates had young stellar populations present but not in sufficient quantities to be detected above the larger surroundings, making them weaker Type Ia candidates. These results suggest a Type Ia fraction between 1 − 29%. While this is not a tight constraint, it is consistent with the ∼15% expected for late-type spirals (Li et al. 2011).

DISCUSSION
Our progenitor age distributions probe the minimum mass at which CCSNe can occur, how SNe are spatially distributed in the disk of M33, and the power-law index of the progenitor mass distribution for the galaxy.

Mass Limits for CCSNe
Using the Bayesian hierarchical analysis developed by DR18, we use our SNR age sample to provide a constraint on the maximum age at which stars undergo CCSNe, t max . Our analysis was sensitive to the assumed minimum age for CCSNe, t min . To account for this, we fit our distribution assuming t min values of 6, 9, 10, 12, 15, and 18 Myr, with each returning similar results. We report the t min = 15 Myr fit since this is the lowest t min value that stabilized the returned progenitor mass distribution slope, finding 54.3 +3.8 −2.0 Myr as the best fit t max which corresponds to a M min of 7.1 +0.1 −0.2 M . Figure 4 shows the distribution of t max (left) and M min (right) returned by the Bayesian analysis for the fit with t min = 15 Myr.
The analysis also attempts to constrain the upper mass limit for CCSNe and the progenitor mass distribution index when fitting a distribution. Our high-mass progenitors, however, have large error bars that prevented the analysis from converging on a rigorous best fit value for the upper mass limit. Since the upper mass is degenerate with the distribution index, it also did not return a reliable index. We estimate the progenitor mass distribution index in Section 4.3 using an alternate technique, since it may be of interest to the community.

Spatial Distribution of Progenitor Masses
To investigate the spatial distribution of SNRs in M33, we plot the locations of our catalog on an Hα image taken with the WIYN 0.9m telescope ( Figure 5). The progenitor mass and most likely SNe type are indicated by the color and symbol, respectively. Locations for which we have mass constraints in column (7) of Table 2, i.e., we were able to measure SF within the last 56 Myr above the background level, are shown as circles. Progenitors with masses <9 M are white, masses of 9 − 12 M are red, masses of 12 − 15M are orange, masses of 15 − 20 M are yellow, and masses >20 M are blue. Our Type Ia SNe candidates are shown as squares, where the color indicates the best fit progenitor mass from the grid SFH that the SNR resides in from Lazzarini et al. (2022). The colors show the mass that could have produced a CCSNe at the location, though these are less likely to be CCSNe than the colored circles due to the lower level of young SF. The coloring depicts the same mass ranges as the circles, with the addition of black indicating the location of LL14-103, our best Type Ia candidate.
Our entire catalog mostly traces the Hα emission and spiral arms of M33. There are many squares (Type Ia candidates) inside star forming regions throughout the galaxy, indicating that young populations are present, just not enough to be detected in fits that include a contamination component. In these cases, fitting without a contamination CMD (i.e., fitting the full population) often finds some massive stars in the region, whereas the fit including a contamination CMD finds no such populations.

Progenitor Mass Distribution Power-Law Index
While the Bayesian hierarchical analysis of DR18 was not able to converge on a power-law index for the progenitor distribution due to the uncertainties at very young ages, it may still be of interest to determine the closest power-law representation of our most likely progenitor masses. To determine this value, we use KS tests to determine the likelihood the locations in our catalog with young populations (<56 Myr) are drawn from various power-law distributions. We compared the data to power-law indices between −6.0 and 0.0, in steps of 0.1, and report the most likely index. To estimate the uncertainties on the index, we employee a bootstrap analysis in which we sample the uncertainties on each mass 1000 times. We then find the indices that return p−values ≥0.05 (∼95% confidence) and report the extremes as our limits.
Performing this analysis on our full catalog of SNRs indicates the progenitor mass distribution is best matched by a power-law with an index of −2.9 +1.2 −1.0 , which does contain the Salpeter index of −2.35 (Salpeter 1955). The best fit index has a p−value of 0.23. Running the progenitor mass distribution from the grid SFHs through this same analysis found that the sample was best matched by an index of −3.8 +1.8 −0.4 , significantly steeper than our SNR-centered catalog though still consistent. This is not surprising given that, as discussed in Section 3.1, L10-043 was the only progenitor found to be more massive than 25 M in this sample.
Our power-law indices were estimated using only the locations that contained SF at ages younger than ∼56 Myr. To estimate the impact that removing locations without young SF has on our indices, we refit the progenitor mass distribution index of our SNR-centered sample while adding in the grid SFH progenitor mass for locations that did not contain young SF in our custom fit. This combined sample was best fit by a −3.2 +1.3 −0.6 index, which is consistent with both the SNR-centered and grid SFH indices. This indicates that removing locations without young SF does not have a large impact on the returned progenitor mass distribution index. Figure 6 shows our ranked progenitor mass distribution. The red points indicate the progenitor mass of each SNR with a constraint in column (7) of Table 2, with uncertainties shown as red lines. Overplotted as gray lines are 50 draws from a power-law distribution with an index of −2.9 +1.2 −1.0 , our best fit index.

Control Sample
As mentioned in Section 2.1, we also performed our analysis on control samples, random locations that did not contain SNRs. We compared our SNR results to these control results to determine if the SNRs are indeed affecting our results. We discuss both the control sample for the mass estimates based on custom SFH measurements and the control samples for mass estimates based on Lazzarini et al. (2022) SFH measurements below.
Our first control sample, containing randomly drawn locations in the PHATTER footprint, returned fewer progenitors with masses >20 M (5 in the control sample and 9 in our catalog). There were also significantly more (33, ∼39% of the sample) Type Ia candidates (i.e., locations with no significant recent SF above what is present in the contamination CMD) than our SNR-centered distribution with a contamination CMD (25, ∼30% of the sample). Both of these suggest that the regions in this control sample contained, on average, older populations than those found near SNRs. The locations in this sample that contained SF at ages younger than 56 Myr were best fit by a power-law index of −4.9 +3.2 −0.2 , which is consistent with our contamination CMD and the grid distributions. While the uncertainties do overlap, this can likely be attributed to the amount of widespread SF within the PHATTER footprint of M33. Comparing this random sample to our SNR-centered sample returned a p−value of 0.08, suggesting these sample are only marginally consistent with being drawn from the same parent distribution.
Of the 2500 random draws in the grid control sample, ∼70% were consistent with the grid power-law index, with the remaining ∼30% resulting in steeper indices. We found the median index to be −4.1 ± 0.5, which includes the grid SFH sample index of −3.8 but excludes the SNR-centered index of −2.9, though the uncertainties do overlap. Of the 2500 draws, 1492 (∼60%) resulted in p−values ≤0.05 when compared to our grid sample. No p−values ≥0.05 were found when compared to the SNR-centered sample. Additionally, our grid sample only contained 1 location without recent SF (LL14-103) whereas only 3 of the grid control draws contained as many or fewer such locations, meaning that >99% of random draws had more locations without young stars present than we find in locations containing SNRs. These results show that the grid SFHs that contain SNRs do differ from those that lack an SNR, with the similar power-law indices likely being explained, again, by the amount of widespread SF in M33.

Comparison to J14 and DR18
J14's catalog contained 33 SNRs in M33, of which 28 are in the PHATTER footprint. Our SNR-centered progenitor mass estimates were consistent with those found by J14 for 16 of these 28 sources. Of the 12 sources that were not consistent between our estimates and those in J14, we identify 8 Type Ia candidates. J14 found their full distribution of 33 SNRs was best fit by an index of −3.8 +0.5 −0.4 . Using our SNR-centered SFHs of the 28 overlapping locations, our analysis pointed to the distribution being well matched by a power-law index of −3.1 +1.2 −1.1 , which is flatter than what J14 found. Their steep index could be from the few progenitors with masses >20 M in their sample, which can likely be partially attributed to the low number of SNRs in their sample. Poisson fluctuations for small numbers can randomly vary to zero quite easily. Additionally, J14 did not include a contamination CMD when fitting, which may have biased their estimates toward older, less massive populations and reduced their number of Type Ia candidates. We have shown in Section 3.1 that CMD-based age dating returns more low-mass progenitors when no background CMD is included in the fitting. DR18 constrained the minimum mass for CCSNe to be 7.32 +0.12 −0.14 M using the J14 measurements, which is consistent with the minimum mass we identified, suggesting that this cutoff may be the most reliable parameter returned from the SNR sample.

Comparison to Other Galaxies
In addition to the SNRs in M33, J14 also constrained the progenitor mass distribution for 82 SNRs in M31 using photometry from the PHAT survey. Their KS test analysis found that this distribution was best matched by a power-law index of −4.4 +0.4 −0.4 , which is not consistent with our distribution. Interestingly, one major difference between J14 and DR18 is that DR18 allow for a uniform background distribution when fitting for the progenitor mass distribution parameters. This should have a similar effect to the use of a contamination CMD during the fitting process in that it accounts for the possibility that some of the measured SF may not be associated with the SN. With the inclusion of the uniform background distribution, DR18 constrained the distribution index to be −2.35 +0.36 −0.48 and found the minimum mass for CCSNe to be 6.5 +0.6 −0.2 M using 62 SNRs in M31. Both of these measurements are consistent with what we have found in M33. Katsuda et al. (2018) gathered the progenitor masses for 40 SNRs in the Milky Way and both Magellanic Clouds from the literature, which were estimated using chemical abundances. They updated many of the measurements using Fe:Si ratios and found a progenitor mass distribution consistent with both a Salpeter index and our measured index for M33. Auchettl et al. (2019) examined 23 SNRs in the Small Magellanic Cloud, finding that 22 were likely core collapse in origin. They report the likelihood that the mass of each progenitor is between 8 − 12.5, 12.5 − 21.5, and >21.5 M assuming single and binary evolution. Regardless of single or binary evolution, 70% of their progenitors had the highest likelihood in the most massive bin, whereas 9% of our SNRs were found to have progenitor masses in this range. The large number of high-mass progenitors led to their distribution being well matched by a power-law index of −1.84, though the uncertainties are consistent with a Salpeter index. A possible reason given by Auchettl et al. (2019) for the top-heavy distribution is that the Small Magellanic Cloud has a lower metallicity than M33. It has been shown that lower metallicity gas is more likely to produce a top-heavy stellar distribution (e.g., Bromm & Larson 2004;Marks et al. 2012). Williams et al. (2019) constrained the progenitor mass of 199 SNRs in M83 using our technique. They found that the progenitor mass distribution was well matched by a power-law index of −2.9 +0.2 −0.7 . A KS test between their M83 distribution and ours resulted in a p−value of 0.04, suggesting their parent distributions may differ, possibly due to the higher star formation intensity or lower metallicity of M33. Koplitz et al. (2021) measured the progenitor mass of 169 SNRs, 8 historically observed SNe, and NGC6946-BH1, the first black hole formation candidate, in the galaxy NGC 6946 using our technique. They found that gas emission impacted their broad V -band photometry, which biased some of their mass estimates. As a result, they only included the 46 sources that were least likely to be biased when constraining the progenitor mass distribution index. In this sample, they found 24% with masses ≥20 M , while we have 11% progenitors in our catalog with similar masses. They found their distribution was best fit by an index of −2.6 +0.5 −0.6 , which is consistent with our measured index. KS tests between their preferred sample and our distribution resulted in a p−value of 0.2. Figure 7 compares our distribution of progenitor masses to those in M83 (Williams et al. 2019) and the preferred sample in NGC 6946 (Koplitz et al. 2021). We normalize each individually so that they integrate to one. Each is dominated by the low-mass progenitors and those less massive than 25 M have similar overall shapes. These led to power-law indices that are consistent with each other. Our distribution is the only one that lacks any progenitors with mass ≥40 M . However, this could be the result of the small number of high-mass progenitors expected combined with the smaller number of SNRs in our sample.

SUMMARY
We constrained the progenitor age and mass of 60 SNRs in the nearby galaxy M33, or the Triangulum Galaxy, using an age-dating technique of the stellar populations near the SNRs. The remaining 25 showed no local SF within the past 56 Myr, making them potential Type Ia candidates. While it is possible that these candidates are binary systems producing delayed CCSNe with ages down to 200 Myr, our analysis is not able to reliably determine the progenitor age beyond ∼56 Myr.
Using the Bayesian hierarchical analysis developed by DR18, we constrained the maximum age for CCSNe to be 54.3 +3.8 −2.0 Myr which, assuming single star evolution, corresponds to a minimum mass of 7.1 +0.1 −0.2 M . A KS test analysis determined that the progenitor mass distribution of our full catalog was best matched by a powerlaw distribution with an index of −2.9 +1.2 −1.0 , which in-cludes the Salpeter index of −2.35. Our distribution is well populated by progenitors with masses 9 − 40 M . When using grid SFHs from Lazzarini et al. (2022), rather than SNR-centered regions with a contamination CMD included, the inferred progenitor mass was biased to lower values, with only 1 progenitor more massive than 25 M . There were also fewer Type Ia candidates when using the grid SFHs, 1 from the grid SFHs and 25 from the SNR-centered SFHs. Additionally, the progenitor mass distribution index that came from the grid SFHs was steeper than the index from the SNR-centered SFHs while KS tests between these samples returned a 0.03 p−value. Without a contamination CMD, our custom SFHs returned a similar distribution to the grid, finding a p−value of 0.14 between the two samples. The grid results differ from a random distribution of grid cells that do not contain an SNR, though not as strongly as our SNR-centered sample. The stronger difference from the overall background suggests that the custom SFHs with a contamination CMD provide a more robust constraint on the age and mass of SNRs than SFHs measured in grids, where the background populations are not taken into account and the SNR may be anywhere in the grid cell.
Previously, J14 used archival HST images to constrain the age and mass of 33 SNRs in M33. We present new age and mass estimates for 28 of these using the deep, uniform photometry from the PHATTER survey. Performing KS test analysis on the SNRs with updated mass estimates pointed to the distribution being well matched by a power-law index of −3.1 +1.2 −1.1 , which is consistent with the index J14 found for their full catalog and our sample of 85 SNRs.
Our normalized progenitor mass distribution is similar to that of M83 (Williams et al. 2019) and NGC 6946 (Koplitz et al. 2021). All the distribution are dominated by low-mass progenitors and have best fit powerlaw indices that are consistent with one another. KS tests between our sample and the preferred sample in NGC 6946 resulted in a p−value of 0.20, suggesting that these are likely drawn from the same parent distribution. A p−value of 0.04 was returned when performing KS tests between our sample and the SNRs in M83. A p−value just below 0.05 and the matching power-law indices means we cannot rule out that the samples are drawn from the same parent distribution.
Each of our distributions shows a sharp drop in the number of progenitors at ∼20 M . Few progenitors are found more massive than this, which coincides with the upper limits found by Smartt (2015) and Davies & Beasor (2020) for Type II SNe. It is possible that the reason we do not see many high mass progenitors is because not all experience a canonical CCSNe and instead collapse directly into a black hole (Pejcha & Thompson 2015). Now that the JWST has launched, similar studies will be able to leverage red supergiants to constrain the age of SN progenitors with higher precision and may resolve the "red supergiant problem". Support for this work was provided by NASA through grants GO-14610 and GO-15216 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555.  Bottom left: the best fit SFH and associated uncertainties produced by hybridMC. The cumulative fraction of stellar mass in each age bin is indicated by the red line. The gray shaded region depicts the 1σ uncertainties on the SFH. The tan region shows the probable age range, which is taken to be median population of the 1σ uncertainties. Bottom right: the differential SFH, showing the SF rates and uncertainties that correspond to the cumulative fraction plot in bottom left. The blue line indicates the SFH for the region, column (3) of  . Histogram comparing the normalized progenitor mass distributions from the SNR-centered SFHs with and without a contamination CMD as well as the distribution from the grid SFHs (Lazzarini et al. 2022).  (7) of Table 2 are shown as colored circles with masses <9 M in white, masses of 9 − 12 M in red, masses of 12 − 15 M in orange, masses of 15 − 20 M in yellow, and masses >20 M in blue. Type Ia candidates (those without an entry in column (7) of Table 2) are shown as colored squares, where the color indicates the grid progenitor mass from column (9) of Table 2. The coloring is the same as the circles with the addition of black indicating the location of LL14-103, our best Type Ia candidate. . Each bin is 2 M wide and each distribution is normalized such that they integrate to one.