# The energy distribution of the first supernovae

I. Koutsouridou, S. Salvadori, Á. Skúladóttir, M. Rossi, I. Vanni, G. Pagnini

► To cite this version:

I. Koutsouridou, S. Salvadori, Á. Skúladóttir, M. Rossi, I. Vanni, et al.. The energy distribution of the first supernovae. Monthly Notices of the Royal Astronomical Society, 2023, 525, pp.190-210. 10.1093/mnras/stad2304. insu-04848714

HAL Id: insu-04848714

https://insu.hal.science/insu-04848714v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The energy distribution of the first supernovae

I. Koutsouridou,<fnref n="1" />,<fnref n="2" />* S. Salvadori,<fnref n="1" />,<fnref n="2" /> Á. Skúladóttir,<fnref n="1" />,<fnref n="2" /> M. Rossi,<fnref n="1" />,<fnref n="2" /> I. Vanni and G. Pagnini<fnref n="3" />

Accepted 2023 July 25. Received 2023 July 25; in original form 2023 May 4

## ABSTRACT

The nature of the first Pop III stars is still a mystery and the energy distribution of the first supernovae is completely unexplored. For the first time we account simultaneously for the unknown initial mass function (IMF), stellar mixing, and energy distribution function (EDF) of Pop III stars in the context of a cosmological model for the formation of a MW-analogue. Our data-calibrated semi-analytic model is based on a N-body simulation and follows the formation and evolution of both Pop III and Pop II/I stars in their proper time-scales. We discover degeneracies between the adopted Pop III unknowns, in the predicted metallicity and carbonicity distribution functions and the fraction of C-enhanced stars. None the less, we are able to provide the first available constraints on the EDF, $dN/dE_* \propto E_*^{-\alpha_e}$ with $1 \le \alpha_e \le 2.5$. In addition, the characteristic mass of the Pop III IMF should be $m_{ch} < 100 \text{ M}_\odot$, assuming a mass range consistent with hydrodynamical simulations (0.1–1000 M}_\odot). Independent of the assumed Pop III properties, we find that all [C/Fe] > +0.7 stars (with [Fe/H] < -2.8) have been enriched by Pop III supernovae at a > 20 per cent level, and all [C/Fe] > +2 stars at a > 95 per cent level. All very metal-poor stars with [C/Fe] < 0 are predicted to be predominantly enriched by Pop III hypernovae and/or pair instability supernovae. To better constrain the primordial EDF, it is absolutely crucial to have a complete and accurate determination of the metallicity distribution function, and the properties of C-enhanced metal-poor stars (frequency and [C/Fe]) in the Galactic halo.

**Key words:** stars: Population III – Galaxy: abundances – Galaxy: formation – Galaxy: halo – galaxies: high-redshift – dark ages, reionization, first stars.

## 1 INTRODUCTION

The formation of the first stars marked a fundamental shift in the history of our universe; from a simple, homogeneous and isotropic state to the structured complexity we observe today. Their light brought an end to the so-called dark ages and initiated a period of reionization and heating of the intergalactic medium (IGM). Formed out of purely metal-free primordial gas, the first stars are commonly referred to as Population III (Pop III) stars to distinguish them from the subsequent generations of Pop II (or metal-poor; $Z \le 0.1 \text{ Z}_\odot$) stars and Pop I (or metal-rich) stars. The Pop III stars were the sources of the first metals (i.e. elements heavier than lithium) and dust grains, which they released in the surrounding medium via supernova (SN) explosions and stellar winds.

Although critical for our understanding of the early Universe, the nature and characteristics of the first stars remain elusive. As of today, no metal-free star has been observed and theoretical models have thus far failed to reach a consensus (e.g. Bromm 2013 and Klessen 2019 for recent reviews). Initially, first stars were thought to be very massive (~100–1000 M$_\odot$), owing to the lack of metals and dust grains, which are more efficient coolants than molecular hydrogen,<fnref n="1" />

and thus facilitate the fragmentation of gas clouds (Omukai & Nishi 1998; Omukai & Palla 2001; Abel, Bryan & Norman 2002; Bromm, Coppi & Larson 2002; O’Shea & Norman 2007). Later, more detailed calculations showed that Pop III stars can have lower masses, of the order of some tens of solar. After the initial spherically distributed gas infall, material falling with non-negligible angular momentum starts building up a rotationally supported disc around the central proto-stellar core. Radiative feedback from the nascent proto-star plays a key role in regulating this accretion process, being able to clear out the accretion disc when the star reaches a mass of ≈ 30–40 M$_\odot$ (McKee & Tan 2008; Hosokawa et al. 2011). Cosmological simulations that include radiative feedback confirm this picture and show that the mass spectrum of Pop III stars is broader (~10–1000 M$_\odot$) than previously thought (Hirano et al. 2014, 2015). Furthermore, detailed 3D simulations investigating the formation of the first proto-stars, show that the accretion disc can be highly susceptible to fragmentation, leading to the formation of sub-solar fragments (see, e.g. Machida, Matsumoto & Inutsuka 2008; Clark et al. 2011; Greif et al. 2011; Dopcke et al. 2013; Stacy, Bromm & Lee 2016; Wollenberg et al. 2020). However, it is still not settled whether most of these fragments migrate inwards and merge together, or are expelled from the system to survive as low mass stars (Hosokawa et al. 2016; Hirano & Bromm 2017). Ultimately,

<footnote n="1">¹Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, Italy</footnote>

<footnote n="2">²INAF/Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="3">³GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, 92195 Meudon, France</footnote>

<footnote>* E-mail: ikoutsouridou@gmail.com</footnote>

<footnote n="1">¹Under certain circumstances, such as the suppression of $H_2$ cooling by an external radiation field (Latif et al. 2013; Agarwal et al. 2016), the critical mass for collapse can increase even further resulting into supermassive stars</footnote>

<footnote>$(m \gtrsim 10^6 \text{ M}_\odot)$, which could in turn be the progenitors of supermassive black holes, powering the first quasars (Mortlock et al. 2011; Bañados et al. 2018).</footnote>
although the mass range of Pop III stars is still largely unknown, it
is probably biased towards massive stars, and most likely extremely
broad, from ≈ 1000 M⊙ down to <1 M⊙.

At the higher mass end, the evolution of Pop III stars is fairly
well-understood. In the absence of rotation, stars with initial masses
$\gtrsim 260$ M$_{\odot}$ collapse directly into black holes swallowing all their
heavy element production (Fryer, Woosley & Heger 2001). Between
~140 and ~260 M$_{\odot}$, stars explode as pair instability supernovae
(PISNe), with explosion energies and ejecta depending on their initial
mass (Heger & Woosley 2002; Takahashi, Yoshida & Umeda 2018).
The mechanism involves electron-positron pair production, which
lowers the internal radiation pressure supporting the star against
gravitational collapse. This pressure drop leads to a rapid contraction,
which in turn ignites a runaway thermonuclear explosion that blows
apart the star completely, leaving no stellar remnant behind (e.g.
Barkat, Rakavy & Sack 1967; Bond, Arnett & Carr 1984). At lower
masses, but higher than ~100 M$_{\odot}$, pair-instability still takes place but
the net reduction in pressure is not sufficient to lead to the complete
disruption of the star. Instead, after a series of pulsations the star
produces a large iron core that likely collapses to a black hole,
swiping most heavy elements inside (Heger & Woosley 2002).

Contrary to PISNe, the explosion of stars with initial masses $m_* \sim 10-100$ M$_{\odot}$ leaves behind a stellar remnant and thus the final ejecta can vary strongly according to the amount of mixing and fallback,<fnref n="2" /> which depend on many factors, including the explosion energy and rotation (e.g. Nomoto et al. 2006; Heger & Woosley 2010; Nomoto, Kobayashi & Tominaga 2013). Therefore, even if the mass distribution of Pop III stars were to be established theoretically, it is not yet settled whether these additional parameters are mass-dependent or follow separate distributions.

Currently there is no 'standard model' that describes the properties
of Pop III stars. However, one can infer indirect constraints from
stellar archaeology, i.e. the study of the oldest, most metal-poor
stars in the Milky Way and its dwarf satellites (e.g. Frebel &
Norris 2015). These long-lived stars preserve in their atmospheres
the chemical abundance patterns of their birth gas clouds, which
were polluted by their ancestors: the first stars. One of the most
interesting populations among them are the carbon-enhanced metal-
poor (CEMP) stars, which are characterized by high relative carbon
abundances [C/Fe] $\ge$ +0.7 (e.g. Aoki et al. 2007). The CEMP stars
are commonly divided into two major sub-classes (Aoki et al. 2002;
Beers & Christlieb 2005): (1) those showing an excess of heavy
elements produced by slow neutron-capture processes (CEMP-s);
and (2) those showing no such enhancement (CEMP-no stars). The
two sub-classes are linked to different formation scenarios. CEMP-s
stars represent $\gtrsim 50$ per cent of all CEMP stars but are extremely
rare at [Fe/H] $< -3$ (Norris et al. 2013; Yoon et al. 2016). Their
abundance pattern and the fact that $\gtrsim 80$ per cent of them are
members of binary systems (Starkenburg et al. 2014; Hansen et al.
2016b) are consistent with them being C-enhanced by mass transfer
from an asymptotic giant branch (AGB) companion star (e.g. Aoki
et al. 2007; Abate et al. 2015). On the contrary, CEMP-no stars are
dominant among the most iron-poor stars: at least 12 out of the 14
known stars with [Fe/H] $< -4.5$ are CEMP-no stars (see Section 3;
Vanni et al. submitted). The CEMP-no stars are not preferentially
found in binary systems (Starkenburg et al. 2014; Hansen et al.
2016a) or, even when they are, show no trace of mass transfer by a

companion (Aguado et al. 2022), and so their chemical abundances
are representative of their birth environment.

CEMP-no stars have been observed in large numbers in the Galactic halo (e.g. Carollo et al. 2012; Norris et al. 2013; Yong et al. 2013b; Placco et al. 2014; Bonifacio et al. 2015; Yoon et al. 2016) and in ultra faint dwarf (UFD) satellites (e.g. Norris et al. 2010; Lai et al. 2011; Gilmore et al. 2013; Frebel, Simon & Kirby 2014; Ji et al. 2016; Spite et al. 2018), but also in more massive and luminous dwarf spheroidal (dSph) galaxies (Skúladóttir et al. 2015, 2023; Susmitha, Koch & Sivariani 2017; Chiti et al. 2018; Yoon et al. 2020) and the Galactic bulge (Howes et al. 2016; Arentsen et al. 2021). Their proposed zero-metallicity progenitors include: (i) massive 'spinstars' with internal mixing and mass loss (Meynet, Ekström & Maeder 2006; Maeder, Meynet & Chiappini 2015; Liu et al. 2021), which are now strongly disfavoured by the measured high values of $^{12}$C/$^{13}$C in several CEMP-no stars (Aguado et al. 2022, 2023); (ii) low energy faint ($E_{51} = E_*/10^{51}$erg < 1) or normal ($E_{51} \sim 1$) core-collapse supernovae (ccSNe) with mixing and fallback that release small amounts of iron and large amounts of carbon and other light elements (e.g. Umeda & Nomoto 2003; Iwamoto et al. 2005; Cooke & Madau 2014; Marassi et al. 2014; Tominaga, Iwamoto & Nomoto 2014; Salvadori, Skúladóttir & Tolstoy 2015; Vanni et al., submitted).

However, higher energy progenitors have also been found con-
sistent with metal-poor stars. Ishigaki et al. (2014) found that the
abundance pattern of the most iron-deficient star observed (SMSS
J031300.36-670839.3; Keller et al. 2014) is well-reproduced both by
a $E_{51} = 1$ ccSN and by a $E_{51} \ge 10$ Pop III hypernova. Ishigaki et al.
(2018) suggested that more than half of a sample of ~200 extremely
metal-poor (EMP; [Fe/H] < -3) literature stars are best fitted by
a $m_* = 20$ M$_{\odot}$ ($E_{51} = 10$) hypernova model. Placco et al. (2015),
however, who analysed a subset of these stars, found systematically
lower energy-progenitors for them (see table 2 in Ishigaki et al.
2018). Ezzeddine et al. (2019) argued that the CEMP-no star HE
1327-2326 shows an imprint of a $E_{51} = 5$ asymmetric hypernova
explosion. More recently, Skúladóttir et al. (2021) and Placco et al.
(2021) discovered two ultra metal poor stars ([Fe/H] < -4) in the
Sculptor dSph galaxy and in the Galactic halo, respectively, with
very low [C/Fe] and abundance patterns indicating that they descend
from $E_{51} = 10$ Pop III hypernovae (see also Skúladóttir et al. 2023).

Besides identifying individual Pop III progenitors, several studies
have employed galaxy formation models of Milky Way (MW)-
analogues to investigate the properties of Pop III stars in a statistical
manner, including their spatial distribution (e.g. White & Springel
2000; Brook et al. 2007; Salvadori et al. 2010; Tumlinson 2010;
Ishiyama et al. 2016; Starkenburg et al. 2017; Hartwig et al. 2022)
and the impact of their initial mass function (IMF) on the abundances
of metal-poor stars (e.g. Salvadori, Schneider & Ferrara 2007;
Komiya et al. 2010; de Bennassuti et al. 2017; Hartwig et al.
2018; Sarmento, Scannapieco & Côté 2019; Tarumi, Hartwig &
Magg 2020). However, none of these studies has examined the full
parameter space of Pop III stars—IMF, mixing, explosion energy—
and none has considered the existence of high $E_{51}$ primordial SNe.

In this work, we aim to fill this gap and explore how varying the
energy distribution of the first SNe affects the properties of the stars
surviving until $z=0$ in the Galactic halo. To this end, we develop a
new semi-analytic model (SAM), named NEFERTITI (NEar-FiEld
cosmology: Re-Tracing Invisible Tlmes), that traces the formation
and evolution of individual Pop III and Pop II/I stars, accounting for
all the unknown related to primordial star-formation. NEFERTITI
can run on any merger-tree or dark matter (DM) simulation to shed
light on the earliest phases of star-formation. Here, we combine it

<footnote n="2">²Fallback refers to the material that collapses to form a neutron star when $m_* \le 20-30$ M$_{\odot}$, or a black hole when $m_* \gtrsim 20-30$ M$_{\odot}$ (Colgate & White 1966; Heger & Woosley 2010; Chan et al. 2018).</footnote>
with a DM simulation of a MW analogue. This allows us to follow
in detail the early chemical evolution of the first star-forming haloes
and link them with the observed properties of Galactic halo stars in
order to constrain the nature of Pop III stars and of the first SNe.

## 2 THE MODEL

This section introduces our newly developed SAM NEFERTITI,
which in this work is coupled with a cosmological N-body simulation
for the formation of a MW analogue. NEFERTITI builds upon our
previous works with the SAM GAMETE (Salvadori et al. 2010;
Graziani et al. 2015, 2017; Pacucci, Loeb & Salvadori 2017; Pagnini
et al. 2023), but represents significant advancements, allowing us to
account for: (i) finite stellar lifetimes of both Pop III and Pop II/I
stars, i.e. relaxing the instantaneous recycling approximation (IRA);
(ii) the incomplete sampling of the IMF for both Pop III and Pop II/I
stars; and (iii) an unknown energy distribution function (EDF) for
Pop III stars exploding as SNe.

In the following, we describe briefly how the *N*-body simulation
traces the hierarchical growth of DM haloes (Section 2.1) and, in
detail, how the evolution of their baryonic content is followed by the
SAM (Section 2.2). Finally, in Section 2.3 we present the calibration
of the free parameters of our model.

## 2.1 The *N*-body simulation

We use a cold DM *N*-body simulation of a MW analogue (Scanna-
pieco et al. 2006; Salvadori et al. 2010) that has been carried out with
the GCD + code (Kawata & Gibson 2003a) using a multiresolution
technique (Kawata & Gibson 2003b). The highest resolution region
has a radius of four times the virial radius<fnref n="3" /> of the system, R<sub>vir</sub> =
239 kpc at z = 0, and a softening length of 540 pc. The system
comprises ~10⁶ DM particles of mass ~7.8 × 10⁵ M<sub>☉</sub>, i.e. a virial
mass M<sub>vir</sub> ≈ 7.8 × 10¹¹ M<sub>☉</sub>, consistent with observational estimates
for the MW (M<sub>vir,MW</sub> ≈ 6–25 × 10¹¹ M<sub>☉</sub>; see Wang et al. 2015
and references therein). A low-resolution simulation including gas
physics and star formation (SF) has been used to confirm that the
initial conditions will lead to the formation of a disc galaxy. The
position and velocities of all DM particles are stored, at each snapshot
of the simulation, and a friend-of-friends algorithm, with a linking
parameter b = 0.15 and a threshold number of particles of 50, is
used to identify the virialized DM haloes. The timestep between the
snapshots is Δt<sub>z</sub> ≈ 22 Myr at 17 < z < 8 and ≈ 110 Myr at z < 8.

## 2.2 The modelling of baryons

The SAM follows the flow of baryons from the IGM into the DM
haloes, the formation of stars and stellar evolution within each galaxy,
and the return of mass and metals into the interstellar medium
(ISM) and the IGM through stellar feedback. In order to resolve
the evolution of the most massive stars, we adopt a shorter sub-
timestep of δt<sub>s</sub> = 1 Myr for the SAM. At the end of each timestep
Δt<sub>z</sub> of the *N*-body simulation, the stellar and gas mass within each
halo is equally distributed among all its DM particles. The baryons
then follow the course of their respective DM particles in the next
integration step of the *N*-body simulation. This way we can extract the

spatial distribution of all stellar populations throughout the Galaxy's
assembly history.

At $z < 3$, when the DM halo of our central galaxy (i.e. the MW) has grown significantly, the assumption that the newly formed stars are equally distributed among its DM particles is no longer valid. Indeed, we know that a disc should form, leading to a more centrally confined SF. However, this approximation is suitable for investigating the ancient very metal-poor (VMP) stellar halo populations we are interested in, which form at $z > 5$ in smaller progenitor haloes. Indeed, we have confirmed that ancient stars (born before $z \sim 5$) inhabiting our mock MW's halo at $z = 0$ have been mostly acquired via mergers, i.e. they form *ex situ*.

### 2.2.1 Gas accretion

According to the traditional view of galaxy formation, all infalling
gas is initially shock heated to the virial temperature, T<sub>vir</sub>, of the
DM halo and forms a quasi-static atmosphere, which then cools
radiatively from the inside out and falls onto the central galaxy (hot
mode accretion; Rees & Ostriker 1977; Silk 1977; White & Rees
1978). More recent simulations have revealed a new paradigm in
which part of the gas enters the halo along cold, dense filaments
and accretes directly onto the galaxy without being shock-heated
(cold-mode accretion; Kereš et al. 2005; Dekel & Birnboim 2006;
Cattaneo et al. 2020). The latter is the dominant accretion mode for
virial masses M<sub>vir</sub> < 10<sup>11</sup> M<sub>☉</sub>, and is therefore more relevant for
our study that focuses on VMP stars ([Fe/H] < -2), which form at
z > 5 in our simulation, when the maximum halo mass is M<sub>vir</sub> ~
5 × 10<sup>10</sup> M<sub>☉</sub>.

In the lowest-mass haloes, baryonic infall may be substantially
reduced due to photodissociating and photoionizing radiation, since
gas cannot cool and accrete onto haloes with virial temperature, T<sub>vir</sub>,
lower than that of the IGM (e.g. Blanchard, Valls-Gabaud & Mamon
1992). To account for this effect, we assume that there is no gas
accretion onto haloes with T<sub>vir</sub> < T<sub>SF</sub>, where T<sub>SF</sub> = 2 × 10<sup>3</sup> K at z
> z<sub>rei</sub> = 6, and T<sub>SF</sub> = 2 × 10<sup>4</sup> K at z ≤ z<sub>rei</sub> when the Milky Way
environment is fully ionized (Salvadori et al. 2014).

For haloes with T<sub>vir</sub> > T<sub>SF</sub>, we assume that gas from the IGM
is continuously added into their cold filaments at a rate that is
proportional to their DM growth:

$$ \dot{M}_{\text{fil,accr}} = f_b \dot{M}_{\text{vir}}, \quad (1) $$

where $f_b = \Omega_b / \Omega_m$ is the universal baryon fraction.

The gas within the filaments is subsequently assumed to stream
onto the central galaxy on a free-fall time-scale:

$$ t_{\text{ff}} = \left( \frac{3\pi}{32G\rho} \right)^{1/2}, \quad (2) $$

where G is the gravitational constant and $\rho = \rho_{200}(z)$ is the total (dark + baryonic) mass density of the halo at redshift $z$. Hence, the gas accretion rate onto a galaxy and the variation of its filaments' gas mass are given, respectively, by

$$ \dot{M}_{\text{gas,accr}} = \frac{M_{\text{fil}}}{t_{\text{ff}}} \quad (3) $$

and

$$ \dot{M}_{\text{fil}} = \dot{M}_{\text{fil,accr}} - \dot{M}_{\text{gas,accr}}. \quad (4) $$

### 2.2.2 SF

SF occurs in a single burst at each sub-timestep $\delta t_s$, at a rate given
by the cold gas mass, $M_{\text{gas}}$, within the galaxy, the free-fall time

<footnote n="3">³For comparison, the virial radius of the MW is estimated at R<sub>vir</sub>, MW = 200–290 kpc (e.g. Dehnen, McLaughlin & Sachania 2006; Posti & Helmi 2019).</footnote>
(equation (2)), and the SF efficiency $\epsilon_{SF}$, which is a free parameter of our model:

$$SFR = \epsilon_{SF} \frac{M_{\text{gas}}}{t_{\text{ff}}}. \qquad (5)$$

Following Salvadori et al. (2015), we assume that the SF efficiency in minihaloes with $T_{vir} \le 2 \times 10^4$ K is reduced by a factor of $\epsilon_{H_2} = 2\epsilon_*[1 + (2 \times 10^4 K/T_{vir})^3]^{-1}$, to account for the ineffective cooling by molecular hydrogen (Salvadori & Ferrara 2012).

At each timestep $\delta t_s$, and for each halo, we compute the stellar mass formed, $M_* = SFR \cdot \delta t_s$, and form a Simple Stellar Population (SSP) only if $M_*$ is greater or equal to the maximum stellar mass $m_*^{\text{max}}$ allowed by the IMF.<fnref n="4" /> This way we ensure that stars throughout the whole mass range of the assumed IMF can be represented (see Section 2.2.3). Each SSP is characterized by its formation time $t_{\text{form}}$, the number and initial masses of its stars and their elemental abundances (which are equal to the ones of the gas in their host halo at $t_{\text{form}}$).

Our SF model is calibrated to act on the total gas content within a galaxy (see Section 2.3); we do not differentiate here between the different phases (e.g. cold, warm, molecular, and atomic) of the ISM. In addition, our model ignores some physical mechanisms, such as merger/instability induced starbursts, mass quenching (including AGN and halo quenching;<fnref n="5" /> Peng et al. 2010) and ram pressure (Gunn & Gott 1972), that can influence the SF rate of our central galaxy and that of its accreted satellite galaxies. Ram pressure can both trigger bursts of SF by compressing the gas within a satellite galaxy and quench its SF by stripping it (e.g. Kapferer et al. 2009; Bekki 2014; Koutsouridou & Cattaneo 2019). These effects might be important for the evolution of Local Group satellites at low $z$ (Simpson et al. 2018; Hausammann, Revaz & Jablonka 2019), but see also Salvadori et al. 2015 for a different view). However, at $z > 5$, the main progenitor (or major branch) of the MW with $M_{vir} < 10^{11} M_\odot$ (and $M_* \le 10^9 M_\odot$), is unlikely to hold a sufficiently massive hot atmosphere to exert high ram pressure (see, e.g. Gelli et al. 2020). The same is true for the AGN and halo quenching mechanisms, which come about at $M_{vir} \sim 10^{12} M_\odot$, or $M_* > 10^{10} M_\odot$ (e.g. Bluck et al. 2020; Cattaneo et al. 2020; Koutsouridou & Cattaneo 2022). Mergers and disc instabilities can drive inward gas flows that provoke bursts of SF in the central galactic regions (Barnes & Hernquist 1991; Teyssier, Chapon & Bournaud 2010; Zolotov et al. 2015; Bournaud 2016). Although the latter can be important for the evolution of high redshift gas-rich galaxies, they can be reasonably neglected for a study focused on the stellar halo such as ours.

### 2.2.3 The IMF of Pop III and Pop II/I stars

Following the critical metallicity scenario (e.g. Omukai 2000; Bromm et al. 2001; Schneider et al. 2002), we form Pop III stars when the total metallicity, $Z_{\text{gas}}$, of the gas within a progenitor halo is below a critical value $Z_{\text{crit}} = 10^{-4.5} Z_\odot$ (de Bennassuti et al. 2017), and PopII/I stars whenever $Z_{\text{gas}} \ge Z_{\text{crit}}$.

We adopt a Larson (1998) IMF for both Pop III and Pop II/I stars:

$$\phi(m_*) = \frac{dN}{dm_*} \propto m_*^{-2.35} \exp\left(-\frac{m_{\text{ch}}}{m_*}\right), \qquad (6)$$

but with very different characteristic mass, $m_{\text{ch}}$, and $m_*$ range.

For Pop II stars we assume $m_* = [0.1, 100] M_\odot$ and $m_{\text{ch}} = 0.35 M_\odot$, which is consistent with observations of present-day forming stars (Krumholz 2015).

For Pop III stars, we consider $m_* = [0.1, 1000] M_\odot$ and an IMF biased towards more massive stars, i.e. with $m_{\text{ch}} \ge 1 M_\odot$. This mass range and characteristic mass are indeed consistent with constraints on the Pop III IMF obtained from UFD galaxies (Rossi, Salvadori & Skúladóttir 2021) and in line with the results of cosmological hydrodynamical simulations for the formation of Pop III stars (Hirano et al. 2014; Susa, Hasegawa & Tominaga 2014; Hirano et al. 2015, see also the Introduction). In Pagnini et al. (2023), we showed that a $m_{\text{ch}} = 10 M_\odot$ is in good agreement with the observed [C/Fe] range within the bulge and can explain the dearth of CEMP-no stars with [C/Fe] > +1 in this environment. Therefore, we adopt $m_{\text{ch}} = 10 M_\odot$ as our starting point and explore the effect of different characteristic masses and a different maximum mass, $m_*^{\text{max}}$, for Pop III stars in Section 4.1.

In some cases, and most commonly in poorly star-forming low-mass haloes, the SF bursts are not strong enough to fully populate the theoretical IMF. This has important consequences for the type of stars formed, stellar feedback, chemical evolution and total stellar mass (Kroupa & Weidner 2003; Weidner & Kroupa 2006; Weidner et al. 2013; de Bennassuti et al. 2017; Applebaum et al. 2020; Rossi et al. 2021). To account for that, we implement in NEFERTITI a Monte Carlo procedure that generates a random sequence of stars, according to the assumed IMF, with total mass equal to the stellar mass formed in each SF burst (see Rossi et al. 2021 for details). In other words, we account for the incomplete sampling of the IMF of both Pop III and Pop II/I stars in poorly star-forming haloes.

### 2.2.4 The energy distribution of Pop III SNe

Currently, there is no theoretical constraint on the explosion energies of Pop III SNe with $m_* = 10-100 M_\odot$, and while observations suggest that they could have spanned almost 2 orders of magnitude, their distribution remains completely unknown (see the section 'Introduction'). Here, we account for the first time for the energy distribution of such primordial SNe in the context of a cosmological galaxy formation model. To this end, we assume a mass-independent EDF of the form:

$$\frac{dN}{dE_*} \propto E_*^{-\alpha_e}, \qquad (7)$$

where $E_*$ is the explosion energy and $\alpha_e$ is a free parameter of the model. Based on this underlying distribution, we assign randomly an energy level to each Pop III SN with $m_* = 10-100 M_\odot$. The top panel of Fig. 1 shows the cumulative probability as a function of $E_*$ for the different $\alpha_e$ values considered hereafter. The bottom panel shows the corresponding probability for a Pop III star to explode as a faint SN ($E_{51} = [0.3, 0.6]$), a core-collapse SN (ccSN; $E_{51} = [0.9, 1.2, 1.5]$), a high energy SN ($E_{51} = [1.8, 2.4, 3]$), or a hypernova ($E_{51} = [5, 10]$), for $\alpha_e = 0.5, 1$, and 2.

### 2.2.5 Relaxing the IRA

In all previous works where the predecessor of our SAM (the GAMETE SAM) was coupled with time consuming N-body simulations

<footnote n="4">⁴Throughout the text, the index (*) is used to refer to the total stellar mass of a galaxy, while the index (★) is used to refer to the mass of individual stars.</footnote>

<footnote n="5">⁵AGN quenching refers to feedback from an active galactic nucleus (AGN) that can eject cold gas from within a galaxy and/or prevent the hot gas surrounding it from cooling. Halo quenching refers to the disruption of cold filamentary flows by a massive hot atmosphere.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top: cumulative probability for a Pop III SN with $m_{\star} = 10-100 M_{\odot}$ to have a given explosion energy, $E_{51}$, for different values of the EDF (equation (7)) exponent, $\alpha_e$. Bottom: the corresponding probability of such star to explode as a faint SN, a ccSN, a high energy SN, or a hypernova, for $\alpha_e = 0.5$ (red), 1 (orange), and 2 (cyan).</figcaption></figure>

(Salvadori et al. 2010; Graziani et al. 2015, 2017; Pacucci et al. 2017; Pagnini et al. 2023), the chemical evolution was computed assuming the IRA, i.e. that all stars that do not survive until $z = 0$, die and return gas and metals into the ISM instantaneously. This approximation overestimates the ISM enrichment rate (steepening the metal-poor tail of the MDF) and describes poorly the abundance of elements produced on long time-scales, such as Fe and C. In addition, it blurs out the chemical signatures of different stellar types, such as primordial faint SNe and PISNe, that in reality explode at different times, by mingling their ejecta instantaneously.

In NEFERTITI, we abandon the IRA and, instead, follow the evolution of each individual star, depending on its initial mass and metallicity. At each timestep and for each halo, we compute the rate, $\dot{R}$, at which gas is restored into the ISM through stellar winds or SN explosions from:

$$ \dot{R} = \int_{m_{\text{turn}}(t)}^{m_{\star}^{\text{max}}} (m_{\star} - w_m(m_{\star})) N(\bar{t}_{\text{form}}, m_{\star}) dm_{\star}, \quad (8) $$

where $N(\bar{t}_{\text{form}}, m_{\star})$ is the number of stars with mass $m_{\star}$ that were formed at time $\bar{t}_{\text{form}} = t - \tau_{\star}$, $w_m$ and $\tau_{\star}$ are the remnant mass and lifetime of a star with initial mass $m_{\star}$, and $m_{\text{turn}}(t)$ is the turn-off mass, i.e. the mass corresponding to $\tau_{\star} = t$.

Similarly, we define the total ejection rate, $\dot{Y}_i$, of an element i that is returned to the ISM without being re-processed (first term in the square brackets) and newly produced (second term):

$$ \dot{Y}_i = \int_{m_{\text{turn}}(t)}^{m_{\star}^{\text{max}}} [(m_{\star} - w_m(m_{\star}) - m_i(m_{\star}, Z_{\star}))Z_i(\bar{t}_{\text{form}}) + m_i(m_{\star}, Z_{\star})] N(\bar{t}_{\text{form}}, m_{\star}) dm_{\star}, \quad (9) $$

where $m_i(m_{\star}, Z_{\star})$ is the mass of element i that is synthesized by a star with initial mass $m_{\star}$ and metallicity $Z_{\star}$, and $Z_i(\bar{t}_{\text{form}})$ is the mass fraction of the element i at the time of formation of each star.

We adopt the stellar lifetimes of Raiteri, Villata & Navarro (1996) for Pop II/I stars and the stellar lifetimes of Schaerer (2002) for Pop III stars.

### 2.2.6 Stellar yields and mixing

The metal yields and remnant masses of Pop III stars, entering equations (8) and (9), are adopted from Heger & Woosley (2002) for PISNe (140 $\le m_{\star}/M_{\odot} \le 260$) and from Heger & Woosley (2010) for less massive Pop III SNe (10 $\le m_{\star}/M_{\odot} \le 100$). The latter are given for 10 different explosion energies in the range 0.3–10 $\times 10^{51}$ erg (see Section 2.2.4). Heger & Woosley (2010) use a 1D code that cannot capture mixing between stellar layers due to its multidimensional nature. They, therefore, implement mixing artificially by moving a running boxcar through the star typically four times (Welsh, Cooke & Fumagalli 2021). For each stellar mass and explosion energy, there are 14 different values for the mixing parameter (in the range $f_{\text{mix}} = 0-0.2512$), which is defined as the width of the boxcar in units of the helium core mass.<fnref n="6" /> The mixing parameter is unknown but it can largely affect the abundance of various chemical elements produced by Pop III SNe, such as carbon, e.g. Vanni et al. (in preparation). In the absence of theoretical yields for the intermediate pulsational PISNe, i.e. stars with $m_{\star} = 100-140 M_{\odot}$, we assume that they collapse into black holes returning no mass into the ISM.

For Pop II/I stars we adopt the yields of Limongi & Chieffi (2018; set R without rotation velocity) for massive stars evolving as core-collapse SNe (ccSNe) and the van den Hoek & Groenewegen (1997) yields for low and intermediate mass ($m_{\star} < 8 M_{\odot}$) AGB stars.

### 2.2.7 Mechanical feedback from SNe

Supernovae drive winds that, if sufficiently energetic, can escape the gravitational potential well of their host halo and expel gas and metals into the surrounding medium. Different kinds of SNe are characterized by different explosion energies. For Pop III PISNe, we adopt the mass-energy relation of Heger & Woosley (2002), while for metal-free stars with $m_{\star} = 10-100 M_{\odot}$ we consider all energy levels ($E_{51} = 0.3-10$) provided by Heger & Woosley (2010), independently of the stellar mass as explained in Section 2.2.4. For Pop II/I ccSNe, we assume an average explosion energy of $10^{51}$ erg. Therefore, at each timestep, the total power output from SNe in a halo is $\sum_i \dot{N}_{\text{SN}}^i \cdot E_{\text{SN}}^i$, where $\dot{N}_{\text{SN}}^i$ is the explosion rate of SNe with energy $E_{\text{SN}}^i$. If a fraction $\epsilon_{\text{wind}}$ of this power is converted into kinetic form, the gas outflow rate from the halo, $\dot{M}_{\text{gas,ej}}$, will satisfy:

$$ \frac{1}{2} \dot{M}_{\text{gas,ej}} u_{\text{esc}}^2 = \epsilon_{\text{wind}} \sum_i \dot{N}_{\text{SN}}^i E_{\text{SN}}^i, \quad (10) $$

where the wind efficiency, $\epsilon_{\text{wind}}$, is the second free parameter of our model and $u_{\text{esc}} = \sqrt{\frac{G M_{\text{vir}}}{r_{\text{vir}}}} = f(M_{\text{vir}}, z)$ (Barkana & Loeb 2001) is the escape speed of the halo.

<footnote n="6">⁶Rotation and mass loss at all stages of stellar evolution are ignored in the Heger & Woosley (2010) models.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Total stellar mass (top row) and gas-to-stellar mass ratio (bottom row) of our MW-analogue at $z=0$, as a function of the two free parameters of our model: the SF efficiency (for $\epsilon_{\text{wind}} = 0.002$; left-hand panel) and the wind efficiency (for $\epsilon_{\text{SF}} = 0.8$; right-hand panel). Grey shaded areas represent the observed global properties of the MW (see text for details).</figcaption></figure>

### 2.2.8 Following the gas evolution

At each sub-timestep $\delta t_s$ of the SAM, we compute the evolution of the total gas mass, $M_{\text{gas}}$ and mass of element $i$, $M_{Z_i}$ in the ISM of each galaxy from the equations:

$$ \dot{M}_{\text{gas}} = \dot{M}_{\text{gas,accr}} - \text{SFR} + \dot{R} - \dot{M}_{\text{gas,ej}} \quad (11) $$

and

$$ \dot{M}_{Z_i} = Z_i^{\text{IGM}} \dot{M}_{\text{gas,accr}} - Z_i \text{SFR} + \dot{Y}_i - Z_i \dot{M}_{\text{gas,ej}} \quad (12) $$

respectively, where $Z_i = M_{Z_i}/M_{\text{gas}}$ is the mass fraction of element $i$ in the ISM and $Z_i^{\text{IGM}}$ is its mass fraction in the IGM. The latter is updated after each sub-timestep by summing the contributions of all haloes $h$:

$$ \dot{M}_{Z_i}^{\text{IGM}} = \sum_h \left( -Z_i^{\text{IGM}} \dot{M}_{\text{gas,accr}}^h + Z_i^h \dot{M}_{\text{gas,ej}}^h \right) \quad (13) $$

Equations (12) and 13 imply that the ejecta of dying stars are instantaneously and homogeneously mixed within both the ISM of each halo and the IGM, which are thus characterized by a unique chemical composition at each timestep. This perfect mixing approximation has significant consequences for the chemical evolution of our system, which are discussed in Sections 4 and 6.

## 2.3 Model calibration

In Fig. 2, we compare the observed global properties of the MW (grey shaded areas) with the results of our model at $z=0$, as a function of the SF efficiency $\epsilon_{\text{SF}}$ (left-hand panel), and the wind efficiency $\epsilon_{\text{wind}}$ (right-hand panels). One can see that the final stellar mass, $M_*$, depends strongly on $\epsilon_{\text{wind}}$ but only weakly on $\epsilon_{\text{SF}}$, while the opposite is true for the gas-to-stellar mass ratio, $M_{\text{gas}}/M_*$. This can also be inferred by solving analytically the equations that govern the evolution of the stellar and gas mass in a galaxy (see equations (4) and (6) in Koutsouridou & Cattaneo 2019). Therefore, these two observables are sufficient to constrain the two free parameters of our model.

We adopt $\epsilon_{\text{SF}} = 0.8$, and $\epsilon_{\text{wind}} = 0.002$, for which we obtain $M_* \approx 4.2 \times 10^{10} M_{\odot}$, and $M_{\text{gas}}/M_* \approx 0.13$ at $z=0$. Bland-Hawthorn & Gerhard (2016) report $M_* = (5.4 \pm 0.57) \times 10^{10} M_{\odot}$, $M_{\text{vir}} = (1.3 \pm 0.3) \times 10^{12} M_{\odot}$; and $M_* = 5.04^{+0.43}_{-0.52} \times 10^{10} M_{\odot}$, $M_{\text{vir}} = 1.08^{+0.20}_{-0.14} \times 10^{12} M_{\odot}$, respectively.<fnref n="7" /> We choose a value for $\epsilon_{\text{wind}}$ that gives a stellar mass at the lower limit of the observational range, since the virial mass of our system is lower than most observational estimates (see Section 2.1). By substituting $\epsilon_{\text{wind}} = 0.002$ in equation (10), we get a mass loading factor<fnref n="8" /> $\eta \equiv \dot{M}_{\text{gas,ej}}/\text{SFR} \approx 0.14$ for our MW-analogue at $z=0$, in agreement with recent observational estimates ($\eta_{\text{MW}} = 0.1 \pm 0.06$, Fox et al. 2019).

Our adopted value for $\epsilon_{\text{SF}}$ results in a gas-to-stellar mass ratio within the typically reported range of 0.1–0.15 (Ferrière 2001; Stahler & Palla 2004). In addition, at $z=0$, $\epsilon_{\text{SF}}$ corresponds to a SF time-scale $t_{\text{SF}} \equiv M_{\text{gas}}/\text{SFR} \sim 1.9$ Gyr, and a SFR $\sim 2.7 M_{\odot} \text{yr}^{-1}$ for our MW-analogue, in agreement with observational estimates ($t_{\text{SF,MW}} = 2$ Gyr, Bigiel et al. 2008; $\text{SFR}_{\text{MW}} = 1-3 M_{\odot} \text{yr}^{-1}$, Chomiuk & Povich 2011; Bland-Hawthorn & Gerhard 2016). At $z=0$, the mean metallicity of the ISM, $Z_{\text{ISM}} = 1.07 Z_{\odot}$, in our MW-analogue and of the IGM, $Z_{\text{IGM}} \sim 0.17 Z_{\odot}$, are in accordance with observations in the Galactic disc and in high-velocity clouds currently accreting onto the disc ($\sim 0.1-0.3 Z_{\odot}$; Tripp et al. 2003; Ganguly et al. 2005; Danforth & Shull 2008). For all the above reasons, we are confident that our model with the selected values for $\epsilon_{\text{SF}}$ and $\epsilon_{\text{wind}}$, is a good representation of the evolution of a MW-like galaxy and its immediate environment.

## 3 STELLAR DATA FOR MODEL COMPARISON

This section describes in detail the available observations of Galactic halo stars that we use to compare with key results of our model— namely the predicted metallicity distribution function (MDF), the fraction of CEMP stars, the carbonicity distribution function (CDF), and the distribution of stars in the [C/Fe]–[Fe/H] diagram.

(i) The MDF: for Galactic halo stars with $-4 < [\text{Fe/H}] < -2$ we adopt the MDF proposed by Bonifacio et al. (2021), which is the largest and most complete (i.e. biased-corrected) MDF for model comparison. It represents the average of three independently derived MDFs, the one by Naidu et al. (2020; H3 Survey), the uncorrected one by Schörck et al. (2009; Hamburg/ESO Survey), and the one determined by Bonifacio et al. (2021) themselves, and corrected for selection biases, from Sloan Digital Sky Survey spectra. The standard deviation in each metallicity bin (shown with the black errorbars in, e.g. the top row of Fig. 3) is provided by Bonifacio et al. (2021) as an error estimate on the MDF. The three MDFs used for computing the average are essentially identical above $[\text{Fe/H}] \sim -3$, explaining the small errors in this metallicity range. However, there are other published Galactic halo MDFs, that appear steeper (e.g. the corrected MDF by Schörck et al. 2009, and the one of Carollo et al.

<footnote n="7">⁷Studies estimating the MW stellar mass from direct integration of starlight usually find higher values, $M_* \sim 6 \times 10^{10} M_{\odot}$ (Bland-Hawthorn & Gerhard 2016), but these studies do not simultaneously constrain its virial mass.</footnote>

<footnote n="8">⁸Note that $u_{\text{esc}}^2 \approx 2.8 \times 10^{47} \text{erg}/M_{\odot}$ for $M_{\text{vir}} = 7.8 \times 10^{11} M_{\odot}$ at $z=0$ (equation (27) in Barkana & Loeb 2001) and that with the assumed Pop II/IMF we form one SN every $100 M_{\odot}$.</footnote>
2010) or shallower (e.g. Youkim et al. 2020) in this [Fe/H] range. At lower [Fe/H], the Naidu et al. (2020) MDF, which is based on high-resolution data, is undefined. The Schörck et al. (2009) MDF and the one determined by Bonifacio et al. (2021) extend down to [Fe/H] ~ -4, but are based on low resolution (~2000) data. Due to the fact that metallicities of [Fe/H] < -3 can only be accurately and precisely determined through high-resolution spectra, and due to the low number statistics at low [Fe/H], the average MDF shows much larger errors at ~4 ≤ [Fe/H] ≤ -3.

(ii) *The lowest-Fe tail of the MDF:* we compute the halo MDF at [Fe/H] < -4 using data from the SAGA<fnref n="9" /> database (2023, April 10 version), which assembles the abundances of all [Fe/H] ≤ -2.5 stars derived from high and medium-resolution follow up observations (Suda et al. 2008, 2011; Yamada et al. 2013). As researchers usually favour detailed studies of the most extreme stars, follow-up observations are biased towards low metallicities. For this reason, we only consider the SAGA MDF at [Fe/H] < -4 (including 43 stars), where we can assume that follow-up is near-complete. Still, it is likely that not all known stars with estimated -4.5 ≤ [Fe/H] ≤ -4 have been followed-up at high resolution. If those are confirmed in the future, the number of stars at [Fe/H] < -4.5 with respect to the number of stars at -4.5 ≤ [Fe/H] ≤ -4 will decline. We represent this possibility qualitatively with down-pointing arrows (shown, for example, in the second row of Fig. 4).

(iii) *The fraction of CEMP-no stars:* we compute the fraction of CEMP-no stars,

$$F_{\text{CEMP}}([\text{Fe/H}]) = \frac{N_{\text{CEMP}}([\text{Fe/H}])}{N_{*}([\text{Fe/H}])}, \quad (14)$$

where $N_{\text{CEMP}}$ is the number of CEMP-no stars, and $N_*$ the total at a given [Fe/H] bin, using the high/medium-resolution sample of Placco et al. (2014) and the high-resolution sample of Yong et al. (2013a). Placco et al. (2014) collected a large sample of VMP literature stars, excluded those with [Ba/Fe] > +0.6 and [Ba/Sr] > 0, which were likely enriched by an AGB companion (CEMP-s stars), and corrected the carbon abundances of the remaining sample to account for evolutionary effects. Their estimated fractions of CEMP-no ([C/Fe] ≥ +0.7) stars are shown, e.g. in the third row of Fig. 5. Yong et al. (2013a) performed a homogeneous chemical abundance analysis of 190 literature and program stars, of which 172 have [Fe/H] < -2. We completed their catalogue at low metallicities by adding the more recently discovered EMP stars shown in Fig. 7 (diamond points) and listed below. We computed $N_{\text{CEMP}}$ using the Yong et al. (2013b) classification of CEMP-no stars (that is based on the Aoki et al. 2007 criterion<fnref n="10" />) and excluding stars with [Ba/Fe] > +0.6. The latter criterion was adopted to be consistent with Placco et al. (2014). In both observational surveys the CEMP-no fraction decreases with increasing [Fe/H]. However, at ~4 ≤ [Fe/H] < -3, the $F_{\text{CEMP}}$ of Yong et al. (2013a) are significantly lower than those of Placco et al. (2014; see Fig. 5) highlighting the uncertainties in the observational estimates of $F_{\text{CEMP}}$ (see also Section 7).

(iv) *The CDF:* we compare our predictions for the [C/Fe] ('carbonicity', Carollo et al. 2012; Lee et al. 2017) distribution function of VMP inner halo stars, to observations from the SAGA database. We do so only at [C/Fe] > +2 where we expect the observational sample to have higher completeness. We compute the SAGA CDF by

including all confirmed CEMP-no stars, CEMP stars with upper limits for barium enhancement at [Ba/Fe] > +0.6 as well as CEMP stars with no measurement of barium abundances (35 stars in total with [C/Fe] > +2 and [Fe/H] ≤ -2, including upper limits for [C/Fe]).

(v) *The [C/Fe] versus [Fe/H] abundances:* in the left-hand panel of Fig. 7, we compare our predicted distribution of halo stars in the [C/Fe]-[Fe/H] diagram with the CEMP-no and C-normal stars in the catalogues of Yong et al. (2013a; X points) and Placco et al. (2014; circles) as well as with individual stars (diamond points) by Christlieb et al. (2002), Norris et al. (2007), Caffau et al. (2011), Hansen et al. (2014), Keller et al. (2014), Bonifacio et al. (2015), Frebel et al. (2015), Li et al. (2015), Caffau et al. (2016), Aguado et al. (2018), Bonifacio et al. (2018), François et al. (2018), Starkenburg et al. (2018), Aguado et al. (2019), Ezzeddine et al. (2019), and Nordlander et al. (2019).<fnref n="11" />

## 4 RESULTS

Using our model (Section 2), we investigate how the MDF, the CDF and the fraction of CEMP-no stars in the Galactic halo depends on the unknown EDF of the first Pop III supernovae. In Section 4.1, we examine the degeneracy between the EDF and the IMF of Pop III stars and in Section 4.2 we determine the key observables to constrain them.

Our findings are impacted by the stochastic sampling of both the masses and the SN explosion energies of Pop III stars (Section 2.2.3). Therefore, for each model presented in this article, we have averaged over 100 realizations as we find that this number is sufficient for our results to converge. We focus on the surviving VMP ([Fe/H] ≤ -2) stars that lie in the inner Galactic halo, i.e. at galactocentric radii 7 kpc ≤ $R_{\text{gal}}$ ≤ 20 kpc, at $z = 0$. Due to the instantaneous mixing of the IGM (that reaches [Fe/H] ≈ -2, by $z = 5$), we find that all VMP stars are formed before $z \sim 5$, regardless of the assumed properties (EDF, IMF, and mixing) of Pop III stars. Since we do not consider binary mass transfer, all CEMP (and C-normal) stars in our models, reflect the abundances of their birth clouds. We adopt the solar abundances from Asplund et al. (2009).

### 4.1 Pop III stars: EDF versus IMF

This section explores how the properties of the surviving inner halo stars in our model change when we vary the EDF of the first SNe. In what follows, we assume a flat distribution for the unknown Pop III mixing parameter and assign randomly one of the 14 mixing levels provided by Heger & Woosley (2010) to each Pop III star with $10 \le m_*/M_\odot \le 100$. We consider four Larson-type IMFs (equation (6)): (i) with characteristic mass $m_{\text{ch}} = 1 M_\odot$ and maximum mass for Pop III stars $m_*^{\text{max}} = 1000 M_\odot$, (ii) with $m_{\text{ch}} = 10 M_\odot$ and $m_*^{\text{max}} = 1000 M_\odot$, (iii) with $m_{\text{ch}} = 100 M_\odot$ and $m_*^{\text{max}} = 1000 M_\odot$, and (iv) with $m_{\text{ch}} = 10 M_\odot$ and $m_*^{\text{max}} = 100 M_\odot$.

Fig. 3 shows the predicted MDFs for the four IMFs assuming different values for the EDF (equation (7)) exponent as denoted by the colour (see Fig. 1). We find that at [Fe/H] > -3, the halo MDF is essentially independent of the assumed Pop III IMF and EDF and is in almost perfect agreement with the one derived by Bonifacio et al. (2021). It is worth noting that the MDF is a genuine prediction of our model since no free parameter was tuned to reproduce it

<footnote n="11">¹¹The Norris et al. (2007) star and the Aguado et al. (2019)/Ezzeddine et al. (2019) star are already included in the sample of Yong et al. (2013a) and were, therefore, not added to it before computing $F_{\text{CEMP}}$.</footnote>

<footnote n="9">⁹http://sagadatabase.jp/</footnote>

<footnote n="10">¹⁰Both the Aoki et al. (2007) criterion for CEMP stars ([C/Fe] ≥ +0.70, for $\log(L/L_\odot) \le 2.3$ and $[C/Fe] \ge +3.0 - \log(L/L_\odot)$, for $\log(L/L_\odot) > 2.3$) and the corrections of Placco et al. (2014), account for the depletion of the surface carbon abundance of stars as they ascent the red-giant branch.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Mean MDFs of VMP inner halo stars, normalized at [Fe/H] = -2.05, in comparison to the observations by Bonifacio et al. (2021; points with errorbars). Four Larson-type (equation (6)) IMFs for Pop III stars are considered: the first three with mass range $m_{\star} = 0.1-1000$ M$_{\odot}$ and characteristic mass $m_{\text{ch}} = 1$ M$_{\odot}$ (top left-hand panel), $m_{\text{ch}} = 10$ M$_{\odot}$ (top right-hand panel), $m_{\text{ch}} = 100$ M$_{\odot}$ (bottom left-hand panel), and the fourth with $m_{\star} = 0.1-100$ M$_{\odot}$ and $m_{\text{ch}} = 10$ M$_{\odot}$ (bottom right-hand panel). For each IMF, results are shown for different values of the Pop III EDF exponent (see equation (7)), $\alpha_e$, as denoted by the colour. All mixing values of the Heger & Woosley 2010 yields are assumed to be equally probable.</figcaption></figure>

(Section 2.3). At [Fe/H] < -3, the MDFs resulting from different models start to differentiate and a clear trend emerges: the MDF steepens as we move both to a higher characteristic mass $m_{\text{ch}}$ (or higher $m_{\star}^{\text{max}}$) at fixed EDF or to a lower $\alpha_e$ at fixed IMF, i.e. as more Pop III stars with high masses and explosion energies are formed. This is to be expected, since high energy SNe, hypernovae, and massive PISNe yield more iron than faint and ccSNe and, thus, accelerate the chemical enrichment in their host haloes resulting in a steeper MDF. Never the less, all MDFs are in agreement with the Bonifacio et al. (2021) one within the (large) errorbars.

In addition, all model MDFs are in agreement with the SAGA MDF at [Fe/H] < -4 within errors (Fig. 4). Yet, the ones predicted by the [$m_{\text{ch}}, m_{\star}^{\text{max}}$] = [1, 1000] M$_{\odot}$ and the [$m_{\text{ch}}, m_{\star}^{\text{max}}$] = [10, 100] M$_{\odot}$ models, which result in less (or zero) PISNe formed, lie above the observational data points in most metallicity bins. Instead the ones with $m_{\text{ch}} = 10$ M$_{\odot}$ and $m_{\text{ch}} = 100$ M$_{\odot}$ (and $m_{\star}^{\text{max}} = 1000$ M$_{\odot}$) lie on top and below the observations, respectively. As explained in Section 3, it is probable that the number of stars at [Fe/H] < -4.5 with respect to the number of stars at [Fe/H] ≈ -4 will decline in the future as more stars in the range -4.5 ≤ [Fe/H] ≤ -4 will be followed up with high resolution observations. Therefore, we could say that the latter two models are preferable here, or else that PISNe are required to steepen the ultra metal poor tail of the MDF. However, this could be a premature conclusion, given that only 14 stars with [Fe/H] < -4.5 have been observed to date.

In Fig. 5, we compare our predicted CEMP fractions to the observations of Yong et al. (2013a) and Placco et al. (2014; see Section 3). All halo stars with [Fe/H] < -5 are predicted to be carbon enhanced, [C/Fe] > +0.7. At higher metallicities and for a

given Pop III IMF, we find that EDFs skewed towards high explosion energies (i.e. with smaller $\alpha_e$, see Fig. 1) result in lower CEMP fractions; naturally, since more energetic Pop III SNe yield less [C/Fe] at fixed mass and mixing level (e.g. see Vanni et al., in preparation).

The dependence of the yielded [C/Fe] on the Pop III stellar mass is not straightforward. At $10 \le m_{\star}/M_{\odot} \le 100$ and at a given explosion energy, the ejected [C/Fe] appears to increase with mass, especially at low mixing levels, but the relation is not monotonous and tends to reverse at the highest explosion energy $E_{\star} = 10 \times 10^{51}$ erg (Heger & Woosley 2010). The opposite is true for PISNe; the yielded [C/Fe] decreases dramatically with stellar mass from $\sim 10^{13}$ at $m_{\star} = 140$ M$_{\odot}$ to $\le 10^{-1}$ at $m_{\star} = 260$ M$_{\odot}$ (Heger & Woosley 2002). In addition, only the lowest explosion energies ($E_{51} \le 1.5$) and the lowest mixing levels produce [C/Fe] that exceed those of the least massive PISNe ($m_{\star} \lesssim 170$ M$_{\odot}$). Yet, all non-PISNe yield higher [C/Fe] than PISNe with $m_{\star} \gtrsim 195$ M$_{\odot}$.

Never the less, Fig. 5 reveals a clear trend. As we increase the characteristic mass from $m_{\text{ch}} = 1$ M$_{\odot}$, to $m_{\text{ch}} = 10$ M$_{\odot}$, and $m_{\text{ch}} = 100$ M$_{\odot}$ (resulting in $M_{\text{PISN}}/M_{\text{PopIII}} \approx 0.04, 0.11$, and $0.22$, respectively) the predicted CEMP fraction for a given EDF decreases. For $m_{\text{ch}} = 1$ M$_{\odot}$ we can reproduce the observed $F_{\text{CEMP}}$ for $\alpha_e > 1$ while for $m_{\text{ch}} = 10$ M$_{\odot}$ we need a higher $\alpha_e \gtrsim 1.5$. That means that as the number of PISNe increases, the number of hypernovae should drop (from < 22 per cent for $m_{\text{ch}} = 1$ M$_{\odot}$ to < 8 per cent for $m_{\text{ch}} = 10$ M$_{\odot}$). For $m_{\text{ch}} = 100$ M$_{\odot}$ even the EDF with $\alpha_e = 4$ (giving $\sim 99$ per cent faint SNe) cannot produce enough CEMP stars to meet the observations. This is because when $m_{\text{ch}} = 100$ M$_{\odot}$, PISNe dominate the ISM enrichment thus washing out the high [C/Fe]
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Fig. 3 but normalized to [Fe/H] = −4.25 and compared to the ultra metal poor MDF from the SAGA database. The downwards pointing arrows indicate the possibility that the SAGA MDF might be biased towards the lowest metallicities (see Section 3).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Differential CEMP fractions of very metal poor inner halo stars, for the four IMFs and different $\alpha_e$ values considered (see Section 4.1 or caption of Fig. 3). Data points show the observational estimates of Yong et al. (2013a; points with errorbars) and Placco et al. (2021; stars) as in the third row of Fig. 9.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cumulative [C/Fe] distribution functions of very metal poor, inner halo stars, normalized at [C/Fe] = +2, in comparison with observations from the SAGA database.</figcaption></figure>

yielded by faint SNe (Pagnini et al. 2023). Indeed, in the case where $m_{\star}^{\text{max}} = 100 \text{ M}_{\odot}$, i.e. when no PISNe are allowed to form, our model fits the observations for $\alpha_e \sim 1.5-2$.

Finally, in Fig. 6, we compare the cumulative CDFs, for stars with $[\text{Fe/H}] \le -2$, predicted by our models with observations from the SAGA database. We find that the CDFs become steeper as the $m_{\text{ch}}$ of the Pop III IMF increases, or else as more PISNe form. Models with $m_{\text{ch}} = 100 \text{ M}_{\odot}$ (bottom left-hand panel) significantly underpredict the number of stars with $[\text{C/Fe}] > +4$, yet all other models are in agreement with the observations within errors. At fixed IMF, we see that when the model CDFs are normalized to $[\text{C/Fe}] = +2$, they show no clear dependence on the assumed EDF. That is not true, however, when we consider the CDFs extending down to lower $[\text{C/Fe}]$. There, our now-familiar trend is evident; the higher the energy of Pop III SNe (or else the lower the $\alpha_e$), the lower the yielded $[\text{C/Fe}]$ and, therefore, the steeper the resulting CDF (see Appendix A).

## 4.2 Metal contribution from Pop III stars

In the previous section, it became clear that there exist degeneracies between the EDF and the IMF of Pop III stars. Furthermore, given the currently large errors in the observational data, it is hard to single out a preferred model. Therefore, it is useful to examine the average behaviour of all our ‘successful’ models, see Table 1, i.e. those that are in better agreement with the observed MDF, CDF and $F_{\text{CEMP}}$.

Fig. 7 (left-hand panel) shows the present-day distribution of VMP inner halo stars in the [C/Fe]–[Fe/H] diagram, averaged among all models of Table 1. Notice that the slope of the [C/Fe]–[Fe/H] relation and the region showing the highest probability (i.e. that of C-normal stars) do not strongly depend on the choice of model. However, the exact value of $P_i$ in each [C/Fe]–[Fe/H] bin is model-dependent. For example, models with a higher number of low energy SNe (higher

$\alpha_e$) and/or a lower characteristic mass of the IMF, produce more stars, i.e. a higher probability, $P_i$, at $[\text{Fe/H}] \lesssim -4$ (see Fig. 3). Therefore, the average probability shown here is only a rough estimate, as $P_i$ varies across models and there is no reason to assume that each of the models in Table 1 should have equal weight in the calculation of the average.

The main bulk of the observed C-normal population is in very good agreement with our model predictions, and coincides with the region predicted to have the highest density of stars. Furthermore, the sparser CEMP-no stars are also well represented by our models. Similar to the observations, our models show a sharply decreasing $[\text{C/Fe}]$ with increasing $[\text{Fe/H}]$. However, our $[\text{C/Fe}]$–$[\text{Fe/H}]$ relation appears shifted towards lower $[\text{Fe/H}]$ compared to the observed one. As a result the CEMP stars with the highest carbonicities in each $[\text{Fe/H}]$ bin are not reproduced by our models. This is a problem faced by several other works (e.g. Cooke & Madau 2014; Komiya et al. 2020; Jeon et al. 2021). In Section 7, we discuss possible solutions to this discrepancy coming both from the modelling and from the observational side.

The right-hand panel of Fig. 7, depicts the minimum metal fraction contributed by Pop III stars, as a function of metallicity and carbon enhancement. In particular, the colours denote the minimum $f_{Z}^{\text{Pop III}} \equiv m_{Z}^{\text{Pop III}}/m_{Z}^{\text{tot}}$ of all stars belonging to each $[\text{C/Fe}]$–$[\text{Fe/H}]$ bin in our models (Table 1), where $m_{Z}^{\text{tot}} = m_{Z}^{\text{Pop III}} + m_{Z}^{\text{Pop II}}$ is the total mass of metals in a star and $m_{Z}^{\text{Pop III}}$ and $m_{Z}^{\text{Pop II}}$ the metals’ mass that it has inherited from Pop III and Pop II progenitors, respectively. We find that all C-enhanced stars at $[\text{Fe/H}] < -3$, are at least $\sim 20$ per cent enriched by Pop III progenitors.<fnref n="12" /> As we go towards higher $[\text{C/Fe}]$, this value increases rapidly, to $> 50$ per cent

<footnote n="12">¹²Note that our model does not include binary transfer, i.e. CEMP-s stars.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left-hand panel: distribution of VMP, inner halo stars on the [Fe/H]–[C/Fe] diagram. The colour denotes the probability $P_i \equiv N_i^i/N_i^{tot}$ of stars to belong in each bin $i$ (averaged among all models of Table 1), where $N_i^i$ is the number of stars in the bin and $N_i^{tot}$ is the total number of [Fe/H] $\le -2$ stars. Right-hand panel: same as the left-hand panel only here the colour denotes the minimum metal fraction inherited by Pop III ancestors for all the stars in each [Fe/H]–[C/Fe] bin. Data points in both panels, show the C-normal and CEMP-no stars from the samples of Yong et al. (2013a; × points), Placco et al. (2014; points), and various authors (diamonds; see Section 3).</figcaption></figure>

<table><caption>**Table 1.** Pop III stellar parameters for the models that successfully reproduce the observed MDF (Bonifacio et al. 2021), the CDF, and the CEMP fractions (Yong et al. 2013a; Placco et al. 2014) of VMP stars in the inner halo. In all models below, all values of stellar mixing given by Heger & Woosley (2010) are assumed to be equally probable.</caption><thead><tr><th></th><th></th><th>Model 1</th><th>Model 2</th><th>Model 3</th></tr></thead><tbody><tr><td rowspan="2">IMF</td><td><i>m</i><sub>ch</sub></td><td>1 M<sub>☉</sub></td><td>10 M<sub>☉</sub></td><td>10 M<sub>☉</sub></td></tr><tr><td><i>m</i><sub>*</sub><sup>max</sup></td><td>1000 M<sub>☉</sub></td><td>1000 M<sub>☉</sub></td><td>100 M<sub>☉</sub></td></tr><tr><td>EDF</td><td><i>α</i><sub>e</sub></td><td>1.0–2.0</td><td>1.5–2.5</td><td>1.5–2.0</td></tr></tbody></table>

for stars with [C/Fe] $> +1$, and to $> 80$ per cent for stars with [C/Fe] $> +1.5$ (at the same [Fe/H] $< -3$). Moreover, all stars with [C/Fe] $\gtrsim +2$ and/or [Fe/H] $\lesssim -4.7$ are pure Pop III descendants; their abundance patterns are less than 5 per cent contaminated by Pop II stars.

In addition to the Pop III enriched CEMP stars, there exists a group of C-enhanced stars that have been entirely enriched by Pop II stars, at [Fe/H] $> -2.8$ and [C/Fe] $< +2$ (dark blue area in Fig. 7). Our adopted Pop II SN yields (Limongi & Chieffi 2018) have a maximum [C/Fe] $= +0.69$ at [Fe/H] $\le -2$, and are, therefore, not able to beget VMP C-enhanced stars. Instead, we find that these (Pop II enriched) CEMP stars are descendants of Pop II AGB stars and can form in minihaloes only after the following conditions are met: (i) SN explosions expel all, or nearly all, gas from the halo, leaving none, or only a small fraction, of the iron they produced, (ii) subsequent accretion of pristine/metal-poor gas from the IGM leads to low [Fe/H] in the ISM, and (iii) previously formed AGB stars release carbon enriching the ISM to high [C/Fe] (see also Rossi et al. 2023).

The stars at $0 < [C/Fe] < +0.7$ are predominantly enriched by Pop II progenitors. They correspond to the highest density region of Fig. 7 (left-hand panel), which implies that most of the observed C-normal stars are not Pop III star descendants. Yet, at low $[C/Fe] < 0$, the Pop III metal contribution starts dominating again. This is a natural consequence of the fact that Pop II stars yield a minimum [C/Fe] $\approx 0.07$ at [Fe/H] $\le -2$ (Limongi & Chieffi 2018), while energetic Pop III SNe can reach down to [C/Fe] $\approx -1.3$ (Heger & Woosley 2002, 2010; see Section 6.1).

Fig. 8 shows the *average* contribution by different Pop III progenitors for stars in each [C/Fe]–[Fe/H] bin, for all models in Table 1 with: $m_{*}^{max} = 1000 M_{☉}$ (left-hand panel); and $m_{*}^{max} = 100 M_{☉}$, i.e. when no PISNe are allowed to form (right-hand panel). We find that at fixed [C/Fe], the surviving CEMP stars with the lowest [Fe/H] have been enriched mostly by faint SNe. Instead, PISNe enrichment dominates the pollution of the most [Fe/H]-rich CEMP stars. Notice that yields from PISNe with $m_{*} = (140-150) M_{☉}$ can reach significantly higher [C/Fe] than many non-PISNe, e.g. Heger & Woosley (2010) and Nomoto et al. (2013). In the absence of PISNe, however, all CEMP stars (except the Pop II AGB-descendants at [Fe/H] $> -2.8$) are on average $> 30$ per cent enriched by faint SNe. Compared to faint SNe and PISNe, the overall contribution of ccSNe, high energy SNe and hypernovae to the surviving stars is not as prominent. This is due to the fact that the EDFs of our preferred models (Table 1) are skewed towards low explosion energies—high $\alpha_e$ exponent—and hence produce much fewer SNe of high energies. None the less, there appears a region in the diagram (at $[C/Fe] < 0$ and $[Fe/H] \lesssim -2.5$) where stars are predominantly enriched by primordial hypernovae. For comparison, in Appendix A we show the results of a model in which all types of Pop III SNe are equally probable.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Distribution of VMP, inner halo stars in the [C/Fe]–[Fe/H] diagram, with the colour denoting the mean fraction of metals inherited from Pop III faint SNe, ccSNe, high energy SNe, hypernovae, and PISNe, based on the models of Table 1, with Pop III mass range $m_{\star} = 0.1–1000 M_{\odot}$ (left-hand panel); and $m_{\star} = 0.1–100 M_{\odot}$, i.e. when no PISNe are allowed to form (right-hand panel). Data points show the observations as in Fig. 7.</figcaption></figure>
The average contribution of any SN type, at a given [C/Fe]–[Fe/H] bin, varies depending on the assumed EDF. However, the qualitative trends described above do not: we find that long-lived descendants of each kind of Pop III SNe always occupy the same regions on the [C/Fe]–[Fe/H] diagram. In particular, certain [C/Fe]–[Fe/H] combinations can only be produced by an enrichment of a specific type of Pop III SN. Hypernovae descendants predominantly populate a well-defined region at [C/Fe] < 0 and [Fe/H] ≲ −2.5, while the [C/Fe] ≲ −0.5 area is only populated with PISNe descendants. Thus, without hypernovae or PISNe these areas are not represented in our models (compare, e.g. the bottom panels in Fig. 8).

# 5 THE IMPACT OF STELLAR MIXING

The convective mixing between stellar layers can influence quite strongly the chemical signature of Pop III SNe. When mixing precedes fallback, heavier nuclei that would not have been ejected otherwise can escape into the ISM. In the previous section, we adopted a uniform distribution for the mixing parameter, $f_{\text{mix}}$ of the Heger & Woosley (2010) yields (see Section 2.2.6). Here, we explore how varying $f_{\text{mix}}$ impacts our model’s predictions.

The left- and right-hand columns of Fig. 9, show our results as a function of $f_{\text{mix}}$, assuming an EDF given by equation (7) with $\alpha_e = 1$ and $\alpha_e = 2$, respectively. Both cases assume a Pop III IMF of the form of Larson (1998; equation (6)) with $m_{\text{ch}} = 10 M_{\odot}$ and $m_{\star}^{\text{max}} = 1000 M_{\odot}$, i.e. the second IMF considered in Section 4.1.

Similar to the case of varying the Pop III EDF and IMF (Section 4.1), the predicted MDFs at [Fe/H] > −3 remain unaffected by the assumed stellar mixing and in agreement with the one by Bonifacio et al. (2021; first row of Fig. 9). At [Fe/H] < −3, the MDFs resulting from different models start to differentiate and this difference becomes pronounced at [Fe/H] ≤ −4; for a given EDF, lower mixing levels produce flatter MDFs. The comparison with the SAGA MDF is crucial to potentially discard models (second row of Fig. 9). Models with $f_{\text{mix}} \le 0.0631$ produce MDFs that lie above the observational values in several metallicity bins. Still we cannot completely discard them since they are in agreement with the observations within errors. On the other hand, higher mixing levels ($f_{\text{mix}} \ge 0.1$) result in more heavy elements, such as iron, being released into the ISM. The chemical enrichment in the first minihaloes proceeds faster giving rise to steeper MDFs that underpredict the number of hyper metal poor stars by more than 1–2 dex.

The third row of Fig. 9 shows the fraction of CEMP-no stars, $F_{\text{CEMP}}$ (equation (14)), in the inner halo, as predicted by the different models. We find that mixing levels $f_{\text{mix}} \le 0.1$ yield almost identical $F_{\text{CEMP}}$ in both models. Those lie well below the observations in the $\alpha_e = 1$ model, but are in agreement with the observations in the $\alpha_e = 2$ model, due to the lower number of high energy SNe and hypernovae there (7 per cent for $\alpha_e = 2$ compared to 43 per cent for $\alpha_e = 1$). The two highest mixing levels, $f_{\text{mix}} = 0.1855$ and 0.2512, underestimate the CEMP fractions in both models and do not produce any stars at [Fe/H] < −6.5. We also notice that the $F_{\text{CEMP}}$ dependence on mixing is weaker for the $\alpha_e = 1$ model, which generates more energetic SNe. Indeed, the higher the explosion energy, or equivalently the smaller the fallback, the weaker the effect of mixing on stellar ejecta<fnref n="13" /> (see Fig. 10).

The $F_{\text{CEMP}}$–[Fe/H] relation can convey only limited information, since CEMP is a binary classification—a star either has [C/Fe] > +0.7 or not. The CDF, instead, can be more informative. The bottom row of Fig. 9 shows the cumulative CDF for all inner halo stars with [Fe/H] < −2 and [C/Fe] > +2 in each model. Here, the effect of mixing is pronounced in both the $\alpha_e = 1$ and the $\alpha_e = 2$ model. In both cases, models with $f_{\text{mix}} \ge 0.0631$ predict too few stars with high carbonicities ([C/Fe] ≳ +3). Instead, the CDFs for lower $f_{\text{mix}}$ are in good agreement with the observations for $\alpha_e = 1$ and even better for $\alpha_e = 2$.

Regardless of the degeneracy between the stellar mixing and the explosion energy in the yielded abundance ratios—increasing either the mixing or the explosion energy yields lower [C/Fe]—we find that when the fraction of energetic SNe is significant (e.g. when $\alpha_e = 1$) even the lowest mixing level is inconsistent with the observations. Similarly, the highest mixing levels fail to reproduce the observations even when the EDF is dominated by faint SNe. More specifically, inspection of the second and fourth row of Fig. 9 reveals that a typical $f_{\text{mix}} \le 0.0631$ is favoured by our model. This result is not at odds with our adoption of a uniform mixing distribution in Sections 4.1 and 4.2, since only 3 out of the 14 mixing levels of (Heger & Woosley 2010) are above $f_{\text{mix}} = 0.0631$. We should note that Heger & Woosley (2010) reached similar conclusions. In particular, they found that regardless of the assumed Pop III IMF and explosion energy, an $f_{\text{mix}} \le 0.0631$ provides the best fit to the abundance patterns of the C-enhanced stars HE1327-2326 (Frebel et al. 2005) and HE0107-5240 (Christlieb et al. 2002).

It would be valuable to generalize this result to other stellar evolutionary models. However, this endeavour is no easy feat, first, due to the fact that the mixing prescription in Heger & Woosley (2010) is artificial and not described using physical principles, and second, due to the many differences in the physical and numerical assumptions employed by different groups. Never the less, one can get a crude idea by comparing the yielded abundance ratios obtained in different studies to their adopted parameter values.

Fig. 10 compares the [C/Fe] yields provided by Heger & Woosley (2010) for Pop III SNe as a function of mixing and explosion energy, to those computed for Pop III SNe by Iwamoto et al. (2005), Tominaga et al. (2007), Marassi et al. (2014), and Limongi & Chieffi (2012). In the first three works, mixing is a free parameter like in the Heger & Woosley (2010) models but is parametrized in terms of a minimum and maximum mass coordinates within which the mass of each element is uniformly mixed. Those two together with the ‘mass cut’ parameter, i.e. the mass coordinate below which all material falls back onto the central remnant, are calibrated<fnref n="14" /> to reproduce the abundance pattern of the stars: (i) HE1327-2326 and HE0107-5240 in the case of Iwamoto et al. (2005), (ii) the CEMP star SMSSJ031300 observed by Keller et al. (2014) in the case of Marassi et al. (2014), and (iii) the average abundance pattern of four C-normal EMP stars Cayrel et al. (2004) in Tominaga et al. (2007). Naturally, a more extended mixing region at fixed mass cut, or a lower mass cut at fixed mixing, would result in lower yielded [C/Fe]. In Limongi & Chieffi (2012), stellar mixing is not artificial but it is coupled to nuclear burning and only the mass cut (or equivalently the explosion energy) is calibrated to reproduce the average abundance pattern of the C-normal stellar sample of Cayrel et al. (2004).

Since the free parameters in the above models have been fitted to reproduce specific stars, they give a quasi-constant [C/Fe], regardless of the assumed stellar mass; the Iwamoto et al. (2005) and Marassi

<footnote n="13">¹³In the limit where the fallback is zero, the stellar yields are independent of mixing.</footnote>

<footnote n="14">¹⁴One can find the adopted parameter values in the respective papers.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Predicted MDF of inner halo stars (7kpc < R<sub>gal</sub> < 20kpc), normalized at [Fe/H] = -2.05 (top row), and [Fe/H] = -4.25 (second row), in comparison with observations from Bonifacio et al. (2021) and the SAGA database (see Section 3 for details). The third row shows the predicted fraction of CEMP stars ([C/Fe] ≥ +0.7) in each [Fe/H] bin, in comparison to the observations by Yong et al. (2013a; points with errorbars) and Placco et al. (2014; stars). The bottom row shows the cumulative [C/Fe] distribution function (CDF) of VMP inner halo stars, normalized at [C/Fe] = +2, in comparison with observations from the SAGA database. Colours in all panels denote the mixing parameter adopted for m<sub>*</sub> = 10–100 M<sub>☉</sub> Pop III stars in each run, as indicated by the colourbar. Lines and shaded areas represent the mean and standard deviation of 100 runs. The left- and right-hand panels show the model results assuming an EDF for Pop III SNe (equation (7)) with exponent α<sub>e</sub> = 1 and α<sub>e</sub> = 2, respectively. In both cases higher mixing for Pop III stars results in steeper MDFs and CDFs and lower CEMP fractions.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of the Heger & Woosley (2010) [C/Fe] yields for Pop III faint SNe (E<sub>51</sub> = 0.3–0.6; top left-hand panel), ccSNe (E<sub>51</sub> = 0.9–1.5; top right-hand panel), high energy SNe (E<sub>51</sub> = 1.8–3.0; bottom left-hand panel), and hypernovae (E<sub>51</sub> = 5–10; bottom right-hand panel) with the [C/Fe] yield of zero-metallicity stars by Iwamoto et al. (2005; stars), Tominaga, Umeda & Nomoto (2007; triangles), Marassi et al. (2014; circles), and Limongi & Chieffi (2012; X symbols). Colours denote the mixing level given by Heger & Woosley (2010), as indicated by the colourbar. For each mixing level, the top (bottom) colouroud solid line corresponds to the lowest (highest) explosion energy shown in each panel, with the area between the two lines shaded. Filled black data points show [C/Fe] yields of Pop III stars with explosion energies within (or close) to the range shown in each panel, while empty data points show SNe yields of different energies.</figcaption></figure>

et al. (2014) yields that have been calibrated to reproduce CEMP stars, lie above our typical $f_{\text{mix}} = 0.0631$ level in all energy bins, while the ones of Tominaga et al. (2007) and Limongi & Chieffi (2012) that have been calibrated to reproduce C-normal stars are consistent with the highest mixing levels and/or the highest explosion energies of Heger & Woosley (2010).

# 6 DISCUSSION

We have developed a new SAM of galaxy formation, named NEFERTITI (NEar-FiEld cosmology: Re-Tracing Invisible Times), and combined it with a cosmological N-body simulation of a MW analogue, to shed light on the properties of the first Pop III SNe, and in particular their EDF parametrized in our model with $\alpha_e$ (equation (7) and Fig. 1). NEFERTITI follows the formation and evolution of individual Pop III and Pop II stars, and considers, for the first time in a SAM for the MW formation, the contribution of Pop III SNe with different masses, stellar mixing, and explosion energies. Subsequently, we have investigated how varying these Pop III stellar parameters affects the properties of the surviving VMP stars in the Galactic halo.

## 6.1 Comparison with previous results and model limitations

For a given Pop III IMF, we find that a higher contribution of low energy SNe (higher $\alpha_e$) results in a greater CEMP-no fraction, and a flatter MDF and CDF for stars with [Fe/H] < −2. In particular, for a Larson (1998) IMF with characteristic mass $m_{\text{ch}} = 1 M_{\odot}$ and range $m_{\star} = [0.1, 1000] M_{\odot}$, we reproduce the observed CEMP fractions by Yong et al. (2013a) and Placco et al. (2021) if 40 per cent–80 per cent of Pop III stars with $10 \le m_{\star}/M_{\odot} \le 100$ explode as faint SNe and only 20 per cent–2 per cent as hypernovae (and the rest with intermediate energies; Fig. 1). When we increase $m_{\text{ch}}$ to $10 M_{\odot}$, the data are best fitted with 63 per cent–90 per cent of faint SNe, while for $m_{\text{ch}} = 100 M_{\odot}$ we always underpredict the fraction of CEMP stars, due to the high number of primordial PISNe (Fig. 5). The effects of the Pop III properties on the halo MDF are only prominent at [Fe/H] $\lesssim -3$, yet all models are in agreement with the current observations within errors (Figs 3 and 4). However, this comparison is somewhat hampered by limitations of the data (see the next subsection).

These results are not easily compared to previous works due to the different assumptions taken (e.g. the Pop III IMF), and the fact that, until now, none have considered simultaneously the contribution of
Pop III SNe of all different energies and Pop II stars. However, we report that Hartwig et al. (2018) find a best matching to the Placco et al. (2014) CEMP fractions by assuming that 40 per cent of Pop III stars explode as faint SNe and the rest as normal ccSNe, but without including high energy SNe, hypernovae and PISNe with $m_*$ > 150 M$_{\odot}$, all of which would lower $F_{\text{CEMP}}$. Contrarily, the simple minihalo model for Pop III star enrichment by Cooke & Madau 2014 manages to reproduce the observed $F_{\text{CEMP}}$ even when 100 per cent of Pop III SNe explode with high energy, when using a flat IMF in the range 10–100 M$_{\odot}$ (hence without including PISNe). This overestimation of the CEMP fraction relative to our model is to be expected, since Cooke & Madau (2014) do not consider the contribution from Pop II stars which, according to our simulation, dominates the chemical enrichment at [C/Fe] < +0.7 (see also Vanni et al. 2023).

One of our key findings is that all VMP stars with subsolar [C/Fe] are predominantly imprinted by Pop III hypernovae and/or PISNe, regardless of our model assumptions (Fig. 8). This result stems from the fact that our adopted Pop II metal yields (Limongi & Chieffi 2018) reach a minimum [C/Fe] = 0.07 at [Fe/H] ≤ −2, whereas Pop III hypernovae and PISNe can reach [C/Fe] as low as ≈ −0.9 and ≈ −1.3, respectively (Heger & Woosley 2002, 2010). The uncertainties associated with Pop II yields are crucial in this regard. Ritter et al. (2018) similarly estimate a yielded [C/Fe] > 0 for massive stars with [Fe/H] < −2, and Kobayashi et al. (2006) and Nomoto et al. (2013) find a yielded [C/Fe] < −0.1 only at [Fe/H] < −3.5 (see fig. 5 of Liang, Gjergo & Fan 2023). Yet, Woosley & Weaver (1995) suggest Pop II [C/Fe] yields reaching down to ≈ −0.5. Therefore, only if we adopted the Woosley & Weaver (1995) yields, we would anticipate a higher contribution of Pop II stars at subsolar [C/Fe]. It should be noted that our model does not include SN type Ia that yield [C/Fe] < −1 (Thielemann, Nomoto & Yokoi 1986; Iwamoto et al. 1999; Seitenzahl et al. 2013). However, we expect their contribution at such low metallicites to be minimal (e.g. Salvadori et al. 2015). Type Ia SNe have a typical delay time of 0.1–1 Gyr (see Chen, Hu & Wang 2021 and references therein), while the majority of our SF minihaloes reach [Fe/H] > −2 within 0.1 Gyr of their formation. In any case, SNe type Ia descendants can be distinguished from those of Pop III hypernovae and PISNe by comparing their complete abundance patterns (see, e.g. fig. 8 in Nomoto et al. 2013 and fig. 12 in Salvadori et al. 2019).

Besides the success of our model in reproducing the MDF, the CDF and the fraction of CEMP stars at [Fe/H] ≤ −2, we find that our predicted [C/Fe]–[Fe/H] relation lies at the lower side of the observations (Fig. 7). This discrepancy has been reported by several previous works, even though the [C/Fe]–[Fe/H] stellar distribution depends on the particularities of each assumed model.<fnref n="15" /> Cooke & Madau (2014), who investigate early chemical enrichment by Pop III stars in isolated minihaloes, find that they cannot reproduce the highest [C/Fe] observed in CEMP-no stars. Using a model calibrated on the UFD Böotes I, Rossi et al. (2023) find that true Pop III descendants have<fnref n="16" /> A(C) < 6, similar to us: the upper envelope of our [C/Fe]–[Fe/H] relation corresponds to A(C) ~ 6 – 6.5 at [Fe/H] ≤ −4.5. Yet, they predict the formation of Pop II AGB descendants with high C-abundances (A(C)~7–7.5) even at [Fe/H] < −4. As explained in Section 4.2, such CEMP stars only form after one or more SN explosions blow out ~all gas from within a halo. This process

removes the iron rich signature of Pop II SNe, allowing previously formed AGB stars to enrich the newly accreted, nearly pristine gas to high [C/Fe]. We find that this condition is satisfied in our model only at $z < 8$ when the IGM has already been enriched to [Fe/H] > −3 (Fig. 7).<fnref n="17" /> However, we must note that our DM simulation does not resolve haloes with $M_{\text{vir}} < 10^7$ M$_{\odot}$. Hydrodynamic simulations, instead, suggest that the minimum mass of the first SF minihaloes can range between $10^{5.5}$ and $10^{7.5}$ M$_{\odot}$ depending on the relative velocity between baryons and DM (see Schauer et al. 2019 and references therein). Including lighter minihaloes in our simulation could allow the formation of AGB-descendant CEMP stars at lower metallicities. In such case, one must make sure that their abundances in s-process elements, inherited by their AGB progenitors, are in accordance with observations.

Another factor that can affect our results is the assumption of the instantaneous mixing approximation.<fnref n="18" /> In reality, SNe ejecta may only mix with a fraction of the available cold gas (see Salvadori et al. 2019 and Magg et al. 2020 for an estimate of the dilution factor). The resulting abundances in the SF clouds are higher, due to less dilution, and may even differ from the gross yields of the SNe (towards higher [C/Fe]; Ritter et al. 2015). However, the minimum dilution mass can not be arbitrarily small but is limited by the mass enclosed within the final size of the SN remnant (Magg et al. 2020). Models that take into account realistic prescriptions for the diffusion of SNe ejecta still find it challenging to reproduce CEMP stars with A(C)~7–7.5 (Chiaki et al. 2020; Jeon et al. 2021; Vanni et al., in preparation). Komiya et al. (2020) find that only an extremely inefficient mixing of SN yields can reproduce the highest [C/Fe] CEMP-no stars, but this results in an inconsistent MDF. They concluded that binary mass transfer from AGB stars is necessary to explain the [C/Fe] abundances of [Fe/H] < −4 stars. Sarmento et al. (2019) do manage to reach A(C) ≳ 7.5 using an Pop III IMF with $m_{\text{ch}} = 60–120$ M$_{\odot}$ and $m_* = [20–120]$ M$_{\odot}$ (the range yielding the highest [C/Fe]), but do not report on their predicted MDF and CEMP fraction. We find that adopting this IMF in our model results in a too flat MDF at [Fe/H] < −4, inconsistent with the SAGA observations.

Finally, we must note that our results could also depend on the adopted merger tree. Chen et al. (2023), for example, find that the predicted MDFs in different MW-like analogues can differ by ~1 dex at [Fe/H] = −4. We plan to explore the level of this dependence in a future work.

## 6.2 Key observables and their intrinsic uncertainties

We have shown that we can constrain the properties of primordial SNe by comparing our model predictions to observations of VMP stars in the Galactic halo. In particular, we find that both the mixing and the explosion energies as well as the IMF of Pop III stars have a strong impact on the present day CEMP-no fraction, the MDF and the CDF.

Stellar mixing affects strongly the halo MDF at [Fe/H] < −4, and the CDF at [C/Fe] > +2 (Fig. 9), where the sample of high-resolution follow-up observations can be deemed unbiased. Instead,

<footnote n="17">¹⁷The mass loading factor $\eta \equiv \dot{M}_{\text{gas,ej}}/\text{SFR} \propto 1/\dot{m}_{\text{esc}}^2$ (equation (10)) is a decreasing function of redshift (Barkana & Loeb 2001), therefore a complete blown-out of the gas at fixed $M_{\text{vir}}$ occurs more easily at low $z$.</footnote>

<footnote n="18">¹⁸Notice that under the instantaneous mixing approximation, including lighter minihaloes would move our [C/Fe]–[Fe/H] relation leftwards, i.e. further away from the observations, since metal enriched gas escapes more easily from less massive haloes, thereby lowering [Fe/H] while keeping [C/Fe] (which depends only on the assumed metal yields) ~constant.</footnote>

<footnote n="15">¹⁵For example one can easily infer from equations (5) and 12 that at fixed [C/Fe], [Fe/H] $\propto \epsilon_{\text{SF}}$.</footnote>

<footnote n="16">¹⁶$A(C) \equiv \log(N_{\text{C}}/N_{\text{H}}) + 12$</footnote>
the effect of adopting different Pop III IMFs and EDFs appears more
prominent when we consider a broader metallicity range, i.e. a MDF
extending from the lowest [Fe/H] up to [Fe/H] > −3 (Figs 3 and
4), and a CDF extending from the highest carbonicities down to
[C/Fe] < +1 (Figs 6 and A1). Unfortunately, the currently available
observations at these abundance ranges are incomplete, and follow-
up is biased towards the lowest metallicities and highest carbonicities.

In addition, all the aforementioned observables suffer from large
observational errors. Arentsen et al. (2022) report that there are
significant systematic differences in the carbon abundances among
various surveys of Galactic halo stars that can translate to more than
50 per cent differences in the estimated CEMP fractions (see their
Fig. 1). Similar uncertainties apply to the determination of [Fe/H] and
the MDF (see e.g. fig. 12 of Youakim et al. 2020). These systematics
can arise from different resolution and pipeline approaches, different
assumptions in the employed synthetic grids, and/or comparison of
stars in different evolutionary phases.

An additional source of large systematic errors comes from the
simplifying assumption of 1D, local thermodynamic equilibrium
(LTE) hydrostatic model atmospheres that is used in standard
spectroscopic abundance analyses. Accounting for 3D non-LTE
effects has been found to lower [C/Fe] estimates by as much as
~1 dex while raising [Fe/H] by < 0.15 dex (Collet, Asplund &
Trampedach 2006; Norris & Yong 2019; Amarsi et al. 2019a; Amarsi,
Nissen & Skúladóttir 2019b). Naturally, this has a dramatic effect on
the fraction of CEMP stars; Norris & Yong (2019) found that after
applying 3D non-LTE corrections, the Yong et al. (2013a) CEMP-
no fraction at -4.5 ≤ [Fe/H] ≤ -3 is reduced by ~ 60 per cent
while the number of CEMP-no stars in the Yoon et al. (2016) sample
decreases by ~ 73 per cent. Correcting for 3D non-LTE effects will
also move the observed [C/Fe]-[Fe/H] stellar distribution downwards
and, perhaps, resolve the discrepancy with our predictions (Fig. 7).

Finally, several CEMP stars with A(C) > 6.5 (and all of them at
[Fe/H] < -4.5) that are not reproduced by our model, have either no
Ba measurements or have only upper limits for barium enhancement
at [Ba/Fe] > 0.6. If a high Ba enhancement is confirmed for those
stars in the future, then their high carbonicities could be explained by
enrichment from a Pop II AGB progenitor (Rossi et al. 2023) or mass
transfer from a Pop III/II AGB companion (Komiya et al. 2020).

7 CONCLUSIONS AND FUTURE OUTLOOK

For the first time, we explore the EDF of the first SNe in the context
of a cosmological galaxy formation model of a MW-analogue. Our
model follows the formation and evolution of individual Pop III stars,
which is uniquely determined by their initial mass, stellar mixing
and explosion energy. Their contribution in the chemical enrichment
of their host minihaloes is imprinted in the present day properties
of VMP galactic halo stars, such as their MDF, CDF and CEMP
fractions. We draw the following main conclusions:

(i) Pop III EDF. The fraction of CEMP stars, $F_{\text{CEMP}}$, is highly sensitive to the primordial EDF, especially at $[\text{Fe/H}] < -3$. Assuming an EDF of the form: $dE/dN \propto E^{-\alpha_e}$, we find that we can reproduce the observed CEMP fractions for $\alpha_e \sim 1-2.5$, depending on the adopted IMF for Pop III stars (Fig. 5 and Table 1). This value corresponds to a $\sim 40-90$ per cent probability for Pop III stars with $m_* = 10-100 M_{\odot}$ to explore as faint SNe, and a 20 per cent–0.5per cent probability for them to explode as hypernovae (intermediate energy SNe have intermediate probabilities; Fig. 1). The effect of the Pop III EDF (and of their IMF) on the halo MDF is only prominent at

$[\text{Fe/H}] \lesssim -3$, but there the observational uncertainties are so large
that they render any comparison inconclusive (Figs 3 and 4).

(ii) Pop III initial mass function. A top-heavy primordial IMF (with characteristic mass $m_{\text{ch}} = 100 M_{\odot}$ in the range 0.1–1000 $M_{\odot}$) is disfavoured, as it underestimates the CEMP fraction and results in a too steep CDF, even if all Pop III stars with $m_* = 10-100 M_{\odot}$ explode as faint SNe (Figs 5 and 6).

(iii) Pop III stellar mixing. At a given EDF and IMF, lower mixing for Pop III stars results in a flatter MDF, higher CEMP fractions, and a CDF skewed towards higher [C/Fe]. We find that typically very low mixing ($f_{\text{mix}} \le 0.0631$ as provided by Heger & Woosley 2010) is required to reproduce the observations (Figs 9 and 10).

(iv) Pop II descendants. The great majority of VMP stars lie at 0 < [C/Fe] < +0.7, i.e. they are C-normal. We predict that these stars have been predominantly polluted by normal Pop II SNe, in agreement with recent studies investigating the abundance patterns of C-normal stars and their small star-to-star scatter (Vanni et al. 2023). In addition, we find a population of CEMP stars at [Fe/H] $\gtrsim -2.8$, which were born from gas enriched by Pop II AGB stars.

(v) Pop III descendants. Regardless of the assumed model, all CEMP stars at $[\text{Fe/H}] \lesssim -2.8$ have been enriched to $> 20$ per cent by Pop III progenitors. This value increases to $> 95$ per cent at $[\text{C/Fe}] \gtrsim +2$ (Fig. 7). At fixed [C/Fe], CEMP stars with the lowest metallicities are faint SNe descendants, while as we move to higher [Fe/H] the contribution of higher energy Pop III SNe prevails. According to our results, VMP stars with $[\text{C/Fe}] \lesssim 0$ are predominantly imprinted by primordial hypernovae (at $[\text{Fe/H}] \lesssim -2.5$) and PISNe (at $[\text{Fe/H}] \gtrsim -2.5$; Fig. 8).

We have demonstrated that the Pop III EDF can be equally
important to their IMF in shaping the abundances of EMP halo stars.
We find that only EDFs that are weighted towards low explosion
energies combined with bottom heavy IMFs (even if they extend to
$1000 M_{\odot}$) can reproduce simultaneously the MDF, the CDF, and
the fraction of CEMP stars in the Galactic halo. However, this
comparison alone does not allow a tighter constraint on the Pop
III IMF, mixing and EDF due to the degeneracies between them
and, most importantly, to the large uncertainties associated with the
observed relations.

We have shown that, regardless of the assumed model, the
descendants of each type of primordial SNe, always appear at specific
regions in the [C/Fe]–[Fe/H] diagram. However, their prevalence
there varies depending on the IMF and EDF of Pop III stars (Figs 8
and B2). In a following study, we will quantify this variation
and compare with the observed fractions of confirmed Pop III-
enriched stars. In addition, we plan to follow additional chemical
elements, or even full abundance patterns (e.g. Vanni et al., in
preparation). This way, we will gain further insight into the prop-
erties of primordial SNe and potentially break the aforementioned
degeneracies.

Additional constraints from forthcoming large spectroscopic sur-
veys, such as WEAVE (Dalton 2016) and 4MOST (Christlieb et al.
2019), as well as surveys dedicated to identifying Pop III descendants
(e.g. Aguado et al. 2023) and complementary studies of high-z
gaseous absorption systems imprinted by the first stellar generations
(e.g. Saccardi et al. 2023), will greatly boost our efforts to unveil the
nature of the first SNe.

ACKNOWLEDGEMENTS

The authors acknowledge support from the ERC starting grant NEFERTITI H2020/808240.
DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

Abate C., Pols O. R., Izzard R. G., Karakas A. I., 2015, A&A, 581, A22

Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93

Agarwal B., Smith B., Glover S., Natarajan P., Khochfar S., 2016, MNRAS, 459, 4209

Aguado D. S. et al., 2022, A&A, 668, A86

Aguado D. S. et al., 2023, MNRAS, 520, 866

Aguado D. S., Allende Prieto C., González Hernández J. I., Rebolo R., 2018, ApJ, 854, L34

Aguado D. S., González Hernández J. I., Allende Prieto C., Rebolo R., 2019, ApJL, 874, L21

Amarsi A. M., Nissen P. E., Asplund M., Lind K., Barklem P. S., 2019a, A&A, 622, L4

Amarsi A. M., Nissen P. E., Skúladóttir Á., 2019b, A&A, 630, A104

Aoki W., Beers T. C., Christlieb N., Norris J. E., Ryan S. G., Tsangarides S., 2007, ApJ, 655, 492

Aoki W., Norris J. E., Ryan S. G., Beers T. C., Ando H., 2002, ApJ, 567, 1166

Applebaum E., Brooks A. M., Quinn T. R., Christensen C. R., 2020, MNRAS, 492, 8

Arentsen A. et al., 2021, MNRAS, 505, 1239

Arentsen A., Placco V. M., Lee Y. S., Aguado D. S., Martin N. F., Starkenburg E., Yoon J., 2022, MNRAS, 515, 4082

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Bañados E. et al., 2018, Nature, 553, 473

Barkana R., Loeb A., 2001, Phys. Rep., 349, 125

Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379

Barnes J. E., Hernquist L. E., 1991, ApJ, 370, L65

Beers T. C., Christlieb N., 2005, ARA&A, 43, 531

Bekki K., 2014, MNRAS, 438, 444

Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

Blanchard A., Valls-Gabaud D., Mamon G. A., 1992, A&A, 264, 365

Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529

Bluck A. F. L., Maiolino R., Sánchez S. F., Ellison S. L., Thorp M. D., Piotrowska J. M., Teimoorinia H., Bundy K. A., 2020, MNRAS, 492, 96

Bond J. R., Arnett W. D., Carr B. J., 1984, ApJ, 280, 825

Bonifacio P. et al., 2015, A&A, 579, A28

Bonifacio P. et al., 2018, A&A, 612, A65

Bonifacio P. et al., 2021, A&A, 651, A79

Bournaud F., 2016, in Laurikainen E., Peletier R., Gadotti D eds, Astrophysics and Space Science Library. Vol. 418, Galactic Bulges. Springer, Cham, p. 355

Bromm V., 2013, Rep. Prog. Phys., 76, 112901

Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23

Bromm V., Ferrara A., Coppi P., Larson R., 2001, MNRAS, 328, 969

Brook C. B., Kawata D., Scannapieco E., Martel H., Gibson B. K., 2007, ApJ, 661, 10

Caffau E. et al., 2011, Nature, 477, 67

Caffau E. et al., 2016, A&A, 595, L6

Carollo D. et al., 2010, ApJ, 712, 692

Carollo D. et al., 2012, ApJ, 744, 195

Cattaneo A., Koutsouridou I., Tollet E., Devriendt J., Dubois Y., 2020, MNRAS, 497, 279

Cautun M. et al., 2020, MNRAS, 494, 4291

Cayrel R. et al., 2004, A&A, 416, 1117

Chan C., Müller B., Heger A., Pakmor R., Springel V., 2018, ApJ, 852, L19

Chen L.-H., Pillepich A., Glover S. C. O., Klessen R. S., 2023, MNRAS, 519, 483

Chen X., Hu L., Wang L., 2021, ApJ, 922, 15

Chiaki G., Wise J. H., Marassi S., Schneider R., Limongi M., Chieffi A., 2020, MNRAS, 497, 3149

Chiti A. et al., 2018, ApJ, 856, 142

Chomiuk L., Povich M. S., 2011, AJ, 142, 197

Christlieb N. et al., 2002, Nature, 419, 904

Christlieb N. et al., 2019, The Messenger, 175, 26

Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011, Science, 331, 1040

Colgate S. A., White R. H., 1966, ApJ, 143, 626

Collet R., Asplund M., Trampedach R., 2006, ApJ, 644, L121

Cooke R. J., Madau P., 2014, ApJ, 791, 116

Dalton G., 2016, ASP Conf. Ser. Vol. 507, Multi-Object Spectroscopy in the Next Decade: Big Questions, Large Surveys, and Wide Fields. Astron. Soc. Pac., San Francisco, p. 97

Danforth C. W., Shull J. M., 2008, ApJ, 679, 194

de Bennassuti M., Salvadori S., Schneider R., Valiante R., Omukai K., 2017, MNRAS, 465, 926

Dehnen W., McLaughlin D. E., Sachania J., 2006, MNRAS, 369, 1688

Dekel A., Birnboim Y., 2006, MNRAS, 368, 2

Dopcke G., Glover S. C. O., Clark P. C., Klessen R. S., 2013, ApJ, 766, 103

Ezzeddine R. et al., 2019, ApJ, 876, 97

Ferrière K. M., 2001, Rev. Mod. Phys., 73, 1031

Fox A. J., Richter P., Ashley T., Heckman T. M., Lehner N., Werk J. K., Bordoloi R., Peeples M. S., 2019, ApJ, 884, 53

François P. et al., 2018, A&A, 619, A10

Frebel A. et al., 2005, Nature, 434, 871

Frebel A., Chiti A., Ji A. P., Jacobson H. R., Placco V. M., 2015, ApJ, 810, L27

Frebel A., Norris J. E., 2015, ARA&A., 53, 631

Frebel A., Simon J. D., Kirby E. N., 2014, ApJ, 786, 74

Fryer C. L., Woosley S. E., Heger A., 2001, ApJ, 550, 372

Ganguly R., Sembach K. R., Tripp T. M., Savage B. D., 2005, ApJS, 157, 251

Gelli V., Salvadori S., Pallottini A., Ferrara A., 2020, MNRAS, 498, 4134

Gilmore G., Norris J. E., Monaco L., Yong D., Wyse R. F. G., Geisler D., 2013, ApJ, 763, 61

Graziani L., de Bennassuti M., Schneider R., Kawata D., Salvadori S., 2017, MNRAS, 469, 1101

Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, MNRAS, 449, 3137

Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ, 737, 75

Gunn J. E., Gott J. R. III., 1972, ApJ, 176, 1

Hansen T. et al., 2014, ApJ, 787, 162

Hansen T. T., Andersen J., Nordström B., Beers T. C., Placco V. M., Yoon J., Buchhave L. A., 2016a, A&A, 586, A160

Hansen T. T., Andersen J., Nordström B., Beers T. C., Placco V. M., Yoon J., Buchhave L. A., 2016b, A&A, 588, A3

Hartwig T. et al., 2018, MNRAS, 478, 1795

Hartwig T. et al., 2022, ApJ, 936, 45

Hausammann L., Revaz Y., Jablonka P., 2019, A&A, 624, A11

Heger A., Woosley S. E., 2002, ApJ, 567, 532

Heger A., Woosley S. E., 2010, ApJ, 724, 341

Hirano S., Bromm V., 2017, MNRAS, 470, 898

Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS, 448, 568

Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60

Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N., 2016, ApJ, 824, 119

Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011, Science, 334, 1250

Howes L. M. et al., 2016, MNRAS, 460, 884

Ishigaki M. N., Tominaga N., Kobayashi C., Nomoto K., 2014, ApJ, 792, L32

Ishigaki M. N., Tominaga N., Kobayashi C., Nomoto K., 2018, ApJ, 857, 46

Ishiyama T., Sudo K., Yokoi S., Hasegawa K., Tominaga N., Susa H., 2016, ApJ, 826, 9

Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439

Iwamoto N., Umeda H., Tominaga N., Nomoto K., Maeda K., 2005, Science, 309, 451
Jeon M., Bromm V., Besla G., Yoon J., Choi Y., 2021, MNRAS, 502, 1

Ji A. P., Frebel A., Ezzeddine R., Casey A. R., 2016, ApJ, 832, L3

Kapferer W., Sluka C., Schindler S., Ferrari C., Ziegler B., 2009, A&A, 499, 87

Kawata D., Gibson B. K., 2003a, MNRAS, 340, 908

Kawata D., Gibson B. K., 2003b, MNRAS, 346, 135

Keller S. C. et al., 2014, Nature, 506, 463

Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2

Klessen R., 2019, in Latif M., Schleicher D. eds., Formation of the First Black Holes. World Scientific Publishing Company, Singapore, p. 67

Kobayashi C., Umeda H., Nomoto K., Tominaga N., Ohkubo T., 2006, ApJ, 653, 1145

Komiya Y., Habe A., Suda T., Fujimoto M. Y., 2010, ApJ, 717, 542

Komiya Y., Suda T., Yamada S., Fujimoto M. Y., 2020, ApJ, 890, 66

Koutsouridou I., Cattaneo A., 2019, MNRAS, 490, 5375

Koutsouridou I., Cattaneo A., 2022, MNRAS, 516, 4194

Kroupa P., Weidner C., 2003, ApJ, 598, 1076

Krumholz M. R., 2015, preprint (arXiv:1511.03457)

Lai D. K., Lee Y. S., Bolte M., Lucatello S., Beers T. C., Johnson J. A., Sivaramani T., Rockosi C. M., 2011, ApJ, 738, 51

Larson R. B., 1998, MNRAS, 301, 569

Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013, MNRAS, 433, 1607

Lee Y. S., Beers T. C., Kim Y. K., Placco V., Yoon J., Carollo D., Masseron T., Jung J., 2017, ApJ, 836, 91

Li H., Aoki W., Zhao G., Honda S., Christlieb N., Suda T., 2015, PASJ, 67, 84

Liang J., Gjergo E., Fan X., 2023, MNRAS, 522, 863

Limongi M., Chieffi A., 2012, ApJS, 199, 38

Limongi M., Chieffi A., 2018, ApJS, 237, 13

Liu B., Sibony Y., Meynet G., Bromm V., 2021, MNRAS, 506, 5247

Machida M. N., Matsumoto T., Inutsuka S.-i., 2008, ApJ, 685, 690

Maeder A., Meynet G., Chiappini C., 2015, A&A, 576, A56

Magg M. et al., 2020, MNRAS, 498, 3703

Marassi S., Chiaki G., Schneider R., Limongi M., Omukai K., Nozawa T., Chieffi A., Yoshida N., 2014, ApJ, 794, 100

McKee C. F., Tan J. C., 2008, ApJ, 681, 771

McMillan P. J., 2017, MNRAS, 465, 76

Meynet G., Ekström S., Maeder A., 2006, A&A, 447, 623

Mortlock D. J. et al., 2011, Nature, 474, 616

Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48

Nomoto K., Kobayashi C., Tominaga N., 2013, ARA&A, 51, 457

Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nucl. Phys. A., 777, 424

Nordlander T. et al., 2019, MNRAS, 488, L109

Norris J. E. et al., 2013, ApJ, 762, 28

Norris J. E., Christlieb N., Korn A. J., Eriksson K., Bessell M. S., Beers T. C., Wisotzki L., Reimers D., 2007, ApJ, 670, 774

Norris J. E., Gilmore G., Wyse R. F. G., Yong D., Frebel A., 2010, ApJ, 722, L104

Norris J. E., Yong D., 2019, ApJ, 879, 37

O'Shea B. W., Norman M. L., 2007, ApJ, 654, 66

Omukai K., 2000, ApJ, 534, 809

Omukai K., Nishi R., 1998, ApJ, 508, 141

Omukai K., Palla F., 2001, ApJ, 561, L55

Pacucci F., Loeb A., Salvadori S., 2017, MNRASL, 471, L72

Pagnini G., Salvadori S., Rossi M., Aguado D., Koutsouridou I., Skúladóttir Á., 2023, MNRAS, 521, 5699

Peng Y.-j. et al., 2010, ApJ, 721, 193

Placco V. M. et al., 2021, ApJ, 912, L32

Placco V. M., Frebel A., Beers T. C., Stancliffe R. J., 2014, ApJ, 797, 21

Placco V. M., Frebel A., Lee Y. S., Jacobson H. R., Beers T. C., Pena J. M., Chan C., Heger A., 2015, ApJ, 809, 136

Posti L., Helmi A., 2019, A&A, 621, A56

Raiteri C., Villata M., Navarro J., 1996, A&A, 315, 105

Rees M. J., Ostriker J. P., 1977, MNRAS, 179, 541

Ritter C., Herwig F., Jones S., Pignatari M., Fryer C., Hirschi R., 2018, MNRAS, 480, 538

Ritter J. S., Sluder A., Safranek-Shrader C., Milosavljević M., Bromm V., 2015, MNRAS, 451, 1190

Rossi M., Salvadori S., Skúladóttir Á., 2021, MNRAS

Rossi M., Salvadori S., Skúladóttir Á., Vanni I., 2023, MNRASL, 522, L1

Saccardi Andrea, 2023, ApJ, 948, 35

Salvadori S., Bonifacio P., Caffau E., Korotin S., Andreevsky S., Spite M., Skúladóttir Á., 2019, MNRAS, 487, 4261

Salvadori S., Ferrara A., 2012, MNRAS, 421, L29

Salvadori S., Ferrara A., Schneider R., Scannapieco E., Kawata D., 2010, MNRASL, 401, L5

Salvadori S., Schneider R., Ferrara A., 2007, MNRAS, 381, 647

Salvadori S., Skúladóttir Á., Tolstoy E., 2015, MNRAS, 454, 1320

Salvadori S., Tolstoy E., Ferrara A., Zaroubi S., 2014, MNRAS, 437, L26

Sarmento R., Scannapieco E., Côté B., 2019, ApJ, 871, 206

Scannapieco E., Kawata D., Brook C. B., Schneider R., Ferrara A., Gibson B. K., 2006, ApJ, 653, 285

Schaerer D., 2002, A&A, 382, 28

Schauer A. T. P., Glover S. C. O., Klessen R. S., Ceverino D., 2019, MNRAS, 484, 3510

Schneider R., Ferrara A., Natarajan P., Omukai K., 2002, ApJ, 571, 30

Schörck T. et al., 2009, A&A, 507, 817

Seitenzahl I. R. et al., 2013, MNRAS, 429, 1156

Silk J., 1977, ApJ, 211, 638

Simpson C. M., Grand R. J. J., Gómez F. A., Marinacci F., Pakmor R., Springel V., Campbell D. J. R., Frenk C. S., 2018, MNRAS, 478, 548

Skúladóttir Á. et al., 2021, ApJ, 915, L30

Skúladóttir A. et al., 2023, preprint (arXiv:2305.02829)

Skúladóttir Á., Tolstoy E., Salvadori S., Hill, V., Pettini, M., Shetrone, M. D., Starkenburg, E., 2015, A&A, 574, A129

Spite M., Spite F., François P., Bonifacio P., Caffau E., Salvadori S., 2018, A&A, 617, A56

Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307

Stahler S. W., Palla F., 2004, The Formation of Stars. Wiley, New York

Starkenburg E. et al., 2017, MNRAS, 471, 2587

Starkenburg E. et al., 2018, MNRAS, 481, 3838

Starkenburg E., Shetrone M. D., McConnachie A. W., Venn K. A., 2014, MNRAS, 441, 1217

Suda T. et al., 2008, PASJ, 60, 1159

Suda T., Yamada S., Katsuta Y., Komiya Y., Ishizuka C., Aoki W., Fujimoto M. Y., 2011, MNRAS, 412, 843

Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32

Susmitha A., Koch A., Sivaramani T., 2017, A&A, 606, A112

Takahashi K., Yoshida T., Umeda H., 2018, ApJ, 857, 111

Tarumi Y., Hartwig T., Magg M., 2020, ApJ, 897, 58

Teyssier R., Chapon D., Bournaud F., 2010, ApJ, 720, L149

Thielemann F. K., Nomoto K., Yokoi K., 1986, A&A, 158, 17

Tominaga N., Iwamoto N., Nomoto K., 2014, ApJ, 785, 98

Tominaga N., Umeda H., Nomoto K., 2007, ApJ, 660, 516

Tripp T. M. et al., 2003, AJ, 125, 3122

Tumlinson J., 2010, ApJ, 708, 1398

Umeda H., Nomoto K., 2003, Nature, 422, 871

van den Hoek L. B., Groenewegen M. A. T., 1997, A&AS, 123, 305

Wang W., Han J., Cooper A. P., Cole S., Frenk C., Lowing B., 2015, MNRAS, 453, 377

Weidner C., Kroupa P., 2006, MNRAS, 365, 1333

Weidner C., Kroupa P., Pflamm-Altenburg J., Vazdekis A., 2013, MNRAS, 436, 3309

Welsh L., Cooke R., Fumagalli M., 2021, MNRAS, 500, 5214

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

White S. D., Springel V., 2000, The First Stars. Springer, Berlin, p. 327

Wollenberg K. M. J., Glover S. C. O., Clark P. C., Klessen R. S., 2020, MNRAS, 494, 1871

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Yamada S., Suda T., Komiya Y., Aoki W., Fujimoto M. Y., 2013, MNRAS, 436, 1362

Yong D. et al., 2013a, ApJ, 762, 26
Yong D. et al., 2013b, ApJ, 762, 27

Yoon J. et al., 2016, ApJ, 833, 20

Yoon J., Whitten D. D., Beers T. C., Lee Y. S., Masseron T., Placco V. M., 2020, ApJ, 894, 7

Youakim K. et al., 2020, MNRAS, 492, 4986

Zolotov A. et al., 2015, MNRAS, 450, 2327

## APPENDIX A: THE EXTENDED CDF

Fig. A1 shows the cumulative CDF of inner halo stars with $[\text{Fe/H}] \le -2$ for the models of Section 4.1, extending down to $[\text{C/Fe}] = -1$. The total number of stars changes slightly between models, from ~9.500.000 for [$m_{\text{ch}}, m_{\star}^{\text{max}}$] = [100, 1000] M$_{\odot}$, i.e. the model predicting the highest number of PISNe, to ~10.000.000 for [$m_{\text{ch}}, m_{\star}^{\text{max}}$] = [10, 100] M$_{\odot}$, i.e. the model with no PISNe. Yet, the slope of the CDF shows significant variation with the IMF of Pop III stars; for $m_{\text{ch}} = 100$ M$_{\odot}$ there are less than 70 stars with $[\text{C/Fe}] > +4$, while for $m_{\star}^{\text{max}} = 100$ M$_{\odot}$ there are more than 600 stars. In addition, one can see a clear dependence of the CDFs on the Pop III EDF that becomes more prominent at high $[\text{C/Fe}]$; the lower the $\alpha_e$ parameter, i.e. the more the high energy SNe and hypernovae, the steeper the CDF at fixed IMF. This trend is obscured when the CDFs are normalized to $[\text{C/Fe}] = +2$ (Fig. 6).

## APPENDIX B: THE CASE OF A UNIFORM EDF

This appendix presents the results of models in which all types of primordial SNe (faint, cc, high energy SNe, and hypernovae) with $m_{\star} = 10-100$ M$_{\odot}$ are assumed to be equally probable. Three different Larson-type IMFs for Pop III stars are considered: two with a mass range $m_{\star} = 0.1-1000$ M$_{\odot}$ and a characteristic mass $m_{\text{ch}} = 1$ M$_{\odot}$ and $m_{\text{ch}} = 10$ M$_{\odot}$ and one with $m_{\text{ch}} = 10$ M$_{\odot}$ and $m_{\star} = 0.1-100$ M$_{\odot}$, i.e. that does not allow the formation of PISNe. The predicted MDFs and CEMP fractions in each case, are shown in Fig. B1. Fig. B2 shows the mean metal contribution from the different types of primordial SNe in each $[\text{C/Fe}]-[\text{Fe/H}]$ bin.

The primordial EDF adopted here predicts a significantly lower CEMP-no fraction than the observations of Yong et al. (2013a) and Placco et al. (2014). The descendants of the different Pop III SN types show similar properties as the ones in the models adopted in Section 4; faint SNe descendants occupy the [Fe/H]-poorest regions of the [C/Fe]-[Fe/H] diagram, while higher energy SNe dominate the metal enrichment at progressively higher [Fe/H]. Stars that have $[\text{C/Fe}] < 0$ are predominantly enriched by primordial hypernovae (at $[\text{Fe/H}] \lesssim -2.5$) and PISNe (at $[\text{Fe/H}] \gtrsim -2.5$).

<figure><img src="image_12.png" /><figcaption>Figure A1. Same as Fig. 6 but extending down to $[\text{C/Fe}] = -1$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure B1.** MDF (top) and CEMP fraction (bottom) of inner halo stars when assuming that all Pop III stars with $m_{\star} = 10-100 M_{\odot}$ have equal probability to explode as faint SNe, ccSNe, high energy SNe, and hypernovae. The colours denote the characteristic mass, $m_{ch}$ and the maximum mass, $m_{\star}^{\text{max}}$ of the Pop III IMF considered in each model, as indicated in the legend. Data points show the observations of Bonifacio et al. (2021), Yong et al. (2013a), and Placco et al. (2014) as in Fig. 9.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure B2.** Same as Fig. 8 but for models in which Pop III stars with $m_{\star} = 10-100 M_{\odot}$ have equal probability to explode as faint SNe, ccSNe, high energy SNe, and hypernovae.</figcaption></figure>

This paper has been typeset from a TeX/LATeX file prepared by the author.