Type IIP Supernova Progenitors. II. Stellar Mass and Obscuration by the Dust in the Circumstellar Medium

Gururaj A. Wagle; Alak Ray

doi:10.3847/1538-4357/ab5d2c

1. Introduction

Core-collapse supernovae (CCSNe) occur as a result of gravitational collapse of massive stars (typically, >8 M_☉ on the zero-age main sequence, ZAMS) at the end of their evolution. CCSNe are classified based on the observations of optical and infrared light curves and spectra. The eponymous plateau in SN IIP visible-band light curves typically lasts for 60–100 rest-frame days, after which exponentially decaying tails follow¹ (for details on classification of SNe, see the review article by Filippenko 1997). A volume-limited sample of CCSNe within 60 Mpc (Smith et al. 2011) shows that nearly 48% of these are supernovae of Type IIP, which have progenitors with a large hydrogen-rich envelope at the time of their explosion. Because of the observational constraints, the progenitors corresponding to supernovae can be identified only for relatively nearby cases (at d ≤ 30 Mpc). SN 2013ej, a Type IIP SN like its predecessor SN 2003gd, occurred in the same host galaxy M74 at a distance of ${9.0}_{-0.6}^{+0.4}$ Mpc (Dhungana et al. 2016) and was followed extensively by many groups in UV, optical, and infrared bands (Fraser et al. 2014; Richmond 2014; Valenti et al. 2014; Bose et al. 2015; Huang et al. 2015; Yuan et al. 2016) and in X-ray bands (Chakraborti et al. 2016). The mass of the progenitor of SN 2013ej estimated in the literature ranges between 11 and 16 M_☉, while the X-ray measurement of Chakraborti et al. (2016) points to a ZAMS mass of 13.7 M_☉ for it. Das & Ray (2017) restrict the ZAMS mass of the progenitor star of SN 2013ej using archival Hubble Space Telescope (HST) observations and their simulations to an upper bound of 14 M_☉. They also show through their 1D simulations of evolution of a 13 M_☉ star and its explosion that the observed light curves are fitted well with their simulations, when they included enhanced mass loss toward the end stages of evolution of the star. On the other hand, Utrobin & Chugai (2017) argue for a significantly higher progenitor mass of 25.5–29.5 M_☉ on the main sequence (and the ejecta mass of 23–26 M_☉) based on arguments of high-velocity ⁵⁶Ni ejecta. They also claim that a sufficient mass-loss rate can produce a circumstellar envelope at a distance about 10¹⁵ cm that would hide the pre-SN light from such a massive progenitor star.

We present here a few cases of an isolated (single), non-rotating, and non-magnetized star as a progenitor of SN 2013ej for two fixed ZAMS masses of 13 M_☉ and 26 M_☉. In the companion Paper I (Wagle et al. 2019), we have studied in detail the effect of convective overshoot on the evolution and explodability of the progenitor for a ZAMS mass of 13 M_☉. Here, we also test the case of a 26 M_☉ star as a possible progenitor of SN 2013ej. We study its luminosity variations in the late stages of stellar evolution and calculate its visual extinction due to dust formation in the circumstellar medium (CSM) formed by mass loss through stellar wind and compare these with archival observations of the progenitor star. We also demonstrate that our models have adequately fine mass resolution. In our future work, we will simulate the explosion of these models through 1D SN explosion codes and compare with the observed light curves and expansion velocity profiles to test the viability of the models.

In Section 2, we describe the methods of computational simulations and the stellar evolution using Modules for Experiments in Stellar Astrophysics (MESA). In Section 3, we discuss the results of variation of different MESA parameters on the stellar structure for our models. We also show that our models have adequate mass resolution. We discuss the extinction due to dust (formed in the CSM) on the observed magnitudes of the progenitor star. In Section 4, we discuss and summarize our conclusions.

2. Methods of Simulations

We use version r-10398 of the 1D stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018) to explore the evolution of the presumed progenitor of SN 2013ej from the ZAMS stage through to the core-collapse (CC) stage. All the MESA inlists are made available publicly at MESA marketplace.² In the following subsections, we discuss the choices of a few important parameters that affect the results discussed in this paper.

2.1. Initial Mass and Abundances

As discussed in the introduction, we investigate an isolated (single) star for two different cases of ZAMS masses of 13 M_☉ and 26 M_☉. The initial metallicity Z = 0.006 (same as in Paper I) is used in our simulations to create a pre-MS model. When provided with the Z value, MESA sets the initial helium abundance (Y) to a value equal to $0.24+2\times {Z}_{{\rm{initial}}}$ (refer to Equations (1)–(3) of Choi et al. 2016). As a result, for our simulations the initial H abundance (X) was set to 0.742 with Y = 0.252, since X + Y + Z = 1. MESA uses abundance values from Grevesse & Sauval (1998) to derive the initial abundances for each of the metals based on the choice of initial Z. The solar values for helium and metal abundances from Grevesse & Sauval (1998) are Y_☉ = 0.2485 and Z_☉ = 0.0169, respectively.

2.2. Mass and Temporal Resolution

In our simulations, we set the values for the mass and temporal resolution controls to those recommended³ in Farmer et al. (2016). We vary the mass resolution by restricting the maximum fraction of a star's mass in a cell, max_dq. This is achieved by setting a minimum mass resolution ${\rm{\Delta }}{M}_{\max }$ or dm, which is defined as max_dq = ${\rm{\Delta }}{M}_{\max }$ /M_*(τ) (Farmer et al. 2016). Adequate mass resolution is required to achieve successful convergence of stellar structure quantities between consecutive mass cells. However, a very high number of cells requires excessive computational resources. Farmer et al. (2016) recommend dm of 0.01 M_☉ for convergence of various quantities. We used this value in our models. We also tested an even finer resolution of 0.007 M_☉. As discussed later in Section 3, all our models with both these dm values are sufficiently resolved.

The time step between consecutive models in a simulation is generally controlled by varcontrol_target. This is achieved by modulating the magnitude of allowed changes in the stellar variables. We set it to its MESA default value of 10⁻⁴.

Farmer et al. (2016) note that the parameters such as delta_lg_XH_cntr_limit, delta_lg_XHe_cntr _limit, etc. offer another useful means of time-step control at critical stages of nuclear burning when one of the major nuclear fuels is depleted in the core, such as the terminal-age main sequence (TAMS). This set of limits ensure the convergence of mass shell locations, smoother transition in the Hertzsprung–Russell diagram (HRD) at the "Henyey Hook" (Kippenhahn et al. 2012), and smoother trajectories in the central temperature–density (T_c–ρ_c) plane. For a few of our model simulations we had to relax some of these limits at steps where we encountered convergence problems (refer to Table 1 for details). The parameter dX_nuc_drop_limit restricts the maximum allowed change in mass fractions between the time steps when the mass fraction is larger than dX_nuc_drop_min_limit. We set these values similar to those in Farmer et al. (2016).

Table 1. Inlist Parameters and MESA-predicted Core Properties for M_ZAMS = 13 M_☉ and 26 M_☉ Models with Z = 0.006 and f = 0.025

dm	Overshoot	delta_lg_X_cntr_limit/hardlimit			Max dt	Dutch Wind	Network	Model ID	${\mathrm{He}}_{\mathrm{core}}$	C_core	O_core	Si_core	Fe_core
(M_⊙)	Parameter, f₀	Ne	O	Si	Change	η	Size
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
M_ZAMS = 13 M_⊙

0.007	0.050	0.02/0.03	0.02/0.03	0.02/0.03	1.15	0.5	22 isotopes	13 M_⊙ model 1	4.271	2.108	2.056	1.576	1.477
0.007	0.005	0.02/0.03	0.02/0.03	0.02/0.03	1.15	0.5	22 isotopes	13 M_⊙ model 2	4.415	2.179	2.120	1.611	1.510
0.007	0.005	0.02/0.03	0.02/0.03	0.02/0.03	1.15	0.5	79 isotopes^a	13 M_⊙ model 3	4.402	2.187	1.898	1.545	1.481
0.01	0.050	0.02/0.03	0.02/0.03	0.02/0.03	1.15	1.0	22 isotopes	13 M_⊙ model 4	4.260	2.102	2.052	1.591	1.495
0.007	0.050	0.02/0.03	0.02/0.03	0.02/0.03	1.2/1.15	0.5	45/204 isotopes	$13\,{M}_{\odot }$ model 5^b	4.265	2.059	1.929	1.497	1.365
⋯	0.005	⋯	⋯	⋯	1.2	0.5	21 isotopes	13 M_⊙ modified test suite	4.381	2.147	1.969	1.593	1.468

M_ZAMS = 26 ${M}_{\odot }$

0.007	0.050	0.02/0.03	0.02/0.03	0.02/0.03	1.2	0.5	22 isotopes	26 M_⊙ model 1	8.406	4.610	1.889	0.000	1.529
0.01	0.050	0.015/0.03	0.015/0.03	0.015/0.03	1.2	0.5	22 isotopes	$26\,{M}_{\odot }$ model 2^c	6.876	3.262	2.948	1.729	1.589
0.01	0.020	0.01/0.02	0.01/0.02	0.01/0.02	1.2	0.5	22 isotopes	26 M_⊙ model 3	7.572	3.952	3.820	1.690	1.538
0.01	0.050	0.01/0.02	0.02/0.04	0.01/0.02	1.2	0.5	79 isotopes	$26\,{M}_{\odot }$ model 4^c	6.876	3.253	2.963	1.811	1.587
⋯	0.050	⋯	⋯	⋯	1.2	0.5	21 isotopes	26 M_⊙ modified test suite	9.209	5.294	2.427	1.660	1.508

Notes. delta_lg_X_cntr_limit and hardlimit values were set to 0.01 and 0.02, respectively, through CC for H, He, and C in all of the above models. See Section 2 for more details on the choices of other parameters listed here.

^aNetwork provided by Farmer et al. (2016) at MESA marketplace. ^bmax_timestep_factor was changed after TAMS from 1.2 to 1.15. The number of isotopes in the network was changed from 45 (mesa_45.net) to 203 (si_burn.net, provided by Renzo et al. 2017), and the maximum numbers of grid points allowed was changed to 5000 from the MESA default value of 8000, after O-depletion (see Renzo et al. 2017 for the definition of O-depletion and the explanation for these changes). ^c26 M_⊙ model 2 run was restarted at model number 7000 by changing delta_lg_X_cntr_limit values for Ne, O, and Si to 0.015 (and corresponding hard limit to 0.03) and model 4 run was restarted at model number 6000 by changing only delta_lg_XNe_cntr_limit to 0.02 (hard limit to 0.04) from the original corresponding limit values of 0.01 (hard limit of 0.02). The model was running into convergence problems with the original lg_XNe_cntr_limit values at those stages. In other models, the values were kept constant through CC.

Download table as: ASCII Typeset image

2.3. Opacities and Equation of State

We have used Type II opacity tables (Iglesias & Rogers 1996) that take into account varying C and O abundances during He-burning and later stages of evolution, beyond that accounted for by Z. MESA uses Type I tables (Cassisi et al. 2007) instead of Type II tables where metallicity is not significantly higher than Zbase. We set Zbase to the same value as the initial Z = 0.006. MESA uses the HELM equation of state to account for an important opacity enhancement during pair production because of the increasing number of electrons and positrons per baryon at late stages of nuclear burning.

MESA uses equation-of-state tables based on OPAL tables (Rogers & Nayfonov 2002), and at lower densities and temperatures SCVH tables (Saumon et al. 1995), with a smooth transition between these tables in the overlapping region defined by MESA.

2.4. Nuclear Reaction Rates and Networks

MESA mainly uses NACRE (Angulo et al. 1999) rates for thermonuclear reactions with updated triple-α (Fynbo et al. 2005) and ¹²C(α, γ)¹⁶O rates (Kunz et al. 2002) among others. (For more details see Section 4.4 of Paxton et al. 2011.)

In our model simulations, we used the softwired network consisting of 79 isotopes (mesa_79.net, Farmer et al. 2016) to optimize between the computational time and convergence of various values at the CC stage. For comparison, we present a few simulations with a smaller hardwired network consisting of 22 isotopes, approx21_cr60_plus_co56.net, and a single simulation with a softwired network consisting of 45 isotopes (mesa_45.net) up to O-depletion and a network consisting of 204 isotopes (si_burn.net, Renzo et al. 2017) thereafter until CC. The so-called "hardwired" networks have the predetermined pathway for each reaction, while the "softwired" networks link all allowed pathways between the isotopes specified in the network. In addition to these simulations, we also present a simulation for each of the two ZAMS masses where we used a slightly modified version of the "test suite" for a pre-CCSN star that accompanies MESA distribution. The test suite uses a network consisting of 21 isotopes, approx21_cr56.net.

2.5. Mixing, Diffusion, and Overshoot

MESA uses the standard mixing length theory (MLT) of convection (Cox & Giuli 1968, chap. 14). We used the Ledoux criterion in our simulations, instead of the Schwarzschild criterion to determine the convective boundaries (Paxton et al. 2013). We used ${\alpha }_{\mathrm{MLT}}\,=$ 2, which is an intermediate value between the values of 1.6 and 2.2 inferred from comparison of observations with stellar evolution models in the literature (see Section 2.3 of Farmer et al. 2016, and references therein). Here the mixing length equals ${\alpha }_{\mathrm{MLT}}$ times the local pressure scale height, ${\lambda }_{{\rm{P}}}=P/g\rho$ . We also enabled two time-dependent diffusive mixing processes, semiconvection and convective overshooting, in our simulations. Semiconvection refers to mixing in regions unstable to the Schwarzschild criterion but stable to Ledoux, ${{\rm{\nabla }}}_{\mathrm{ad}}\lt {{\rm{\nabla }}}_{{\rm{T}}}\lt {{\rm{\nabla }}}_{{\rm{L}}}$ . Semiconvection only applies when the Ledoux criterion is used in MESA. We set the dimensionless efficiency parameter (α_SC) to 0.1 for efficient mixing (same as in Das & Ray 2017). MESA's treatment of overshoot mixing is described in our Paper I in more detail. The parameter f₀ controls the degree of mixing in the overshoot region, whereas the parameter f determines the extent of the overshoot. We have only applied the overshoot for H-burning and non-burning cores and shells. For the simulations presented in this paper, we keep the value of the overshoot parameter, f, fixed at 0.025.

2.6. Mass Loss by Stellar Winds

For the "hot" phase of evolution we used the "Vink" scheme (Vink et al. 2001) with η_Vink = 1.0 (Vink_scaling_ factor), and for the "cool" phase we used the "Dutch" wind scheme for both the asymptotic giant branch and red giant branch phases, with η_Dutch = 0.5 (Dutch_scaling_factor) in the simulations presented here. This particular combination is used to confine the mass-loss rate to moderate levels during the red supergiant (RSG) phase. We note that the same combination was used by Das & Ray (2017).⁴ The particular combination used for the "Dutch" scheme in MESA is based on Glebbeek et al. (2009). The MESA mass-loss rates are switched between the rates of Vink et al. (2000, 2001) and that of de Jager et al. (1988) at temperatures below 10,000 K.

2.7. Stopping Criterion

The model is stopped very close to the collapse of the star when the infall velocity at any location in its interior reaches fe_core_infall_limit = 10⁸ $\mathrm{cm}\,{{\rm{s}}}^{-1}$ . We turned on the velocity variable "v" defined at cell boundaries, by setting v_flag=.true. when the central electron fraction Y_e (electrons per baryon, $\bar{Z}/\bar{A}$ ) drops below 0.47 (center_ye_limit_for_v_flag). This enables MESA to compute the hydrodynamic radial velocity, which is used for evaluating the stopping criterion.

3. Results

3.1. HRD and Mass Resolution

Figure 1 shows the evolutionary tracks in the HRD for (a) the 13 M_☉ star and (b) the 26 M_☉ star. Various 13 M_☉ models closely track each other throughout the evolution of the star, with minor differences stemming from differences in choices of parameters. Model 4 with η_Dutch = 1.0 remains less luminous than other models with η_Dutch = 0.5. The variations in the models with same η_Dutch, which are due to differences in the choices of one or two parameters (see Table 1), are not significant enough to be distinguished through observations.

On the other hand, we notice that the effects of choices of parameters are more prominent in the 26 M_☉ model, especially with dm, which controls the minimum level of mass resolution. We see in panel (b) for the HRD that models 2, 3, and 4, which use a coarser mass resolution (dm = 0.01 M_☉), exhibit spirals in their evolutionary tracks around core Ne-ignition. This feature is independent of the choices for other parameters. The greatest variation in luminosity takes place in model 2 (red) between 1.3 and 0.6 yr before collapse from log(L/L_☉) = 5.4 to 5.5 (see the inset in panel (b)). The corresponding variation in V- and I-band magnitudes reported by MESA is 0.3 mag, which can be easily observed. Unfortunately, observations of the progenitor of SN 2013ej are available only at about 8 and 10 yr before its collapse. The variation in magnitudes between 10 to 8 yr before collapse in the models above is consistent with the HST observations. Note, however, that the spirals are not exhibited in model 1, which uses a finer mass resolution (dm = 0.007 M_☉). The modified test suite models for both the ZAMS masses explored here are presented for comparison purposes only. The test suite models suffer from a lack of strict temporal and mass resolution, the effect of which is prominent in the inset of panel (b) for 26 M_☉.

According to Sukhbold et al. (2018) one of the criteria for adequate zoning is that the key variables such as density and temperature do not vary significantly between consecutive zones. In Figure 2, we plot density and pressure scale heights along with the mass resolution through the interior of the 13 M_☉ and 26 M_☉ stars. The 13 M_☉ model 3 with dm = 0.007 is shown in panel (a) and the 26 M_☉ models 2 and 3 with dm = 0.007 M_☉ and 0.01 M_☉ are shown in panels (b) and (c), respectively. A scale height for a quantity "x" is defined as M_x = (d ln x/dm)⁻¹, where the absolute value of M_x is considered. The large upward spikes in Figure 2 represent regions of nearly constant temperature and density. The discontinuities in mean molecular weight at the edge of convective shells result in quantities being artificially small. The edge of the helium convective shell (mass coordinate ∼4 M_☉ in the 13 M_☉ model and 8 M_☉ in the 26 M_☉ model) shows a steep gradient where the density changes by a few orders of magnitude, where the zoning is not well resolved. A finer resolution comes at a cost of increased computational time. Other than the edge of the He shell, the zoning is 2 dex finer than the scale heights everywhere inside both the stars, indicating that the quantities are over-resolved for both the dm values.

**Figure 2.** The density (gray) and temperature (red) scale heights along with the mass resolutions (zoning, blue dotted line) for the M_ZAMS = 13 M_☉ model (a) and 26 M_☉ models (b), (c) with Z = 0.006 at CC. The ratio of scale height M_x = (d ln x/dm)⁻¹ (Sukhbold et al. 2018) for a quantity "x" to mass resolution is an indicator of zoning efficiency. Panel (a) shows the 13 M_☉ model with dm = 0.007 M_☉, and panels (b) and (c) show 26 M_☉ models 1 and 2 with dm = 0.007 M_☉ and 0.01 M_☉, respectively. The fluctuations in density and pressure scale heights are over-resolved in our simulations, which is apparent from the 2 dex fine zoning in most places.
Download figure:
Standard image High-resolution image

3.2. Dust Extinction for a 26 M_⊙ Progenitor

Table 1 of Fraser et al. (2014) lists the observed magnitudes from archival pre-explosion images of the progenitor for SN 2013ej. The observations made in 2003 November and 2005 June are available, about 10 and 8 yr before the collapse of the star. We converted the HST-observed magnitudes in the F555 and F814 bands to Johnson's photometric magnitudes V and I using a formalism explained by Sirianni et al. (2005). To compare the MESA-calculated V- and I-band magnitudes to the observations, we evaluated the dust extinction due to the stellar wind formed by CSM, using the formalism provided in the Appendix. These dust extinction values are listed in Table 2 along with the HST-observed magnitudes and the MESA-calculated magnitudes for a 26 M_☉ fiducial star. We notice that the V-band magnitude is reduced by about 0.3 mag in ∼2 yr from 2003 to 2005. Our MESA calculations do not show such variations for any of our 13 M_☉ or 26 M_☉ simulations between 10 and 8 yr before collapse. The dust extinction values calculated based on the stellar wind model do not account for this variation either.

Table 2. Magnitudes and Dust Extinction for M_ZAMS = 26 M_☉, Z = 0.006 Model 1

HST	Johnson	Synthetic	MESA-	Min.	Calculated A(ν) for Graphite					Calculated
	Band	Observed	calculated	A(ν)	Minimum	Heat Balance		Adiabatic Cooling		A(ν) for
Magnitude		Magnitude^a	Magnitude	Required	Grain Size	${T}_{d,\max }$	${T}_{d,\max }$	${T}_{d,\max }$	${T}_{d,\max }$	Silicates^c
		( ${m}_{\mathrm{obs}}$ )	( ${M}_{\mathrm{MESA}}+{\mu }_{d}$ )^b		(Å)	1500 K	2000 K	1500 K	2000 K
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
${\tau }_{{\rm{CC}}}-\tau \approx 10$ yr, ${R}_{\mathrm{star}}=1017.73\,{R}_{\odot }$

25.16 ± 0.09	V	24.99 ± 0.13	21.91	3.08	800	1.35	3.89	4.51	5.28	0.74
					1500	1.13	3.27	3.80	4.44
22.66 ± 0.03	I	22.70 ± 0.09	20.40	2.30	800	0.65	1.86	2.16	2.53	0.36
					1500	0.54	1.57	1.82	2.13

${\tau }_{{\rm{CC}}}-\tau \approx 8$ yr, ${R}_{\mathrm{star}}=1017.69\,{R}_{\odot }$

24.84 ± 0.05	V	24.69 ± 0.08	21.91	2.78	800	1.35	3.90	4.51	5.32	0.74
					1500	1.14	3.28	3.79	4.47
22.66 ± 0.03	I	22.69 ± 0.08	20.40	2.29	800	0.65	1.87	2.16	2.55	0.36
					1500	0.54	1.57	1.82	2.14

Notes. The calculated magnitudes from our MESA simulation are listed along with the observed HST WFC magnitudes in F555W and F814W bands from Fraser et al. (2014). The dust extinction values calculated using the formalism explained in the Appendix are listed here.

^aThe synthetic magnitudes in the V and I bands are calculated from HST observations using the algorithm explained in Sirianni et al. (2005). ^bA distance modulus μ_d = ${29.77}_{-0.1}^{+0.09}$ ${29.77}_{-0.1}^{+0.09}$ is calculated considering the distance of M74 to be ${9}_{-0.6}^{+0.4}$ ${9}_{-0.6}^{+0.4}$ Mpc. ^cA single value of ${T}_{d,\max }$ ${T}_{d,\max }$ = 1500 K was used to find R_min for silicates, because silicates do not withstand a higher temperature.

Download table as: ASCII Typeset image

We see from Table 2 that a minimum A_V ≈ 2.8–3.1 is required to shield our fiducial 26 M_☉ progenitor. Although we see such high numbers among our calculated A_V values (columns (6)–(10) in the table), it is significantly higher than the value of ${0.45}_{-0.04}^{+0.04}$ calculated by Maund (2017) in the host for the progenitor of SN 2013ej through Bayesian analysis of the stellar population, using archival pre-SN data. As the host galaxy M74 is observed face-on, most of this extinction would be due to the CSM around the star. For another supernova of Type IIP, SN 2003gd in the same host galaxy M74 as SN 2013ej, Maund (2017) calculates A_V = ${0.46}_{-0.03}^{+0.04}$ .

3.3. Central Temperature versus Density

Figure 3 shows the evolution of central temperature and density from ZAMS through CC. The 13 M_☉ models shown in panel (a) are in good agreement with each other until about when C-burning starts in the core ( ${\mathrm{log}}_{10}{T}_{c}\,\approx$ 8.9 and ${\mathrm{log}}_{10}{\rho }_{c}\,\approx$ 5.5). The 26 M_☉ models shown in panel (b) diverge after the TAMS stage ( ${\mathrm{log}}_{10}{T}_{c}\,\approx$ 8 and ${\mathrm{log}}_{10}{\rho }_{c}\,\approx$ 2). We notice that models with smaller networks tend to be hotter and less dense at the center through various burning stages. The insets in both the panels show evolutionary tracks during core O-burning and core Si-burning stages. We note that the tracks in this part show rapid variations, more so during core Si-burning, indicating that these phases of nuclear burning are challenging due to very high and nearly balancing reaction rates, as also pointed out in Renzo et al. (2017).

3.4. The Kippenhahn Diagram

The Kippenhahn (KH) diagram is a useful representation of the succession of various convective burning zones inside the star throughout its evolution. We have discussed the KH diagram for a 13 M_☉ star in Paper I. In Figure 4 we compare the 26 M_☉ model 1 with dm = 0.007 M_☉ and f₀ = 0.050, model 3 with dm = 0.01 M_☉ and f₀ = 0.020, and model 4 with dm = 0.01 M_☉ and f₀ = 0.050. Models 1 and 3 both use a network consisting of 22 isotopes while model 4 uses one consisting of 79 isotopes. Note that the 26 M_☉ models 3 and 4 (Figures 4(b) and (c)), which both exhibit spirals in the HRD, have a distinct intermediate convective zone (ICZ) at the bottom of the convective hydrogen shell (at about 16 M_☉) starting at about a few thousand years before the collapse, although the ICZ is very limited in mass extent. This ICZ is missing in model 1 (Figure 4(a)). The mass resolution and the network size affect the internal structure of the star. The various core boundaries in models 1 and 4 with the coarser mass resolution (dm = 0.01 M_☉) are systematically lower than those in model 1 with the finer mass resolution (dm = 0.007 M_☉). The star is more massive at CC in model 1 than in models 3 and 4. The larger network in model 4 results in a more compact C–O core (yellow and green solid lines) than model 3. However, the final sizes of the Si and Fe cores remain somewhat less affected by differences in network size.

**Figure 4.** Kippenhahn (KH) diagram for the M_ZAMS = 26 M_☉ and Z = 0.006 models described in Table 1. Panel (a) shows model 1 with dm = 0.007 M_☉ and f₀ = 0.050, while panel (b) shows model 3 with dm = 0.01 M_☉ and f₀ = 0.020, both using a network consisting of 22 isotopes. Panel (c) shows model 4 with dm = 0.01 M_☉ and f₀ = 0.050, with a network consisting of 79 isotopes. Various core boundaries reported by MESA are marked in the diagram. By default, the core boundaries are determined by the point at which the mass fraction of the previous isotope decreases below a threshold value of 0.01 and where the same for the current isotope is above 0.1.
Download figure:
Standard image High-resolution image

3.5. Entropy from Oxygen-burning onward

The entropy profiles at various stages are shown in Figure 5, from ZAMS through C-ignition in panels (a) and (c) and C-ignition through CC in panels (b) and (d) for the two ZAMS masses. The profiles show sharp increases at the boundaries of various fossil shell-burning zones. Due to high core temperatures, the core is cooled predominantly by neutrinos emitted by thermal processes rather than by radiation from the carbon-burning stage onward. However, from oxygen-burning onward, the non-thermal emission of neutrinos as a result of electron capture on nuclei and the beta decay of others starts contributing. Neutrinos are efficient cooling agents because they escape from the core without any further interaction once produced. This causes entropy to decrease in the central parts of the core. At the edges of the core (or in regions of shell-burning), sharp entropy gradients develop that extend up to the region that is well mixed by convective transport and where the temperature is high. These sharp entropy gradients that prevent the mixing and penetration of burning products lead to the onion-skin structure in the core with distinct elements (e.g., C, Ne, O, Si) . Panel (a) shows that as the star's entropy in the core decreases, the entropy in the envelope increases with time, even though during the various early nuclear burning stages the entropy profile remains flat for a large part of the star except near its surface. Panel (c) shows a similar trend for the 26 M_☉ star.

3.6. Mass-loss Rate as the Star Evolves

As discussed in Section 2, we use a combination of mass-loss rates calculated for hot stars as in Vink et al. (2001) appropriate for OB stars near the main sequence and a combination of rates calculated under the "Dutch" scheme for cooler stars. The rates for RSGs have been tested recently by Mauron & Josselin (2011), who found them consistent with the rates of de Jager et al. (1988). van Loon et al. (2005), however, have considered RSGs that are believed to be dust-enshrouded and have argued for a much higher mass-loss rate. Meynet et al. (2015) pointed out that these mass-loss rates at a given luminosity can differ by more than two orders of magnitude.

In Figure 6 we show the mass-loss rate as a function of luminosity for our fiducial 13 M_☉ models. The star in the ZAMS stage has a low mass-loss rate of 10⁻⁹ M_☉ yr⁻¹, which climbs to a value roughly three times higher as hydrogen in its core is depleted and its color evolves toward the red. At a luminosity of about 3 × 10⁴ L_☉ (log(L/L_☉) = 4.47) and a surface temperature of T_eff = 25,000 K, the mass-loss rate undergoes a sharp increase by a factor of nearly 18, due to the so-called "bi-stability jump"⁵ (Vink et al. 2001). The 26 M_☉ star (plot not shown here) also goes through the "bi-stability jumps" after the TAMS loop and about 10⁶ yr before CC, but because it remains an RSG, its mass-loss evolution is somewhat simpler than that of a 13 M_☉ star. The 13 M_☉ star crosses this transition temperature several times near the TAMS stage. As a result the mass loss rate undergoes both sudden upward and downward transitions within a narrow range of luminosities.

When surface temperatures fall below 5000 K in massive stars, dust begins to form in the stellar wind as the gas cools while it moves away from the star. As the wind mass loss is driven by stellar luminosity, and the mass-loss rate is one of the factors that determines the amount of dust formed, the dust production rate was found to correlate with the bolometric magnitude (Massey et al. 2005) as

$\begin{eqnarray*}&&\mathrm{log}{\dot{M}}_{\mathrm{dust}}=-0.43{M}_{\mathrm{bol}}-12.0.\end{eqnarray*}$

This relation is valid for stars with M_bol < −5, which corresponds to stars more massive than >10M_⊙. This therefore gives a direct handle on how much dust formation is expected from an RSG of given luminosity at the end stage. The dust affects stellar luminosity as well as color, and its effects must be taken into account when comparing with observational color–magnitude data. Dust formation and extinction due to dust are calculated in the Appendix, and the results for a 26 M_☉ star and the corresponding calculated V- and I-band magnitudes are reported in Table 2 and discussed in the following section.

4. Summary and Discussion

We have explored two ZAMS masses claimed in the literature for the progenitor star of SN 2013ej. We simulated the evolutionary stages up to the core-collapse of the progenitor star using MESA code. We used strict mass and temporal resolution controls recommended by Farmer et al. (2016). Our ZAMS 13 M_☉ and 26 M_☉ progenitor models have sufficient mass resolution and have successfully converged. We compared these models with the pre-SN observations by the HST about a decade before collapse (Fraser et al. 2014) and the visual extinction derived for the line of sight of the progenitors of SN 2013ej and SN 2003gd (Maund 2017).

Utrobin & Chugai (2017) claim that a progenitor star with a ZAMS mass of 27.5 ± 2 M_☉, an ejecta mass of 23.1–26.1 M_☉, and a radius of 1500 R_☉ at the pre-SN stage explains better the initial peak of the light curve. In addition, they argue that the bipolar ⁵⁶Ni distribution gives a better fit to the observed photospheric velocity profiles in their models. They prefer the 26 M_☉ ejecta model over a 23 M_☉ ejecta model. However, they do not give any explicit details about their pre-SN models, especially about the various physical parameters that would lead to such a progenitor star. We have run several MESA simulations to show that a star with ZAMS mass of 26 M_☉ will end up with a much smaller pre-SN mass of 22–23.5 M_☉ (including the collapsing core mass), even with a moderate mass-loss rate (η_Dutch = 0.5) and reasonable choices for other parameters for the evolution of the star. Our 26 M_☉ progenitor has a pre-SN radius of ∼1000 R_☉, substantially smaller than what is taken by Utrobin & Chugai (2017), while our 13 M_☉ progenitor has a pre-SN radius of ∼660 R_☉. We also find that a star with a higher ZAMS mass and the same mass-loss scheme produces an even smaller pre-SN progenitor (see also Figure 3 of Ugliano et al. 2012, which shows a similar trend). In addition, a larger ZAMS mass (typically above 30 M_☉) would enter the Wolf–Rayet regime, reducing its chances to explode as a Type II supernova.

As seen in Figure 1, such a massive star would be highly luminous (L_star/L_☉ ∼ a few times 10⁵) for most of its life before the SN explosion. Such a star would require a considerable amount of dust in the CSM to obscure its high luminosity. Even so, an observation in the K band would not be affected by dust extinction. Unfortunately, we do not have any pre-explosion K-band observations for the progenitor of SN 2013ej. As discussed in Section 2.3, our 26 M_☉ model 2 (red curve in Figure 1(b)) shows a large variation in luminosity and spirals in the HRD in the final year before collapse, which results in both V- and I-band variations as large as 0.3 mag. These spirals, if observed, may indicate that the progenitor star is approaching the CC stage. This variation would not be affected by dust extinction. However, such observations are not available for the progenitor of SN 2013ej. Nevertheless, there exist HST observations for the progenitor star taken at about 8 and 10 yr before the explosion. We calculate the extinction values using the CSM profiles calculated from the output of our MESA simulation for our fiducial ZAMS 26 M_☉ model 1 as discussed in the Appendix. These values are tabulated in Table 2. To dim a star as bright as a 26 M_☉ star, the amount of extinction required, A_V, would be of the order of 3 mag. Even though our calculations of the A_V values based on the CSM profiles allow such a high extinction, such a value of A_V is excluded by the observed value of 0.45 ± 0.04 determined by Maund (2017) in the host for SN 2013ej using archival pre-explosion observations through Bayesian analysis of stellar population. Maund also determined a value of A_V = ${0.46}_{-0.03}^{+0.04}$ for another Type IIP SN 2003gd, which occurred in the same host galaxy. As this galaxy is observed face-on, most of this extinction would be due to the dust formed in the CSM around the progenitor star. This poses a difficulty for a ZAMS mass of 26 M_☉ for the progenitor star of SN 2013ej, as advocated by Utrobin & Chugai (2017).

The multi-wavelength observations of the SN, observationally derived pre-explosion values of dust extinction, as well as our model results point toward a ZAMS mass of 13 M_☉ for the progenitor of SN 2013ej rather than a star of about twice that mass. In a future paper, we will discuss detailed simulations of post-explosion multi-waveband light curves and the evolution of Fe line velocity of the two stars considered here to compare with available observational data for additional scrutiny of the mass of the progenitor star.

We thank the directors and the staff of the Tata Institute of Fundamental Research (TIFR) and the Homi Bhabha Center for Science Education (HBCSE-TIFR) for access to their computational resources. This research was supported by a Raja Ramanna Fellowship of the Department of Atomic Energy (DAE), Govt. of India to Alak Ray and a DAE postdoctoral research associateship to G.W. G.W. thanks Rob Farmer for valuable feedback through private communication. The authors thank the anonymous referee for his/her constructive comments. The authors thank Ajay Dev and Viraj Meruliya, students of the National Initiative of Undergraduate Science (NIUS) program at HBCSE (TIFR). The authors acknowledge the use of NASA's Astrophysics Data System.

Software: MESA r-10398 (Paxton et al. 2011, 2013, 2015, 2018), Anaconda Spyder (Python 3.6).

Note added in proof.

After we posted our paper to the e-print (arXiv:1911.12831v1), Matteo Cantiello drew our attention to Figure 45 of Paxton et al. (2013, ApJS, 208, 4) where the authors describe unstable radial pulsations in RSGs of a 25 M_sun star driven by the kappa-mechanism in the hydrogen ionization zone in the envelope (opacity due to H-recombination). This is due to the high L/M ratio leading to a pulsation instability. For non-rotating RSGs with solar metallicity, their pulsation period is reported to be in the range 1–8 yr. We note again here that the spirals seen in our HRD (Figure 1(b)) occur only when dm = 0.01 M_sun, but not when dm = 0.007 M_sun for 26 M_sun star, irrespective of other parameters. The model time-steps are comparable (<0.01 yr) for both these mass resolutions at the evolutionary stages in consideration.

Appendix: Extinction by Circumstellar Material

Here, we discuss our calculations for the extinction by dust particles (and the underlying assumptions) in the CSM formed by the mass lost in stellar winds (values listed in Table 2). The dust particles can form only when the wind moves sufficiently far from the star to avoid dust particles being eradicated by radiation near the star, and cools down to a temperature below the dust destruction temperature. The typical value used in the literature for the sublimation temperature at which the dust particles will be destroyed is ${T}_{d,\max }$ = 1500 K (Kochanek et al. 2012; Das & Ray 2017); however, graphite dust particles can withstand temperatures as high as 2000 K (see, e.g., Fox et al. 2010). We calculate the minimum distance (R_min) from the stellar surface at which the dust can form, where the temperature of the CSM falls below ${T}_{d,\max }$ . Das & Ray (2017) determined this distance by simple adiabatic cooling of the wind ( ${{TV}}^{\gamma -1}={\rm{constant}}$ ) due to expansion with an adiabatic monotonic gas ( $\gamma =5/3$ ). In this case, the minimum distance the gas has to travel to cool down below the graphite destruction temperature, supposing that the dust temperature is same as that of the surrounding gas, is given by

$\begin{eqnarray}&&{R}_{\min }={R}_{\mathrm{star}}\sqrt{\displaystyle \frac{{T}_{\mathrm{star}}}{{T}_{d,\max }}}.\end{eqnarray} \tag{ 1 }$

We use a more robust way of calculating the minimum distance using the heat balance equation for dust grains in thermal equilibrium, where the rate of dust heating due to stellar radiation is equal to its rate of cooling due to thermal emission, written as (see, e.g., Draine & Lee 1984; Kruegel 2003)

$\begin{eqnarray}&&{\int }_{0}^{\infty }{Q}_{\mathrm{abs}}(a,\nu ){J}_{\nu ,\mathrm{star}}d\nu ={\int }_{0}^{\infty }{Q}_{\mathrm{abs}}(a,\nu ){B}_{\nu }({T}_{d})d\nu .\end{eqnarray} \tag{ 2 }$

The absorption efficiency Q_abs in the mid-infrared where the grain emission takes place can be approximated by a power law using a simple spherical grain model (Kruegel 2003):

$\begin{eqnarray}&&{Q}_{\mathrm{abs}}(a,\nu )={{aQ}}_{0}{\nu }^{\alpha }\end{eqnarray} \tag{ 3 }$

where Q₀ is a constant. The value of α, the emissivity index, ranges from 1 to 2; however, a value of 2 is favored in the literature. The right-hand side of Equation (2), which is the cooling rate, ${ \mathcal W }$ ( $\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ ), then becomes

$\begin{eqnarray}\begin{array}{rcl}{ \mathcal W } & = & {\displaystyle \int }_{0}^{\infty }{Q}_{\mathrm{abs}}(a,\nu ){B}_{\nu }(T)d\nu \\ & = & \displaystyle \frac{2h}{{c}^{2}}\ {{aQ}}_{0}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\nu }^{3+\alpha }}{{e}^{h\nu /{kT}}-1}d\nu \\ & = & \displaystyle \frac{2h}{{c}^{2}}\ {{aQ}}_{0}\ {\left(\displaystyle \frac{{kT}}{h}\right)}^{4+\alpha }{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{y}^{3+\alpha }}{{e}^{y}-1}{dy}.\end{array}\end{eqnarray} \tag{ 4 }$

The integral in the above equation yields approximate values (Kruegel 2003)

$\begin{eqnarray}{\int }_{0}^{\infty }\displaystyle \frac{{y}^{3+\alpha }}{{e}^{y}-1}{dy}\simeq \left\{\begin{array}{ll}6.49, & \,\mathrm{if}\ \alpha =0\\ 24.89, & \,\mathrm{if}\ \alpha =1\\ 122.08, & \,\mathrm{if}\ \alpha =2.\end{array}\right.\end{eqnarray} \tag{ 5 }$

Then with the favored value of α = 2, Equation (4) becomes (Lequeux 2005)

$\begin{eqnarray}&&{ \mathcal W }=4.6\times {10}^{-11}\displaystyle \frac{a}{0.1\,\mu {\rm{m}}}{T}^{6}\,({\rm{erg}}\,{{\rm{cm}}}^{-2}\,{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 6 }$

If we assume that the emitting dust is highly absorbing in UV and the visible, which are the only relevant wavelengths, then ${Q}_{\mathrm{abs}}(a,\nu (\mathrm{UV}))\,\simeq$ 1. In this case, the left-hand side of Equation (2) will be simply equal to ${L}_{* }/4\pi {R}^{2}$ , where R is the distance of the dust grain from the star. Combining this with Equations (2) and (6), we get the minimum distance at which a dust grain of size a with temperature ${T}_{d,\max }$ could exist as

$\begin{eqnarray}&&{R}_{min}[{\rm{cm}}]=\sqrt{\displaystyle \frac{{L}_{\ast }[{\rm{erg}}\,{{\rm{s}}}^{-1}]}{4\pi \times 4.6\times {10}^{-11}\,{T}_{d,max}^{6}(a/0.1\,\mu {\rm{m}})}}.\end{eqnarray} \tag{ 7 }$

Any ejected mass closer to the star than this will not be able to form dust grains.

On the other hand, mass that has moved too far away from the star will be too dilute to contribute to the extinction. We calculate the visual extinction assuming only the contribution from dust particles existing in the CSM between R_min and a maximum distance of R_max ≈ 10¹⁵ cm (similar to Utrobin & Chugai 2017), beyond which the dust will be part of the interstellar medium (ISM), and perhaps too dilute to contribute significantly to the extinction for our models. The visual extinction can be calculated for a distribution of grains between minimum and maximum grain sizes (a_min and a_max) using the equation (Perna et al. 2003)

$\begin{eqnarray}&&{A}_{V}={\int }_{{R}_{min}}^{{R}_{max}}{dr}\displaystyle \sum _{i}{\int }_{{a}_{min}}^{{a}_{max}}{da}\,\pi {a}^{2}\displaystyle \frac{{{dn}}_{i}}{{da}}{Q}_{{\rm{abs}},i}(a,V)\end{eqnarray} \tag{ 8 }$

where we neglect the angle-dependent term involving Q_scat(a, ν) because MESA is a 1D code. The size distribution dn_i/da of dust grains of type i, per unit volume per grain size a per H-atom (in units of cm⁻⁴), typically follows an MRN (Mathis et al. 1977) power law given as

$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{i}}{{da}}={A}_{i}{n}_{{\rm{H}}}{a}^{-\beta },\,{a}_{min}\leqslant a\leqslant {a}_{max}\end{eqnarray} \tag{ 9 }$

where the typical values for the ISM for a_min, a_max, and the index β are 0.005 μm, 0.25 μm, and 3.5, respectively. This distribution would change for extreme conditions close to a highly luminous RSG star. Radiation from a star fragments larger grains into smaller grains, while much smaller grains are easily sublimated away. This leads to a flatter grain size distribution over a narrower range of grain sizes. Following Perna et al. (2003, Figure 3), we see that after a considerable amount of time, there would exist only graphite grains distributed between grain sizes a_min and a_max of 0.15 μm and 0.22 μm, respectively, and with β ≈ 1.4. The coefficients A_i in Equation (9) are related to the dust-to-gas ratio (by mass), f_d (see Laor & Draine 1993; Perna et al. 2003). For dust consisting of only graphite, we have

$\begin{eqnarray}&&{f}_{d}=\displaystyle \frac{{m}_{\mathrm{dust}}}{{m}_{\mathrm{gas}}}=\displaystyle \frac{4\pi {a}_{\max }^{4-\beta }}{3{m}_{{\rm{H}}}(4-\beta )}\left[1-{\left(\displaystyle \frac{{a}_{\min }}{{a}_{\max }}\right)}^{4-\beta }\right]{A}_{\mathrm{gra}}{\rho }_{\mathrm{gra}}\end{eqnarray} \tag{ 10 }$

where ${\rho }_{\mathrm{gra}}\approx 2.26\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ is the density of graphite grains (Draine & Lee 1984) and m_H is mass of a hydrogen atom. We conservatively assume that only about 50% of the available ¹²C has formed dust. We determine f_d to be a few times 10⁻⁴, using the surface composition of the star. We find that the values of Q_abs(a, V) are somewhat independent of a for the above range of grain sizes. We used an averaged value for Q_abs(V) = 1.5 in this range using data provided online by Laor & Draine (1993). Within this formalism, we can now reduce Equation (8) to the form

$\begin{eqnarray}&&{A}_{V}=\displaystyle \frac{\pi {A}_{{\rm{gra}}}{Q}_{{\rm{abs}}}(V)}{{m}_{{\rm{H}}}}{\left(\displaystyle \frac{{a}^{3-\beta }}{3-\beta }\right)}_{{a}_{min}}^{{a}_{max}}{\int }_{{R}_{min}}^{{R}_{max}}{dr}\,\rho (r)\end{eqnarray} \tag{ 11 }$

where ρ(r) is the CSM density profile calculated by using the mass lost by stellar winds with a constant speed of 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (see Equation (4) of Das & Ray 2017). We calculated A_V using R_min determined from the adiabatic cooling of the gas as in Equation (1) or the heat balance as in Equation (7) up to R_max. While choosing R_min for the heat balance equation, which depends on the grain size, we chose the distance corresponding to a_min. The values of the extinction in the I band, A_I, were calculated using the relation given by Cardelli et al. (1989), A_I/A_V = 0.479.

We also calculated the visual extinction values using Equation (4) of Walmswell & Eldridge (2012), where they assume that the CSM dust consists of low-density silicates ( ${\rho }_{\mathrm{sil}}=2.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ), for comparison, using our CSM profiles. We used the surface composition of the star to find the mass of the dust. We considered the elements that constitute silicate dust grains. To find the R_min value we used the heat balance formalism explained above, with ${T}_{d,\max }$ = 1500 K, which is typical for silicates. We used the same value of R_max (10¹⁵ cm) from the surface of the star as before. These extinction values for the V and I bands are listed in Table 2 and discussed in the Section 3.

Type IIP Supernova Progenitors. II. Stellar Mass and Obscuration by the Dust in the Circumstellar Medium

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction