© ESO 2022

1. Introduction

Core-collapse supernovae (CCSNe) are produced by the explosion of massive stars (> 8 − 10 M_⊙). The collapse is initiated when the iron core achieves the Chandrasekhar mass (Chandrasekhar 1939). When the collapsing core reaches nuclear densities, the implosion rebounds, generating a shock wave that stalls into accretion due to the flat density profile around the iron core. However, neutrino heating from the proto-neutron star revives the shock and produces the supernova (SN) explosion. This is the so-called delayed neutrino-driven mechanism, and it is the most favoured explosion mechanism for CCSNe (see e.g. Burrows & Vartanyan 2021, for a recent review).

Type II supernovae (SNe II) are the most common type of CCSN in nature (Li et al. 2011; Shivvers et al. 2017). In the context of single stellar evolution, they are believed to arise from the least massive stars among the core-collapse mechanism. Type II SNe are classified by the presence of hydrogen lines in their optical spectra (Minkowski 1941). Previous theoretical studies have shown that an extensive hydrogen-rich envelope, typical of red supergiant (RSG) stars, is required to reproduce SN II light-curve (LC) morphologies (e.g. Grassberg et al. 1971; Chevalier 1976). More recently, this picture was confirmed by the detection of a significant number of RSG stars at the position of SN II explosion sites on archival images prior to explosion, demonstrating that they are the progenitors of most SNe II (e.g. Van Dyk et al. 2003; Smartt 2009). Type II SNe were historically grouped into SNe IIP and SNe IIL categories based on the shape of their LCs (Barbon et al. 1979). However, in this paper we use ‘SNe II’ to refer to both groups together since recent studies indicate that they come from a continuous population (Anderson et al. 2014a; Sanders et al. 2015; Galbany et al. 2016; Valenti et al. 2016; Rubin & Gal-Yam 2016; de Jaeger et al. 2019, although see Davis et al. 2019 for distinct populations in near-infrared spectral features). Other hydrogen-rich SNe (SNe IIb, SNe IIn, and SN 1987A-like events) are not analysed in this study, and they will not be discussed.

During the last few decades, several works have focused on the analysis of large samples of SNe II to examine their photometric and spectroscopic diversity (e.g. Patat et al. 1994; Hamuy 2003; Bersten & Hamuy 2009; Arcavi et al. 2012; Anderson et al. 2014a,b; Faran et al. 2014a,b; Spiro et al. 2014; González-Gaitán et al. 2015; Gutiérrez et al. 2014, 2017a,b; de Jaeger et al. 2018; Davis et al. 2019). The great diversity observed in SNe II is indicative of a large variety of progenitor and explosion properties. A key aim of SN research is to determine the full range of parameters that produce these events and constrain the predominant physical properties that yield the observed diversity. Obtaining such knowledge is critical to furthering our understanding of how massive stars evolve through their lives to produce hydrogen-rich events, along with the properties of the explosion mechanism.

There are different approaches to connecting the characteristics of SNe with their progenitors (e.g. pre-SN imaging, LC modelling, and nebular-phase spectral modelling), but here we concentrate on those methods that use LCs and spectral information to extract the physical properties of SNe II. Since the pioneering works of Litvinova & Nadezhin (1983, 1985), who presented a set of analytical relations connecting SN observables (the absolute V magnitude and photospheric velocity at mid-plateau and the plateau duration) with progenitor and explosion properties (ejecta mass, progenitor radius, and explosion energy), significant efforts have been made to improve the estimation of confident physical parameters. While the concept of using simple relations to derive progenitor properties is in principle powerful, atypical values have sometimes been obtained (see e.g. Hamuy 2003), motivating advances in this methodology. The hydrodynamical models utilised to calibrate these analytical relations were based on simplified physical assumptions (e.g. simple polytropic models as initial configurations, old opacity tables, and the omission of ⁵⁶Ni heating in the calculations). Similar relations were presented in Popov (1993), although they were built from a two-zone analytical model for hydrogen recombination in SN II ejecta (see also Arnett 1980, for an analytic analysis for SN II LCs). New analytical relations based on modern hydrodynamical models (including the effects of ⁵⁶Ni heating and contemporary opacity tables) were presented by Kasen & Woosley (2009) and more recently by Goldberg et al. (2019). While these methods enable the user to extract physical parameters of SNe II in a relatively simple manner, it is not clear how accurate such constraints are because they are obtained using only three observables.

An alternative method for deriving the physical parameters of SNe II consists of detailed modelling of the complete LC evolution – sometimes modelled together with spectral information. In such techniques, the explosion is often simulated by artificially adding internal energy near the centre of the star, known as a ‘thermal bomb’. Here, the structure of the progenitor prior to core collapse has to be assumed, together with the explosion energy, ⁵⁶Ni yields, and ⁵⁶Ni distribution within the ejecta (e.g. Utrobin 2007; Bersten et al. 2011; Morozova et al. 2015). These parameters are freely chosen to reproduce the observations. Most studies use hydrodynamical codes that solve the radiation transport in the diffusion approximation, thus producing bolometric LCs. More sophisticated codes incorporate radiation-hydrodynamical modelling with LC and spectral information (Pumo & Zampieri 2011) or multi-group radiative transfer that produces multi-band LC simulations (Blinnikov et al. 1998). Additionally, radiative transfer codes have been developed to calculate spectra in rapidly expanding SN atmospheres, assuming homologous expansion (e.g. Dessart & Hillier 2005; Kasen et al. 2006). Further sophistication can be attained by calculating more realistic explosions using CCSN simulations, parameterising the delayed neutrino-driven mechanism, and coupling this to explosive nucleosynthesis calculations. In this context, the explosion energy, ⁵⁶Ni mass, and distributions are no longer free parameters. These calculations simulate the collapse and explosion, but they need alternative codes to reproduce SN observables. Recent efforts have been made in this direction (Sukhbold et al. 2016; Curtis et al. 2021; Barker et al. 2021), where multi-dimensional explosions are sometimes mimicked in 1D models and mapped to different codes to further obtain SN observables. While this methodology is more consistent (given that it only needs the progenitor structure as input), sometimes it cannot reproduce the observations and is restricted to a limited parameter space.

Progenitor models can also be constructed in different ways, either as an ad hoc configuration or with detailed stellar evolution from main sequence to core collapse. In general, non-rotating single-star models are used, whereas binarity or stellar evolution with enhanced mass loss could have a significant effect on the structure of the progenitor. Stellar evolution produces a final structure for each star, for which the pre-explosion mass and the progenitor size are not independent. To treat the progenitor final mass and radius as independent parameters, non-evolutionary models – polytropic models – are used (e.g. Utrobin 2007; Bersten et al. 2011). These models can also mimic Rayleigh-Taylor (RT) mixing during shock propagation by smoothing the transitions between zones of different chemical abundances (Utrobin & Chugai 2008). It is important to note that different approaches assumed in the construction of pre-SN models and the subsequent explosion and transient modelling may produce different results because of their assumed physics and the assumed free parameters.

Many studies have been published that show hydrodynamical modelling of individual SNe II. Recently, the number of multiple studies of SNe II has increased (e.g. Pumo et al. 2017; Ricks & Dwarkadas 2019; Utrobin & Chugai 2019; Eldridge et al. 2019; Martinez & Bersten 2019; Martinez et al. 2020). However, there are still only a few analyses of statistically significant samples (González-Gaitán et al. 2015; Morozova et al. 2018; Förster et al. 2018). In Martinez et al. (2020, hereafter M20), we analysed a sample of eight SNe II with observed progenitors by fitting their bolometric LC and photospheric velocity evolution. For this purpose, we presented a large grid of explosion models and a fitting procedure based on Markov chain Monte Carlo (MCMC) methods. We found that our results are consistent with those from pre-SN imaging and late-time spectral modelling. The analysis by M20 served as the basis for the present study, where we derive physical parameters for a large sample of SNe II from a homogeneous dataset.

The present work is the second in a series of three papers based on the study of a sample of 74 SNe II observed by the Carnegie Supernova Project-I (CSP-I; Hamuy et al. 2006). CSP-I was a SN follow-up programme based at the Las Campanas Observatory that obtained high-quality optical and near-infrared (NIR) LCs and optical spectra with high observational cadence. In Martinez et al. (2022a, hereafter Paper I), we calculated the bolometric LCs for this sample of SNe II, discussing our methodology in detail and presenting an analysis of the observed bolometric parameters. An additional conclusion of Paper I is the importance of having NIR observations to derive reliable bolometric luminosities. Here, in Paper II, we assess the progenitor and explosion properties for this sample of events and characterise the physical parameter distributions. Finally, in Martinez et al. (2022b, hereafter Paper III), we study correlations between physical parameters and different LC and spectroscopic measurements, and analyse SN II diversity, tying this to the physics of massive star explosions and their progenitors.

The current paper is organised as follows. In the following section we briefly outline the physical processes involved in SN II explosions. Section 3 describes the data sample. In Sect. 4 we present the hydrodynamical simulations and the fitting procedure used to derive the physical parameters of the sample. Section 5 presents the distributions of the physical parameters for our sample of SNe II. In Sect. 6 we compare our findings with those presented in previous works. In Sect. 7 we discuss possible explanations for the initial mass distribution found. We provide our concluding remarks in Sect. 8. Further analysis and figures not included in the main body of the manuscript are presented in the appendices.

2. SN II physics

The observed behaviour of SN II LCs and expansion velocities is related to physical properties of the progenitor star and the explosion; such as the ejecta mass (M_ej), the progenitor radius (R), the explosion energy (E), the amount of ⁵⁶Ni (M_Ni), and the distribution of ⁵⁶Ni in the ejecta (⁵⁶Ni mixing). Below, different phases of SNe II evolution are discussed in the context of their physical processes. Detailed reviews can be found in the literature (e.g. Arnett 1996; Woosley et al. 2002; Zampieri 2017).

During the core collapse, a powerful shock wave is generated near the centre of the star. The shock wave heats and accelerates the progenitor envelope until it arrives at the stellar surface, where photons begin to diffuse outwards (a phase known as the ‘shock breakout’). After the shock breakout, the outermost layers of the ejecta expand at great speed leading to a fast decrease in the temperature and, therefore, causing a rapid decline in the bolometric luminosity. This phase is commonly known as the ‘cooling phase’, and it is directly related to the progenitor size. However, this phase is also closely related to any material near the stellar surface.

Type II SNe enter the recombination phase when the temperature drops to ∼6000 K. This phase is commonly referred to in the literature as the ‘plateau’, although it is not necessarily a phase of constant luminosity. During this phase, hydrogen recombination takes place at different layers of the ejecta as a recombination wave recedes (in mass coordinate) through the expanding ejecta (e.g. Grassberg et al. 1971; Bersten et al. 2011). Therefore, the hydrogen-rich envelope mass of the progenitor star at the time of explosion (M_H, env) is related to the duration of the plateau phase, although other physical parameters also play a role (see Paper III). In addition, the explosion energy significantly drives SN II observational properties (Kasen & Woosley 2009; Dessart et al. 2013; Bersten 2013). More energetic explosions produce more luminous SNe II expanding at higher velocities, thus cooling and recombining the ejecta more rapidly, and producing shorter plateau phases. On the other hand, more extended progenitors produce more luminous and longer plateau phases. The additional heating of the SN II ejecta at late times by the ⁵⁶Ni decay chain extends the duration of the plateau (Kasen & Woosley 2009) and increases the luminosity in the late-plateau phase (Bersten 2013). The mixing of ⁵⁶Ni determines when energy deposition from radioactive decay starts to influence the LC, affecting the duration of the plateau and its shape (Bersten et al. 2011; Kozyreva et al. 2019). An extensively mixed ⁵⁶Ni impacts the LC earlier (see Bersten et al. 2011, their Fig. 12), that is, ⁵⁶Ni powers the LC sooner, producing a slowly declining plateau (typical SN IIP). The cooling and plateau phases together are also known as the optically thick or photospheric phase.

When the hydrogen-rich ejecta is recombined, SNe II enter a transition phase, which is marked by a rapid decline in luminosity. After this, the luminosity is mainly powered by the ⁵⁶Co decay. This phase is known as the radioactive tail and its luminosity primarily depends on the amount of ⁵⁶Ni in the ejecta.

3. Supernova sample

The sample of SNe II used in this study is the same as that analysed in Paper I where we present bolometric LCs for 74 SNe II observed by the CSP-I (Hamuy et al. 2006, PIs: Phillips and Hamuy; 2004−2009). The sample is characterised by a high cadence and quality of the observations. In addition, CSP-I LCs cover a wide wavelength range from optical (uBgVri) to NIR (YJH) bands. The data were homogeneously observed and processed (see Hamuy et al. 2006; Contreras et al. 2010; Stritzinger et al. 2011; Folatelli et al. 2013; Krisciunas et al. 2017, for a detailed description of the data reduction and photometric calibration). CSP-I SN II photometry is presented in Anderson et al. (in prep.), while the optical spectra were published by Gutiérrez et al. (2017a,b).

CSP-I was a SN follow-up programme that generally obtained observations for any SN that (a) was sufficiently bright to observe for a number of weeks; (b) had reasonable sky visibility to enable such observations; and (c) where the classification spectrum indicated a reasonably young explosion. This resulted in a magnitude-limited sample of SNe II, which we use in the current study. Thus our sample will be biased towards the inclusion of intrinsically brighter SNe II. We discuss how this may affect our results and conclusions in Sect. 7.2.

In the present work, we use the bolometric LCs from Paper I and the Fe II 5169 Å line velocities from Gutiérrez et al. (2017a) to infer the physical properties of the SNe II in the sample through hydrodynamic modelling (Sect. 5). The explosion epochs are taken from Gutiérrez et al. (2017b, see also Table 1 from Paper I for details).

4. Methods

4.1. Hydrodynamical simulations

The determination of the physical properties of SNe II is based on comparing models constructed with different physical parameters with observations. In the current work, we used the grid of hydrodynamic models presented in M20. These models were calculated using a 1D Lagrangian code that simulates the explosion of the SN and produces bolometric LCs and expansion velocities at the photospheric layers (see Bersten et al. 2011, for details). The grid of explosion models comprises a wide range of parameters (Table 1 and M20, for details). Pre-SN models in hydrostatic equilibrium are necessary to initialise the explosion. The public stellar evolution code MESA¹ (Paxton et al. 2011, 2013, 2015, 2018, 2019) was used to obtain non-rotating solar-metallicity pre-SN RSG models with initial masses (M_ZAMS) in the 9−25 M_⊙ range². Massive-star models are affected by the uncertainties in stellar modelling, such as the treatment of convection, rotation, and mass loss. In this context, we utilised the prescriptions adopted by the community to model these physical processes, and adopted standard values for the parameters involved to test standard single-star evolution against SN II LC and velocity modelling.

Following the above, for convection we adopted the Ledoux criterion and a mixing-length parameter of α_mlt = 2.0. We used exponential overshooting parameters of f_ov = 0.004 and f_ov, D = 0.001, a semiconvection efficiency of α_sc = 0.01, and thermohaline mixing with the coefficient α_th = 2 according to Farmer et al. (2016). We adopted the ‘Dutch’ prescriptions for the wind mass loss defined in the MESA code with an efficiency η = 1. For each value of M_ZAMS (9−25 M_⊙ with 1 M_⊙ increment), there is a corresponding value of M_ej and R. In this grid of simulations, M_ej cover a range of 7.9−15.7 M_⊙, while R is found in the range of 445−1085 R_⊙, similar to the values derived for RSGs (Levesque et al. 2005).

4.2. Fitting procedure

A fitting procedure based on MCMC methods was employed to find the posterior probability of the model parameters given the observations. This technique was first used for SNe in Förster et al. (2018). A similar procedure to that used in the current work was presented in M20. The MCMC is implemented via the emcee package (Foreman-Mackey et al. 2013). Since the models presented in M20 may not be sufficient when fitting SN observations using statistical inference techniques, we used the interpolation method presented in Förster et al. (2018). This is a robust and quick method that allows for irregular grids of models in the space of parameters to be used. For the MCMC runs, we used 400 parallel samplers (or walkers) and 10 000 steps per sampler, with a burn-in period of 5000 steps. The walkers were randomly initialised covering the entire parameter space (see Sect. 4 of M20, for details).

As in M20, there are six parameters in our model. Four of them are physical parameters of the explosions and their progenitors: M_ZAMS, E, M_Ni, and ⁵⁶Ni mixing. We included two additional parameters: the explosion epoch (t_exp) and a parameter named ‘scale’. The scale parameter multiplies the bolometric luminosity by a constant dimensionless factor to allow for errors in the bolometric LC because of the uncertainties in the distance and host-galaxy extinction. We used uniform priors for the following parameters: t_exp, M_ZAMS, E, M_Ni, and ⁵⁶Ni mixing. We allowed the sampler to run within the observational uncertainty of t_exp (see Table 1 from Paper I) and within the ranges of the physical parameters in our grid of explosion models (Table 1). We used a Gaussian prior for scale, where the Gaussian is centred at one. The uncertainty in the distance allows the SN to yield a more or less luminous bolometric LC. For this effect, we allowed variations to the bolometric LC produced by scale values within ±1σ_d, where σ_d is the relative difference in luminosity because of the uncertainty in the distance. The effect of unaccounted host-galaxy extinction allows the SN to be only intrinsically more luminous, making the scale prior asymmetric. To mimic this effect, we let the scale vary up to +1σ_AVhost (i.e. only in the direction that would make the SN brighter), where σ_AVhost is the relative difference in bolometric luminosity produced by taking a host-galaxy extinction equal to two standard deviations of the measurements of the host extinction from Anderson et al. (2014a)³. To summarise, the Gaussian prior for the scale is asymmetric. If the scale takes values below unity, the standard deviation is equal to σ_d, and σ_d + σ_AVhost everywhere else.

Extinction also changes the shape of the bolometric LC⁴. The scale prior does not consider different bolometric LC shapes. However, we tested how much the shape of the bolometric LC changes when different host-galaxy extinctions are considered and found that the shape is only significantly affected at early times (≲30 days after explosion) and for large extinction values. Given that we removed the first 30 days of evolution from the fitting (see Sect. 4.3), this effect is not considered in our prior.

4.3. Modelling background

We used the models and the fitting procedure described in the previous subsections to derive the physical properties of the SNe II in our sample. Expansion velocities are important as they provide restrictions to the ejecta expansion rate. In M20, we showed that in some cases, fits to the expansion velocities are crucial in breaking the LC degeneracies. In addition, the physical parameters derived from fits to the LC alone are significantly different to those that also include velocity data.

The bolometric LCs from the CSP-I sample were already presented in Paper I. The photospheric velocity is typically estimated by measuring the velocity at maximum absorption of optically thin lines (Leonard et al. 2002). Consequently, we used the Fe II 5169 Å line – measured by Gutiérrez et al. (2017b) – as a photospheric velocity indicator. This assumption is extensively used in the literature given that Dessart & Hillier (2005) showed that Fe II line delivers high accuracy in reproducing the photospheric velocity. However, these results are restricted to a minimum velocity of ∼4000 km s⁻¹. Recently, Paxton et al. (2018) proposed an alternative approximation to model the ejecta velocities. They calculated the Fe II line velocity in the Sobolev approximation and compared with the observed Fe II velocities. M20 analysed the variation on the physical parameters if different model velocities are used to compare to observations. They found a tendency to larger M_ej and E values of, on average, ∼0.2 M_⊙ and ∼0.2 foe, respectively. While differences exist, they do not significantly alter the results. The reader is referred to Sect. 5.2 of M20 for details.

Recently, it has been proposed that some SNe II show signatures of interaction with a dense circumstellar material (CSM) shell surrounding the star (e.g. González-Gaitán et al. 2015; Khazov et al. 2016; Yaron et al. 2017; Förster et al. 2018; Bruch et al. 2021). The interaction of the ejecta with a low-mass CSM is thought to only significantly affect the early evolution, with negligible effect at later epochs where the evolution is dominated by the hydrogen recombination and radioactive decay (Morozova et al. 2018; Hillier & Dessart 2019). Given that we focus on inferring the intrinsic properties of SN II progenitors instead of characterising the CSM properties, our explosion models were calculated without including CSM. Therefore, our fitting is not valid during the early evolution of SNe II. For this reason, we did not consider the first 30 days of evolution of the observed LC in our fitting procedure. This value is similar to the median of the distribution of the cooling phase duration determined in Paper I. Differences between our models and observations are therefore to be expected during the cooling phase, which may be strongly affected if CSM interaction exists. The results of LC and ejecta velocity modelling are presented in Sect. 5, but here we note that around 60% of the SNe II in the CSP-I sample with photometric observations earlier than 30 days after explosion show differences in the early bolometric LC. Those SNe II that show differences between our models and observations during the first 30 days are good candidates to study the general properties of the CSM. In spite of this, we did not remove the early data from the velocity evolution for two reasons. On the one hand, the effect of moderate CSM on the velocity evolution, as expected for normal SNe II, is less important than for the LC (Englert Urrutia et al. 2020). On the other hand, its effect is to reduce the velocities only at early times (Moriya et al. 2018). When modelling SNe II with interaction signatures using models that do not account for CSM interaction, the models should display higher or similar early velocities than observed, given that these models may serve as the basis for future models considering interaction with a CSM. Then, these new models would show lower early velocities than our models, possibly consistent with observations. If early velocities are not taken into account, our best-fitting model could display higher or lower velocities than early-time observations. From the above discussion, the former case is not an issue. However, the latter case would be incompatible within the above scenario, since any future models including CSM interaction would have early velocities even lower than observed. Therefore, to avoid this incompatibility, we decided not to remove the early velocities from the fitting.

5. Results

We searched for the probability distributions of the parameters for each SN II in the sample employing the set of explosion models and the MCMC procedure presented in Sect. 4. We define a ‘gold sample’ of our SNe II, selecting those events that followed our selection criteria: (a) bolometric LCs covering the photospheric phase and at least the beginning of the transition to the radioactive tail phase, since the latter is crucial to constraining principally the M_ej and E; (b) at least two Fe II velocity measurements during the photospheric phase and separated in time by more than ten days; and (c) the maximum a posteriori model reproduces the observations⁵. The assessment of the quality of the fits to the models was achieved visually and independently by the first three authors of this paper (LM, MB, and JA). A small number of contentious cases were discussed; however, their inclusion or exclusion in the gold sample does not affect our results or conclusions (see Appendix A). In total, 24 SNe II fall into the gold sample (see Table B.1).

The full sample of SNe II comprises 53 events and is formed by the SNe in the gold sample plus 29 additional objects for which progenitor and/or explosion parameters were determined, but insufficient data coverage or fitting quality prevent these SNe II from being considered gold events. SN 2008bk was considered a gold event even though the conditions were not fulfilled with our dataset since there are no observations during the transition from the plateau to the tail phase. SN 2008bk is a well-studied SN and for this reason, we used already published data to restrict the end of the plateau (Van Dyk et al. 2012).

Unfortunately, the following SNe do not fall into either of the previous groups and were excluded from the rest of this work. The models in our set of explosions cannot reproduce the observed behaviour of five SNe II⁶: 2006Y, 2008bu, 2008bm, 2009aj, and 2009au. The first two SNe show atypical optically thick phase durations (‘optd’⁷) of only 64 ± 4 days and 52 ± 7 days for SN 2006Y and SN 2008bu, respectively (Paper I). No model within our grid presents such a short optd (see also Sect. 6.7). We are able to model short-plateau SNe II by increasing the mass loss during the evolution of their progenitors, which reduces the extent of the hydrogen-rich envelope at the time of collapse. These results are presented in Paper III (see also Hiramatsu et al. 2021, for short-plateau SN II modelling). SNe 2008bm, 2009aj, and 2009au were already analysed in Rodríguez et al. (2020). Those authors showed that these events share the following common characteristics: low expansion velocities, absolute V-band LCs much brighter compared to normal SNe II with such velocities, and signs of interaction of the ejecta with CSM, among others. Moreover, based on hydrodynamic simulations, Rodríguez et al. (2020) found that a massive CSM of ∼3.6 M_⊙ was needed to reproduce the entire LCs and expansion velocities of SN 2009aj. Such a large CSM mass is expected to influence SNe II properties at epochs much later than the 30 day limit we assumed for the rest of the sample. As we mentioned in Sect. 4.3, our explosion models were calculated without including any CSM. Therefore, it is understandable that we are not able to find a set of parameters that can represent the full observations of these three events.

The explosion epoch was not estimated for SN 2005gk, SN 2005hd, or SN 2005kh, which makes it impossible to infer reliable results. Finally, a number of SNe II have insufficient data for constraining their physical properties from LC modelling (SNe 2004dy, 2005K, 2005es, 2006bc, 2006it, 2006ms, 2008F, 2008bh, 2008bp, 2008hg, 2008ho, 2008il, and 2009A) and were also excluded from this work.

Figure 1 gives an example of our fits for one specific SN (SN 2006ai), where models drawn from the posterior distribution of the parameters are shown in comparison with the observed bolometric LC (top panel) and Fe II line velocities (bottom panel). Additionally, Appendix B compares models and observations for the entire CSP-I SN II sample, and shows an example corner plot with the posterior probability distribution of the parameters of SN 2006ai. We characterise the results of the sample by the median of the marginal distributions. The 16th and 84th percentiles measure the width of the distribution of the sample, which are adopted as the lower and upper uncertainties. The results are reported in Table B.1 where we also include estimates of M_ej, M_H, env, the total mass of hydrogen at the time of collapse (M_H), and R. We emphasise that these progenitor properties are not model parameters and, therefore, they were not fitted. These values were linearly interpolated from the M_ZAMS derived from the fitting. We note that most of our estimations have small error bars. This arises from the high cadence and quality of the observations. However, these errors do not take into account systematics such as the uncertainties in the hydrodynamic simulations and stellar evolution modelling. In the latter, we used ‘standard’ values for various parameters (wind efficiency, convection, metallicity). Changes in these parameters give different pre-SN configurations, and therefore, different results. As a consequence, the errors on the physical parameters are likely to be underestimated. We reiterate that the results presented in the following sections are achieved by using standard pre-SN models (see details in Sect. 4.1), similar to those used by the studies that determine initial masses from progenitor detection in pre-explosion images. While this brings several caveats (that are discussed in later sections), it affords a consistent comparison with various other literature results.

Fig. 1. Observed bolometric LC (top panel) and Fe IIλ5169 line velocities (bottom panel) of SN 2006ai (filled markers) with models randomly sampled from the posterior distribution of the parameters (solid lines). The median value and 68% confidence range for every physical parameter are also shown. The grey shaded region shows the early data removed from the fitting (first 30 days after explosion).

Having inferred the progenitor and explosion properties for a large sample of SNe II (24 SNe in the gold sample and 53 in total), we now analyse the distributions of the physical parameters. This is the largest set of physical properties of SNe II derived to date. For two of the 53 SNe II (SN 2005af and SN 2007it) only M_Ni was derived. SN 2005af was only observed at late times, from the late recombination phase to the radioactive tail phase. SN 2007it was also observed at late phases and, additionally, it has early observations (< 30 days). These data are insufficient to determine all the physical parameters, and only M_Ni was inferred.

Figures 2 and 3 display histograms of progenitor (M_ZAMS, M_ej, and R) and explosion parameters (E, M_Ni, and ⁵⁶Ni mixing), respectively. In each panel, the yellow bars represent SNe II in the gold sample, while the histogram of the full sample is indicated in red bars. We characterise the distributions by the 16th, 50th, and 84th percentiles. Additionally, Figs. 2 and 3 show the probability density function of the parameters. These functions were computed using a kernel density estimation of the sum of the posterior distributions marginalised over the parameters for the SNe II in the gold (solid lines) and full (dashed lines) samples. We note that some values exceed unity. This is not erroneous since the figure shows the probability density function. The integral of the probability density function along the entire range of values equals unity.

Fig. 2. Histograms of three progenitor parameters: MZAMS (top panel), Mej (middle panel), and R (bottom panel). The gold sample is represented by yellow bars, and the red bars are the histograms for the full sample. In each panel, the number of SNe II is listed, together with the median and the 16th and 84th percentiles. Probability density functions of the physical parameters for the gold and full samples are represented by solid and dashed lines, respectively.

Fig. 3. Same as in Fig. 2 but for E (top panel), MNi (middle panel), and 56Ni mixing (bottom panel).

The M_ZAMS distribution for the SNe II in the gold sample (Fig. 2, top panel) is characterised by a median value of 10.4 M_⊙, with most progenitors (21 of 24) being less massive than 13 M_⊙. When the full sample is considered the median is slightly larger, M_ZAMS = 11.7 M_⊙. For both the gold and full samples, the lowest and highest M_ZAMS values are found to be 9.2 and 20.9 M_⊙ for SN 2009N and SN 2008ag, respectively. SN 2008ag presents one of the longest and brightest optically thick phases in the sample (Paper I), and thus it is not surprising to find a massive progenitor with large pre-SN mass and radius to reproduce its observations. None of the SNe II in the sample is consistent with explosion models for stars more massive than 21 M_⊙. This is in accordance with several studies that analyse SN II progenitors in pre-explosion images, that also find a lack of high-mass RSG progenitors (Smartt 2015; Davies & Beasor 2018). On the other hand, the lower boundary of our M_ZAMS distribution is consistent with the lowest-M_ZAMS progenitor model. In Sect. 6.1, a full analysis of the M_ZAMS distribution is given.

The histograms of M_ej (Fig. 2, middle panel) and R (Fig. 2, bottom panel) can be explained from the M_ZAMS distributions. The M_ej and R are not model parameters, and therefore these values were directly inferred from the M_ZAMS value derived from the fitting. As more low-M_ZAMS progenitors are recovered, it is expected that most of the SNe II have low M_ej and R (under the assumption of standard single-star evolution used in this study). The median M_ej is 8.4 M_⊙ for the gold sample, and 9.2 M_⊙ when all SNe II are analysed. As for the M_ZAMS, the lowest and highest M_ej and R are obtained for SN 2009N and SN 2008ag, respectively (for both gold and full samples). The values of M_ej range from 7.9 to 14.8 M_⊙, and from 450 to 1077 R_⊙ for the pre-explosion radius. The median R is 495 R_⊙ for the gold sample, and 582 R_⊙ for the full sample.

The top panel of Fig. 3 shows the histograms for the explosion energy. The gold sample ranges from 0.15 foe (SN 2008bk) to 1.25 foe (SN 2006ai) with a median value of 0.61 foe. The full sample spans a larger range of E from 0.15 foe (SN 2008bk) to 1.40 foe (SN 2006bl), with a median value of 0.63 foe. The M_Ni distribution is displayed in the middle panel of Fig. 3. We could determine the ⁵⁶Ni mass for 17 SNe II as only these present observations during the radioactive tail. Uncertainties in the explosion epoch, host-galaxy extinction, and distance impact on the determination of M_Ni. While some SNe II have lower uncertainties in the parameters mentioned above, these SNe are not necessarily the same as those in the gold sample above. Thus, to avoid confusion with the previously defined sub-samples, we did not separate SNe II into different samples for the derivation of their ⁵⁶Ni masses. The median of the distribution is 0.036 M_⊙, and it ranges from 0.006 M_⊙ for SN 2008bk, to 0.069 M_⊙ for SN 2007X. In Paper I, we found that SN 2008bk is the lowest-luminosity event at the radioactive tail phase, and therefore it is unsurprising that SN 2008bk has the lowest estimated M_Ni. The bottom panel of Fig. 3 shows the ⁵⁶Ni mixing distributions. ⁵⁶Ni mixing values range from the inner ejecta (0.2) to the outer ejecta (0.8) with a median of 0.53 for the gold sample, similar to that for the full sample. Only a few SN II explosions are consistent with extended ⁵⁶Ni mixing (∼0.8) into the ejecta.

6. Analysis

6.1. The mass function of SN II progenitors

Most studies of the initial mass function (IMF) of SN II progenitors use information from pre-explosion images (e.g. Smartt et al. 2009; Smartt 2015; Davies & Beasor 2018, although see Morozova et al. 2018). A discussion of their results is presented in Sect. 7.1. In the current work, we study the M_ZAMS distribution of SN II progenitors using the values derived for the CSP-I sample. Restricting ourselves to the gold sample we have 24 M_ZAMS estimations. This number increases to 51 if we consider the full sample. To date, this is the most homogeneous and largest sample of progenitors of SNe II analysed from hydrodynamical modelling. Our analysis is based on fitting a power law to the inferred M_ZAMS distribution employing a Monte Carlo (MC) method to randomly sample the masses from a master population. This technique is similar to that presented by Davies & Beasor (2020), where the reader is referred for details.

We first constructed the cumulative distribution (CD) of the M_ZAMS for the progenitors. We performed 10 000 MC simulations to determine the probability distribution of the CD. For each simulation, M_ZAMS was randomly sampled from the probability distribution of each progenitor, which was obtained from our MCMC procedure (Sect. 4.2), and then the progenitors were ordered in increasing M_ZAMS. As noted in Davies & Beasor (2020), when sampling from the posterior distribution, the M_ZAMS of a progenitor can take different values producing changes in the order of the progenitors. For this reason, it is necessary to re-order the progenitor masses at each MC trial. Figures 4 and 5 show the progenitor masses in increasing order (filled squares) and the confidence regions of the CD of the masses (filled contours) for the gold and full samples, respectively.

Fig. 4. Cumulative distribution of MZAMS for the SN II progenitors in the gold sample. The derived masses are shown in yellow squares. The shaded contours show the confidence regions of the CD.

Fig. 5. Same as in Fig. 4 but for the SN II progenitors in the full sample. Yellow squares represent the SNe II in the gold sample.

Then, the CDs derived above (i.e. the shaded contours in Figs. 4 and 5) were fitted with models of CD based on power laws. These models have the following input parameters: the lower and upper M_ZAMS limits of the distribution (M_low and M_high, respectively), and the slope of the power law (Γ). For each set of parameters, we employed 10 000 MC trials to determine the posterior probability distribution for the model CD of the masses. We found the probability distributions of the parameters via MCMC methods. Results are characterised by the median of the marginal distributions, adopting the 16th and 84th percentiles as the lower and upper uncertainties, with one exception (see below).

In the following, we focus on the analysis of the gold sample. The full sample is analysed later. Figure 6 shows the CD of M_ZAMS obtained for the gold sample in comparison with the theoretical CD based on power laws. The yellow triangles and error bars represent the median and the 68% confidence level of the distribution presented in Fig. 4 as shaded contours. The shaded regions in Fig. 6 correspond to the theoretical CD constructed with the median values of the marginal distributions of the parameters. As previously indicated, 10 000 MC simulations were performed to construct this theoretical CD. A corner plot of the joint posterior probability distribution of the parameters is presented in the right panel of Fig. 6. We find M_low = 9.3${}_{-0.1}^{+0.1}$ M_⊙, M_high = 24.7${}_{-3.2}^{+4.2}$ M_⊙, and Γ = −6.35${}_{-0.57}^{+0.52}$. These values are also listed in Table 2. We note that the median value of M_high is offset to larger masses from the peak of the distribution, that is, the distribution is skewed to the right. For this reason, we used the mode to characterise the M_high distribution and the 68% confidence interval for the uncertainties, finding $M_{\mathrm{high}}^{\mathrm{mode}}$ = 21.3${}_{-0.4}^{+3.8}$ M_⊙. The right-skewed distribution is because of the steepness of the CD of the progenitors. The mass function of SN II progenitors we derive is much steeper than that for a Salpeter massive-star IMF with Γ = −2.35 (Salpeter 1955). With a steep power law of Γ = −6.35 we do not expect high-mass progenitors to be within our sample. The vast majority of the stars are found near the value of M_low, and changes in M_high to larger values do not significantly alter the distribution. For a standard Salpeter IMF, 90% of the stars between 9−25 M_⊙ (the range of our pre-SN models) are in the 9−21 M_⊙ range. Then, in our sample of 24 progenitors, we would expect to find two progenitors with masses above 21 M_⊙. This is not the case for our defined Γ. Not only we do not have progenitors with M_ZAMS above 21 M_⊙, we also find a large number of low-mass progenitors. In the gold sample, 87% of the progenitors have M_ZAMS < 13 M_⊙, while the expected value for a Salpeter IMF is ∼50%.

Fig. 6. Fits to the CD of MZAMS of SN II progenitors in the gold sample. Left panel: CD of MZAMS derived in our work for the gold sample (yellow triangles and error bars), corresponding to the median and the 68% confidence limit of the distribution presented in Fig. 4 as shaded contours, in comparison with the model CD constructed with the median values of the marginal distributions of the parameters (shaded contours). Right panel: corner plot of the joint posterior probability distribution of the parameters. Dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions.

Following the above results when leaving the power law as a free parameter, we now compare the CD of M_ZAMS to models with the power-law slope set to −2.35 (i.e. a standard Salpeter IMF). In this present case, we find M_low = 8.5${}_{-0.1}^{+0.1}$ M_⊙ and M_high = 19.1${}_{-0.2}^{+0.2}$ M_⊙ (Fig. 7, right panel). The CD model drawn from the posterior distribution of the parameters shows discrepancies with the CD of progenitor masses derived from our modelling (Fig. 7, left panel). The differences are due to the large number of low-M_ZAMS progenitors found.

Fig. 7. Same as in Fig. 6 but when the power-law slope is fixed to Γ = −2.35, that is, the slope of a Salpeter massive-star IMF.

The same analysis was performed for the full sample of SN II progenitors. Here, we find the following parameter values: M_low = 9.3${}_{-0.1}^{+0.1}$ M_⊙, M_high = 21.5${}_{-0.8}^{+1.2}$ M_⊙, and Γ = −4.07${}_{-0.29}^{+0.29}$. In Fig. 8, we show the CD of M_ZAMS obtained for the full sample (red triangles and error bars corresponding to the median and 68% confidence limit of the distribution presented in Fig. 5 as shaded contours) in comparison with the model CD constructed with the median values of the marginal distributions (shaded contours). In addition, Fig. 8 shows a corner plot of the posterior probability distribution of the parameters. The values of M_low and M_high are similar for both the gold and full samples. The largest differences are found for the power-law slope. This is shallower than that obtained for the gold sample, although it is still steeper than a Salpeter IMF. In Sect. 7.2.1 we analyse whether the different sample sizes influence the power-law slope obtained.

Fig. 8. Fits to the CD of MZAMS of SN II progenitors in the full sample. Top left: CD of MZAMS derived in our work for the full sample (red triangles and the error bars), corresponding to the median and the 68% confidence limit of the distribution presented in Fig. 5 as shaded contours, in comparison with the model CD constructed with the median values of the marginal distributions of the parameters (shaded contours). Bottom left: corner plot of the joint posterior probability distribution of the parameters. Dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions. The right panels show the same as the left panels but when the power-law slope is fixed to Γ = −2.35.

Additionally, we computed the lower and upper mass limits if we constrain Γ = −2.35 (top-right and bottom-right panels of Fig. 8). We find M_low = 9.0${}_{-0.1}^{+0.1}$ M_⊙ and M_high = 19.0${}_{-0.2}^{+0.2}$ M_⊙. Here, the model drawn from the median values of the parameter distributions can mostly reproduce the behaviour of the M_ZAMS distribution. At first glance, M_high is smaller than the upper mass boundary of the derived cumulative mass distribution (red triangle in the top-right panel of Fig. 8). However, the upper mass boundary has a value of M_ZAMS = 21.0${}_{-2.2}^{+1.0}$ M_⊙ (99.7% confidence), that is, they are coincident at the 99.7% significance level. In both the gold and full samples, we find steeper power laws than that of a Salpeter massive-star IMF. In Sect. 7 we discuss possible reasons of this discrepancy, which we term ‘the IMF incompatibility’.

6.2. Ejecta masses and explosion energies

Here we discuss the range of ejecta masses and explosion energies obtained in our work. Figure 9 compares our results to those from other studies involving large samples of SNe II. For this purpose, we used the results published in Utrobin & Chugai (2019), which summarises previous estimates from a series of papers by those authors and the results of Pumo et al. (2017) on under-luminous SNe II. Pumo et al. (2017) also complies the results of several ‘normal’ SNe II obtained with the same code. We also compare with the ejecta masses and explosion energies from Ricks & Dwarkadas (2019). All the above-mentioned studies present detailed modelling of LCs and expansion velocities for several SNe II. Morozova et al. (2018) presented LC modelling for 20 SNe II. Given that our distributions were derived from the modelling of the LC and expansion velocities simultaneously, a direct comparison is difficult. However, we include the results from Morozova et al. (2018) in the comparison. Unfortunately, we cannot compare with the findings from Eldridge et al. (2019) as only the initial masses were published.

Fig. 9. Comparison of ejecta masses (top panel) and explosion energies (bottom panel) derived in the current work to those from previous studies containing large samples of SNe II.

We observe that our M_ej and E are generally consistent with those from previous studies as our results fall within the range of previous estimates. However, we note a spike at ∼8−9 M_⊙ that is not found in the other studies⁸. This spike is the reason why we found so many low-M_ZAMS progenitors, causing the IMF incompatibility (see Sect. 7.2).

The highest estimated M_ej and E for the CSP-I SN II sample are 14.8 M_⊙ and 1.25 foe (1.40 foe for the full sample), respectively. A number of higher M_ej and E are found in the other works, especially in those from Utrobin & Chugai (2019). This could be because these authors use non-evolutionary models as pre-SN configurations that are capable of producing a large variety of pre-explosion structures, each with different mass, radius, chemical composition, and density profile. In the present work, we used standard single-star evolution calculations as pre-SN models. Therefore, the initial density profiles were not freely chosen, and the range of pre-SN parameters is confined. In our set of progenitor models, the highest M_ej is 15.7 M_⊙, coming from a 24 M_⊙ progenitor in the ZAMS (see M20, their Fig. 2). Therefore, it is impossible for us to find larger M_ej given that we only consider M_ZAMS up to 25 M_⊙. Ricks & Dwarkadas (2019) obtained a M_ej of ∼20 M_⊙ for SN 2015ba using evolutionary models to initialise the explosion. However, for this particular progenitor model, the authors used a lower scaling factor for the efficiency of mass loss via winds during the stellar evolution producing a more massive star at explosion time. We again stress that our pre-SN models were calculated using the standard assumptions and parameters during the evolution. Different stellar evolution assumptions can lead to different results.

Pumo et al. (2017) found, on average, lower M_ej and E than other studies probably because most of their sample corresponds to low-luminosity SNe II. While the range of parameters found is in good agreement with that from Morozova et al. (2018), we find some differences in the distributions. Morozova et al. (2018) found more SNe II with lower energies, while at the same time, a spike at M_ej ∼ 14 M_⊙.

Finally, we compare the explosion energies in this work with those expected from neutrino-driven explosion models. 1D core-collapse models found explosion energies ranging from 0.1 to 2.0 foe (Ugliano et al. 2012; Ertl et al. 2016; Sukhbold et al. 2016). The SNe II modelled in this work do not cover the entire range of explosion energies predicted by the explosion models, particularly because our maximum value is somewhat lower (1.25 foe for the gold sample and 1.40 foe for the full sample). More energetic explosions produce more luminous SNe II, which are not found in our sample. For example, the highest explosion energies found in Sukhbold et al. (2016) produce plateau luminosities higher than 10^42.6 erg s⁻¹ (see their Fig. 31), while in the CSP-I SN II sample, only one object – SN 2009aj – has plateau luminosities above this value. In addition, as previously discussed in Sect. 5, SN 2009aj may be powered by interaction between the ejecta and a massive CSM. More recently, Burrows et al. (2020) conducted 3D core-collapse simulations in a wide range of progenitor masses. The lowest-M_ZAMS model (9 M_⊙) in their grid achieves an explosion energy of ∼0.1 foe, which is compatible with the lowest explosion energy found in our study. Unfortunately, the other core-collapse models have not reached the asymptotic explosion energy by the end of the simulations. However, these authors do find that progenitors with higher M_ZAMS attain higher explosion energies. We find the same trend in Paper III.

6.3. Progenitor radii

Levesque et al. (2005, 2006) studied Galactic and Magellanic Cloud RSGs and obtained their radii using the Stefan-Boltzman law from effective temperatures and bolometric luminosities. These studies found that most RSGs have radii between 100 and 1500 R_⊙. In Sect. 5, most progenitor radii inferred from our LC and velocity modelling are found within the range of 450−600 R_⊙. However, we emphasise that progenitor radius is not an independent parameter in our modelling. Our fitting routine derives M_ZAMS, which relates to pre-SN structures with different characteristics (e.g. ejecta mass, radius, among other physical parameters). Additionally, while early-time optical LCs are most sensitive to progenitor radius, in the current work fits were performed only to observations later than 30 days from explosion (see Sect. 4.3).

6.4. ⁵⁶Ni masses

The ⁵⁶Ni masses for the CSP-I SN II sample range from 0.006 M_⊙ through 0.069 M_⊙ with a median of 0.036 M_⊙. Our values were determined for 17 SNe II which do not cover the entire ⁵⁶Ni mass range found in the literature. Some SNe II have been found with more and less ⁵⁶Ni (e.g. Hamuy 2003; Pastorello et al. 2004; Pumo & Zampieri 2011). M_Ni was estimated for all CSP-I SNe II with bolometric data during the radioactive tail phase⁹, suggesting that the M_Ni range of our models (see Table 1) was sufficient and the inclusion of models with M_Ni larger than 0.08 M_⊙ into our grid of simulations was not necessary.

Pejcha & Prieto (2015) estimated M_Ni for 21 SNe II and found a minimum of 0.0045 ± 0.0008 M_⊙ for SN 2001dc¹⁰ and a median of 0.030 M_⊙. Both the median and minimum values are consistent with our determination. However, a large discrepancy is found for the maximum value. Pejcha & Prieto (2015) inferred a ⁵⁶Ni mass yield of 0.28 M_⊙ for SN 1992H, which is significantly higher than the 0.069 M_⊙ estimated in our study for SN 2007X. Pejcha & Prieto (2015) inferred a larger distance to SN 1992H than previous estimates (Schmidt et al. 1994; Clocchiatti et al. 1996) making this SN more luminous compared to normal SNe II, and therefore, with a larger M_Ni.

Müller et al. (2017) determined M_Ni for 19 SNe II and combined their sample with that from Pejcha & Prieto (2015) making a larger sample of 38 objects. Their M_Ni distribution is described by a median value of 0.031 M_⊙ – consistent with ours – with values between 0.005 and 0.280 M_⊙. The latter value corresponds to SN 1992H as described above. More recently, from a compilation of 115 SN II ⁵⁶Ni mass estimates from the literature, Anderson (2019) estimated a median value of 0.032 M_⊙, in agreement with our results. The minimum and maximum values in Anderson (2019) are 0.001 and 0.360 M_⊙, respectively. This range is larger than that determined in our study, particularly because the high ⁵⁶Ni mass estimate of 0.360 M_⊙ (for SN 1992am; Nadyozhin 2003) is significantly larger than the maximum value estimated for the CSP-I sample. The larger sample in Anderson (2019) may indicate that values above the adopted maximum in our models of 0.08 M_⊙ are exceptional. Around 10% of the SNe II in the literature have estimated ⁵⁶Ni masses above 0.08 M_⊙ (see Fig. 1 from Anderson 2019). However, some of these SNe II are classified as 1987A-like events, which have relatively high ⁵⁶Ni masses compared to the ‘normal’ SNe II studied in the current work. Additionally, a few SNe II have more than one ⁵⁶Ni mass estimate in the literature, and some of these are consistent with our adopted maximum ⁵⁶Ni mass. Differences are attributed to distinct distance and reddening values used in the literature for the same SN.

Sukhbold et al. (2016) report ⁵⁶Ni masses based on neutrino-powered explosions for numerous pre-SN models within a large range of initial masses. The minimum M_Ni values in Sukhbold et al. (2016) are 0.003 and 0.006 M_⊙ for the 9.25 and 9.0 M_⊙ models, respectively. The different values correspond to the different codes used for their calculation. These values are similar to the lowest M_Ni derived for our SN II sample. For stars between 10 and 25 M_⊙ (the maximum initial mass in our pre-SN models), most M_Ni values are found within the range of 0.01−0.08 M_⊙, which is consistent with our findings. Only a few explosion models produce larger M_Ni in Sukhbold et al. (2016), with a maximum value of ∼0.10 M_⊙.

6.5. ⁵⁶Ni mixing

During a SN explosion, a shock wave emerges heating the stellar material. If the shock temperature is sufficiently high, explosive nucleosynthesis takes place. Explosive nucleosynthesis efficiently produces heavy elements from Si to Zn (Umeda & Nomoto 2002). ⁵⁶Ni dominates the production of nuclear species in the inner regions of the star, but Co, Zn, and additional Ni isotopes are also mostly produced here (Woosley & Weaver 1995; Thielemann et al. 1996; Umeda & Nomoto 2002). Once the shock reaches the composition interfaces between the carbon-oxygen core and the helium core, and the helium core and the hydrogen-rich envelope, RT instabilities appear causing a large-scale spatial mixing. This is the mechanism by which the heavy elements synthesised during explosive burning can reach the outer regions of the envelope while hydrogen can be mixed inwards.

The first evidence of chemical mixing during the explosion was based on the observations of SN 1987A. A substantial amount of ⁵⁶Ni mixed into the hydrogen-rich envelope is required to explain the rise time to peak of SN 1987A (Shigeyama et al. 1988; Shigeyama & Nomoto 1990), and the LC and spectral evolution of X-rays and γ-rays (Kumagai et al. 1989). Furthermore, the presence of significant Fe emission at high velocities during the nebular phase demonstrates that a significant amount of ⁵⁶Ni was mixed into the hydrogen-rich envelope during the explosion (Haas et al. 1990). While the progenitor of SN 1987A is a blue supergiant (BSG) star, SN II progenitors are known to be RSGs with significantly distinct pre-SN structures. Recent 3D hydrodynamical simulations of the evolution of the SN shock for both BSG and RSG explosions have shown that RSG models achieve higher maximum velocities for the ⁵⁶Ni into the ejecta (Wongwathanarat et al. 2015) implying more extended mixing of ⁵⁶Ni.

In Sect. 5 we showed that the ⁵⁶Ni distribution in the ejecta covers the range of initial values, from 0.2 (inner ejecta) to 0.8 (outer ejecta). However, in both the gold and full samples, the vast majority of the events are found within low (0.2−0.4) and moderate (0.4−0.6) values and only a low number of SNe II are found with a large extent of ⁵⁶Ni mixing (Fig. 3, bottom panel). An extreme mixing of ⁵⁶Ni is hard to obtain in RSGs because of their large envelopes. Our results are consistent with this concept given that an extended ⁵⁶Ni mixing was only found for one SN II in the gold sample (SN 2008ag, ⁵⁶Ni mixing = 0.8). We note that the ⁵⁶Ni mixing is characterised as a fraction of the pre-SN mass. This means that the same value of ⁵⁶Ni mixing for different progenitors may represent different extents of ⁵⁶Ni in mass and radial coordinates.

We stress that in this study the ⁵⁶Ni mixing is treated as a free parameter. Recently, Wongwathanarat et al. (2015) computed 3D hydrodynamical simulations and found that the extent of the mixing of the metal-rich ejecta depends on the early-time asymmetries generated by the neutrino-driven mechanism, explosion energy, and progenitor structure. In the progenitor structure, the most important features are the width of the carbon-oxygen core, the density structure of the helium core, and the density gradient at the composition interface between the helium core and hydrogen-rich envelope. However, only a few hydrodynamical simulations for RSG models have been carried out in the literature. Further studies are crucial for better understanding of the ⁵⁶Ni mixing process in SNe II.

6.6. Explosion epochs

We compared the explosion epochs derived in our analysis with those from Gutiérrez et al. (2017b) for each SN. Our estimations are always inside the range of values derived by Gutiérrez et al. (2017b), which is by design since we did not allow our fitting procedure to sample values outside that range (see Sect. 4.2). The mean difference between both estimates is 0.1 days with a standard deviation of 4.8 days. There are a few cases where our explosion epoch estimates are very close to the observational limits, possibly suggesting explosion dates beyond the uncertainties of Gutiérrez et al. (2017b). In M20, we tested this by relaxing the prior for the explosion epoch and obtained a similar marginal distribution of the physical parameters, indicating no changes in our results.

That the mean offset is so close to zero gives significant support to both methodologies of estimating the explosion epoch. While we did not allow our modelling estimates to be outside the range of the observational errors, the modelling explosion epochs could have been biased towards – on average – higher or lower values than the observational values. In the case of observational explosion epochs, the constraints from the non-detection that the SN has not exploded at that date are only as strong as the survey depth providing the non detection. One could therefore imagine that ‘true’ explosion epochs could be biased towards values close to the non detection. However, the above analysis shows this not to be the case, and therefore gives support to the observation epochs as estimated by Gutiérrez et al. (2017b). At the same time, one could imagine that a bias in the modelling technique could mean that the fitting could attempt to use the prior on the explosion epoch to systematically offset the plateau duration – and therefore the ejecta mass – to lower or higher values by constraining the explosion epoch to be later or earlier, respectively (with respect to the observational value). Again, this is not seen in our data, giving support to the robustness of our model fitting methods.

6.7. Non-standard SNe II

In this study, we utilise a large grid of explosion models assuming standard values for several stellar evolution parameters (metallicity, wind efficiency, mixing length, and overshooting) and explosion parameters (M20). Although most of the well-sampled SNe II in the CSP-I sample can be adequately reproduced by the models in this grid, some are not.

SN 2006Y and SN 2008bu reside in the lower end of the optically thick phase duration distribution with 64 ± 4 days and 52 ± 7 days, respectively (Paper I). Unfortunately, no model in our grid has such a short optd. In the case of SN 2006Y, the closest fitting models can reproduce the velocities and the plateau luminosity reasonably well. However, the optd is overestimated by ∼10 days. At first sight, this difference is small, but the models also largely underestimate the radioactive tail luminosity by ∼0.6 dex. Additional ⁵⁶Ni can easily fit the tail, but it is known that more ⁵⁶Ni also produces a longer plateau phase (e.g. Kasen & Woosley 2009; Bersten 2013), which increases the discrepancy with the observations. On the other hand, for SN 2008bu, the situation is even worse as the optd lasts ∼10 days less than that of SN 2006Y.

Theoretical studies suggest that progenitors with smaller hydrogen-rich envelope masses at time of collapse produce shorter optd (e.g. Litvinova & Nadezhin 1983; Kasen & Woosley 2009; Bersten 2013). Following the standard mass-loss rate in the stellar modelling of single stars (as we assume in this study) it is not possible to find such stripped progenitors with M_ZAMS smaller than 25 M_⊙, which is the highest-mass progenitor in our grid. More massive stars can finish with smaller envelopes but they are more difficult to explode (Sukhbold & Adams 2020). Therefore, lower ejecta masses can be achieved by increasing the mass-loss rate via stellar winds, interacting binaries, or rotation, among other possibilities (as discussed in Sect. 7). We calculated new pre-SN models with enhanced mass loss providing much better fits to SN 2006Y and SN 2008bu. These results are presented in Paper III. Our results agree with those from Hiramatsu et al. (2021). These authors presented models for SN 2006Y, as well as other two short-plateau SNe II, from high-mass progenitors that experience more mass loss than standard stellar evolution.

SN 2004er also shows some discrepancies between models and observations, especially with the velocities. SN 2004er shows the longest optd with 146 ± 2 days (Paper I). Its bolometric LC is mostly well-reproduced with the exception that the fitting models show shorter optd by ∼10 days. The largest discrepancy is found in the velocity evolution, where the observed Fe II velocities are greater than the models by ∼2000 km s⁻¹. In fact, SN 2004er displays higher Fe II velocities than most SNe II. For example, the Fe II velocity at 55 days post-explosion is 5033 ± 751 km s⁻¹, while the mean Fe II velocity of SNe II at 53 days is 3537 ± 851 km s⁻¹ (Gutiérrez et al. 2017b). Such large velocities can be achieved through higher explosion energies, although this would produce more luminous and shorter plateau phases. However, larger ejecta masses have the opposite effect, that is, dimmer SNe II with longer plateaus. Therefore, SN 2004er was probably a high-energy explosion of a star with a massive pre-SN envelope – a model that is beyond the parameter space sampled in the current study (and therefore beyond stellar evolution with standard assumptions and/or values of the input parameters). Something similar happens for SN 2007sq. Its bolometric LC is well fit by our standard models, but our model velocities fall below observations. The SNe II above analysed (SNe 2004er, 2006Y, 2007sq, and 2008bu) may be examples of high-mass RSG progenitors that could explain the lack of more massive progenitors in the sample, although a detailed analysis is necessary.

Finally, we briefly describe the fitting models of SN 2008K. These models show offsets both in the bolometric LC and velocities with respect to the observations. The bolometric LC is well fit in general, but when looking at the details, it is seen that the models do not resemble the linear behaviour of SN 2008K, which has a plateau decline rate of 1.70 ± 0.31 mag per 100 days (Paper I). SN 2008K also has higher Fe II velocities than most SNe II. In principle, more energy is needed to reproduce the velocities. This implies more mass to conserve the plateau duration (more energy reduces the plateau length). However, we compared SN 2008K with higher-energy and more-massive explosion models and found even larger discrepancies given that the new parameters produce LCs similar to typical SNe IIP. The detailed pre-SN and explosion properties of SN 2008K are still uncertain. The examples above reflect that some of the standard assumptions in stellar evolution cannot cover all the physical properties of SN II progenitors. Non-standard stellar evolution is required to explain the properties of some individual SNe II.

7. Discussion

7.1. The RSG problem

During recent years, the number of works studying the M_ZAMS distribution of SN II progenitors has increased. Smartt et al. (2009) were the first to analyse the M_ZAMS distribution from observed progenitors in pre-explosion images, founding a maximum M_ZAMS of 16.5 ± 1.5 M_⊙. This is in contrast with model predictions and that more massive RSGs are observed in the Local Group (Levesque et al. 2006; Neugent et al. 2020). The lack of higher-mass progenitors is called the ‘RSG problem’. Since then, special attention has been given to the upper mass boundary of the M_ZAMS distribution. This limit is crucial to understanding the evolutionary pathways of massive stars. A larger set of SN II progenitors consisting in 13 detections and the same number of upper limits was studied by Smartt (2015), but no progenitors were found with M_ZAMS above 18 M_⊙.

Inferring the progenitor M_ZAMS from a direct detection requires a previous estimation of the progenitor luminosity. This is used to compare with the luminosity predicted by stellar evolution theory. Thus, the physical models provide an estimate of the mass from a luminosity measurement. In this context, Straniero et al. (2019) studied the role of convection, rotation, and binarity in SN II progenitor evolution finding that neither of these uncertain processes in stellar modelling can mitigate the RSG problem.

The progenitor luminosity can be achieved by fits to the spectral energy distribution if the progenitor is detected in several bands, or by using bolometric corrections to convert single-band flux into luminosity. The uncertainties in the latter case may be large if only single-band or two-band detections are available. Davies & Beasor (2018) investigated this source of systematic error and derived new bolometric corrections and extinction values. Those authors found an increased upper mass limit of M_high = 19${}_{-1.3}^{+2.5}{M}_{\odot }$ (68% confidence), with a 95% upper confidence limit of < 27 M_⊙. Additionally, Davies & Beasor (2018) analysed the effects of a small sample size on M_high, concluding that this causes a systematic error of ∼2 M_⊙ that shifts the upper mass limit to larger values. Recently, Davies & Beasor (2020) inferred the properties of the progenitor distribution but working directly with the observed luminosities, thus eliminating the uncertainties introduced in stellar modelling when converting the observed luminosities into M_ZAMS. Finally, they compared their maximum luminosity with stellar models and found a M_high of 18−20 M_⊙. However, they remark that the sample size should be at least doubled to enable a reduction in the large uncertainties on M_high.

The direct detection of progenitors in pre-explosion images is the most powerful tool to determine the nature of progenitor stars as it can be directly linked to a progenitor system. However, the analysis of direct detections can only be applied – when pre-SN images exist – to the nearest SNe (d ≲ 30 Mpc) due to the lack of resolution for more distant objects. The study of the stellar populations in the immediate SN environments is an alternative indirect method that is also useful for nearby objects where individual stars or clusters are resolved. This information can yield constraints on the ages, and therefore, on the progenitor masses. These measurements are mostly consistent with initial masses of SN II progenitors lower than 20 M_⊙, although they present larger uncertainties than the direct identification of the progenitor star (e.g. Van Dyk et al. 1999; Williams et al. 2014; Maund 2017).

Nebular spectra can also be used to constrain the progenitor mass. After the optically thick phase, when all the hydrogen is recombined, the core of the progenitor star becomes visible and nucleosynthesis yields can be analysed (Jerkstrand et al. 2012; Dessart & Hillier 2020). Using this method, most of the mass estimates agree with the lack of high-mass progenitors (see e.g. Jerkstrand et al. 2014; Valenti et al. 2016). However, Anderson et al. (2018) infer an initial mass of 17−25 M_⊙ for SN 2015bs through comparison to nebular-phase spectral models suggesting an explosion of a higher-mass progenitor than previously observed for a SN II.

Different alternatives have been proposed to solve the RSG problem. One possibility is that high-mass progenitors do not explode and instead collapse into black holes. Massive stars in the range of the missing RSG progenitors could have pre-SN core structures that are generally more difficult to explode (see e.g. O’Connor & Ott 2011; Sukhbold & Adams 2020). The search for failed SNe has found possible candidates corresponding to high-mass stars (Adams et al. 2017; Allan et al. 2020; Neustadt et al. 2021), lending support to theoretical predictions. In addition, it has been suggested that more massive RSGs produce weak explosions leading to faint SNe II (Fryer 1999; Heger et al. 2003). Zampieri et al. (2003) investigated the possibility that two low-luminosity SNe II (SN 1997D and SN 1999br) could be produced in low-energy explosions of massive progenitors. An alternative solution posits that stellar evolution models could be modified in order to increase the mass loss and reduce the maximum theoretical mass of stars exploding as a RSG. With this assumption, massive stars in the range of the missing RSG progenitors will explode as stripped-envelope SNe instead of SN II events (see e.g. Yoon & Cantiello 2010; Ekström et al. 2012).

Unfortunately, most of the previously mentioned methods to constrain progenitor parameters of SNe are restricted to nearby events because they require either images of high enough resolution and sensitivity or spectra at late-time epochs when the SN is significantly dim. Hydrodynamical modelling of LCs (and sometimes spectra) is one of the most commonly applied indirect methods to infer physical properties of SNe II as it can be applied to large distances, and therefore statistically significant samples. Morozova et al. (2018) inferred progenitor and explosion properties of a sample of 20 SNe II from LC modelling and used the derived M_ZAMS to explore the mass distribution. They found an upper mass cut-off of M_high = 22.9 M_⊙ (∼30 M_⊙ at 95% confidence limit). However, Morozova et al. (2018) did not consider expansion velocities in their modelling, the exclusion of which has been shown to bias progenitor masses to larger values (see M20).

In Sect. 6.1 we determined the IMF parameters (M_low, M_high, and Γ) for the gold and full samples. For each sample, we also fitted the IMF leaving Γ as a free parameter or fixing it to −2.35. Thus, in total, we have four different estimations of the IMF parameters. All these estimations find M_high ≲ 21.5 M_⊙. In principle, this is consistent with the RSG problem. However, different levels of significance are found for each estimation, which are examined below. We first analyse the cases when Γ is left as a free parameter. We find M_high = 21.3${}_{-0.4}^{+3.8}$ M_⊙ for the gold sample, that is, M_high = 25 M_⊙ is within the 68% confidence limit. The study of the full sample indicates M_high = 21.5${}_{-0.8}^{+1.2}$ M_⊙ (68% confidence) with a 95% upper confidence limit of < 27.1 M_⊙. While this increases the significance of our results with respect to the RSG problem, we note that the results obtained from the full sample are not as robust as those from the gold sample. When Γ is fixed to the value of a Salpeter massive-star IMF, we find M_high < 21 M_⊙ for both the gold and full samples. In these two cases, the RSG problem is recovered with a significance above the 99% significance level.

7.2. The IMF incompatibility

The SN II progenitor IMF derived in Sect. 6.1 is inconsistent with a Salpeter IMF due to the large difference found in the slope of the distributions – which we named the IMF incompatibility. However, the massive-star IMF is very well established. Several studies of OB associations and clusters suggest that the majority of the stars above a few solar masses are drawn from a Salpeter power-law (Bastian et al. 2010). In addition, progenitor mass constraints from the stellar populations surrounding a SN event and SN II progenitor mass estimations from pre-SN imaging are also drawn from a Salpeter IMF (Smartt 2015; Koplitz et al. 2021). The apparent incompatibility of our results with a standard IMF is driven by the large number of low-M_ZAMS progenitors found with respect to more massive progenitors. In Sect. 3 we pointed out that CSP-I built a magnitude-limited sample of SNe II. The luminosity function of a magnitude-limited distribution could be biased against lower-luminosity events that generally come from low-M_ZAMS progenitors and low-energy explosions (see Paper III). This is the opposite of the volume-limited sample of SNe II with progenitor detection in pre-explosion images (Smartt et al. 2009). Given that most SNe II in our sample are consistent with such low-mass events, there does not seem to be significant evidence that this bias is affecting our results or conclusions. We therefore believe that this incompatibility reflects a lack of understanding of some physical ingredients in our study or the completeness of our sample. Here we discuss possible reasons, but a deep and systematic analysis should be carried out in order to identify the ultimate reason for the IMF incompatibility.

7.2.1. Sample size

First, we explore if the IMF incompatibility is consequence of the sample size. We performed a test similar to that presented in Sect. 6.1 for different sample sizes (N). We randomly sampled the initial masses of N progenitors assuming a Salpeter massive-star IMF with M_low = 9 M_⊙ and M_high = 25 M_⊙ (the M_ZAMS range of our pre-SN models), and ordered the progenitors in increasing mass. We performed 10 000 MC simulations to determine the CD of the masses for each value of N, and then we fit power-law models to these CDs. Given that we are interested in testing the sample size effect in the estimation of Γ, we constrained M_low and M_high to 9 and 25 M_⊙, respectively, which means that Γ was the only parameter to determine. The posterior distribution of the Γ parameter was found via an MCMC procedure. We employed 10 000 MC trials to determine the posterior probability distribution for the progenitors at each Γ value. We performed the same experiment for a range of sample sizes.

The results of this test are presented in Fig. 10. The larger the sample size, the best-fitting Γ tends to the value of the Salpeter IMF slope (−2.35, black solid line). For 24 objects (i.e. the number of SNe II in the gold sample), the best-fitting Γ is −2.50${}_{-0.35}^{+0.25}$, which means that Γ = −2.35 is within the 68% confidence interval. Something similar happens for N = 51 (i.e. the number of SNe II in the full sample with M_ZAMS estimates), where we find −2.45${}_{-0.15}^{+0.15}$. This indicates that our sample size do not produce the bias we found in the IMF. As the IMF incompatibility cannot be explained by the sample size, in the following sections we discuss some possible physical reasons for this discrepancy.

Fig. 10. Best-fitting power-law slope of the CD of the progenitor masses for a range of sample sizes. The error bars correspond to the 68% confidence limit. The larger the sample, the higher the probability that the best-fitting power-law slope is that of a Salpeter IMF. The solid black line indicates the Salpeter IMF slope for massive stars.

From Fig. 10, we also note the best-fitting Γ differs significantly from that of a Salpeter IMF for samples smaller than 20 objects. This restricts future studies to work with larger samples for a reliable estimate of the power-law slope.

7.2.2. Pre-SN models

The IMF incompatibility is produced by the large relative number of low- to high-M_ZAMS progenitors. In Sect. 5 we found a large number of explosions with progenitor structures compatible with M_ej in the range of 8−10 M_⊙ and R between 450 and 650 R_⊙ (Fig. 2, middle and bottom panels). In the context of standard single-star evolution adopted in the current study, these pre-SN structures are consistent with low-M_ZAMS progenitors. However, if the standard assumptions in stellar evolution change (which is possible given the uncertainties in massive-star modelling), different pre-SN configurations could be obtained for the same M_ZAMS value. Below we discuss how different factors such as the envelope stripping processes, metallicity, mixing length, and overshooting can affect the pre-SN structure, and therefore our results.

Mass-loss rate. Similar M_ej (for higher M_ZAMS) can be achieved by increasing the mass-loss rate via stellar winds, eruptive mass-loss events, rotation, binary interaction, or a combination of these phenomena. If more mass were lost by progenitors, then the estimated M_ej could be reduced for more massive M_ZAMS, and thus this would reduce the inconsistency of our results with a standard IMF. While higher mass loss may change the M_ZAMS distribution derived from hydrodynamical modelling, additional mass stripping does not change the progenitor M_ZAMS estimations from pre-SN imaging. This is because the progenitor final luminosity is closely related to its helium core mass and it is almost independent of the final envelope mass (direct progenitor observations do not provide ejecta mass estimates). The sensitivity of the initial mass–final luminosity relationship on the adopted mass-loss rate is small (Straniero et al. 2019; Farrell et al. 2020).

Recent studies infer that most massive stars have rotational velocities in the range of 0−300 km s⁻¹ (Ramírez-Agudelo et al. 2013; Dufton et al. 2013). Since rotation increases mass loss, the final mass of rotating stars is smaller than that of non-rotating stars (Hirschi et al. 2004). Following this, if we consider pre-SN models of rotating massive stars we could find lower-M_ej for higher-M_ZAMS progenitors, that is, in the direction to explain the IMF incompatibility.

Beasor et al. (2020) and Humphreys et al. (2020) derived new mass-loss prescriptions based on direct measurements and did not find larger mass-loss rates in the RSG phase than previous estimations. In addition, Beasor et al. (2020) conclude that quiescent mass-loss during the RSG phase (when most of the mass loss occurs) has little effect in stripping the envelope. In principle, this excludes the RSG winds as a mechanism to strip a larger part of the hydrogen-rich envelope (see also Beasor et al. 2021). However, Ekström et al. (2012) proposed new mass-loss prescriptions for RSGs considering that they may undergo episodic mass-loss events (i.e. more mass is lost). Neugent et al. (2020) compared the observed luminosity function of RSGs in M31 to that predicted by the stellar models from Ekström et al. (2012) and concluded that the new prescription for RSG mass-loss rate reproduces the observations. Given the inconsistencies in previous estimations, it is not clear whether mass-loss processes on isolated RSGs are responsible for stripping a significant part of the envelope.

Very early spectra have shown the presence of CSM around SN II progenitors (e.g. Khazov et al. 2016; Yaron et al. 2017; Bruch et al. 2021) indicating that enhanced mass loss might occur prior to core collapse. This is also inferred from early-time LC modelling (e.g. Morozova et al. 2018; Förster et al. 2018). Such high mass-loss rates are not consistent with direct measurements of RSG winds possibly implying another mechanism to induce mass loss. Fuller (2017) suggests that pre-SN outbursts via waves driven by convection during the late stages of evolution may reproduce the CSM properties. However, despite early-time observations of many SNe II indicate very high mass-loss rates shortly before the explosion, the amount of mass involved is small, which means that this process cannot explain the low pre-SN masses.

An efficient way to remove a large fraction of (if not all) the outer envelope of massive stars is through mass transfer to a companion star in a binary system (e.g. Podsiadlowski et al. 1992; Benvenuto et al. 2013). Therefore, this mechanism can produce SN II events with lower ejecta masses than predicted by single-star evolution. According to recent studies, most massive stars belong to binary systems (or higher-order multiples; Moe & Di Stefano 2017). Moreover, it has been shown that a large fraction of massive stars will exchange mass with a companion prior to the RSG phase (Sana et al. 2012, 2013). Recently, from analytical estimations and population synthesis simulations, Zapartas et al. (2019) showed that a large fraction (from one-third to half) of SN II progenitors will interact with a companion star during its evolution. Additionally, Eldridge et al. (2018) showed that the diversity of CCSN LCs can be interpreted as the result of binary interaction. Previously, it has generally been assumed that SN II arise from progenitors that explode with most of their hydrogen envelopes intact. SN II modelling assuming such progenitors has generally been consistent with observations and thus it has not been necessary to invoke binarity. However, this does not mean that the effects of binary evolution are unimportant in explaining SN II observations and diversity. Our results and the above mentioned studies support the hypothesis that many SN II progenitors may experience larger mass loss than expected assuming standard single-star evolution. In addition, a different treatment of mass loss is required in the cases of SN 2006Y and SN 2008bu (see Sect. 6.7 and Paper III). Thus, it may well be that binarity plays a larger role in shaping SN II behaviour than previously assumed. We now discuss whether other stellar parameters can explain the discrepancy.

Metallicity. The influence of metallicity on the winds of hot-massive stars is well known (Vink et al. 2001). However, for single-massive stars, most of the mass lost happens in the RSG phase, where the rate of mass loss is assumed to be nearly independent of metallicity (van Loon 2006; Goldman et al. 2017) at least within the metallicity range of the CSP-I SN II sample (see below). The influence of metallicity on the winds can change the pre-SN structure at the time of collapse. Lower metallicities produce more compact and more massive progenitors (Dessart et al. 2013). However, we do not expect large variations in the metallicity of the CSP-I sample. Anderson et al. (2016) inferred environment metallicities for the CSP-I SN II sample and found a narrow range in metallicities, which is not expected to significantly alter the stellar evolution. In addition, Anderson et al. (2016) found that there is no relationship between progenitor metallicity and SN II LC and spectral properties (except for the metal line equivalent widths). Gutiérrez et al. (2018) also concluded that the hydrogen-rich envelope mass and the explosion energy are not correlated with the metallicity. Therefore, we do not expect that metallicity significantly modifies our results obtained from solar-metallicity stellar models.

Mixing-length parameter. Pre-SN models with different mixing-length parameters differ primarily in the final radius (larger mixing-length parameters produce smaller radii). In our study, we find that most SN II progenitors are consistent with RSG radii between 450 and 600 R_⊙ (Fig. 2, bottom panel). The minimum progenitor radius in our grid of pre-SN models is 450 R_⊙ (using a mixing-length parameter of two). Smaller radii can be obtained with a mixing-length parameter of three, as previously suggested by Dessart et al. (2013) based on the analysis on SN 1999em. Adopting a larger mixing-length parameter, we would obtain pre-SN structures with similar pre-SN masses but smaller radii. A smaller pre-SN radius for the same initial mass might help to solve the IMF incompatibility as more massive progenitors are needed to derive our radius distribution. However, a direct comparison is not possible given that larger M_ZAMS progenitors also imply larger ejecta masses and the LC models will dramatically change. Therefore, it would be necessary to carry out a series of simulations before being able to conclude whether or not changes in this parameter could be a determining factor in the IMF incompatibility.

Overshooting. A more extended convective zone produces a more massive helium core, and therefore, higher luminosities. Because of the higher luminosity, the mass-loss rate is larger producing smaller final hydrogen-rich envelope masses. In principle, this can help to solve the discrepancy, but larger overshooting parameters produce much more extended RSGs with dense envelopes (Wagle et al. 2020), acting in the opposite direction to our results. A detailed future analysis should be carried out to fully understand these effects.

7.2.3. Failed SNe II

In the previous sections, we give possible explanations for the IMF incompatibility, but always assuming that the mass function of SN II progenitors should be consistent with a Salpeter massive-star IMF. However, these two distributions may be intrinsically different. Here, we test whether the lack of massive progenitors might be caused by high-M_ZAMS progenitors failing to produce a SN explosion, where their collapse leads to an implosion into a black hole. If the shock produced by the core collapse fails to explode the RSG star, it may produce a low-luminosity red transient lasting for about a year (Lovegrove & Woosley 2013). Such weak explosions would be difficult to detect and may contribute to the lack of high-M_ZAMS progenitors.

We constructed a SN II progenitor IMF using the results from Sukhbold et al. (2016), where the authors calculated explosion models for a large and fine grid of pre-SN models and distinguish between exploding and non-exploding cases. The procedure is similar to that presented in Sect. 6.1, but with one exception. Here, at each MC simulation, we randomly sampled a certain number of stars in the range from 9 to 25 M_⊙ (the lower and upper M_ZAMS of our pre-SN models) assuming a Salpeter massive-star IMF with the condition that the sample has 24 objects (i.e. the number of SNe II in the gold sample) after removing the non-exploding stars. We determined the probability distribution of the CD of the M_ZAMS for the SN II progenitors. We fitted this CD in the same manner as in Sect. 6.1 and found the following parameters: M_low = 9.1${}_{-0.2}^{+0.1}$ M_⊙, M_high = 20.8${}_{-2.3}^{+3.0}$ M_⊙, and Γ = −3.32${}_{-0.68}^{+0.71}$. These results imply that if failed SNe II are taken into account, the upper M_ZAMS limit and the slope of the SN II progenitor IMF are modified. A steeper power-law than that of a Salpeter IMF is found, although not as steep as the slope inferred for the gold sample. In addition, for the case of M_high, a lower value than that assumed for the Salpeter IMF is obtained.

7.2.4. Modelling

Here we briefly discuss if the differences found in the IMF can be attributed to the simplifications in our hydrodynamical simulations, or to a bias in our fitting procedure.

Hydrodynamical simulations. While the details of the explosion mechanism are still under debate (see Burrows et al. 2020, for a recent review), observations (polarisation, SN remnants, and nebular spectroscopy) show that the explosion is asymmetric in the inner regions (e.g. Hoflich 1991; DeLaney et al. 2010; Taubenberger et al. 2009; Tanaka et al. 2012). However, our hydrodynamic code assumes spherically symmetric explosions, which might produce systematic effects (that would be inherent to all 1D codes).

In our explosion models, the explosion energy is deposited at the inner layers by artificially adding internal energy in form of a thermal bomb. The explosion energy is used to overcome gravity and lifts the envelope out of its potential. The high initial explosion energy produces a strong shock front and pushes the material above resulting in the ejection of the envelope. Inner low-velocity and low-density regions expanding with a certain minimum velocity appear. However, such regions are artefacts of the 1D nature of the modelling. In a spherical explosion, a hot low-density and fast region below a cold and denser stellar structure is formed. Thus, this region is highly RT unstable, with instabilities growing on timescales of seconds to minutes (Janka & Müller 1993). In multi-dimensional calculations, RT instabilities mix material filling the inner low-density region. The inner low-density and low-velocity region present in 1D codes might produce systematic errors in the physical parameter estimations, although a detailed analysis of this effect in our results is beyond the scope of this work.

In addition, RT instabilities cause a large-scale mixing of heavy elements into the hydrogen-rich envelope when the shock wave reaches the composition interfaces between the carbon-oxygen core and the helium core, and the helium core and the hydrogen-rich envelope (Mueller et al. 1991; Herant & Benz 1991; Kifonidis et al. 2003), and also the mixing of hydrogen inside the helium core (see Utrobin & Chugai 2008, and references therein). Such mixing of hydrogen inwards is not taken into account in our simulations. However, as shown in previous studies, this causes a minor effect on the LC morphology (Utrobin & Chugai 2008). Additionally, it has been proposed that RT instabilities could also modify the pre-SN density profile (Utrobin & Chugai 2008), leading to a smoother transition between zones of different chemical abundances. As this process happens before the arrival of the shock wave to the surface, it could be taken into account simply by modifying the pre-SN model assumed in the 1D hydrodynamic simulation. The above discussion is the reason why some authors have decided not to use stellar evolution models as initial configurations in 1D simulations. However, it is difficult to conclude whether this effect could be the reason for the IMF incompatibility.

Another important factor that can impact our results is found in the opacity calculation. Our hydrodynamical code adopts a minimum value of the opacity sometimes referred to as the ‘opacity floor’ to account for the actual contribution of the line opacity and the non-thermal excitation or ionisation of electrons that are not included in our opacity calculation (see more details in Bersten et al. 2011). Different values of the opacity floor produce differences in our LC and photospheric velocity models that could modify our results. This could be a systematic source of error that needs to be studied. However, the adopted value for the opacity floor was selected based on several tests (see Bersten et al. 2011; Bersten 2013) with the STELLA code (Blinnikov et al. 1998) that uses multi-group radiative transfer and includes the effect of line opacities (Sorokina & Blinnikov 2002).

Fitting technique. Another possibility is that our fitting procedure could bias our results to low-M_ZAMS progenitors because of the degeneracy between physical parameters, meaning that explosions with different physical properties can produce similar photometric and spectroscopic properties (Dessart & Hillier 2019; Goldberg et al. 2019; Martinez & Bersten 2019). However, in M20, we derived the progenitor and explosion parameters for eight SNe II using the same fitting method and found very good agreement with the progenitor initial mass estimations from pre-SN imaging and nebular spectral modelling. Furthermore, in M20, we compared SNe II observations with several high-mass models, finding that no such model could reproduce the observations (see M20, their Fig. 11). We note that the eight SNe II from M20 have progenitor and observed properties within the ranges derived in this study and Paper I, respectively. We also propose several high-mass models as observations to test our procedure, recovering the high-mass models in each case (see Appendix C). We conclude that our fitting procedure is not biased against high-mass solutions.

In this section we have discussed several possible explanations for the IMF incompatibility. It is clear that the uncertainties in stellar evolution can have a significant impact on the final structure of progenitors, and therefore, on LCs and spectra (see e.g. Dessart et al. 2013). From the present discussion, we conclude that higher mass-loss rates have a clear effect going in the direction needed to explain the IMF incompatibility, although we cannot rule out other possible factors. Additional mass stripping during evolution may be required to obtain similar ejecta masses as those from this study, but for higher-M_ZAMS progenitors. Lower masses at the time of explosion can be reached by a higher wind efficiency (although see Beasor et al. 2020), unstable and episodic mass-loss events, stellar rotation, and/or mass transfer in binary systems. Such mass transfer may be the main mechanism for stripping larger amounts of envelope masses since most massive stars are found in interacting binary systems (Sana et al. 2012). Higher M_ZAMS also imply larger progenitor radii, which modifies the pre-SN structure. Therefore, a larger mixing-length parameter may be more appropriate to use, as suggested by Dessart et al. (2013) and González-Gaitán et al. (2015). To test this conclusion, we carried out a smaller number of additional simulations of high-M_ZAMS progenitors evolved with a larger mixing-length parameter that experienced an increased mass loss. The analysis of these models is presented in Appendix D, and suggests that indeed such models may also provide reasonable fits to some of our SNe II. A detailed study in this direction is clearly warranted in the future.

8. Summary and conclusions

In this, the second paper of a series, we have derived the physical properties (M_ZAMS, M_ej, R, E, M_Ni, and ⁵⁶Ni mixing) from the hydrodynamical modelling of bolometric LCs and expansion velocities for a sample of 53 SNe II. This is the largest sample of SNe II analysed so far that comes from the same follow-up programme where the SNe were all observed and processed in the same manner. We used a grid of explosion models applied to stellar evolution progenitors and compared them to SN observations using a quantified fitting technique based on MCMC methods. The pre-SN models were constructed assuming non-rotating solar-metallicity single-star evolution with standard values for various evolutionary parameters (efficiency of mass loss through winds, mixing length, overshooting, and semiconvection). Therefore, our conclusions are valid for standard pre-SN models, and variations in stellar modelling may provide different results.

Using progenitors produced by standard stellar evolution, we were able to model the majority of the objects in the sample, with the exception of a few SNe II that will be analysed in Paper III. The following ranges of physical parameters were estimated, M_ej = 7.9−14.8 M_⊙, R = 450−1077 R_⊙, E = 0.15−1.40 foe, and M_Ni = 0.006−0.069 M_⊙, with different degrees of ⁵⁶Ni mixing within the ejecta, from a more centrally concentrated ⁵⁶Ni to a more extended spatial distribution (with a median value of ∼0.5 in mass fraction). An extreme mixing of ⁵⁶Ni was only found for one SN in the gold sample: SN 2008ag. Values of M_ZAMS are found between 9.2 and 20.9 M_⊙. In summary, we find ranges of progenitor and explosion properties of SNe II that qualitatively agree with previous observational and theoretical estimates. Five SNe II in the CSP-I sample with well-sampled bolometric LCs (and sometimes well-sampled velocity curves) could not be fully reproduced with our grid of hydrodynamical simulations applied to standard single-star models, implying that non-standard stellar evolution and/or binary evolution models are required¹¹.

We analysed the distributions for all the physical parameters and found that most SNe II in our sample are consistent with ejecta masses in the range of 8−10 M_⊙ and RSG radii between 450 and 600 R_⊙, which are the result of the evolution of low-M_ZAMS progenitors between 9 and 12 M_⊙. The study of the M_ZAMS distribution for both the gold and full samples, assuming a Salpeter massive-star IMF or a free power-law slope, suggests a lack of SN II progenitors more massive than 21.5 M_⊙, consistent with the RSG problem at different levels of significance. This is the first time that the RSG problem is recovered from LC and photospheric velocity modelling. The gold sample gives an upper mass boundary of M_high = 21.3${}_{-0.4}^{+3.8}$ M_⊙ and a power-law slope of Γ = −6.35${}_{-0.57}^{+0.52}$. This power law is not consistent with previous IMF estimations for massive stars. The steeper IMF found for the SN II progenitors is caused by the large relative number of low- to high-ejecta mass that biases the M_ZAMS distribution to lower masses. We discussed several possible explanations for this discrepancy, suggesting that the bias in our results towards low M_ZAMS may be caused by the pre-SN structures utilised. Our analysis indicates that similar pre-SN structures, for higher M_ZAMS, may be obtained by considering a higher rate of mass loss during massive-star evolution, probably as a result of binary interaction prior to the RSG phase, and larger mixing-length parameters. However, a systematic analysis of the different factors needs to be carried out to properly identify the reason for the IMF incompatibility.

The progenitor and explosion properties inferred in this work are used in Paper III to search for correlations between physical parameters and observed bolometric and spectroscopic properties for the same sample of SNe II. This provides an extensive analysis of SN II diversity dependence on the physics of the progenitors and their explosions.

Acknowledgments

We thank the referee for the useful comments that improved the manuscript. The work of the Carnegie Supernova Project was supported by the National Science Foundation under grants AST-0306969, AST-0607438, AST-1008343, AST-1613426, AST-1613472, and AST-1613455. L.M. acknowledges support from a CONICET fellowship. L.M. and M.O. acknowledge support from UNRN PI2018 40B885 grant. M.H. acknowledges support from the Hagler Institute of Advanced Study at Texas A&M University. S.G.G. acknowledges support by FCT under Project CRISP PTDC/FIS-AST-31546/2017 and Project No. UIDB/00099/2020. M.S. is supported by grants from the VILLUM FONDEN (grant number 28021) and the Independent Research Fund Denmark (IRFD; 8021-00170B). F.F. acknowledges support from the National Agency for Research and Development (ANID) grants: BASAL Center of Mathematical Modelling AFB-170001, Ministry of Economy, Development, and Tourism’s Millennium Science Initiative through grant IC12009, awarded to the Millennium Institute of Astrophysics, and FONDECYT Regular #1200710. L.G. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under the 2019 Ramón y Cajal program RYC2019-027683 and from the Spanish MICIU project PID2020-115253GA-I00. P.H. acknowledges the support by National Science Foundation (NSF) grant AST-1715133. Software: corner.py (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), NumPy (Oliphant 2006; Van Der Walt et al. 2011), matplotlib (Hunter 2007), MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), SciPy (Virtanen et al. 2020), Pandas (McKinney 2010), ipython/jupyter (Perez & Granger 2007).

References

Appendix A: Influence of the gold sample selection on the results

As described in detail in Sect. 5, a gold SN II requires both good data coverage and good fits to our models. The classification of the gold sample was achieved visually by the first three authors of this paper. In most cases those authors agreed in the classification, but a small number of events were controversial. In this section, we analyse whether the inclusion or exclusion of some debatable cases affect our results on the IMF parameters.

We started by including three additional SNe II to the gold sample reaching a total of 27 events. These SNe II were not initially included in the gold sample for the following reasons: SN 2005dw presents two Fe II line velocity measurements but very close in time; SN 2007hm does not show a clear end of the plateau phase; and SN 2008K shows offsets in the fits to both bolometric LC and photospheric velocities. The inclusion of these SNe II in the gold sample leads to the following IMF parameters: M_low = 9.3${}_{-0.1}^{+0.1}$ M_⊙, M_high = 24.5${}_{-3.0}^{+4.3}$ M_⊙ ($M_{\mathrm{high}}^{\mathrm{mode}}$ = 21.0${}_{-0.6}^{+3.4}$ M_⊙), and Γ = −6.48${}_{-0.52}^{+0.51}$. We next excluded three events from the authentic gold sample defined in Sect 5. SN 2006be and SN 2007W have observations at the beginning of the transition from the plateau to the tail phase, although it might not be clear. In addition, we excluded SN 2008bk, which does not fulfil the conditions to be considered a gold event since there is no CSP-I photometry during the transition from the plateau to the radioactive tail phase (published data are used to restrict the end of the plateau). With the exclusion of these three SNe II from the gold sample, the results are: M_low = 9.3${}_{-0.1}^{+0.1}$ M_⊙, M_high = 25.1${}_{-3.5}^{+4.0}$ M_⊙ ($M_{\mathrm{high}}^{\mathrm{mode}}$ = 22.7${}_{-1.7}^{+2.8}$ M_⊙), and Γ = −6.39${}_{-0.64}^{+0.65}$. In conclusion, this analysis suggests that even if those limiting cases are included/excluded from the gold sample, the parameters derived for the IMF are not statistical different with the results presented in Sect. 6.1.

Appendix B: Physical parameters of SNe II

Figures B.1, B.2, B.3, and B.4 present bolometric LCs and Fe II velocities for the 53 SNe II included in this study in comparison with the models drawn from the posterior distribution of the parameters. Fe II velocities are not available for SNe 2005af, 2005lw, and 2007it. Table B.1 presents the physical parameters derived from the hydrodynamical modelling of SN II bolometric LCs and expansion velocities. In Table B.1, we also include the estimates of M_ej, M_H, env, M_H, and R. We note that these progenitor properties are not model parameters and, therefore, they are not fitted. These values are linearly interpolated from the M_ZAMS derived from the fitting. Additionally, an example corner plot with the posterior probability distribution of the parameters is presented in Fig. B.5.

Fig. B.1. Comparison between models and observations. The open circles show the observed bolometric LCs (left) and Fe II velocities (right). Solid lines represent 30 models randomly chosen from the probability distribution. The panels present SNe in order of their discovery dates, from SN 2004ej to SN 2005dt. Grey shaded regions show the early data we removed from the fitting. The errors in the observed bolometric LCs are not plotted, for better visualisation.

Fig. B.2. Same as in Fig. B.1, but from SN 2005dw to SN 2008K.

Fig. B.3. Same as in Fig. B.1, but from SN 2008M to SN 2009N.

Fig. B.4. Same as in Fig. B.1, but from SN 2009ao to SN 2009bz.

Fig. B.5. Corner plot of the joint posterior probability distribution of the parameters for SN 2006ai.

Appendix C: High-mass progenitor fits

We mapped LCs and expansion velocities from high-mass models as observations to test our fitting technique. Figure C.1 shows a corner plot of the joint posterior probability distribution of the parameters for a model of M_ZAMS = 24 M_⊙, E = 1.4 foe, M_Ni = 0.03 M_⊙, and ⁵⁶Ni mixing = 0.5. The blue squares represent the input physical parameters of the model.

Fig. C.1. Corner plot of the joint posterior probability distribution of the parameters when a high-mass progenitor model is considered as observations for our fitting technique. Blue squares represent the input physical parameters of the model.

Appendix D: High-mass progenitor models evolved with enhanced mass loss

In Sect. 7.2 we proposed that one possible solution for the IMF incompatibility could be to consider a higher rate of mass loss, in addition to a larger mixing-length parameter. A detailed investigation is beyond the scope of this work; however, we give an example to show that our suggestion allows high-M_ZAMS progenitors to fit some SNe II in the CSP-I sample.

We used MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) to obtain pre-SN models employing the same configurations detailed in Sect. 4, with two exceptions. We adopted a wind scaling factor of η = 1.5 (i.e. 50% higher than standard) and a mixing-length parameter of three. Then, we used the code presented in Bersten et al. (2011) to obtain bolometric LC and photospheric velocity models. We generated a small grid of models for M_ZAMS = 20 and 22 M_⊙, E = 0.8, 1.0, and 1.2 foe, M_Ni = 0.03 M_⊙, and ⁵⁶Ni mixing = 0.8 (given in fraction of the pre-SN mass). Thus, we obtained a total of six hydrodynamical models.

We used these six models to compare with all CSP-I SNe II with bolometric luminosities during the photospheric phase and – at least – the transition to the radioactive tail and velocity measurements. The observations of a small number of SNe II in the CSP-I sample can be reasonably reproduced by some of these models. We discuss the case of SN 2006ai as an example. The 20 M_⊙ models reproduce the plateau phase and the transition to the radioactive phase. The models display lower luminosities during the radioactive tail phase, indicative of a higher ⁵⁶Ni mass in the ejecta of SN 2006ai than that assumed in the models. For this reason, we carried out new models with higher ⁵⁶Ni masses until the observed radioactive tail phase is reproduced by our models. Figure D.1 shows the models that reproduce the observations of SN 2006ai¹². Both models assume M_ZAMS = 20 M_⊙, M_Ni = 0.049 M_⊙, and ⁵⁶Ni mixing = 0.8. Model 1 and Model 2 were exploded with E = 1.0 and 1.2 foe, respectively. We stress that this test seeks reasonable agreement between models and observations and not detailed fits to observations. A fine tuning of the model parameters is required to find the best-fitted model.

Fig. D.1. Observed bolometric LC (top panel) and Fe IIλ5169 line velocities of SN 2006ai (black markers). Solid lines represent the high-MZAMS models with increased mass loss that best reproduce the observations. Dashed lines represent the maximum a posteriori model using our standard grid of models. This model is corrected by the scale factor and the explosion epoch found in Sect. 5. Since our estimation suggests an earlier explosion epoch, the model shown as dashed lines starts before zero on the x axis (where zero corresponds to the explosion epoch derived by Gutiérrez et al. 2017b).

This test shows that a high-M_ZAMS progenitor evolved with enhanced mass loss and larger mixing-length parameter can reproduce the observations of a SN II (i.e. in the direction to solve the IMF incompatibility). However, a deeper analysis is necessary for stronger conclusions.

All Tables

All Figures

Fig. 1. Observed bolometric LC (top panel) and Fe IIλ5169 line velocities (bottom panel) of SN 2006ai (filled markers) with models randomly sampled from the posterior distribution of the parameters (solid lines). The median value and 68% confidence range for every physical parameter are also shown. The grey shaded region shows the early data removed from the fitting (first 30 days after explosion).

In the text

Fig. 2. Histograms of three progenitor parameters: MZAMS (top panel), Mej (middle panel), and R (bottom panel). The gold sample is represented by yellow bars, and the red bars are the histograms for the full sample. In each panel, the number of SNe II is listed, together with the median and the 16th and 84th percentiles. Probability density functions of the physical parameters for the gold and full samples are represented by solid and dashed lines, respectively.

In the text

	Fig. 3. Same as in Fig. 2 but for E (top panel), MNi (middle panel), and 56Ni mixing (bottom panel).
In the text

	Fig. 4. Cumulative distribution of MZAMS for the SN II progenitors in the gold sample. The derived masses are shown in yellow squares. The shaded contours show the confidence regions of the CD.
In the text

	Fig. 5. Same as in Fig. 4 but for the SN II progenitors in the full sample. Yellow squares represent the SNe II in the gold sample.
In the text

Fig. 6. Fits to the CD of MZAMS of SN II progenitors in the gold sample. Left panel: CD of MZAMS derived in our work for the gold sample (yellow triangles and error bars), corresponding to the median and the 68% confidence limit of the distribution presented in Fig. 4 as shaded contours, in comparison with the model CD constructed with the median values of the marginal distributions of the parameters (shaded contours). Right panel: corner plot of the joint posterior probability distribution of the parameters. Dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions.

In the text

	Fig. 7. Same as in Fig. 6 but when the power-law slope is fixed to Γ = −2.35, that is, the slope of a Salpeter massive-star IMF.
In the text

Fig. 8. Fits to the CD of MZAMS of SN II progenitors in the full sample. Top left: CD of MZAMS derived in our work for the full sample (red triangles and the error bars), corresponding to the median and the 68% confidence limit of the distribution presented in Fig. 5 as shaded contours, in comparison with the model CD constructed with the median values of the marginal distributions of the parameters (shaded contours). Bottom left: corner plot of the joint posterior probability distribution of the parameters. Dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions. The right panels show the same as the left panels but when the power-law slope is fixed to Γ = −2.35.

In the text

	Fig. 9. Comparison of ejecta masses (top panel) and explosion energies (bottom panel) derived in the current work to those from previous studies containing large samples of SNe II.
In the text

	Fig. 10. Best-fitting power-law slope of the CD of the progenitor masses for a range of sample sizes. The error bars correspond to the 68% confidence limit. The larger the sample, the higher the probability that the best-fitting power-law slope is that of a Salpeter IMF. The solid black line indicates the Salpeter IMF slope for massive stars.
In the text

Fig. B.1. Comparison between models and observations. The open circles show the observed bolometric LCs (left) and Fe II velocities (right). Solid lines represent 30 models randomly chosen from the probability distribution. The panels present SNe in order of their discovery dates, from SN 2004ej to SN 2005dt. Grey shaded regions show the early data we removed from the fitting. The errors in the observed bolometric LCs are not plotted, for better visualisation.

In the text

	Fig. B.2. Same as in Fig. B.1, but from SN 2005dw to SN 2008K.
In the text

	Fig. B.3. Same as in Fig. B.1, but from SN 2008M to SN 2009N.
In the text

	Fig. B.4. Same as in Fig. B.1, but from SN 2009ao to SN 2009bz.
In the text

	Fig. B.5. Corner plot of the joint posterior probability distribution of the parameters for SN 2006ai.
In the text

	Fig. C.1. Corner plot of the joint posterior probability distribution of the parameters when a high-mass progenitor model is considered as observations for our fitting technique. Blue squares represent the input physical parameters of the model.
In the text

Fig. D.1. Observed bolometric LC (top panel) and Fe IIλ5169 line velocities of SN 2006ai (black markers). Solid lines represent the high-MZAMS models with increased mass loss that best reproduce the observations. Dashed lines represent the maximum a posteriori model using our standard grid of models. This model is corrected by the scale factor and the explosion epoch found in Sect. 5. Since our estimation suggests an earlier explosion epoch, the model shown as dashed lines starts before zero on the x axis (where zero corresponds to the explosion epoch derived by Gutiérrez et al. 2017b).

In the text