The first fireworks: A roadmap to Population III stars during the Epoch of Reionization through Pair-Instability Supernovae

With the launch of JWST and other scheduled missions aimed at probing the distant Universe, we are entering a new promising era for high-z astronomy. One of our main goals is the detection of the first population of stars (Population III or Pop III stars), and models suggest that Pop III star formation is allowed well into the Epoch of Reion-ization (EoR), rendering this an attainable achievement. In this paper, we focus on our chance of detecting massive Pop IIIs at the moment of their death as Pair-Instability Supernovae (PISNe). We estimate the probability of discovering PISNe during the EoR in galaxies with different stellar masses ( 7 . 5 ⩽ Log( M ⋆ / M ⊙ ) ⩽ 10 . 5 ) from six dustyGadget simulations of 50 h − 1 cMpc per side. We further assess the expected number of PISNe in surveys with JWST/NIRCam and Roman/WFI. On average, less than one PISN is expected in all examined JWST fields, while ≃ 1 . 5 η III PISNe may be found in a ∼ 1 deg 2 Roman field, with potential for increased discoveries when considering more optimistic choices for the Pop III initial mass function and formation efficiency η III , and when including the contribution of coarsely-resolved environments. JWST/NIRCam and Roman/WFI allow the detection of massive-progenitor ( ∼ 250 M ⊙ ) PISNe throughout all the optimal F200W-F356W, F277W-F444W, and F158-F213 colors. PISNe are also predominantly located at the outskirts of their hosting haloes, facilitating the disentangling of underlying stellar emission thanks to the spatial-resolution capabilities of the instruments.


INTRODUCTION
Pair-instability supernovae (PISNe) are a type of supernova (SN) that occurs when a massive star, with a mass between 140 M⊙ and 260 M⊙, undergoes a thermonuclear explosion due to the production and subsequent annihilation of electron-positron pairs in its hot core.These reactions lead to a rapid loss of radiation pressure support causing the core to collapse, and thus igniting explosive oxygen and silicon ⋆ E-mail:alessandra.venditti@inaf.itburning that ultimately results in the complete disruption of the star (Rakavy & Shaviv 1967;Barkat et al. 1967;Fraley 1968;Bond et al. 1984;Fryer et al. 2001).PISNe are particularly interesting because they provide valuable insights on the first generation of stars, known as Population III (Pop III) stars.
Pop III stars formed from pristine gas, almost entirely composed of hydrogen and helium.Because of the absence of heavy elements, they are expected to have very different properties compared to second-generation (Pop II) stars, that formed from the gas pre-enriched by Pop III.Most no-tably, they are expected to be predominantly massive (Abel et al. 2002;Bromm et al. 2002), with masses ranging from ∼10s M⊙ to ∼100s M⊙ (Hosokawa et al. 2011;Hirano et al. 2014;Stacy et al. 2016), or even ∼1000s M⊙ (Hirano et al. 2015a,b;Susa et al. 2014;Hosokawa et al. 2016;Sugimura et al. 2020;Latif et al. 2022), implying that some of the most massive Pop III stars may be able to produce PISNe.When a PISN explodes, large amounts of metals are released into the surrounding medium (Heger & Woosley 2002).Hence, PISNe arising from Pop III stars are thought to be one of the main sources of heavy elements in the early Universe, and are therefore key to understand the process of early chemical enrichment (e.g.Salvadori et al. 2007Salvadori et al. , 2008;;Karlsson et al. 2013;de Bennassuti et al. 2017;Aguado et al. 2023a).
Although a number of potential Pop III systems at z > 6 have been recently proposed in the literature (Vanzella et al. 2020(Vanzella et al. , 2023;;Wang et al. 2022;Maiolino et al. 2023), an unambiguous detection of a Pop III-hosting galaxy is currently still missing.The identification of Pop III stars in high-z galaxies through different spectroscopic line diagnostics has been extensively discussed over the last decades (Bromm et al. 2001a;Inoue 2011;Zackrisson et al. 2011;Mas-Ribas et al. 2016;Nakajima & Maiolino 2022;Trussler et al. 2023;Katz et al. 2023;Cleri et al. 2023).However, the properties of PISNe, such as their explosive energy output and light curve, can also provide important constraints on the characteristics of their progenitor stars (Mackey et al. 2003;Scannapieco et al. 2005;Whalen et al. 2013;Pan & Loeb 2013;de Souza et al. 2014;Wang et al. 2017), and thus shed light on the properties of Pop III stars.Several simulations have been carried out to model the distinctive light curve and afterglow emission of PISNe arising from Pop IIIs (Scannapieco et al. 2005;Woosley et al. 2007;Pan et al. 2012;Chen et al. 2014;Whalen et al. 2014;Jerkstrand et al. 2016;Kozyreva et al. 2017;Gilmer et al. 2017;Hartwig et al. 2018).Simulations suggest that they are very energetic, with total energy produced up to 10 53 erg, and have an extended light curve duration, ∼ 1 yr in the source frame (Kasen et al. 2011).
The detection of PISNe is especially challenging due to their rarity.Pop III stars are expected to have formed at high redshifts (z ∼ 20 − 30; Bromm 2013; Klessen & Glover 2023), when the Universe was only a few hundred million years old, and their rate of formation is expected to have declined rapidly as the Universe became more metal-enriched -although recent simulations suggest a late Pop III star formation, down to the Epoch of Reionization (EoR, z ∼ 6), is also possible (Xu et al. 2016;Jaacks et al. 2019;Liu & Bromm 2020;Sarmento et al. 2018;Sarmento & Scannapieco 2022;Skinner & Wise 2020;Visbal et al. 2020;Venditti et al. 2023).The exact rate of PISNe across cosmic time is however highly uncertain (Miralda-Escudé & Rees 1997;Wise & Abel 2005;Hummel et al. 2012;Pan et al. 2012;Johnson et al. 2013;Tanaka et al. 2013;Magg et al. 2016), as it depends on the details of the metal enrichment process, as well as other factors such as the Pop III initial mass function (IMF).A more top-heavy Pop III IMF, for example, would yield a higher rate of PISNe, as a large fraction of the stars would be massive enough to undergo pair instability (Lazar & Bromm 2022).However, the exact shape of the Pop III IMF is unconstrained.The typical metallicity of the gas in which Pop III stars formed is also uncertain, as it depends on the details of the primordial star-formation process and the efficiency of metal mixing.Simulations suggest that Pop III stars formed in very low-metallicity environments, with metallicity not higher than 10 −4 − 10 −6 Z⊙ (Bromm et al. 2001a;Omukai et al. 2005;Maio et al. 2010;Schneider et al. 2012a,b;Chiaki et al. 2016;Chiaki & Yoshida 2022), although a massive stellar component in excess of what is predicted by a present-day IMF may persist up to metallicity ≲ 10 −2 Z⊙ (Chon et al. 2021(Chon et al. , 2022)).
Despite the challenges, several efforts have been made to detect PISNe in the past.One of the most promising candidates proposed is SN 2007bi, a super-luminous SN (SLSN), discovered in 2007 in a nearby dwarf galaxy (Gal-Yam et al. 2009).SN 2007bi has several properties that are consistent with a PISN, including its high energy output, long duration, measured yield of radioactive nickel and elemental composition of the ejecta.However, its exact nature is still a matter of debate, and other explanations have been proposed in the literature (see e.g.Moriya et al. 2010;Dessart et al. 2012).Other examples of SLSNe with slowly evolving light curves at high redshift (z = 2.05 and z = 3.90) have been found with Keck I follow-up, late-time spectroscopy of two objects first identified in archival data from the Canada-France-Hawaii Telescope Legacy Survey Deep Fields (Cooke et al. 2012).Possible survey strategies dedicated to the search of PISNe during the EoR with JWST and with the Nancy Grace Roman Space Telescope have been discussed in Hartwig et al. (2018), Regős et al. (2020) and Moriya et al. (2022b), while the possibility of detecting PISNe in the Euclid Deep Survey (EDS, Laureijs et al. 2011;Euclid Collaboration et al. 2022) at z ≲ 2.5 is studied in Moriya et al. (2022a).
Stellar archaeology can provide a complementary probe to constrain the properties of the first stars -particularly their IMF -and their SN outcomes, by looking at the atmospheres of old and metal-poor stars in our Galaxy.Notably, the detection of a peculiar chemical abundance pattern in a star with [Fe/H] = −2.5, possibly indicative of a PISN nucleosynthetic yield (Aoki et al. 2014), presents evidence that very high masses are indeed allowed for Pop III stars (de Bennassuti et al. 2017).Recently, Xing et al. (2023) have reported the identification of a halo star, with [Fe/H] = −2.42,whose abundance pattern is remarkably well fit by a PISN imprint from a 260 M⊙ progenitor, although a core-collapse SN-origin is at least equally plausible with the currently measured elements (Jeena et al. 2023).Other instances of metal-poor stars with suggested PISN signatures in their atmosphere have been discussed in Salvadori et al. (2019) and Aguado et al. (2023a), and a contribution of PISNe from massive Pop III stars was also suggested to explain the unusual abundance feature of Fe and Mg in the broad-line region of a quasar at z = 7.54 (Yoshii et al. 2022).Simulations of early chemical enrichment (e.g.Jeon et al. 2015) have suggested that PISN-produced material predominantly resides in the intergalactic medium (IGM), raising the possibility of using the ultra-luminous afterglows of high-redshift gamma ray bursts (GRBs) to probe this IGM signature with deep absorption spectroscopy (Wang et al. 2012).
An important factor in the detectability of PISNe is their host environment.Pop III stars are expected to have formed predominantly in low-mass, faint haloes.However, some models also suggest that Pop III could still form in more massive and luminous haloes (Liu & Bromm 2020;Bennett & Sijacki 2020;Riaz et al. 2022;Venditti et al. 2023).Although massive galaxies are generally easier to detect, their higher luminosity might even hinder our ability to clearly identify a PISN in such hosts.This is especially true if the SN is observed at late times, when the emission has faded considerably (Kasen et al. 2011).An extensive study of the expected rate of PISNe as a function of their host environments is currently missing.
This work is aimed at assessing the probability of observing PISNe arising from Pop III stars during the EoR, in galaxies with different stellar masses.We rely on a suite of six 50h −1 cMpc simulations already introduced in Di Cesare et al. ( 2023) and Venditti et al. (2023), carried out with the hydrodynamical code dustyGadget (Graziani et al. 2020).Our findings can guide the design of optimal strategies to search for potential PISN-hosting candidates, both from archival data and from dedicated surveys.While achieving the first clear direct detection of a PISN event would be a remarkable accomplishment in and of itself, finding PISNe at high redshifts would also have deeper implications for the search of the first stars.As PISNe are in fact the natural outcome of very massive stars, hardly attainable for normal Pop II/I stars, they can be used as markers to identify potential Pop III hosts.Hence, the synergy between wide-field PISNe searches and spatially resolved observations of their host galaxies (possibly enhanced by gravitational lensing) may offer one of our finest tools to directly probe Pop III stars and test our predictions on primordial star-formation.
The paper is organised as follows.Section 2 presents our methodology: Section 2.1 describes our simulation suite, while Section 2.2 outlines our strategy for evaluating the probability of finding PISNe during the EoR from the simulations.Our results are presented in Section 3. Specifically, in Section 3.1 we compute the expected number of PISNe in galaxies of different stellar mass per unit halo and volume, and we predict the expected number of PISNe within JWST and Roman fields.In Section 3.2 we discuss our capability of observing PISNe and discriminating from the stellar emission of their host galaxies, by considering their integrated emission (Section 3.2.1)and possibly resolved observations (Section 3.2.2).A critical discussion is presented in Section 4, including a number of caveats in our modelling (Section 4.1), and a discussion on the implication of either future detections or non-detections of PISNe (Section 4.2), and on the advantages and disadvantages of many possible detection strategies (Section 4.3).Finally, we summarize and offer conclusions in Section 5.

Simulating the cosmological environment
Our study is based on a suite of eight cosmological simulations carried out with the hydrodynamical code dustyGadget.The code is extensively described in Graziani et al. (2020), especially regarding its most novel feature, i.e. the implementation of a self-consistent model for dust production and evolution, while our simulations (hereafter named as U6 -U13) are introduced in Di Cesare et al. (2023) and Venditti et al. (2023).More specifically, Di Cesare et al. (2023) studied the main scaling relations in the context of available observations and model predictions, while Venditti et al. (2023) focused on the properties of galaxies hosting Pop III star formation during the EoR, which are of particular interest for the present work.We briefly describe relevant aspects of the simulations here, and refer to the aforementioned papers for further details.
Pop II/I stars with masses ⩾ 40 M⊙ and Pop III stars outside the PISN mass range do not contribute to the metal enrichment, as they are assumed to directly collapse into black holes (not followed explicitly here).The impact of our choice of the Pop III IMF, and hence of neglecting the contribution of traditional core-collapse SNe from Pop IIIs, is discussed in Section 4.
The gas chemical evolution model is adopted from Tornatore et al. (2007b).Dust and metals are spread in the ISM through a spline kernel, and galactic winds are also modelled as in Springel & Hernquist (2003) with a constant velocity of 500 km s −1 , as indicated by observations of normal galaxies in the ALMA Large Program to INvestigate [CII] at Early times (ALPINE) survey (Ginolfi et al. 2020).
The identification of dark matter haloes and their substructures is performed in post-processing with the AMIGA halo finder (Knollmann & Knebe 2009).Finally, the simulations are carried out from z ≃ 100 down to z ≃ 4, while this work only focuses on the redshift range 6 ≲ z ≲ 8. Indeed, the number of well-resolved Pop III-hosting haloes (with a stellar mass log(M⋆/M⊙) ≳ 7.5, corresponding to a number of stellar particles ≳ 20) decreases significantly at higher redshifts, while a drop in Pop III star formation is expected at lower redshifts due to the joint effect of cosmic metal enrichment and ultraviolet (UV)/Lyman-Werner (LW) radiative feedback.A more extensive discussion can be found in Venditti et al. (2023).
The simulations described in Di Cesare et al. ( 2023) and Venditti et al. (2023) have been chosen for the present study because, given the rarity of PISNe, having a big simulated volume is key to capture the statistics of these events.The convergence of our model has been tested in Graziani et al. (2020) 2 .Given our limited mass resolution, we focus on a much higher mass regime for Pop III-hosting halos (log(M⋆/M⊙) ≳ 7.5) than that of the first mini-halos.In a future work, zoom-in simulations will be performed to investigate the contribution of low-mass halos to Pop III star formation across cosmic times, including the very first starforming regions.

Computing the probability of finding PISNe
The average number of PISNe, N PISN, produced by a Pop III stellar population per unit Pop III mass MIII, with our assumed IMF ϕ(m) (whose lower and upper limits are m low and mup, respectively), is given by: indicating an average of one PISN event per 460 M⊙.For Pop III stellar populations with mass ∼ 2 × 10 6 M⊙, this would result in more than 4000 PISNe produced on average by each stellar population (but see the discussion below on more realistic outcomes).We note that N PISN/MIII can vary depending on the adopted IMF.As both the shape and mass range of the Pop III IMF are still largely uncertain, in Table 1 we evaluate this quantity for a compilation of IMFs that have been adopted in the literature to describe Pop III stellar populations (see Figure 1).The shape of all these IMFs can be expressed in the general form (Lazar & Bromm 2022): i.e. a modified-Larson IMF (Larson 1998), with various parameters describing the slope (α) and the exponential cut-off at low stellar masses (mcut, β).We note that N PISN/MIII can vary in the range ∼ [0.2 − 4] × 10 −3 M −1 ⊙ over all the included IMFs.
Another source of uncertainty comes from our Pop III stellar mass resolution element MIII,res ∼ 2 × 10 6 M⊙.Indeed, since the expected number of PISNe is directly proportional to the total Pop III mass MIII, it is important to carefully estimate the amount of stellar mass produced in a single star-formation event.In principle, this could be lower than MIII,res, which should more accurately be interpreted  2018) (pink ), and they can all be expressed through the general expression at the bottom (Equation 2), with parameters specified in Table 1.Our Salpeterlike IMF (Graziani et al. 2020;Di Cesare et al. 2023;Venditti et al. 2023) is shown as a black, solid line.We also include IMFs with a broader mass range, extending up to 1000 M ⊙ , as allowed from ab-initio simulations of star-forming clouds (e.g.Susa et al. 2014;Hirano et al. 2015a;Hosokawa et al. 2016, grey).The gray shaded region indicates the stellar progenitor mass range of PISNe.
Table 1.Values of the average number of PISNe produced per unit Pop III mass N PISN /M III computed for a compilation of Pop III IMFs, expressed through the general form of Equation 2, with parameters α, β, mcut and mass range mrange as specified.The IMF adopted in this work is highlighted in bold.See Figure 1 for references.
The mass range for these IMFs has been extended up to 260 M ⊙ with respect to the original upper limit (150 M ⊙ ) of Jaacks et al. (2018), to account for the proposed constraint on the Pop III IMF upper limit from Xing et al. (2023).We note that higher values for N PISN /M III (in the range 9.32 − 12.13) are found when this constraint is removed.
as the amount of extremely metal-poor gas above our density threshold that is available for star formation.We can formally express this through an efficiency factor ηIII < 1: and place a lower limit on ηIII ∼ 0.01 from a wealth of simulation data describing the first sites of Pop III star formation in mini-haloes (see Bromm 2013 and references therein).
However, mini-haloes at the Rees-Ostriker cooling threshold (Rees & Ostriker 1977) are known to be very inefficient hosts for star formation, given their shallow potential well.We can thus argue that ηIII could be larger in more massive haloes at later times, and explore our results for a range of different values of ηIII ⩾ 0.01.As no available simulation constrains the mass regime we are currently exploring, we allow this parameter to vary up to about an order of magnitude more than our lower limit ηIII = 0.01.An empirical hint that higher values of ηIII (even up to ηIII ≃ 0.3) might be allowed at least in some cases is offered by the tentative detection of HeII emission at ≃ 2.5 kpc from an exceptionally luminous galaxy at z = 10.6, that might indicate the presence of a Pop III cluster with a top-heavy IMF and a total Pop III mass of 6 − 7 × 105 M⊙ (Maiolino et al. 2023).With this caveat in mind, the average number of PISNe produced by a single Pop III stellar population is reduced to ∼ 4000 ηIII, for our adopted IMF.We emphasize that, although both the Pop III IMF and star-formation efficiency are still very uncertain, they only affect the expected number of PISNe through the scaling parameters N PISN/MIII and ηIII.Hence, our results can be easily adjusted if further constraints on either quantities will become available.In fact, our goal is to provide a framework to estimate the probability of achieving a direct detection of PISNe and using these PISN discoveries to tag sites of late Pop III star formation close to reionization, rather than compute an exact rate.In the following, the case N PISN/MIII = 2.17 × 10 −3 M −1 ⊙ (corresponding to our Salpeter-like IMF in the range [100, 500] M⊙) and ηIII = 0.1 will be referred to as our "reference model", as baseline for our discussion.We emphasize that this choice of reference values is for convenience of presentation only, and does not imply that they are the most realistic ones.
Since Pop III stellar particles in our simulations are formed in an instantaneous burst, and Pop IIIs with initial masses m in the range [140,260] M⊙ explode as PISNe over a time span ∆tPISN = τ (140 M⊙) − τ (260 M⊙), where τ (m) is the lifetime of a star with initial mass m, we can approximately compute the average number of PISNe expected to be seen at a given time t within ∆tPISN as: where ∆tprompt ≃ 1 yr is the rest-frame time interval over which the prompt PISN emission is expected to be bright (Kasen et al. 2011).
The Pop III lifetimes τ (m) as a function of initial stellar mass m are taken from the models of Schaerer (2002) 3 .Given the uncertainties for existing Pop III evolution models, they reported results for models assuming strong mass loss4 arising from high-mass Pop III stars and models assuming no mass loss at all (see their tables 5 and 4 respectively).Schaerer (2002), in the mass range of our assumed IMF.The black, solid line refers to the model assuming no mass loss (see table 4 of the original paper), while the black, dotted line refers to the model assuming strong mass loss (see their table 5).The interpolated lifetimes at the limits of the PISN progenitors mass range [140, 260] M ⊙ (grey, shaded area) are indicated for the two models.The average lifetimes τ of Pop III stellar populations for the two models with our assumed IMF are also shown, on top of the horizontal, solid/dotted, black lines.
In Figure 2, we show τ (m) for both models over the mass range of our assumed Pop III IMF, to emphasize the existing differences, indicating that mass loss can result in longer lifetimes by up to a factor of ≃ 1.3, and even in an inversion of the trend for masses m ≳ 300 M⊙, with τ increasing rather than decreasing for larger stellar masses.The lifetimes of stars with masses 140 M⊙ and 260 M⊙ are also shown in the plots.With these lifetimes, ∆tPISN ≃ 0.65/0.33Myr respectively for the models assuming strong/no mass loss.This results in an average number of PISNe from our stellar populations of the order of ∼ 10 −2 ηIII at a given time 5 (Equation 4).
To compute the average number of PISNe that we expect to find in haloes with stellar mass in the range [M⋆, M⋆ + ∆M⋆] of a given simulated volume V and simulation snapshot, i.e. at a given time t, we sum over all Pop III particles with age τ (260 M⊙) ⩽ tIII ⩽ τ (140 M⊙) in these haloes to recover the total mass MIII,PISN(M⋆) of Pop III stellar populations that can produce PISNe within a time ∆tPISN.The average number of PISNe as a function of M⋆ is then: If we assume the actual number of PISNe k found in a fixed volume V at time t follows a Poisson distribution et al. ( 2020) considered models for Pop IIIs that have lost their outer envelope due to rotationally-induced mixing before exploding as PISNe.P (k; λ) = λ k e −λ /k! with parameter λ ≡ N PISN,t(M⋆), we can estimate the probability of finding at least 1 PISN in haloes with stellar mass within [M⋆, M⋆ + ∆M⋆] as: and if the average number of PISNe that we find in the volume is indeed very small -as we would expect -we will have P (⩾ 1 PISN) ≃ λ ≡ N PISN,t(M⋆).

Expected number of PISNe in galaxies of different mass
We compute the average number of PISNe at a given redshift z as a function of the stellar mass M⋆ of the halo in which they are found, using the method outlined in Section 2.2 (see Equation 4).We consider Pop III stellar particles within all the simulated volumes U6, U7, U8, U10, U12, and U13 to provide the largest possible statistics6 , meaning we are virtually probing a total volume We also focused on six bins of M⋆ spaced by 0.5 dex in the range 7.5 ⩽ LogM⋆/M⊙ < 10.5, for three different redshift points, z = 8.1, 7.3 and 6.7.
The results are shown in the top and middle panels of Figure 3, for the models assuming no/strong mass loss (solid/dotted, black lines).In the plots, the average number of PISNe is normalised by the total number of haloes N haloes in each bin (top panels) and by the total volume V eff (middle panels), to recover the expected number of PISNe per halo/per unit volume; these numbers provide an indication respectively on the number of galaxies one should be looking at to have a reasonable chance of observing PISNe at a given redshift and on the average comoving number density of PISNe expected in a blind survey.Different assumptions on the efficiency ηIII (see Equation 3) and on the IMF ϕ(m) (see Equation 1) are explored.The total number and comoving number density of haloes are also shown in the bottom panels to demonstrate our explored statistics, and we refer to Di Cesare et al. (2023) for a detailed comparison of the dustyGadget galaxy stellar mass functions at various redshifts with observations and with other models and simulations.
We see that the probability of finding PISNe in our simulated volumes is small but non-negligible: in our reference model for Pop III (Salpeter-like IMF in the range [100, 500] M⊙, Pop III star formation efficiency of ηIII = 0.1 and no mass loss), the comoving number density of PISNe can reach values ∼ 10 −1 cMpc −3 , while the average number of PISNe per halo can be up to ∼ 5 × 10 −6 , meaning we would expect about 1 PISN every two hundred thousand haloes, on average (if we take into account the contributions of all the considered stellar mass bins).These numbers depend however on the considered redshift, on the assumed Pop III model and on the halo stellar mass.Particularly, high-mass haloes (109.5 M⊙ ≲ M⋆ ≲ 1010 M⊙) at z = 6.7 seem to be the most favourable for the search of PISNe at high redshifts (N PISN/Nhaloes ≳ 3 × 10 −5 in our reference model, see top panels).In general, the probability of finding PISNe in a given halo increases with M⋆, dropping to zero at the highest masses (M⋆ ≳ 10 9.5 M⊙ at z = 8.1, M⋆ ≳ 10 10 M⊙ at z = 7.3/6.7).On the other hand, given the higher number of haloes at the low-mass end, the expected number density of PISNe in lower-mass galaxies is generally higher, even though the fraction of low-mass galaxies hosting PISNe is lower; the highest number density is found in haloes with M⋆ ≲ 108 M⊙ at z = 6.7 (nPISN ∼ 2 × 10 −8 cMpc −3 in our reference model, see middle panels).We emphasize that the statistics for the highest stellar mass bins is very limited, as no more than a few/a few tens of haloes are found in these bins (see bottom panels), and hence our results are very sensitive to the adopted selection criterion, see e.g. the difference between the no/strong mass loss cases.
While the discrepancies due to the assumption of no/strong mass loss in our estimate of Pop III lifetimes are mostly sub-dominant (apart from the rare massive haloes), a more significant variation is expected when considering different IMFs and/or different values of ηIII.In particular, a factor of ten in the assumed ηIII yields a difference of ∼ 1 dex in the number of PISNe.The scatter resulting from different assumptions on the IMF (covering all the cases explored in Table 1 and shown as a shaded region for the reference case ηIII = 0.1), is instead slightly higher than 1 dex.
Figure 4 shows the average number of PISNe at z = 8.1 that we expect as a function of effective survey volume in our reference model, normalised by the efficiency ηIII; the expected scatter due to different choices of the IMF is also shown as a grey, shaded area.We multiply the values of nPISN in our simulations by the comoving volume of selected JWST surveys at z ≃ 8, with ∆z = 1: (i) the Next Generation Deep Extragalactic Exploratory Public (NGDEEP) Survey (Finkelstein et al. 2021;Pirzkal et al. 2023;Bagley et al. 2023); (ii) the Grism Lens-Amplified Survey from Space (GLASS7 , Treu et al. 2017Treu et al. , 2022;;Castellano et al. 2022 In our reference model, we expect less than 1 PISN on average in all the examined JWST surveys, even when considering all fields together.We also consider a possible ∼ 1 deg 2 survey with the Wide Field Instrument (WFI) of the Nancy Grace Roman Space Telescope12 in the same redshift range: we see that the expected number of PISNe in this volume in our reference model is ≃ 1.5 ηIII.The most pessimistic ) that we expect to find at a given redshift z in haloes within a given range of stellar mass M⋆ as a function of M⋆, computed in six bins of M⋆ (with a spacing of 0.5 dex in the range 7.5 ⩽ LogM⋆/M ⊙ < 10.5).The various colours in the top panels refer to different values of the efficiency η III of a single star-formation episode; η III = 0.01 (gold), η III = 0.1 (black ), and η III = 0.3 (brown).The solid lines refer to the model assuming no mass loss, while the dotted line refers to the model assuming strong mass loss for the case η III = 0.1 (Figure 2).The shaded area also demonstrates the expected scatter around the model assuming no mass loss and η III = 0.1 for different assumed IMFs (see Figure 2 and Table 1).Bottom panels: comoving number density n haloes of haloes found in each M⋆ bin.The total number of haloes found in each bin is also indicated on top of the bin.Results are shown for the combined simulated volumes U6, U7, U8, U10, U12, and U13 at redshifts z = 8.1 (left panels), z = 7.3 (middle panels) and z = 6.7 (right panels).didates during the EoR (also see Section 3.2.1).This would further increase the expected number of PISNe by a factor ten.
We emphasize that in this paper we investigate the possibility of detecting PISNe during the EoR and using them as a marker for Pop III stars.Hence, we focus on galaxies that are potentially observable with JWST and Roman, so that we can also study the underlying stellar populations in details.However, mini-haloes are the most favourable environments for Pop III star formation at high redshifts, meaning we expect a higher number of PISNe occurring in low-mass halos closer to Cosmic Dawn (also refer to the first caveat of Section 4.1).We can compute the cumulative observed The black line refers to the average number of PISNe found in all haloes with 7.5 < Log(M⋆/M ⊙ ) ⩽ 9.5, while the coloured lines refer to the average number found in haloes of different stellar mass bins (see Figure 3).Solid lines refer to our reference model for Pop III stars, i.e. a Salpeter-like IMF in the range [100, 500] M ⊙ with no mass loss, while dotted lines refer to the corresponding model assuming strong mass loss; the shaded area shows the expected scatter for different assumed IMFs (see Figure 1 1).
PISN rate per unit time per unit solid angle at z > 8 as: where shows the brief, "shock-breakout" phase (from figure 3 of Kasen et al. 2011), expanded to display the time evolution in units of hours, while the right-hand side shows the long-term emission in day units (from figure 5 of Kasen et al. 2011).The L bol of galaxies with given M⋆ is computed assuming all the stellar mass consists of Pop II/I stellar populations following a Salpeter IMF in the range [0.1, 100] M ⊙ , and that all stars obey standard massluminosity relations (Equation 10, valid for stars at solar metallicity on the ZAMS, with C = 10 2.68 L ⊙ M −1 ⊙ ).The magenta, shaded region also shows a comparison with the average stellar luminosity of haloes with virial mass ∼ 10 11 M ⊙ at 6 ≲ z ≲ 30, corresponding to M⋆ ∼ 10 8.2 − 10 8.9 M ⊙ , from Riaz et al. (2022, see their figures 1 and 2).
their model15 , reduced to ≃ 10 when accounting for a more extended IMF up to ≃ 260 M⊙, the inferred mass of the PISN progenitor found in Xing et al. (2023) (mod-Larson IMF in Table 1).

Integrated emission
To estimate the chance of observing PISNe against the stellar emission of their hosting galaxies, in Figure 5 we compare the predicted bolometric light curves of representative PISNe with the bolometric luminosity arising from the stellar populations in galaxies of given stellar mass.
We consider PISNe with five different progenitor stars from Kasen et al. (2011), particularly two zero-metallicity stars dying as compact blue super-giants (BSGs, i.e.B250 and B200, respectively with initial masses of 250 M⊙ and 200 M⊙) and three 10 −4 Z⊙ stars dying as red super-giants16 (RSGs, i.e.R250, R200 and R150, with 250 M⊙, 200 M⊙ and 150 M⊙ respectively), with radii 10-50 times larger.The highest luminosities are reached during the "shock-breakout" phase of the explosion, when the radiation-dominated shock, resulting from the expansion of the exploded He core into the hydrogen envelope, approaches the surface of the star and escapes in a luminous X-ray/UV burst.The peak luminosity for each SN can reach values higher than 10 12 L⊙ (see table 2 of Kasen et al. 2011).However, this phase only lasts a few hours.We are instead more interested in the longer-lasting emission following the breakout, as radiation continues to diffuse out of the expanding, cooling ejecta and it may be visible at rest-frame optical wavelengths for several weeks, and even longer in the IR (see e.g.figure 7 of Kasen et al. 2011).The initial phase of the light-curve is powered by the diffusion of thermal energy deposited by the shock and it is dimmer for the BSG models due to the relatively small radii of the progenitors; the secondary peak results from the radioactive decay of synthesized 56 Ni.
To roughly estimate the luminosity arising from the stellar populations in galaxies of given stellar mass, we assume the luminosity of a given star of mass m in the range [mi, mi+1] can be inferred from a simple mass-luminosity relation of the form (see, e.g.Riaz et al. 2022) 17 .The luminosity arising from all stars within this mass range in a given stellar population can then be obtained by integrating over our assumed IMF: Hence, we find a simple relation between the total stellar mass of the haloes M⋆ and their luminosity L⋆: with mϕ(m)dm .Since we are interested in the total emission, as a first approximation we assume all the stellar mass consists of Pop II/I stars and neglect the contribution of Pop III18 , making up less than 5% of the total stellar mass (see e.g.figure 3 of Venditti et al. 2023 and the lower curves of figure 3 in Riaz et al. 2022).We use standard mass-luminosity relations 19for Pop II/I stars (in Eq. 8): (11) where L50 and L0.4 are, respectively, the luminosity of a star with 50 M⊙ and 0.4 M⊙.By integrating these relations over our assumed Salpeter IMF in the range [0.1, 100] M⊙, we find C = 10 2.68 L⊙ M −1 ⊙ (Equation 10).Note that the relations in Equation 11are only valid for the bolometric luminosities of stars on the Zero Age Main Sequence (ZAMS), while older stars would be found at progressively lower luminosities.Therefore, these results should be strictly interpreted as upper limits.Keeping this caveat in mind, we see from Figure 5 that the peak emission from most PISNe would easily outshine the stellar emission of their hosting haloes, at least for M⋆ ≲ 10 9 M⊙, especially when considering the most likely RSG-scenarios.However, as previously emphasised, this emission is shortlived with respect to the total light curve, and our chance of catching the event close to its peak is of the order of ∼ 1 hour/1 yr∼ 10 −4 .Considering the long-term evolution, when the luminosity can drop by more than two orders of magnitude, we have a much less favourable scenario, with only the extreme R250 PISN clearly rising above the stellar emission from M⋆ ≲ 10 8 M⊙ haloes.
Figure 6 shows the observed PISN light curves for the RSG-progenitor models at z ≃ 6, in four JWST/NIRCam filters (F444W, F356W, F277W and F200W), as well as in two filters of the Nancy Grace Roman Space Telescope (F213 and F158).The combination of F356W and F200W was suggested by Hartwig et al. (2018) as the optimal 2-filter diagnostic to detect the prompt emission of PISNe at 6 ≲ z ≲ 12 with JWST, thus identifying possible candidates to be confirmed as transients through follow-up observations.This choice of filters maximizes the visibility time (i.e. the fraction of the light curve that we are able to observe) and thus the probability of finding a PISN, with an optimal exposure time of ∼ 600 s (see figure 6 of Hartwig et al. 2018).Similarly, Moriya et al. (2022b) determined that F158-F213 is the best two-filter combination to discriminate between the classical SNeIa/SNeII at z > 1.5 and PISNe/SLSNe at z > 6 with Roman, by comparing the color-magnitude diagrams of SNe at all phases.
The sensitivity limits 20 in the reference filters of CEERS and of the shallowest/deepest JWST surveys (COSMOS-Web/NGDEEP) included in Figure 4  optimal filters F356W and F200W are not available, we report the limits of the closest filters in terms of wavelength coverage, i.e.F444W and F277W21 .We see that the PISNe arising from the most massive progenitors would be visible from a time ≲ 2 months (in the case of R200 in the COSMOS-Web/F277W filter) up to a time ∼ 9 yr (in the case of R250 in the NGDEEP/F356W filter).The leastmassive-progenitor PISN (R150) would also be visible for almost one year in NGDEEP/F200W and almost two years in NGDEEP/F356W.The limiting sensitivity necessary to separate SNeIa and SNeII at z > 1.5 from PISNe and SLSNe at z > 6 with Roman (see Moriya et al. 2022b) is also indicated: with these sensitivities, the most massive PISN would be observable for about three years in F213 and for about a year and a half in F158.

Resolved observations
In Figure 5 we considered the integrated emission arising from all stellar populations in the galaxies, but in reality it is unlikely that this emission would be concentrated in a small region near their centre.In Venditti et al. (2023) we discussed the complex and irregular morphology of early galaxies in dustyGadget simulations in terms of their gas, dust and stellar components, as they are mostly still in an assembly phase.Moreover, given the low metallicity required for their formation, many Pop IIIs will be found either at the periphery of Pop II clusters or in pristine satellites of their hosting haloes.This could result in lower levels of contamination from higher-metallicity stellar populations to the signal of Pop IIIs and of their PISNe, provided that we are able to observe these galaxies with enough spatial resolution.
In Figure 7 we show the ratio between the distance22 of Pop III particles23 that may produce PISNe at a given redshift (z = 8.1, 7.3 and 6.7) from the centre of mass of their hosting halo and the stellar-mass-weighted radius r⋆, as a function of the halo stellar mass M⋆.We see that there is a large scatter in the values, and that Pop IIIs can be found at a distance up to twenty-five times r⋆, in the most extreme case.To highlight the trend of rIII/r⋆ as a function of M⋆, we also plot the mean and standard deviation in the same M⋆ bins of Figure 3.Most particles have an rIII slightly higher than r⋆, although if we consider a stricter criterion to indicate "external" PISNe (i.e.rIII > 2r⋆), only particles in haloes higher than 10 8.5 M⊙ satisfy this requirement on average, and only particles in haloes higher than 10 9.5 M⊙ at z = 6.7.Note, however, that low-mass galaxies also have a generally lower luminosity, and therefore they are less likely to outshine their hosted PISN (see Figure 5).Despite the large scatter, there is indeed a slightly growing trend of rIII/r⋆ with increasing M⋆.
We compare the physical size of our galaxies to the resolution element of the selected JWST/NIRCam and Roman/WFI filters discussed in Section 3. Pop IIIs with r III above these thresholds (depending on the strictness of the adopted criterion) are considered "external".The mean and standard deviation of r⋆ normalized to the PSF of selected JWST/NIRCam and Roman/WFI filters (F200W/F356W and F158/F213 respectively) are also shown in the plots as thin, coloured lines, with same colour legend of Figure 6.The segments at the top-left corners indicate the relative size of the filters PSF (coloured segments) compared to the average value of r⋆ over the considered sample (gold segments); the absolute size is also indicated next to each segment.
size of the Point Spread Function (PSF) at given redshift is s(z) = θDA(z), with DA(z) the angular diameter distance at redshift z and θ the angular size24 of the PSF.In Figure 7, we show the mean and standard deviation of r⋆ normalized to the PSF in the same bins of M⋆.For reference, we also show the size of the various PSFs relative to the average r⋆ of our PISN-hosting galaxies at each redshift.We see that, on average, the bulk of our galaxies is resolved with more than one pixel in all the considered filters.Particularly, they are best resolved in NIRCam/F200W, while WFI/F213 has the coarsest spatial resolution.High-mass galaxies up to 10 9.5 M⊙ are also generally better resolved than their lower-mass counterparts.

Main caveats
The following caveats should be kept in mind: (i) the statistics for high-mass, well-resolved haloes in our simulations is quite scarce.More reliable results for the highest-stellar-mass bins in Figure 3 would require even larger simulated volumes.We also did not take into account the contribution of galaxies with M⋆ < 10 7.5 M⊙, that are not well resolved in our simulations due to our mass resolution limits, and that would also be the most affected by radiative feedback (see point (v)).Neglecting this halo population, that hosts a significant portion of the total Pop III mass (≳ 60% at the considered redshifts, see Venditti et al. 2023), will result in an underestimation of the total number of expected PISNe, at least by a factor 1.5.However, we remark that we are mainly interested in PISN hosts in which we can potentially study the underlying stellar populations in details, to look for candidate Pop III stars during the EoR.A full picture on all the PISNe produced in our cubes regardless of their environment is beyond the scopes of the present work.
(ii) In Equation 4, we are assuming for simplicity that PISNe have uniform probability of occurring within ∆tPISN for a given stellar population.However, in reality this probability depends on the assumed IMF; in particular, with our IMF, lower-mass Pop III stars are more likely to form, therefore we should expect more PISNe arising from stars closer to the lower-mass threshold for PISN progenitors, and hence having longer lifetimes.Moreover, a random sampling of the IMF would be required for a more realistic estimate of the average number of PISNe produced per stellar population.
(iii) We have explored two different models with different assumptions for the stellar mass loss arising from Pop IIIs, however, many other uncertainties exist.For example, we have not explored different choices for the value of Zcrit, although the chosen assumptions may change the results for Pop III star formation at cosmological scales considerably (see e.g.Maio et al. 2010).We did consider the impact of choosing different IMFs in the estimate of the average number of PISNe produced per unit Pop III mass (Equation 1).Nonetheless, we emphasize that a full exploration of the impact of IMF variations would require a self-consistent approach.
(iv) Our chosen Pop III IMF covers the mass range [100, 500] M⊙.On the other hand, studies on the properties of carbon-enhanced extremely metal-poor stars show that they are best matched by the nucleosynthetic yield of traditional core-collapse SNe25 (Iwamoto et al. 2005;Keller et al. 2014;Ishigaki et al. 2014;de Bennassuti et al. 2014de Bennassuti et al. , 2017;;Fraser et al. 2017;Magg et al. 2022;Aguado et al. 2023b), suggesting that Pop III progenitors with masses ∼ 10−40 M⊙ are indeed more prevalent.An even lower minimum value of ∼ 1 M⊙ has been suggested by several works of stellar archaeology at a high confidence level (Hartwig et al. 2015;Rossi et al. 2021).Here, we are not interested in the detailed nucleosynthetic pattern, hence we do not require a high precision on the IMF as in classical archaeology studies and we do not expect our conclusions to be significantly affected by this assumption.Indeed, the main impact of the IMF on our results is through the parameter N PISN/MIII, and the value of 2.17 × 10 −3 M⊙ that we find with our assumed IMF is within the allowed range spanned by IMFs extended down to ∼ 1 M⊙.We also note that, while considering an IMF more biased towards massive stars might cause in principle an overproduction of PISNe, the reality is more complex.In fact, given the exceptionally high metal yield of PISNe, this choice may actually lead in the opposite direction, as a faster enrichment might induce an earlier transition to the Pop II mode of star formation.The effect of the Pop III IMF on cosmic star formation history is debated (Maio et al. 2010(Maio et al. , 2016;;Pallottini et al. 2014).In Venditti et al. (2023), we found that our Pop III SFRD is in line with other models and simulations, although large uncertainties are involved.
(v) Even though the simulation includes an homogeneous UV background as in Haardt & Madau (1996), this does not have feedback on cosmic star formation at z > 6, i.e. in the epoch considered in the present work.The lack of a proper treatment of radiative feedback, -especially in the UV and LW bands -certainly affects the late Pop III star formation.This limitation is extensively discussed in Venditti et al. (2023), in which our Pop III star formation rate density is also carefully compared with the results of other small-scale simulations modelling the local/global effects of UV and LW radiation.
(vi) At the current stage, dustyGadget does not include a model for metal mixing and turbulent metal diffusion below our gas mass resolution.While allowing for incomplete inhomogeneous mixing at the sub-grid level may allow to enhance Pop III star formation by a factor 2-3 (Sarmento et al. 2016(Sarmento et al. , 2017(Sarmento et al. , 2018;;Sarmento & Scannapieco 2022), the study by Su et al. (2017) showed that sub-grid turbulent metal diffusion has a negligible impact on general star formation and ISM properties.The effect of diffusion has also been shown to be weaker at lower metallicity and relatively unimportant for Pop III star formation (Jeon et al. 2017).
We refer to Venditti et al. (2023) for further discussion and for a comparison with the models of Sarmento et al. (2018) and Sarmento & Scannapieco (2022).
(vii) We have not discussed our ability to discriminate PISNe from other theoretical red sources such as dark stars (Spolyar et al. 2008;Freese et al. 2008Freese et al. , 2010;;Yoon et al. 2008;Natarajan et al. 2009;Hirano et al. 2011;Banik et al. 2019) and direct-collapse black holes (Haemmerlé et al. 2020).We note, however, that we should be able to distinguish these sources through follow-up observations after a few years, due to the time variability of PISNe.

Implications of future detections
In light of our results, we wish to discuss the implications of future detections/non-detections of PISNe with both JWST and Roman.In our model, we predict that we should be able to see some PISNe at high-redshift at least through wide surveys with the Roman Space Telescope.Hence, the following questions arise: what if we do not find any PISN at high redshift?What if, conversely, we find a surprisingly high number of PISNe?
A higher number of PISNe than allowed by our models may be explained if a substantial fraction of Pop II stars is also able to reach sufficiently high masses.The JWST capabilities for discovering PISNe from 0.14-0.43Z⊙ stars at z < 10 have been discussed e.g. in Regős et al. (2020).Although the general consensus is that Pop II/I stars usually follow a standard Salpeter IMF, with masses ∼ 0.1−100 M⊙, our constraints on the IMF are mainly derived from observations of the local Universe.On the other hand, more topheavy IMFs might actually be more common at high redshifts (see e.g.Chon et al. 2022).The prevalence of a topheavy IMF has been suggested to alleviate the tension between JWST observations and existing models and simulations, which seem to systematically under-produce massive, bright galaxies at high redshifts (z ≳ 10, Yung et al. 2023;Harikane et al. 2023;Trinca et al. 2023).Olivier et al. (2022) also demonstrated that a high-temperature black-body spectrum (T eff ∼ 80000 K, associated to very massive stars) is necessary to explain the ionization level of low-z, extremeemission-line galaxies that are considered analogues to EoR galaxies.We further note that we focused on our chance of observing PISNe through their prompt emission over a time ∆tprompt ∼ 1 yr, hence neglecting the possibility that PISNe may be observable through their afterglow on much longer timescales.Hartwig et al. (2018) predict that the PISN afterglow is not accessible through current observational facilities, although observatories and telescopes with stronger capabilities may be able to detect it in the future.However, more accurate modelling may find that this signal is actually brighter than these predictions, implying that the ∆tprompt −→ ∆t afterglow in Equation 4 -and hence the expected number of PISNe -can be boosted up to a hundred times.
Not finding any PISNe during the EoR with either JWST or Roman would also be surprising, and it would require a revision of our models for Pop III stars.Maybe Pop IIIs do not reach high enough masses to produce a significant number of PISNe, either because the physics of star-forming clouds only allows stellar masses lower than ∼ 140 M⊙, or because the stars lose most of their mass be-fore dying.An alternative explanation is that Pop III star formation is suppressed much earlier than z ∼ 10.As most models and simulation currently predict that the inhomogeneous metal enrichment should allow a late Pop III star formation, our models for chemical and radiative feedback may also need to be reconsidered in this scenario.

Detection strategies: archaeology vs. direct detections at high z
In the introduction of this paper we listed many possible strategies to study PISNe and their Pop III progenitors.These include direct detection, especially at high redshift where massive Pop III stars may still be actively forming, and cosmic archaeology (Aoki et al. 2014;Salvadori et al. 2007;Yoshii et al. 2022;Aguado et al. 2023a;Xing et al. 2023).Both techniques present unique advantages and challenges.
The first obvious advantage of cosmic archaeology is that it makes the study of PISNe accessible through observations in the local Universe.However, the nucleosynthetic footprint of low-metallicity gas may result from a complex stratification of multiple sources of metal pollution, that can be hard to disentangle (e.g.Ji et al. 2015;Vanni et al. 2023).Moreover, an inherent issue lies in the selection of candidate PISN descendants: given its significant metal yields, a single PISN may cause its environment to "overshoot" and promptly reach a relatively high metallicity of ≳ 10 −3 Z⊙ (Karlsson et al. 2008;Greif et al. 2010;Wise et al. 2012), and all PISN-descendant candidates to date (Aoki et al. 2014;Salvadori et al. 2019;Xing et al. 2023;Aguado et al. 2023a) exhibit in fact [Fe/H] > −2.5.This means we should not only look at the most metal-poor stars/environments, rendering the identification and interpretation of potential candidates even more difficult (see also de Bennassuti et al. 2017).We certainly need improved selection strategies (Aguado et al. 2023a), possibly relying on photometry rather than costly spectroscopy on each source.
Throughout this paper, we demonstrated that a direct detection of PISNe during the EoR is not beyond hope.The identification of PISNe in this case is probably more straightforward and it has been already discussed in the literature (Kasen et al. 2011;Hartwig et al. 2018;Moriya et al. 2022b).We also note that no SN of any kind has yet been found at z ≳ 3 (e.g.Cooke et al. 2012;Jones et al. 2013).Due to our sensitivity limits, it is actually more likely that if we do find SNe at high redshifts, these will be either PISNe or other kinds of SLSNe.The main obstacle is the fundamentally transient nature of SNe: although PISNe light curves can be stretched up to ∼ 10 yr in the observer frame, this transient nature, combined with their rarity, make their detection less likely.A longer visibility time might be attained through their afterglow (see e.g. the discussion in Section 4.2).The direct imprints of SN feedback may also be studied as a fossil record for longer time scales after the explosion.If we consider for example a galaxy with M⋆ ∼ 10 7 M⊙, of which ∼ 10 4 M⊙ in Pop IIIs, ∼ 40% of the Pop III mass would end up as PISNe with our assumed IMF, while ∼ 10% of the Pop II mass would end up as a traditional core-collapse SNe.By assuming a return fraction of ∼ 0.1 for core-collapse progenitors and ∼ 0.5 for PISN progenitors (Karlsson et al. 2013), the ratio of metals from PISNe vs core-collapse would be of the order of 2 × 10 −2 , i.e. non negligible.This means we may be able to reliably tell apart the chemical signature of PISNe outflows.A thorough investigation of such matters is left to future works.
Looking for PISNe at even higher redshifts (z ≳ 10) may also be feasible thanks to the efficient Pop III formation in low-mass, pristine halos towards Cosmic Dawn (see e.g. the discussion at the end of Section 3.1).We did not investigate such possibility in depth as the galaxies hosting these PISNe are less likely to be observable, while in the current work we were mainly interested in using PISNe as tracers for Pop III stars.Nonetheless, the tentative detection of HeII emission in the vicinity of a luminous z = 10.6 galaxy (Maiolino et al. 2023) might already indicate the presence of observable Pop III clusters at these redshifts.Dedicated surveys aimed at finding PISNe at very high redshifts may reveal even more of these fascinating events.
In conclusion, many promising strategies to look for PISNe at high redshifts are now on the horizon.So far, stellar archaeology has been our best tool.While this technique will continue to provide invaluable information, ultimately, the combination of all the proposed diagnostics of PISNe will serve our purposes, as they can provide independent constraints on the rate of PISNe across time and on the properties of their Pop III progenitors.We stress that a synergistic approach, exploiting the different strengths of all these methods, will give us the best chance to gain a complete picture on massive Pop III stars.

SUMMARY AND CONCLUSIONS
In this paper we studied the probability of finding PISNe arising from Pop III stars during the EoR.Starting from a suite of six 50h −1 cMpc dustyGadget simulations, we provided indications on the most promising galaxy candidates with M⋆ > 10 7.5 M⊙ to look for PISNe at 6 ≲ z ≲ 8.We also provided predictions on the expected number of PISNe in selected JWST surveys and in a future ∼ 1 deg 2 survey with the Nancy Grace Roman Space Telescope.Finally, we discussed the observability of PISNe within their hosting galaxies in these surveys, in terms of both their integrated and resolved stellar emission.We considered many possible sources of uncertainties in our results, including the effect of mass loss on Pop III lifetimes, the impact of choosing different IMFs on the expected number of PISNe per unit Pop III mass, and the possibility of a reduced Pop III star formation efficiency with respect to the nominal efficiency of the simulation.
We find that, although the probability of observing PISNe is indeed small, it is non-negliglible.In our reference model for Pop III 500] M⊙, Pop III star formation efficiency of ηIII = 0.1 and no mass loss), the comoving number density of PISNe and the average number of PISNe per halo can be up to ∼ 10 −1 cMpc −3 and ∼ 5 × 10 −6 (i.e. about 1 PISN every two hundred thousand haloes, on average) respectively.A higher number density of PISNe is generally found in haloes with decreasing stellar mass, given the higher absolute number density of low-mass haloes at all the considered redshifts.However, the relative fraction of low-mass haloes hosting PISNe is smaller with respect to higher-mass haloes, and hence the probability of finding PISNe in a given halo increases with M⋆, and it is highest in haloes with 10 9.5 M⊙ ≲ M⋆ ≲ 10 10 M⊙ at z = 6.7.Therefore, different observing strategies may be considered: either targeted follow-up observations of candidate high-mass galaxies -with a higher PISNe-hosting fraction -, or blind wide-field surveys -also including the numerous low-mass galaxies that can also potentially host PISNe.
PISNe from progenitors with masses higher than 200 M⊙ would be observable for a variable time range (from ∼ 2 months up to ∼ 9 yr) through all the considered JW-ST/NIRCam filters, particularly the F356W and F200W filters that have been indicated by Hartwig et al. 2018 as the best two-filter diagnostic to identify possible PISNe candidates, and the F444W and F277W filters available for the largest considered JWST survey (COSMOS-Web).The 150 M⊙-progenitor PISN would also be visible for almost one/two years in the F200W/F356W filters of the deepest considered JWST survey (NGDEEP).However, detecting PISNe with JWST may still be challenging due to their rarity: in our reference model, we expect on average less than 1 PISN in all the examined JWST surveys.
A higher potential might be obtained through the WFI instrument on board of Roman, thanks to its large field-ofview.By considering an example ∼ 1 deg 2 survey, we see that ≃ 1.5 ηIII PISNe are expected in our reference model.An even higher survey volume would further increase the number of expected PISNe.For example, the ∼ 10 deg 2 survey suggested by Moriya et al. 2022b with limiting magnitudes of 27.0 mag and 26.5 mag in the F158 and F213 filters (low enough to see at least the 250 M⊙-progenitor PISN for about a year and a half/three years in F158/F213) would further increase this number by a factor ten.More favourable scenarios are also obtained when considering different Pop III IMFs and/or higher Pop III formation efficiencies, and when including the contribution of coarsely-resolved environments that have not been specifically targeted in this study.
While the integrated flux of the underlying galaxies might exceed the flux of PISNe in the considered filters, Pop III stars -and hence PISNe -are usually found at the periphery of their hosting haloes.This should mitigate the contamination of the signal arising from stars, especially for the most massive galaxies in which the stellar emission is more likely to outshine PISNe.The bulk of most PISNhosting galaxies would also be resolved in JWST/NIRCam and Roman/WFI at these redshifts.
We remark once again that if we do find a galaxy hosting a PISN at high redshifts, this incredible accomplishment may also pave the way for an even more incredible discovery.Indeed, by observing the target galaxy again after some years, when the SN has faded considerably, we might even be able to discern an active underlying Pop III stellar population.This would possibly lead to the first clear detection of Pop III stars in the history of astronomy.

Figure 1 .
Figure 1.Compilation of various theoretical models for the Pop III IMF, ϕ(m), normalised to a total mass of 2 × 10 6 M ⊙ (i.e.approximately our Pop III stellar particle mass resolution).The models shown here are from Schauer et al. (2022) (blue), Hartwig et al. (2018) (orange) and Jaacks et al. (2018) (pink ), and they can all be expressed through the general expression at the bottom (Equation2), with parameters specified in Table1.Our Salpeterlike IMF(Graziani et al. 2020;Di Cesare et al. 2023;Venditti et al. 2023) is shown as a black, solid line.We also include IMFs with a broader mass range, extending up to 1000 M ⊙ , as allowed from ab-initio simulations of star-forming clouds (e.g.Susa et al. 2014;Hirano et al. 2015a; Hosokawa et al. 2016, grey).The gray shaded region indicates the stellar progenitor mass range of PISNe.

Figure 2 .
Figure 2. Pop III main-sequence lifetimes τ as a function of the initial mass m of the stars, fromSchaerer (2002), in the mass range of our assumed IMF.The black, solid line refers to the model assuming no mass loss (see table 4 of the original paper), while the black, dotted line refers to the model assuming strong mass loss (see their table5).The interpolated lifetimes at the limits of the PISN progenitors mass range[140, 260]  M ⊙ (grey, shaded area) are indicated for the two models.The average lifetimes τ of Pop III stellar populations for the two models with our assumed IMF are also shown, on top of the horizontal, solid/dotted, black lines.

Figure 3 .
Figure3.Top/middle panels: Average number of PISNe per halo (N PISN,t /N haloes )/per unit volume (n PISN,t ) that we expect to find at a given redshift z in haloes within a given range of stellar mass M⋆ as a function of M⋆, computed in six bins of M⋆ (with a spacing of 0.5 dex in the range 7.5 ⩽ LogM⋆/M ⊙ < 10.5).The various colours in the top panels refer to different values of the efficiency η III of a single star-formation episode; η III = 0.01 (gold), η III = 0.1 (black ), and η III = 0.3 (brown).The solid lines refer to the model assuming no mass loss, while the dotted line refers to the model assuming strong mass loss for the case η III = 0.1 (Figure2).The shaded area also demonstrates the expected scatter around the model assuming no mass loss and η III = 0.1 for different assumed IMFs (see Figure2and Table1).Bottom panels: comoving number density n haloes of haloes found in each M⋆ bin.The total number of haloes found in each bin is also indicated on top of the bin.Results are shown for the combined simulated volumes U6, U7, U8, U10, U12, and U13 at redshifts z = 8.1 (left panels), z = 7.3 (middle panels) and z = 6.7 (right panels).

Figure 4 .
Figure 4. Average number of PISNe N PISN,t (normalised to the Pop III formation efficiency η III ) as a function of the effective survey volume V eff at z = 8.1.Vertical, dashed lines indicate the effective volume of selected JWST surveys and their cumulative volume at z ≃ 8, with ∆z = 1 (filled circles, see text for details), and the comoving volume of an equivalent ∼ 1 deg 2 survey with the Roman Space Telescope (RST, filled triangles); the average number of PISNe found in our six dustyGadget cubes is also indicated in the plot (DG, empty squares).A horizontal, dashed line further indicates the reference value of 1 PISN per volume.The black line refers to the average number of PISNe found in all haloes with 7.5 < Log(M⋆/M ⊙ ) ⩽ 9.5, while the coloured lines refer to the average number found in haloes of different stellar mass bins (see Figure3).Solid lines refer to our reference model for Pop III stars, i.e. a Salpeter-like IMF in the range [100, 500] M ⊙ with no mass loss, while dotted lines refer to the corresponding model assuming strong mass loss; the shaded area shows the expected scatter for different assumed IMFs (see Figure1and Table1).As a reference, the purple pentagon indicates the cumulative number of PISNe at z > 8 expected in a survey of ∼1 deg 2 in one year, computed by integrating over the observed PISN rate as a function of redshift from Jaacks et al. (2019, see their figure15) in this redshift range, and re-normalizing by the considered survey area and by the extended IMF (see text and Table1).

Figure 7 .
Figure7.Ratio between the Pop III distance from the centre of mass of the hosting halo r III and the stellar mass-weighted radius r⋆ of the halo, for each Pop III stellar particle that could potentially produce PISNe (gold stars, i.e.Pop III particles with age τ noML (260 M ⊙ ) ⩽ t III ⩽ τ strongML (140 M ⊙ )).The mean and standard deviation computed in five bins of M⋆ (with a spacing of 0.5 dex in the range 7.5 ⩽ LogM⋆/M ⊙ < 10) is also shown as a gold, thick, solid line.Results are reported at z = 8.1 (top panel), z = 7.3 (middle panel) and z = 6.7 (bottom panel).The black, dashed lines indicate where r III = r⋆ or r III = 2r⋆; Pop IIIs with r III above these thresholds (depending on the strictness of the adopted criterion) are considered "external".The mean and standard deviation of r⋆ normalized to the PSF of selected JWST/NIRCam and Roman/WFI filters (F200W/F356W and F158/F213 respectively) are also shown in the plots as thin, coloured lines, with same colour legend of Figure6.The segments at the top-left corners indicate the relative size of the filters PSF (coloured segments) compared to the average value of r⋆ over the considered sample (gold segments); the absolute size is also indicated next to each segment.
Moriya et al. (2022b)g a Salpeter-like IMF in the range [1, 1000] M⊙, see Table1) yields 0.1 ηIII PISNe, while this number can grow up to ∼ 2.3 ηIII in the most optimistic case 13 (i.e.assuming a power-law IMF with α = −0.17 in the range [1, 260] M⊙).Even higher survey volumes with Roman have been proposed in the literature:Moriya et al. (2022b), for example, suggested a 10 deg 2 transient survey with a limiting magnitude of 27.0 mag and 26.5 mag in the F158 and F213 filters respectively, to be conducted for five years with a cadence of one year with the aim of looking for PISNe can- 13When the original IMFs in the range [1, 150] M ⊙ from Jaacks et al. (2018) are considered, the upper limit on the total number of PISNe can grow up to ∼ 7.6 η III .
Hartwig et al. 2018 and curve of PISNe from the different RSG progenitors of Figure5in the same F356W and F200W filters (thick, yellow lines and thick, cyan lines), in the F444W and F277W filters (thin, red lines and thin, dark-green lines), and in the F213 and F158 filters of the Roman Space Telescope (thick, orange lines and thick, aquamarine lines), at z ≃ 6, as a function of observer-frame time (from figure4ofHartwig et al. 2018 and figure 6 of Moriya et al. 2022b).The sensitivity limits of selected JWST surveys and the limiting sensitivity necessary to separate SNeIa and SNeII at z > 1.5 from PISNe and SLSNe at z > 6 with Roman (see text for details) are indicated as horizontal, dashed lines, with the same colour legend used for the filters.