Pairinstability supernovae PISNe have crucial implications for many astrophysical topics including the search for very massive stars the black hole mass spectrum and galaxy chemical enrichment. To this end we need to understand where PISNe are across cosmic time and what are their favourable galactic environments. We present a new determination of the PISN rate as a function of redshift obtained by combining uptodate stellar evolution tracks from the PARSEC and FRANEC codes with an uptodate semiempirical determination of the star formation rate and metallicity evolution of starforming galaxies throughout cosmic history. We find the PISN rate to exhibit a huge dependence on the model assumptions including the criterion to identify stars unstable to pair production and the upper limit of the stellar initial mass function. Remarkably the interplay between the maximum metallicity at which stars explode as PISNe and the dispersion of the galaxy metallicity distribution dominates the uncertainties causing a  sevenordersofmagnitude PISN rate range. Furthermore we show a comparison with the corecollapse supernova rate and study the properties of the favourable PISN host galaxies. According to our results the main contribution to the PISN rate comes from metallicities between inlineformulatexmath idTM0001 notationLaTeXsim 103texmathinlineformula and inlineformulatexmath idTM0002 notationLaTeX102texmathinlineformula against the common assumption that views very low metallicity Population III stars as exclusive or dominant PISN progenitors. The strong dependencies we find offer the opportunity to constrain stellar and galaxy evolution models based on possible future or the lack of PISN observations.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
The cosmic rate of pairinstability supernovae
{'The cosmic rate of Pair-Instability Supernovae': "Francesco Gabrielli 1 , 2 ★ , Andrea Lapi 1 , 2 , 3 , 4 , Lumen Boco 1 , 2 , 3 , Cristiano Ugolini 1 , 8 , Guglielmo Costa 7 , Cecilia Sgalletta 1 , 5 , Kendall Shepherd 1 , 6 , Ugo N. Di Carlo 1 , Alessandro Bressan 1 , Marco Limongi 8 , 9 , 10 , Mario Spera 1 , 2 , 8  \n- 1 SISSA, via Bonomea 265, 34136 Trieste, Italy\n- 2 National Institute for Nuclear Physics - INFN, Sezione di Trieste, I-34127 Trieste, Italy\n- 3 Institute for fundamental physics of the Universe - IFPU, Via Beirut 2, 34014 Trieste, Italy\n- 4 Istituto di Radioastronomia - INAF/IRA, Via Piero Gobetti 101, 40129 Bologna, Italy\n- 5 INFN-Padova, Via Marzolo 8, I-35131 Padova, Italy\n- 6 Osservatorio Astronomico di Padova-INAF, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy\n- 7 Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France\n- 8 Istituto Nazionale di Astrofisica - Osservatorio Astronomico di Roma, Via Frascati 33, I-00040, Monteporzio Catone, Italy\n- 9 Kavli Institute for the Physics and Mathematics of the Universe, Todai Institutes for Advanced Study, University of Tokyo, Kashiwa, 277-8583, Japan 10 INFN. Sezione di Perugia, via A. Pascoli s/n, I-06125 Perugia, Italy  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'Pair-instability supernovae (PISNe) have crucial implications for many astrophysical topics, including the search for very massive stars, the black hole mass spectrum, and galaxy chemical enrichment. To this end, we need to understand where PISNe are across cosmic time, and what are their favourable galactic environments. We present a new determination of the PISN rate as a function of redshift, obtained by combining up-to-date stellar evolution tracks from the PARSEC and FRANEC codes, with an up-to-date semi-empirical determination of the star formation rate and metallicity evolution of star-forming galaxies throughout cosmic history. We find the PISN rate to exhibit a huge dependence on the model assumptions, including the criterion to identify stars unstable to pair production, and the upper limit of the stellar initial mass function. Remarkably, the interplay between the maximum metallicity at which stars explode as PISNe, and the dispersion of the galaxy metallicity distribution, dominates the uncertainties, causing a ∼ seven-orders-of-magnitude PISN rate range. Furthermore, we show a comparison with the corecollapse supernova rate, and study the properties of the favourable PISN host galaxies. According to our results, the main contribution to the PISN rate comes from metallicities between ∼ 10 -3 and 10 -2 , against the common assumption that views very-low-metallicity, Population III stars as exclusive or dominant PISN progenitors. The strong dependencies we find offer the opportunity to constrain stellar and galaxy evolution models based on possible future (or the lack of) PISN observations.  \nKey words: galaxies: general - supernovae: general - stars: evolution', '1 INTRODUCTION': "Pair-instability supernovae (PISNe) are explosions that develop inside the cores of very massive stars (VMSs) at the end of their evolution, leading to the complete disruption of the progenitor. The physical mechanism behind PISNe has been well understood ever since its discovery (Fowler & Hoyle 1964; Bisnovatyi-Kogan & Kazhdan 1967; Rakavy & Shaviv 1967; Barkat et al. 1967; Fraley 1968). At the end of C burning in the core, where temperatures approach 10 9 𝐾 and densities are greater than ∼ 100 𝑔 𝑐𝑚 -3 , photons become energetic enough to create electron-positron pairs. The pair-production process removes radiation pressure, which counteracts the gravitational pull from the inner layers, and lowers the adiabatic index Γ below 4/3. As a result, the star becomes unstable and begins to collapse in a runaway fashion. The onset of explosive O and Si burning releases  \nenergies of ∼ 10 52 -10 53 ergs (Heger & Woosley 2002), high enough to eject all the star's material, without leaving any remnant behind. The explosion produces high amounts of 56 Ni, up to ≳ 50 𝑀 ⊙ (Heger & Woosley 2002). Its radioactive decay is held responsible for luminosities up to 10 2 times those of typical core-collapse supernovae (CCSNe, e.g. Scannapieco et al. 2003; Kasen et al. 2011; Dessart et al. 2012; Whalen et al. 2013; Kozyreva et al. 2014a,b; Jerkstrand et al. 2015; Smidt et al. 2015; Kozyreva et al. 2016; Gilmer et al. 2017; Hartwig et al. 2018; Chatzopoulos et al. 2019). Despite PISNe being so luminous, and the several hundreds of CCSN observations achieved so far (e.g. Cooke et al. 2012; Yaron & Gal-Yam 2012; Gal-Yam et al. 2013; Guillochon et al. 2017), no PISN has ever been confidently discovered. Several candidate detections have been reported, including super-luminous supernovae (SLSNe), but none has been confirmed as a PISN (Woosley et al. 2007; Gal-Yam et al. 2009; Quimby et al. 2011; Gal-Yam 2012; Cooke et al. 2012; Kozyreva & Blinnikov 2015; Lunnan et al. 2016; Kozyreva et al. 2018; Mazzali", 'et al. 2019; Nicholl et al. 2020; Gomez et al. 2019; Schulze et al. 2023).': "PISNe are predicted to occur in stars with masses on the zeroage main sequence (ZAMS) in the range 140 ≲ 𝑀 ZAMS / 𝑀 ⊙ ≲ 260, and metallicities ( 𝑍 ) below some threshold (Heger & Woosley 2002; Heger et al. 2003). Stars undergo pair instability starting from 𝑀 ZAMS ∼ 100 𝑀 ⊙ , but below ∼ 140 𝑀 ⊙ they experience a series of pulsations, accompanied by the ejection of the most external layers, in a pulsational pair-instability supernova (PPISN, Heger & Woosley 2002; Heger et al. 2003). In this case, the core stays mostly intact, and the star finally collapses into a black hole (BH). Above ∼ 140 𝑀 ⊙ , the first pulsation is so energetic that it completely disrupts the star. For masses higher than ∼ 260 𝑀 ⊙ , the star directly collapses into an intermediate-mass BH (Heger & Woosley 2002; Heger et al. 2003). If 𝑍 is too high, mass loss due to stellar winds prevents the star from forming cores massive enough to become unstable (Heger et al. 2003). The maximum metallicity at which stars explode as PISNe is uncertain. Stellar evolution simulations succeed in producing PISNe up to some fraction of the solar metallicity, generally below ∼ 0 . 5 𝑍 ⊙ (Langer et al. 2007; Kozyreva et al. 2014b; Spera & Mapelli 2017; Langer 2012; Costa et al. 2021; Higgins et al. 2021; Sabhahit et al. 2023; Martinet et al. 2023).  \nStellar evolution codes simulate the evolution of stars from the ZAMS throughout the nuclear burning stages, providing a link between 𝑀 ZAMS and the final core mass, at the pre-SN stage. In order to identify stars undergoing pair instability, it is common to adopt a criterion based on the final mass of the core. As described above, the physical processes behind the onset of instability are much more complex. Nonetheless, the core mass criterion represents a good and useful proxy. By assuming that stars with He or CO core mass in a certain range end their life as PISN, it is possible to obtain a range of PISN progenitor masses. Depending on the details of the adopted stellar evolution code, the 𝑀 ZAMS range where PISNe occur can fluctuate (e.g. Heger & Woosley 2002; Heger et al. 2003; Langer et al. 2007; Takahashi et al. 2015; Spera & Mapelli 2017; Woosley 2017; Marchant et al. 2019; Iorio et al. 2022; Tanikawa et al. 2023). Moreover, detailed stellar evolution calculations show that the PISN range tends to shift to higher 𝑀 ZAMS at increasing metallicity (e.g. Heger & Woosley 2002; Heger et al. 2003; Spera & Mapelli 2017).  \nSince PISNe are expected to occur only in low-metallicity stars, pristine, very-low-metallicity Population III stars (Pop III) have been traditionally considered as main PISN progenitors (e.g. Ober et al. 1983; El Eid et al. 1983; Baraffe et al. 2001; Umeda & Nomoto 2001; Heger & Woosley 2002; Scannapieco et al. 2005; Wise & Abel 2005; Langer et al. 2007; Kasen et al. 2011; Pan et al. 2012; Dessart et al. 2012; Yoon, S.-C. et al. 2012; de Souza et al. 2013; de Souza et al. 2014; Whalen et al. 2013, 2014; Smidt et al. 2015; Magg et al. 2016; Regős et al. 2020; Venditti et al. 2023; Bovill et al. 2024; Wiggins et al. 2024). However, the fact that stellar evolution simulations allow for PISNe up to ∼ 0 . 5 𝑍 ⊙ , suggests that also higher𝑍 , Population II/I (Pop II/I) stars might give a significant contribution to the PISN rate.  \nComprehending the physics behind and the occurrence of PISNe, particularly in relation to their surrounding environments, holds a myriad of astrophysical implications.  \nFor instance, the study of PISNe is strongly linked to the debate on the upper limit of the stellar initial mass function (IMF). Stars with masses consistent with ≳ 200 -300 𝑀 ⊙ have indeed been observed in local galaxies (Crowther et al. 2010; Evans et al. 2010; Walborn et al. 2014; Crowther et al. 2016; Schneider et al. 2018; Crowther  \n2019; Bestenlehner et al. 2020; Kalari et al. 2022; Brands et al. 2022), challenging the previous consensus that placed the maximum stellar mass at ∼ 150 𝑀 ⊙ (Vink 2015). The existence of VMSs is also supported by chemical abundances studies (e.g. Romano et al. 2017, 2020a,b; Goswami et al. 2021, 2022). The observation of massive stellar BH binaries in gravitational waves (GWs) offers a new opportunity to investigate VMSs, even though all confident detections achieved so far do not necessarily require such massive progenitors (e.g. Abbot 2016, 2020; Spera et al. 2015; Spera & Mapelli 2017; Vink et al. 2021). Due to the very-high masses of PISN progenitors, the location of the IMF upper limit can critically determine the PISN rate.  \nMoreover, PISNe can help to shed light on many uncertain aspects of galaxy evolution. Indeed, PISN occurrence strongly depends on the properties of the galactic environments in which they take place. Therefore, their study requires a determination of the 𝑍 -dependent SFR density (SFRD) across cosmic history. This quantity defines the amount of mass available to form stars in the Universe, per unit time, comoving volume and metallicity. To estimate it, two main approaches are usually adopted, relying either on cosmological simulations (e.g. Mapelli et al. 2017; O'Shaughnessy et al. 2016; Lamberts et al. 2018; Mapelli & Giacobbo 2018; Artale et al. 2019), or on empirical prescriptions for the SFRD and galaxy metallicity distribution, derived from observations (e.g. Belczynski et al. 2016a; Lamberts et al. 2016; Cao et al. 2017; Elbert et al. 2017; Li et al. 2018; Boco et al. 2019, 2021; Chruslinska & Nelemans 2019; Neijssel et al. 2019; Santoliquido et al. 2021). Many uncertainties still exist around this subject (see e.g. Chruslinska & Nelemans 2019; Neĳssel et al. 2019; Boco et al. 2021; Santoliquido et al. 2021 for a comprehensive overview). In particular, the metallicity distribution of galaxies is still considerably unknown. Moreover, the low-mass end of the galaxy stellar mass functions (GSMFs), describing the mass distribution of galaxies, is still poorly constrained, since lowmass galaxies are very faint and thus difficult to observe. Specifically, it is not clear whether the GSMF slope in this mass range is constant, or redshift dependent (see e.g. Navarro-Carrera et al. 2023, where the latter case is supported by JWST data at redshifts 4 ≲ 𝑧 ≲ 8).  \nPISNe are also commonly invoked in studies about the chemical enrichment of the Universe (e.g. Ricotti & Ostriker 2004; Matteucci & Pipino 2005; Ballero et al. 2006; Cherchneff & Dwek 2009; Rollinde et al. 2009; Cherchneff & Dwek 2010), and in particular to explain the chemical abundance patterns observed in the Milky Way and local galaxies (e.g. Kojima et al. 2021; Goswami et al. 2021, 2022). The detection of PISN descendants, i.e. stars with chemical abundances compatible with at least partial enrichment by a PISN, represents another avenue to find out about PISN occurrence, besides direct observation (e.g. Heger & Woosley 2002; Aoki et al. 2014; Takahashi et al. 2018; Salvadori et al. 2019; Caffau et al. 2022; Aguado et al. 2023; Xing et al. 2023; Koutsouridou et al. 2023). Hence, obtaining insights into the existence and rate of PISNe across cosmic time would have strong implications for numerous unresolved inquiries in both astrophysics and cosmology.  \nIn this work, we compute the PISN rate as a function of redshift, by combining up-to-date stellar evolution tracks from the PARSEC and FRANEC codes, with an up-to-date semi-empirical determination of the 𝑍 -dependent SFRD across cosmic history. The aim is to study the dependence of the PISN rate on stellar and galactic prescriptions, an aspect which has not been explored extensively in the past. We also compute the ratio between PISN and CCSN rate, in order to provide a comparison with these observed transients. Finally, we explore the  \nproperties of the favourable PISN host galaxies. This work represents the premise to a follow-up work, where we will employ the theoretical framework presented here to study PISN observability with the James Webb Space Telescope (JWST), and address the question why these transients have never been observed. These works offer the opportunity to put constraints on stellar and galaxy evolution models, via the comparison with possible future PISN observations, or with the absence of observations in the eventuality that they are never discovered.  \nThe paper is structured as follows. In Section 2, we present our theoretical framework. We start by describing our semi-empirical, 𝑍 -dependent SFRD determination, and how we use stellar evolution tracks to compute the number of PISNe produced per unit star forming mass. We also show the considered variations on stellar and galactic prescriptions, and how we finally compute the PISN rate. Our results are presented in Section 3, and discussed in Section 4, where we also show additional variations. We draw our conclusions in Section 5.  \nThroughout this work, we assume the flat Λ CDMcosmology from Planck Collaboration 2020, with parameters Ω 𝑀 = 0 . 32, Ω 𝑏 = 0 . 05, 𝐻 0 = 67 𝑘𝑚𝑠 -1 𝑀𝑝𝑐 -1 . A standard Kroupa IMF is adopted (Kroupa 2001), defined from 0.1 𝑀 ⊙ . The IMF upper limit is among the parameters we decide to vary. Following Caffau et al. 2010, we use the value of 𝑍 ⊙ = 0 . 0153 for the Solar metallicity, and 12 + log(O/H) ⊙ = 8 . 76 for the Solar oxygen abundance.", '2 METHODS': 'We take into account how galaxies evolve throughout cosmic history by constructing a detailed semi-empirical determination of the 𝑍 -dependent SFRD, 𝑑 3 𝑀 SFR / 𝑑𝑡𝑑𝑉𝑑 log 𝑍 , directly based on observations, following Boco et al. 2021. This provides us with the amount of mass available for star formation at a certain redshift, per unit time, comoving volume and metallicity. In order to compute the number of PISNe produced per unit star forming mass, 𝑑𝑁 PISN / 𝑑𝑀 SFR , we make use of stellar evolution tracks computed with the PARSEC code (Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014; Chen et al. 2015; Costa et al. 2019, 2021; Nguyen et al. 2022). Finally, we compute the PISN rate density as a function of redshift by convolving these two quantities, according to the following formula:  \n𝑑 2 𝑁 PISN 𝑑𝑡𝑑𝑉 ( 𝑧 ) = ∫ 𝑑 log 𝑍 𝑑 3 𝑀 SFR 𝑑𝑡𝑑𝑉𝑑 log 𝑍 ( 𝑍, 𝑧 ) × 𝑑𝑁 PISN 𝑑𝑀 SFR ( 𝑍 ) . (1)  \nWedescribe in detail how we compute each of these quantities in the following Sections. In order to account for the uncertainties in stellar and galaxy evolution models, we decide to follow a variational approach, considering alternative stellar evolution codes, and different values for the relevant stellar and galactic parameters.', '2.1 Galaxy evolution': 'Wefollow the semi-empirical approach presented in Boco et al. 2021 in order to compute the 𝑍 -dependent SFRD, with some updates informed by more recent works (Chruślińska et al. 2021; Popesso et al. 2022). This approach is built on galaxy observations up to 𝑧 = 6. First, wecomputegalaxystellar mass functions (SMFs), providing a galaxy statistics based on their stellar mass, Φ ( 𝑀 ★ ) = 𝑑 2 𝑁 / 𝑑𝑉𝑑 log 𝑀 ★ . Here, 𝑁 indicates the number of galaxies, 𝑀 ★ their stellar mass, and  \n𝑉 the comoving volume. We then convolve the SMFs with the galaxy main sequence (MS), 𝜓 ( 𝑀 ★ ) , relating the galaxy stellar mass to their SFR, 𝜓 . We also implement a distribution of galaxies as a function of SFR, 𝑑𝑝 / 𝑑 log 𝜓 , as indicated by observations (see below). In this way we are able to obtain SFR functions, i.e. a galaxy statistics based on SFR. For the galaxy metallicity, we define a log-normal distribution, 𝑑𝑝 / 𝑑 log 𝑍 ( 𝑍, 𝑍 FMR ) , with a constant dispersion 𝜎 Z . The mean value, 𝑍 FMR , is given by the Fundamental Metallicity Relation (FMR) linking galaxy mass, metallicity and SFR. We finally convolve all these factors according to the following formula, where we integrate over 𝑀 ★ and 𝜓 :  \n𝑑 3 𝑀 SFR 𝑑𝑡𝑑𝑉𝑑 log 𝑍 ( 𝑍, 𝑧 ) = ∫ 𝑑 log 𝑀 ★ 𝑑 2 𝑁 𝑑𝑉𝑑 log 𝑀 ★ ( 𝑀 ★ , 𝑧 ) × ∫ 𝑑 log 𝜓 𝜓 𝑑𝑝 𝑑 log 𝜓 ( 𝜓, 𝑀 ★ , 𝑧 ) × 𝑑𝑝 𝑑 log 𝑍 ( 𝑍, 𝑍 FMR ( 𝑀 ★ , 𝜓 )) (2)  \nEach ingredient in Equation 2 is computed as described in the following.', 'SMF+MS': 'We follow the work by Chruslinska & Nelemans 2019 in order to compute the galactic SMFs. They perform a comprehensive determination by gathering several previous results from the literature, consisting in Schechter and/or double Schechter analytical fits to observational data. In particular, they consider discrete redshift bins, and in each of them they average among previous results in order to obtain the SMFs at the corresponding redshifts. The SMFs at arbitrary redshifts are computed by interpolating between these curves. In order to take into account the uncertainties around the low-mass end of the SMFs, they consider two variations, defining a constant low-mass end slope, 𝛼 SMF = -1 . 45, or prescribing a redshift dependence, 𝛼 SMF ( 𝑧 ) = -0 . 1 𝑧 -1 . 34. In Figure 1 we show the SMFs we obtain by following the same procedure, for both 𝛼 SMF cases (solid and dashed lines respectively), in the redshift range relevant to this work, 0 ≤ 𝑧 ≤ 6.  \nIt is to be noted that Chruslinska & Nelemans 2019 adopt a Kroupa IMF up to 100 𝑀 ⊙ , while we extend it up to 300 𝑀 ⊙ , as we will see later. This should not affect our results significantly, given that the Kroupa IMF predicts a relatively low number of massive stars between 100 and 300 𝑀 ⊙ . Therefore the SMFs are not expected to differ appreciably.  \nWe define our galactic MS following Popesso et al. 2022, which is the most complete determination up to date, taking into account all previous works in the literature and covering an unprecedented redshift and mass range, 0 < 𝑧 < 6 and 10 8 . 5 -10 11 . 5 𝑀 ⊙ . Additionally, we implement a double-Gaussian distribution in SFR, following Sargent et al. 2012:  \n𝑑𝑝 𝑑 log 𝜓 ( 𝜓, 𝑀 ★ , 𝑧 ) = 𝑓 MS 𝜎 MS √ 2 𝜋 exp " -( log 𝜓 - ⟨ log 𝜓 ⟩ MS ) 2 2 𝜎 2 MS # + + 𝑓 SB 𝜎 SB √ 2 𝜋 exp " -( log 𝜓 - ⟨ log 𝜓 ⟩ SB ) 2 2 𝜎 2 SB # , (3)  \nwhere ⟨ log 𝜓 ⟩ MS is the MS relation, and 𝜎 MS = 0 . 188. As one can note by looking at the second term, we also consider the presence  \nFigure 1. SMFs computed following Chruslinska & Nelemans 2019, at redshift 𝑧 = 0 to 6. We show both variations with 𝛼 SMF = -1 . 45 and 𝛼 SMF ( 𝑧 ) = -0 . 1 𝑧 -1 . 34, in solid and dashed lines respectively.  \n<!-- image -->  \nof starbursts (SBs), i.e. galaxies experiencing intense star formation, located in a separate region above the MS in the 𝜓 -𝑀 ★ plane. The meanofthe SB distribution is given by ⟨ log 𝜓 ⟩ SB = ⟨ log 𝜓 ⟩ MS + 0 . 59, while 𝜎 SB = 0 . 243. 𝑓 MS and 𝑓 SB are the fractions of MS and SB galaxies, respectively, such that 𝑓 MS + 𝑓 SB = 1. Differently from Boco et al. 2021, where a fixed SB fraction 𝑓 SB = 0 . 03 is assumed, we implement a dependence of 𝑓 SB on galaxy stellar mass and redshift, as done in Chruślińska et al. 2021. In Figure 2 we show the 𝑓 SB we compute following their work. This prescription enhances 𝑓 SB at low masses and increasing redshifts, bringing it from 0.03 up to 0.35. After 𝑧 = 4 . 4, 𝑓 SB saturates at values constant in 𝑧 . All in all, the 𝑀 ★ and 𝑧 dependencies enter in Equation 3 through ⟨ log 𝜓 ⟩ MS , ⟨ log 𝜓 ⟩ SB , 𝑓 MS and 𝑓 SB .', 'Z distribution': 'The FMR prescribes a correlation between galaxy stellar mass, SFR and metallicity, 𝑍 FMR ( 𝑀 ★ , 𝜓 ) . We define it following Curti et al. 2020:  \n𝑍 FMR ( 𝑀 ★ , 𝜓 ) = 8 . 779 - ( 0 . 31 / 2 . 1 ) × log GLYPH<16> 1 + ( 𝑀 ★ / 𝑀 0 ( 𝜓 )) -2 . 1 GLYPH<17> , (4)  \nwhere 𝑀 0 ( 𝜓 ) = 10 10 . 11 × 𝜓 0 . 56 .  \nFurthermore, we assume 𝑍 to follow a log-normal distribution around the FMR,  \n𝑑𝑝 𝑑 log 𝑍 ( 𝑍, 𝑍 FMR ( 𝑀 ★ , 𝜓 )) ∝ exp " -( log 𝑍 -log 𝑍 FMR ( 𝑀 ★ , 𝜓 )) 2 2 𝜎 2 Z # , (5)  \nWe consider variations 𝜎 Z = [ 0 . 15 , 0 . 35 , 0 . 70 ] , to study the effect  \nFigure 2. SB fractions 𝑓 SB as a function of 𝑀 ★ , for different 𝑧 , computed following Chruślińska et al. 2021. Thicker lines represent the 𝑓 SB computed for five initial redshifts (shown in the legend), while the 𝑓 SB at intermediate redshifts are computed by interpolation. At 𝑧 = 0, 𝑓 SB remains constant at 0.03, while above 𝑧 = 4 . 4 all curves saturate at 0.35, due to the dearth of observational data at those high redshifts.  \n<!-- image -->  \nof this parameter on the PISN rate.  \nWe report the final outcome of our galactic model in Figure 3, where we show the 𝑍 -dependent SFRD for all three variations on 𝜎 Z , fixing 𝛼 SMF = -1 . 45 (Equation 2). We also show the metallicity corresponding to the SFRD peak, as a function of redshift, for both variations on 𝛼 SMF (solid and dashed white lines), as well as the position of the overall SFRD peak (white star). In Figure 4 we instead report the SFRD as a function of redshift, obtained by integrating the previous quantity over 𝑍 , according to:  \n𝑑 2 𝑀 SFR 𝑑𝑡𝑑𝑉 ( 𝑧 ) = ∫ 𝑑 log 𝑍 𝑑 3 𝑀 SFR 𝑑𝑡𝑑𝑉𝑑 log 𝑍 ( 𝑍, 𝑧 ) . (6)  \nIn order to show the agreement of our SFRD with observations, we also plot several empirical determinations from galaxy surveys in different bands (Schiminovich et al. 2005; Gruppioni et al. 2013; Dunlop et al. 2016; Rowan-Robinson et al. 2016; Novak, M. et al. 2017; Casey et al. 2018; Liu et al. 2018; Gruppioni et al. 2020; Bouwens et al. 2021), as well as the Kistler et al. 2009, 2013 data obtained from long gamma-ray burst observations.  \nDifferent ways to compute this quantity exist in the literature. Among empirical or semi-empirical approaches, one of the main alternatives is to rely on SFR functions (SFRFs), describing the number density of galaxies with a given SFR. They can be computed from galaxy UV and IR luminosity functions, using the conversion between luminosity and SFR. Moreover, it is possible to implement a 𝑍 -dependence on the SFRD using a Mass Metallicity Relation (MZR), linking galaxy mass and metallicity, instead of a FMR. Another common approach is to combine an existing SFRD determination with a standalone galaxy 𝑍 distribution in redshift. We refer to Boco et al. 2021 for a comparison between these different methods. Finally, we stress that one of the main advantages of empirical determinations like ours, with respect to those relying on cosmological simulations, is that they are free from theoretical assumptions, being directly informed by observations.  \nFigure 3. SFRD as a function of redshift and metallicity, for 𝜎 Z = [ 0 . 15 , 0 . 35 , 0 . 70 ] (from left to right), and 𝛼 SMF = -1 . 45. The colorbar shows the SFRD values. White lines indicate the metallicity at which the SFRD peaks, as a function of redshift, for 𝛼 SMF = -1 . 45 (solid lines), and also 𝛼 SMF = 𝛼 SMF ( 𝑧 ) (dashed lines), for comparison. White stars indicate the redshift corresponding to the SFRD peak.  \n<!-- image -->  \nFigure 4. Cosmic SFR density as a function of redshift, for 𝛼 SMF = -1 . 45 and 𝛼 SMF ( 𝑧 ) variations (solid and dashed black lines respectively). The data points and bands show the observational determinations available at different redshifts, obtained from galaxy surveys in different bands (see the text for references). We also show the Madau & Dickinson 2014 plot, corrected for our Kroupa IMF, for comparison (dotted grey line).  \n<!-- image -->', '2.2 Stellar evolution': 'The second part of our method consists in computing the number of PISNe produced per unit star forming mass available, 𝑑𝑁 PISN / 𝑑𝑀 SFR . We do so by integrating an assumed IMF, 𝜙 ( 𝑀 ) , over the mass range of PISN progenitors, according to the following formula:  \n𝑑𝑁 PISN 𝑑𝑀 SFR ( 𝑍 ) = ∫ 𝑀 exit Mentry 𝜙 ( 𝑀 ) 𝑑𝑀 ∫ 𝑀 up 0 . 1 𝑀𝜙 ( 𝑀 ) 𝑑𝑀 , (7)  \nwhere 𝑀 entry and 𝑀 exit are the ZAMS masses for entering and exiting the PISN range, respectively. We adopt a Kroupa IMF (Kroupa 2001) defined in the mass range [0 . 1 -𝑀 up] 𝑀 ⊙ , where the lower limit is the minimum stellar mass to ignite core H burning. We choose to vary 𝑀 up from 150 to 300 𝑀 ⊙ , in order to study the effects of  \ndifferent IMF upper limits on the PISN rate (see the Appendix B for an additional variation with 𝑀 up = 600 𝑀 ⊙ ).  \nIn order to compute 𝑀 entry and 𝑀 exit , we adopt the criterion on the final core mass, according to which a star will develop pair instability leading to a PISN explosion if its He or CO core mass ( 𝑀 He / CO ) in the pre-SN stage lies in a certain range. Then we employ stellar evolution tracks to link the initial, ZAMS mass to the final core mass, and obtain a 𝑀 ZAMS range for PISN progenitors. We gather several works studying the evolution of VMSs from the literature, where the authors consider a range of He and/or CO core masses for which stars end their life as PISN (Heger & Woosley 2002; Heger et al. 2003; Langer et al. 2007; Chatzopoulos & Wheeler 2012; Takahashi et al. 2015; Spera et al. 2016; Belczynski et al. 2016b; Spera & Mapelli 2017; Woosley 2017; Eldridge et al. 2018; Marchant et al. 2019; Farmer et al. 2019; Stevenson et al. 2019; Kinugawa et al. 2020; Belczynski 2020; Marchant & Moriya 2020; du Buisson et al. 2020; Tanikawa et al. 2021; Briel et al. 2022; Woosley & Heger 2021; Olejak et al. 2022; Briel et al. 2023; Tanikawa et al. 2023). Specifically, we decide to adopt the criterion on 𝑀 CO , instead of 𝑀 He , motivated by the fact that it is in the CO core that PI sets in. Moreover, varying C burning reaction rate, and consequently the C/O ratio at He core depletion, impacts significantly on the position of the star in the core density-temperature plane, i.e. on the onset of PI (e.g. Farmer et al. 2019). In any case, as shown in Appendix A, varying between CO and He core criterion does not significantly affect our results. Due to the uncertainties existing on these ranges, we consider an optimistic and pessimistic variation, featuring the widest and shortest intervals based on the literature, i.e. 𝑀 CO ∈ [ 60 -105 ] 𝑀 ⊙ and 𝑀 CO ∈ [ 45 -120 ] 𝑀 ⊙ , respectively. We also consider an intermediate case with 𝑀 CO ∈ [ 55 -110 ] 𝑀 ⊙ .  \nWe consider two sets of stellar evolution tracks computed with the PARSEC code, that follows the evolution of single stars from the pre-main sequence up to the most advanced burning phases. The first set (Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014; Chen et al. 2015) was used in Spera & Mapelli 2017, implemented in the binary population synthesis code SEVN , and has been employed in several works in the last years. We will address it as PARSEC-I . The second and more recent set, which we will call PARSEC-II (Costa  \net al. 2019, 2021), was implemented in the recently published version of SEVN (Iorio et al. 2022).  \nIn PARSEC-II , the nuclear reactions and elemental mixing are coupled and solved at the same time in a diffusive scheme (Marigo et al. 2013; Costa et al. 2019). Additionally, PARSEC-II has updated prescriptions for the mass loss of massive stars, including Wolf-Rayet (WR) type wind (Sander et al. 2019, see Costa et al. 2021 for more details). Moreover, these new models include an updated equation of state (EOS), also accounting for the effects of e -e + pair creation.  \nIn order to see the effects of varying stellar evolution code, we also consider the results obtained with the FRANEC stellar evolution code (Chieffi & Limongi 2013; Limongi & Chieffi 2018). Comparing the results obtained using different sets of evolutionary tracks, especially when they are calculated by different groups that use different codes, is not a trivial procedure because, besides differences already arising from the different numerical procedures adopted, there are also variations due to different assumptions in the model calculations, and even in the description of physical complex phenomena that still require a calibration against observations. For example, different assumptions could refer to the criterion (basically Schwarzschild or Ledoux) adopted to establish whether a region is stable or not to convective motions. Instead, the calibration process may refer to the calibration against the Solar model to fix the value of the Mixing Length Theory (MLT), or to the size of the overshooting adopted in the convective core or at the bottom of the convective envelope. It will be impossible to trace back all these differences (and their impact on the results), and we will list below only the major differences that we believe must be considered when performing a sound comparison, leaving the details to the individual papers that describe the model calculations.  \nPARSEC models use the Caffau et al. 2010 solar partition of heavy elements (with 𝑍 ⊙ = 0 . 0153), even at very-low metallicity, while FRANEC adopts the Asplund et al. 2009 solar partition (with 𝑍 ⊙ = 0 . 0134), and consider 𝛼 -element enhancement at low metallicity (see e.g. Limongi & Chieffi 2018). This small discrepancy can cause models to evolve at slightly different luminosities already on the main sequence (Sibony et al. 2024). Furthermore, since the metallicity of the Sun is used as a reference in the scaling law of the mass-loss process with metallicity, this small difference can induce differences in the mass-loss rates at different absolute metallicities, even adopting the same mass-loss prescriptions.  \nPARSEC-I provides tracks for 12 different metallicities from 1 × 10 -4 to 3 × 10 -2 , while PARSEC-II considers 13 metallicities from 1 × 10 -4 to 4 × 10 -2 . The FRANEC tracks are instead available for 4 metallicities, from 3 × 10 -5 to ∼ 1 . 35 × 10 -2 . We linearly interpolate to obtain stellar tracks at metallicities in between, following Iorio et al. 2022. Moreover, while the PARSEC-II tracks are computed for masses up to 600 𝑀 ⊙ , the PARSEC-I and FRANEC ones extend up to ≳ 300 and 120 𝑀 ⊙ respectively. We linearly extrapolate at higher masses. We caution that the actual trend of the stellar tracks at these high masses might deviate quite significantly from this approximation, which should be taken into account in examining our results.  \nFor the convective stability, PARSEC models adopt the Schwarzschild criterion, while FRANEC models adopt the Ledoux one. The latter is more restrictive than the former and, if we consider that the core overshooting prescriptions are slightly different between PARSEC and FRANEC , and that PARSEC also accounts for overshooting from the bottom of the convective envelope, it is clear the the evolution in the low-mass range of massive stars may be different  \n(but given the high non-linearity of the solutions, it is not easy to trace back the nature of the differences).  \nOn the other hand, the evolution of more massive stars is strongly regulated by mass-loss (Smith 2014), and the effect that may have the strongest impact on the results is the adopted prescription for the mass-loss rates in the different evolutionary phases. Concerning massive stars, we remind that PARSEC includes radiative winds depending on the mass and evolutionary phase as described in Chen et al. 2015 and Costa et al. 2021. In particular, in hot stars ( 𝑇 eff ≥ 10000 𝐾 ), it uses the mass-loss prescriptions by Vink et al. 2000, 2001, including a surface iron abundance dependence; it also includes a dependence on the Eddington ratio (Grafener & Hamann 2008; Vink et al. 2011), which becomes important for the most luminous stars. For WR stars, i.e. when 𝑇 eff ≥ 20000 𝐾 and the hydrogen surface abundance is less than 0.3 in mass fraction, PARSEC-II uses the luminosity-dependent prescription for the mass loss from Sander et al. 2019, while PARSEC-I uses a combination of literature massloss rates. Finally both PARSEC-I and II use the prescription by de Jager et al. 1988 for cold massive stars (i.e. red super giants, RSGs).  \nIn FRANEC models, the Eddington ratio is checked within each structure, and if a region where it is larger than unity is found, this region and the overlying layers are removed from the star. Furthermore, for RSGs a dust-driven mass-loss rate is used (van Loon et al. 2005). Both of these latter recipes for the mass-loss rates may cause strong differences in the evolution of the most massive stars, even for those with initial mass as low as 𝑀 ZAMS = 15 𝑀 ⊙ or 20 𝑀 ⊙ at near-solar metallicity.  \nAll stellar tracks employed here are calculated without rotation. The main effects of including rotation would be increased mass loss, and bigger cores due to enhanced chemical mixing. Overall, we expect the 𝑀 ZAMS ranges of PISN progenitors to shift to lower masses (see also the recent work by Winch et al. 2024), resulting in higher PISN rates. It would be interesting to evaluate the entity of these effects by employing evolution tracks of rotating stars, which we do not explore here.  \nIt is important to note that, as one can see, in this work we restrict to metallicities down to 10 -4 , and redshifts up to 𝑧 = 6, effectively only considering Pop II/I stars as PISN progenitors. This is motivated by the fact that huge uncertainties still exist around quantities such as the SFRD at very-high redshifts, 𝑧 > 6, and the Pop III IMF, which prevents from extending the study to Pop III stars in a comparably robust way. We discuss this issue in more depth in Section 4.6.  \nTable 1 shows the 𝑀 ZAMS ranges obtained with each stellar evolution code, at some representative metallicities. We note that 𝑀 up cuts these ranges above its value. Masses greater than 300 𝑀 ⊙ were fixed to 300 𝑀 ⊙ , since we only consider IMFs up to this value in this work (but see also the Appendix B). As one can see, these ranges can differ quite significantly from the canonical [ 140 -260 ] 𝑀 ⊙ . Moreover, they generally shift to higher masses at higher metallicities. Indeed, mass loss due to stellar winds is enhanced, thus the star must be initially more massive in order to produce a core satisfying the PISN criterion. Above a certain metallicity, which depends on the considered variation, the core cannot reach this mass threshold. Werefer to the maximum metallicity for a star to explode as PISN, as 𝑍 max. As we will later show, 𝑍 max turns out to be a crucial quantity in determining the PISN rate.  \nIn Figure 5, we report the 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) obtained using each stellar code, combined with different variations on 𝑀 CO criterion and 𝑀 up. To avoid redundancy, among the eighteen cases that we  \nTable 1. 𝑀 ZAMS ranges of PISN progenitors, for representative metallicities, obtained with the PARSEC-I , PARSEC-II , and FRANEC stellar evolution tracks. The ranges at metallicities in between the original ones from the stellar tracks are obtained via interpolation, as described in the text. All values are in solar units. Blanks indicate cases with no PISN, at 𝑍 > 𝑍 max. Masses > 300 𝑀 ⊙ are fixed to 300 𝑀 ⊙ , i.e. the highest 𝑀 up considered in this work (see Appendix B for an additional variation with 𝑀 up = 600 𝑀 ⊙ ). The double intervals at 𝑍 = 1 × 10 -4 of PARSEC-II are due to the non-monotonicity of the corresponding stellar track, which exits the PISN mass range and then re-enters again (see Figure 8 of Iorio et al. 2022).  \ncompute, we decide to select six representative combinations, described in Table 2. F stands for our fiducial variation. P represents our pessimistic case, producing the smallest amount of PISNe, while O is our optimistic case. All variations in the middle are marked with an 𝑀 . One can clearly see the 𝑍 max resulting from each combination, as the metallicity at which 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) vanishes. The exact values are reported in Table 2. The effects of combining different stellar codes, 𝑀 CO range and 𝑀 up, can be boiled down to 𝑍 max, and secondly to the height of the 𝑑𝑁 PISN / 𝑑𝑀 SFR curve. For this reason, when presenting our results in Section 3, every stellar variation will be represented by the corresponding 𝑍 max.  \nWe note that some different combinations would lead to results superimposing to the ones we show, due to the degeneracy of the PISN rate on stellar prescriptions, but falling in any case between our P and O variations.', '2.3 Rate computation': 'The galactic and stellar parts of our model provide us with the 𝑍 -dependent SFRD, 𝑑 3 𝑀 SFR / 𝑑𝑡𝑑𝑉𝑑 log 𝑍 ( 𝑍, 𝑧 ) , and the number of PISNe produced per unit star forming mass, 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) , respectively. In order to finally compute the PISN event rate as a function of redshift, 𝑑 2 𝑁 PISN / 𝑑𝑡𝑑𝑉 ( 𝑧 ) , we need to convolve these two quantities by integrating over 𝑍 , as shown in Equation 1.  \nIn Section 3.4, we study the dependence of the PISN rate on galaxy stellar mass and metallicity, 𝑀 ★ and 𝑍 . In order to do that, we need to start from the SFRD defined also per unit 𝑀 ★ , which we obtain by simply avoiding integrating over this quantity in Equation 2. We then obtain the dependence on 𝑀 ★ and 𝑍 by integrating over 𝑧 ,  \nTable 2. Considered variations on stellar evolution code, CO core mass criterion and IMF upper limit, as reference for Figure 5. The maximum metallicity to have PISN according to each variation, 𝑍 max, is also shown.  \nafter converting the unit volume into unit redshift via the comoving volume element per unit redshift and steradian, 𝑑𝑉 / 𝑑𝑧 :  \n𝑑 4 𝑀 SFR 𝑑𝑧𝑑𝑡𝑑 log 𝑍𝑑 log 𝑀 ★ ( 𝑍, 𝑀 ★ , 𝑧 ) = = 𝑑 4 𝑀 SFR 𝑑𝑡𝑑𝑉𝑑 log 𝑍𝑑 log 𝑀 ★ ( 𝑍, 𝑀 ★ , 𝑧 ) × 4 𝜋 𝑑𝑉 𝑑𝑧 ( 𝑧 ) , (8)  \n𝑑 3 𝑁 PISN 𝑑𝑡𝑑 log 𝑍𝑑 log 𝑀 ★ ( 𝑍, 𝑀 ★ ) = = ∫ 𝑑𝑧 𝑑 4 𝑀 SFR 𝑑𝑧𝑑𝑡𝑑 log 𝑍𝑑 log 𝑀 ★ ( 𝑍, 𝑀 ★ , 𝑧 ) × 𝑑𝑁 PISN 𝑑𝑀 SFR ( 𝑍 ) . (9)  \nBy further integrating, one can get the PISN rate as a function of 𝑀 ★ or 𝑍 . The possibility to handle the masses of single galaxies  \nFigure 5. Number of PISN produced per unit star forming mass as a function of metallicity, for different combinations of stellar evolution code, CO core mass criterion and IMF upper limit (see Table 2). This figure clearly shows the maximum metallicity at which a star can explode as PISN, 𝑍 max, according to each variation.  \n<!-- image -->  \nrepresents one of the main advantages of SFRD semi-empirical determinations like the one employed in this work. In particular, here it allows us to identify the masses of the favourable PISN host galaxies, besides their metallicity, as we show in Section 3.4.', '3 RESULTS': 'In the following Section, we present the results obtained in this work. First, we show how our variations on stellar evolution code, CO core mass criterion and IMF upper limit, affect the PISN rate (Section 3.1). Then we focus on the interplay between 𝑍 max, resulting from each stellar variation, and 𝜎 Z in determining the PISN rate, and the effect of changing 𝛼 SMF (Section 3.2). The aim is to study to what extent the uncertainties in both stellar and galactic models affect the PISN rate, and obtain a range of results spanning from a pessimistic to an optimistic case.  \nDifferently from PISNe, several hundreds of CCSNe have been observed so far. Thus it can be useful to compare the rates of these transients. We compute the ratio between PISN and CCSN rate, and show the results in Section 3.3, for all of our variations.  \nFinally, in Section 3.4 we study the dependence of the PISN rate on galaxy 𝑀 ★ and 𝑍 , in order to understand which galactic environments are favourable for PISN production, and how much galaxies with given properties contribute to the PISN rate.', '3.1 Stellar variations': 'Figure 6 shows the PISN rate density as a function of redshift, computed with Equation 1 for the stellar variations described in Section 2.2. The parameters of the galactic model are fixed to 𝜎 Z = 0 . 35 and 𝛼 SMF = -1 . 45, which can be considered our fiducial values. The rate closely resembles the trend of the SFRD (Figure 4), peaking at around 𝑧 = 2 and smoothly declining towards higher redshifts.  \nWe can see that our whole range of stellar variations, from pessimistic ( P ) to optimistic ( O ), leads to roughly three orders of magni-  \ntude in the PISN rate, with values at 𝑧 = 0 from ∼ 10 0 to 10 3 𝐺𝑝𝑐 -3 𝑦𝑟 -1 , and values at peak from ∼ 3 × 10 1 to 2 × 10 4 𝐺𝑝𝑐 -3 𝑦𝑟 -1 . As can be easily understood, a larger 𝑀 CO interval produces larger 𝑀 ZAMS ranges for the progenitors, i.e. more stars going into PISN. Analogously, extending the IMF from 𝑀 up = 150 to 300 𝑀 ⊙ adds the contribution from more massive stars. This translates into higher 𝑑𝑁 PISN / 𝑑𝑀 SFR , with increased 𝑍 max, allowing for the production of more PISNe, at higher metallicities (Figure 5). As a result, one gets a higher PISN rate. E.g. if we fix FRANEC as stellar evolution code and 𝑀 up = 150 𝑀 ⊙ , changing the 𝑀 CO criterion from [ 60 -105 ] to [ 45 -120 ] (i.e. variation 𝑃 to 𝑀 2) increases the rate by about two orders of magnitude. Selecting instead PARSEC-I , 𝑀 CO ∈ [ 55 -110 ] , and bringing 𝑀 up from 150 to 300 𝑀 ⊙ ( 𝑀 1 to 𝐹 ), leads to between one and two orders of magnitude increase in the rate.  \nThe effects described above are weaker going towards higher PISN rates. E.g. if we consider PARSEC-II with 𝑀 CO ∈ [ 45 -120 ] , bringing 𝑀 up from 150 to 300 𝑀 ⊙ ( 𝑀 3 to 𝑂 ) only leads to a factor ∼ 2 difference in the rate. This can be understood by looking at Figure 7, where we present the PISN rate distribution in the 𝑧 -𝑍 plane for variations 𝑃 , 𝐹 and 𝑂 (lower, middle and upper panels respectively). We decide to only select these variations for simplicity, being representative of the whole set. In particular, the central column shows results for 𝜎 Z = 0 . 35, considered here (other variations on 𝜎 Z are discussed in the next Section). These plots exhibit a sharp cut in the metallicity distribution, which is a direct effect of 𝑍 max, completely suppressing PISNe at higher metallicities. Variation 𝑂 displays 𝑍 max ≳ 10 -2 , thus including the majority of the contribution from the SFRD peak, located at metallicities just above or below 10 -2 depending on redshift (see the white lines in Figure 3). This can be best appreciated in Figure 7 by comparing the PISN rate distribution with the SFRD contour levels, indicated by violet dashed lines. Also variation 𝑀 3 includes a fair portion of the SFRD peak contribution, with only a sligthly smaller 𝑍 max. Together with the fact that 𝑑𝑁 PISN / 𝑑𝑀 SFR changes just by a factor less than 2, this causes the rate to increase only by a relatively small amount going from 𝑀 3 to 𝑂 . On the contrary, the PISN rate distribution for the more pessimistic variations ( 𝑃 and 𝑀 1) extends to metallicities far below the SFRD peak. As a consequence, moving to variations 𝑀 2 and 𝐹 adds a significant contribution to the PISN rate, from metallicities closer to the SFRD peak. This effect is further enhanced by the stronger increase in 𝑑𝑁 PISN / 𝑑𝑀 SFR , especially from 𝑃 to 𝑀 2, with respect to the more optimistic cases.  \nAs indicated by the stars in Figure 6, also the peak of the PISN rate distribution in redshift appears to shift with stellar variations, an effect which is enhanced when combining with galactic variations (see the following Section). It is due to the fact that, on average, metallicity tends to decrease with redshift, according to our 𝑍 -evolution recipe. Variations with lower 𝑍 max thus favour the contribution to the PISN rate coming from high redshift. As a result, the position of the peak tends to shift towards higher redshifts, always below 𝑧 = 3 considering only stellar variations.', '3.2 Galactic variations': 'After showing how different stellar evolution prescriptions, and IMF upper limits, influence the PISN rate, we now focus on the variations on parameters 𝜎 Z and 𝛼 SMF of our galactic model (see Section 2.1).  \nAs anticipated above, we find a strong interplay between 𝑍 max, resulting from the combination of stellar evolution code, 𝑀 CO criterion and 𝑀 up, and 𝜎 Z . Among the stellar variations discussed in  \nthe previous Section, for clarity we decide to select only the most pessimistic and optimistic cases ( 𝑃 and 𝑂 ), as well as the fiducial one ( 𝐹 ). Figure 8 shows the PISN rates we obtain for the combination of each of these variations, described by the corresponding 𝑍 max (see Table 2), with different 𝜎 Z . As we can see, this results in PISN rates spanning roughly seven orders of magnitude considering values at 𝑧 = 0, or five orders of magnitude considering values at 𝑧 = 6. Most of this range is due to variations with lowest 𝑍 max = 1 . 5 × 10 -3 . Indeed, as already discussed in the previous Section and shown in Figure 7, such low 𝑍 max removes the main contribution to the PISN rate, coming from metallicities close to the SFRD peak. Selecting a low 𝜎 Z means considering a SFRD which does not extend to the lowest metallicities (Figure 3, left panel), and as a consequence it strongly suppresses the PISN rate contribution from the tail of the galaxy metallicity distribution, 𝑍 < 𝑍 max. This is the reason why variation 𝑃 is so dependent on 𝜎 Z . On the other hand, variation F yields a 𝑍 max = 6 . 6 × 10 -3 which is closer to the SFRD peak, and even in the lowest 𝜎 Z case it includes part of that contribution (middle panels in Figure 7). Therefore the rate suffers less dramatically from changes in 𝜎 Z . Variation 𝑂 , with 𝑍 max = 1 . 7 × 10 -2 (top panels in Figure 7), includes most of the contribution from the SFRD peak, making the dependence of the PISN rate on 𝜎 Z even fainter. For all 𝑍 max, this dependence appears more enhanced at lower redshifts. This is due to the fact that, as shown in Figure 3, our SFRD experiences a substantial decrease below 𝑧 ∼ 1, resulting in a drop in the rate which is more and more significant going to lower 𝜎 Z .  \nIn variations 𝑃 and 𝐹 , the PISN rate increases with 𝜎 Z . Indeed, in these cases 𝑍 max lies below the peak of the SFRD, therefore a higher dispersion for the galaxy metallicity distribution increases the PISN rate below 𝑍 max. On the contrary, variation 𝑂 exhibits a rate increasing with decreasing 𝜎 Z (namely, the dotted line is above the solid one). This is because in this case the peak of the SFRD is already included below 𝑍 max. A higher 𝜎 Z thus spreads a fraction of the PISN rate distribution above 𝑍 max, where it gets lost (Figure 7).  \nAll in all, the PISN rate dependence on 𝜎 Z varies dramatically with 𝑍 max, increasing significantly going to lower values of this parameter (see the color bands in Figure 8, which get larger for lower 𝑍 max). This shows the strong interplay between these parameters, which turns out to be crucial in determining the PISN rate. It is to be noted that, if we applied our galactic variations to the other stellar variations ( 𝑀 1, 𝑀 2 and 𝑀 3), we would obtain ranges of results partially superimposing to the ones we show for variations 𝑃 , 𝐹 and 𝑂 , as expressed by the grey band in Figure 8. This reveals the degeneracy of the PISN rate on stellar and galactic prescriptions.  \nOne can also notice a change in the peak position, from 𝑧 ∼ 2 (reproducing the SFRD peak) for 𝜎 Z = 0 . 70, to 𝑧 ≳ 3 for 𝜎 Z = 0 . 15, in the lowest 𝑍 max = 1 . 5 × 10 -3 variation. As explained in the previous section, this is an effect of the average metallicity decrease with redshift, prescribed by our Z evolution recipe. Here this effect is enhanced by the 𝜎 Z = 0 . 15 variation, which further favours the PISN rate contribution from higher 𝑧 .  \nFigure 9 shows the effect of varying 𝛼 SMF prescription. The grey band represents the difference between the constant 𝛼 SMF case, 𝛼 SMF = -1 . 45, and that with 𝛼 SMF defined as a function of redshift, 𝛼 SMF = 𝛼 SMF ( 𝑧 ) . We consider only some of the variations reported in Figure 8, for clarity purposes, since the trend is analogous. We can see that prescribing a redshift dependence for 𝛼 SMF has an appreciable effect only at high redshifts, 𝑧 > 3 -4, where it increases the rate by a factor in any case smaller than one order of magnitude. Indeed, as shown in Figure 1, this variation increases  \nthe number of galaxies at low masses, an effect which is stronger going to higher redshift. Moreover, the FMR prescribes that lowmass galaxies are also metal poor. Since this relation produces an average metallicity decreasing with redshift, these low-mass galaxies end up increasing the PISN rate at high redshift. This can also be appreciated by looking at Figure 3, where one can see that the 𝑍 corresponding to the peak of the SFRD rate experiences a steeper decrease with redshift for variation 𝛼 SMF ( 𝑧 ) , with respect to constant 𝛼 SMF (an effect already shown in Chruslinska & Nelemans 2019). For the same reason, in variation 𝑃 the peak of the 𝛼 SMF ( 𝑧 ) rate is shifted to between 𝑧 = 4 and 5. Indeed, the effects of the 𝑍 max -𝜎 Z dependence are exacerbated, and in particular the PISN rate peak is brought to even higher redshifts than the constant𝛼 SMF case.  \nWe warn that the huge range we obtain for the PISN rate, by accounting for both stellar and galactic variations, is actually strongly dependent on the adopted 𝑍 -evolution prescriptions. We show this point in Section 4.5, where we explore an additional variation on FMR.  \nFinally, Figure 10 shows the uncertainty intervals on the local PISN rate, computed at 𝑧 = 0, due to each individual stellar and galactic variation (red and blue bands respectively). Uncertainties due to all stellar variations, as well as all possible 𝑍 max -𝜎 Z combinations, are also shown (purple band). Every bar has been computed by fixing all other prescriptions to the fiducial case, and varying only the quantity of interest over the whole range. Notice how the 𝛼 SMF variation does not lead to significant uncertainties on the PISN rate at 𝑧 = 0 since, as explained in this Section, it produces appreciable effects only at 𝑧 ∼ 3 -4. This Figure clearly shows how crucial the interplay between 𝑍 max and 𝜎 Z is in determining the PISN rate, extending the possible values downwards by several orders of magnitude, with respect to considering stellar and galactic variations separately.', '3.3 PISN over CCSN ratio': 'In this Section, we compute the ratio between the PISN and CCSN rate as a function of redshift. We follow the same procedure outlined in Section 2, fixing 𝛼 SMF = -1 . 45. We modify the mass range for IMF integration in Equation 7 to [8-50] 𝑀 ⊙ . The upper ZAMS-mass limit of CCSN progenitors is highly uncertain, possibly ranging from ∼ 20 -25 𝑀 ⊙ to 125 𝑀 ⊙ (Heger & Woosley 2002; Heger et al. 2003; Dahlen et al. 2004; Cappellaro et al. 2005; Botticella et al. 2007; Smartt et al. 2009; Botticella et al. 2012; Dahlen et al. 2012; Mattila et al. 2012; Melinder et al. 2012; Strolger et al. 2015; Petrushevska et al. 2016; Ziegler et al. 2022). For this reason, we select an intermediate, fiducial value of 50 𝑀 ⊙ . We note that varying this value along the credible range does not change our results significantly, due to our bottom-heavy IMF. It is also to be stressed that, while for PISNe we obtain a dependence of the 𝑀 ZAMS range on metallicity, here we keep it constant, since the CCSN mass range is not expected to exhibit such a crucial dependence (namely CCSNe are expected to take place at all metallicities).  \nThe results are shown in Table 3, for the same variations as Figure 8. Analogously to the PISN rate, the PI/CC ratio spans from ∼ five to less than seven orders of magnitude, depending on redshift. It ranges from 2 . 5 × 10 -9 to 1 . 5 × 10 -2 at 𝑧 = 0, and from 2 . 3 × 10 -7 to 2 . 8 × 10 -2 at 𝑧 = 6. Considering the redshift at which the PISN rate peaks for each variation, 𝑧 PI peak , we get a PI/CC ratio ranging from 1 . 4 × 10 -7 to 2 . 2 × 10 -2 . The PISN rate dependence on 𝑍 max -𝜎 Z combinations is closely reproduced. This is easily understood, since  \nFigure 6. PISN rate as a function of 𝑧 for our stellar variations on stellar evolution code, CO core mass criterion and IMF upper limit. See Table 2 for a description of each variation. Stars indicate the peak of the PISN rate.  \n<!-- image -->  \nFigure 7. PISN rate as a function 𝑧 and 𝑍 for variations 𝑃 , 𝐹 and 𝑂 on stellar evolution code, 𝑀 CO criterion and 𝑀 up (lower, middle and upper panels respectively). See Table 2 for a description of each variation, and the corresponding 𝑍 max. Different columns correspond to variations 𝜎 Z = 0 . 15, 0.35, 0.70 (left to right). We also show the SFRD contour levels for comparison (violet dashed lines), corresponding to fractions of 10 -10 , 10 -2 , 0.1 and 0.9 times its peak value.  \n<!-- image -->  \nFigure 8. PISN rate as a function of 𝑧 for different 𝜎 Z, represented by different linestyles, for 𝑍 max corresponding to variations 𝑃 , 𝐹 and 𝑂 (blue, orange and green lines, see Table 2). Bands of the corresponding colors show the range of results produced by varying 𝜎 Z, for every fixed 𝑍 max. Stars indicate the peak of the PISN rate. We also plot a grey band to stress that, for clarity reasons, here we do not show all possible combinations of our model variations, but the ones we left out would fill in the gaps between, and be included within, the representative ones that we selected.  \n<!-- image -->  \nFigure 9. PISN rate as a function of 𝑧 for 𝛼 SMF = -1 . 45 and 𝛼 SMF = 𝛼 SMF ( 𝑧 ) variations. We show the difference between the two variations with a grey band. Lines are as in Figure 8 (note that we only reported some of those variations, for clarity). Stars indicate the peak of the PISN rate for each variation.  \n<!-- image -->  \nthe CCSN rate consists simply in the multiplication of the SFRD by the constant factor 𝑑𝑁 CCSN / 𝑑𝑀 SFR . Therefore, all the features of the PISN rate, which are ultimately due to its metallicity dependence, are again found in the PI/CC ratio. Values at 𝑧 = 6 are always higher than at 𝑧 = 0 and even 𝑧 PI peak . This is because PISNe occur in environments with 𝑍 < 𝑍 max, which are found preferentially at high redshift, according to our 𝑍 -evolution recipe. Since CCSNe are independent from metallicity, this favours PISNe over CCSNe, and causes their ratio to increase with redshift. The PISN rate drop at low redshifts, which is again due to the metallicity dependence and is thus absent in the CCSN rate, further accentuates this feature.', '3.4 Host galaxy properties': 'In this Section we explore the dependencies of the PISN rate on galaxy mass and metallicity, and the interplay between them. Among our set of variations, we select three cases which best serve this  \npurpose. The rates have been computed with Equations 8 and 9, further integrating over 𝑀 ★ or 𝑍 . We stress that, due to the integration over 𝑧 in Equation 9, these rates represent the contribution coming from all redshifts up to 𝑧 = 6.  \nIn Figure 11 we present a corner plot showing the individual and joint dependencies of the PISN rate on 𝑀 ★ and 𝑍 , for our fiducial stellar variation 𝐹 , with galactic parameter 𝜎 Z = 0 . 15. The difference between the two GSMF low-mass end slope cases, 𝛼 SMF = -1 . 45 and 𝛼 SMF = 𝛼 SMF ( 𝑧 ) , is indicated by a grey band. We find it informative to also show the SFR for the corresponding variations, with 𝛼 SMF = -1 . 45, telling us about the total star forming mass available in the first place. In particular, we show both the contour levels for the SFR as a function of 𝑀 ★ and 𝑍 , as well as the individual dependencies of the SFR on each of these quantities (violet dashed lines). We find again the sharp cut in the metallicity distribution due to 𝑍 max, already discussed for Figure 7. In this case 𝑍 max = 6 . 6 × 10 -3 , which lies below the peak of the SFRD, as indicated by the violet dashed contours in the bottom left panel (see also Figure 3). As a  \nFigure 10. Uncertainties on the local PISN rate, computed at 𝑧 = 0, due to each considered variation on stellar and galactic prescriptions, indicated with a red and blue band respectively (see the text). Global uncertainties due to all stellar variations combined are also shown, as well as those due to all possible combinations of 𝑍 max and 𝜎 Z (purple band). Indeed, as highlighted in the text, the interplay between these two parameters plays a crucial role in determining the PISN rate. This figure was inspired by Figure 5 in Farmer et al. 2019.  \n<!-- image -->  \nTable 3. PISN over CCSN rate ratio for our set of 𝑍 max -𝜎 Z combinations (same as Figure 8). We report values at 𝑧 = 0 and 𝑧 = 6, as well as the redshift at which the PISN rate peaks ( 𝑧 PI peak ), in order to show the whole range across redshift.  \nconsequence, the main contribution to the PISN rate, coming from metallicities where the SFRD peaks, gets partially cut out. As discussed in Sections 3.1 and 3.2, this is the reason why the PISN rate turns out to be lower with respect to other variations with higher 𝑍 max.  \nAs one can see from the right panel of Figure 11, the PISN rate peaks at metallicities just below 𝑍 max, as a result of combining 𝑑𝑁 PISN / 𝑑𝑀 SFR (variation 𝐹 in Figure 5) with the SFRD. In other  \nwords, PISN production is favoured in low𝑍 environments, but stars formedat 𝑍 closer to the SFRD peak are more abundant, and therefore provide a larger contribution to the PISN rate.  \nIn the top panel of Figure 11 we show the PISN rate dependence on 𝑀 ★ . This variation favours host galaxies with stellar masses between 10 9 and 10 10 𝑀 ⊙ , from which comes the main contribution to the PISN rate. These masses are somewhat lower than those at which the SFR peaks, around ∼ 10 10 𝑀 ⊙ , as indicated by the violet dashed line in the top panel. Indeed, partially cutting the peak of the SFRD distribution over metallicity, also stops the rise of the PISN rate with 𝑀 ★ , as can be clearly seen in the bottom left panel by comparing the PISN rate distribution with the SFR violet dashed contours. This is due to the correlation between 𝑀 ★ and 𝑍 prescribed by our FMR.  \nFigure 12 shows the results obtained for variation 𝑃 , with 𝜎 Z = 0 . 15. This is the most pessimistic case considered in this work. Variation 𝑃 features the lowest 𝑍 max = 1 . 5 × 10 -3 , leading to an even more dramatic cut in the metallicity distribution, with respect to the previous case shown. As a consequence, the rate ends up being more than four orders of magnitude lower than in the previous case. Moreover, the rate peak shifts to a lower metallicity, just above 10 -3 , reflecting the trend of the 𝑑𝑁 PISN / 𝑑𝑀 SFR over metallicity (variation 𝑃 in Figure 5). The mass distribution gets shifted towards lower values, favouring galaxies with 𝑀 ★ between 10 8 and 10 9 𝑀 ⊙ . This is again due to the correlation between 𝑀 ★ and 𝑍 , given the lower 𝑍 max cut. Since 𝑍 max is far below the SFRD peak, a major contribution to the PISN rate is taken out of the game. This is the reason why in this variation the rate experiences such a dramatic decrease, as discussed in Sections 3.1 and 3.2.  \nFinally, in Figure 13 we present the host galaxy properties obtained for variation 𝑂 , combined with 𝜎 Z = 0 . 35. As already discussed, here 𝑍 max is greater than the SFRD peak, thus the rate is increased with respect to the fiducial case, by a factor less than one order of magnitude. Moreover, the peaks in both the 𝑍 and 𝑀 ★ distributions resemble the SFR ones since, differently from the previous cases shown, here the main contribution to the SFR is almost-fully included in the PISN rate distribution. In particular, the peak lies at 𝑍 between 10 -2 and 10 -2 . 5 , and 𝑀 ★ ≲ 10 10 𝑀 ⊙ .  \nIn all cases, prescribing a redshift dependence of 𝛼 SMF increases the rate at 𝑀 ★ ≲ 10 10 𝑀 ⊙ , as can be seen by the grey band in the top panel. Indeed, the 𝛼 SMF ( 𝑧 ) variation produces more galaxies with low mass with respect to the fixed 𝛼 SMF case (Figure 1). Because of the correlation between 𝑀 ★ and 𝑍 , these galaxies will also be metal poor, explaining why the PISN rate slightly increases at metallicities below the SFRD peak (grey band in the lower right panel).', '4.1 PISN rate and variations': "Our results indicate that using different stellar evolution prescriptions, and varying the relevant parameters of the galactic model, hugely affects the PISN rate. Different stellar evolution codes, 𝑀 CO criteria and IMF upper limits produce a ∼ three orders of magnitude range from the most pessimistic to optimistic case (variations 𝑃 and 𝑂 in Figure 6). These variations play a crucial role through the quantity 𝑍 max, and secondly through the height of the 𝑑𝑁 PISN / 𝑑𝑀 SFR curve. As shown in Figure 8, the interplay between 𝑍 max and 𝜎 Z  \nFigure 11. PISN rate as a function of 𝑀 ★ , 𝑍 and both, for our fiducial stellar variation ( PARSEC-I , 𝑀 CO ∈ [ 55 -110 ] 𝑀 ⊙ , 𝑀 up IMF = 300 𝑀 ⊙ ), with 𝜎 Z = 0 . 15 and 𝛼 SMF = -1 . 45. Grey areas indicate the transition from the constant 𝛼 SMF to the 𝛼 SMF ( 𝑧 ) case. For comparison, we also show the SFR contour levels as a function of 𝑀 ★ and 𝑍 for the corresponding variations, with fixed 𝛼 SMF = -1 . 45, as well as the individual dependencies of the SFR on 𝑀 ★ and 𝑍 (violet dashed lines). In particular, we show the contours at 10 -10 , 10 -4 , 0.1 and 0.9 the SFR peak. As stressed in the text, these rates represent the contribution coming from all redshifts up to 𝑧 = 6 (see Equation 9).  \n<!-- image -->  \nbroadens the PISN rate range to between ∼ five and seven orders of magnitude, throughout the whole redshift range considered. Indeed, 𝑍 max completely cuts the SFRD distribution over metallicity, thus 𝜎 Z strongly regulates the amount of SFRD which gets saved or lost (Figure 7). This effect becomes huge for low 𝑍 max. E.g., in variation 𝑃 , the peak of the PISN rate drops from ∼ 10 1 to 10 -1 𝐺𝑝𝑐 -3 𝑦𝑟 -1 , by decreasing 𝜎 Z from 0.35 to 0.15. Moreover, the PISN rate peak is moved from 𝑧 ≲ 2 to 𝑧 ≳ 3, for the most pessimistic variation. Finally, varying 𝛼 SMF increases the PISN rate at 𝑧 > 3 -4 by a factor less than one order of magnitude, and shifts its peak to even higher redshifts than the constant𝛼 SMF case, up to 𝑧 ∼ 4 -5 (Figure 9).  \nThese results show how the existing uncertainties in stellar and galaxy evolution hamper the determination of the cosmic PISN rate, with uncertainties dominated by the interplay between 𝑍 max and 𝜎 Z . On the other hand, this offers the opportunity to constrain stellar and galactic models, thanks to possible PISN observations in the future. Also in the eventuality that PISNe are never discovered, their lack of observations would help pose limits on the models' parameters.  \nAccording to our PISN rate distributions over galactic 𝑀 ★ and 𝑍 , the stellar mass of the favourable PISN host galaxies ranges from  \n10 8 -10 9 𝑀 ⊙ , to around 10 10 𝑀 ⊙ , going from pessimistic to optimistic variation. Favoured galaxy metallicities range from ∼ 10 -3 to just below 10 -2 , following 𝑍 max. These results tell us about the galactic environments which are most favourable to host PISNe, and can be taken as indication for current and future surveys aimed at observing these transients.", '4.2 Comparison with previous works': 'It is useful to compare our results with previous works that, similarly to us, compute the PISN rate by combining a stellar evolution code with a cosmic star formation history and 𝑍 -evolution recipe (Stevenson et al. 2019; du Buisson et al. 2020; Briel et al. 2022; Hendriks et al. 2023; Tanikawa et al. 2023; Tanikawa 2024). We note that, in all these comparisons, stellar and galactic prescriptions were selected in order to be as compatible as possible. Moreover, these works also consider PISNe arising from binary stars, while we only consider single stars. As we will discuss in the next Section, we do not expect taking binaries into account to influence the PISN rate significantly. The 𝑍 -dependent SFRDs adopted in the following works exhibit 𝜎 Z  \nFigure 12. As in Figure 11, for our pessimistic variation ( FRANEC , 𝑀 CO ∈ [ 60 -105 ] 𝑀 ⊙ , 𝑀 up IMF = 150 𝑀 ⊙ ), with 𝜎 Z = 0 . 15. Notice the different axis ranges for the PISN rate, due to the low values.  \n<!-- image -->  \nequal, or roughly compatible, to our 𝜎 Z = 0 . 35 choice. Therefore, we compare with the results shown in Figure 6, where we apply stellar variations keeping 𝜎 Z = 0 . 35 fixed.  \nOverall, these previous determinations tend to fall in the lower end of our Figure 6 range. This is partially due to the employed stellar evolution codes, producing different 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) with respect to our work. Moreover, the interplay between 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) and the adopted 𝑍 -evolution prescription can greatly enhance the differences with our results. Remarkably, the PISN rate obtained by Stevenson et al. 2019 extends up to three orders of magnitude downwards our pessimistic case (variation 𝑃 in Figure 6). Indeed, on one hand their 𝑑𝑁 PISN / 𝑑𝑀 SFR ( 𝑍 ) is between one and two orders of magnitude lower than in our variation 𝑃 , at metallicities crucial for the PISN rate, close to 𝑍 max ∼ 2 × 10 -3 (compare our Figure 5 with their Figure 2). This is also due to the fact that they consider an IMF upper limit of 150 𝑀 ⊙ . Secondly, their galaxy metallicity distribution declines much more steeply in redshift with respect to ours. As explained in Sections 3.1 and 3.2 for our results, this can dramatically lower the PISN rate, besides favouring the contribution from higher redshifts. Interestingly, this is the same reason why e.g. in Tanikawa et al. 2023 the PISN rate turns out to peak at much higher redshifts than ours, around 𝑧 ∼ 6 -7. As already discussed, we obtain the same effect, albeit of smaller entity, by adopting a GSMF low-  \na steeper decrease of the galaxy metallicity with redshift, and causing the PISN rate to peak up to 𝑧 ∼ 4 -5(see Section 3.2). If we change the Curti et al. 2020 FMR prescription adopted here to that of Mannucci et al. 2010, we are even able to move the peak to 𝑧 > 6, as discussed below (Section 4.5). Therefore, a redshift dependence of 𝛼 SMF might motivate the employment of steeper 𝑍 -evolution prescriptions.  \nNoticeably, while like us most of previous works focus on 𝑍 ≥ 10 -5 -10 -4 , Tanikawa et al. 2023; Tanikawa 2024 also consider very-low-metallicity, Pop III stars. As a result, e.g. in Tanikawa 2024 the PISN rate features two distinct peaks, one around 𝑧 = 2 -3 due to Pop II/I stars, and the other at 𝑧 ≳ 12 due to Pop III. Depending on the considered variations, including different IMF upper limits, Tanikawa et al. 2023; Tanikawa 2024 find these two peak contributions to be roughly compatible, or the Pop III one to be dominant. We further discuss the contribution to the PISN rate from Pop III stars below.', '4.3 PI/CC ratio': 'In Section 3.3, we computed the ratio between PISN and CCSN rate, for our 𝑍 max -𝜎 Z combinations (see Table 3). Our range of results  \nFigure 13. As in Figure 11, for our optimistic variation ( PARSEC-II , 𝑀 CO ∈ [ 45 -120 ] 𝑀 ⊙ , 𝑀 up IMF = 300 𝑀 ⊙ ), with 𝜎 Z = 0 . 35.  \n<!-- image -->  \nexhibits an upper limit ≳ 10 -2 , which is close to the value typically considered in the literature, coming from the simple integration of the IMF over the expected PISN and CCSN progenitor mass ranges. Adopting our Kroupa IMF, and integrating over the typical [140 -260] 𝑀 ⊙ range for PISN, and [8 -50] 𝑀 ⊙ for CCSN, gives a value of ∼ 1 . 5 × 10 -2 for the PI/CC ratio. Considering the hundreds of observed CCSNe, and neglecting the different observational biases linked to these transients for the sake of simplicity, this would lead to expect at least 1 PISN detection so far, which is not the case. On the other hand, the lowest part of our range would tell us that PISNe are simply too rare with respect to CCSNe, making it trivial to understand why we never observed one. We stress that the differences in PISN and CCSN observation can be crucial and must be taken into account in order to make more robust statements.  \nUsing the missed PISN observation as a constraint, it is possible to infer an upper limit on the PI/CC ratio, as done in Nicholl et al. 2013. They obtain a PI/CC ratio < 6 × 10 -6 within 𝑧 < 0 . 6 (at 3 𝜎 ), which is compatible with the lowest part of our range. Studies of this kind can be extremely valuable in constraining the parameters of stellar and galaxy evolution models, by excluding all parameter values that produce too high PI/CC ratios. However, these studies would inevitably suffer from important degeneracies. We leave a more accurate assessment regarding these matters to the follow-up work.', '4.4 PISNe in binaries': 'As already outlined, so far we only considered PISNe arising from single stars, under the assumption that interactions in binaries (e.g. mass transfer) do not affect significantly the rate of PISN occurrence. In order to test this assumption, we run the most recent version of the SEVN code (Iorio et al. 2022) to generate a synthetic population of binary stars. We then generate a population of single stars, and compute the 𝑑𝑁 PISN / 𝑑𝑀 SFR for both populations. By assuming that a fraction 𝑓 bin = 0 . 5 of stars lie in binaries, we are able to compute the fraction of PISNe arising in binaries with respect to single stars, 𝑓 PISN bin .  \nWe get 𝑓 PISN bin = 0 . 56, meaning that PISN events are distributed roughly equally among single stars and binaries. Indeed, we find that stellar mergers are the dominant process varying the masses of stars, while other processes such as common envelope and Roche-lobe overflow only play a secondary role. The number of stars entering and exiting the PISN mass range, by merging with another star, roughly balance each other. 𝑓 bin is slightly in favour of binaries due to our bottom-heavy IMF, making it so that stars with lower mass which merge into the PISN range are more abundant than higher mass stars merging out of that range.  \nWe stress that this computation requires an assumption on 𝑓 bin ,  \nwhich is still quite uncertain (e.g. Sana et al. 2012). Nonetheless, we regard 𝑓 bin = 0 . 5 as a fiducial value, based on what established so far in the literature.', '4.5 FMR variation': 'In this work, we adopted the FMR from Curti et al. 2020, one of the most recent determinations (see Section 2.1). As we saw, this prescription leads to a drop in the 𝑍 -dependent SFRD at 𝑧 ≲ 1, at metallicities crucial for PISN production (Figure 3), causing the PISN rate to drop at 𝑧 ≲ 1 (e.g. Figure 8). It is useful to consider an alternative FMR, in order to see how dependent our results are on this relation. We choose to follow Mannucci et al. 2010, who present the following fit:  \n12 + log ( 𝑂 / 𝐻 ) = 8 . 90 + 0 . 37 𝑚 -0 . 14 𝑠 -0 . 19 𝑚 2 + 0 . 12 𝑚𝑠 -0 . 054 𝑠 2 , (10)  \nwhere 𝑚 = log ( 𝑀 ★ ) -10 and 𝑠 = log ( 𝑆𝐹𝑅 ) . Here, log 𝑍 = 12 + log ( 𝑂 / 𝐻 ) -10 . 58, as in Boco et al. 2021. In the following, we will address Curti et al. 2020 as C20, and Mannucci et al. 2010 as M10.  \nFigure 14 shows the 𝑍 -dependent SFRD obtained by adopting the M10 FMR, for 𝜎 Z = [ 0 . 15 , 0 . 35 , 0 . 70 ] and 𝛼 SMF = -1 . 45, to be compared with Figure 3 for the C20 one. We show in Figure 15 the PISN rate we obtain for the same 𝑍 max -𝜎 Z combinations as Figure 8. One can clearly see how this new FMR strongly reduces our range of results, which now sums up to just two/three orders of magnitude, as opposed to the five/seven orders of magnitude of Figure 8. This is due to the fact that, for 𝜎 Z = 0 . 15 and 0.35, the SFRD computed with the M10 FMR extends to lower metallicities, past 𝑍 = 10 -3 , with respect to the one computed with the C20 FMR. This is due to the steeper dependence of galaxy metallicity with stellar mass found in M10 (see also Figure 3 in C20). As a consequence, the PISN rate contribution coming from those metallicities is enhanced, as can be seen in Figure 16. Since this is where 𝑍 max lies for our most pessimistic variations (Figure 5), the corresponding PISN rate is especially enhanced. Remarkably, in the case with 𝑍 max = 1 . 5 × 10 -3 and 𝜎 Z = 0 . 15 (lowest curve in Figures 8 and 15), the difference in rate due to the FMR variation amounts to two/three orders of magnitude, reaching more than four orders of magnitude at 𝑧 = 0. On the other hand, in the case with 𝜎 Z = 0 . 70 the PISN rate does not vary significantly since, as can be seen in the right panels of Figures 3 and 14, the SFRD distribution is quite similar in the two FMRvariations. Moreover, the 𝛼 SMF = 𝛼 SMF ( 𝑧 ) variation produces an even steeper decrease of galaxy metallicity with redshift, with respect to the Curti et al. (2020) FMR case. As a consequence, the PISN rate peak is brought to even higher redshifts, past 𝑧 = 6, as shown in Figure 17.  \nThis warns us about how dependent the PISN rate can be on the adopted FMR. Our results must thus be taken with caution. We stress how it is the lowest end of this range to be affected, while its upper end appears to be immune. It again becomes clear how crucial the interplay between 𝑍 max and 𝜎 Z is in determining the PISN rate, it being the cause of these delicate dependencies.', '4.6 PISNe from Pop III stars': 'It is traditionally believed that only low𝑍 stars explode as PISNe. For this reason, due to their extremely low metallicities 𝑍 ≲ 10 -10  \n(e.g. Cassisi & Castellani 1993), very massive Population III stars are usually considered as main PISN progenitors (e.g. Ober et al. 1983; El Eid et al. 1983; Baraffe et al. 2001; Fryer et al. 2001; Umeda & Nomoto 2001; Heger & Woosley 2002; Scannapieco et al. 2005; Wise & Abel 2005; Langer et al. 2007; Kasen et al. 2011; Pan et al. 2012; Dessart et al. 2012; Yoon, S.-C. et al. 2012; de Souza et al. 2013; de Souza et al. 2014; Whalen et al. 2013, 2014; Smidt et al. 2015; Magg et al. 2016; Venditti et al. 2023; Bovill et al. 2024; Wiggins et al. 2024). On the other hand, our stellar evolution tracks and 𝑀 CO variations allow for PISNe up to ∼ solar metallicity, ∼ 1 . 7 × 10 -2 for variation 𝑂 . Combining with the SFRD, we find that the main PISN rate contribution comes from 𝑍 ∼ 10 -3 -10 -2 , typical of Pop II/I stars. This result identifies Pop II/I stars as prominent PISN candidates. However, we stress that in this work we only consider 𝑍 ≥ 1 × 10 -4 and 𝑧 ≤ 6. In order to study the contribution to the PISN rate coming from Pop III stars, and make a comparison with that from Pop II/I, one would need to extend the treatment to lower metallicities and higher redshifts. This task is hindered by the uncertainties surrounding relevant quantities like the Pop III SFRD and IMF, due to the lack of Pop III observations and the overall challenges of searching those very high redshifts. Nonetheless, we attempt at estimating the PISN rate contribution from Pop III stars, and compare it with our results for Pop II/I.  \nWe compute the Pop III PISN rate analogously as in Equation 1. In particular, we consider the 𝑍 = 10 -11 PARSEC stellar evolution tracks, introduced in Costa et al. 2023 inside the framework of SEVN . For the Pop III SFRD, we consider the determination by Hartwig et al. 2022 (see also Santoliquido et al. 2023), obtained with the semi-analitical model A-SLOTH. Pop III stars are believed to follow a more top-heavy IMF with respect to Pop II/I, that might extend up to 1000 𝑀 ⊙ (e.g. Hirano et al. 2015). This motivates us to adopt a Larson IMF (Larson 1998) of the form ∝ 𝑀 𝜉 𝑒 -𝑀 𝑐 / 𝑀 (where 𝑀 ≡ 𝑀 ZAMS ), defined up to 1000 𝑀 ⊙ , with Kroupa-like slope 𝜉 = -2 . 3, and characteristic mass 𝑀 𝑐 = 30 𝑀 ⊙ . We note that varying the upper mass limit from 1000 to 300 𝑀 ⊙ does not affect the results, since the 𝑍 = 10 -11 tracks we use predict a PISN progenitor mass range of ∼ [ 107 -223 ] 𝑀 ⊙ . Finally, we vary 𝑀 CO criterion for PISN as done in the rest of the work.  \nIn Figure 18 we show the 𝑑𝑁 PISN / 𝑑𝑀 SFR obtained with our Larson IMF, compared to the case of a Kroupa IMF defined up to 1000 𝑀 ⊙ . As one can see, the Kroupa IMF produces a Pop III 𝑑𝑁 PISN / 𝑑𝑀 SFR comparable to the Pop II/I one for the most optimistic stellar variations. On the contrary, the more top-heavy, Larson IMF elevates this quantity by almost an order of magnitude. Indeed, the Larson IMF increases the relative abundance of stars with masses in the PISN range, compared to the Kroupa IMF. Since this result depends only slightly on the chosen 𝑀 CO criterion, we fix it to the fiducial case in computing the PISN rate, that we show in Figure 19 for the Larson IMF. As one can see, the Pop III PISN rate lies around the same values of the most pessimistic results obtained for Pop II/I. Indeed, even though the Hartwig et al. 2022 Pop III SFRD is orders of magnitude lower than the one employed in this work for Pop II/I stars, the Larson IMF produces more PISNe per unit mass compared to the Kroupa IMF (Figure 18). Moreover, the inclusion of the whole SFRD contribution in 𝑍 = 10 -11 , allows to avoid the metallicity cut due to 𝑍 max which, as discussed in depth in the paper, can lower significantly the PISN rate, especially in the lowest 𝜎 Z cases.  \nPrevious Pop III PISN rate determinations in the literature tend  \nFigure 14. As Figure 3, for the Mannucci et al. 2010 FMR.  \n<!-- image -->  \nFigure 15. As Figure 8, for the Mannucci et al. 2010 FMR.  \n<!-- image -->  \nto distribute in the middle/upper range of our results for Pop II/I stars (Figure 8), with some even reaching our most optimistic results (Scannapieco et al. 2005; Wise & Abel 2005; Pan et al. 2012; de Souza et al. 2013; de Souza et al. 2014; Magg et al. 2016; Regős et al. 2020; Tanikawa et al. 2023; Venditti et al. 2023; Wiggins et al. 2024; Tanikawa 2024). All in all, despite the uncertainties on the Pop III SFRD and IMF, these considerations suggest that Pop III could somewhat contribute to the PISN rate at 𝑧 ≳ 6.', '5 CONCLUSIONS': 'In this work, we compute the PISN rate as a function of redshift up to 𝑧 = 6, by combining stellar evolution tracks with a semi-empirical determination of the metallicity-dependent cosmic star formation history. The aim is to study how the uncertainties in both stellar and galaxy evolution theory affect our knowledge about the rate of PISN occurrence throughout cosmic history. We vary stellar evolution code, 𝑀 CO criterion to have PISN, and IMF upper limit, as well as the dispersion of the galaxy metallicity distribution and the low-mass end slope of the GSMFs. We find these variations to have  \na huge effect on our results. The PISN rate spans roughly three orders of magnitude under stellar variations. Their effect can be boiled down to the dependence on the maximum metallicity to have PISN, 𝑍 max, and the height of the 𝑑𝑁 PISN / 𝑑𝑀 SFR curve. Remarkably, the interplay between 𝑍 max and 𝜎 Z can extend the PISN rate range up to ∼ seven orders of magnitude, depending on redshift. The local, 𝑧 = 0 rate ranges from 4 × 10 -4 to 2 × 10 3 𝑦𝑟 -1 𝐺𝑝𝑐 -3 , while the values at peak range from ∼ 10 -1 to 3 × 10 4 𝑦𝑟 -1 𝐺𝑝𝑐 -3 . Moreover, prescribing a redshift dependence of 𝛼 SMF increases the PISN rate at 𝑧 > 3 -4. Also the position of the peak gets affected, shifting from 𝑧 ≲ 2 up to 𝑧 ∼ 4 -5 throughout all variations. This warns us about the delicate link between these two metallicity parameters, in choosing stellar evolution prescriptions and IMF leading to a certain 𝑍 max, and in adopting recipes for the cosmic metallicity distribution of galaxies.  \nAll in all, our limited knowledge makes the PISN rate very uncertain. On the other hand, the strong dependencies we find offer the chance to constrain these parameters based on possible PISN observations in the future, or the lack of PISN observations in the eventuality that these transients are never discovered. Even accounting for the significant degeneracies, this would help dispel the uncertainties around massive star evolution, such as the criterion for a star to undergo pair instability and the maximum stellar mass, and the evolution of star formation and galaxy metallicity throughout cosmic history. We will delve into these aspects in a follow-up work, dedicated to studying PISN detection with telescopes such as JWST, in which we will attempt at answering the question why PISN explosions have never been observed.  \nOur fiducial model indicates galaxies with stellar masses between ∼ 10 9 and 10 10 𝑀 ⊙ , and metallicities ≲ 10 -2 , to be the favourable PISN hosts, providing the main contribution to the PISN rate. These values shift to 𝑀 ★ ∼ 10 8 -10 9 𝑀 ⊙ and 𝑍 ≳ 10 -3 in our most pessimistic variation, and to 𝑀 ★ ∼ 10 10 𝑀 ⊙ and 𝑍 ∼ 10 -2 in our most optimistic one. This can represent useful information for the design of current and future campaigns aimed at observing these elusive transients, with facilities such as JWST, the Vera Rubin Observatory , Euclid , the Nancy Grace Roman Space Telescope , the Zwicky Transient Facility , and ULTIMATE-Subaru (e.g. Weinmann & Lilly 2005; Whalen et al. 2012; Smidt et al. 2015; Kozyreva et al. 2014a; Hartwig et al. 2018; Regős et al. 2020; Moriya et al. 2019, 2022a,b; Tanikawa et al. 2023; Tanikawa 2024).  \nFigure 16. As Figure 7, for the Mannucci et al. 2010 FMR.  \n<!-- image -->  \nFigure 17. As in Figure 9, for the Mannucci et al. 2010 FMR.  \n<!-- image -->  \nFinally, our stellar evolution tracks and 𝑀 CO variations allow stars to explode as PISNe up to ∼ solar metallicity. As a consequence, we find peak metallicities for the PISN rate that suggest Pop II/I stars to be prominent PISN progenitors. This goes against the traditional expectation that PISNe come exclusively or preferentially from verylow-metallicity, Pop III stars, a point which is worth exploring in the future.  \nFigure 18. Number of PISNe per unit star forming mass produced by Pop III stars with 𝑍 = 10 -11 , adopting a Larson IMF up to 1000 𝑀 ⊙ (filled markers). We also show the case of a Kroupa IMF defined up to 1000 𝑀 ⊙ (empty markers). We indicate the 𝑀 CO variations with different shapes. For comparison, on the right is shown the same plot as in Figure 5.  \n<!-- image -->  \nZ', 'ACKNOWLEDGEMENTS': 'The authors thank the anonymous referee for constructive comments, and Max Briel and David Hendriks for helpful discussions.  \nFigure 19. PISN rate from Pop III stars as a function of redshift, computed by adopting a Larson IMF and the Pop III SFRD determination by Hartwig et al. 2022. The red band indicates the uncertainty coming from the 95% credible interval of the Hartwig et al. 2022 SFRD, as computed in Santoliquido et al. 2023. For clarity, we fix the fiducial, 𝑀 CO = [ 55 -110 ] 𝑀 ⊙ variation. On the left we report the results obtained for Pop II/I stars, as in Figure 8, for comparison.  \n<!-- image -->  \nThis work was partially funded from the projects: \'Data Science methods for MultiMessenger Astrophysics & Multi-Survey Cosmology\' funded by the Italian Ministry of University and Research, Programmazione triennale 2021/2023 (DM n.2503 dd. 9 December 2019), Programma Congiunto Scuole; EU H2020-MSCA-ITN-2019 n. 860744 BiD4BESt: Big Data applications for black hole Evolution STudies ; Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations); European Union - NextGenerationEU under the PRIN MUR 2022 project n. 20224JR28W "Charting unexplored avenues in Dark Matter"; INAF Large Grant 2022 funding scheme with the project "MeerKAT and LOFAR Team up: a Unique Radio Window on Galaxy/AGN coEvolution; INAF GO-GTO Normal 2023 funding scheme with the project "Serendipitous H-ATLAS-fields Observations of Radio Extragalactic Sources (SHORES)". The authors declare no conflict of interest.', 'APPENDIX A: HE CORE MASS CRITERION': 'In this Section, we explore how adopting the criterion on 𝑀 He , instead of 𝑀 CO , affects our results. We impose that only stars with 𝑀 He in the range [64-135] 𝑀 ⊙ explode as PISNe, as adopted in SEVN (Spera & Mapelli 2017). We consider [70-120] and [60-140] 𝑀 ⊙ as pessimistic and optimistic ranges respectively, based on previous works (see references in Section 2.2). Table A1 shows the stellar variations we considered, namely the same as in Table 2 but with 𝑀 He ranges. Note that 𝑍 max is modified as a consequence.  \nAs shown in Figure A1, the shape of 𝑑𝑁 PISN / 𝑑𝑀 SFR is similar to the 𝑀 CO case (Figure 5), except for variation 𝑀 3, which exhibits a harsh drop at metallicities between ∼ 4 and 9 × 10 -3 . This is due to envelope undershooting, reducing the mass of the core and thus preventing stars in the corresponding mass range from going into  \nTable A1. As in Table 2, with 𝑀 He mass ranges and, consequently, different 𝑍 max.  \nFigure A1. As in Figure 5, using the 𝑀 He criterion (see Table A1).  \n<!-- image -->  \nZ  \nPISN. The effect on 𝑀 He is stronger than on 𝑀 CO , where it manifests as a softer drop. See Costa et al. 2021 for a detailed explanation of this process in PARSEC . We note that this drop is not particularly relevant for the PISN rate, the main contribution coming from the bulge between 𝑍 ∼ 2 × 10 -3 and 7 × 10 -2 . 𝑍 max is affected in a nonnegligible way in most variations, as shown in Table A1. The values of 𝑑𝑁 PISN / 𝑑𝑀 SFR change by some small factor. Note that we chose core mass ranges which are not correspondent to each other, between the 𝑀 CO and 𝑀 He cases, which is mostly causing these differences.  \nIn Figure A2, we show the PISN rate resulting from each stellar variation. As one can see, the rates only change by some small factor, in any case less than one order of magnitude, mainly due to the different 𝑍 max (compare with Figure 6). In particular, variations 𝑃 and 𝑂 are not subject to dramatic changes, so that the range of results is left substantially unaffected.  \nAnalogously, we find that also applying our galactic variations does not lead to significant differences, with the PISN rate ranges being only slightly shrinked with respect to the 𝑀 CO case, as shown in Figures A3 and A4.  \nOverall, we find that varying between 𝑀 CO and 𝑀 He criterion does not change our results significantly.  \nFigure A2. As in Figure 6, using the 𝑀 He criterion.  \n<!-- image -->  \n<!-- image -->  \nz  \nFigure A3. As in Figure 8, using the 𝑀 He criterion.Figure A4. As in Figure 9, using the 𝑀 He criterion.  \n<!-- image -->  \nFigure B1. As in Figure 5, for 𝑀 up = 300 and 600 𝑀 ⊙ (solid and dotted lines respectively), and for 𝑀 CO ∈ [ 60 -105 ] and [ 45 , 120 ] 𝑀 ⊙ (blue and orange lines). The PARSEC-II stellar evolution tracks have been employed.  \n<!-- image -->', 'APPENDIX B: ADDITIONAL IMF UPPER LIMIT VARIATION': 'Throughout this work, for simplicity we assumed a Kroupa IMF, and explored variations on its upper limit, 𝑀 up = 150 and 300 𝑀 ⊙ . Extending this IMF to higher masses does not affect our results significantly. In order to prove this point, we consider an additional variation with 𝑀 up = 600 𝑀 ⊙ .  \nIndeed, considering the typical 𝑀 ZAMS of ∼ [ 140 , 260 ] 𝑀 ⊙ for PISN progenitors, one would not expect IMFs with 𝑀 up > 260 𝑀 ⊙ to have any effect. However, the stellar evolution codes we consider here produce 𝑀 ZAMS ranges varying significantly with metallicity, as shown in Table 1. In particular, they predict even VMSs beyond 300 𝑀 ⊙ to end their life as PISN, at high metallicity. Therefore, in principle extending the IMF to 600 𝑀 ⊙ could play a role.  \nAmong our set of stellar variations, we select 𝐹 and 𝑂 , since they have respectively the lowest and highest 𝑍 max among all possible variations with 𝑀 up = 300 𝑀 ⊙ . For simplicity, we do not consider different variations on 𝑀 CO , which do not play a significant role in this case. Figure B1 shows the 𝑑𝑁 PISN / 𝑑𝑀 SFR we obtain by imposing 𝑀 up = 300 and 600 𝑀 ⊙ . In variation 𝑂 , setting 𝑀 up to 600 𝑀 ⊙ has an appreciable effect only at 𝑍 ≳ 10 -2 , where it slightly increases 𝑑𝑁 PISN / 𝑑𝑀 SFR . In variation 𝐹 , 𝑑𝑁 PISN / 𝑑𝑀 SFR gets increased starting from 𝑍 ≲ 10 -3 , and also 𝑍 max is slightly higher. As one can see in Figure B2, the effects on the PISN rate are negligible for variation 𝑂 , while in variation 𝐹 the rate is increased only by a small factor, between 1 and 2. This is because in the first case 𝑍 max is greater than the SFRD peak (see Section 2.1), therefore the main contribution to the PISN rate is already included. In variation 𝐹 , 𝑍 max is lower than the SFRD peak, so the rate is more sensible to this variation, even though it is affected only in a minor way.  \nAll in all, we find that varying 𝑀 up from 300 to 600 𝑀 ⊙ does not have a significant effect on the PISN rate. This is due to the fact that our Kroupa IMF prescribes a low number of stars with mass between 300 and 600 𝑀 ⊙ . Adopting a more top-heavy IMF might lead to different results.  \nFigure B2. PISN rate as a function of redshift, obtained with the PARSEC-II tracks for the stellar variations shown in Figure B1. Galactic parameters are fixed to 𝜎 Z = 0 . 35 and 𝛼 SMF = -1 . 45.  \n<!-- image -->', 'DATA AVAILABILITY': 'The data used in this work can be made available upon request.', 'REFERENCES': "Abbot e. a., 2016, Physical Review Letters, 116  \nAbbot e. a., 2020, Physical Review Letters, 125  \n- Aguado D. S., et al., 2023, Monthly Notices of the Royal Astronomical Society, 520, 866\n- Aoki W., Tominaga N., Beers T. C., Honda S., Lee Y. S., 2014, Science, 345, 912\n- Artale M. C., Mapelli M., Giacobbo N., Sabha N. B., Spera M., Santoliquido F., Bressan A., 2019, Monthly Notices of the Royal Astronomical Society, 487, 1675  \nAsplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481  \nBallero S. K., Matteucci F., Chiappini C., 2006, New Astronomy, 11, 306  \n- Baraffe I., Heger A., Woosley S. E., 2001, The Astrophysical Journal, 550, 890\n- Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379\n- Belczynski K., 2020, The Astrophysical Journal, 905, L15\n- Belczynski K., Holz D. E., Bulik T., O'Shaughnessy R., 2016a, Nature, 534, 512\n- Belczynski K., et al., 2016b, Astronomy & Astrophysics, 594, A97 Bestenlehner J. M., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 1918\n- Bisnovatyi-Kogan G. S., Kazhdan Y. M., 1967, Soviet Astronomy, 10, 604 Boco L., Lapi A., Goswami S., Perrotta F., Baccigalupi C., Danese L., 2019, The Astrophysical Journal, 881, 157\n- Boco L., Lapi A., Chruslinska M., Donevski D., Sicilia A., Danese L., 2021, The Astrophysical Journal, 907, 110  \nBotticella M. T., et al., 2007, Astronomy andamp; Astrophysics, 479, 49-66  \nBotticella M. T., Smartt S. J., Kennicutt R. C., Cappellaro E., Sereno M., Lee  \nJ. C., 2012, Astronomy andamp; Astrophysics, 537, A132  \n- Bouwens R. J., et al., 2021, The Astronomical Journal, 162, 47\n- Bovill M. S., Stiavelli M., Wiggins A. I., Ricotti M., Trenti M., 2024, The Astrophysical Journal, 962, 49\n- Brands S. A., et al., 2022, Astronomy and Astrophysics, 663, A36  \nBressan A., Marigo P., Girardi L., Salasnich B., Cero C. D., Rubele S., Nanni  \nA., 2012, Monthly Notices of the Royal Astronomical Society, 427, 127  \n- Briel M. M., Eldridge J. J., Stanway E. R., Stevance H. F., Chrimes A. A., 2022, Monthly Notices of the Royal Astronomical Society, 514, 1315  \n- Whalen D. J., et al., 2013, The Astrophysical Journal, 777, 110\n- Whalen D. J., Smidt J., Even W., Woosley S. E., Heger A., Stiavelli M., Fryer C. L., 2014, The Astrophysical Journal, 781, 106\n- Wiggins A. I., Bovill M. S., Strolger L.-G., Stiavelli M., Bowling C., 2024, Kindling the First Stars II: Dependence of the Predicted PISN Rate on the Pop III Initial Mass Function ( arXiv:2402.17076 )\n- Winch E. R. J., Vink J. S., Higgins E. R., Sabhahit G. N., 2024, Monthly Notices of the Royal Astronomical Society, p. stae393\n- Wise J. H., Abel T., 2005, The Astrophysical Journal, 629, 615\n- Woosley S. E., 2017, The Astrophysical Journal, 836, 244\n- Woosley S. E., Heger A., 2021, The Astrophysical Journal Letters, 912, L31\n- Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390 Xing Q.-F., et al., 2023, Nature, 618, 712\n- Yaron O., Gal-Yam A., 2012, Publications of the Astronomical Society of the Pacific, 124, 668\n- Yoon, S.-C. Dierks, A. Langer, N. 2012, A&A, 542, A113\n- Ziegler J. J., Edwards T. D. P., Suliga A. M., Tamborra I., Horiuchi S., Ando S., Freese K., 2022, Monthly Notices of the Royal Astronomical Society, 517, 2471-2484\n- de Jager C., Nieuwenhuĳzen H., van der Hucht K. A., 1988, A&AS, 72, 259 de Souza R. S., Ishida E. E. O., Johnson J. L., Whalen D. J., Mesinger A., 2013, MNRAS, 436, 1555\n- de Souza R. S., Ishida E. E. O., Whalen D. J., Johnson J. L., Ferrara A., 2014, Monthly Notices of the Royal Astronomical Society, 442, 1640-1655\n- du Buisson L., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5941\n- van Loon J. T., Cioni M. R. L., Zĳlstra A. A., Loup C., 2005, A&A, 438, 273  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author."}