EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 703 (November 2025) A&A, 703 (2025) A215 Full HTML
Open Access
Issue
A&A
Volume 703, November 2025
Article Number A215
Number of page(s) 12
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202555490
Published online 18 November 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Very massive stars (VMSs), with an initial mass >100 M⊙, represent a challenge to our understanding of stellar evolution due to their extreme physical properties. Several studies have explored their evolutionary paths (Belkus et al. 2007; Yungelson et al. 2008; Gräfener et al. 2011; Brott et al. 2011; Georgy et al. 2013; Yusof et al. 2013; Hirschi 2015; Limongi & Chieffi 2018) and their observational properties in the local Universe (Crowther et al. 2010; Bestenlehner et al. 2014).

As VMSs evolve, they lose mass via stellar winds, including nitrogen and other products of hot hydrogen burning. Hence, they might be one of the main sources of the N enhancement observed in some compact high-redshift galaxies (Vink 2023; Cameron et al. 2023; Charbonnel et al. 2023; Gieles et al. 2025). Several works show that, if they directly collapse onto a black hole, VMSs can form intermediate-mass black holes (black holes with masses >102 M⊙; Ebisuzaki et al. 2001; Heger & Woosley 2002; Portegies Zwart & McMillan 2002; Portegies Zwart et al. 2004, 2005; Giersz et al. 2015; Mapelli 2016; Spera & Mapelli 2017; Costa et al. 2023).

If VMSs do not undergo direct collapse at the end of their life but instead explode as supernovae, they release enormous amounts of material into the interstellar medium (Ohkubo et al. 2006; Yokoyama & Tominaga 2012; Kojima et al. 2021; Goswami et al. 2022; Vink 2023). Specifically, one of the possible final fates of VMSs is a pair-instability supernova (PISN; Heger & Woosley 2002; Woosley 2017; Spera & Mapelli 2017). PISNe occur when the efficient formation of electron-positron pairs softens the equation of state of the carbon-oxygen core (Barkat et al. 1967; Fraley 1968; Bond et al. 1984; Heger & Woosley 2002). The resulting contraction of the core triggers explosive burning of oxygen and silicon. As a consequence, the star may undergo a single, powerful explosion that leaves no compact remnant (PISN), or it may experience multiple pulsations before ultimately collapsing into a black hole; i.e., pulsational pair-instability supernova (PPISN; Woosley et al. 2002; Chen et al. 2014; Yoshida et al. 2016; Woosley 2017). Since the onset of the instability depends on the central density and temperature, the final mass of the He core (MHe,f) is a good proxy to infer whether a star enters the pair-instability regime or not: according to Woosley (2017), if 32 M⊙ ≤ MHe,f < 64 M⊙, it undergoes a PPISN; if 64 M⊙ ≤ MHe,f < 135 M⊙, it explodes in a PISN; otherwise, if MHe,f > 135 M⊙, it collapses directly onto an intermediate-mass black hole.

Pair-instability supernovae (PISNe) and pulsational pair-instability supernovae (PPISNe) may have a deep impact on the mass function of black holes, as they are expected to open a gap between ≈60 and ≈120 M⊙ (Woosley 2017; Spera & Mapelli 2017). Our knowledge of these boundaries is, however, affected by multiple sources of uncertainties, concerning nuclear reaction rates, stellar rotation, convection, and the energetics of failed supernovae (e.g., Leung et al. 2019; Farmer et al. 2019, 2020; Mapelli et al. 2020; Marchant & Moriya 2020; Renzo et al. 2020; Song et al. 2020; Tanikawa et al. 2021; Farrell et al. 2021; Vink et al. 2021; Woosley & Heger 2021; Farag et al. 2022; Hendriks et al. 2023; Winch et al. 2025).

As for observations, we still lack an uncontroversial detection of a (P)PISN, even if several studies reported possible candidates (Woosley et al. 2007; Gal-Yam et al. 2009; Quimby et al. 2011; Cooke et al. 2012; Kozyreva & Blinnikov 2015; Kozyreva et al. 2018; Mazzali et al. 2019). The supernova SN2018ibb is possibly the strongest PISN candidate reported to date (Schulze et al. 2024). The same authors estimated that the rate density of PISNe ℛPISN should lie between 0.009 yr−1 Gpc−3 and 0.7 yr−1 Gpc−3 at 2σ confidence, based on the spectroscopically complete ZTF bright transient survey of superluminous supernovae (Fremling et al. 2020; Perley et al. 2020).

Theoretical models estimate the PISN rate to be between 10−5 and 10−2 times the core-collapse supernova rate (Langer et al. 2007; Briel et al. 2022), but such estimates are affected by uncertainties about the initial mass function and metallicity of the progenitors (Langer et al. 2007), the cosmic star formation history (Briel et al. 2022; Gabrielli et al. 2024), models of stellar evolution (Takahashi 2018; du Buisson et al. 2020; Farmer et al. 2020; Tanikawa et al. 2023; Gabrielli et al. 2024), and whether the progenitor evolves as an isolated star or in a binary system (Briel et al. 2022; Tanikawa et al. 2023).

Briel et al. (2022, hereafter B22) investigated how star-forming environments influence electromagnetic and gravitational-wave transients by employing empirical models (Madau & Dickinson 2014) and cosmological simulations (EAGLE, Schaye et al. 2015; Crain et al. 2015; ILLUSTRISTNG, Springel et al. 2005; Nelson et al. 2018; Pillepich et al. 2018; Naiman et al. 2018; Marinacci et al. 2018; MILLENNIUM, Springel et al. 2005). Among their models, the only one consistent with the upper limit found by Schulze et al. (2024) is based on the EAGLE simulation (Schaye et al. 2015), combined with the B-PASS population synthesis models (Eldridge et al. 2017). In a follow-up study, Briel et al. (2024) focused on the ILLUSTRISTNG simulation and obtain an expected PISN rate density of ℛPISN = 2–29 Gpc−3 yr−1 at redshift ɀ = 0 by quantifying the uncertainties from their binary evolutionary models.

Tanikawa et al. (2023, hereafter T23) explored the influence of the 12C(α,γ)16O reaction rate on the PISN rate in models of binary populations. By adopting the standard reaction rate and assuming a maximum zero-age-main-sequence (ZAMS) mass of 300 M⊙, they found ℛPISN ≥ 1 yr−1Gpc−3. In contrast, using a 3σ-lower reaction rate, they found ℛPISN ⪡ 1 yr−1Gpc−3. The critical role of the 12C(α,γ)16O reaction rate has been largely debated in the context of the pair-instability mass gap of black holes, with lower rates leading to higher black-hole masses below the gap (e.g., Farmer et al. 2019, 2020; Costa et al. 2021 ; Woosley & Heger 2021).

Gabrielli et al. (2024, hereafter G24) explored various criteria for the onset of pair instability, including the impact of the metallicity threshold for PISNe, the maximum ZAMS mass, and the dispersion in the metallicity distribution of the host galaxies. They found that all these quantities affect the PISN rate, resulting in a seven orders of magnitude variation up to ℛPISN ~ 2000 yr−1Gpc−3. According to their analysis (Fig. 8 from Gabrielli et al. 2024), only models with a metallicity threshold for the onset of pair instability Z ≤ 1.5 × 10−3, a maximum ZAMS star mass MZAMS ≤ 150 M⊙, and a narrow metallicity spread yield rates ℛPISN < 1 Gpc−3 yr−1. This result confirms that the PISN rate is very sensitive to the assumed metallicity threshold for the onset of PISNe. Moreover, assuming a narrow metallicity dispersion suppresses pockets of low-metallicity star formation in the nearby Universe (Boco et al. 2019; Santoliquido et al. 2022).

In summary, most theoretical estimates of the PISN rate density exceed the upper limit reported by Schulze et al. (2024). The only models that yield PISN rates below this upper limit assume low values for the 12C(α,γ)16O reaction rate, enforce a low metallicity threshold for the onset of PISNe, suppress the low-metallicity tails of star formation in the nearby Universe, and/or put a cap on the maximum ZAMS mass. The upper limit to the maximum ZAMS mass is particularly critical, because we observed stars with MZAMS > 150 M⊙ even in the local Universe (Crowther et al. 2010; Bestenlehner et al. 2014, 2020).

Here, we explore an alternative scenario inspired by recent models of optically thick winds for VMSs (Sabhahit et al. 2023, hereafter S23). This scenario does not require fine-tuning of the 12C(α,γ)16O reaction rate, works for a realistic evolution of the metallicity spread in galaxies, and does not require a cap on the maximum star mass, while it naturally leads to a low metallicity threshold Zth for the activation of PISNe (Zth ~ Z⊙/10, S23).

Although the stellar winds of massive and very massive stars remain an active topic of research (Yang et al. 2023; Cheng et al. 2024; Pauli et al. 2025), recent models of mass loss by stellar winds (Vink et al. 2011; Sabhahit et al. 2022, 2023) indicate that VMSs already develop optically thick winds on the main sequence, becoming WNh stars (Bestenlehner et al. 2020). The transition between optically thin winds of O-type stars and optically thick ones is a function of the electron-scattering Eddington ratio, ΓEdd. At the metallicity of the Large Magellanic Cloud and even the Small Magellanic Cloud, this model strips the H-rich envelope of VMSs completely, without the need for mass loss through binary interaction. Here, we investigate the impact of such new stellar wind models (S23) on the expected PISN and PPISN rate density. We show that a larger mass-loss rate for VMSs at intermediate metallicity (down to Z ~ 0.002, S23) seems to be the key to predicting PISN rates in agreement with the observational constraints.

This paper is organized as follows. We briefly describe the physics and the prescription adopted by the (P)PISN rate density in Sect. 2. In Sect. 3, we show the rate densities obtained with the two different stellar models. We discuss the impact of an enhanced mass-loss rate in Sect. 4. Finally, our results are summarized in Sect. 5.

2 Method

We used the stellar evolution code MESA (version r12115; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) to model massive stars and VMSs and produce stellar tables for the population synthesis code SEVN (Spera & Mapelli 2017; Spera et al. 2019; Mapelli et al. 2020; Iorio et al. 2023). Then, we combine our population synthesis models with a metal-dependent star formation rate density history by means of the semi-analytic code COSMOℛATE (Santoliquido et al. 2021, 2022). Below, we describe our methodology in detail.

2.1 Stellar models with MESA

We adopted the MESA mass-loss scheme and parameter setup introduced by S231, with some changes to enforce convergence until the formation of the CO core. In the following, we give an overview of our main parameters and initial conditions.

During the main sequence, we used the mass-loss framework described by Sabhahit et al. (2022) and extended to metal-poor stars by S23. This model describes the transition from an optically thin to an optically thick wind in VMSs (Vink et al. 2011; Bestenlehner et al. 2014), exploiting the concept of transition mass-loss rate (Vink & Gräfener 2012). The transition happens when ΓEdd ≥ ΓEdd,tr, where ΓEdd is the Eddington ratio of a star of luminosity L*, defined as the ratio between L* and the Eddington luminosity of the star, whereas ΓEdd,tr is the value of the Eddington ratio of the same star computed at the transition point. When stellar winds become optically thick, the mass-loss rate is enhanced according to Vink et al. (2011): M˙=M˙tr(LLtr)4.77(MMtr)−3.99,$\dot M = {\dot M_{{\rm{tr}}}}{\left( {{L \over {{L_{{\rm{tr}}}}}}} \right)^{4.77}}{\left( {{M \over {{M_{{\rm{tr}}}}}}} \right)^{ - 3.99}},$(1)

where Ṁtr, Ltr, and Mtr are the mass-loss rate, the luminosity, and the mass at the transition point. To obtain ΓEdd,tr, the following condition is iteratively verified: ηVink≡M˙Vink v∞L/c=ητF, s(v∞vesc)~(1+vesc2v∞2)−1,${\eta _{{\rm{Vink}}}} \equiv {{{{\dot M}_{{\rm{Vink}}}}\,{v _\infty }} \over {L/c}} = {\eta \over {{\tau _{F,\,s}}}}\left( {{{{v _\infty }} \over {{v _{{\rm{esc}}}}}}} \right)\~{\left( {1 + {{v _{{\rm{esc}}}^2} \over {v _\infty ^2}}} \right)^{ - 1}},$(2)

where ṀVink is the mass-loss rate computed following Vink et al. (2001), v∞ the asymptotical velocity, vesc the escape velocity, η the wind efficiency of the model, and τF the flux-weighted mean optical depth. This condition allows us to evaluate Ltr at each iteration, which is necessary to compute ΓEdd,tr.

During the core He-burning phase, we only used Eq. (1) for temperatures between 4 × 103 K and 105 K, while for lower temperatures (T < 4 × 103 K), we followed de Jager et al. (1988). For T > 105 K and the hydrogen-mass fraction XC < 0.01, we adopted the Wolf–Rayet (WR) mass-loss models by Sander & Vink (2020). Recently, Boco et al. (2025) showed that this formalism, when combined with fast rotation, can strip single metal-poor O-type stars and convert them into WNh stars even at low metallicity (Z ~ 0.0002). Thus, this mechanism offers a potential explanation for the presence of WR stars in metal-poor galaxies, such as the Small Magellanic Cloud (Schootemeijer et al. 2021).

We treat convective mixing with the mixing-length theory (MLT) by Cox & Giuli (1968), with a constant mixing-length ratio of αMLT = 1.5. The convective boundaries were set by applying the Ledoux criterion (Ledoux 1947). We used the semi-convective diffusion parameter αSC = 1. We chose the exponential overshooting by Herwig (2000) to describe the over-shooting above the convective regions with efficiency parameter ƒOV = 0.03. We only considered non-rotating models. Finally, we only turned on the MLT++ option (Paxton et al. 2013) during the core He-burning phase to ensure convergence. MLT++ is an option that allows us to reduce superadiabaticity (∇T > ∇ad, where ∇T is the temperature gradient and ∇ad the adiabatic temperature gradient) in radiation-dominated convective regions with a posteriori calculations. We refer the reader to S23 for a more in-depth description of the parameters and the model.

We made exactly the same assumptions as S23 with only two exceptions: the definition of core boundaries and the nuclear reaction network. As for the core boundaries, we define the radius of the He core as the radius inside which the mass fraction of hydrogen drops below ƒX = 10−4. Similarly, we define the radius of the helium (carbon) core as the radius inside which the mass fraction of 4He (12C) drops below ƒY = 10−4 ( ƒC = 10−4). The default of MESA instead sets the boundary of the H, He, and C core where the mass fraction of H, 4He, and 12C goes below 0.1, respectively. Our definition yields a more precise description of the core and is the same as that assumed by the MIST libraries (Paxton et al. 2011, 2013, 2015; Dotter 2016; Choi et al. 2016).

We adopted the default nuclear reaction network in MESA with one exception: after the main sequence, we enabled the co_burn net (parameter auto_extend_net = true), which includes 28Si, in order to have a more complete and extended network for C- and O-burning and α-chains with respect to S23.

We calculated stellar models with initial masses ranging from 50 to 500 M⊙ by intervals of 25 M⊙. The initial metallicities are Z = 0.008, 0.006, 0.004, 0.002, 0.001, 0.0008, 0.0006, 0.0004, 0.0002, and 0.0001. Thus, our grid consists of 190 models (19 initial masses and 10 initial metallicities). The initial He-mass fraction, Y, was computed as Y = Yprim + (∆Y/∆Z), Z, where Yprim is the primordial He abundance and ∆Y/∆Z describes the helium-to-metal ratio. We adopted the default values from MESA (Pols et al. 1998) with Yprim = 0.24 and ∆Y/∆Z = 2. This formulation spans the range from a nearly primordial composition (X = 0.76, Y = 0.24, Z = 0) to a quasisolar one (X = 0.70, Y = 0.28, Z = 0.02). The hydrogen mass fraction is given by X = 1 − Y − Z. We evolved all the stellar models until the luminosity due to carbon burning, LC, is less than the 20% of the total luminosity LTOT (LC/LTOT < 0.2) and the central He abundance is YC < 10−8.

Figure 1 shows the mass-loss rate of our models for four selected metallicities, until the end of core He burning. During the main sequence, the initial mass required to enter the optically thick wind regime grows with decreasing metallicity: at Z = 0.0008, models with MZAMS ≤ 125 M⊙ do not experience optically thick winds, while at Z = 0.008 even the less massive models experience enhanced winds.

When stars that begin their life with ΓEdd < ΓEdd,tr enter the optically thick wind regime, they experience a boost in mass loss. Instead, the mass-loss rate of stars that begin their life with ΓEdd > ΓEdd,tr reaches a peak of ~10−3.5 M⊙ yr−1 during the main sequence and then starts to decrease. We discuss other properties of our MESA stellar models in Appendix A.

2.2 Population synthesis with SEVN

We used the population-synthesis code SEVN (Spera & Mapelli 2017; Spera et al. 2019; Mapelli et al. 2020; Iorio et al. 2023)2 to interpolate among our MESA models and estimate the (P)PISN rate density. SEVN evolves populations of single and binary stars by interpolating pre-computed stellar tracks and incorporating semi-analytic models for mass transfer, tides, supernova explosions, and gravitational-wave decay. We refer to Iorio et al. (2023) for a complete description of the code.

Here, we used the outputs from our new MESA models to obtain new evolutionary tables for stars with masses of M★ ≥ 50 M⊙. We considered the same initial metallicities as in the MESA models. We simulated 105 stars for each metallicity. The initial masses are distributed following a Kroupa initial mass function (IMF; Kroupa 2001) with a slope of α = 2.3 between 50 M⊙ and 500 M⊙. To model PISNe and PPISNe, we used the fitting formulas described by Mapelli et al. (2020) based on the models of Woosley (2017). Namely, a star undergoes a PPISN if the final mass of the He core is around 32 M⊙ ≤ MHe,ſ < 64 M⊙. A PISN takes place if 64 M⊙ ≤ MHe,f ≤ 135 M⊙. For higher masses of the He core, the star directly collapses onto an intermediate-mass black hole.

As a comparison, we also considered the SEVN tables produced with the stellar evolution code PARSEC (Bressan et al. 2012; Chen et al. 2015; Costa et al. 2019a,b; Nguyen et al. 2022; Costa et al. 2025) and described by Iorio et al. (2023). The main difference between the PARSEC and MESA models is that the former adopts a different stellar wind model from the one described by S23. The PARSEC wind model introduced by Chen et al. (2015, hereafter C15) is based on the fitting formulas by Vink et al. (2000) and Vink et al. (2001) for the winds of O-type stars, with additional corrections for electron scattering, and a dependence of the mass loss on the Eddington ratio (Gräfener & Hamann 2008; Vink et al. 2011). We provide more details about the PARSEC models in Appendix B.

thumbnail Fig. 1

Mass-loss rate of our MESA models as function of time, from the ZAMS until the end of core He burning. The models adopt the winds by Sabhahit et al. (2023). From top to bottom: Z = 0.008, 0.004, 0.001, and 0.0008. The solid lines correspond to the evolutionary phase in which the stellar models are in the optically thick regime (ΓEdd > ΓEdd,tr), while the dashed lines mark the optically thin regime.

2.3 (P)PISN rate-density evolution across cosmic time

We modeled the evolution of the (P)PISN rate density as a function of redshift with the code COSMOℛATE3 (Santoliquido et al. 2020, 2021) based on the observed metal-dependent cosmic star formation rate density (Madau & Fragos 2017). The default version of COSMOℛATE estimates the merger rate density of compact binaries as a function of the redshift, but it can be adapted to trace the rate density of other transients of stellar origin, including (P)PISNe. Specifically, we calculated the rate of (P)PISNe as follows: R(ɀ)=∫ɀmaxɀ[ ∫ZminZmaxS(ɀ′, Z)F(ɀ′, ɀ, Z)dZ ]dt(ɀ′)dɀ′dɀ′,$R\left( z \right) = \int_{{z_{\max }}}^z {\left[ {\int_{{Z_{\min }}}^{{Z_{\max }}} {S\left( {z',\,Z} \right)F\left( {z',\,z,\,Z} \right){\rm{d}}Z} } \right]{{{\rm{d}}t\left( {z'} \right)} \over {{\rm{d}}z'}}{\rm{d}}z',} $(3)

where 𝒮(ɀ′, Z) = ψ(ɀ′) p(Z|ɀ′).

Here, ψ(ɀ′) is the cosmic star formation rate density at red-shift ɀ′, and p(Z|ɀ′) is the distribution of metallicity Z at fixed formation redshift ɀ′. In Equation (3), dt(ɀ′)/dɀ′=H0−1(1+ɀ′)−1[ (1+ɀ′)3ΩM+ΩΛ ]−1/2${\rm{d}}t\left( {z'} \right)/{\rm{d}}z' = H_0^{ - 1}{\left( {1 + z'} \right)^{ - 1}}{\left[ {{{\left( {1 + z'} \right)}^3}{\Omega _{\rm{M}}} + {\Omega _\Lambda }} \right]^{ - 1/2}}$, where H0 is the Hubble parameter and ΩM and ΩΛ are the matter and energy density, respectively. We adopted the values from Planck Collaboration (2020).

The term ℱ(ɀ′,ɀ, Z) in Equation (3) is given by F(ɀ′, ɀ, Z)=1M*(Z)dN(ɀ′, ɀ, Z)dt( ɀ),$F\left( {z',\,z,\,Z} \right) = {1 \over {{M_*}\left( Z \right)}}{{{\rm{d}}N\left( {z',\,z,\,Z} \right)} \over {{\rm{d}}t\left( {\,z} \right)}},$(4)

where d𝒩(ɀ′, ɀ, Z)/dt(ɀ) is the rate of (P)PISNe from stars with a metallicity of Z that formed at redshift ɀ′ and undergo a (P)PISN at redshift ɀ, and M*(ɀ′, Z) is the initial total stellar mass that forms at redshift ɀ′ with a metallicity of Z. We obtain d𝒩(ɀ', ɀ, Z)/dt(ɀ) directly from our simulations. We note that 𝒩(ɀ', ɀ, Z) also depends on the stellar lifetime (and, hence, on the stellar mass) and on the assumed initial mass function. While we do not explicitly state these dependencies for notation simplicity, our catalogs fully account for them.

We modeled the star formation rate density through the fitting formula by Harikane et al. (2022): ψ(ɀ)=[ 61.7(1−ɀ)−3.13+100.22(1+ɀ)+2.40.5(1+ɀ)−3 ]−1M⊙Mpc3 yr.$\psi \left( z \right) = {\left[ {61.7{{\left( {1 - z} \right)}^{ - 3.13}} + {{10}^{0.22\left( {1 + z} \right)}} + {{2.4}^{0.5\left( {1 + z} \right) - 3}}} \right]^{ - 1}}{{{{\rm{M}}_ \odot }} \over {{\rm{Mp}}{{\rm{c}}^3}\,{\rm{yr}}}}.$(5)

They obtained this formula through the analysis of new measurements of the rest-frame UV luminosity at ɀ ~ 4–7 based on wide and deep optical images obtained in the Hyper Supreme Cam (HSC) Subaru Strategic Program (SSP) survey (Aihara et al. 2018) and the CFHT Large Area U-band Deep Survey (CLAUDS; Sawicki et al. 2019). This formula is a better match to recent James Webb Space Telescope (JWST) data (Donnan et al. 2023; Harikane et al. 2023, 2024) than the traditionally adopted fit by Madau & Dickinson (2014). The star formation rate density is then multiplied by a correction factor of ƒKroupa = 0.66 to convert the SFR from a Salpeter IMF to a Kroupa IMF (Madau & Dickinson 2014).

The metallicity-evolution term p(Z|z′) is the probability density function that a star that formed at redshift ɀ′ has a metallicity of Z, modeled as a log-normal distribution function (Chruslinska et al. 2019; Santoliquido et al. 2023): p(Z|ɀ′)=12 π σZ2exp{ −[ log(Z/Z⊙)−〈 log(Z/Z⊙) 〉 ]22 σZ2 }.$p\left( {Z|z'} \right) = {1 \over {\sqrt {2\,\pi \,\sigma _{\rm{Z}}^2} }}\exp \left\{ { - {{{{\left[ {\log \left( {Z/{Z_ \odot }} \right) - \left\langle {\log \left( {Z/{Z_ \odot }} \right)} \right\rangle } \right]}^2}} \over {2\,\sigma _{\rm{Z}}^2}}} \right\}.$(6)

We derived 〈 log(Z/Z⊙)〉 from the fitting formula by Madau & Fragos (2017): log〈 Z/Z⊙ 〉=0.153−0.074 ɀ1.34,$\log \left\langle {Z/{Z_ \odot }} \right\rangle = 0.153 - 0.074\,{z^{1.34}},$(7)

knowing that 〈 logZ/Z⊙ 〉=log〈 Z/Z⊙ 〉−ln(10)σZ2/2$\langle \log Z/{Z_ \odot }\rangle = \log \langle Z/{Z_ \odot }\rangle - \ln \left( {10} \right)\sigma _Z^2/2$. We assumed σZ = 0.4, based on the results by Sgalletta et al. (2025).

thumbnail Fig. 2

Contour plot showing He core masses at the end of He burning for our MESA tracks adopting S23 winds as a function of the metallicity and ZAMS mass. The white lines highlight the levels at 32 M⊙. 64 M⊙, and 135 M⊙, which set the boundaries to have PPISNe, PISNe, and direct collapse via photodisintegration.

3 Results

3.1 The metallicity threshold for (P)PISNe

Figure 2 shows the final He core mass for our MESA models with S23 winds as a function of metallicity and ZAMS mass. We highlight the levels corresponding to MHe,f = 32, 64, and 135 M⊙, because these correspond to the threshold for PPISN, PISN, and direct collapse according to Woosley (2017). In the S23 model, the maximum He core mass of stars with MZAMS up to 500 M⊙ is below 65 M⊙ for metallicity Z ≳ 0.003 and below 135 M⊙ for Z≳ 0.0015.

Assuming that the He core mass is a good proxy for the onset of pair instability and adopting the models by Woosley (2017), the results presented in Figs. 2 and 3 allow us to estimate the threshold metallicity for the onset of (P)PISNe for both the S23 and C15 models. With our new models (S23), the maximum metallicity for PISNe and PPISNe are Z ~ 2 × 10−3 and Z ~ 4.5 × 10−3, respectively. For comparison, the model adopting the PARSEC stellar tracks (C15) yields Z ~ 1.7 × 10−2 and Z ≳ 2 × 10−2 as thresholds for PISNe and PPISNe, respectively (Fig. 3). Such metallicities are extremely high, above the solar value (Z ~ 0.015; e.g., Caffau et al. 2011).

In the S23 model (Fig. 2), even stars whose initial masses are higher than 250 M⊙ explode as PISNe if their metallicity lies between 0.001 and 0.002. At lower metallicities, stars with ZAMS masses >250 M⊙ collapse directly into black holes, because photodisintegration prevents the explosion, while stars with ZAMS mass between 110 and 230 M⊙ undergo PISNe. Stars with 64 M⊙ ≤ MZAMS ≤ 110 M⊙ undergo PPISN for Z ≤ 0.004 and this mass range extends up to MZAMS = 500 M⊙ for 0.002 ≤ Z ≤ 0.004.

Figure 4 shows the evolution of the He core mass at the end of core He burning for some selected metallicities, comparing S23 and C15. From this Figure, it is apparent how the onset of optically thick winds during the main sequence progressively “erodes” the core of the most massive stars for Z ≥ 0.0004 in our MESA models. In contrast, the PARSEC models grow large cores even at relatively high metallicity (Z ≥ 0.0004) because of the different wind model.

However, at low metallicities (Z ≲ 0.0004), the PARSEC models tend to develop smaller cores than the MESA models for ZAMS masses MZAMS ≥ 80 M⊙. This difference arises from the different envelope undershooting and convection models adopted by MESA and PARSEC (see Appendices A and B), rather than by the wind model. Specifically, the PARSEC tracks develop effective dredge-up episodes that exchange mass between the envelope and the He core, reducing the overall He core mass at the end of core He burning (Costa et al. 2021).

thumbnail Fig. 3

Same as Fig. 2, but for models C15.

3.2 (P)PISN rate density

Figure 5 shows the evolution of the PISN and PPISN rate density across cosmic time. We compare our new MESA models adopting the mass-loss rate by S23 to the PARSEC models with the mass-loss rate presented by C15. We also compare the models with the observational constraints reported by Schulze et al. (2024).

At low redshift, the (P)PISN rate density of the S23 model is lower by two orders of magnitude compared to the C15 model (see Table 1). The PISN rate of the S23 model is below the observational upper limit set by Schulze et al. (2024), whereas the C15 model yields a much higher rate than the upper limit.

The main reason for this difference is that the metallicity threshold to activate the PISNe is lower in the S23 model (Z ≤ 0.002) than in the C15 model (Z ≤ 0.017), as shown in Figs. 2 and 3. According to the C15 model, we expect PISNe to take place even at the metallicity of the Small (Z ≈ 0.003) and Large Magellanic Clouds (Z ≈ 0.008), whereas the S23 model predicts that PISNe only take place at metallicities lower than that of the Small Magellanic Cloud.

At high redshift, the main contribution to the PISN rate comes from very metal-poor stars (Z < 0.002), for which the mass-loss rates adopted in the S23 and C15 models behave almost in the same way (Fig. 6). However, the (P)PISN rates obtained with models S23 and C15 show some differences even at high redshift. Specifically, at redshift ɀ ≥ 8, the rate of the S23 model becomes slightly higher than the rate of the C15 models. This happens because of other differences between the MESA and PARSEC stellar tracks that we already discussed in Section 3.1. Specifically, Fig. 4 shows that the final He core mass of the PARSEC tracks (C15) at metallicity 0.0001 ≤ Z ≤ 0.0006 tends to be lower than the one of the MESA tracks (S23), reducing the overall PISN rate in the C15 models.

thumbnail Fig. 4

Final He core mass as a function of the ZAMS mass. The orange lines are our new MESA tracks (models S23), while the dashed blue lines correspond to results from PARSEC tracks (models C15). Dark (light) gray shaded areas: regime for PISNe (PPISNe). From top left to bottom right, we show metallicities Z = 0.0001, 0.0004, 0.0008, 0.001, 0.004, and 0.008.
thumbnail Fig. 5

PISN (orange lines) and PPISN (blue lines) rate densities as function of redshiſt. The solid and dashed lines are obtained using the models by S23 and C15, respectively. The green shaded area represents the observational constraints set by Schulze et al. (2024). The black arrow highlights the upper limit from Schulze et al. (2024).
Table 1

PISN and PPISN rate densities at redshiſt ɀ = 0.

thumbnail Fig. 6

PISN (upper panels) and PPISN (lower panels) rate density evolution as function of redshift. The left-hand and right-hand plots show the results obtained with the new wind model by S23 and the one by C15, respectively. The thick black line represents the total rate density, whereas the colored lines show the contribution of individual metallicities from Z = 1 × 10−4 (blue dashed line) to Z = 2 × 10−2 (solid yellow line).

3.3 The impact of metallicity evolution

Figure 6 shows the contribution of various metallicities to the (P)PISN rate for our two models. Below redshift ɀ = 5, the two models behave in a completely different way: if we assume the winds by S23, the main contribution to PISNe (PPISNe) comes from metallicity Z = 0.001–0.002 (Z = 0.002), whereas if we adopt the winds by C15 the dominant metallicity is Z = 0.01 for both PISNe and PPISNe. At ɀ > 5, the two models evolve in a more similar way, and lower metallicities also become important.

Figure 6 assumes our fiducial metallicity spread σZ= 0.4 (see Eq. (6)). In Figure 7, we show the impact of the parameter σZ on our results. Specifically, we change the metallicity spread σZ between 0.1 and 0.6. Our fiducial value of σZ = 0.4 is compatible with observational calibrations (see, e.g., Chruslinska et al. 2019; Santoliquido et al. 2021). However, in Figure 7 we treat σZ as a parameter, given the significant uncertainties and observational biases (see, e.g., Sgalletta et al. 2025, for a discussion).

In the results based on the C15 models, the rates are mainly influenced by medium-to-high metallicities with a peak at ɀ ~ 2 determined by metallicity Z = 0.01 (see also Fig. 6). As a consequence, changing σZ has little effect, introducing only minor fluctuations at redshift ɀ = 7 for (P)PISN rate density, and leading to PPISN rate densities at ɀ = 0 that range from 2 to 8.5 yr−1 Gpc−3. This occurs because of the high metallicity threshold for PISN in C15 models. In this case, the rate density is dominated by the peak of the SFRD, with the spread of the metallicity playing only a minor role (Gabrielli et al. 2024). As a result, the rate is almost insensitive to σZ.

In contrast, when enhanced winds are considered, there is no contribution from metallicities Z > 2 × 10−3 for PISNe (Z > 4 × 10−3 for PPISNe). This leads to a strong dependence on metallicity evolution, allowing the PISN rate to be always below the upper limit estimated by Schulze et al. (2024), even when σZ ~ 0.5. The PISN rate density at redshift ɀ = 0 even falls below the lower limit set by Schulze et al. (2024) when assuming σZ < 0.3. In this case, a different value of σZ shifts also the expected peak from ɀ ~ 7 (assuming σZ = 0.1) to ɀ ~ 2 (assuming σZ = 0.6). These results depend on our simplifying assumption that the distribution of metallicities at a given red-shift follows a log-normal distribution around the mean value. We will study more detailed distributions in a follow-up study. Model C15 is almost unaffected by the spread (σZ) because the threshold metallicity for (P)PISNe is so high that there is always a subpopulation of VMSs that end their life with a (P)PISN, even at low redshift. The rate then becomes more sensitive to the IMF (maximum stellar mass and slope of the mass function) than to the metallicity spread.

Gabrielli et al. (2024) found a similar behavior. They modeled the metallicity as a log-normal distribution centered on the fundamental metallicity relation (FMR; Curti et al. 2020), with a dispersion of σ˜Z${\tilde \sigma _{\rm{Z}}}$. When adopting the metallicity threshold Zth = 0.017, varying σ˜Z${\tilde \sigma _{\rm{Z}}}$ between 0.15 and 0.70 leads to nearly identical PISN rates at redshift ɀ = 0. In contrast, for a lower threshold (Zth = 0.0015), the same variation in σ˜Z${\tilde \sigma _{\rm{Z}}}$ results in a discrepancy of about five orders of magnitude (see Table 1). This indicates that the PISN rate is determined not only by σZ, but also by the interplay between σZ and the metallicity threshold Zth.

thumbnail Fig. 7

PISN (upper row) and PPISN (lower row) as function of redshift, varying the metallicity spread from σZ = 0.1 (black dashed line) to σZ = 0.6 (yellow dashed line). Left hand (right hand) panels: Models with stellar winds from S23 (from C15).

4 Discussion

The stellar-wind model proposed by S23 has a strong impact on the evolution of VMSs. The high mass-loss rate in the optically thick regime causes the most massive stars at metallicity Z ≥ 0.004 to end their life with an He core mass lower than 32 M⊙, without entering the (pulsational) pair-instability regime. In contrast, at metallicity Z = 0.001, even stars with masses between 250 and 500 M⊙ can undergo PISNe. This has important consequences on the expected mass function of black holes (Törniamenti et al. 2025). Here, we find that the enhanced mass-loss rate lowers the maximum metallicity for PISNe from Z ~ 1.7 × 10−2 to Z ~ 2 × 10−3 as also shown by S23. The absence of contributions by higher metallicities suppresses the predicted peak of the PISN rate at redshift ɀ ~ 2 and reduces the expected rate at redshift ɀ = 0 by two orders of magnitude compared to the PARSEC stellar tracks (C15). The PISN rate density of our MESA models falls within the observational constraints estimated by Schulze et al. (2024): 0.009 Gpc−3 ≤·ℛPISN ≤ 0.7 Gpc−3−yr−1.

To better comprehend the critical impact of stellar winds and metallicity on the (P)PISN rate, we compare our results with those of previous works. A summary of all PISN rate densities at ɀ = 0 discussed here is provided in Table 1.

B22 demonstrated that the variation of the star formation rate model widely affects the PISN rate. Their empirical model aligns with our results from model C15, showing a peak around ɀ ~ 2 and a similar value at ɀ = 0 of approximately ℛPISN ~ 6 yr−1 Gyr−3. In the other models they considered the PISN rate peak is at ɀ ~ 3, resulting in a lower rate at ɀ = 0, but still exceeding the more recent observational constraints.

T23 explored the impact of variations in the 12C(α,γ)16O reaction rate, predicting a peak of the PISN rate around ɀ ~ 6. This higher redshift may be attributed to the different metallicity distribution they assumed. As already discussed by G24, binary interactions are not expected to significantly impact the predicted rate. Moreover, the inclusion of enhanced stellar winds can lead to wider binaries due to the large amount of mass ejected, reducing the effects of binary evolution.

G24 also investigated how metallicity evolution influences the expected PISN rate. They assumed the metallicity distributed as a log-normal distribution around the fundamental metallicity relation (FMR; Curti et al. 2020) with a metallicity spread of σ˜Z${\tilde \sigma _{\rm{Z}}}$. Changing σ˜Z${\tilde \sigma _{\rm{Z}}}$ from 0.15 to 0.7 and the maximum metallicity to have PISN Zth from 1.5 × 10−3 to 1.7 × 10−2, they observed a change of several orders of magnitude in the PISN rate, largely exceeding the upper limits by Schulze et al. (2024) when (σ˜Z≳0.35${\tilde \sigma _{\rm{Z}}} \mathbin{\lower.3ex\hbox{$\buildrel>\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 0.35$ and Zth > 1.5 × 10−3. In our work, the metallicity spread σZ is defined as the standard deviation of a log-normal distribution around the average metallicity evolution of stars (Eq. (6)). Hence, it is not directly comparable with the parameter σ˜Z${\tilde \sigma _{\rm{Z}}}$ adopted by G24, even if both parameters quantify the dispersion of the stellar metallicity distribution (see the discussion by Sgalletta et al. 2025).

As shown in the left hand panels of Fig. 7, the PISN rates of our MESA models are also very sensitive to metallicity, but, considering σZ ≤ 0.5, the results are consistent with the observational upper limit. We do not expect that a more detailed modeling of metallicity evolution would affect our results. Indeed, since the S23 models do not include contributions from high metallicities, the high-metallicity tail predicted by the log-normal distribution of Z has no impact on the outcome.

Overall, previous studies able to reproduce the observed PISN rate (see Table 1) either adopt a value that falls in the tail of the observational distribution for some of the parameters, or use a metal-dependent SFRD derived from cosmological simulations. For instance, T23 adopted a value of the 12C(α, γ)16O reaction rate that is 3σ lower than the measured one, whereas B22 assumed the metal-dependent SFRD derived from the EAGLE cosmological simulation, which generally yields a higher median stellar metallicity than empirical models (e.g., Artale et al. 2020). Here, we show that our approach, based on a stellarwind model calibrated against observational data (S23) and data-driven prescription for the SFRD (Harikane et al. 2022), results in a PISN rate density in agreement with the recent observational constraints (Schulze et al. 2024) without assuming extreme values for nuclear reaction rates or other parameters. Moreover, due to optically thick winds, the S23 model yields a metallicity threshold (Z ~ 0.003), above which the PISN is completely avoided, regardless of the maximum stellar mass considered.

In a companion paper (Boco et al. 2025), the authors argue that the S23 model, coupled with fast rotation or a higher over-shooting parameter, can also explain the existence of companionless WR stars in the Small Magellanic Cloud (Schootemeijer et al. 2021). Additional tests and comparisons are needed in order to address whether such models are compatible with other observational benchmarks; for example, the Humphreys-Davidson limit (Humphreys & Davidson 1979).

5 Summary

Here, we studied the impact of the new stellar wind model by S23 on the rate of PISNe and PPISNe, by integrating new VMS tracks with MESA (Paxton et al. 2019) and by convolving our results with a metal-dependent star formation history (Santoliquido et al. 2021). The model by S23 describes the transition from optically thin to optically thick winds already resulting in high mass-loss rates (Ṁ ≳ 10−4 M⊙/yr) during the main sequence of VMSs for metallicities of Z ≳ 0.0004. This quenches the growth of VMSs’ cores (Fig. 4) and prevents the onset of PISNe at metallicity Z ≥ 2 × 10−3, even if we assume a maximum ZAMS mass as high as 500 M⊙ (S23).

We find that the PISN (PPISN) rate density at ɀ = 0 is ℛPISN = 0.08 Gpc−3 yr−1 (ℛPPISN = 0.8 Gpc−3 yr−1 ) for the S23 wind model, if we assume a metallicity dispersion of σZ = 0.4 (Fig. 6 and Table 1). By changing the metallicity dispersion, σZ, from 0.1 to 0.6, the local (P)PISN rate density changes by several orders of magnitude (Fig. 7), because it is extremely sensitive to the fraction of star formation that occurs above the metallicity threshold for (P)PISNe.

For σZ = 0.4–0.5, our expected local PISN rate density (ℛPISN, z=0 = 0.08–0.4 Gpc−3 yr−1) is consistent with the observational constraints (0.009 < ℛPISN/Gpc−3 yr−1 < 0.7) inferred by Schulze et al. (2024). Notably, this result does not require to cap the maximum stellar mass or to assume an unrealistically thin metallicity dispersion.

In our S23 models, the PISN rate density grows up to a few × Gpc−3 yr−1 at redshift ɀ = 2–6. The main contribution to PISNe and PPISNe comes from stars with metallicity Z = 0.001–0.002.

We compared our new MESA models with the PARSEC stellar tracks reported by Costa et al. (2025), adopting the fitting formulas by Vink et al. (2001) for O-type stars, with corrections for electron scattering introduced by C15. The PARSEC tracks lead to a higher maximum metallicity (Z ≥ 1.7 × 10−2) for PISNe if we assume the same maximum ZAMS mass. As a consequence, the PISN rate density at ɀ = 0 is much higher (~6 Gpc−3 yr−1), above the observational upper limit. Our results confirm the importance of stellar winds to describe the evolution of VMSs, the rate of (P)PISNe, and the expected mass function of black holes (Törniamenti et al. 2025).

Acknowledgements

We thank the anonymous referee for the constructive comments and suggestions that helped improve the manuscript. We whole-heartedly thank Gautham Sabhahit and collaborators for sharing their MESA inlists (https://github.com/Apophis-1/VMS_Paper1, https://github.com/Apophis-1/VMS_Paper2). Our work would not have been possible without this invaluable open-science input. We thank Guglielmo Costa, Marco Dall’Amico, Varsha Ramachandran, Andreas Sander, Raffaele Scala, Jorick Vink, Gautham Sabhahit, Daniel Pauli and Amedeo Romagnolo for useful discussions. MM, GI, and ST acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. MM and ST also acknowledge financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) STRUCTURES. ST acknowledges financial support from the Alexander von Humboldt Foundation for the Humboldt Research Fellowship. GI acknowledges financial support from the La Caixa Foundation for the La Caixa Junior Leader fellowship 2024. GI also acknowledges financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.4, – Call for tender No. 3138 of 18/12/2021 of Italian Ministry of University and Research funded by the European Union – NextGenerationEU. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG. We use the MESA software (https://docs.mesastar.org/en/latest/) version r12115; (Paxton et al. 2011, 2013, 2015, 2018, 2019). We used SEVN (https://gitlab.com/sevncodes/sevn) to generate our BBHs catalogs (Spera et al. 2019; Mapelli et al. 2020; Iorio et al. 2023). We used TRACKCRUNCHER (https://gitlab.com/sevncodes/trackcruncher) (Iorio et al. 2023) to produce the tables needed for the interpolation in SEVN. We used COSMORATE (https://gitlab.com/Filippo.santoliquido/cosmo_rate_public) (Santoliquido et al. 2020, 2021) to compute the rate densities. This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (The pandas development Team 2024). For the plots we used MATPLOTLIB (Hunter 2007).

References

  1. Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4 [NASA ADS] [Google Scholar]
  2. Artale, M. C., Mapelli, M., Bouffanais, Y., et al. 2020, MNRAS, 491, 3419 [Google Scholar]
  3. Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379 [Google Scholar]
  4. Belkus, H., Van Bever, J., & Vanbeveren, D. 2007, ApJ, 659, 1576 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  7. Boco, L., Lapi, A., Goswami, S., et al. 2019, ApJ, 881, 157 [Google Scholar]
  8. Boco, L., Mapelli, M., Sander, A. A. C., et al. 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202556187 [Google Scholar]
  9. Böhm-Vitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
  10. Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bressan, A. G., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25 [Google Scholar]
  12. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
  13. Briel, M. M., Eldridge, J. J., Stanway, E. R., Stevance, H. F., & Chrimes, A. A. 2022, MNRAS, 514, 1315 [NASA ADS] [CrossRef] [Google Scholar]
  14. Briel, M. M., Metha, B., Eldridge, J. J., Moriya, T. J., & Trenti, M. 2024, MNRAS, 533, 3907 [Google Scholar]
  15. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [Google Scholar]
  17. Cameron, A. J., Katz, H., Rey, M. P., & Saxena, A. 2023, MNRAS, 523, 3516 [NASA ADS] [CrossRef] [Google Scholar]
  18. Charbonnel, C., Schaerer, D., Prantzos, N., et al. 2023, A&A, 673, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Chen, K.-J., Woosley, S., Heger, A., Almgren, A., & Whalen, D. J. 2014, ApJ, 792, 28 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 [Google Scholar]
  21. Cheng, S. J., Goldberg, J. A., Cantiello, M., et al. 2024, ApJ, 974, 270 [NASA ADS] [CrossRef] [Google Scholar]
  22. Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
  23. Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012 [Google Scholar]
  24. Cooke, J., Sullivan, M., Gal-Yam, A., et al. 2012, Nature, 491, 228 [Google Scholar]
  25. Costa, G., Girardi, L., Bressan, A., et al. 2019a, A&A, 631, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Costa, G., Girardi, L., Bressan, A., et al. 2019b, MNRAS, 485, 4641 [NASA ADS] [CrossRef] [Google Scholar]
  27. Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514 [NASA ADS] [CrossRef] [Google Scholar]
  28. Costa, G., Ballone, A., Mapelli, M., & Bressan, A. 2022, MNRAS, 516, 1072 [NASA ADS] [CrossRef] [Google Scholar]
  29. Costa, G., Mapelli, M., Iorio, G., et al. 2023, MNRAS, 525, 2891 [NASA ADS] [CrossRef] [Google Scholar]
  30. Costa, G., Shepherd, K. G., Bressan, A., et al. 2025, A&A, 694, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure [Google Scholar]
  32. Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]
  33. Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [Google Scholar]
  34. Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]
  35. de Jager, C., Nieuwenhuijden, H., & van der Hucht, K. A. 1988, Bull. Inform. Centre Donnees Stellaires, 35, 141 [Google Scholar]
  36. Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011 [Google Scholar]
  37. Dotter, A. 2016, ApJS, 222, 8 [Google Scholar]
  38. du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941 [Google Scholar]
  39. Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19 [NASA ADS] [CrossRef] [Google Scholar]
  40. Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058 [Google Scholar]
  41. Farag, E., Renzo, M., Farmer, R., Chidester, M. T., & Timmes, F. X. 2022, ApJ, 937, 112 [NASA ADS] [CrossRef] [Google Scholar]
  42. Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53 [Google Scholar]
  43. Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36 [CrossRef] [Google Scholar]
  44. Farrell, E., Groh, J. H., Hirschi, R., et al. 2021, MNRAS, 502, L40 [NASA ADS] [CrossRef] [Google Scholar]
  45. Fraley, G. S. 1968, Ap&SS, 2, 96 [NASA ADS] [CrossRef] [Google Scholar]
  46. Fremling, C., Miller, A. A., Sharma, Y., et al. 2020, ApJ, 895, 32 [NASA ADS] [CrossRef] [Google Scholar]
  47. Gabrielli, F., Lapi, A., Boco, L., et al. 2024, MNRAS, 534, 151 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gal-Yam, A., Mazzali, P., Ofek, E. O., et al. 2009, Nature, 462, 624 [NASA ADS] [CrossRef] [Google Scholar]
  49. Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Giersz, M., Leigh, N., Hypki, A., Lützgendorf, N., & Askar, A. 2015, MNRAS, 454, 3150 [Google Scholar]
  51. Gieles, M., Padoan, P., Charbonnel, C., Vink, J. S., & Ramírez-Galeano, L. 2025, MNRAS [arXiv:2501.12138] [Google Scholar]
  52. Goswami, S., Silva, L., Bressan, A., et al. 2022, A&A, 663, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945 [Google Scholar]
  54. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Harikane, Y., Ono, Y., Ouchi, M., et al. 2022, ApJS, 259, 20 [NASA ADS] [CrossRef] [Google Scholar]
  56. Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5 [NASA ADS] [CrossRef] [Google Scholar]
  57. Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56 [NASA ADS] [CrossRef] [Google Scholar]
  58. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  59. Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532 [Google Scholar]
  60. Hendriks, D. D., van Son, L. A. C., Renzo, M., Izzard, R. G., & Farmer, R. 2023, MNRAS, 526, 4130 [NASA ADS] [CrossRef] [Google Scholar]
  61. Herwig, F. 2000, A&A, 360, 952 [NASA ADS] [Google Scholar]
  62. Hirschi, R. 2015, in Astrophysics and Space Science Library, 412, Very Massive Stars in the Local Universe, ed. J. S. Vink, 157 [Google Scholar]
  63. Humphreys, R. M., & Davidson, K. 1979, ApJ, 232, 409 [Google Scholar]
  64. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  65. Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426 [NASA ADS] [CrossRef] [Google Scholar]
  66. Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kojima, T., Ouchi, M., Rauch, M., et al. 2021, ApJ, 913, 22 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kozyreva, A., & Blinnikov, S. 2015, MNRAS, 454, 4357 [NASA ADS] [CrossRef] [Google Scholar]
  69. Kozyreva, A., Kromer, M., Noebauer, U. M., & Hirschi, R. 2018, MNRAS, 479, 3106 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  71. Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Ledoux, P. 1947, ApJ, 105, 305 [NASA ADS] [CrossRef] [Google Scholar]
  73. Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72 [Google Scholar]
  74. Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13 [NASA ADS] [CrossRef] [Google Scholar]
  75. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  76. Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [Google Scholar]
  77. Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  78. Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [NASA ADS] [CrossRef] [Google Scholar]
  79. Marchant, P., & Moriya, T. J. 2020, A&A, 640, L18 [EDP Sciences] [Google Scholar]
  80. Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 [NASA ADS] [Google Scholar]
  81. Mazzali, P. A., Moriya, T. J., Tanaka, M., & Woosley, S. E. 2019, MNRAS, 484, 3451 [NASA ADS] [CrossRef] [Google Scholar]
  82. Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 [Google Scholar]
  83. Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 [Google Scholar]
  84. Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Ohkubo, T., Umeda, H., Maeda, K., et al. 2006, ApJ, 645, 1352 [NASA ADS] [CrossRef] [Google Scholar]
  86. Pauli, D., Oskinova, L. M., Hamann, W. R., et al. 2025, A&A, 697, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  88. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  89. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  90. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  91. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  92. Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, ApJ, 904, 35 [NASA ADS] [CrossRef] [Google Scholar]
  93. Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648 [Google Scholar]
  94. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525 [Google Scholar]
  96. Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899 [Google Scholar]
  97. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [Google Scholar]
  98. Portegies Zwart, S. F., Dewi, J., & Maccarone, T. 2005, Ap&SS, 300, 247 [Google Scholar]
  99. Quimby, R. M., Gal-Yam, A., Arcavi, I., et al. 2011, Astronomer’s Telegram, 3841, 1 [Google Scholar]
  100. Renzo, M., Cantiello, M., Metzger, B. D., & Jiang, Y. F. 2020, ApJ, 904, L13 [NASA ADS] [CrossRef] [Google Scholar]
  101. Sabhahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736 [NASA ADS] [CrossRef] [Google Scholar]
  102. Sabhahit, G. N., Vink, J. S., Sander, A. A. C., & Higgins, E. R. 2023, MNRAS, 524, 1529 [NASA ADS] [CrossRef] [Google Scholar]
  103. Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873 [Google Scholar]
  104. Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Santoliquido, F., Mapelli, M., Bouffanais, Y., et al. 2020, ApJ, 898, 152 [Google Scholar]
  106. Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877 [NASA ADS] [CrossRef] [Google Scholar]
  107. Santoliquido, F., Mapelli, M., Artale, M. C., & Boco, L. 2022, MNRAS, 516, 3297 [NASA ADS] [CrossRef] [Google Scholar]
  108. Santoliquido, F., Mapelli, M., Iorio, G., et al. 2023, MNRAS, 524, 307 [NASA ADS] [CrossRef] [Google Scholar]
  109. Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202 [NASA ADS] [Google Scholar]
  110. Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
  111. Schootemeijer, A., Langer, N., Lennon, D., et al. 2021, A&A, 646, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Schulze, S., Fransson, C., Kozyreva, A., et al. 2024, A&A, 683, A223 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Schwarzschild, M. 1958, Structure and evolution of the stars [Google Scholar]
  114. Sgalletta, C., Mapelli, M., Boco, L., et al. 2025, A&A, 698, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Song, H., Meynet, G., Li, Z., et al. 2020, ApJ, 892, 41 [Google Scholar]
  116. Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739 [Google Scholar]
  117. Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [Google Scholar]
  118. Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]
  119. Takahashi, K. 2018, ApJ, 863, 153 [Google Scholar]
  120. Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021, ApJ, 910, 30 [NASA ADS] [CrossRef] [Google Scholar]
  121. Tanikawa, A., Moriya, T. J., Tominaga, N., & Yoshida, N. 2023, MNRAS, 519, L32 [Google Scholar]
  122. The pandas development Team 2024, pandas-dev/pandas: Pandas [Google Scholar]
  123. Törniamenti, S., Mapelli, M., Boco, L., et al. 2025, arXiv e-prints [arXiv:2510.12465] [Google Scholar]
  124. Vink, J. S. 2023, A&A, 679, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Vink, J. S., & Gräfener, G. 2012, ApJ, 751, L34 [Google Scholar]
  126. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] [Google Scholar]
  127. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146 [NASA ADS] [CrossRef] [Google Scholar]
  130. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  131. Winch, E. R. J., Sabhahit, G. N., Vink, J. S., & Higgins, E. R. 2025, MNRAS, 540, 90 [Google Scholar]
  132. Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
  133. Woosley, S. E., & Heger, A. 2021, ApJ, 912, L31 [NASA ADS] [CrossRef] [Google Scholar]
  134. Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015 [NASA ADS] [CrossRef] [Google Scholar]
  135. Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  136. Yang, M., Bonanos, A. Z., Jiang, B., et al. 2023, A&A, 676, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Yokoyama, T., & Tominaga, N. 2012, in Astronomical Society of the Pacific Conference Series, 458, Galactic Archaeology: Near-Field Cosmology and the Formation of the Milky Way, eds. W. Aoki, M. Ishigaki, T. Suda, T. Tsujimoto, & N. Arimoto, 51 [Google Scholar]
  138. Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351 [Google Scholar]
  139. Yungelson, L. R., van den Heuvel, E. P. J., Vink, J. S., Portegies Zwart, S. F., & de Koter, A. 2008, A&A, 477, 223 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114 [NASA ADS] [CrossRef] [Google Scholar]
1

We refer to the MESA inlist and run_star_extras files by S23, which are publicly available at this gitlab repository.

2

We used version 2.10.1 of SEVN, updated to commit hash c9dc8e4578990fd85cba91e082ebff6068388f56. SEVN is publicly available at this gitlab repository.

3

COSMOℛATE is publicly available at this gitlab repository.

Appendix A MESA models

Here, we summarize the main properties of our MESA models, including the evolution of luminosity (L) and effective temperature (Teff) in the Hertzsprung-Russell (HR) diagram (Fig. A.1), the radius (Fig. A.2) and the He core mass (Fig. A.3) until the end of core-He burning.

Figure A.1 shows the tracks in the HR diagram. All the models at low metallicity (Z ≤ 0.001) do not develop optically thick winds during the main sequence. These models evolve horizontally in the HR diagram towards cooler temperatures and the mass loss is not high enough to cause a drop in luminosity. Certain models at higher metallicity, e.g. the MZAMS = 100 M⊙ model at Z = 0.008, go back to the blue part of the HR diagram after entering the high-ΓEdd phase. At these metallicities, the more massive stars enter the high-ΓEdd regime just after the ZAMS, delineating the boundary between two distinct evolutionary behaviors: these stars evolve vertically in the HR diagram, exhibiting a steep decline in luminosity while remaining within a narrow range of temperature.

The radii of stars that enter the enhanced wind regime immediately after the ZAMS rapidly decrease as a consequence of the large amount of mass loss: indeed, in these cases, the H-rich envelope is quickly lost and also mass from the inner regions of the star is expelled (Fig. A.2). In the models that switch to high-ΓEdd regime during the main-sequence, the radius reaches the peak at ~ 2 Myr and then decreases. When the models enter this regime in the very last part of the main sequence, the mass loss is not able to contrast the inflation and the radius grows up to log(R/R⊙) ~ 3.5 in a short timescale. In all the models, the radii decrease rapidly after the ignition of the He core.

Figure A.3 shows the evolution of the He-core mass, while Table A.1 reports the final values of the He core mass, at the end of core-He burning. At Z = 0.008, the mass of the He core decreases quickly because of wind mass loss, and all models end their life with 13.7 M⊙ ≤ MHe,f < 20 M⊙. These systems cannot enter the pair-instability regime and undergo a core-collapse supernova. Lowering the metallicity, at Z = 0.004, only stars with 150 M⊙ ≤ MZAMS ≤ 200 M⊙ can experience PPISNe, while in all the other cases we find MHe,f < 32 M⊙. At Z = 0.001, stars can experience different fates: only the MZAMS = 50 M⊙ model undergoes a core-collapse supernova; models with MZAMS = 75, 100 M⊙ go through PPISN, while models with 125 M⊙ ≤ MZAMS ≤ 450 M⊙ can experience PISNe. Finally, if MZAMS = 475, 500 M⊙ the model directly collapses onto a black hole. At even lower metallicity (Z = 0.0008), the initial stellar mass range of systems that can undergo a PISN is reduced: 125 M⊙ ≤ MZAMS ≤ 350 M⊙.

thumbnail Fig. A.1

HR diagram of the simulated models. We report four selected metallicities: Z = 0.008, 0.004, 0.001, and 0.0008.
thumbnail Fig. A.2

Time evolution of the radius. We report four selected metallicities: Z = 0.008, 0.004, 0.001, and 0.0008. The solid lines show the evolutionary phase in which the stars enter the optically thick regime (ΓEdd > ΓEdd,tr), while the dashed lines represent the optically thin phase of stellar winds.
thumbnail Fig. A.3

He core mass evolution as a function of time. From top-left to bottom-right: Z = 0.008, 0.004, 0.001, and 0.0008.
Table A.1

Final He core mass of single stellar models.

Appendix B PARSEC models

Our comparison sample is based on PARSEC stellartracks, which have been extensively discussed by Costa et al. (2025); see also Bressan et al. (2012); Chen et al. (2015); Costa et al. (2022); Nguyen et al. (2022); Iorio et al. (2023); Costa et al. (2023) for additional references. Here, we briefly summarize their main properties.

For the winds of massive O-type stars, PARSEC adopts the fitting formulas by Vink et al. (2000) and Vink et al. (2001), with a correction introduced by Chen et al. (2015) to account for the effects of electron scattering, based on the models of Gräfener & Hamann (2008). In this framework, the mass-loss rate of an O-type star scales as Ṁ ∝ Zβ where Z is the metallicity and the exponent β depends on the Eddington ratio ΓEdd: β = 0.85 for ΓEdd < 2/3, β = 2.45–2.4 × ΓEdd for 2/3 ≤ ΓEdd < 1, and β = 0.05 for ΓEdd ≥ 1 (Chen et al. 2015). This formalism provides a fit to the models by Gräfener & Hamann (2008) and Vink et al. (2011), incorporating the dependence of mass loss on the Eddington ratio.

For Wolf–Rayet (WR) stars, PARSEC adopts the prescriptions by Sander et al. (2019), which successfully reproduce the observed properties of Galactic WR stars of type C (WC) and type O (WO).

PARSEC models include internal mixing by adopting the mixing-length theory (MLT; Böhm-Vitense 1958), calibrated on the Sun with a mixing-length parameter αMLT = 1.74 (Bressan et al. 2012). Convective regions are identified using the Schwarzschild criterion (Schwarzschild 1958), as implemented in Bressan et al. (1981). The PARSEC tracks used in this work assume an overshooting parameter of λov = 0.5, where λov represents the mean free path of convective elements beyond the convective boundary, expressed in units of the pressure scale height. In addition, an envelope undershooting distance of Λenv = 0.7 pressure scale heights is adopted.

All Tables

Table 1

PISN and PPISN rate densities at redshiſt ɀ = 0.

In the text
Table A.1

Final He core mass of single stellar models.

In the text

All Figures

thumbnail Fig. 1

Mass-loss rate of our MESA models as function of time, from the ZAMS until the end of core He burning. The models adopt the winds by Sabhahit et al. (2023). From top to bottom: Z = 0.008, 0.004, 0.001, and 0.0008. The solid lines correspond to the evolutionary phase in which the stellar models are in the optically thick regime (ΓEdd > ΓEdd,tr), while the dashed lines mark the optically thin regime.
In the text
thumbnail Fig. 2

Contour plot showing He core masses at the end of He burning for our MESA tracks adopting S23 winds as a function of the metallicity and ZAMS mass. The white lines highlight the levels at 32 M⊙. 64 M⊙, and 135 M⊙, which set the boundaries to have PPISNe, PISNe, and direct collapse via photodisintegration.
In the text
thumbnail Fig. 3

Same as Fig. 2, but for models C15.
In the text
thumbnail Fig. 4

Final He core mass as a function of the ZAMS mass. The orange lines are our new MESA tracks (models S23), while the dashed blue lines correspond to results from PARSEC tracks (models C15). Dark (light) gray shaded areas: regime for PISNe (PPISNe). From top left to bottom right, we show metallicities Z = 0.0001, 0.0004, 0.0008, 0.001, 0.004, and 0.008.
In the text
thumbnail Fig. 5

PISN (orange lines) and PPISN (blue lines) rate densities as function of redshiſt. The solid and dashed lines are obtained using the models by S23 and C15, respectively. The green shaded area represents the observational constraints set by Schulze et al. (2024). The black arrow highlights the upper limit from Schulze et al. (2024).
In the text
thumbnail Fig. 6

PISN (upper panels) and PPISN (lower panels) rate density evolution as function of redshift. The left-hand and right-hand plots show the results obtained with the new wind model by S23 and the one by C15, respectively. The thick black line represents the total rate density, whereas the colored lines show the contribution of individual metallicities from Z = 1 × 10−4 (blue dashed line) to Z = 2 × 10−2 (solid yellow line).
In the text
thumbnail Fig. 7

PISN (upper row) and PPISN (lower row) as function of redshift, varying the metallicity spread from σZ = 0.1 (black dashed line) to σZ = 0.6 (yellow dashed line). Left hand (right hand) panels: Models with stellar winds from S23 (from C15).
In the text
thumbnail Fig. A.1

HR diagram of the simulated models. We report four selected metallicities: Z = 0.008, 0.004, 0.001, and 0.0008.
In the text
thumbnail Fig. A.2

Time evolution of the radius. We report four selected metallicities: Z = 0.008, 0.004, 0.001, and 0.0008. The solid lines show the evolutionary phase in which the stars enter the optically thick regime (ΓEdd > ΓEdd,tr), while the dashed lines represent the optically thin phase of stellar winds.
In the text
thumbnail Fig. A.3

He core mass evolution as a function of time. From top-left to bottom-right: Z = 0.008, 0.004, 0.001, and 0.0008.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.