Introduction

Massive stars (>8 M) have a significant impact on the interstellar medium (ISM) by regulating it through stellar winds, ionising radiation, and supernova (SN) explosions. SNe contribute to the origin of heavy elements, a still poorly understood aspect in cosmology1. The details of SN feedback and metal production depend primarily on which stars explode as which type of SNe because the mechanism of explosion and element production yield is different for each of them. SN feedback models that take into consideration progenitor stars and mechanisms of explosions are essential to improve simulations.

The connection between SN types and their progenitors is of great importance for galaxy evolution and cosmology. This can be directly studied by identifying progenitor stars on pre-explosion images, with follow-up observations that confirm their disappearance. However, to date, only 23 (18 Type II, four Type IIb, and one Type Ib) core-collapse SN (CCSN) progenitor stars have been confirmed to disappear in post-explosion images2. No Type Ic SN (without hydrogen or helium in the spectrum) progenitor has been confirmed in this way, with SN 2017ein as the only candidate, but with a wide range of derived progenitor masses3,4,5,6. Hence, most of our knowledge on the nature of Type Ic SN progenitors is based on photometric and spectroscopic observations after the explosion.

Type Ic SNe likely originates either from core-collapses of very massive stars or from less massive stars in binary systems. The mechanism responsible for the lack of hydrogen or helium lines is still a subject of debate. In the case of the very massive star model, the hydrogen and helium envelope is lost due to stellar winds from the star itself7,8 (e.g., radiatively driven winds, episodic mass-loss, or rapid rotation). In the binary model, the companion radically affects the evolution of the progenitor due to the mass exchange9,10,11. Non-detections of Type Ic SN progenitors in pre-explosion images and relative rates of different types of SNe suggest that most of Type Ic SN progenitors are binaries with initial masses  <  20 M12,13.

Another way to address the relation between progenitor stars and resulting SNe is to investigate the molecular gas properties at the explosion ___location. Molecular gas at the locations of SNe of different types was recently investigated at a spatial resolution comparable to giant molecular clouds (GMCs)14, using the Atacama Large Millimetre Array (ALMA) carbon monoxide 2-1 line transition [CO(2-1)] observations from the Physics at High Angular resolution in Nearby GalaxieS (PHANGS)15,16 survey. Their sample consisted of a total of 59 SNe: 12 thermonuclear (Type Ia SNe), 32 Type II SNe, eight stripped-envelope SNe (SESNe, hereafter, Type Ib, Ic or Ib/c), and seven unclassified. They found that Type Ia and II SNe are associated with little or no molecular gas emission, while SESNe and unclassified SNe mostly show strong molecular gas emission. They concluded that there is a clear dependence of the type of SN and the molecular gas environment, however, their conclusions are drawn based on a low sample size for SESNe and, thus are not statistically significant.

In this work, our goal is to constrain lifetimes and initial masses of Type Ic SN progenitors. To this end, we compare the molecular gas densities at the positions of Type II and Ic SNe. By targeting a large sample of SNe we aim to uncover their nature. This statistical approach offers strong constraints on the overall progenitor characteristics of different SN populations but does not provide a strong constraint on individual SN progenitor properties. We report a statistically significant study to do so with spatial resolution comparable to the GMC sizes. This is an important factor because molecular hydrogen column surface densities and lifetimes of GMCs can only be measured accurately if the resolution at least matches the cloud sizes17,18.

Results

In order to investigate the environments of a significant number of SNe at high resolution, we initiated the ALMA CO SN (ACOS) survey, obtaining CO(2-1) observations of the locations of 16 Type Ic SNe. Together with the PHANGS survey this results in a sample of 63 SNe: 12 Type Ia, 30 Type II, and 21 Type Ic SNe. These CO(2-1) observations have a spatial resolution of ~ 100 pc, similar to sizes of GMCs. The spatial resolution and the large sample allow us to study the immediate environments in which the SNe exploded. Our main conclusion is from the comparison of Type II and Ic SNe, whereas Type Ia SNe are shown only to contrast the different progenitor natures. See 'Methods', subsections ALMA CO SN survey, PHANGS–ALMA data, and Supernova sample for detailed information about the SNe and their host galaxies. Supplementary Data 1 lists the information for the SNe used in this work.

As an example, Fig. 1 shows the molecular hydrogen gas surface density (Σmol) map of NGC 4254 (M 99, a typical PHANGS–ALMA galaxy) with its four CCSNe; three Type II (SN 1967H, SN 1972Q, and SN 1986I) and one Type Ic (SN 2014L) SNe, plus the ___location of SN 2009bb (Type Ic SN) hosted in NGC 3278 (observed by ACOS). Σmol is computed from the CO(2-1) line intensity, see methods, subsection CO-to-Σmol for the description used.

Fig. 1: Distribution of molecular gas in galaxy host and SN locations.
figure 1

SN 1967H, SN 1972Q and SN 1986l hosted in NGC 4254 (a), and SN 2009bb hosted in NGC 3278 (b). Colour-coded Σmol intensity is represented in logarithmic scale. Pixels without signal are masked and shown as white. Red dots represent SN locations. North is up and East is left.

The Σmol value for a given SN was calculated in two ways. First, we measured it at the exact pixel of the SN ___location and this is denoted as “SN ___location”. However, the SN ___location might not be the exact site where the progenitor star formed. The true ___location of the formation of the SN progenitor star could be shifted due to astrometric displacement and/or peculiar motion of the progenitor system with respect to the parent GMC. To take into account these effects, together with the maximum size of GMCs, we also measured the maximum value within a radius of 200 pc centred at the SN position (see methods, subsection 200 pc regions). To assess an average level of Σmol for our host galaxies, we also measured Σmol for 103 random pixels in each PHANGS galaxy and the respective maximum value within a radius of 200 pc centred at these random pixels. The ACOS sample was not included in the random pixel calculations because these observations cover the SN positions but not the whole galaxy, so it is not possible to obtain a representative measurement of Σmol for these galaxies.

Figure 2 shows the empirical cumulative distribution function (eCDF) of Σmol in the environments of SNe with their respective confidence intervals calculated using the Dvoretzky–Kiefer–Wolfowitz inequality19 at a significance level of 68%. If a given Σmol measurement has a significance lower than 2σ (non-detection), where σ is the noise at the pixel ___location, then the 2σ upper limit is shown on the eCDF. This happens for 66% of the random pixels (15077 out of 23000), 67% of Type Ia (eight out of 12), 37% Type II (11 out of 30) and 43% Type Ic SNe (nine out of 21). These upper limits are not used for any calculations, only for the visualisation of the eCDFs. Supplementary Data 2 and 3 lists the Σmol environments of SNe and random locations in galaxy hosts, respectively. Random pixel locations extend towards lower values when compared to the SN positions, indicating that the SN locations have a strong connection with their GMC parents or larger molecular structures. This is an expected result for CCSNe, since they are associated with massive progenitor stars and therefore they are expected to explode close to their birthplaces in spiral arms where molecular gas densities are higher20.

Fig. 2: Σmol eCDFs for SN locations.
figure 2

Random locations, Type Ia, Type II, and Type Ic SNe are represented by black, blue, red, and grey lines, respectively. Upper limits (2σ) are used in case of non-detections. The shaded areas represent confidence intervals at 1σ.

To quantify if our samples are drawn from different parent GMC populations, we performed a Kolmogorov—Smirnov (KS) test for each combination. KS statistics and p-values are shown in Table 1 (and for the 200 pc regions in Supplementary Table 1). Comparing the locations of CCSNe and random pixels, the low p-values indicate that their distributions are different, suggesting high densities of molecular gas environments for the sites where SNe was observed. On the other hand, high p-values (> 0.05) for SN combinations indicate that it is not possible to reject the null hypothesis that the Σmol at the SN positions are from the same distribution.

Table 1 KS statistics (and p-value in parenthesis) for different SN groups (SN locations)

In order to obtain confidence interval ranges of Σmol for each SN type and random locations, we performed Monte Carlo simulations. Each measurement of an individual SN was perturbed according to its uncertainty, which was assumed to have a Normal distribution. We created 104 of such perturbed sets, each time calculating its median value. The uncertainty was adopted to be at 16% and 84% of the simulated distribution. We show the results of these simulations in Fig. 3. The medians and 1σ confidence intervals for SN locations and 200 pc regions are summarised in Supplementary Table 2.

Fig. 3: Median values and 1σ confidence intervals (using 104 Monte Carlo simulations) of Σmol for SN locations.
figure 3

Random locations, Type Ia, Type II, and Type Ic SNe are represented by black, blue, red, and grey error bars, respectively.

Median molecular gas densities increase from the values measured for the random pixels (\(4.4{7}_{-0.04}^{+0.05}\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2}\)), through Type Ia SNe (\(6.9{3}_{-2.36}^{+3.70}\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2}\)), to Type II and Ic SNe (\(20.1{5}_{-2.46}^{+3.38}\) and \(20.6{2}_{-4.88}^{+4.28}\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2}\), respectively). The results are shown in Fig. 3. The random locations and Type Ia SNe have much lower median molecular gas densities than CCSNe. At the positions of Type II and Ic SNe we obtained similar molecular gas densities within the 1σ confidence levels.

Under the assumption of single very massive star progenitors for Type Ic SNe (for masses of 30–100 M with lifetimes of 7–3 Myr, respectively21), it is expected that the respective parent GMC would not have been dispersed before the SN explosion due to a short progenitor lifetime, and the progenitors would not have enough time to shift away from their birthplaces significantly. Therefore, if progenitors of Type Ic SNe were very massive stars, then the distribution of Σmol at their positions should be shifted toward higher values compared to that of Type II SNe. This is because the lower masses of progenitors of Type II SNe imply longer lifetimes, and therefore more time for the parent clouds to disperse and for the progenitors to move away. This scenario is not supported by our results. In the alternative scenario, the progenitors of Type II SNe evolve as single stars (or in wide binaries in which their hydrogen layers are not affected) and those of Type Ic SNe are similarly massive stars that evolve in binary systems with a companion being responsible for removing the external layers of hydrogen and helium21. Then the progenitor masses, and therefore lifetimes, of both types are similar. Thereby, their distributions of Σmol should be comparable, which is consistent with our data. See methods, subsection Timing the SN progenitor lifetime with molecular gas information for the justification of using the molecular gas densities to constrain the stellar population age.

We note that at lower gas densities the distribution of Type Ic SNe is slightly lower than that of Type II SNe (see Fig. 2). This may indicate that the lifetimes of some of the Type Ic progenitors were increased by the binary interaction21. However, the statistical significance of this difference is too low to draw any definitive conclusions.

In order to assess the maximum difference between the lifetimes of progenitors of Type II and Ic SNe, we assumed that the molecular gas density at the SN progenitor positions, \({\Sigma }_{{{{{\rm{mol,SN}}}}}}\), decreases exponentially with time as

$${\Sigma }_{{{{{\rm{mol,SN}}}}}}={\Sigma }_{{{{{\rm{0}}}}}}\,{e}^{-{t}_{{{{{\rm{SN}}}}}}/{\tau }_{{{{{\rm{GMC}}}}}}},$$
(1)

where Σ0 is a normalisation constant, \({t}_{{{{{\rm{SN}}}}}}\) is the lifetime of the SN progenitor and τGMC is the lifetime of the GMC. In order to constrain the lifetimes of Type Ic SN progenitors, we made use of our measured ratio of \({\Sigma }_{{{{{\rm{mol,SN}}}}}}\) for Type Ic and II SNe (\({\Sigma }_{{{{{\rm{mol,Ic}}}}}}/{\Sigma }_{{{{{\rm{mol}}}}},{{{{\rm{II}}}}}}=1.0{2}_{-0.27}^{+0.27}\)) and a measured value of the GMC lifetime of \({\tau }_{{{{{\rm{GMC}}}}}}=1{6}_{-5}^{+6}\,{{{{\rm{Myr}}}}}\)22.

This characteristic cloud evolution lifetime is in agreement with theoretical works and simulations18,22,23,24,25,26,27 and while the exact value of the assumed average GMC lifetime influences these calculations, it does not change the interpretation when lifetimes of two SN types are compared. Assuming that progenitors of Type II and Ic SNe are born in GMCs with similar average initial conditions, i.e. that average Σ0 is the same for both, from eq. (1) it is possible to calculate the difference between the SN progenitor lifetimes as

$${t}_{{{{{\rm{II}}}}}}-{t}_{{{{{\rm{Ic}}}}}}={\tau }_{{{{{\rm{GMC}}}}}}\,\ln \left({\Sigma }_{{{{{\rm{mol,Ic}}}}}}/{\Sigma }_{{{{{\rm{mol}}}}},{{{{\rm{II}}}}}}\right)=0.3{7}_{-4.26}^{+4.27}\,{{{{\rm{Myr}}}}}.$$
(2)

A zero-age main sequence (ZAMS) mass \({M}_{{{{{\rm{init,SN}}}}}}\) for Type II SN progenitors of \({M}_{{{{{\rm{init,II}}}}}}=10.6{6}_{-0.20}^{+0.20}\,{{{{{\rm{M}}}}}}_{\odot }\) (obtained by averaging nine SNe with pre-explosion detection28 and confirmed by the disappearance in post-explosion images2) yields a lifetime of tII = \(25.2{2}_{-0.80}^{+0.80}\,{{{{\rm{Myr}}}}}\). Using this tII in eq. (2) results in a lifetime for Type Ic SN progenitors of \({t}_{{{{{\rm{Ic}}}}}}=24.{9}_{-4.3}^{+4.3}\,{{{{\rm{Myr}}}}}\) and a ZAMS mass of Minit,Ic = \(10.9{0}_{-1.20}^{+1.20}\,{{{{{\rm{M}}}}}}_{\odot }\). On the other hand, if we assume a typical mass for red supergiant progenitors for type II SNe29 of \({M}_{{{{{\rm{init,II}}}}}}=1{5}_{-1}^{+1}\,{{{{{\rm{M}}}}}}_{\odot }\), then we obtain Minit,Ic = \(15.{3}_{-3.2}^{+3.2}\,{{{{{\rm{M}}}}}}_{\odot }\). This also means that if Type II SN progenitors include rare examples of very massive stars, so can Type Ic SN progenitors30. To account for signal unrelated to SNe, as the first step we also subtracted the random values of Σmol from those of SNe, which resulted in \({t}_{{{{{\rm{II}}}}}}-{t}_{{{{{\rm{Ic}}}}}}=0.4{7}_{-5.45}^{+5.47}\,{{{{\rm{Myr}}}}}\), \({t}_{{{{{\rm{Ic}}}}}}=24.{8}_{-5.5}^{+5.5}\,{{{{\rm{Myr}}}}}\), and Minit,Ic = \(10.{9}_{-1.5}^{+1.5}\,{{{{{\rm{M}}}}}}_{\odot }\), indistinguishable from the original results.

Another effect to take into account is that SN progenitors may be runaway stars, which are moving away from their parent clusters with significant velocities. Maximum velocities for OB runaway stars are ~ 30 km s−1 or ~ 30 pc/Myr31, so they would need ~ 3 Myr to cross a GMC. Replacing τGMC with the effective crossing timescale τcro, Eq. (1) will then be smaller, including both the movement of the progenitors and the cloud dispersal. Using τcro = 3 Myr, Eq. (2) yields an even smaller lifetime difference between Type II and Ic SNe of ~ 0.06 Myr, what makes our conclusions even stronger.

If most of Type Ic SN progenitors were very massive stars with masses around 30 M, then for their lifetime of 7 Myr, Eq. (1) results in molecular gas densities a factor of 4 higher than those at the positions of Type II SNe, which is not in agreement with our results.

Discussion

Our findings indicate that the binary interaction model (mass transfer due to a companion) is the main mechanism extracting outer layers for most Type Ic SN progenitors. However, we remark that we do not reject the possibility that strong stellar winds of a high-mass star can blow away the layers of hydrogen and helium leading to an explosion as a Type Ic SN, individual Type Ic SNe can be due to very massive star progenitors32. Indeed our calculations show that the 1σ range of the accepted fraction of low-mass stars in the Type Ic SNe population is 66–100% (see methods, subsection Statistical significance of the sample). Another possible binary progenitor scenario includes mergers or accretion of significant mass from a companion so that an initially low-mass star becomes massive and luminous enough to launch strong winds shedding the outer envelope before explosion.

Our results are consistent with low measurements of Type Ic SN progenitor masses from light curve modelling33, the comparison of the SN rates12, the modelling of emission lines at the positions of SNe Type Ib/c29, and the direct observational evidence of a binary system for Type Ic SN 2022jli34. Moreover, in our Galaxy, ~ 70% of massive O-type stars are formed in close binary systems and are expected to experience mass transfer during their lifetime35. All these works support our conclusion about binary systems as progenitors of Type Ic SNe.

On the other hand, our low progenitor masses for Type Ic SNe are inconsistent with some other measurements. The 56Ni yields of Type Ic SNe are five times higher than those of Type II SNe36. However, this may not necessarily imply higher progenitor masses, but different explosion mechanisms or energy sources37. If Type Ic SN explosions are asymmetric or not predominantly powered by the radioactive decay of 56Ni, then the nickel mass will be overestimated for them.

Moreover, using PHANGS data, ref. 14 reported that SESNe (Type IIb, Ib, Ib/c, and Ic SNe) are associated with stronger CO(2-1) emission than those for Type Ia and II SNe. Their sample included 59 SNe in 28 galaxies, with the Type Ic sample size being the main difference compared to this work (they included only five of them, see methods, subsection Statistical significance of the sample for a statistical comparison).

SESNe (including Type Ic SNe) are associated with higher Hα intensities than Type II SNe, suggesting younger star-forming regions and higher progenitor masses38,39,40,41,42,43,44. However, these results are complicated by the fact that the Hα emission disappears on a timescale of several Myr only for isolated star-forming regions, whereas for larger complexes (several hundreds of pc, as probed by these observations), this timescale may be as long as 20 Myr45. Low ages for Type Ic SN locations were also obtained from population fitting of stars in the SN vicinity46,47, but the ages may be underestimated due to blending of stars47. Moreover, stronger associations of Type Ic SNe with Hα and UV-bright stars may also be explained within a scenario in which their lower-mass progenitors prefer regions of high stellar density or more top-heavy initial mass function (IMF), which increase the binarity fraction, but also the Hα and UV emission. High-resolution multi-wavelength observations are needed to alleviate this tension.

Broad-lined Type Ic (Type Ic-BL) SNe and gamma-ray bursts, both known to be connected with very massive progenitor stars48, show twice the amount of synthesised 56Ni than regular Type Ic SNe49. Our sample contains only four Type Ic-BL SNe and when they are excluded from the Type Ic SN sample, the results remain the same (at the 1σ level). Similarly, when Type IIn and IIn/LBV (2 in total) are removed from the Type II SN sample, the results do not change.

The lack of statistically significant difference between Type Ia and II SN positions could be due to a low number of statistics for the former. Subtracting the median molecular density at the random positions from that of Type Ia SNe yields a non-detection: \({\Sigma }_{{{{{\rm{mol}}}}}}=1.9{3}_{-2.36}^{+3.70}\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2}\). Using a corresponding 2σ upper limit in Eq. (2) gives a lower limit on the lifetime difference between the Type Ia and II SN progenitors of tIa − tII > 12 Myr, consistent with the expected higher lifetimes of Type Ia SN progenitors.

Implementation of our findings into cosmological simulations will have important impact on our understanding of the feedback processes. The progenitor type can be implemented in sub-grid prescriptions in numerical cosmological simulations50,51,52 with regards to SN feedback and chemical yields into the ISM. Moreover, whether the progenitors of Type Ic SNe are binary systems or very massive stars changes their contribution to the formation of heavy elements, one of the key aspects of stellar and galaxy evolution1. This includes carbon, one of the most fundamental elements in the Universe and building block of life. At solar metallicity, within the Local Universe (including the Milky Way), stellar winds and SN explosions from binary–stripped stars are found to produce twice more 12C than similar single stars53. The fourth most abundant element in the Solar system is carbon (after hydrogen, helium, and oxygen) and a significant contribution could be produced from Type Ic SNe. This method can also be applied to larger samples divided into different properties (e.g., explosion characteristics, host galaxy types, environmental metallicities, etc) and more rare events to learn about their nature.

Methods

ALMA CO SN survey

ACOS (ALMA ID 2021.1.00099.S, P.I. M.J.M.) consists in observations of the J = 2 → 1 transition of the 12CO line in the environment of 16 Type Ic SNe with an angular spatial resolution range of 0.4–1.1 so that the physical resolution is around 50–100 pc. The selection criteria were the observability with ALMA, i.e. declination < 20°, and distances < 55 Mpc (redshift z < 0.013) to allow the detection of individual GMCs in a reasonable observing time. We excluded seven hosts which are edge-on, for which projection effects would make it difficult to measure the gas surface density at the SN position. The distances, masses and luminosities of these galaxies are comparable to PHANGS–ALMA galaxies.

PHANGS–ALMA data

PHANGS survey provides CO(2-1) line observations using ALMA for 74 galaxies in the Local Universe (<20 Mpc), which mostly are face-on (i < 75o). The typical resolution was ~ 2, corresponding to ~ 100 pc comparable to the sizes of GMCs. We used calibrated data from ref. 15 for the PHANGS–ALMA galaxy sample16. Further information about these procedures can be found in refs. 15,16.

Supernova sample

Our SN sample was compiled from the Open Catalogue for Supernova (https://github.com/astrocatalogs/supernovae/) in April 202254. The SN compilation, as designated in the catalogue, consists of Type Ia (Ia, Ia-02cx, and IaPec), Type II (II, IIP, IIn and IIn/LBV), and Type Ic (Ic and Ic-BL). We did not consider Type I, Ib/c, IIb or SNe without classification because they cannot be categorised as any of our well-defined thermonuclear and CCSNe. Only two Type Ib SNe exploded in PHANGS galaxies, so they were also excluded from the analysis. We obtained a total of 63 SNe hosted in 39 galaxies (either single or multiple SNe located in one galaxy): 16 Type Ic SNe within ACOS and 47 SNe in 23 PHANGS galaxies. Out of these 47 SNe, 12 were thermonuclear Type Ia, 30 Type II and five Type Ic. We note that the sample may be biased by observational limitations, so it can miss dust-embedded SNe and those in the outskirts of the PHANGS galaxies not covered by the ALMA data. We also do not constrain SNe that do emerge from very massive stars that have not been able to escape from their parent GMCs but are so extinguished that they do not make it into an optical sample. This can only be addressed by an infra-red-selected SN sample. We also do not constrain the properties of stars which collapse directly to black holes (BHs)55, because they do not feature in our SN samples. However, their exclusion does not affect significantly our results (see methods, subsection Statistical significance of the sample).

CO-to-Σ mol conversion

The CO(2-1) velocity-integrated intensities from moment-0 maps were converted to Σmol using the following equation (eq. 10 from ref. 16):

$${\Sigma }_{{{{{\rm{mol}}}}}}={\alpha }_{{{{{\rm{CO}}}}}}^{1-0}\,{R}_{{{{{\rm{21}}}}}}^{-1}\,{I}_{{{{{\rm{CO(2-1)}}}}}}\,\cos i,$$
(3)

where \({\alpha }_{{{{{\rm{CO}}}}}}^{1-0}\) is the CO(1-0) conversion factor, R21 is the CO(2-1)-to-CO(1-0) line ratio, i is the inclination angle of the galaxy, and ICO(2-1) is the line-integrated CO(2-1) intensity. We adopt a Galactic CO-to-H2 conversion factor of \({\alpha }_{{{{{\rm{CO}}}}}}^{1-0}\,=\,5\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2}\,{({{{{\rm{K}}}}}{{{{\rm{km}}}}}{{{{{\rm{s}}}}}}^{-1})}^{-1}\), the same as in ref. 56, and a line ratio of R21 = 0.5, from ref. 57. Inclination angles were taken from ref. 58 for the PHANGS sample, and from the Hyperleda (http://leda.univ-lyon1.fr/) galaxy database59 for ACOS galaxy hosts.

200 pc regions

An SN explosion could be located away from the centre of a cloud. In order to have a better understanding of the parent GMCs, the maximum value Σmol in a circumference within a radius of 200 pc centred on the SN position was also calculated and denoted as “200 pc region”. In the Milky Way, this radius is comparable to the maximum GMC size60. Supplementary Fig. 1 shows the Σmol eCDFs in such 200 pc regions, with a clear shift to higher densities compared with SN locations, as expected. The two-sample KS test from Supplementary Table 1 shows high probabilities that each of the ___location pairs is drawn from the same distribution. The median and 1σ values obtained via Monte Carlo simulations are shown in Supplementary Fig. 2. The fact that Type II SNe reach molecular gas densities higher than Type Ic SNe strengthens our conclusions.

Cosmological model

We use the nine-year Wilkinson Microwave Anisotropy Probe cosmological model61 with parameters H0 = 69.32 km s−1 Mpc−1, ΩΛ = 0.7134, and Ωm = 0.2865. Redshift values were used only to compute the 200 pc region for each SN and have no influence on the physical interpretation of the results.

Timing the SN progenitor lifetime with molecular gas observations

The use of molecular gas density as a tracer of stellar age is justified by the strong correlation between the age distribution and the cluster-GMC distance62,63,64. Moreover, there is a close association between the birth environment (i.e. GMC separation) and age of the cluster, measured by the equivalent width (EW) of Hα [EW(Hα)]65. Finally, the analysis of stellar associations and the CO(2-1) emission revealed that the percentage of overlap between the regions of stellar associations and GMCs is ~ 60%66.

In order to test if there is a correlation between molecular gas densities and stellar ages in the PHANGS sample, a pixel-to-pixel comparison was computed for Σmol and EW(Hα), a proxy for age. The Hα maps were obtained from the Multi Unit Spectroscopic Explorer (MUSE)67. The continuum maps were collected from the Wide Field Imager (La Silla’s 2.2m MPG/ESO telescope)68 and also available in the PHANGS–MUSE dataset.

Supplementary Fig. 3 shows the relation of Σmol and EW(Hα) for every pixel of our galaxy sample, i.e. 11 PHANGS galaxies (NGC 628, NGC 1087, NGC 1365, NGC 1385, NGC 1433, NGC 1566, NGC 1672, NGC 3627, NGC 4254, NGC 4303, and NGC 4321), with both ALMA and MUSE data, and hosting at least one SN from our sample. There is a clear correlation with pixels with lower EW(Hα) (older) having lower molecular gas density, The scatter is significant, but we take the scatter of this magnitude into account in our significance test below (and this justifies the need of a sample of the order of a few tens of SNe).

Statistical significance of the sample

To assess the statistical significance of our results with respect to the sample size, we generated 104 sets of synthetic parent GMC densities for 30 Type II SNe (as in our data) and a variable number of very massive stars in order to test if we can distinguish them. We have done it in three ways, first starting from the measured gas densities of Type II SNe (method 1), second starting from the measured gas density distribution in PHANGS galaxies (method 2), and last from lifetimes of binary systems from a numerical model (method 3).

For the former case, in order to have a realistic distribution of GMC densities we need to remove the outliers of Σmol,II data because their high values do not correspond to densities of single GMCs (as we intend to probe), but the accumulation of GMCs along the line-of-sight towards to galaxy centres, where indeed, all identified outliers are located. We obtained the first, second, and third quartiles of Σmol,II (Q1, Q2, and Q3, respectively) and considered outliers as values lower than Q1 − 1.5 IQR or higher than Q3 + 1.5 IQR, where IQR = Q3 − Q1 is the interquartile range. After removing outliers (292, 486, 515, 1442, and 5599 Mpc−2, higher than Q3 + 1.5 IQR = 157 M pc−2), we found an analytical function which best reproduces the distribution of Type II SN Σmol locations by fitting around ~80 different distributions69. The best function was a generalised normal continuous random distribution \(f(x,\,\beta )=\frac{\beta }{2\Gamma (1/\beta )}{e}^{-| x{| }^{\beta }}\), where x is a real number, β > 0 is the shape parameter, and Γ is a gamma function. The fitted parameters were β = 0.51, centred at 11.7 with a scale of 3.52. From this distribution, we constructed two different synthetic distributions corresponding to Type II SNe and very massive stars to assess our ability to distinguish them. For Type II SNe we randomly drew from the function fitted above. For the very massive stars, we made use of Eq. (2) to derive their median \({\Sigma }_{{{{{\rm{mol}}}}},{{{{\rm{massive}}}}}}={\Sigma }_{{{{{\rm{mol}}}}},{{{{\rm{II}}}}}}\,{e}^{({t}_{{{{{\rm{II}}}}}}-{t}_{{{{{\rm{massive}}}}}})}/{\tau }_{{{{{\rm{GMC}}}}}}=4\,{\Sigma }_{{{{{\rm{mol}}}}},{{{{\rm{II}}}}}}\) and drew from a similar function scaled by this factor. In this calculation we assumed the initial mass for Type II SN progenitor of Minit,II = 11 M, corresponding to a lifetime of tII = 25 Myr and an initial mass of Minit,massive = 30 M, corresponding to a lifetime of tmassive = 3 Myr. Finally, we assumed τGMC = 16 Myr.

In the second method, for each SN we drew a random progenitor age from a normal distribution of 25 ± 5 Myr and 3 ± 1 Myr for Type II SNe and very massive stars, respectively, and the lifetime of the GMC of τGMC = 16 ± 5 Myr. We also drew an initial GMC gas density from a distribution with a mean value 0.5 dex higher than the observed distribution70 (because the initial densities were higher for all of the observed clouds) and the same width, so that \(\log ({\Sigma }_{{{{{\rm{0}}}}}}\,/\,{{{{{\rm{M}}}}}}_{\odot }\,{{{{{\rm{pc}}}}}}^{-2})=2.0\pm 0.5\). The value of this parameter has no influence on the results, as this is only a normalisation and was chosen so that the median of the simulated distribution for Type II SNe is consistent with the observed value. Then we evolved the clouds as an exponential decay to calculate the surface densities at the time of the SN explosions (Eq. (1)).

For the third test, in order to take into account the effect of binarity in a simplified way, we drew samples of Type II SN progenitors and massive stars from the ranges of 6 M – Mthresh and Mthresh – 100 M, respectively, for Mthresh = 15, 20, 25, and 30 M, weighting with the Kroupa IMF. We randomly assigned an age according to the age probability distribution of SN progenitors for a given initial mass according to the models of ref. 21. Prior mass exchange of the progenitor with its companion leads in general to a longer lifetime. The two mass ranges represent tentative progenitors of Type II SNe and stripped-envelope SNe, although various binary scenarios violate this threshold. In a way, this test takes into account the change of lifetimes due to binarity, without accounting for a possible change in the SN type due to it.

For each method and for each simulated pair of sets (Type II SNe and very massive stars), we performed the KS test in order to check if we could reject the incorrect-by-design null hypothesis that they are drawn from the same distribution. Supplementary Fig. 4 shows the percentages of p-values below 0.05 (to reject the null hypothesis) and 0.37 (measured value from Table 1) for 104 Monte Carlo simulations from a KS two-sample test between the distributions of the 30 random values of Type II SNe and the massive stars constructed above, as a function of sample size for such massive stars. With the sample size of 21, as in our sample of Type Ic SNe, in these simulations, in  ~96% of the cases we obtained the p-value lower than 0.05 (and in 99.9% of cases lower than the measured value of 0.37). This means that we have statistical significance to correctly reject the null hypothesis and if Type Ic SNe were very massive stars, then we would obtain a lower p-value than we measured for virtually all the cases, so our data have enough statistical significance to rule out the very massive star hypothesis.

We also tested how the data can constrain a mixed Type Ic SN population, by analysing the fraction of the simulations with higher p-values than measured when we replaced some of the massive stars by lower mass progenitors in the same range as we assumed for Type II SNe. The 1σ range (68% of the simulated samples having a p-value higher than measured) of the accepted fraction of low-mass stars in the Type Ic SN population is 66–100%. Hence, only a third of the Type Ic SNe could be very massive stars, so that we could still measure the high p-value.

Moreover, in method 2, instead of drawing ages from normal distributions, we also drew masses according to the Kroupa IMFs71 and calculated their lifetimes according to the relationship of ref. 21. In this case, for all the values of Mthresh listed above, the number of the simulated samples having p-values lower than the measured value decreased from 99.9% to 97–98%. Finally, none of these calculations was significantly affected by the exclusion of the mass ranges for which no SNe are expected, due to a direct collapse into BHs, i.e. within the ranges 22–25 and 27–60 M55. If this is taken into account the significance increases by 1–2% due to making the difference between the Type II SNe and massive stars more pronounced.

Lifetime–initial mass relation

We converted the ZAMS masses to lifetimes using the lifetime–initial mass relation for single stars from ref. 21 (see their Fig. 1).