© ESO 2020

1. Introduction

It has been established that most of the stars with initial masses greater than 8 M_⊙ finish their evolution in a violent explosion (Woosley et al. 2002; Heger et al. 2003), known as a core-collapse supernova (CCSN). These CCSNe are observationally classified according to their spectral characteristics (Filippenko 1997). Type II supernovae (SNe II) show strong and prominent P-Cygni hydrogen lines. A subsequent class division was introduced based on light curve (LC) decline rates after maximum into II-Plateau (IIP) and II-Linear (IIL) events. However, recent studies have questioned this subdivision and propose the existence of a continuous sequence of LC slopes among SNe II (Anderson et al. 2014; Sanders et al. 2015; Galbany et al. 2016; Rubin & Gal-Yam 2016). Therefore, throughout this paper we use “SNe II” to refer to these two groups together. A further classification exists for SNe that show hydrogen lines: the SN 1987A-like events, displaying unusually long-rising LCs (e.g. Taddia et al. 2012, 2016); the type IIn, showing narrow emission features in their spectra (e.g. Schlegel 1990; Taddia et al. 2013); and the type IIb, which show hydrogen features at early times while later such lines disappear (e.g. Filippenko et al. 1993). A recent study analysed the possible existence of a continuum between the SNe II and IIb in terms of their photometric properties (Pessi et al. 2019). However, these latter authors found clear differences between the two subgroups. These three groups show characteristics sufficiently distinct from SNe II (as defined above) that we do not consider them in the present work. Only SNe II are studied in this paper.

It is generally assumed that the progenitors of SNe II are massive stars that have retained a significant fraction of their hydrogen-rich envelopes before explosion. In addition to this, it has been shown that these assumptions are consistent with detections of progenitor stars in pre-explosion images. These detections have constrained the progenitors to be red supergiant (RSG) stars in the zero-age main sequence (ZAMS) mass range of ∼8–20 M_⊙ (Smartt 2015; Davies & Beasor 2020).

Although SNe II are the most common type of SN in nature (Arcavi et al. 2010; Li et al. 2011), significant gaps remain in our knowledge of the different processes involved. Moreover, important discrepancies can be found in the literature regarding the masses of the progenitors depending on the different methods used for the analysis. Archival images provide an opportunity to detect progenitor stars in images taken prior to SN explosions. The luminosity and effective temperature of these progenitors can be obtained from fits to the spectral energy distribution or using bolometric corrections to convert single-band flux into luminosity (see Van Dyk et al. 2012; Davies & Beasor 2018, among others). Once stars are located in a Hertzsprung-Russell (HR) diagram, their mass in the ZAMS is estimated by comparison with stellar evolutionary tracks. The acquisition of late-time imaging is the next step in the analysis of pre-explosion observations to confirm a progenitor identification through its disappearance (Maund et al. 2014b; Folatelli et al. 2016, among others). Despite being the most immediate method to determine the types of stars that produce the different types of SNe, this method can only be applied to the most nearby SNe (out to distances of d ≲ 30 Mpc) because it requires images of high enough resolution and sensitivity.

Nebular-phase spectral modelling can also be used to constrain the progenitor masses of SNe. During the photospheric phase, the LC is mostly powered by reemission of the energy deposited by the shock wave. As time goes on, hydrogen recombination occurs in different layers of the object as a recombination wave moves inward. After the hydrogen has recombined, the envelope becomes transparent, the inner ejecta become visible and the nucleosynthesis yields can be analysed. Therefore, using the dependency of oxygen production on progenitor initial mass, it is possible to distinguish between different progenitors (see, e.g. Jerkstrand et al. 2012, 2014).

However, one of the most used methods to analyse the progenitor properties is the hydrodynamical modelling of SN LCs. It is well known that LCs are extremely sensitive to the physical properties of their progenitors (final masses and radii), as well as the properties of the explosion itself (released energy, amount of synthesised radioactive nickel and its distribution; see e.g. Shigeyama & Nomoto 1990; Bersten et al. 2012, among others). The main problem in using this method is that LC modelling cannot always provide a unique solution for the ejecta mass of SNe II. There is a degeneracy among some progenitor properties when reproducing the observations. This means that progenitors with different physical properties can have similar photometric and spectroscopic properties (Bersten et al. 2011; Dessart & Hillier 2019; Goldberg et al. 2019; Goldberg & Bildsten 2020).

Following the degeneracies involved in constraining SN II properties, the masses inferred by hydrodynamic simulations are usually much larger than those estimated from pre-SN imaging (see Utrobin & Chugai 2009, 2017, among others). Recently, Morozova et al. (2018) and Eldridge et al. (2019) inferred initial masses for a group of SNe II from LC modelling, and found that their results are mostly consistent with those from pre-explosion data. However, they did not take into account the ejecta velocities in their modelling and therefore the parameters derived are not unequivocally determined. Additionally, Martinez & Bersten (2019, MB19 hereafter) presented hydrodynamic modelling of LCs and photospheric velocities of six objects with confirmation of the progenitor star, and found that in most cases masses inferred by both methods were compatible. They also noted that the degeneracy in some physical parameters may be the reason for the differences found in the literature. Hydrodynamical modelling can be applied to large distances and large samples as well, contrary to other methods which are more difficult to employ or present limitations when used for sources at great distance. For this reason, efforts must be made in order to solve these issues and enable us to accurately constrain the physical and explosion properties of SNe II.

In order to acquire a detailed outlook, we generated a grid of hydrodynamic models in the parameter space and a quantified fitting procedure. Similar grids and techniques were published in the recent works of Morozova et al. (2018) and Eldridge et al. (2019). The aim of the current work is to develop a procedure based on Markov chain Monte Carlo (MCMC) methods to construct the posterior distribution of the parameters involved given the LC and the photospheric velocity evolution together, and then to constrain the progenitor properties of SNe II. In this paper we propose to verify the robustness of this method with several tests. First, we derive the progenitor parameters of the same sample as in MB19 for comparison. In this latter work, the authors present detailed modelling of six objects using double polytropic progenitor structures and find their optimal models by fitting to the observations by eye. The current work is different in that we use pre-SN structures based on stellar evolution calculations and a more robust statistical analysis. The same hydrodynamic code is used in both cases. The sample includes SN 2004A, SN 2004et, SN 2005cs, SN 2008bk, SN 2012aw, and SN 2012ec. We chose these objects as they represent some of the best-observed SNe II with enough photometric and spectroscopic monitoring during all the phases of evolution, detection of the progenitor star in pre-explosion images, and confirmation of the progenitor through its disappearance in late-time images.

We also test the validity of our method by comparing the initial masses derived from our fitting to those based on the analysis of the progenitor star in pre-explosion images for the same sample described above. For this purpose, we also include SN 2017eaw and SN 2018aoq as these are the last SNe II to be discovered and analysed with this method. It is important to mention that these progenitor detections have not yet been confirmed. However, despite the lack of confirmation, a nebular spectral analysis is available for SN 2017eaw and SN 2018aoq has not yet been analysed by hydrodynamical modelling, which make these two SNe relevant in our study for further validation between different methods. Finally, whenever possible, we compare with the results from late-time spectral modelling found in the literature as well.

This paper can be considered as a companion of a forthcoming paper which analyses a large sample of SNe II using the same grid of simulations and fitting procedure. In this work we present the stellar evolution calculations we use to obtain the structure of stars at core collapse, the grid of explosion models, the technical part of the fitting procedure, and a sanity check.

Our paper is organised as follows. We first present a description of our hydrodynamic code and pre-SN models used in Sect. 2. We then provide a brief description of the sample of SNe (Sect. 3). In Sect. 4, we present the fitting procedure used and its characteristics. In Sect. 5 we present comparisons between our results and those of previous studies using different methods. Section 6 provides a discussion on LC degeneracies and the limitations of our models, and finally, we summarise our conclusions in Sect. 7.

2. Hydrodynamic simulations

Theoretical LCs are calculated using a 1D Lagrangian hydrodynamical code that simulates the explosion of the SN and produces bolometric LCs and photospheric velocities of SNe (Bersten et al. 2011). The explosion is simulated by injecting a certain amount of energy near the centre of the progenitor object, which produces a powerful shock wave that propagates through the star transforming the thermal and kinetic energy of the matter into energy that can be radiated from the stellar surface.

The code assumes that the fluid motion can be described as a 1D, radially symmetric flow and that radiation and matter are strongly coupled, that is to say that local thermodynamical equilibrium (LTE) describes the radiative transfer. The code uses opacity tables calculated assuming LTE and a medium at rest (see Bersten et al. 2011, for details).

A pre-supernova model in hydrostatic equilibrium that simulates the conditions of the star before exploding is necessary to initialise the explosion. We use the public stellar evolution code MESA¹ version 10398 (Paxton et al. 2011, 2013, 2015, 2018, 2019) to obtain non-rotating solar-metallicity RSG models. Each stellar model is evolved from the pre-main sequence until core collapse, which we take as the time when any location inside the stellar model reaches an infall velocity of 1000 km s⁻¹. We use Ledoux criterion for convection and set a mixing-length parameter of α_mlt = 2.0, exponential overshooting parameters f_ov = 0.004 and f_ov, D = 0.001, a semiconvection efficiency α_sc = 0.01 according to Farmer et al. (2016), and thermohaline mixing with coefficient α_th = 2 (Kippenhahn et al. 1980). For every model, we use the “Dutch” wind scheme (de Jager et al. 1988; Vink et al. 2001; Glebbeek et al. 2009) defined in the MESA code with an efficiency η = 1. Figure 1 shows the density profiles for the 17 pre-SN models used in this work. Some additional properties of the pre-SN models are given in Table 1. The 9, 10, and 11 M_⊙ progenitor models were calculated up to the end of core carbon burning since the evolution to core collapse for these stars is computationally expensive. The absence of later stages of evolution for these three pre-SN models is noted in Fig. 1 as the inner core density is about three orders of magnitude lower than that of the other stellar models. However, only the inner core changes from this part of evolution until core collapse and, additionally, we remove this part for the hydrodynamic calculation as we assume it will collapse and leave a compact remnant.

Fig. 1. Density profiles for the pre-SN models used in this work. The 9, 10, and 11 M⊙ models were calculated up to the end of core carbon burning since the evolution to core collapse for these stars is computationally expensive.

The determination of the physical properties of SNe II is based on describing the bolometric LC and the expansion velocity at the photospheric layers by means of comparing hydrodynamical models with observations. The morphology of the LC and the evolution of the photospheric velocity are related to physical properties of the progenitor star and the explosion itself, such as the mass (M_presn) and radius (R) prior to explosion, the energy that is transferred to the envelope after core-collapse (denoted “explosion energy”; E), the amount of ⁵⁶Ni synthesised in the explosion (M_Ni) and its degree of mixing into the outer layers of the ejecta (⁵⁶Ni mixing). Additionally, M_presn and R depend on the evolution of the star and are directly connected to the progenitor initial mass (M_ZAMS). Figure 2 shows the ejecta mass (pre-SN mass minus the compact remnant mass) and final radius for every model as a function of M_ZAMS. Ejecta masses cover a range of 7.9–15.7 M_⊙, while progenitor radii are found in the range of 445–1085 R_⊙.

Fig. 2. Ejecta masses (top panel) and progenitor radii (bottom panel) for the pre-SN models as a function of the initial mass.

We computed a grid of explosion models in the parameter space covering a M_ZAMS range of 9–25 M_⊙ in intervals of 1 M_⊙ (which represent the ranges of ejecta masses and final radius as described above) and explosion energies between 0.1 and 1.5 foe (1 foe ≡ 10⁵¹ erg) in steps of 0.1 foe with the exception of the largest masses and lowest energies due to numerical difficulties. For the 20 M_⊙ and 21 M_⊙ models, the lowest explosion energies are 0.2 foe and 0.3 foe, respectively. Models of 22 M_⊙ and 23 M_⊙ were calculated for explosion energies higher than 0.4 foe, and for the 24 M_⊙ and 25 M_⊙ models, only explosion models with energies higher than 0.5 foe are available. We also consider M_Ni in the range of 0.01–0.08 M_⊙ in intervals of 0.01 M_⊙, together with nickel masses of 0.0001 and 0.005 M_⊙ to be consistent with the lowest SN II estimated ⁵⁶Ni masses in the literature (Müller et al. 2017; Anderson 2019). To account for the effect of the spatial distribution of ⁵⁶Ni within the ejecta we consider three degrees of ⁵⁶Ni mixing for each model: out to the 20%, 50% and 80% of the pre-SN structure in mass coordinate.

Recent studies have shown that the interaction of the ejecta with a circumstellar material (CSM) shell surrounding the progenitor star can affect the early LC of SNe (González-Gaitán et al. 2015; Moriya et al. 2017; Förster et al. 2018; Morozova et al. 2018). In the current work, we do not attempt to characterise the CSM. In this way, we do not include any CSM surrounding the star in our explosion models and focus on deriving intrinsic properties of the progenitors. Therefore, we restrict the analysis of the observed LC to times later than 30 days after explosion (see Sect. 5).

3. Data sample

In this work, we aim to determine whether or not our method to infer physical properties of SNe II is consistent with other methods in the literature. For this purpose we use some of the most well-studied SNe II. In particular, we use the same sample analysed in MB19 defined as SNe II that: (a) have well-constrained pre-explosion progenitor detections, (b) have post-explosion images confirming the disappearance of the progenitor, and (c) have sufficiently well-sampled photospheric-phase observations to enable accurate modelling fits. In addition, we include two further objects in the sample, SN 2017eaw and SN 2018aoq, as these are the most recently published SNe II with progenitor identification in pre-explosion images. However, these identifications have not yet been confirmed by post-explosion images. Information for these two SNe is presented in Sects 3.1 and 3.2. Details of the other objects can be found in MB19 and references therein.

Following our sample definition, we calculate their bolometric LCs using the correlation between bolometric correction and colours inferred in Bersten & Hamuy (2009), which allows us to calculate bolometric luminosities using only two optical filters. We use the dispersion values listed in Table 1 of Bersten & Hamuy (2009) and the uncertainties in colours to estimate the uncertainty in the bolometric luminosities via error propagation. Neither the uncertainties in the distance nor extinction were propagated to the bolometric LCs as we include an additional parameter in the MCMC procedure that models these uncertainties (Sect. 4). The values of host extinctions, distances, and explosion times are the same as those presented in Table 1 of MB19 with the exception of the distance assumed for SN 2004et. Here we recalculate the bolometric LC of SN 2004et using the most recent estimation of the distance to its host galaxy NGC 6946 (see Sect. 3.1).

In addition to the bolometric luminosity, an estimation of the photospheric velocity is also needed to compare with the models. This velocity can be estimated through the measurement of certain spectroscopic lines. We use the Fe IIλ5169 Å line because this line is formed in internal regions of SNe and has been proposed as a good estimator of the photospheric velocity (see Dessart & Hillier 2005, and Sect. 6.2 for discussion).

With the aim of assessing the nature of progenitors and comparing our constraints with results obtained from pre-explosion information and late-time spectral modelling, we set out to perform a detailed modelling of the available observations. Figure 3 shows the bolometric and absolute V-band LCs of SN 2017eaw and SN 2018aoq in comparison to the sample from MB19. It seems like SN 2017eaw is intermediate in luminosity between, for example, SN 2012aw and the more luminous SN 2004et, while SN 2018aoq appears to be an intermediate case between the normal and the low-luminosity SNe II (O’Neill et al. 2019). From this comparison we also note that the last two V-band observations of SN 2018aoq constrain the end of the plateau phase.

Fig. 3. Bolometric LCs (top panel) and absolute V-band LCs (bottom panel) of SN 2017eaw and SN 2018aoq in comparison to well-studied SN 2004A, SN 2004et, SN 2005cs, SN 2008bk, SN 2012ec, and SN 2012aw. Absolute curves were computed using the distances and reddenings in Table 1 of MB19, except in the case of SN 2004et for which we use the new estimates of the distance (see text).

3.1. SN 2017eaw

SN 2017eaw was discovered on 2017 May 14.238 UT in NGC 6946 at an unfiltered magnitude of 12.8 mag (Wiggins 2017). This object was classified as a young SN II by Cheng et al. (2017), Xiang et al. (2017), and Tomasella et al. (2017). Wiggins (2017) also observed the site of explosion on 2017 May 12.20 UT but nothing was visible. This limitation in the detection restricts the uncertainty in the explosion epoch (t_exp) to only one day. We adopt t_exp as JD 2457886.72 ± 1.01 (Rui et al. 2019).

Pre-explosion images of the SN location were obtained with the Hubble Space Telescope (HST) and the Spitzer Space Telescope covering the last ∼13 yr before explosion. The progenitor can be detected in eight optical and infrared bands making it one of the most-characterised progenitors to date. Several works studying the progenitor properties have been published. Kilpatrick & Foley (2018) found a RSG progenitor compatible with an initial mass of ∼13 M_⊙. They detected an increase in its 4.5 μm luminosity over the final 3 yr before the explosion and argue that it is a signature of circumstellar dust near the progenitor star. Moreover, Rui et al. (2019) found a narrow and blueshifted Hα emission component in a spectrum taken a few hours after discovery that disappeared in less than two days, suggesting the presence of a CSM shell. The authors propose that the progenitor could have experienced dramatically enhanced mass-loss during the last 1–2 yr before explosion. They also found a RSG progenitor with an initial mass of 12 ± 2 M_⊙. On the other hand, Van Dyk et al. (2019) established that the progenitor was a dusty, luminous RSG consistent with an initial mass of ∼15 M_⊙. Unfortunately, these three studies assumed different values for the distance to the object. Recently, Murphy et al. (2018), Anand et al. (2018), and Van Dyk et al. (2019) used archival HST data taken in the outer regions of NGC 6946 to measure the tip of the red giant branch and infer the distance to the galaxy. All these studies arrive at the same value for the distance within the uncertainties. In addition, Eldridge & Xiao (2019) used these new measurements to re-evaluate the final luminosity for some SN progenitors in NGC 6946 (SN 2004et, among others). With the new distance, these latter authors estimate that the initial mass of the progenitor of SN 2017eaw is 14${}_{-3.5}^{+3.0}$ M_⊙, consistent with the results from Van Dyk et al. (2019) for the same distance. In the present paper we set the distance at 7.73 ± 0.78 Mpc and the total extinction at 0.941 mag (see discussion in Van Dyk et al. 2019).

Additionally, using nebular-phase spectral modelling, Van Dyk et al. (2019) and Szalai et al. (2019) analysed late-time spectra of SN 2017eaw and found that a progenitor with initial mass near 15 M_⊙ is most consistent with observations.

In this study, we take photometric data and Fe IIλ5169 Å line velocities from Szalai et al. (2019). Additional optical photometry can also be found in Tsvetkov et al. (2018), Van Dyk et al. (2019), Rui et al. (2019), and Buta & Keel (2019).

3.2. SN 2018aoq

SN 2018aoq was discovered on 2018 April 01.43 in the galaxy NGC 4151 by the Lick Observatory Supernova Search (LOSS) at the unfiltered magnitude of 15.3 mag (Nazarov et al. 2018). Yamanaka (2018) carried out spectroscopic observations on 2018 April 02 using the Hiroshima One-shot Wide-field Polarimetry (HOWPol) installed on the 1.5-m Kanata telescope and found a spectrum dominated by a blue continuum and the Hα line with P-Cygni profile consistent with a young SN II. O’Neill et al. (2019) presented optical imaging data and spectra using a combination of Asteroid Terrestrial-impact Last Alert System (ATLAS, Tonry et al. 2018; Smith et al. 2020), the 2.0 m Liverpool Telescope (LT), and the 2.5 m Nordic Optical Telescope (NOT) as part of the NOT Unbiased Transient Survey (NUTS). They also presented an estimation of the total reddening and a new estimation of the distance that is in good agreement with the value based on geometric methods. In addition, from the ATLAS non-detection on 2018 March 28, the explosion epoch is well constrained to within four days. We adopt the same values of distance (d = 18.2 ± 1.2 Mpc), reddening (E(B − V)_tot = 0.04 mag), and explosion epoch (JD 2458208.5) as O’Neill et al. (2019) in what follows.

Archival pre-explosion images of the SN site are available. These images were taken with the Wide Field Camera 3 (WFC3) on board the HST approximately 2 yr before explosion. From these observations, O’Neill et al. (2019) detected a source at the SN location in four bands: F350LP, F555W, F814W and F160W. From fits to the spectral energy distribution of the progenitor candidate, these latter authors found a luminosity range of log(L/L_⊙) = 4.56–4.83 and an effective temperature of T = 3500 ± 150 K, implying an M-type red supergiant progenitor. Using single and binary star models, they conclude that the explosion of a star with a ZAMS mass of 10 ± 2 M_⊙ is the most favoured scenario.

4. Fitting procedure

In MB19, the authors derive physical parameters for the sample presented in Sect. 3, with the exception of SN 2017eaw and SN 2018aoq, based on visual comparisons between observations and hydrodynamic simulations. In the present work we aim to develop and test a robust fitting procedure that automatically obtains optimal solutions for the fitted parameters in a statistically sound manner, following Förster et al. (2018). To this end, we computed a large grid of models of bolometric LCs and photospheric velocity evolution. The range of physical parameters considered is described in Sect. 2. However, these models may not sufficient when trying to fit SN observations using statistical inference techniques. We need to be able to interpolate between models with different physical parameters. First we define the set of parameters for which we want to compute an interpolated model, that is, M_ZAMS, E, M_Ni, and ⁵⁶Ni mixing. Then, in our grid of simulations we find the models with the closest values in all the physical parameters. The interpolated LC and velocity evolution are calculated using Eqs. (2) and (3) from Förster et al. (2018). This is a robust and quick method which attempts to provide a scale–free interpolation, relevant when combining variables with different physical dimensions. It also allows for irregular grids of models in the space of parameters to be used. We show examples of interpolated models in Fig. A.1.

Having a powerful interpolation method, we can attempt to infer the physical parameters using Bayesian statistics, that is, by computing the posterior probability of the model parameters given the observations and assuming prior distributions. To do this we use a MCMC sampler which uses an affine invariant approach (Goodman & Weare 2010). This method estimates the properties of a distribution by examining random samples from the distribution which are generated using parallel Markov chains. The characteristic of the Markov chains is that each random sample is used to generate the next random sample. While each new sample depends on the one before it, new samples do not depend on any samples before the previous one (van Ravenzwaaij et al. 2018). This method is implemented via emcee in python (Foreman-Mackey et al. 2013).

We run the MCMC sampler using flat distributions as priors for the following parameters: t_exp, M_ZAMS, E, M_Ni, and ⁵⁶Ni mixing. We allow the sampler to run within the observational uncertainty of the t_exp. In addition, we define a variable named scale for which we use a Gaussian prior. The scale parameter multiplies the bolometric luminosity by a constant factor to allow for errors in the bolometric LC due to the uncertainties on the distance and extinction. For this parameter we use a Gaussian prior centred at 1.0 with a standard deviation equal to the uncertainties in the distance estimation. It is worth emphasising that the uncertainties on the distance generally dominate over those on extinction. In this way, we assume that the distance errors include the extinction errors. Additionally, we constrained the scale parameter to ±1σ. The reason for this choice is as follows. Progenitor initial mass is not an observational parameter. For example, in pre-SN imaging, the detection of the progenitor star gives the observed luminosity of the object close to the explosion epoch, which is then converted to an initial mass using stellar evolution calculations. Thus, the initial mass depends on the luminosity of the progenitor, and this latter depends on the distance to the object and the host galaxy extinction, among others. Constraining the scale parameter to ±1σ we make sure that the physical properties we derive are consistent with the distance and extinction estimates. In this context we can compare our mass estimation with that from pre-SN imaging (or any other method) as we assume the same range of distances and total extinction.

We use 400 parallel samplers (or walkers) and 10 000 steps per sampler, with a burn-in period of 1000 steps. These numbers were set via trial and error through checking randomness and stationarity of the chains. The walkers are randomly initialised covering the entire parameter space.

An example corner plot with the posterior probability distributions can be visualised in Fig. 4. It is seen that the marginal distributions for the t_exp, scale, and ⁵⁶Ni mixing are limited. The posterior distribution marginalised over the ⁵⁶Ni mixing parameter is restricted to values above 0.2 because of limitations in our model. Additionally, the t_exp and the scale are constrained to the uncertainties in the explosion epoch and distance, respectively. If we relax the priors for these two parameters we obtain a similar marginal distribution of the physical parameters (see Fig. A.9). The corner plots for the entire sample are shown in Appendix A.2, and the models drawn from the posterior distribution for the SN sample are shown in Figs. 5 and 6. Examples of autocorrelation plots and trace plots are presented in Figs. A.10 and A.11, respectively.

Fig. 4. Corner plot of the joint posterior probability distribution of the parameters for SN 2012ec.

Fig. 5. Comparison between models (solid lines) and observations (filled dots) for our SN sample. We show 50 models randomly chosen from the posterior probability distribution. Left: bolometric LC. Right: evolution of the photospheric velocity. From top to bottom: SN 2004A, SN 2004et, SN 2005cs, and SN 2008bk. The grey shaded region shows the early data we removed from the fitting. For SN 2004et we show the results using d = 5.9 ± 0.4 Mpc to calculate the bolometric LC.

Fig. 6. Comparison between models (solid lines) and observations (filled dots) for our SN sample. We show 50 models randomly chosen from the posterior probability distribution. Left: bolometric LC. Right: evolution of the photospheric velocity. From top to bottom: SN 2012aw, SN 2012ec, SN 2017eaw, and SN 2018aoq. The grey shaded region shows the early data we removed from the fitting.

5. Results

We have developed a fitting procedure to derive physical parameters of SNe II from the hydrodynamic modelling of LCs and photospheric velocities using the models described in Sect. 2 and the method from Sect. 4.

It is known that the presence of a dense CSM affects the early evolution of SNe II with almost no effect at later epochs (t ≳ 30 days after explosion; see Figs. 3 and 4 of Morozova et al. 2018) where the evolution is dominated by the hydrogen recombination and radioactive decay. This is correct if the ejecta interacts with a low-mass CSM (Englert, in prep.). The effect of the ejecta-CSM interaction can dominate the behaviour of the early LC suggesting that the general characteristics of the CSM have a significant role in this early phase, and not the progenitor properties as mentioned by Utrobin & Chugai (2008). As we have noticed in Sect. 2, we did not include CSM in our set of progenitor models as we are interested in deriving global properties of the SN progenitor instead of analysing the CSM properties. Therefore, we do not consider the first 30 days of evolution of the observed LC in our fitting procedure. As a consequence, differences are expected between the models and the observations during the cooling phase, which is strongly affected by CSM interaction (where it exists). Despite this, we decided not to remove the early data from the velocity evolution. If we do not take into account these observations, the discrepancy between the fitted models and the observed velocity evolution at early times could be large. This could result in incompatibilities as the interaction of the ejecta with a CSM produces a decrease in the photospheric velocities at early times (see, e.g. Fig. 13 of Rodríguez et al. 2020). Our velocity models should therefore be of the same order of magnitude as or display higher velocities than early observations.

Figures 5 and 6 show models drawn from the posterior distribution for each SN in our sample, together with median and 1σ confidence range for every parameter. The ejecta mass and the progenitor radius were interpolated linearly to the M_ZAMS we derived from the fitting. We also report the results in Table 2. It is important to note that the published errors on our estimated progenitor and explosion parameters are statistical in nature. The size of the errors indicates that our fitting technique is robust. However, these errors do not take into account systematic errors such as the uncertainties in stellar evolution modelling. There are a number of additional parameters which one could change in MESA that would give different pre-SN configurations (for the same M_ZAMS). A full exploration of these effects is beyond the scope of this work and will be the focus of future efforts. Here we use “standard” values for various stellar evolution parameters, which may not cover the full parameter space. As a consequence, the errors on the physical parameters are likely to be underestimated.

We can see good quality fits for the whole sample. However, we note some issues. First, it is important to note that SN 2018aoq presents only a few observations for the LC which includes the photospheric phase and the transition to the radioactive tail phase (see Fig. 3), and only one measurement of the Fe II 5169 Å velocity. Since there are no observations during the radioactive tail phase, the amount of ⁵⁶Ni synthesised cannot be derived (Hamuy 2003).

We were not able to obtain results for SN 2004et. Light curves drawn from the posterior distribution do not represent the observed LC. While the evolution of the photospheric velocity and the radioactive tail phase are well reproduced, the maximum a posteriori (MAP) model is fainter by an average of 0.27 dex during the photospheric phase. With the most recent estimation for the distance, the bolometric LC becomes remarkably luminous. Additionally, this SN presents an extensive plateau phase (∼125 days). The combination of both properties makes this SN an outlier in our sample and, unfortunately, we could not find any model in our grid that matches these properties (see Fig. A.12). As we mention in Sect. 4, we constrain the scale parameter to the uncertainties in the distance estimates. As with this restriction we could not find any set of parameters that fit the observations of this object, we tested whether or not solutions were possible in which the scale parameter is unconstrained. In this context we find a solution, but with a scale factor that represents a shorter distance to the SN which is ∼0.4 times the current distance estimate, or a combination of shorter distance and lower extinction. In both cases this distance is outside the range of allowed distances by the latest estimates. These results may indicate that the shorter distance is in fact correct or that our modelling cannot reproduce the observations because our range of initial models and explosion energies does not include any model compatible with this object. Recently, three studies arrive at similar distance estimates for its host galaxy (Murphy et al. 2018; Anand et al. 2018; Van Dyk et al. 2019). This may indicate a necessity to change some parameters in the evolutionary calculations to obtain a different progenitor structure that reproduces the properties of this SN. Therefore, we discard this solution for SN 2004et as it is not compatible with the new estimates of the distance. However, as we want to compare the results of this work using hydrodynamic simulations and MCMC fitting with those using other methods to look for compatibility, we modelled SN 2004et again but now assuming a shorter distance of 5.6 ± 0.4 Mpc (Smartt et al. 2009) as this is the previous estimate of the distance. The best fitting models for SN 2004et with this distance are shown in the second panel of Fig. 5.

Additionally, we find poor agreement in the radioactive tail phase for SN 2005cs. It seems that the decline rate of our model does not follow the observations. Pastorello et al. (2009) already pointed out that the decline rates of SN 2005cs in different bands are significantly smaller than the decline rate expected from the ⁵⁶Co decay and describe this phenomenon as a residual contribution from radiation energy, as first suggested by Utrobin (2007). According to Utrobin (2007), by the end of the optically thick phase the total radiation energy is not exhausted completely. A radiation flow generated in the warmer inner layers propagates throughout the optically thin layers and results in an additional source of energy. This transitional phase (labelled as plateau tail phase by Utrobin 2007) was also observed in SN 1999em and SN 1999eu.

Finally, we also find discrepancies in the velocities of SN 2005cs and SN 2012aw at early times. The largest difference found is about 1000 km s⁻¹.

5.1. Comparison with the analysis of pre-explosion imaging

Here we compare our estimation of the progenitor initial mass to those derived from the direct analysis of the progenitor star in pre-explosion images. We use the results from Smartt (2015) and Davies & Beasor (2018) for SN 2004A, SN 2004et, SN 2005cs, and SN 2012ec. It is important to note that these two papers use the old estimate for the distance to SN 2004et and so we can compare our results. For SN 2008bk we use a combination of the results from Mattila et al. (2008), Smartt et al. (2009), Van Dyk et al. (2012), and Davies & Beasor (2018). We discard the results from Maund et al. (2014a) as their mass determination is subject to a foreground extinction almost ten times higher than the later estimation made by Maund (2017). The solutions of Kochanek et al. (2012), Smartt (2015), and Davies & Beasor (2018) were used for SN 2012aw (see discussion in Smartt 2015). For SN 2017eaw we make use of the results from Van Dyk et al. (2019) and Eldridge & Xiao (2019). We do not consider the conclusions of Kilpatrick & Foley (2018) and Rui et al. (2019) with respect to SN 2017eaw as both studies assumed a shorter distance to the host galaxy (see Sect. 3.1). Finally, we use the mass estimation of O’Neill et al. (2019) for SN 2018aoq. In cases where more than one estimation of M_ZAMS is available, we take the complete range of values predicted in the different works instead of using some specific value. It should be noted that in almost all of these cases, the most accurate values are within the error bars of other solutions.

Figure 7 compares our results with those mentioned above. We find a good agreement between the masses estimated by both methods for almost every SN in the sample. SN 2004et and SN 2012ec are the only ones that escape from the trend. Our analysis suggests a low-mass progenitor for SN 2012ec while pre-SN imaging propose a more massive one, and the opposite for SN 2004et.

Fig. 7. Comparison between MZAMS obtained in this work using hydrodynamic simulations and MCMC methods and those based on the analysis of progenitors in pre-explosion images. For cases where multiple values of MZAMS from pre-SN imaging exist, we take the complete range of values predicted in the literature instead of using some specific value and its uncertainties. In these cases, the central point of the initial mass is the midpoint of the range.

Similar works to that presented here were published by Morozova et al. (2018) and Eldridge et al. (2019). These latter authors use large grids of hydrodynamic simulations and a non-Bayesian fitting procedure based on χ² minimisation. Morozova et al. (2018) use a large sample of SNe II which includes eight objects with observed progenitors. At first glance, there is a clear discrepancy between both quantities. They attribute it to the fact that they use pre-SN simulations from a different stellar evolution code to that used to connect the luminosity of the progenitor with its initial mass. Eldridge et al. (2019) only consider the SNe II with observed progenitors of Smartt (2015). Eldridge et al. (2019) claim that their results are consistent with those from pre-explosion imaging, although they note that their results have a tendency towards higher masses. In summary, both above papers present results that are consistent with those from pre-SN imaging but with a tendency to higher progenitor masses. It is important to mention that while both works use a larger sample of SNe II than we do here, some progenitor candidates have not yet been confirmed. Post-explosion images when the SN has faded sufficiently are needed in order to confirm the progenitor through its disappearance. Until this happens, results should be taken with caution as they can lead to incorrect determination of the progenitor star. Additionally, the mass discrepancy may be due to the fact that the authors only use the LC to obtain the progenitor properties without using any spectral information such as the expansion velocity. One of the differences with our work is that we fit the photospheric velocity simultaneously with the LC (see Sect. 6.1 for discussion).

In our work, the root-mean-square (RMS) deviation of our progenitor mass estimates against those from direct detections is found to be 2.8 M_⊙. We now compare with the RMS values from Morozova et al. (2018) and Eldridge et al. (2019) to quantify the differences. It must be stressed that when comparing with other works, we recalculate the RMS of our results taking into account only those SNe within both works. We find a RMS of 3.5 M_⊙ in our results and 5.8 M_⊙ in those by Morozova et al. (2018). For the results in Eldridge et al. (2019) we estimate a RMS value of 3.4 M_⊙, while 3.1 M_⊙ is found in our work for the same objects. In conclusion, we find that our mass estimations are more consistent with those from pre-explosion imaging than the estimations of previous works.

During the last two decades, tension has emerged between hydrodynamic modelling and pre-SN imaging in the sense that the progenitor mass estimated by the former was usually larger than the estimated or upper limits given by the direct analysis of progenitors. However, we do not find such discrepancy in our analysis.

5.2. Comparison with results from late-time spectral modelling

We have already compared our results with those which come from the detection and analysis of the progenitor star in pre-explosion images. Similar analysis can be performed by comparing with the results of progenitor M_ZAMS determined through nebular spectral modelling. Late-time spectra of SNe allow examination of the nucleosynthesis yields, especially the emission lines of [OI]λλ6300, 6364 as these lines characterise the core mass of the progenitor. Then, by comparison with synthetic nebular spectra available for different M_ZAMS, it is possible to distinguish between different progenitors (Jerkstrand et al. 2014).

We use the results from Jerkstrand et al. (2012, 2014, 2015, 2018) for SNe 2004et, 2012aw, 2012ec, and 2008bk, respectively, and Silverman et al. (2017) for SN 2004A, and Van Dyk et al. (2019) for SN 2017eaw to contrast with our values. A comparison of the progenitor masses obtained with both methods is shown in Fig. 8. Excellent agreement is found between the two methods for all objects with the exception of SN 2012ec. This is the only object that displays a different solution (see Sect. 6.3).

Fig. 8. Comparison between MZAMS obtained from hydrodynamic and late-time spectral modelling. Excellent agreement can be seen for all SNe II except SN 2012ec.

5.3. Ejecta masses

In this section we compare the ejecta masses derived in the present study with those presented by MB19. Even though both works use the same hydrodynamic code, the progenitor models and the selection of the preferred model was done in a very different way. In MB19, double polytropic models were used as pre-SN structures, while here we use stellar evolution calculations. Polytropic models allow to obtain the final structure of the star dependent only on the pre-SN properties, a priori without information on the initial mass. Therefore, we only provide a comparison to the ejecta masses instead of main-sequence masses. In addition, MB19 chose their preferred models by visual comparison while here we use a robust statistical method.

Figure 9 shows poor agreement between the ejecta masses from the two studies, with those of MB19 being systematically larger than ours. SN 2004et and SN 2012aw present the largest differences, being of the order of 3.5 M_⊙ and 10 M_⊙, respectively. As we mentioned above, the only physical difference between both works is in the calculation of the pre-SN models. Double polytropic calculations allow a large variety of pre-SN structures to be produced, each with different mass, radius, chemical composition, and density profile. The pre-SN models in MB19 may indicate that different solutions could be found if the standard assumptions in stellar evolution change, for example with respect to mixing processes and mass-loss rates. This is an option given the uncertainties still present in stellar modelling, especially in massive stars. Additionally, in this work we use a large grid of simulations and a fitting method with statistical support that initially covers the entire parameter space. While visual comparison can find solutions that reproduce the observations, it does not consider whether other solutions are possible, some of which can even be more probable.

Fig. 9. Ejecta masses estimated in this work in comparison with those from MB19.

5.4. Progenitor radii

Direct detections provide a unique opportunity to place the progenitor in an HR diagram. The observations determine the progenitor luminosity and effective temperature. Additionally, assuming a black body, the progenitor radius can be estimated. On the other hand, we can recover the final radius and other properties from the M_ZAMS derived from the hydrodynamical modelling as we use progenitor models from stellar evolution calculations.

In the following, we compare the progenitor radius from pre-SN imaging to our estimations. We use the values of the luminosity and effective temperature from Smartt (2015) to estimate the progenitor radius for SN 2004A, SN 2004et, SN 2005cs, SN 2008bk, SN 2012aw, and SN 2012ec, and O’Neill et al. (2019) for SN 2018aoq. Finally, we use the radius of the progenitor of SN 2017eaw from Van Dyk et al. (2019). Progenitor radii for each pre-SN model are listed in Table 1. The progenitor radius was interpolated linearly to the M_ZAMS we derived from the fitting.

Figure 10 shows a comparison of progenitor radii from our modelling with those derived from observed progenitor properties before explosion. We note that our results are consistent for five objects in our sample, while our modelling underpredicts the remaining three radii. In general, there is a tendency to yield lower progenitor final radius by an average of 175 R_⊙. Nevertheless, the agreement is quite good considering that the progenitor radius is not a free parameter in our modelling.

Fig. 10. Comparison between progenitor radii from our modelling with that derived from observed progenitor properties before explosion.

6. Discussion

6.1. Light-curve degeneracies

One of the main problems of inferring progenitor properties from LC modelling is that there is a degeneracy among some progenitor properties when reproducing the observations. Sometimes, similar photometric properties can be achieved with different progenitor and explosion characteristics. In Dessart & Hillier (2019), the authors pointed out that different progenitors can finish with a comparable H-rich envelope mass and produce similar photospheric phases. Moreover, Goldberg et al. (2019) argue that the ejecta mass, the explosion energy, and the progenitor radius cannot be constrained from the LC and velocity measurements. These latter authors solve it using two of these explosion properties as a function of the third. This requires independent knowledge of one parameter.

During recent years, hydrodynamic modelling of SNe II LCs and velocity measurements has suggested a significant discrepancy between the SN ejected masses and the initial masses of the observed progenitors. For example, Utrobin & Chugai (2008) present detailed modelling of SN 2005cs and a summary for another three SNe II (1987A, 1999em, and 2003Z). They argue that the hydrodynamic progenitor masses are systematically higher than if SNe II had originated from the range of 9–25 M_⊙, assuming a Salpeter initial mass function. This differs markedly from the direct detection of progenitors in pre-explosion images. However, MB19 analysed a sample of six SNe II with confirmed progenitors in post-explosion images and find that hydrodynamic masses are not systematically larger that those from pre-SN imaging. However, these latter authors do note that using similar pre-SN models and explosion parameters to Utrobin & Chugai (2008) and Utrobin & Chugai (2009), they arrive at similar LCs. Once again, this shows the high degree of degeneracy present in this problem. We consider that a detailed inspection of the degeneracy and the discrepancy between hydrodynamic masses and the masses inferred by the direct detection of the progenitors can be achieved with a large grid of hydrodynamic models in parameter space consisting of a considerable variety of LC and velocity models and a robust fitting procedure with statistical support, as presented in this work. With the above-mentioned considerations and from the analysis presented in Sect. 5.1 we conclude that we do not find such discrepancy between the progenitor initial masses inferred by hydrodynamical modelling and pre-SN imaging.

Furthermore, we also compare high-mass models with the observed LCs and photospheric velocities in order to look for possible differences in their appearance. We show the case of SN 2012ec as an example, but we performed the same analysis for the complete sample. Figure 11 shows hydrodynamic models for different initial masses (dashed lines) compared to the MAP model for SN 2012ec (solid line). In addition, models for several explosion energies are plotted. Every model has the same M_Ni and ⁵⁶Ni mixing as the MAP model. The top panel of Fig. 11 shows the case for a 15 M_⊙ progenitor. Despite the large range of explosion energies plotted we note that none of the models reproduce the LC. We do not show a larger range of energies since more energetic models will produce brighter and shorter plateau phases, and the opposite for models with lower explosion energy. The plateau luminosity and the photospheric velocity are well represented by the model with an explosion energy of 0.7 foe, but it predicts a longer plateau duration (∼40 days longer). Something similar is seen in the middle panel of Fig. 11 for the 20 M_⊙ model. The main difference is that in this example the model with 0.7 foe of energy also fails to reproduce the expansion velocities at early times. The bottom panel shows models computed for a 25 M_⊙ progenitor. Here, the model with 1.1 foe reproduces the plateau length but fails to reproduce its shape. Additionally, this model predicts low expansion velocities which differ from the observations. We conclude that we are not missing high-mass solutions. Therefore, we feel confident that our fitting procedure finds the best solutions and, at least within the grid of models we are using, other models are much less probable.

Fig. 11. Hydrodynamic models for different MZAMS (dashed lines) compared to the MAP model (solid line) and observations of SN 2012ec (green dots). Dashed lines are colour-coded according to different explosion energy. The models plotted correspond to the same values of MNi and 56Ni mixing as the MAP model. The MAP model is reproduced with MZAMS = 9.9 M⊙, E = 0.57 foe, MNi = 0.034 M⊙, and 56Ni mixing = 28%. Left panels: bolometric LCs. Right panels: evolution of the photospheric velocity. Top panel: 15 M⊙ model. Middle panel: 20 M⊙ model. Bottom panel: 25 M⊙ model.

In the case of the 25 M_⊙ progenitor, the contribution of the observed velocities when discarding the models with explosion energies between 1.0 and 1.1 foe is considerable. In the following we show with a simple test how, in some cases, the photospheric velocity helps to solve the dichotomy between different solutions. For this purpose we ran the MCMC sampler for SN 2017eaw again but this time considering only the bolometric LC. Results can be seen in Fig. 12. As expected, we obtain good agreement with the LC while photospheric velocities are not well reproduced. In this case the estimated M_ZAMS is 21.2 M_⊙, that is ∼6 M_⊙ larger than that estimated through LC and photosperic velocity modelling.

Fig. 12. Bolometric LC (left panel) and photospheric velocity evolution (right panel) for SN 2017eaw in comparison with fifty hydrodynamical models randomly chosen from the posterior probability distribution. In this example we did not take into account the photospheric velocity for the fitting. We note the importance of fitting both observables simultaneously.

Clearly, fits to the LC alone are not a good method to estimate the physical properties of explosions. At least in some cases the photospheric velocity evolution is essential in breaking the degeneracy (see also Ricks & Dwarkadas 2019). However, this is the methodology used recently by Morozova et al. (2018) and Eldridge et al. (2019). In both of these latter papers, there is a tendency to higher progenitor masses with respect to the stellar masses based on direct analysis of progenitors in pre-explosion imaging. With this example we emphasise that care has to be taken when deriving progenitor properties without using any spectral information. If this observable is not taken into account, it could lead to incorrect determination of the mass and energy.

6.2. Limitations and caveats

The stellar evolution simulations presented in this work require a large number of assumptions. In the calculation of our progenitor models we assume non-rotating stars and standard values for the mixing-length parameter and overshooting. Dessart et al. (2013) explore how variations of these parameters affect the final structure of a SN II progenitor. Different values produce changes in the progenitor radius, the H-rich envelope mass, and the helium-core mass, among others, which significantly influence the LCs, although the initial mass is the same. Studying all the existing possibilities according to the different evolutionary parameters that can be used is difficult and beyond the scope of this work.

As we mention in Sect. 2, the hydrodynamic code we use assumes a radially symmetric flow and adopts LTE to describe the radiative transfer. Although these approximations might not be entirely correct, the approach appears to be good. The very extended and massive hydrogen envelopes that characterise SNe II are expected to smooth the asymmetries of the explosion mechanism which makes spherical symmetry a good approximation for the bulk of the ejecta, though the ⁵⁶Ni distribution is more likely to be in some preferred direction (Wongwathanarat et al. 2015). On the other hand, LTE assumes that radiation and matter are strongly coupled. This is not valid at shock breakout and during and after the transition phase to a completely recombined ejecta.

Another approximation in the code is found in the opacity calculation. The code uses opacity tables calculated assuming LTE and a medium at rest. These calculations underestimate the true line opacity when considering rapidly expanding envelopes where large velocity gradients are present (Karp et al. 1977). In addition, the effect of the non-thermal excitation or ionisation of electrons that are created by Compton scattering of γ-rays emitted by radioactive decay of ⁵⁶Ni and ⁵⁶Co is not included in the calculation of the opacity. Our assumptions in the calculation of the opacity lead to considerable underestimation of the true ionisation. To partially solve this underestimation, the code adopts a minimum value of the opacity sometimes referred to as the “opacity floor” (see more details in Bersten et al. 2011). This approach has been extensively used in the literature (see, e.g. Young 2004; Morozova et al. 2015).

The analytical study of Popov (1993) shows the dependence of the bolometric luminosity on the opacity as it also shows the subsequent dependence of the explosion energy, mass, and radius on the opacity. Qualitatively, a larger value of the opacity decreases the plateau luminosity while increasing the duration of the plateau. As a consequence, this leads to different progenitor and explosion parameters. From Eq. (27) of Popov (1993), if the bolometric luminosity is fixed, an increase in the opacity leads to a higher explosion energy and a lower mass. Thus, given the opacity has an important effect on the calculation of the bolometric LC models, this gives rise to a corresponding uncertainty in the progenitor and explosion parameters.

Finally, in this section we discuss the uncertainties on our results that arise from the assumed photospheric velocities from both the models and the observations. As mentioned in Sect. 3, our hydrodynamical modelling requires measurement of the ejecta photospheric velocity. One of the typical procedures used to estimate the photospheric velocity is the measurement of the velocity at maximum absorption of optically thin lines as it is assumed that these lines are formed near the photosphere (Leonard et al. 2002). Dessart & Hillier (2005) analysed several synthetic line velocities and determined that the Fe II 5169 Å line delivers high accuracy in matching the photospheric velocity. This assumption is extensively used in the literature. Consequently, we use this line velocity as photospheric velocity indicator. However, the results achieved by Dessart & Hillier (2005) are restricted to a minimum velocity of ∼4000 km s⁻¹. Since some of the Fe II velocities in our sample are below that limit, we now discuss and analyse how the use of different techniques for estimating this velocity may affect our results.

Jones et al. (2009) used synthetic spectra from Eastman et al. (1996) and Dessart & Hillier (2005), and found polynomial relations to convert the observed H_β velocities into photospheric velocities. Moreover, Jones et al. (2009) suggest that models from Eastman et al. (1996) predict more realistic line profiles in the SN ejecta than those of Dessart & Hillier (2005), and therefore should provide a better photospheric velocity estimation. For this reason, we use the relation found by Jones et al. (2009) with models from Eastman et al. (1996) to derive the photospheric velocity. Figure 13 compares this photospheric velocity with the observed Fe II velocities. Within the uncertainties both velocities are generally consistent for the full sample. However, there are some caveats: (a) despite the good agreement for SN 2008bk, only one photospheric velocity is available for comparison as only one measurement of the H_β velocity is within the range of validity of the polynomial relation, (b) H_β velocities are not available for SN 2018aoq, and (c) for SN 2004et, there are differences between Fe II and photospheric velocities of the order of 1000 km s⁻¹ at 30–40 days after explosion. These differences decrease with time. This comparison provides additional support to the use of Fe II velocities as photospheric velocities. However, an additionally caveat is that the photospheric position in our models can be different from those in atmospheric models. This is due to the differences in the opacities involved in determining this location. In our models, the photosphere is defined where the Rosseland mean optical depth is two-thirds, while in Eastman et al. (1996) the photosphere is located where the Thomson scattering optical depth is two-thirds.

Fig. 13. Comparison between observed Fe II velocities (filled markers, solid lines) and photospheric velocities (open markers, dashed lines) inferred from Hβ velocities through the polynomial relation from Jones et al. (2009) for each SN in our sample. Errors are not plotted for better visualisation.

Recent studies have used another approximation to model ejecta velocities. Instead of using observed line velocities as photospheric velocity indicators, these works calculate the Fe II line velocity in the Sobolev approximation where the Sobolev optical depth is equal to one (Paxton et al. 2018; Goldberg et al. 2019; Ricks & Dwarkadas 2019; Bostroem et al. 2019). However, the precise Sobolev optical depth where the line is formed is not known, and different values translate into different velocities (see Fig. 36 of Paxton et al. 2018). We compare our model photospheric velocity with the Fe II line velocity as defined by these authors to check how our results could be influenced by this issue. We use Eq. (53) of Paxton et al. (2018) to calculate the Sobolev optical depth for the Fe II line. In order to do so, the ionisation fraction of iron atoms is needed. This information is provided in a table as a function of density and temperature and was obtained from the public version of MESA. The expression for the Sobolev optical depth (Eq. (53) of Paxton et al. 2018) is valid in a homologously expanding atmosphere. Therefore, the Fe II velocities are calculated only for times later than 25 days after explosion.

Figure 14 compares models of photospheric velocity and Fe II velocity for a 12 M_⊙ progenitor and three values of the explosion energy: 0.5, 1.0, and 1.5 foe. Differences of about 1000 km s⁻¹ are found at early times. As time goes by, the differences decrease. Therefore, if this way of comparing model Fe II velocities with observations is precise, we expect some changes in our results. According to Fig. 14, our photospheric velocity models underestimate Fe II velocities. As the expansion velocities are mostly affected by the energy of the explosion, our results could overestimate the explosion energy. Changes in progenitor masses are also possible.

Fig. 14. Comparison of the photospheric velocity and the Fe II line velocity for a 12 M⊙ progenitor and different explosion energies.

The following analysis estimates how our results could change if observed Fe II velocities are fitted using these model Fe II velocities. We compute bolometric LC and Fe II velocity models for 12, 15, 18, and 20 M_⊙ progenitors, with explosion energies of 0.5, 0.8, 1.0, and 1.2 foe for each mass value. We then fit these LCs and Fe II velocity models using the grid of models described in Sect. 2, that is, using the photospheric velocity models. This analysis gives an estimation of the differences we can expect in our results if Fe II velocity models are used to compare with the observations. We find a tendency to yield larger ejecta masses and explosion energies by an average of ∼0.2 M_⊙ and ∼0.2 foe, respectively. Changes in our results would go in the opposite direction as photospheric velocities are slower than Fe II line velocities (Fig. 14). Therefore, how one defines model velocities for comparison to observations can be associated to a systematic error in the ejecta masses and explosion energies. However, while our results on individual SNe II would change moderately, this would not significantly alter our conclusions.

In conclusion to this analysis, we have discussed different procedures to estimate the ejecta photospheric velocity, as well as other techniques to model the ejecta velocities. While small differences in best-fit physical parameters emerge, it is not completely clear which model velocities one should use, as both have uncertainties. The differences found can introduce a possible small bias in our results, but it does not affect our conclusions.

6.3. SN 2012ec

In Sects. 5.1 and 5.2 we compare the progenitor masses we obtained from the hydrodynamic modelling with the estimations from the direct detection of the progenitor star in pre-explosion images and nebular spectral modelling recovering a strong agreement. However, SN 2012ec is the only one that shows different solutions in both cases. From the hydrodynamic modelling point of view, we conclude an ejecta mass of ∼8 M_⊙, related to a main-sequence star of ∼10 M_⊙. On the other hand, pre-SN imaging infers a luminous and red progenitor. It is worth emphasising that pre-SN imaging provides a luminosity range and not the mass. The luminosity is then converted to an initial mass after comparison with evolutionary tracks of single-star models of different initial masses that terminates in the RSG phase within that range in luminosity. This analysis shows a progenitor of 16 ± 5 M_⊙. We observe that despite being more massive than our value, both estimates are not statistically distinct, given the large uncertainty on the pre-explosion image analysis. Furthermore, nebular spectral modelling suggests a main-sequence mass range of 13–15 M_⊙ due to the core mass of the progenitor. Collating this information, we obtain a luminous progenitor with a core mass corresponding to those stars of ∼14 M_⊙ in the ZAMS, but with a typical ejecta mass of a single 10 M_⊙ star. A detailed analysis of the progenitor star of SN 2012ec is beyond the scope of this work. Nevertheless, we note that the above combination can be obtained if we assume an enhanced rotation or binarity for SN 2012ec, as proposed by Straniero et al. (2019). Rotation produces higher mass helium cores and lower mass H-rich envelopes. A binary system with a primary star of initial mass estimated by nebular spectral modelling that experiences mass transfer episodes could also explain the disparity. In addition, from analytical estimations and by performing population synthesis simulations, Zapartas et al. (2019) conclude that a significant fraction (from one-third to half) of SN II progenitors are expected to interact with a companion before exploding, which supports this idea.

7. Summary and conclusions

We calculated a large grid of hydrodynamic models applied to stellar evolution progenitors in order to study the nature of SNe II. Light-curve modelling can provide constraints on progenitor and explosion properties although there is not always a unique solution. Therefore, we developed a robust method to derive physical properties based on MCMC methods using the observed bolometric LC and the expansion velocity simultaneously.

We applied this method to the observations of a well-studied set of SNe II (SNe 2004A, 2004et, 2005cs, 2008bk, 2012aw, and 2012ec) in order to compare with previous results from the analysis of pre-SN imaging. Progenitor identification was confirmed for these SNe via its disappearance in post-explosion images. In addition, we also include SN 2017eaw and SN 2018aoq in the sample, as these are the last SNe to be discovered and analysed with this method. We find that our results are entirely consistent between the initial masses estimated by both methods for almost every SN in the sample. Moreover, some works have questioned the ability of the hydrodynamic modelling to recover progenitor and explosion parameters, in the sense that progenitor masses from hydrodynamic modelling are usually larger. With this analysis we discard such a discrepancy and find a robust method to recover the progenitor mass, among other progenitor and explosion properties. We note that future high-resolution observations of the explosion site of SN 2017eaw and SN 2018aoq will be required in order to confirm their progenitor candidates. An additional comparison between our progenitor mass estimations and those from nebular spectral modelling was also carried out, showing very good agreement between these methods.

From the proposed analysis we conclude that we have developed a robust method to infer progenitor and explosion properties of SN II progenitors which is in complete agreement with results derived from other methods. We are now confident in our method and are able to move on to analysing a larger sample of SNe II.

Acknowledgments

We thank the anonymous referee for providing constructive comments that improved the content of this paper. SGG acknowledges support by FCT under Project CRISP PTDC/FIS-AST-31546 and UIDB/00099/2020. L. M. acknowledges support by UNRN under Project PI2018-40B696. Software: emcee (Foreman-Mackey et al. 2013), corner.py (Foreman-Mackey 2016), matplotlib (Hunter 2007), MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), MesaScript (Wolf et al. 2017), jupyter (Kluyver et al. 2016), pandas (McKinney 2010).

References

Appendix A: Additional plots

A.1. Interpolation method

In this section, some examples of interpolated hydrodynamical models of bolometric LCs and photospheric velocities are presented. Figure A.1 shows interpolated models after varying M_ZAMS, E, M_Ni, and ⁵⁶Ni mixing.

Fig. A.1. Interpolated synthetic bolometric LCs (left panels) and photospheric velocities (right panels) after varying the progenitor initial mass between 9 and 25 M⊙ (panel a), the explosion energy between 0.1 and 1.5 foe (panel b), the 56Ni mass between 0.005 and 0.08 M⊙ (panel c), and the 56Ni mixing between the 20% and 80% of the final structure in mass coordinate (panel d). The values increase from purple to yellow. Eighteen models are shown for each parameter being varied. The parameters not being varied are fixed at an initial mass of 10 M⊙, explosion energy of 1.3 foe, MNi of 0.01 M⊙, and 50% of 56Ni mixing, with the exception of the panel showing the 56Ni mixing effect for which a larger MNi of 0.06 M⊙ is used to enable better visualisation of this effect.

A.2. Corner plots of the posterior probability distributions

Additional corner plots of the joint posterior probability distribution of the parameters are presented in Figs. A.2–A.9.

Fig. A.2. Same as in Fig. 4 but for SN 2004A.

Fig. A.3. Same as in Fig. 4 but for SN 2004et.

Fig. A.4. Same as in Fig. 4 but for SN 2005cs.

Fig. A.5. Same as in Fig. 4 but for SN 2008bk.

Fig. A.6. Same as in Fig. 4 but for SN 2012aw.

Fig. A.7. Same as in Fig. 4 but for SN 2017eaw.

Fig. A.8. Same as in Fig. 4 but for SN 2018aoq.

Fig. A.9. Corner plot of the joint posterior probability distribution of the parameters for SN 2012ec when the priors for texp and the scale are relaxed.

A.3. Autocorrelation and trace plots

Here, examples of autocorrelation plots are shown in Fig. A.10. These were performed using the autocorrelation_plot tool implemented in the Python library pandas (McKinney 2010). Additionally, Fig. A.11 shows the trace plots of the MCMC samples. In both cases, we use SN 2017eaw as an example.

Fig. A.10. Autocorrelation plots for twenty chains randomly chosen using SN 2017eaw as an example. Each panel shows the autocorrelation for a different parameter.

Fig. A.11. Trace plots of the parameters for twenty chains randomly chosen using SN 2017eaw as an example. These plots show the evolution of the chains over time.

A.4. SN 2004et

Figure A.12 shows models drawn from the posterior distribution of the parameters for SN 2004et when using a distance to the host galaxy of 7.73 ± 0.78 Mpc. The lack of agreement between models and observations is easily seen. The LC models present large discrepancies during the photospheric phase. The models are on average ∼0.27 dex fainter and evolve more rapidly.

Fig. A.12. Comparison between models and observations for SN 2004et assuming d = 7.73 ± 0.78 Mpc. We show 50 models randomly chosen from the posterior probability distribution. The lack of agreement during the photospheric phase is clearly seen. Left: bolometric LC. Right: evolution of the photospheric velocity. The grey shaded region shows the early data we removed from the fitting.

All Tables

All Figures

	Fig. 1. Density profiles for the pre-SN models used in this work. The 9, 10, and 11 M⊙ models were calculated up to the end of core carbon burning since the evolution to core collapse for these stars is computationally expensive.
In the text

	Fig. 2. Ejecta masses (top panel) and progenitor radii (bottom panel) for the pre-SN models as a function of the initial mass.
In the text

	Fig. 3. Bolometric LCs (top panel) and absolute V-band LCs (bottom panel) of SN 2017eaw and SN 2018aoq in comparison to well-studied SN 2004A, SN 2004et, SN 2005cs, SN 2008bk, SN 2012ec, and SN 2012aw. Absolute curves were computed using the distances and reddenings in Table 1 of MB19, except in the case of SN 2004et for which we use the new estimates of the distance (see text).
In the text

	Fig. 4. Corner plot of the joint posterior probability distribution of the parameters for SN 2012ec.
In the text

Fig. 5. Comparison between models (solid lines) and observations (filled dots) for our SN sample. We show 50 models randomly chosen from the posterior probability distribution. Left: bolometric LC. Right: evolution of the photospheric velocity. From top to bottom: SN 2004A, SN 2004et, SN 2005cs, and SN 2008bk. The grey shaded region shows the early data we removed from the fitting. For SN 2004et we show the results using d = 5.9 ± 0.4 Mpc to calculate the bolometric LC.

In the text

	Fig. 6. Comparison between models (solid lines) and observations (filled dots) for our SN sample. We show 50 models randomly chosen from the posterior probability distribution. Left: bolometric LC. Right: evolution of the photospheric velocity. From top to bottom: SN 2012aw, SN 2012ec, SN 2017eaw, and SN 2018aoq. The grey shaded region shows the early data we removed from the fitting.
In the text

Fig. 7. Comparison between MZAMS obtained in this work using hydrodynamic simulations and MCMC methods and those based on the analysis of progenitors in pre-explosion images. For cases where multiple values of MZAMS from pre-SN imaging exist, we take the complete range of values predicted in the literature instead of using some specific value and its uncertainties. In these cases, the central point of the initial mass is the midpoint of the range.

In the text

	Fig. 8. Comparison between MZAMS obtained from hydrodynamic and late-time spectral modelling. Excellent agreement can be seen for all SNe II except SN 2012ec.
In the text

	Fig. 9. Ejecta masses estimated in this work in comparison with those from MB19.
In the text

	Fig. 10. Comparison between progenitor radii from our modelling with that derived from observed progenitor properties before explosion.
In the text

Fig. 11. Hydrodynamic models for different MZAMS (dashed lines) compared to the MAP model (solid line) and observations of SN 2012ec (green dots). Dashed lines are colour-coded according to different explosion energy. The models plotted correspond to the same values of MNi and 56Ni mixing as the MAP model. The MAP model is reproduced with MZAMS = 9.9 M⊙, E = 0.57 foe, MNi = 0.034 M⊙, and 56Ni mixing = 28%. Left panels: bolometric LCs. Right panels: evolution of the photospheric velocity. Top panel: 15 M⊙ model. Middle panel: 20 M⊙ model. Bottom panel: 25 M⊙ model.

In the text

	Fig. 12. Bolometric LC (left panel) and photospheric velocity evolution (right panel) for SN 2017eaw in comparison with fifty hydrodynamical models randomly chosen from the posterior probability distribution. In this example we did not take into account the photospheric velocity for the fitting. We note the importance of fitting both observables simultaneously.
In the text

	Fig. 13. Comparison between observed Fe II velocities (filled markers, solid lines) and photospheric velocities (open markers, dashed lines) inferred from Hβ velocities through the polynomial relation from Jones et al. (2009) for each SN in our sample. Errors are not plotted for better visualisation.
In the text

	Fig. 14. Comparison of the photospheric velocity and the Fe II line velocity for a 12 M⊙ progenitor and different explosion energies.
In the text

Fig. A.1. Interpolated synthetic bolometric LCs (left panels) and photospheric velocities (right panels) after varying the progenitor initial mass between 9 and 25 M⊙ (panel a), the explosion energy between 0.1 and 1.5 foe (panel b), the 56Ni mass between 0.005 and 0.08 M⊙ (panel c), and the 56Ni mixing between the 20% and 80% of the final structure in mass coordinate (panel d). The values increase from purple to yellow. Eighteen models are shown for each parameter being varied. The parameters not being varied are fixed at an initial mass of 10 M⊙, explosion energy of 1.3 foe, MNi of 0.01 M⊙, and 50% of 56Ni mixing, with the exception of the panel showing the 56Ni mixing effect for which a larger MNi of 0.06 M⊙ is used to enable better visualisation of this effect.

In the text

	Fig. A.2. Same as in Fig. 4 but for SN 2004A.
In the text

	Fig. A.3. Same as in Fig. 4 but for SN 2004et.
In the text

	Fig. A.4. Same as in Fig. 4 but for SN 2005cs.
In the text

	Fig. A.5. Same as in Fig. 4 but for SN 2008bk.
In the text

	Fig. A.6. Same as in Fig. 4 but for SN 2012aw.
In the text

	Fig. A.7. Same as in Fig. 4 but for SN 2017eaw.
In the text

	Fig. A.8. Same as in Fig. 4 but for SN 2018aoq.
In the text

	Fig. A.9. Corner plot of the joint posterior probability distribution of the parameters for SN 2012ec when the priors for texp and the scale are relaxed.
In the text

	Fig. A.10. Autocorrelation plots for twenty chains randomly chosen using SN 2017eaw as an example. Each panel shows the autocorrelation for a different parameter.
In the text

	Fig. A.11. Trace plots of the parameters for twenty chains randomly chosen using SN 2017eaw as an example. These plots show the evolution of the chains over time.
In the text

	Fig. A.12. Comparison between models and observations for SN 2004et assuming d = 7.73 ± 0.78 Mpc. We show 50 models randomly chosen from the posterior probability distribution. The lack of agreement during the photospheric phase is clearly seen. Left: bolometric LC. Right: evolution of the photospheric velocity. The grey shaded region shows the early data we removed from the fitting.
In the text