Numerically Modeling the First Peak of the Type IIb SN 2016gkg

We first consider fitting envelope models with a convective density profile to SN 2016gkg. We include 20 different radii with R_e = 12.6, 15.8, 20.0, 25.1, 31.6, 39.8, 50.1, 63.1, 79.4, 100.0, 125.9, 158.5, 199.5, 251.2, 316.2, 398.1, 501.2, 631.0, 794.4, and $1000.0\,{R}_{\odot }$, and 15 different density scalings with ${K}_{c}=4.0\times {10}^{9}$, $6.3\times {10}^{9}$, $1.0\times {10}^{10}$, $1.6\times {10}^{10}$, $2.5\times {10}^{10}$, $4.0\times {10}^{10}$, $6.3\times {10}^{10}$, $1.0\times {10}^{11}$, $1.6\times {10}^{11}$, $2.5\times {10}^{11}$, $4.0\times {10}^{11}$, $6.3\times {10}^{11}$, $1.0\times {10}^{12}$, $1.6\times {10}^{12}$, and $2.5\times {10}^{12}\,{\rm{g}}\,{\mathrm{cm}}^{-3/2}$. These are exploded at eight different energies of ${E}_{51}=0.6$, 0.8, 1.0, 1.3, 1.6, 2.0, 2.3, and 2.6. Thus, in total we ran 2400 convective envelope explosion models.

The models were compared to the data only over the first $3.5\,\mathrm{days}$ so that we could focus on the first peak of the SN. In this way we are not sensitive to the details of the amount, location, or mixing of ⁵⁶Ni, which would influence the rise to the second peak. This comparison was evaluated with a simple ${\chi }^{2}$ calculation, where we take

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}{({M}_{i,\mathrm{observed}}-{M}_{i,\mathrm{model}})}^{2}/{\rm{\Delta }}{M}_{i,\mathrm{observed}}^{2},\end{eqnarray} \tag{ 6 }$$

where ${M}_{i,\mathrm{observed}}$ and ${M}_{i,\mathrm{model}}$ are the observed and model absolute magnitudes, respectively, ${\rm{\Delta }}{M}_{i,\mathrm{observed}}$ is the observed magnitude error, and the index i runs over all of the data points in all photometric bands. In addition, because of their constraining nature, we require the model to go through the first two data points by A. Buso and S. Otero. For each calculated model, we search through various potential explosion times using a bisectional algorithm until we find an explosion time that minimizes ${\chi }^{2}$. This is taken to be the ${\chi }^{2}$ for that particular model. Once the full grid of models is run, we can identify the minimum ${\chi }_{\min }^{2}$ for the entire grid. The next step is interpreting these values of ${\chi }_{\min }^{2}.$ The models will never be an exact fit to nature, and in addition there are systematic uncertainties in both the modeling and the data, as well as the fact that the grid spacing is not infinitely small. Therefore, we should not necessarily find ${\chi }_{\min }^{2}=1$. Our approach to this problem is to consider ${\chi }_{\min }^{2}$ as the best we can do, and therefore effectively treat ${\chi }_{\min }^{2}=1$. From this we can then estimate $1\sigma$, $2\sigma$, and $3\sigma$ uncertainties as models with ${\chi }^{2}\lt 2{\chi }_{\min }^{2},$ ${\chi }^{2}\lt 5{\chi }_{\min }^{2},$ and ${\chi }^{2}\lt 10{\chi }_{\min }^{2},$ respectively.

The results of applying this procedure are shown for the particular explosion energy ${E}_{51}=1.3$ in Figure 2. There are 300 models considered across this panel. The best-fit model corresponds to the red region with a value ${R}_{e}=199.5\,{R}_{\odot }$ and ${K}_{c}=1.0\times {10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-3/2}$, so that ${M}_{c}\approx 0.02\,{M}_{\odot }$. However, we note that given the coarseness of our grid (see the values listed above) and uncertainties in the modeling, there still remains at least a 20% error in these quantities.

Zoom In Zoom Out Reset image size

Also plotted in Figure 2 are lines of constant M_c using Equation (4). From this we see that there are two degeneracies acting in Figure 2. The first is at roughly constant R_e. This is simply set by the maximum luminosity of the first peak. The second runs along at roughly constant M_c, which is set by the width of the first peak. Such degeneracies are expected from the scalings described in more detail by Piro (2015), but it is reassuring for our fitting routine here that they naturally appear. This lends some robustness to our derived parameters for the envelope material, even if there are uncertainties in the modeling in detail.

In Figure 3, we consider the ${\chi }^{2}$ contours across all of the considered energies. This shows that although the energy we considered in Figure 2 of ${E}_{51}=1.3$ gives the best fit overall, an energy of ${E}_{51}=1.0$ can give a similarly good fit with a slightly larger radius. One can also see that as the energy is varied, different degeneracies gain or lose strength. At low energy, most of the degeneracy is for fixed R_e because these models have difficulty matching the width of the first peak. As the energy increases, the width is better matched and a degeneracy at constant M_c grows stronger. Nevertheless, independent of these issues, this demonstrates how robustly M_c and R_e can be constrained from this modeling with values of $\sim 0.02\,{M}_{\odot }$ and $\approx 180\mbox{--}260\,{R}_{\odot }$, respectively, where we take the spread for the value of R_e as roughly our $3\sigma$ uncertainties. Distance is an additional uncertainty not included in this estimate. Since R_e scales are proportional to the luminosity at peak, for a distance D the scaling is roughly ${R}_{e}\propto {D}^{2}$.

Zoom In Zoom Out Reset image size

It should be noted that there is also a degeneracy between the explosion energy and the mass of the core, because the larger the core, the less energy there is available for the envelope. Therefore, one should really think about the energy provided by the shock to the envelope, which scales as (Nakar & Piro 2014)

$$\begin{eqnarray}&&{E}_{e}\approx 2\times {10}^{49}{E}_{51}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{3{M}_{\odot }}\right)}^{-0.7}{\left(\displaystyle \frac{{M}_{c}}{0.01{M}_{\odot }}\right)}^{0.7}\,\mathrm{erg}.\end{eqnarray} \tag{ 7 }$$

This is an important distinction because the early time evolution and the associated velocities will better reflect E_e rather than E₅₁. For our best-fit energy of ${E}_{51}=1.3$ and ${M}_{\mathrm{core}}=3.55\,{M}_{\odot }$ (once the neutron star mass is subtracted off), we estimate from Equation (7) that ${E}_{e}\approx 3.8\times {10}^{49}\,\mathrm{erg}$. Thus, a larger or smaller E₅₁ would be expected for a smaller or larger ${M}_{\mathrm{core}}$, respectively, to keep this value of E_e roughly fixed.

The overall best-fit model is shown in comparison to the multi-band data in Figure 4. In comparison to previous semi-analytic fits to the data, the numerical models do a better job of following the changes from short to long wavelengths. In particular, the model by Piro (2015) predicts a much more symmetrical peak, which does an adequate job of representing the data at optical wavelength, but becomes increasingly worse at shorter wavelengths. This is because the accelerating shock velocity in the decreasing density near the surface of the envelope causes a stronger temperature evolution than predicted in the one-zone model by Piro (2015), which does not follow this velocity gradient.

Zoom In Zoom Out Reset image size

Furthermore, this comparison demonstrates how crucial the early data from A. Buso and S. Otero are for constraining the rise of the first peak and thus the envelope model. In fact, we infer from our modeling that the explosion must have occurred within ≈2–3 hr of this first data point! This is very close to the moment of explosion; for example, SN 2011fe was also observed very early as well and this was at roughly 4 hr (Bloom et al. 2012). As wide-field, transient surveys grow in the future, this work demonstrates the powerful constraints that we will be able to provide once more early time data is available.

Another aspect to note is that the best-fit extended mass of $0.02\,{M}_{\odot }$ is much less than the values around $\sim 0.1\,{M}_{\odot }$ that are typically presented by Woosley et al. (1994) and Bersten et al. (2012) for SNe IIb. This is because the first peak is most sensitive to only the mass near the maximum radius and not the total amount of hydrogen present (see the more detailed discussion in Nakar & Piro 2014, and in particular their Figure 2, which explicitly shows how the mass measured by the first peak compares with the total hydrogen shell mass). So while the total hydrogen mass can indeed be $\sim 0.1\,{M}_{\odot }$ to produce a realistic, hydrostatic model of a convective envelope, the first peak itself only can be utilized to measure the outer $0.02\,{M}_{\odot }$ of material.

We next consider a wind-like density profile for the envelope material. Traditionally SN IIb progenitors are thought to have extended, convective envelopes due to a recent mass transfer event that stripped the majority of its hydrogen-rich envelope (e.g., Benvenuto et al. 2013; Yoon et al. 2017). Nevertheless, if a star is in the midst of a mass transfer event, and the mass loss rate is great enough, the progenitor could in principle look like a yellow supergiant just from the optically thick wind. This is not traditionally considered for SNe IIb, but we explore such a density profile here in case assumptions about the density profile introduce uncertainties in the parameter estimations.

For our grid of wind models we consider the same 20 values for R_e and 8 values for E₅₁ as for the convective models discussed in Section 3.1. For the mass loading factor, we again consider 15 values with ${K}_{w}=6.3\times {10}^{15}$, $1.0\times {10}^{16}$, $1.6\times {10}^{16}$, $2.5\times {10}^{16}$, $4.0\times {10}^{16}$, $6.3\times {10}^{16}$, $1.0\times {10}^{17}$, $1.6\times {10}^{17}$, $2.5\times {10}^{17}$, $4.0\times {10}^{17}$, $6.3\times {10}^{17}$, $1.0\times {10}^{18}$, $1.6\times {10}^{18}$, $2.5\times {10}^{18}$, and $4.0\times {10}^{18}\,{\rm{g}}\,{\mathrm{cm}}^{-1}$. This again constitutes 2400 different wind models.

The comparison of the wind models with the data is summarized in Figure 5, plotted in the same way as before in Figure 3. The lines of constant M_w use Equation (5), and one can see that the degeneracy in M_e (especially at larger E₅₁) follow the slope of these lines. Overall, the best-fit model is at ${E}_{51}=1.6$, with ${R}_{e}=251.2\,{R}_{\odot }$ and ${K}_{w}=2.5\times {10}^{17}\,{\rm{g}}\,{\mathrm{cm}}^{-1}$, which corresponds to ${M}_{w}\approx 0.03\,{M}_{\odot }$. This is remarkably close to the result using a convective density profile, especially considering the coarseness of the grid. This demonstrates a strong constraint on these quantities, independent of the density profile.

Zoom In Zoom Out Reset image size

Even though the best-fit convective and wind density profile models have essentially the same radii and mass associated with them, the best-fit light curve for a wind profile is compared to the data in Figure 6. This shows that the wind profile gives a noticeably worse fit that the convective profile. The decline from the peak is just too steep in comparison to the wind profile model, perhaps arguing that such a profile cannot explain the data. In addition, it is important to ask whether such a wind model is even physically plausible. Using Equation (3) and an estimated wind velocity of $v\approx 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, the corresponding mass loss rate would be $\dot{M}\approx 0.8\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. This appears very high compared to what is normally expected for massive stars, although maybe it represents a star in the midst of mass transfer to a close binary companion, as is expected to take place for the progenitors of SNe IIb. Furthermore, the length of time implied by the radius of the extended material is ∼ weeks. With such a short timescale, it is implausible that the presence of enhanced wind generation should randomly occur so close to explosion. Either the two are casually linked, or perhaps the wind model does not really happen in nature. Future work should be done to better understand whether a wind environment as inferred here could actually be present in these binary scenarios.

Zoom In Zoom Out Reset image size

Besides having early photometry, SN 2016gkg was also exceptional because it had especially early spectra and velocity measurements in comparison to other SNe IIb. This provides a unique opportunity to use early velocities as another constraint on these models.

In Figure 7, we plot velocity data for SN 2016gkg from Tartaglia et al. (2017), which were computed by measuring the positions of the minima of the P-Cygni absorption components. This is not all of the lines that were measured, rather this is meant to be a representation of the fastest (${\rm{H}}\alpha$) and slowest (Ca ii H&K and Fe ii 5169) features observed during these early phases. Nominally, one would expect that these velocities represent an upper limit to the photospheric velocity, where the photospheric velocity roughly corresponds to where the continuum emission is mostly being generated.

Zoom In Zoom Out Reset image size

We also plot in Figure 7 the photospheric velocities inferred from our explosion modeling. These are defined as the velocity at the radius where the optical depth satisfies

$$\begin{eqnarray}&&\tau ={\int }_{r}^{\infty }\kappa \rho {dr}=2/3,\end{eqnarray} \tag{ 8 }$$

where κ is the specific opacity. Each color line in Figure 7 represents a best-fit model to the photometry for a given value of E₅₁. The specific value of ${E}_{51}=1.3$ is represented with a thicker line to highlight that this was our best-fit model overall. The numerical models generically show a large gradient in velocity as the photosphere transitions from the low-density extended material into the higher-density helium core. A comparison between the numerical results and the data show that none of the calculations are consistent with the observations. Even the lowest energy explosion that we consider overpredicts the velocities at the earliest times. The best-fit explosion energy does even worse. We make a similar comparison in Figure 8 for our wind models to emphasize that this problem cannot be reconciled by using a different density profile. At the earliest times, high velocities can cause the lines to become diluted. This can mean that the velocity measurements are more difficult to make, but it does not appear to explain the discrepancy we find here.

Zoom In Zoom Out Reset image size

How can these differences between the theory and observations be reconciled when the photometric fits seem reasonable? One resolution would be if the hydrogen-rich material of this SN is not distributed spherically symmetrically. In such a case, the material we are modeling to produce the first peak is indeed moving faster and is hotter than the material represented by the ${\rm{H}}\alpha$ absorption, but the two components have different angular distributions. This is interesting because it might provide an important clue about the nature of hydrogen-rich material surrounding the helium core, which is related to the progenitor scenarios. For example, one could expect that for a wind or strong mass transfer event the material would be denser (and thus able to sustain the shock generated by the SN) in some regions rather than others, depending on the binary configuration.

Alternatively, the asymmetries that we are inferring here may point to the generation of deeper asymmetries from the explosion itself. Light echoes from the SN that generated the Cas A remnant show that it was an SN IIb as is being studied here (Krause et al. 2008; Rest et al. 2008, 2011; Finn et al. 2016). Detailed studies of the distribution of material and kinematics in the remnant show that the explosion must have been very asymmetric, with some indication that it could have even had a jet-like component (Milisavljevic & Fesen 2013, 2015; Fesen & Milisavljevic 2016).

Cas A is additionally interesting because the region within its remnant does not show evidence of a massive ($\gtrsim 1\,{M}_{\odot }$) companion at the moment of explosion (Kochanek 2017). This is surprising both because models of SN IIb progenitors typically invoke a binary origin (e.g., Benvenuto et al. 2013; Yoon et al. 2017), but also because the great majority of massive stars are in binaries (Sana et al. 2012). One way to reconcile this is if Cas A was instead the result of a merger (a solution also suggested by Kochanek 2017), so that there was a companion, but it was lost before the SN. A merger might also explain the asymmetries that we infer for SN 2016gkg, as well as those seen in the remnant of Cas A. In the future, it will be important to continue looking for companions to nearby SNe IIb to get better statistics on how many were in binaries and how many may be due to mergers.

Finally, we note that Folatelli et al. (2014) also pointed out that some SNe IIb have strangely low velocities and discussed whether this can be produced by asymmetries. In that case the features in question are He i 5876 and He i 7065, and they are being measured well after the first peak and close to the second peak. Thus, that work is probing different material than the very shallow material that we are focusing on here. At larger depths, the low velocities are more natural, and may potentially be attributed to the high ionization energy of helium (Piro & Morozova 2014), but further work should be done to probe how deep asymmetries may be present in SN IIb ejecta.