The role of Pop III stars and early black holes in the 21cm signal from Cosmic Dawn

The semi-analytical code adopted throughout this work is the Cosmic Archaeology Tool (cat, Trinca et al. 2022). A thorough description of the model can be found in the original paper, where it has been used to investigate the evolution of galaxies and their nuclear black holes at $z\geq 4$. cat describes the formation of the first stars and black hole seeds in a self-consistent way, following the feedback-regulated co-evolution of nuclear BHs and their host galaxies from cosmic dawn down to the post-reionization era (Trinca et al. 2022). As such, cat provides us with the information required to estimate the strength of the 21cm global signal generated from the first stars and BHs.

cat uses the model galform to describe the redshift evolution of dark matter halos (Cole et al. 2002; Parkinson et al. 2008), while the description of star formation, chemical evolution, BH seeding and growth is imported from the semi-analytical model gqd (Valiante et al., 2011, 2016; Sassano et al., 2021). cat simulations have been run assuming a $\Lambda$CDM model and taking the cosmological parameters according to Planck Collaboration et al. (2020): $\Omega_{\Lambda}=0.685$, $\Omega_{m}=0.315$, $\Omega_{b}=0.05$, h=0.674 and $\sigma_{8}=0.81$.

The amount of UV photons produced by the first stars and BHs is crucial to compute the reionization history and it has a direct impact on the evolution of the 21cm global signal. Following Barkana & Loeb (2001), we estimate the ionized hydrogen fraction as:

where $f_{\rm esc}$ is the escape fraction of ionizing photons and $\dot{n}_{\rm ion}(z)$ is the ionizing photon rate density (number of photons cm^-3 s^-1). We compute the latter quantity from the luminosities of the stellar populations and accreting BHs predicted by cat. The term $F(z_{\rm obs},z)$ accounts for the recombination rate of ionized hydrogen. It depends on the clumping factor $C$, which quantifies the patchiness of the IGM, on the present-day hydrogen number density $\bar{n}^{0}_{\rm H}$ and on the case-B recombination factor $\alpha_{\rm B}=2.6\times 10^{-13}$cm³  s^-1 (valid for hydrogen at $T=10^{4}$ K, Barkana & Loeb 2001). The clumping factor and $f_{\rm esc}$ are free parameters, that we tune in order to obtain ionization histories consistent with the latest Planck Collaboration et al. (2020) measurements and with the ionizing photon emissivities constrained from observations at $z<6$ (Bolton & Haehnelt, 2007; Kuhlen & Faucher-Giguère, 2012; Becker & Bolton, 2013; Becker et al., 2021).

A thorough description of the reionization scenarios predicted by cat will be presented in Trinca et al. (in prep). In the following, we adopt a redshift-dependent escape fraction $f_{\rm esc}(z)=0.02\times[(1+z)/5)]^{3/2}$, in agreement with what has been done in other recent studies (e.g Mutch et al., 2016).

At each redshift, cat is able to compute for each galaxy its SFR ($M_{\odot}$/yr), stellar metallicity (this is crucial to distinguis between Pop III and Pop II/I stars), nuclear BH mass and accretion rate. The stellar contribution to the ionizing (E > 13.6 eV) photon rate can be recovered using time- and metallicity-dependent UV luminosities from Bruzual & Charlot (2003) for Pop II/I stars and from Schaerer (2002) for Pop III stars (see also Fig. 3 in de Bennassuti et al. 2017). Following Valiante et al. (2016) and Trinca et al. (2022), the ionizing UV photon rate produced by each accreting BH depends on the Eddington ratio and on the BH mass. We then account for obscuration in the UV band, assuming an obscured fraction of (Merloni et al., 2013):

where $L_{\rm X,2-10\,keV}$ is the logarithmic X-ray AGN luminosity in the 2 - 10 keV band and the best-fit parameters are: $A=0.56$, $l_{0}=43.89$ and $\sigma_{\rm X}=0.46$. The procedure we adopt to compute $L_{\rm X,2-10\,keV}$ is described in Section 2.4. In our model, we assume that only a fraction equal to $1-f_{\rm obs}$ of all AGNs is unobscured and contributes to the ionizing photon emissivity.

The resulting $\dot{n}_{\rm ion}$ multiplied by the above $f_{\rm esc}$ is shown in the left panel of Fig. 1, with the separate contributions coming from Pop III and Pop II/I stars (these are the two dominant contributions). We compare our models with some observational constraints at $z<6$ as indicated in the Figure. Overall, the dominant contribution is provided by Pop II/I emission (thin black line), with Pop III stars (dashed red line) dominating during the first phase of the evolution, at $z+1>20$ and disappearing at $z+1<10$ as a consequence of the increased level of metal enrichment in the galaxies and the IGM. The Pop III contribution to the ionizing photon rate has a less smooth evolution compared to the Pop II/I especially in the redshift interval $z+1\sim 10-16$. The right panel of Fig. 1 shows the same quantities but eliminating the contributions of galaxies hosted by mini-halos (systems with $T_{\rm vir}<10^{4}$ K). Both Pop III and Pop II contribution at $z+1<15$ are not affected with the Pop II contribution largely dominating over the Pop III one. However, at high-$z$ both contributions are weaker. The impact on the ionizing emission from Pop II stars is only mildly reduced, while the Pop III one suffers of a much stronger suppression ($\sim$ 1 order of magnitude at $z+1>22$). In addition, this latter contribution is more "bursty". This evolution reflects the fact that at high-$z$ a large fraction of Pop III (and a smaller fraction of Pop II) star formation occurs in mini-halos. Conversely, at lower redshift, most of the star formation occurs in metal enriched atomic cooling halos. At $z+1<15$ the emission of ionizing photons is halted as a consequence of the effect of the LW background that suppress molecular cooling in these less massive systems.

Overall, our model predicts a production of ionizing photons at $z+1=5-7$ consistent with observations reported by Bolton & Haehnelt (2007), Kuhlen & Faucher-Giguère (2012) Becker & Bolton (2013) and Becker et al. (2021). In addition, depending on the adopted clumping factor²²2In Trinca et al. (in prep) we explore different reionization histories depending on the adopted evolution of the clumping factor $C(z)$. The values quoted above refer to a model obtained using the maximum value between a constant clumping factor of 3 and the analytic form presented in Iliev et al. (2007) ($z_{\rm reio}=4.67$, $\tau_{e}=0.050$), and a model obtained using a constant clumping factor of 3 ($z_{\rm reio}=6$, $\tau_{e}=0.065$)., we obtain a total Thomson scattering optical depth of $\tau_{e}=0.050-0.065$, consistent with the latest measurements of Planck Collaboration et al. (2020), and a complete reionization redshift of $z_{\rm reio}=4.67-6$. All the models predict the same high-redshift reionization history with an optical depth in the redshift range $15\leq z\leq 30$ of $\tau_{e}(15-30)=0.050-0.065$. We refer to Trinca et al. (in prep) for a more in-depth discussion of the reionization histories predicted by cat.

The Lyman-$\alpha$ background impacts the evolution of the 21cm global signal as it determines the coupling between the kinetic and spin temperature through the coupling coefficient $x_{\alpha}$, which in turn depends on the Lyman-$\alpha$ flux $J_{\alpha}$ and on its evolution in time.

Following Dayal et al. (2008), we simply assume the Lyman-$\alpha$ photon rate to be set by the amount of UV ionizing photons produced by stars and accreting BHs which do not escape from the galaxies where they are generated, hence $\dot{n}_{\alpha}=\dot{n}_{\rm{ion}}\,(1-f_{\rm{esc}})$. Thus, the Lyman-$\alpha$ flux is computed as:

with $\Delta\nu=\nu_{\rm LL}-\nu_{\alpha}$ the frequency interval between the Lyman limit and the Lyman-$\alpha$ line and $f_{\alpha}$ the fraction of Lyman-$\alpha$ photons that escape the galaxy without being destroyed by dust, that we assume to be equal to 1³³3Here we are mostly interested in the 21cm global signal at cosmic dawn, hence we do not expect significant dust enrichment in the interstellar medium of galaxies at $z>15$.. The main advantage of this formalism is that we do not need any new quantity to compute $J_{\alpha}$ since we use the ionizing UV photon rate already discussed in the previous section.

Fig. 2 shows the results including all stars and accreting BHs predicted by cat (solid line), and without the contribution of mini-halos (dashed line). The evolution of the Lyman-$\alpha$ flux shows the same trend of the UV ionizing emission (see Fig. 1). At high redshift, sources hosted in mini-halos strongly contribute to the Lyman-$\alpha$ background while for $z+1\leq 15$ the evolution of ${\rm J_{\alpha}}$ is determined by sources in atomic cooling halos. Thus, the contribution of mini-halos is relevant in determining a coupling between $T_{\rm S}$ and $T_{\rm K}$ at early times, while at later times it becomes negligible.

The X-ray background is responsible for gas heating and determines the evolution of the kinetic temperature of the gas, as shown by Eq. (4). In particular, we need to estimate the X-ray emissivity $\epsilon_{\rm X}$ (erg s^-1 Mpc^-3) produced by X-ray binaries and accreting BHs.

Following Grimm et al. (2003), the first contribution is modeled according to the locally observed correlation between the star formation rate and the bolometric X-ray luminosity, $L_{\rm X}$ (erg s^-1):

This relation is based on CHANDRA and ASCA observations of nearby star-forming galaxies. In particular, the value of the constant correctly matches the proportionality between the X-ray luminosity of a galaxy arising from high-mass X-ray binaries (HMXBs) and the SFR for high values of SFR. The contribution of low-mass X-ray binaries (LMXBs) is not considered as it has been found to be subdominant with respect to the one coming from HMXBs at $z\geq 2.5$ (Fragos et al. 2013; Madau & Fragos 2017). Our ignorance of the accuracy of this relation at high-redshift (in particular the proportionality constant) is factorized in the parameter $f_{\rm X}$⁴⁴4This parameter likely will evolve with redshift and/or metallicity but for simplicity we assume it to be a constant. The dependence on metallicity though has not a very strong impact as found in Kaur et al. (2022).. We can thus express $\epsilon_{\rm X}$ as (Chatterjee et al., 2019):

where $f_{\rm X,h}=0.2$ is the fraction of the X-rays that contribute to heating (Shull & van Steenberg, 1985).

In this work we assume $f_{X}=1$. This choice is consistent with the recent results of the HERA collaboration Abdurashidova et al. (2022), which imply that galaxies at higher redshift were more efficient at producing X-rays than local ones, constraining the ratio between $L_{\rm X}$ and SFR between 1.6 $\times 10^{40}$ and 8 $\times 10^{41}$. However, it is worth noting that HERA observations are taken at z = 8 and z = 10 and do not provide constraints at earlier times.

Finally, we need to account for the fact that X-ray photons will not be equally efficient at heating the IGM. The higher the energy of the photon, the higher its mean free path $\lambda_{\rm X}$, and lower the probability of interacting with the IGM. As shown by Fialkov et al. (2014) and more recently by Madau & Fragos (2017) a photon will not interact if $\lambda_{\rm X}>c/H(z)$, yielding an upper limit on the energy of X-ray photons that can be absorbed by the IGM of:

which gives $E_{\rm lim}\in$ [1.23, 3.29] keV in the redshift interval of our simulation. Following Muñoz et al. (2022), we adopt a power-law spectral energy distribution (SED) for HMXBs with a spectral index $\alpha_{\rm X}=-1$ and a low-energy cut-off $\rm{E}_{0}=0.5$. This SED is normalized to the bolometric X-ray emissivity computed from Eq. (21). This choice is broadly consistent with emission spectra of HMXBs in low-metallicity environments (e.g Fragos et al., 2013; Madau & Fragos, 2017; Das et al., 2017; Qin et al., 2020). We highlight that we are considering a simple SED that we apply to Pop III and Pop II/I stars. This is a crude approximation, and we defer to a future work a more in-depth analysis of the X-ray SED. The X-ray emissivity that enters in Eq. (4) is then $\epsilon_{\rm X,E<E_{lim}}=\int_{E_{0}}^{E_{\rm lim}}\epsilon_{\rm X}(E)dE$.

The X-ray contribution from BHs is calculated by applying a correction factor $K_{\rm X}$ to the bolometric luminosity of the AGN so that $L_{\rm X,2-10\,keV}=L_{\rm bol}/K_{\rm X}$. We take the correction factor from Duras et al. (2020):

where $a=10.96,b=11.93$, and $c=17.79$. This relation is valid for both type 1 and type 2 AGNs (the sampled analysed included $\sim$1000 type 1 and type 2 AGNs). The hard X-ray spectrum of AGNs is important to determine the AGN obscuration fraction in the UV band (see Eq. (18)). As stressed above, the contribution which is relevant to heat up the neutral IGM comes from the soft X-ray band (0.5-2 keV). In this case we applied a correction factor $K_{\rm X}$ from Shen et al. (2020) (which has been found also by Graziani et al. 2018):

with $c_{1}=5.712$ , $k_{1}=-0.026$ , $c_{2}=17.67$ and $k_{2}=0.278$. As we have already done for UV photons, we must account for the fact that some X-ray photons (probably a smaller fraction than the UV) will not escape from the galaxy where they are produced. Following Trinca et al. (2022) our reference model for X-ray obscuration is the one proposed by Ueda et al. (2014). This correction is expressed in terms of the $\psi(L_{\rm X},z)$ parameter which represents the fraction of absorbed AGNs. It is expressed as a linear function of $\log L_{\rm X}$ within a range $\psi_{\rm min}=0.2$ and $\psi_{\rm max}=0.84$:

where $\beta=0.24$ and $\psi_{43.75}(z)$ represents the absorption fraction of AGNs with $\log L_{\rm X}=43.75$ located at $z$. This redshift dependence disappears for $z\geq 2$ (which is our case) and $\psi_{43.75}\simeq 0.73$. The X-ray luminosity from AGNs is then $L_{\rm X,unobs}=L_{\rm X}\,[1-\psi(L_{\rm X})]$. It is worth mentioning here that, all the fitting equations for the AGNs spectra adopted in this section, have been calibrated to observations at low redshifts. Extrapolating those up to the cosmic dawn it may not be correct and lead to an overestimation of the contribution of AGNs to the radiative background. This will be better discussed in section 4.2.

Fig. 3 shows the total soft X-ray emissivity (the one contributing to IGM heating), disentangling between the various components: Pop II/I stars (solid thin line), Pop III stars (red dashed line) and unobscured AGNs (light blue thin line). The left panel shows the results of the entire sample of cat galaxies and the right panel illustrates how the results are affected by removing the contribution of mini-halos. Differently from the ionizing background, the soft X-ray emissivity from AGNs is dominant already at very high redshift ($z+1\sim 34$) in both panels. Indeed, a significant contribution of the bolometric X-ray emission from HMXBs falls in the hard regime (E $\geq 2$keV), while AGNs have softer X-ray spectra. Removing the contribution of mini-halos mildly affects the global evolution of the soft X-ray emission. As for the ionizing emission, the two panels in Fig. (3) at $z+1>15$ differ mostly for the red dashed curve (Pop III contribution), which is largely reduced when we are not accounting for sources hosted in mini-halos. Another interesting feature of the Pop III HMXBs emission is its "burstiness". Since we compute the X-ray luminosity to be proportional to the SFR, the scattered evolution of the dashed line in Fig. (3) reflects the "burstiness" of Pop III star formation (Furlanetto & Mirocha, 2022). The properties of Pop III star formation predicted by cat will be discussed in a future work (Trinca et al. in prep.). Here we just highlight that the latest episodes of Pop III star formation occur at $z+1\sim 10$. The Pop II contribution (solid thin line) instead shows a smoother evolution and it becomes dominant over the AGN’s contribution at $z+1<10$.

In Eq. (1) we have implicitly assumed that the only 21cm background against which we expect to observe the signal is the CMB. If a strong radio background with a brightness temperature $T_{\rm rad}>T_{\gamma}$ is already present at early times, the background brightness temperature would increase by a multiplicative factor of $1+T_{\rm rad}/T_{\gamma}$ (Feng & Holder, 2018). This means that in Eq. (1) we would need to substitute $T_{\rm rad}+T_{\gamma}$ in place of $T_{\gamma}$, which leads to a decrease in the differential brightness temperature (DBT). The background radio temperature determines the 21cm global signal also through Equations 2 and 3 changing the spin temperature evolution and the Lyman-$\alpha$ coupling coefficient.

As pointed out by Mittal & Kulkarni (2022), it is extremely challenging to explain the peculiar shape of the EDGES signal with just Lyman-$\alpha$ and X-ray emission from stars. The presence of a radio background might be considered in order to explain the large amplitude of the absorption signal detected by EDGES (Bowman et al. 2018a). Having a deep absorption signal means that $T_{\rm s}$ is much lower than $T_{\rm rad}$, therefore, an extra radio background would make the HI spin temperature much lower than $T_{\rm rad}$ before the X-ray heating starts to be effective.

The EDGES signal could also be explained if there is a (weak and non-gravitational) interaction between baryons and dark matter (Barkana 2018). In virtue of their earlier decoupling, dark matter particles are colder than the early cosmic gas and, if such an interaction between baryons and dark matter exists, it would provide a cooling mechanism decreasing $T_{\rm K}$ (and $T_{\rm S}$) and therefore increasing the difference between $T_{\rm S}$ and $T_{\rm rad}$. This second possibility will not be considered in this work, and in order to reproduce the EDGES signal we will estimate the radio background produced by the population of sources predicted by cat.

Fig. 4 shows the specific intensity of the radio background as a function of redshift assuming $f_{\rm R}=100$ (the motivation for this choice will be discussed in Section 4.3). The solid line illustrates the result when all the BHs predicted by cat are assumed to contribute, while the dashed and the dash-dotted lines show the radio background intensity considering only BHs that are accreting with an Eddington ratio $f_{\rm edd}\geq 0.01$ and 0.1, respectively. The motivation to disentangle between BHs accreting at different paces (consistently with Ewall-Wice et al. 2019) is that the radio emissivity is expected to depend strongly on the BH accretion rate.

The radio background predicted by our model starts to be built up as soon as the first BH seeds are formed, and rapidly increases up to $z+1\sim 16$. At lower redshift the growth is smoother. There is a very little difference between the solid and the dashed curve meaning that, in our model, a large fraction of BHs accretes gas with an Eddington ratio larger than 0.01. If we assume a more extreme cut at $f_{\rm edd}\geq 0.1$, the radio background drops by almost one order of magnitude at $z+1\leq 16$.

The thermal history predicted by the model when both stars and accreting BHs are considered is shown in Fig. 5. Here we show the redshift evolution of the spin temperature (thick dashed line), of the CMB temperature (dotted line) and of the gas kinetic temperature (thick solid line). We also show the contribution of unobscured AGN to the kinetic temperature of the gas (blue solid line) as it is the dominant one. The two panels of Fig. 5 show the same quantities computed considering all the galaxies predicted by cat (left panel) and without the sources formed in mini-halos (right panel).

At the beginning of the evolution, the IGM cools adiabatically. However, after the first sources start to form, the kinetic temperature starts to increase ($z+1\sim 27$, left panel) due to X-ray heating until it becomes hotter than the CMB temperature ($z+1\sim 20$, left panel). Thereafter $T_{\rm K}$ continues to increase reaching the maximum allowed value of $10^{4}$ K when hydrogen is fully ionized. We highlight that the evolution of $T_{\rm K}$ is consistent with the constraint obtained by Abdurashidova et al. (2022) at $z+1=8.9$ (95% confidence interval) represented by the vertical line. If we neglect the contribution of mini-halos, the evolution of $T_{\rm K}$ at low redshift is almost unaffected, while some differences can be seen at $z+1\geq 15$. The heating starts to be effective at $z+1\sim 25$ and $T_{\rm K}>T_{\gamma}$ at $z+1\sim 18$. The thin blue line in both panels is close to the thick black one (at $z+1\geq 14$ the two lines overlap), showing that the heating is largely dominated by accreting BHs at all the redshifts with a minor contribution from stars.

The spin temperature is initially coupled to the CMB temperature. Once the first structures are formed, $T_{\rm S}$ is driven toward the kinetic temperature by the Lyman-$\alpha$ photons until, at $z+1\sim 22$ (left panel), a tight coupling between $T_{\rm K}$ and $T_{\rm S}$ is reached and $T_{\rm S}\simeq T_{\rm K}$ thereafter. Again, the effect of not considering sources formed in molecular cooling halos is to delay the thermal evolution of the IGM; the tight coupling between $T_{\rm S}$ and $T_{\rm K}$ is reached at $z+1\sim 19$ when the IGM is already hotter than the CMB. While the contribution of unobscured AGN was dominant for the evolution of $T_{\rm K}$, it is negligible for the spin temperature. This is a result of the fact that AGN are highly obscured in the Lyman-$\alpha$ band. Thus they are not able to provide Lyman-$\alpha$ photons in order to couple $T_{\rm S}$ to $T_{\rm K}$. Ultimately, this tells us that this coupling must be provided by stars.

The corresponding 21cm global signal is shown by the thick solid line in the upper left panel of Fig. 6. We find that between $z+1\sim 20.5$ and $z+1\sim 32$ (44 $\leq\nu\leq$ 69 MHz) the signal is observed in absorption, with a maximum depth of $\delta T_{\rm b}\simeq-95$ mK at $z+1\sim 26.5$, which corresponds to an observed frequency of $\nu\sim 54$ MHz. This feature reflects the early coupling between the kinetic and spin temperature (see lower left panel). The early heating of the IGM due to accreting BHs (as shown in the previous section) plays an important role in determining the shape and the depth of the absorption feature. Indeed, if we neglect the contribution of AGN and leave only the stellar one (red dashed line) we predict that (i) the absorption feature becomes wider ($\delta T_{\rm b}\leq-50$ mK until $z+1\sim 17$), (ii) deeper ($\delta T_{\rm b}\sim-110$ mK at $20\leq z+1\leq 26$) and (iii) the transition from an absorption to an emission signal (which corresponds to the spin temperature to become larger than the CMB one) occurs later ($z+1\sim 14.5$, 97 MHz). If we subtract also Pop III stars and we let only Pop II stars contribute (solid thin line) we obtain a 21cm global signal history consistent with the "classical" evolution found by many reference studies in the literature (e.g. Furlanetto, 2006; Visbal et al., 2012; Fialkov et al., 2013b; Cohen et al., 2017; Muñoz et al., 2022). In our computation, the 21cm global signal driven by only Pop II stars reaches a maximum amplitude in absorption of $\delta T_{\rm b}\sim-65$ mK at $z+1\sim 18$ (78 MHz). The difference between the red dashed line and the solid thin one highlights the importance of the contribution of Pop III stars to the absorption feature as it introduces an earlier coupling between $T_{\rm S}$ and $T_{\rm K}$ (see lower left panel) which anticipates the timing and increase in the depth of the absorption feature compared to the case where only Pop II stars are present. Another effect of Pop III stars is to make the signal more "noisy". If only Pop II stars were present, the evolution of the signal would have been smoother, reflecting the "burstiness" of Pop III star formation.

Thus, as shown by the left panels of Fig. 6, the 21cm global signal encodes the properties of the sources that are present in the Universe at high redshift and their evolution from cosmic dawn to the epoch of reionization, as the various classes of sources contribute in different ways to the kinetic and spin temperature evolution of the IGM. Pop III stars are very effective at driving $T_{\rm S}$ toward $T_{\rm K}$ at high redshift, but they provide a negligible contribution to gas heating. Pop II stars are very important for the coupling and have a small impact on the heating. Finally, accreting BHs are crucial to heat the neutral hydrogen erasing most of the contribution to the signal coming from Pop II stars and leaving only the strong absorption feature caused mainly by Pop III stars.

In this work we are mostly interested in the 21cm absorption feature originating at cosmic dawn, since this provides direct evidence for the onset of Pop III star formation and may encode important information on the interplay between physical processes that regulate early star formation at $z>17$. In Fig. 7 we show a zoom-in of the 21cm global signal between $z+1=32$ and $z+1=14$, highlighting separately the contributions of Pop III (red dashed line) and Pop II (thin black line) stars (we note that, differently from Fig. 6 where the red dashed line accounts for the entire stellar contribution Pop III + Pop II, here the red dashed line refers to the Pop III contribution alone). When all sources are considered (left panel), the absorption feature reaches a maximum depth $\delta T_{\rm b}\sim-95$ mK at $z+1\sim 26.5$ (54 MHz), with a smooth decrease starting at ($z+1=34$). As we can see from the thin black and red dashed lines, the absorption signal at $z+1\geq 26$ is entirely determined by the early coupling caused by Pop III stars. At such high-$z$ the heating is still not effective and the IGM is not chemical enriched yet so that Pop II star formation episodes are rarer than Pop III ones. As a consequence, the coupling due to Lyman-$\alpha$ photon emission from Pop III stars is predicted to be the main contribution to the shape of the 21cm global signal at these redshifts. At lower redshifts, even if we do not consider the heating from unobscured AGNs, the stellar contribution to the heating of the IGM (coming from HMXBs as described in Section 2.4) is sufficient to heat up the neutral hydrogen above the CMB temperature, suppressing the 21cm global signal (see red dashed line in Fig. 6). The zoom-in of the signal (Fig. 7) shows that Pop II stars are much more effective in heating up the IGM compared to Pop III stars. The red dashed line in Fig. 7 is always negative ($\delta T_{b}\leq-70$ mK) meaning that, if we consider only Pop III stars, we are able to drive the spin temperature toward the kinetic temperature, but the latter will always be much smaller than the background radiation and thus $T_{\rm S}<T_{\gamma}$ making the signal always be in absorption. This does not happen when we compute the evolution of the 21cm global signal only from Pop II stars. The absorption feature is delayed compared to the red dashed line, but it turns into an emission signal at $z+1\sim 14.5$ ($\nu\sim$ 98 MHz). Ultimately this implies that Pop III stars are mostly important for the early coupling between $T_{\rm S}$ and $T_{\rm K}$ but not for the heating of the IGM, while Pop II stars have an impact on both.

The previous results have been obtained considering all sources that form inside both molecular (or mini-halos) and atomic cooling halos. Here we investigate the results without the contribution of mini-halos, in order to see how the 21cm global signal changes. In practice, we compute the 21cm global signal considering only sources hosted by dark matter halos with $T_{\rm{vir}}\geq 10^{4}$  K. It is important to stress, however, that this procedure is applied by post-processing cat simulations output. Hence the properties of the sources hosted in atomic cooling halos are still sensitive to the physical processes occurring in the less massive mini-halos.

The 21cm signal evolution is shown in the upper right panel of Fig. 6 (thick line). The overall profile is now quite different compared to the one obtained from the entire halo population and shown in the upper left panel of the same plot. Without the contribution of mini-halos, the absorption feature is strongly suppressed, showing a much smaller amplitude compared to the previous case ($\delta T_{\rm b}\simeq-40$ mK for $\nu\sim 55-63$). At higher frequencies, $\delta T_{\rm b}$ starts to increase until the signal vanishes at $z+1\sim 18.5$. This indicates that sources hosted in molecular cooling halos are crucial in achieving an early tight coupling between the spin and kinetic temperature. However, their contribution is soon cancelled out due to the effect of the Lyman-Werner feedback that stops the molecular cooling inside mini-halos preventing star formation. However, the Pop III star contribution, even if smaller, still dominates the shape of the absorption feature. This is visible when we compare the 21cm signal after removing AGN (red dashed line), with only Pop II stars (black thin line) and with both stars and AGN (black thick line). The red dashed line and the black thick one have a very similar evolution for $z+1>25$ indicating that the role of accreting BHs is quite visible even when we do not consider sources in mini-halos. The signal computed only from Pop II stars instead is almost unaffected, reflecting the fact that Pop II stars mostly form in more massive atomic cooling halos. As in the previous section, the zoom-in of the absorption feature in the 21cm global signal obtained after removing mini-halos (Fig. 7, right panel) shows that Pop III stars efficiently drive the coupling between $T_{\rm S}$ and $T_{\rm K}$ but do not heat up the IGM. Overall, neglecting mini-halos leads to important changes in the signal at high-$z$ when the absorption feature is predicted while at lower redshift ($z+1<15$ in our simulation) star formation in mini-halos is suppressed by radiative feedback from the Lyman-Werner background, and no longer contributes to the signal.

The results shown so far consider the redshifted CMB photons as the only radio background against which the DBT is measured ($T_{\rm R}=T_{\gamma}$). In this subsection we show the impact of an additional radio background to the 21cm signal computed assuming the whole halo population as explained in Section 2.5. We stress that, using cat outputs, we are able to compute this additional background consistently with the other radiation backgrounds that we have considered so far. The only free parameter that enters here is $f_{\rm R}$ and this can be tuned to match the amplitude of the EDGES claimed detection (Bowman et al. 2018a).

The resulting 21cm signal is shown in Fig. 8 adopting a value of $f_{\rm R}=100$ and considering BHs that are accreting at different rates, as discussed in Section 2.5 (solid, dashed and dash-dotted lines refer to all BHs, BHs with $f_{\rm edd}\geq 0.01$ and 0.1, respectively). The solid thick line shows that the signal is strongly altered by the presence of a radio background. The absorption feature between $20\leq z+1\leq 24$ has now an increased amplitude ($\delta T_{\rm b}\sim-600$ mK at $\nu\sim 59-71$ MHz). In addition, the 21cm signal is in absorption at all frequencies, with an amplitude $\delta T_{\rm b}\leq-50$ mK at $z+1\geq 11$. The evolution of the signal is now smoother across all the redshift steps, since the leading contribution is now the radio background from accreting BHs shown in Fig. 4.

In order to match the amplitude of EDGES detection Bowman et al. (2018a), a large value of $f_{\rm R}=100$ is required. If we consider only BHs that are accreting with $f_{\rm edd}\geq 0.01$ (dashed line), the absorption feature at $\nu\sim 63-74$ MHz is mildly shallower ($\delta T_{\rm b}\sim-500$ mK) while the remaining evolution is almost identical. Restricting the contribution to systems with $f_{\rm edd}\geq 0.1$ (dash-dotted line) leads to an absorption signal which is much weaker and potentially turning into an emission signal at $z+1\sim 14$ ($\nu\sim 101$).

From now on, we will take the dashed curve as the reference 21cm signal evolution when an additional radio background is considered, as this appears to be the model that best fits the amplitude (although not the timing) of the EDGES detection (see the red bar in Fig. 8). It is evident that, even considering this additional radio background we are not able to reproduce the EDGES results with our reference model. This point will be further discussed in Section 4.2 and 4.3. This choice of $f_{\rm R}$ leads to a radio background temperature consistent with the current observational constraint from Abdurashidova et al. (2022): $\log_{10}(T_{\rm{rad}}/T_{\rm K})<1.7$ at $z+1=8.9$ (95% confidence).

The main novelty of the 21cm signal model presented in this work is the inclusion of Pop III stars and BH evolution, starting from their seeding. The vast majority of current theoretical models for the 21cm global signal, identify Pop III stars with stars that are formed inside mini-halos (e.g. Qin et al., 2020; Muñoz et al., 2022) without allowing Pop III star formation in more massive atomic cooling halos. Among the few studies on the impact of Pop III stars to the 21cm global signal we mention Magg et al. (2022) and Gessey-Jones et al. (2022) (who focus respectively on the Pop III/II chemical transition and on the IMF of Pop III stars.) Other previous works considered the impact of Pop III stars on the radio and X-ray background in order to match the EDGES detection reported by Bowman et al. (2018a) (e.g. Schauer et al., 2019; Chatterjee et al., 2020; Mebane et al., 2020).

Given the large uncertainties on the parameters that determine Pop III star formation, in this work we choose not to explore such a huge parameter space, but pick a single model of Pop III star formation (see Section 2.1 and Trinca et al. in prep.) and examine its impact on the 21cm signal. The choice of a Pop III star formation model also determines the formation and the evolution of the first BHs: Pop III remnants are one of the two channels through which we form BH seeds in cat, and therefore, depending on the characteristic IMF of Pop III stars we will end up with a different number of BH seeds. Thus, for the first time, we are able to provide a self-consistent way to estimate both the impact of Pop III stars and BHs to the 21cm signal.

In Section 3 we already assessed the role of Pop III star formation in the early coupling between $T_{\rm S}$ and $T_{\rm K}$ and in the narrow 21cm absorption feature (solid thick line in Fig. 6). A detailed description of Pop III star formation in cat will be presented in a forthcoming paper (Trinca et al. in prep). Here we discuss those factors that have a direct impact on the 21cm global signal.

Firstly, we notice that the absorption feature presented in this work differs from many other previous works. This is mostly caused by the fact that cat is able to resolve all the minihalos down to $M_{\rm halo}\sim 5\times 10^{5}M_{\odot}$ and the SF episodes within those. This leads to a higher SFRD at $z>25$ compared to other works (e.g. we obtain a Pop III SFRD of 0.5 dex higher than Visbal et al. 2018 in the redshift interval 25 - 35). This small difference can also be attributed to the different MTs used. As described in Section 2.1, our semi-analytical model uses galform which is based on the Extended Press-Schechter theory and it is tuned in order to obtain merger histories consistent with the N-body Millennium simulation (Springel et al. 2005). This choice leads to a higher halo number density at $z>25$ compared to other N-body simulations ultimately impacting on the number density of star formation episodes and on the earlier coupling between $T_{\rm S}$ and $T_{\rm K}$. As the cosmic dawn is a poorly constrained epoch of the evolution of our Universe, we believe that our anticipated evolution of the 21cm global signal is not a major concern as the model that we present is self-consistent and well-anchored both to the observations at lower redshifts and, as shown in this section, to the results of the hydrodynamical simulations for Pop III star formation. Lastly, a comparison between different models of the global signal must account for the differences on how mini-halos and their relative contribution are modeled.

When Pop III star formation occurs in mini-halos, we enhance the star formation efficiency to $\rm\epsilon_{SF,PopIII}=0.15$ in order to achieve a total stellar mass formed in each star formation episode of $M_{\rm PopIII}\geq M_{\rm min,PopIII}=150M_{\odot}$, in agreement with state-of-the-art hydro-dynamical simulations (Chon et al., 2021). We then randomly sample the assumed top-heavy Pop III IMF until we saturate $M_{\rm PopIII}$ (see Valiante et al. 2016, de Bennassuti et al. 2017 and Trinca et al. 2022 for the random sampling in cat). Hence, the value of $M_{\rm PopIII}$ (and thus the adopted Pop III star formation efficiency) has an impact on the masses of newly formed Pop III stars, with larger values of $M_{\rm PopIII}$ favouring a larger frequency of Pop III stars at the high-mass end of the distribution.

Their high masses and large number are the main reason why we obtain a very early coupling between $T_{\rm S}$ and $T_{\rm K}$ and thus the absorption feature at $z=21-32$ as shown in the top-left panel of Fig. 7. Moreover, in cat a significant fraction of the total number of Pop III stars are formed inside mini-halos (this is a realistic scenario as mini-halos are the first to virialize at high redshift and are more abundant than the more massive halos). For this reason, when we consider only the sources that form inside atomic cooling halos, the impact of Pop III stars on the 21cm signal is strongly reduced for $z+1>21$. Since in our model this is the main contribution to the absorption feature at high-$z$, removing mini-halos strongly suppresses the absorption feature so that this signal would not be detectable even by the new generation of radio telescopes ($\delta T_{\rm b}>-45$ mK over the entire evolution).

The modeling of the typical masses of Pop III stars also affects the evolution of high-$z$ BHs. Pop III stellar remnants (if the initial masses of Pop III stars are large enough) provide the light BH seed population which then grow via both gas accretion and mergers. Gas accretion onto these BHs provides an important source of heating of the Universe already at high redshift, causing a quick suppression of the absorption signal that around $z+1\sim 20.5$ and then an emission signal in our reference model. As already stressed, in our model it is the large number, rather than the masses, of BH seeds that makes the accreting BHs the main contribution to the total X-ray background.

In conclusion, a reliable model of the 21cm global signal, and in particular of the high-redshift 21cm absorption feature, requires accurate modelling of Pop III star formation and their associated BH remnants. If Pop III stars are numerous, the 21cm signal will be characterised by a narrower and earlier absorption feature; the first effect is due to the increased X-ray emissivity of unobscured AGNs while the second is determined by the larger Lyman-$\alpha$ emissivity of Pop III stars.

One surprising result discussed in this paper, is the role of accreting BHs in heating the IGM already at $z+1\sim 27$, leading to a quick suppression of the absorption signal. This is a consequence of (i) the large number of light seeds formed after the death of Pop III stars and (ii) the SED of accreting BHs adopted in this work (see section 2.4). The second point comes with several caveats. The fitting functions adopted to compute the soft X-ray emission from BHs as a function of the bolometric luminosity are calibrated to observations at much lower redshift. In addition, the correction for the obscuration is calibrated on AGNs with much larger luminosities compared to the ones we have at high-$z$ (Ueda et al. 2014, Shen et al. 2020). Finally, we are implicitly assuming that all the light seeds that we are forming in cat, are located at the center of their hosting halo during their entire life, and so will be able to accrete the surrounding gas very efficiently. These three approximations likely lead to an overestimation of the AGN contribution to the total X-ray background. For this reason, we also compute the 21cm signal reducing the X-ray luminosity of AGNs to 10% of the value adopted by our reference model discussed above.

Fig. 9 shows the 21cm signal computed under the approximation described above compared with our reference signal when all sources (black thin line) and only stellar sources (red dashed line) are considered. Reducing the amount of X-rays produced by unobscured AGN, leaves an imprint on the signal of a wider and slightly deeper absorption feature ($\delta T_{\rm b}\sim-100$ mK at $z+1\sim 22-26$ ($\nu\sim 54-64$ MHz) and delays the transition to an emission signal to $z+1\sim 15.5$ ($\nu\sim 91$ MHz). Nevertheless, accreting BHs would still have an impact on the signal as the black thick line still shows a faster evolution than the red dashed line. This result may suggest that the approximations done within cat for the estimation of the X-ray background produced by the accreting BH seeds lead to an overestimation of the total heating of the neutral hydrogen when predicting the timing of the absorption feature. Finally, we notice that the model with the reduced X-rays from accreting BHs, is consistent also with the updated constraints recently provided by HERA (The HERA Collaboration et al. 2022).

In the previous section we did not consider the presence of a radio background (Fig. 8); here we will discuss the consistency between our signal and the claimed detection by the EDGES collaboration Bowman et al. (2018a).

The three key features of the absorption profile detected by EDGES are: (i) timing: $17\leq z+1\leq 20$ centered at $z+1\sim 18$, (ii) depth: $\delta T_{\rm b}=-500^{+200}_{-500}$ mK and (iii) shape: flat profile (all these features are reported in Fig. 8 with a red bar and a red shaded region representing the uncertainty). Our prediction shown in Fig. 8 only match the depth of the absorption profile with a maximum of $-500$ mK at $\nu\sim 63-74$ MHz. The most difficult of the EDGES constraints to match is the flat shape of the profile. The dashed line in Fig. 8 drops quite fast at high-$z$ (but not as fast as EDGES detection), it has a broad and almost constant peak between $19\leq z+1\leq 22.5$ and then the differential brightness temperature increases gently. This slow increase in the DBT is not consistent with EDGES claimed detection which instead shows a sharp increase. Finally, we note that, the effect of a radio background is to "wash out" all the "bursty" evolution originated from Pop III stars. Nevertheless, Pop III stars have a key role also in the computation of the radio background as they provide the BH seeds that, as soon as they start accreting gas, emit radiation in the radio band. Finally, we noticed that the timing of our absorption feature is predicted to be $\Delta z\sim 2$ earlier compared to EDGES. In the previous section we assessed that this is probably due to an overestimation of the X-ray background that erases the absorption feature faster.

In Fig. 10 we recomputed the 21cm signal in presence of a radio background with the reduced X-ray background from accreting BHs (dashed line). In this case, in order to achieve the depth of -500 mK we required $f_{\rm R}=30$. Once we account for a smaller X-ray background from unobscured AGNs, the resulting signal is more consistent with the timing of EDGES. However, the flat profile is still not fully achieved.

In conclusion, our model is not able to correctly reproduce all features of the EDGES signal. The main issue is with regards to the peculiar flat profile of EDGES that suggests a strong X-ray heating is responsible for a quick rise in the DBT. In our model BHs are important both for the X-ray and the radio emission. However the two contributions act in the opposite direction, so, while to match the depth of the signal, we need a strong radio background, this tends to wash out the X-ray contribution so we cannot achieve a sharp drop in the DBT. These results are also consistent with what has been found in Mebane et al. (2020).