The Progenitor Age and Mass of the Black Hole Formation Candidate N6946-BH1

We used pipeline-calibrated HST Wide Field Camera 3 (WFC3) UVIS images from several general observing campaigns. See Table 1 for the proposal IDs, PI, filters, and total exposure times. Figure 1 shows an image for the region surrounding N6946-BH1. The location of the progenitor for N6946-BH1 was at the center of the red circle, R.A. 20:35:27.56 and decl. +60:08:08.29. The GO-13392 images were taken on 2014 February 21, the GO-14266 images were observed on 2015 October 8, and the GO-14786 images were taken on 2016 October 26, which is 5–7 yr after the observations in which N6946-BH1 had disappeared (2009).

Zoom In Zoom Out Reset image size

We derived photometry from the images using a modified version of the PHAT pipeline, which is built around DOLPHOT (Dolphin 2002; Williams et al. 2014b). We selected stars with with signal-to-noise ratios S/Ns > 4, sharpness squared <0.15, and crowding <1.3. We combined the photometry from all epochs of imaging; the resulting data reach a depth of 27.64, 28.23, and 26.74 in F438W, F606W, and F814W, respectively. These correspond to the brightness of a solar-metallicity MS star with mass 12, 9, and 15 M_⊙, respectively (Marigo et al. 2017). To convert these magnitude depths into MS mass, we assume a distance to NGC 6946 of 7.83 ± 0.29 Mpc (see the derivation of the distance in Section 2.2) and total foreground extinction of A(V) = 0.94 (Schlafly & Finkbeiner 2011; Adams et al. 2017).

Gerke et al. (2015) measured the photometry for the progenitor of N6946-BH1 from HST archival images. The resulting magnitudes are 23.09 ± 0.01 in F606W and 20.77 ± 0.01 in F814W.

Figure 1 also shows regions for which we characterize the stellar populations. The red circle is centered on the vanisher. We adopt a radius of 100 pc for these regions. Early experiments on the size of the region (Gogarten et al. 2009) suggested that ∼50 pc is big enough to include a significant number of MS stars for age-dating but not too big to include too many stars from neighboring populations. However, that exploration was for one region of one galaxy, and the optimum size will likely differ depending upon the local stellar density. For now, we adopt 100 pc radii and leave a more thorough analysis of the dependence of size for future work. In addition to the region centered on the vanisher, we also infer the age of eight nearby 100 pc radius regions (blue circles), which allows us to consider scenarios in which the vanisher was ejected from a neighboring star-forming region. For a characteristic ejected velocity of ∼10 km s⁻¹ and a lifetime of ∼10 Myr, the vanisher could have traveled 100 pc before undergoing core collapse.

In Figure 2, we plot the F814W versus F606W–F814W (top panel) and F606W versus F438W–F606W (bottom panel) CMDs for the central region colocated with the vanisher. The star in the top plot represents the photometry of N6946-BH1 before it vanished in the optical bands (Gerke et al. 2015). Before N6946-BH1 vanished, it was the brightest star in the F814W filter in this region. In fact, it is the brightest star within 300 pc. Later, we find the SFH both with and without the progenitor; it made no real difference in the resulting SFH. For reference, we plot isochrones for six ages (7, 10, 14, 19, 25, and 50 Myr) assuming solar metallicity and a foreground reddening of A(V) = 0.94 and R(V) = 3.1 (Marigo et al. 2017).

Zoom In Zoom Out Reset image size

We infer the SFH for both F438W versus F438W–F606W and F606W versus F606W–F814W CMDs independently, but we only report the more accurate SFH for F606W versus F606W–F814W. In principle, one may use MATCH to infer an SFH based upon all three bands. However, MATCH is not able to calculate the uncertainty in the SFH based upon three bands. It calculates the uncertainty for two or four bands. We find that the SFH is a little more accurate when using the F606W and F814W bands. Furthermore, the progenitor is extremely red and only has upper limits for F438W (Adams et al. 2017). To give a broader sense of the CMD, we report both CMDs, but for inferring SFHs and ages, we primarily report ages using the F606W and F814W CMD.

When inferring the SFH, one must convert model luminosities into model fluxes; this requires an accurate distance to NGC 6946. NGC 6946 has had roughly 10 observed SNe in the last 50 yr, so deriving an accurate distance to this one galaxy will greatly improve the statistics of many local SNe. Here we describe our method to calculate the distance to NGC 6946 using the tip of the red giant branch (TRGB).

There are 32 distance estimates in the NASA/IPAC Extragalactic Database with a mean of 5.5 Mpc and a standard deviation of 1.5 Mpc. A variety of techniques were used to determine these distances; they include the brightest blue stars, Tully–Fischer, TRGB, and the expanding photosphere for II SN methods. The reported uncertainties on most of the distance moduli are of order 0.1. Gerke et al. (2015) and Adams et al. (2017) used a distance of 5.9 Mpc derived from the brightest blue stars method (Karachentsev et al. 2000). However, we note that Karachentsev et al. (2000) actually reported a distance of 6.8 Mpc for NGC 6946. The confusion in distance from Karachentsev et al. (2000) is understandable. The text of that manuscript only refers to the average distance to the NGC 6946 group (5.9 Mpc), but Table 4 of Karachentsev et al. (2000) explicitly gives the distance to NGC 6946 as 6.8 Mpc. More recently, Tikhonov (2014) used the expanding photosphere method and inferred a distance of 6.72 ± 0.15 Mpc. In any case, these distance measures are fairly uncertain compared to a well-calibrated TRGB distance estimate.

McQuinn et al. (2016) described a precise method for determining the distance from the TRGB. Our first step is to use HST images to calculate accurate photometry for a large number of stars in NGC 6946. The images were taken in the F606W and F814W bands by program HST-GO-14786 (PI: B. Williams). They contain 5470 s of exposure in F606W and 5538 s of exposure in F814W. They reach 27.5 in F814W and 28.6 in F606W. They were observed on 2016 November 02 pointed at R.A. = 308.8530225, decl. = 60.040649634. We derived F606W and F814W photometry from the images using a modified version of the PHAT pipeline, which is built around DOLPHOT (Dolphin 2002; Williams et al. 2014b). The TRGB location was measured using a Bayesian maximum-likelihood technique described in McQuinn et al. (2016), which is based on Makarov et al. (2006).

The measured TRGB magnitude is 26.004 ± 0.035 at a color of F606W − F814W = 1.71. The fit is clean; in particular, there are no secondary peaks that might indicate issues. With a total foreground extinction of A(V) = 0.94 (Schlafly & Finkbeiner 2011; Adams et al. 2017), the corresponding reddenings are A(F606W) = 0.88 and A(F814W) = 0.57. Therefore, the reddened-corrected TRGB is 25.44 with a color of 1.40. Rizzi et al. (2007) provided a calibration for the absolute magnitude for the TRGB: M(ACS F814W) = −4.06 + 0.20[(F606W − F814W) − 1.23]. Given the reddening-corrected color, the F814W absolute magnitude is M(ACS F814W) = −4.03. Sources of uncertainty in the Rizzi et al. (2007) calibration and its underlying horizontal-branch calibration (Carretta et al. 2000) add 0.05 mag of uncertainty, while an adopted 10% uncertainty in the foreground extinction adds another 0.05 mag of uncertainty (Schlafly & Finkbeiner 2011). Combined with the measurement uncertainty, we calculate a distance modulus of m − M = 29.47 ± 0.079. The corresponding metric distance is 7.83 ± 0.29 Mpc.

This new distance for NGC 6946 suggests that the vanisher was intrinsically 1.3 times more luminous than assumed in Adams et al. (2017). According to the Marigo et al. (2017) models, the larger distance would shift their mass estimate from ∼25 to ∼27 M_⊙. As a rough measure of systematics in the stellar evolution models, Adams et al. (2017) also estimated the mass using rotating models. Their SED-derived mass from the rotating models was ∼22 M_⊙. With the new distance, this would shift to ∼24 M_⊙.

We derive SFHs from the F606W–F814W CMD using the program MATCH (Dolphin 2002, 2012, 2013), which generates model CMDs that include the effects of observational errors (as characterizing the artificial star tests), foreground and internal dust extinction, and distance and then adjusts the SFHs to maximize the likelihood of the observed CMD. We run MATCH using logarithmically spaced age bins, Δ log₁₀(t yr⁻¹) = 0.05, and the left edge of the minimum age bin is log₁₀(t yr⁻¹) = 6.6). Using the TRGB distance, we fix the distance to be 7.83 Mpc.

We assume an extinction law with R(V) = A(V)/E(B − V) = 3.1 (Cardelli et al. 1989; O'Donnell 1994) and use the same galactic foreground reddening as Adams et al. (2017) of E(B − V) = 0.303 (Schlafly & Finkbeiner 2011) to fix the foreground extinction in MATCH to be A(V) = 0.94. However, we allow MATCH to fit for a distribution of internal extinctions; the default model in MATCH uses a top-hat distribution of extinction, with width dA_V. This is designed to approximate the behavior of young stars enshrouded in a layer of patchy dust. We infer dA_V for each region in Figure 1. Starting from the top left, the values for dA_V are 0.00, 0.19, and 0.00 (top row); 0.70, 0.39, and 0.00 (middle row); and 0.05, 0.17, and 0.00 (bottom row). Taking the internal extinction as a tracer of the presence of molecular gas, one might expect the youngest population to be in the region containing the vanisher and the region to its left in Figure 1.

The SFH provides the stellar mass formed as a function of age. To estimate the age probability distribution function (PDF), we assume that the probability is proportional to the mass of stars formed. If one knows the exact SFH, then the age PDF would simply be proportional to the best-fit SFH. However, MATCH provides a hybrid Markov chain Monte Carlo (MCMC) algorithm, which returns a distribution of SFHs given the photometric data of the surrounding stars. Therefore, one must derive the age distribution given the distribution of accepted SFHs.

The goal is to derive the marginalized age distribution P(τ), given the distribution of SFHs ({S_j(τ)}) that MATCH returns; j indexes each MCMC sample. When converting an SFH to a PDF, one must choose a maximum age for the normalization. A reasonable choice for this maximum age is the maximum age for core collapse, T_max. Given these assumptions and dependencies, the joint probability for observing a burst of star formation is P(τ, S, T_max). Marginalizing this distribution gives the age PDF

$$\begin{eqnarray}&&P(\tau )=\int P(\tau ,S,{T}_{\max })\,{dS}\,{{dT}}_{\max }.\end{eqnarray} \tag{ 1 }$$

To derive the joint probability distribution, we use the conditional probability theorem, which states that

$$\begin{eqnarray}&&P(\tau ,S,{T}_{\max })=P(\tau | S,{T}_{\max })P(S)P({T}_{\max }).\end{eqnarray} \tag{ 2 }$$

Here $P(\tau | S,{T}_{\max })$ is the likelihood for τ given the SFH, S, and maximum age in converting the SFH into a PDF, T_max:

$$\begin{eqnarray}&&P(\tau | S,{T}_{\max })=\displaystyle \frac{S(\tau )}{{M}_{\star }({T}_{\max })},\end{eqnarray} \tag{ 3 }$$

where ${M}_{\star }={\int }_{0}^{{T}_{\max }}S(\tau )\,d\tau$ is the total stellar mass formed in the last T_max years. Here P(S) represents the distribution of SFHs from MATCH hybrid MCMC runs, and P(T_max) represents a prior on T_max. Díaz-Rodríguez et al. (2018) recently found the maximum age for core collapse to be ${50.3}_{-0.5}^{+2.5}$ Myr, and one could use this distribution as the prior. However, this distribution is preliminary in that the Díaz-Rodríguez et al. (2018) analysis does not include the uncertainties in the SFH in the conversion from SFH to PDF. The work of this section will provide the foundation for a more thorough analysis. Therefore, for the purposes of this manuscript, we will employ the same prior that was used in previous studies (Gogarten et al. 2009; Murphy et al. 2011; Williams et al. 2014a); we will only consider ages below ${\widetilde{T}}_{\max }$ in the conversion from an SFH to a PDF. Hence, the prior for T_max is a delta function at ${\widetilde{T}}_{\max }$,

$$\begin{eqnarray}&&P({T}_{\max })=\delta ({T}_{\max }-{\widetilde{T}}_{\max }),\end{eqnarray} \tag{ 4 }$$

where ${\widetilde{T}}_{\max }=50\,\mathrm{Myr}$. This prior simplifies the marginalization in Equation (1) to

$$\begin{eqnarray}&&P(\tau )=\int P(\tau | S,{\widetilde{T}}_{\max })P(S)\,{dS}.\end{eqnarray} \tag{ 5 }$$

MATCH returns discreet samples of SFHs with star formation rates (SFRs) at discreet ages; hence, the age PDF will be discreet. In general, samples returned by MCMC algorithms are designed so that the density of the samples is proportional to the probability density function. Therefore, the probability of each SFH sample is

$$\begin{eqnarray}&&P(S)=\displaystyle \frac{1}{N},\end{eqnarray} \tag{ 6 }$$

where N is the number of samples. The discreet marginalized PDF for the age then becomes

$$\begin{eqnarray}&&P({\tau }_{i})=\displaystyle \frac{1}{N}\displaystyle \sum _{j=1}^{N}\displaystyle \frac{{S}_{j}({\tau }_{i})}{{M}_{j}({\widetilde{T}}_{\max })}.\end{eqnarray} \tag{ 7 }$$

The index i specifies the age bin, j specifies the MCMC sample, and ${M}_{j}({\widetilde{T}}_{\max })$ is the total stellar mass formed in the last ${\widetilde{T}}_{\max }$ yr. If the total mass formed is the same for each sample, then the discreet marginalized PDF is equal to the average SFH divided by the total stellar mass formed in the last ${\widetilde{T}}_{\max }$ yr.

In this manuscript, we report the most likely age and the narrowest 68% confidence interval (CI). Alternatively, one might report the 68% CI bounded by the 16th and 84th percentiles. However, this often-used technique is prone to several biases when applied to real data. If the PDF has exceptionally large wings, as real data often do, then the 16th and 84th percentiles tend to be centered around the median, not the mode. The age PDF in Figure 3 has a relatively narrow mode but very large wings. Therefore, if we simply report the 16th and 84th percentiles, then the CI will be strongly biased to the middle of the age range that we consider. Effectively, this CI will be centered around ${\widetilde{T}}_{\max }/2$, making the CI strongly dependent upon the maximum age chosen for normalization. Reporting the narrowest 68% CI mitigates these problems. The narrowest 68% CI is less sensitive to the large wings in the PDF, is much less dependent on ${\widetilde{T}}_{\max }$, and naturally centers on the mode of the distribution.

Zoom In Zoom Out Reset image size

N6946-BH1 offers a unique opportunity to compare the SED-fitting and age-dating methods.

From SED fitting, Adams et al. (2017) inferred a progenitor mass of ${25}_{-1}^{+2}$ M_⊙. With the more accurate TRGB-derived distance, the SED-derived mass becomes ${27}_{-1}^{+2}$ M_⊙ (the vertical black bar in Figure 3). Formally, they derived an uncertainty that is quite small, but this uncertainty in mass only accounts for the relatively small random uncertainty in the photometry of the progenitor. However, the last stages of stellar evolution are the most difficult to model, so the systematic uncertainties associated with modeling the SED are likely much larger. For example, Beasor & Davies (2016, 2018) found that a combination of significant mass loss, dust formation, and/or extreme bolometric correction in the last stages of red supergiants (RSGs) can dramatically reduce the flux in the most easily accessible visible and near-IR bands. Hence, the SED–to–progenitor mass mapping can be quite uncertain.

It is also not entirely clear that the progenitor would be in hydrostatic equilibrium just before collapse. The relative steadiness of the progenitor's luminosity (Adams et al. 2017) supports that the progenitor was likely in a hydrostatic state, but this is not entirely certain. Theoretical models suggest that there is a great deal of convective energy in the cores, enough to unbind the outer envelope if that energy can couple efficiently (Quataert & Shiode 2012; Fuller 2017).

The observed CMD is consistent with a population with an age of ${10.6}_{-5.9}^{+14.5}$ Myr. Within the assumption of single-star evolution, this age corresponds to an initial mass of ${17.9}_{-7.6}^{+29.9}$ M_⊙. By contrast, the SED-derived mass, ∼27 M_⊙, experiences core collapse at 7 Myr (Marigo et al. 2017; 7.6 Myr for the SED-derived mass, ∼24, from the rotating models). Formally, these ages and masses are consistent within our uncertainties. If the progenitor is truly associated with the burst at ∼10 Myr, then the most likely age either requires an ∼27 M_⊙ progenitor to live 40% longer than standard evolution theories predict or an ∼18 M_⊙ progenitor must have a luminosity 2.5 times brighter than predicted just before collapse. This discrepancy could be explained by either binary evolution or episodic eruptions at the end of the progenitor's life.

Unfortunately, in the central 100 pc, there are only a handful of bright MS stars to constrain the age of the stellar population. As a result, while the best-fit SFH does not show any populations younger than 10 Myr, 91% of the age PDF (Figure 3) has ages older than 7 Myr, making the brightness of the progenitor marginally consistent with the age PDF of the surrounding stars. In this case, standard single-star evolutionary models may be appropriate in describing the fate of the progenitor.

In general, theory predicts that the explodability of a star depends upon these parameters: the neutron star mass (M_NS), neutron star radius (R_NS), neutrino luminosity (L_ν), and amount of mass accreting during collapse $(\dot{{ \mathcal M }})$ (Murphy & Dolence 2017). In the context of these important parameters, there is a critical hypersurface that divides nonexploding and explosive solutions. During collapse, these parameters evolve over time, and if the evolution crosses the critical threshold, then the star explodes. The progenitor's structure prior to collapse likely determines the evolution in these parameters and whether the evolution crosses the critical hypersurface for explosion. To date, there is no direct mapping between the progenitor structure and whether the evolution crosses the critical hypersurface. However, there are many studies exploring more approximate explodability conditions (Burrows & Goshy 1993; Janka 2001; Murphy & Burrows 2008; O'Connor & Ott 2011; Ugliano et al. 2012; Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016).

In one such study, Sukhbold et al. (2016) explored the explodability of progenitors from 9 to 120 M_⊙; they found a general trend that the least massive stars explode more easily, but they also found that the explodability may not be entirely monotonic with mass. More specifically, they found islands of SN production. Generically, their simulations indicate that all masses below 15 M_⊙ explode. Above this mass, there are islands of failed explosions and BH formation. The existence of these islands seems to be robust to the details of stellar evolution assumptions, but the size and frequency of these islands depends upon the progenitor prescription. For typical stellar evolution assumptions, the region between 15 and 22 M_⊙ shows both explosions and failed explosions. Between 22 and 25 M_⊙, there seems to be a robust region of failed SNe. Above 28 M_⊙, most models fail to explode. Considering all results, the minimum progenitor to exhibit a failed SN is 15 M_⊙. This prediction is consistent with our best-fit progenitor mass for N6946-BH1, ${17.9}_{-7.6}^{+29.9}$ M_⊙.

If the explodability mostly depends upon the progenitor's initial mass, then constraining the progenitor distribution of SNe could constrain the progenitor distribution of failed SNe. For example, Smartt (2015) constrained the progenitor distribution for 18 SNe IIP (plus 27 upper limits on progenitors). Assuming a Salpeter mass function, they found a minimum mass of ${9.5}_{-2}^{+0.5}$ M_⊙ and a maximum mass of ${16.5}_{-2.5}^{2.5}$ M_⊙. This maximum mass for SNe IIP is consistent with the theoretical minimum mass for failed SNe. It is also consistent with our best-fit mass for N6946-BH1. However, one should note that the maximum mass for SNe IIP may not be the maximum mass for explosion. Progenitors above this mass may explode as other types of SNe. Complicating this interpretation, Davies & Beasor (2018) recently reanalyzed the progenitor masses of the SNe IIP using bolometric corrections that are derived from very red supergiants near death. With these new bolometric corrections, they found a higher upper mass (${19.0}_{-1.3}^{+2.5}$ M_⊙) for the progenitor masses of SNe IIP.

The distribution of masses from the age-dating technique does not show a maximum mass for explosion. However, the distributions are consistent with the most massive stars failing to explode more often than the least massive stars. Williams et al. (2018) age-dated the stellar populations for 25 SNe and found a distribution of median masses that is consistent with a wide range of masses. However, all 25 SNe have uncertainties that extend below 18 M_⊙, so the masses are also consistent with lower masses exploding and higher masses failing to explode. Díaz-Rodríguez et al. (2018) inferred the progenitor mass distribution for 94 S/Ns. They found a very well-constrained minimum mass (${7.33}_{-0.16}^{+0.02}$ M_⊙) and a power-law slope that is steeper than Salpeter ($-{2.96}_{-0.25}^{+0.45}$) but a maximum mass that is >59 M_⊙. While the maximum mass is well above the best-fit mass for N6946-BH1, the steeper-than-Salpeter (−2.35) slope could suggest that the most massive stars explode less frequently.

Both observations and theory of SN progenitors suggest that the least massive stars tend to explode more readily than the most massive stars. This implies that the most massive stars are more likely to form BHs by a failed SN. Our best-fit mass of N6946-BH1 is consistent with this general trend.