Model Light Curves for Type Ib and Ic Supernovae

The presupernova evolution of mass-losing helium stars was studied by Woosley (2019) for star masses, at helium ignition, from 1.6 $\mathrm{M}_{\odot}$ to 120 $\mathrm{M}_{\odot}$. This range corresponds to main sequence masses from approximately 10 to 250 $\mathrm{M}_{\odot}$. For standard assumptions regarding mass loss (Yoon, 2017, with $f_{\rm WR}=1$), helium stars with initial masses below 2.5 $\mathrm{M}_{\odot}$ were found to produce white dwarfs or electron-capture supernova and their final evolution was not followed. For the other stars, mass loss resulted in presupernova masses in the range 2 to 60 $\mathrm{M}_{\odot}$. The effect of a larger mass loss rate, closer to that favored by Yoon, $f_{\rm WR}$ = 1.5, was also examined.

The explosion of a large subset of these models, those that did not end up as as pulsational-pair instability supernovae or collapse directly to black holes, was simulated in a one-dimensional (1D) neutrino-transport code, PHOTB, by Ertl et al. (2020). The collapse to a black hole was determined by following the accretion onto the proto-neutron star long enough to assure that no outgoing shock emerged. In some cases not considered here, a black hole was formed by fall back even though the outer part of the star was ejected. This reduced model set consisted of 133 helium stars with initial masses from 2.5 to 40 $\mathrm{M}_{\odot}$, corresponding to presupernova masses from 2.1 to 19.6 $\mathrm{M}_{\odot}$. The efficiency of the neutrino-powered explosions was calibrated to SN 1987A and the Crab. The ejected mass, fall-back mass, final kinetic energy, and remnant mass were determined for each star in a self-consistent way that depended upon the progenitor structure. Considerable attention was paid to limits on the amount of ⁵⁶Ni each explosion produced, and lower and upper bounds consistent with the (1D) neutrino-powered assumption and nucleosynthesis were established. An additional set of models with mass loss rates 1.5 times the standard Yoon values ($f_{\rm WR}$ = 1.5) was also considered. More recently, Woosley, Sukhbold, & Janka (2020) used the remnant mass distributions calculated by Ertl et al. (2020) to estimate the initial mass-function averaged birth function and mean masses for black holes and neutron stars in close binary systems and found good agreement with present observables.

Ertl et al. (2020) used a variety of central engines to simulate their explosions. Their standard model was W18, but four other cases were also considered that used different presupernova models for SN 1987A. The average outcomes for these such as explosion energy, mass cut, and ⁵⁶Ni synthesis did not vary greatly for the five different central engines. See Figs. 2 and 3 and Table 2 of Ertl et al. (2020) and Figs. 8, 9, 13, and 17 and Table 4 of Sukhbold et al. (2016). Since we will find (§2.3; Fig. 2) that bolometric light curves calculated using SEDONA and KEPLER agree quite well, we will, with one exception, use only the models from Ertl et al. (2020) with their W18 central engine for the present study. The exception is an 8.00 $\mathrm{M}_{\odot}$ helium star evolved with twice the standard mass loss rate ($f_{\rm WR}$ = 2) and exploded manually without the benefit of a neutrino-transport calculation. This was necessary to examine the effects of the composition on the outcome of a star with presupernova mass low enough to be common yet composed mostly of carbon and oxygen and not helium (see §3.2 and Yoon et al., 2019). The evolution of this star gave a presupernova mass of 3.63 $\mathrm{M}_{\odot}$ and an iron core of 1.48 $\mathrm{M}_{\odot}$, similar to a standard model ($f_{\rm WR}=1$) with initial mass 4.7 $\mathrm{M}_{\odot}$. An explosion with final kinetic energy $1.28\times 10^{51}$ erg was generated using a piston situated near the edge of the silicon core at 1.64 $\mathrm{M}_{\odot}$.

Except for this single model, the results of Ertl et al. (2020) for the bolometric properties using the other central engines remain valid and sufficient. In addition to the standard model set with $f_{\rm WR}=1$, we continued to carry the parallel set with $f_{\rm WR}=1.5$ which was also exploded by Ertl et al. (2020).

The supernova models used are summarized in Table 1 and Table 2. Except for the 8.0 $\mathrm{M}_{\odot}$ model mentioned above, quantities in the tables were adopted without modification from Ertl et al. (2020). “${\rm Ni+Tr}$” and “$3/4\times{M_{\rm NSE}}$” refer to a typical and a maximum ⁵⁶Ni mass and their meanings are discussed extensively in that publication. Here, in keeping with the goal of producing the brightest possible supernovae this model set allows, the “$3/4\times{M_{\rm NSE}}$” value is adopted here as standard. These values can reasonably be reduced in any light curve simulation, but not increased. $M_{\rm presn}$ and $M_{\rm ej}$ are the masses of the presupernova star and the ejecta, after any fallback is over, and $E_{\rm exp}$ is its terminal kinetic energy. For the exceptional 8.00 $\mathrm{M}_{\odot}$ model that was manually exploded, a ⁵⁶Ni mass of 0.099 $\mathrm{M}_{\odot}$ was enforced by adjusting the mass cut. The reason for this will become apparent in §3.2 when the colors of this model are compared with another (4.50 $\mathrm{M}_{\odot}$; $f_{\rm WR}$ = 1) model with this mass of ⁵⁶Ni.

The standard models that used $f_{\rm WR}$ = 1 are labeled by their initial helium star mass and the prefix “He”. Model He5.00 thus began its computational lifetime as a 5 $\mathrm{M}_{\odot}$ helium star, lost mass using $f_{\rm WR}=1$, and exploded with a final mass of 3.81 $\mathrm{M}_{\odot}$ using the W18 central engine. The other explosion models in Table 1 that used $f_{\rm WR}=1.5$ will have an “x1.5” attached to their name. Model He5.00x1.5 thus also began as a 5 $\mathrm{M}_{\odot}$ helium star, but ended with a mass of 3.43 $\mathrm{M}_{\odot}$ and was exploded using the W18 central engine. All of these models in Table 1 are candidates for normal Type Ib and Ic supernovae (§3). The files are available to those wanting to use these explosion models for their own studies.

An additional 4 low and high mass models using the same explosion and mass-loss assumptions, but not considered “normal” are given in Table 2. Model He8.00x2 had an initial mass of 8.00 $\mathrm{M}_{\odot}$ but twice the standard mass loss rate.

SEDONA (Kasen et al., 2006, 2008; Roth & Kasen, 2015) is a multi-dimensional implicit Monte Carlo code, used here in one dimension. Given the ejecta properties i.e., the density, composition, and velocity profile of the freely expanding supernova material, SEDONA does a detailed treatment of the gamma-ray transport and self-consistently determines the emergent broadband light curves and spectra. The temperature structure is solved assuming radiative equilibrium. Lines are treated in the expansion opacity formalism and assumed to be purely absorbing. Otherwise no free parameters need be adjusted in the radiative transfer calculation. It is assumed that the excited states of the various ions are in thermal equilibrium. Opacities include the effects of over 20 million iron-group lines and about 800,000 lines of hydrogen through calcium (http://kurucz.harvard.edu/vitabib.html; Kurucz, 1994, 1995, 2009). The ions included here were He, C, N, O, Ne, Mg, Si, Ar, Ca, Ti, Cr, Fe, Co, and Ni. The only energy source was the decay of ⁵⁶Ni and ⁵⁶Co. Each ion except helium carried 7 stages of ionization and from 300 to 2000 levels. The initial thermal energy from the KEPLER explosion model was included.

To make a link with the explosion simulations, KEPLER models were evolved to a coasting configuration, typically 1 day, but 2 days in some of the more massive models with slower expansion speeds, and mapped into SEDONA. The original zoning was preserved except that few bottom zones with coasting speeds less than 500 km s^-1 and surface zones moving faster than 30,000 km s^-1 were excised. The ⁵⁶Ni mass and velocity of the remaining zones was then adjusted slightly so as to preserve the final total ⁵⁶Ni and kinetic energy given by Ertl et al. (2020, see also Table 1 and Table 2). In a few cases where sensitivity to ⁵⁶Ni mass was studied, the ⁵⁶Ni mass was multiplied by a constant and the other mass fractions renormalized so as to keep the sum of the mass fractions equal to one. Typical spatial resolution varied from 200 zones for the models with small ejected mass to 1220 zones for the heaviest model. The models were then mixed as discussed in §2.3.

Most runs were initialized with 1 million particles and added 140,000 particles per time step. Particles that escaped were allowed to leave the grid. Typical runs took about 350 time steps and terminated after 100 days. Finer time steps were used early in the calculation. The typical time step near maximum light was 0.3 days. The typical maximum number of particles on the grid at that time was 15 million. Several runs were carried out with four times the number of particles in order to obtain well converged spectra. These showed that the broad-band light curves were well converged for the smaller value. Calculations were run in the LUX supercomputer at UCSC using 4 to 8 nodes (160 to 320 CPU). Run times were shorter for the lighter models with higher velocity and fewer zones, but ranged from several hours to 15 hours.

The principal output of SEDONA consists of a time series of spectra from the far UV to infrared (10¹⁴ to $2\times 10^{16}$ Hz; 100 to 20,000 Å). In order to generate multiband light curves, these spectra were then folded with different filter response functions (§3.2).

One of the greatest uncertainties in any modeling of the multi-color light curves or spectra of Type Ib and Ic supernovae is the extent to which the ejecta, especially ⁵⁶Ni and ⁴He, have been mixed. See Dessart et al. (2012), Dessart et al. (2015), and Yoon et al. (2019) for discussions. The mixing of iron-group elements to large velocity provides a source of opacity that blankets the blue and ultraviolet emission making the supernova redder at peak. The decay of ⁵⁶Ni provides heat that maintains ionization to a larger radius in a mixed model hence supporting a larger photosphere and cooler effective temperature. The early color evolution of the supernova is thus especially sensitive to mixing Yoon et al. (2019). If ⁵⁶Ni and ⁴He are mixed together, the decay of the former can aid in exciting lines of the latter (Lucy, 1991; Dessart et al., 2012; Woosley & Eastman, 1997), dramatically affecting the spectrum and even the classification of the supernova as Type Ib or Ic. Above some mass, the small residual helium may be ejected with such high velocity that, depending on mixing, no ⁵⁶Ni and ⁴He commingle. Hence no strong, non-thermal lines of helium will be present in the spectrum and it will be Type Ic. But the supernova can also be Type Ic if there is a thick buffer of carbon, oxygen, and heavier elements between the ⁵⁶Ni and even a thick shell of helium. When ⁵⁶Ni is mixed to larger radii, the light curve rises earlier and may decline earlier. The peak bolometric luminosity is thus affected, though the change is generally not large.

Despite its importance, mixing in Type Ib and Ic supernova has not received nearly as much attention in multi-dimensional simulation as Type IIp. Focus in the latter has frequently been on the Rayleigh-Taylor instability created as the outgoing shock decelerates in regions of increasing $\rho r^{3}$, where $\rho$ is the local density and $r$ the shock radius. Major mixing occurs when the helium core encounters the hydrogen envelope and is forced to decelerate. This kind of mixing is present in exploding helium stars as well, but lacking a hydrogen envelope, the consequences are not as dramatic.

A second kind of mixing occurs in all varieties of core-collapse supernovae and is driven by the essential asymmetry of the central engine. In neutrino-powered explosions, the “hot bubble” is unstable and expands at different rates for different angles. See for example, Figs 10 and 11 of Wongwathanarat et al. (2017). If a magnetar powers the explosion, similar anisotropies will exist (Chen et al., 2017). For explosions that impart their energy in a time short compared with the shock crossing time for the presupernova star, about a minute here, the final velocity of the ejected plumes depends on their initial speed and how much matter they interact with on the way out. There is usually an upper bound to their speed (Wongwathanarat et al., 2017).

As is often noted, this mixing is not microscopic. It does not lead to the homogenization of the composition. Rather material is ejected in plumes and clumps that extend to higher final velocity than some of the other ejecta that initially was at smaller radius. Mixing is thus intrinsically a multi-dimensional process that cannot be properly replicated in 1D. For a clear demonstration in nature, see the supernova remnant for Cas-A and its well studied “fast moving knots”. Whatever is done here in 1D will be a gross approximation.

Many artificial prescriptions exist for mixing in one-dimensional simulations of supernovae. Most have, at their root, early attempts to model mixing in SN 1987A, and involve moving a “boxcar average” through the supernova ejecta multiple times. A boxcar with an interval of specified mass, is moved outwards through the supernova model, zone by zone until either some maximum mass or the surface of the ejecta is reached. At each zone the overlying interval is completely mixed while conserving mass. The mixing operation is often performed through the supernova several times with a decreasing mixing interval so that the final curve for abundance as a function of mass is smooth. In our previous studies (e.g. Ertl et al., 2020), an initial boxcar width of 0.15 $M_{\rm ej}$ was employed , where $M_{\rm ej}$ is the mass ejected in the supernova (presupernova mass minus remnant mass). This mixing region was moved through all the ejecta three times and then a final fourth mixing was applied using an interval half as great, i.e., 0.075 $M_{\rm ej}$. Importantly, this mixing was continued to the surface of the star, so some ⁵⁶Ni was mixed to arbitrarily high velocity.

A similar prescription was employed by Dessart et al. (2015), but instead of mixing in mass, they mixed in velocity intervals. Specifically, for the study of Type Ib and Ic supernovae, they recommended a 1000 km s^-1 mixing interval as characteristic of weak mixing and 2000 km s^-1 as representative of strong mixing and explored the consequences. The Dessart et al mixing and Ertl et al mixing are shown for Model He6.00, which ejected 2.32 $\mathrm{M}_{\odot}$ (see Table 1) in Fig. 1. Note the tail extending to very high velocity for the Dessart major mixing case and the Ertl et al prescription. Since the abundance of ⁵⁶Ni in these high velocity zones exceeds the abundance of solar iron in the initial model, 0.0014, this will have a major effect on the spectrum.

While we could find no published three-dimensional studies of post-explosive mixing in completely stripped massive stars exploded with neutrino transport, the results should resemble closely what Wongwathanarat et al. (2017) calculated for a Type IIb supernova with a very low mass hydrogenic envelope. Their initial model had a helium and heavy element core of 4.4 $\mathrm{M}_{\odot}$ capped by 0.3 $\mathrm{M}_{\odot}$ of low density hydrogen envelope. The ejected mass was 3.0 $\mathrm{M}_{\odot}$ of core material plus the 0.3 $\mathrm{M}_{\odot}$ of envelope, and the explosion energy was $1.47\times 10^{51}$ erg. Though intended as a model for Cas-A, the helium core mass and structure of their Model W15-2-cw-IIb is very similar to Model He6.0 here (Table 1). Wongwathanarat et al. (2017) found that mixing in their ejecta was characterized by a maximum speed for ⁵⁶Ni of about 7000 km s^-1 with only a small amount, about 1%, being mixed out to 7000 - 9000 km s^-1. They point out that this is consistent with the maximum speed seen for ⁴⁴Ti in Cas-A of 6300 $\pm$ 1250 km s^-1 (Grefenstette et al. 2017). They also point out that this speed is consistent with dimensional arguments for the bulk speed $(2E/M_{\rm ej})^{1/2}\approx 7000$ km s^-1. Dessart et al. (2016) say that this bulk speed is characteristic of the Doppler velocity measured for He I (5875 Å) at maximum in supernovae of Type IIb and Ib. It thus seems quite doubtful that substantial ⁵⁶Ni will be mixed to speeds in excess of 10,000 km s^-1 in models like He6.00. Dessart’s high mixing and Ertl et al. standard mixing are thus unrealistic at high speeds where the ejected material has been accelerated by a shock going down a density gradient.

With admittedly large residual uncertainty, we adopt a mixing formalism where the bulk of the ⁵⁶Ni and other heavy elements is mixed out to a speed given by the bulk velocity. That is the mixing is characterized by a maximum speed of

where $E_{51}$ is the kinetic energy of the ejecta in units of 10⁵¹ erg. To account for a small amount of material at greater speeds, a small amount of matter is mixed out to a value 30% greater. This formula would not describe mixing in a red or blue supergiant and is only for stripped helium stars.

This mixing is accomplished by the usual boxcar averaging technique. Two zones, $n_{1}$ and $n_{2}$ are determined that have velocities equal to $v_{\rm mix}$ and 1.3 $v_{\rm mix}$. Two mixing masses (boxcar sizes) are taken to be $m_{1}=0.1M_{\rm ej}$ and $m_{2}=0.02M_{\rm ej}$. The boxcar $m_{1}$ is passed 7 times from the center of the star to zone $n_{1}$ and then the boxcar $m_{2}$ is passed 4 times from the center of the star to zone $n_{2}$. The results for three different choices for $v_{\rm mix}$ are shown for Model He6.00 in Fig. 1. The choice of 6400 km s^-1 corresponds to the characteristic value in eq. (1).

Calculations using SEDONA show that neither the bolometric luminosity nor $V$ magnitude at peak is very sensitive to the mixing prescription, but the $B$ and especially $U$ magnitudes are. Were the mixing assumed in Ertl et al. (2020) to be employed here, the colors of the supernova at peak, as characterized by ($B-V$), would be too red compared with observations.

The peak bolometric luminosities and times for the standard models are given in Appendix B. All used the new prescription for mixing (§2.3). Also given in the tables are rise times ($t_{-1/2}$), decline times ($t_{+1/2}$), and decline rates ($\Delta m_{15}$) for the bolometric luminosity and similar information for the filtered luminosities. From this information one can approximately reconstruct any model light curve near its peak.

In Fig. 3, the light curves of typical models, He4.10 (based on a Salpeter IMF average of supernovae for models He2.70 to He8.00) and He4.60 (based on average between He3.30 and He8.00), are compared with the Carnegie sample of stripped-envelope supernovae (Taddia et al., 2018). Also given are the light curves of the faintest (He2.80) and the most luminous models (He5.38). Other models would fill in the space between these two curves. The agreement with most of the observations, including the faintest, is encouraging. There is a substantial set though, perhaps as many as one-third of the observed events, for which the brightest models are just too faint. This is essentially the same dilemma noted by Ertl et al. (2020) restated with better calculations of radiation transport and new observations. Is a substantial fraction of the stripped supernovae called “normal” by the observers outside of the reach of standard models based on a neutrino-powered explosion and a radioactively illuminated light curve? For this model set, it seems so.

We also calculated the “quasi-bolometric” light curves of the same models for direct comparison with the tabulations of Prentice et al. (2019). This entailed integrating the SEDONA spectral histories assuming a complete measurement within the 4000 to 10000 Å wavelength range and no detection outside. Fig. 4 shows the comparison. When peak luminosity is plotted against rise time, the models rise too late and are too faint to explain the brightest events. Some of these bright events are Type Ic-BL supernovae and almost certainly involve other physics not included in this study. Two other events are unusual. SN 2013bb (off scale in Fig. 4) ejected a very large mass and is slowly moving. Its broad, faint light curve is better explained by a more massive model than the common events. The discussion of such events is deferred to §5. At the other other extreme, SN 2018ie is a faint, rapidly evolving, high velocity, high temperature event. This might have been a low mass explosion related to the fast blue transients discussed in §4.

A comparison that excludes these anomalous cases is given in the second panel of Fig. 4. Since more supernovae had their decay time measured than their rise time, there is slightly more data in this figure. The agreement is significantly improved. One outlier, SN 2013F, is a fast rising, very luminous supernova. The host galaxy extinction, $E(B-V)_{\rm host}=1.4\pm 0.2$ mag, is much larger than for any other supernova in the sample though and its colors are uncertain. Perhaps its luminosity might be correspondingly uncertain.

Prentice et al. (2019) also give median luminosities, time scales, and ejected masses for a larger set of several dozen supernovae with data extracted from several sources. All luminosities are only for the 4000 - 10000 Å range. Table 3 compares these median values to the IMF-averaged models of this paper for the two different assumptions about the progenitor mass range.

Now the agreement for all measured quantities is quite good. The median rise time is a day or so longer in the models and the decline time two days shorter. The sum, i.e., the light curve width agrees quite well. Even in the models, it is difficult to determine the exact time of peak to better than a day (more so in heavier models) because the luminosity varies so slowly there. Though the data is for a limited number of supernovae, the agreement is somewhat better than that of the same models with the Carnegie data set (Fig. 3, Taddia et al., 2018, see also Table 3). Perhaps it is only the statistics of small numbers and the definition of what constitutes a “normal” supernova. We see no obvious explanation for a discrepancy.

Many others have extracted average ejected masses, ⁵⁶Ni masses, and kinetic energies from surveys of Type Ib and Ic supernovae (Richardson et al., 2006; Drout et al., 2011; Cano, 2013; Taddia et al., 2015; Lyman et al., 2016; Prentice et al., 2016; Taddia et al., 2018). For a summary table of these other results see Taddia et al. (2018). Most are consistent with Prentice et al. (2019), though often with larger error bars. Caution should be exercised when using the ⁵⁶Ni masses derived using Arnett’s Rule, as some of these surveys did. These are usually overestimates (§A.2). Using numerical modeling, Taddia et al. (2018) concluded, for Type Ib and Ic respectively, that $M_{\rm ej}$ = $3.8\pm 2.1$ and $2.1\pm 1.0$ $\mathrm{M}_{\odot}$; $E_{\rm exp}$ = $1.4\pm 0.9$ and $1.2\pm 0.7\times 10^{51}$ erg; and $M_{\rm Ni}$ = $0.14\pm 0.09$ and $0.13\pm 0.04$ $\mathrm{M}_{\odot}$.

Table 4 compares the IMF-averaged characteristics of the two model sets that used different mass loss rates ($f_{\rm WR}$ = 1 and $f_{\rm WR}$ = 1.5; Table 1). Unfortunately, no low mass explosions with $M_{\rm preSN}$ less than 3.43 $\mathrm{M}_{\odot}$ were computed for $f_{\rm WR}=1.5$ by Ertl et al. (2020), so one can only compare the averages for a limited range of masses whose boundaries do not reflect the full range of Type Ib and Ic progenitors. Table 4 should thus not be compared with surveys that do not select against low mass (typically low luminosity) events.

The use of a greater mass loss rate does not cause any major differences in the bulk characteristics of the supernovae. Even though the models to be averaged in Table 4 were selected on the basis of a common presupernova mass, 3.42 $\mathrm{M}_{\odot}$ to 5.63 $\mathrm{M}_{\odot}$ for $f_{\rm WR}$ = 1 and 3.43 $\mathrm{M}_{\odot}$ to 5.64 $\mathrm{M}_{\odot}$ for $f_{\rm WR}$ = 1.5, the average presupernova mass still differs by about 0.2 $\mathrm{M}_{\odot}$ because the weighting factors in the averaging process depend on the estimated main sequence masses of the stars. The larger range of higher masses for $f_{\rm WR}$ = 1.5 result in a more massive average. Even so, the ⁵⁶Ni masses and explosion energies are virtually identical. Because of the greater emphasis on more massive explosions, the ejected masses are a bit larger and that results in a slight lengthening of the time scales and decrease in the peak luminosity. From this limited study, one may conclude that the bulk characteristics of the supernovae are insensitive to the mass loss rate used to produce a given range of presupernova masses. This insensitivity does not carry over to the colors and spectra, however (§3.2).

Absolute magnitudes were computed by passing the SEDONA spectral histories through appropriate filters and using the prescribed zero points for various magnitude systems. In order to accommodate both convention and the data of several sets of observers, three different magnitude systems were used: Vega for Johnson–Cousins $UBVRI$, AB for SDSS $ugriz$, and “natural” for CSP $uBgVri$. All filters, and also zero points for Vega magnitudes, are given at http://svo2.cab.inta-csic.es/svo/theory/fps3. The zero points for the AB system was taken to be 3631 Jy. For natural magnitudes, they were calculated according to Krisciunas et al. (2017). Offsets at peak for a typical model are given in Table 1 in the Appendix. Unless otherwise stated, $UBVRI$ magnitudes in this paper are on the Vega scale and $ugriz$ magnitudes are on the AB scale. Exceptions are when models are compared with data exclusively from the Carnegie group. Then natural magnitudes are used.

The photometry and spectrum of a typical model calculated with SEDONA is shown in Fig. 5. As noted previously, He4.50 is close to the median presupernova mass, 4.6 $\mathrm{M}_{\odot}$, that one calculates for an IMF-weighted distribution of helium cores in the mass range 3.3 to 8 using the equivalent ZAMS masses. The mean ejected mass, kinetic energy, and ⁵⁶Ni mass for the same model set (Table 3) are also close to those of Model He4.50 (Table 1). The spectrum in Fig. 5 has been color coded to show the wavelengths to which the $ugriz$ filters are most sensitive.

Fig. 6 shows that the bolometric light curve depends chiefly on the bulk properties of the supernova - ejected mass, kinetic energy, ⁵⁶Ni mass - and not much on the composition of the presupernova star. Models He4.50 with $f_{\rm WR}=1$ and He8.00x2 with $f_{\rm WR}=2$ end up ejecting similar masses with the kinetic energies (Table 1 and Table 2). They also eject, by design, the same masses of ⁵⁶Ni. It is then perhaps not too surprising that their luminosities at peak are very similar, 10^42.39 erg s^-1 at day 19.5 for He4.50 and 10^42.37 erg s^-1 at day 20.2 for He8.00x2. The spectrum and colors are different though, especially at early times (see also similar results by Yoon et al., 2019). The carbon-rich star (He8.00x2) is redder than the helium-rich star (He4.50). Even with large mass loss rates the surface is never devoid of helium (e.g. Woosley, 2019) due to the recession of the helium convective core in the presence of mass loss, but Model He8.00x2 has a lot less total helium than He4.50. The helium, carbon, and oxygen masses in the two explosions are 0.13 $\mathrm{M}_{\odot}$, 0.47 $\mathrm{M}_{\odot}$, and 0.84 $\mathrm{M}_{\odot}$ for He8.00x2, but 0.96 $\mathrm{M}_{\odot}$, 0.13 $\mathrm{M}_{\odot}$, and 0.41 $\mathrm{M}_{\odot}$ for He4.50. The little bit of helium in He8.00x2 is in the outer layers and is ejected with high velocity. Removing an electron from carbon or oxygen is easier than taking one from helium, so the matter remains ionized farther out for a longer time in Model He8.00x2. This sustains a photosphere with a larger radius and may be why He8.00x2 is redder at a given luminosity, especially at times so early that something like a photosphere still exists. Another factor that could potentially affect the color is the abundance of iron group elements at the photosphere, but the abundance of iron in the unmixed outer layers was the same (solar), and the mixing and iron-group synthesis was also the same in the two models. This is therefore not an important effect.

While the rest of the paper will focus on the results using the standard mass loss rate, $f_{\rm WR}$ = 1, one should keep in mind these possible mass-loss-rate dependent variations in color indices. The mass loss rate used in He8.00x2, $f_{\rm WR}$ = 2, is probably an upper bound to the actual mass loss ate (Yoon, 2017; Woosley, Sukhbold, & Janka, 2020) so the variations shown in Fig. 6 may also be upper bounds.

Not all wavebands are treated equally well in SEDONA. $U$ (and $u$) are probably the least accurate and will not be emphasized in this study. The $U$-band is more sensitive to iron-line blanketing compared to other bands and to the strength of the calcium H and K lines. That means the $U$-band is quite sensitive to the distribution of iron-group elements with velocity, hence to mixing. It is also sensitive to temperature and luminosity, since small changes in temperature can change the ionization state (e.g., the ratio of Fe III to Fe II) in some zones, which has a big effect on the line blanketing. Uncertainties due to non-LTE effects can also play a role. A small error in, say, the LTE computation of the ionization state may have a significant effect on the $U$-band while only a modest effect on on the $R$-band.

Uncertainties in the input atomic line data can also play a role. The calculations for the main models here used a large line list (§2.2). Calculations with a reduced set of 370,000 lines showed a modest increase in the ultraviolet brightness and blue magnitudes, but not much change for other filters.

IMF-averaged color characteristics are given in Table 5 for the mass range considered as typical for normal Type Ib and Ic supernovae. Compared with observations, the absolute magnitudes of our models in various filters agree quite well with the Carnegie tabulation for 34 Type Ib, Ic and IIb supernovae (Fig. 7). Since the data in Stritzinger et al. (2018a) and Stritzinger et al. (2020) are given in CSP natural magnitudes, model light curves were computed accordingly. The faintest observed curve in Fig. 7 is SN 2007Y, a brief Type Ib supernova (Stritzinger et al., 2009, 2018a; Taddia et al., 2018). This supernova is treated in greater detail in §3.6.1. Explaining faint events generally poses no problems since the ⁵⁶Ni masses used here are upper bounds (Table 1).

The distributions of peak magnitudes in $V$ and $R$ magnitudes are given for our model set in Fig. 8. These are averages calculated using a Salpeter IMF (Salpeter, 1955) applied on the zero age main sequence mass distribution as described in Ertl et al. (2020). The time of maximum is defined for each color, i.e., the brightest $V$ magnitude is evaluated at the time of $V$-band peak, not bolometric maximum. The brightest $V$ magnitude for any of the standard models is $-17.83$. The brightest $R$ magnitude is $-17.99$. This is for Model He5.38. Slightly brighter models exist for other choices of explosion energy and mass loss. For Model He5.25S, exploded with slightly more powerful S19.8 engine, the magnitudes at peak are $B=-17.48$, $V=-17.96$, and $R=-18.03$. For He5.00 with 1.5 times the standard mass loss (Table 1) and with W18 engine the peak magnitudes were $B=-17.52$ , $V=-17.97$, and $R=-18.00$. Interestingly, the peak luminosities reach a maximum value around He5.00 to He6.00 and decline for heavier masses, especially above He8.00 (Table 1). This is because the ⁵⁶Ni masses cease rising (Table 1), but the slower expansion results in later peak times, more ⁵⁶Ni decay before peak, and thus lower luminosities. Our distribution of peak brightnesses is not being influenced on the upper end by the neglect of more massive models, at least for the neutrino-transport models considered here.

These distributions agree reasonably well with observations for about 2/3 of the ordinary Type Ib and Ic supernovae given in Fig 19 of Drout et al. (2011) which have been corrected for host galaxy extinction. Fainter, briefer events could come from lower mass stars (§4) and fainter broader events from higher mass ones (§5). Fainter values for peak magnitude could also be obtained by acceptable downwards adjustments of the ⁵⁶Ni yield in any of our models. The distributions in Fig. 8 lack, however, the brighter events in the observed distribution of Drout et al. (2011), even after broad-lined Type Ic supernovae are removed and we see no obvious reason for the discrepancy. A comparison with the more recent data of Taddia et al. (2018, see their Fig 8 and Fig. 3 here) shows somewhat better agreement and some support for a cut off around $V=-18$.

Fig. 9 shows the evolution of several frequently used color indices for the standard model set. The first three days of the explosion are not plotted because the KEPLER model is linked to SEDONA on day one, after homologous expansion has become a good approximation, and it takes another day or so for the radiation transport to adjust to the new code. The evolution near and shortly after peak has been best studied observationally and there is considerable interest in “pinches” in the various color indices about 10 days post-maximum (Drout et al., 2011; Taddia et al., 2015; Stritzinger et al., 2018b) since they might be used to estimate reddening.

Model results for ($B-V$), ($g-r$), ($g-i$), and ($V-r$) are plotted in CSP natural magnitudes and compared, at 10 days post-peak of $V$ or $g$, with measurements by Stritzinger et al. (2018b). Model results for ($V-R$) and ($V-I$) are given in Johnson-Cousins Vega magnitudes and are compared, for ($V-R$), with data from Drout et al. (2011). Approximate offsets for other magnitude systems are given in Table 1 and are less than 0.1 magnitudes. All data has been de-reddened and reflects the intrinsic values at the host. Model results for ($V-I$) have previously been published by Dessart et al. (2015). Overall the agreement between observations and the models at 10 days post-peak in Fig. 9 is reasonably good, though there is a tendency for the models to be a little too red at early times and far too blue at very late times. Compare our model curves with, for example, the intrinsic color templates in Fig. 9 of Stritzinger et al. (2018b).

There are several reasons for the disparity. One is an inherent uncertainty in the time of maximum light in the various bands. The visual brightness is slowly varying in the models near maximum with a typical full width at 97% maximum of about 6 days. Small changes in the spectrum would affect our estimated time of visual maximum appreciably and the observed dates of maximum brightness might also not be precise. At the same time, the color indices are rapidly evolving. An offset of 3 days would significantly affect the comparison between the models and observations. The model colors are also intrinsically sensitive to mixing, ⁵⁶Ni production (Fig. 10), the LTE approximation used in our studies, and the geometry of the explosion. Less line blanketing in the outer layers would make the model bluer near peak. Greater ionization in the outer layers, possibly due to non-LTE effects at late times, would make them redder.

Larger values for the color indices in Fig. 9 are indicative of cooler emission, so a positive slope indicates a time when the supernova is becoming redder. Abrupt changes of slope in a color index reflect evolving physical circumstances in the supernova. During the first several days (not shown), the color evolves sharply to the red as shock deposited energy diffuses out of the expanding outer layers. The luminosity declines or is nearly constant while the photospheric radius increases. Helium recombines at about 10000 K and a recombination front starts to move inwards in mass while still moving outwards in radius. This brief phase was not accurately captured in the SEDONA calculations and is not plotted, but is clearly seen in the KEPLER light curves and effective temperatures.

Once the outward moving diffusion front carrying energy from radioactive decay reaches the photosphere, the luminosity rises rapidly while the radius continues a slow increase. The combination causes the color to become blue and this is the cause of the first downward sloping part seen near the origin in some of the color index plots.

Shortly before peak, the luminosity rises more slowly while the radius rapidly increases. The color thus begins a long gradual ascent into the red that is prominent in all the panels of Fig. 9. By maximum luminosity roughly half the supernova has recombined. The photospheric radius is still increasing though and, as the luminosity declines, the reddening of the color indices continues. Other effects also come into play and the idea of a well-defined photosphere with a color governed by recombination like in a Type IIp supernova is increasingly overly simplistic. The opacity is not just due to electrons, but to lines that are sensitive to the ionization state and generally more important at lower temperature. The deeper the photosphere, the more heavy elements lie outside it, especially in mixed models, and this line blanketing also drives the color to the red.

About 10 to 20 days after maximum luminosity, with the longer time scale being appropriate to the more massive cases with longer diffusion times and larger fractions of carbon and oxygen, recombination reaches the center of the explosion and the photospheric phase ends. The concept of a photosphere is only an approximation anyway since different wavelengths originate from different depths and recombination is never complete, but from this point on, the supernova is increasingly just a cloud of ashes pumped by radiative decay. The power deposited by this decay (that which is not escaping) balances the luminosity. This is the onset of the nebular stage, though the full transition will take some time. Colors after this point vary according to the lines present in the filter, the mass of ⁵⁶Ni, and the density. The LTE model becomes an increasingly poor approximation. Among other things, it underestimates the ionization in a radioactively energized plasma. Underestimating the ionization means one sees to a deeper depth and hotter matter. For these reasons, the colors of our models are increasingly unreliable more than 20 days after peak luminosity.

As previously noted (Drout et al., 2011; Dessart et al., 2015; Stritzinger et al., 2018b), a “pinch” is seen in the distribution of colors measured in many wavebands about 10 days after peak. This reflects the regularity of the model set selected as Type Ib and Ic candidates and basic physics. The ratio of explosion energy to ejected mass and the ⁵⁶Ni mass do not vary greatly in the standard set (Table 1). One day after the explosion, all the models, except a few low mass ones with presupernova radius expansion, have a similar distribution of internal temperature (Fig. 11). The entropy in most of the ejecta at this time is about half due to radiation; the other half is electrons and ions. For adiabatic expansion, temperature thus declines more slowly than (radius)^-1, and hence slower that 1/t. This common thermodynamic configuration means that the light curves for the models without radius expansion will be quite similar during the first week. Nevertheless, the high density cores of all models would be fully recombined by 10 days after the explosion, well before peak, were it not for radioactive decay.

Helium recombines at $\sim$10,000 K. Carbon and oxygen combine with their last electrons at a lower temperature, typically $\sim$5000 K. For a blackbody, these two temperatures correspond to ($B-V$) = 0 and 1 which is very roughly the range of ($B-V$) in Fig. 9. These are approximate values because of the variable density, the non-uniformity of the temperature, and gradients in the composition after mixing, but for temperatures cooler than about 5000 K, the supernova is effectively recombined. Prior to recombination, the photosphere continues to move out in radius, forcing the color to the red. This is the epoch corresponding to the well-defined upward slopes in Fig. 9. The lower mass models have faster expansion speeds in their deep interiors and greater helium fractions. They experience recombination earlier, so their color-index curves (e.g., for $B-V$) are steeper. These differences in slope result in a convergence, a “pinching” of the distribution as time passes. The convergence is truncated when the supernova recombines. After that the sudden collapse of the photosphere sends things to the blue, though what happens later is not well determined in the present study.

As Fig. 9 shows, there is observational evidence for these pinches. Drout et al. (2011) have found a convergence of ($V-R$) evaluated 10 days after $V$-band maximum. For their “gold” set of well-studied Type Ib and Ic supernovae, ($V-R$)${}_{\rm 10\ days}=0.26\pm 0.06$ mag. In a study similar to ours, Dessart et al. (2016) found a similar convergence, but with a value $0.33\pm 0.035$ magnitudes, if they removed two high energy outliers. Stritzinger et al. (2018b) also observed similar pinches in a number of different filters, especially ($V-r$). These are shown as vertical error bars on the plot and agree reasonably well with the models.

It is encouraging that four studies, two observational and two theoretical agree on ($V-R$) at 10 days. As Dessart et al. (2016) note, weaker mixing favors higher temperatures in the inner ejecta, causing redder colors early on, but bluer colors around maximum (see also Dessart et al., 2012), so this color index might be useful for constraining mixing. As discussed by Drout et al. (2011) and Stritzinger et al. (2018b), the narrow range of several color indices in Fig. 9 might also be useful for better estimating the host galaxy extinction.

Fig. 12 shows the time-dependent color indices for the Carnegie collection of of 34 stripped supernovae (Stritzinger et al., 2018a). Here both data and models have been displayed in the “natural” magnitudes in which they were measured. The agreement for ($B-V$) is quite good. Both the average supernova models (between He4.10 and He4.50) and the band defined by the brightest and faintest models show near congruence with the observations. Interestingly the narrowing in ($B-V$) for the models at 10 days post-peak is more pronounced than in the observations. Typical error bars on the observations are 0.03 magnitudes.

The agreement is not as good for ($V-r$) and ($V-i$) showing again that our models emit inadequate radiation above about 5500 Å at late times. See also the spectroscopic comparisons to individual supernovae in §3.5. This could be due to the LTE approximation used in the models. Inadequate blue radiation is being processed into the red. The photosphere is collapsing and becoming quasi-nebular prematurely.

Fig. 13 shows the spectra for our models. Similar to Shivvers et al. (2019), the spectra within 10 days of peak are rather uniform, similar to the uniformity in colors in Fig. 9. This is especially true if the very light and very heavy models, those with initial masses over 7 $\mathrm{M}_{\odot}$ and below 2.9 $\mathrm{M}_{\odot}$ are excluded from the sample. The lighter models have experienced radius expansion and are unusual because of it. The heavier models are fainter, rarer, and may be unusually red. The differences that exist in the models are mostly in the ultraviolet short of 4000 Å where there is a deficiency in emission, perhaps due to line blanketing. This same deficiency is also seen in the similar (but non-LTE) models of Dessart et al. (2015). The 5880 Å line of He I here is not particularly strong, even 10 days post-maximum, suggesting that all of these supernovae could be identified as Type Ic. Helium is certainly present in the models, ranging from 83% of the ejected mass in He2.90 to 21% in He8.00. Probably the weak line reflects the use of the LTE approximation which does a poor job of reproducing the level populations in a mixture of ⁴He irradiated by the radiation from ⁵⁶Ni decay (Lucy, 1991; Dessart et al., 2012; Woosley & Eastman, 1997). Other line identifications can be made by comparison with Fig. 9 of Shivvers et al. (2019) and Fig. 14 of Dessart et al. (2015). The O I 7740 Å line is apparent, but not strong, and the 6200 Å line may be due to Si II. No hydrogen was included in the present study.

To illustrate the strengths and shortcomings of our models and to gain a better perspective on their differences in peak luminosity, five models were compared with five well-analyzed Type Ib and Ic supernovae. The events chosen were selected by observers (see “Acknowledgments”) to be representative cases. No effort was made to seek out the observations with which we agreed best. The five cases span a range of luminosities from very faint (SN 2007Y) to quite bright (SN 2009jf) and low mass to high mass. They include two events with double peaks (LSQ13abf and SN 2008D) and one low mass Type Ic event (SN 2007gr).

Spectra from the UV to the near IR for the 100 days following explosion were calculated for all models using SEDONA and photometry was obtained by numerically passing those spectra through various filters. As discussed in §3.4, the supernova undergoes a major adjustment about 20 to 30 days after the explosion (Fig. 9) when it ceases to be substantially ionized. In reality, gamma-rays from radioactive decay may keep the gas at least partly ionized longer than calculated in SEDONA, which assumes thermal equilibrium. This may have the effect of maintaining something resembling a photosphere for a longer time and make the supernova generally redder than calculated here. We thus trust the SEDONA calculations of spectrum and photometry only to about 20 - 30 days past peak, or about 50 days after the explosion.