Semi-analytical modelling of Pop. III star formation and metallicity evolution - II. Impact on 21 cm power spectrum

MNRAS 540, 483–497 (2025) https://doi.org/10.1093/mnras/staf699 Advance Access publication 2025 April 30 Semi-analytical modelling of Pop. III star formation and metallicity evolution – II. Impact on 21 cm power spectrum Emanuele M. Ventura , 1 , 2 ‹Yuxiang Qin , 2 , 3 Sreedhar Balu 1 , 2 , 4 and J. Stuart B. Wyithe 2 , 3 1 School of Physics, University of Melbourne, Parkville, VIC 3010, Australia 2 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Canberra, ACT 2611, Australia 3 Researc h Sc hool of Astr onomy and Astr ophysics, Austr alian National University, Canberr a, ACT 2611, Austr alia 4 Facultad de F ´ ısicas, Multidisciplinary Unit for Energy Science, Universidad de Sevilla, 41012 Seville, Spain Accepted 2025 April 27. Received 2025 April 11; in original form 2025 February 12 A B S T R A C T Simulating Population (Pop.) III star formation in mini-haloes in a large cosmological simulation is an extremely challenging task but it is crucial to estimate its impact on the 21 cm power spectrum. In this work, we develop a framework within the semi-analytical code MERAXES to estimate the radiative backgrounds from Pop. III stars needed for the computation of the 21 cm signal. We computed the 21 cm global signal and power spectrum for different Pop. III models varying star formation efficiency, initial mass function and specific X-ray luminosity per unit of star formation (L X /SFR). In all the models considered, we find Pop. III stars have little to no impact on the reionization history but significantly affect the thermal state of the intergalactic medium (IGM) due to the strong injection of X-ray photons from their remnants that heat the neutral IGM at z ≥15. This is reflected not only on the 21 cm sk y-av eraged global signal during the Cosmic Dawn but also on the 21 cm power spectrum at z ≤ 10 where models with strong Pop. III X-ray emission have larger power than models with no or mild Pop. III X-ray emission. We estimate observational uncertainties on the power spectrum using 21CMSENSE and find that models where Pop. III stars have a stronger X-ray emission than Pop. II are distinguishable from models with no or mild Pop. III X-ray emission with 1000 h observations of the upcoming SKA1-low. Key words: stars: Population III – galaxies: high-redshift – cosmology: dark ages, reionization, first stars. 1 INTRODUCTION When and where did Population III (Pop. III) stars form? What role did they play in the Cosmic Dawn and Epoch of Reionization (EoR)? And what is the best way to detect them? These questions remain open as no definitive observation of a mini-halo or Pop. III star has been reported. To gain insight, small size and high-resolution hydrodynamical simulations have been performed (e.g. Greif et al. 2011 ; Hirano et al. 2018 ; Chon, Omukai & Schneider 2021 ; Chon et al. 2022 ; Toyouchi et al. 2023 ; Sadanari et al. 2024 ) that suggest that metal-free (or poor) mini-haloes fa v our the formation of Pop. III stars with a more top-heavy initial mass function (IMF) and with lower star formation efficiencies than observed today. It is also thought that Pop. III star formation might occur down to the end of the EoR at z  6 in pristine metal free pockets of gas (e.g. Venditti et al. 2023 ). The variety and complexity of the processes involved in Pop. III star formation and the resolution required to keep track of the evolution of the gas particles, limits the size of these hydrodynamical simulations to ∼100 kpc. In order to mitigate this problem, semianalytical models that account for Pop. III star formation have been developed (e.g. Visbal, Bryan & Haiman 2020 ; Hegde & Furlanetto  E-mail: e[email protected].edu.au 2023 ; Liu et al. 2024 ). These models allow a statistical study of Pop. III star formation in mini-haloes out to scales of ∼10 Mpc. While these volumes start to investigate the chemical enrichment of the intergalactic medium (IGM) and the Pop. III/II transition, they are still too small to study the EoR as volumes of at least ∼200 Mpc are required (Iliev et al. 2014 ; Kaur, Gillet & Mesinger 2020 ; Balu et al. 2023a ). Observations and models are converging on a scenario where the Universe was completely ionized by z ∼5 . 3 (e.g. Fan, Carilli & Keating 2006 ; Ouchi et al. 2010 ; McGreer, Mesinger & D’Odorico 2015 ; Qin et al. 2021a , 2024 ; Bosman et al. 2022 ) with reionization likely driven by low-mass haloes (e.g. Kuhlen & Faucher-Gigu ` ere 2012 ; Qin et al. 2021b ; Mutch et al. 2024 ; Saxena et al. 2024 ). Ho we ver, the impact of Pop. III stars and mini-haloes on the EoR is unclear. Pop. III stars are likely to be the dominant contribution to the total star formation rate density (SFRD) at z > 15 −20 and, if their IMF is more top-heavy than the present day one, Pop. III could significantly contribute to the heating and the ionization of the IGM which determines the evolution and shape of the 21 cm signal (e.g. Qin et al. 2020 ; Gessey-Jones et al. 2022 ; Sartorio et al. 2023 ). The 21 cm signal represents our most promising tool to put constraints on the thermal state of the IGM during the Cosmic Dawn and EoR. Even though no confirmed detection has been reported so far, the first upper limits on the 21 cm power spectrum obtained © 2025 The Author(s). Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https:// creativecommons.org/ licenses/ by/ 4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 484 E. M. Ventura et al. MNRAS 540, 483–497 (2025) Table 1. Main free parameters for galaxy formation. Parameter Description Fiducial value αSF , II Pop. II Star formation efficiency 0.1 αSF , III Pop. III Star formation efficiency see Table 4 η0 Mass loading normalization 7.0 0 Supernova energy coupling normalization 1.5 Z crit Critical metallicity for Pop III star formation 10 −4 Z   norm Critical surface density of cold gas for star formation 0.37 M pc −2 Pop. III IMF Shape of Pop. III IMF Sal [1, 500] M  E PISN Energy from pair instability SN 10 52 erg E CCSN Energy from core collapse SN 10 51 erg Table 2. Main free parameters for reionization. Parameter Description Fiducial value f 0 esc , III Pop. III escape fraction normalization 0.14 f 0 esc , II as abo v e for Pop. II 0.14 αesc , III Pop. III escape fraction redshift scaling 0.2 αesc , II as abo v e for Pop. II 0.2 L X < 2keV , III / SFR Specific Pop. III X-ray luminosity per unit star formation see Table 4 L X < 2keV , II / SFR as abo v e for Pop. II 3.16 ×10 40 erg s −1 M −1 yr with HERA phase I strongly disfa v our cold reionization scenarios (HERA Collaboration 2023 ). In the last few years, the impact of Pop. III stars on the 21 cm signal has been studied using both analytical and semi-analytical models (e.g. Cohen et al. 2017 ; Chatterjee et al. 2020 ; Mebane, Mirocha & Furlanetto 2020 ). Ho we ver these models either did not compute reionization (e.g. Magg et al. 2022 ; Hegde & Furlanetto 2023 ; Cruz et al. 2024a ), focusing only on the absorption trough of the 21 cm signal occurring at z ∼13 −20, or used a very simple analytical approach to compute reionization (e.g. Ventura et al. 2023 ). On the other hand, Cohen et al. ( 2017 ), Qin et al. ( 2021a ), and Mu ˜ noz et al. ( 2022 ) used a simple analytical model for modelling Pop. III star formation but computed the reionization self-consistently. In this work, we o v ercome these challenges using a realistic Pop. III star formation and mini-halo model (Ventura et al. 2024 ) developed within the semi-analytical model MERAXES designed to self-consistently couple galaxy formation and reionization. While in Ventura et al. ( 2024 ) we ran this model on a small ( L = 10 h −1 cMpc) and high-resolution box, here we extend it to a significantly larger volume simulation ( L = 210 h −1 cMpc) enabling the study of cosmic reionization. Since at such large volumes we cannot directly resolve mini-haloes, we implemented scaling relations between the SFRD and the dark matter density field calibrated on the results from the small and high-resolution box discussed in Ventura et al. ( 2024 ). With this new model we are able to accurately follow the evolution of the radiative backgrounds relevant to the EoR and 21 cm signal (X-rays, Lymanα, ionizing UV, Lyman–Werner) and to disentangle the contribution of Pop. III star formation to the 21 cm global signal and power spectrum. Pop. III stars are expected to have a stronger impact at z ≥15 where they dominate star formation and ionization. Dif ferently from pre vious works who explored the differences in the 21 cm signal at Cosmic Dawn due to various Pop. III models, here we focus our attention on the residual signature of Pop. III on the 21 cm power spectrum at z ≤10 where the sensitivity of the Square Kilometer Array (SKA) is expected to be significantly better and a detection is more plausible. To achieve this, it is crucial to model both Pop. III star formation and reionization in a self-consistent framework. This study allows us to assess under which conditions an early heating of the IGM from Pop. III stars leaves a detectable imprint on the 21cm power spectrum at z ≤10. This paper is structured as follows: In Section 2 , we give a brief o v erview of Pop. III star formation in MERAXES . In Section 3 , we present the scaling relation between the SFRD in mini-haloes and the density field calibrated from the small and high-resolution box which is implemented in the large (210 h −1 cMpc) 3 box. In Section 4 , we discuss the impact of different Pop. III star formation models on the 21 cm power spectrum and in Section 5 we make forecasts on the observability of these power spectra with SKA. Finally, we summarize our main results and conclusions in Section 6 . Our simulations use the best-fitting parameters from the Planck Collaboration XIII ( 2016 ): h = 0.6751, m = 0.3121, b = 0.0490,  = 0.6879, σ8 = 0.8150, and n s = 0.9653. 2 POP. III GALAXIES IN MERAXES MERAXES 1 is a semi-analytical model designed to study the interplay between galaxy formation and reionization (Mutch et al. 2016 ; Qin et al. 2017 ; Qiu et al. 2019 ; Ventura et al. 2024 ). MERAXES includes a number of free parameters that are calibrated against observations (see Tables 1 and 2 ). Values in Table 1 are calibrated against observed luminosity functions and stellar mass functions at z ∼5 −8 while those in Table 2 are calibrated against constraints on the neutral hydrogen fraction, ionizing emissivity and the Thomson scattering optical depth τe from Planck Collaboration VI ( 2020 ). The most recent version of MERAXES (Ventura et al. 2024 , V24) includes Pop. III star formation and mini-halo physics. As shown in V24, the Pop. III parameters with the largest impact are the star formation efficiency αSF , III and the shape of the IMF. The latter has a large impact on galaxy evolution as it determines the strength of the feedback and the emission properties of the Pop. III stellar population. In the following sections, we quickly summarize the main features of MERAXES rele v ant for this work. 1 https:// github.com/ meraxes-devs/ meraxes Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 21 cm PS with Pop. III star formation 485 MNRAS 540, 483–497 (2025) 2.1 Galaxy formation MERAXES post-processes the output of an N -body dark matter only simulation, reading the spatial and physical information of dark matter haloes and computing the baryonic physics of galaxy formation. In particular processes included are: (i) gas infall on to dark matter haloes, (ii) radiative cooling of the infalling gas, (iii) star formation, and (iv) supernova and active galactic nuclei feedback. In V24, the cooling prescriptions were updated to account for H 2 cooling (the main cooling channel in mini-haloes) and a more detailed metal enrichment model to keep track of the metallicity of each gas reservoir in a halo (which is crucial to distinguish between Pop. III and Pop. II star formation episodes). We also account for the effects of both baryon-dark matter streaming velocity and H 2 photodissociation by the Lyman–Werner background which increases the minimum mass of a mini-halo capable of hosting stars (e.g. Schauer et al. 2021 ). Our model accounts for spatial variations only for the LW background, while for the relativ e v elocity we assume a mean value throughout the entire box. For this work, we updated MERAXES by adding the effect of H 2 self-shielding which counteracts the H 2 photo-dissociation, increasing the Pop. III SFRD by up to one order of magnitude at z ∼10 (see e.g. Feathers et al. 2024 ). We discuss the details of the implementation and impact of H 2 self-shielding in MERAXES in the Appendix A . We refer the reader to Mutch et al. ( 2016 ) for a more detailed explanation of the main architecture of MERAXES , Qiu et al. ( 2019 ) for the supernova model and V24 for the mini-halo model. The main free parameters that regulate the galaxy formation in MERAXES are summarized in Table 1 . The Pop. II related ones are taken from Balu et al. ( 2023a ) where MERAXES was run on a cosmological volume of L = 210 h −1 cMpc resolving all atomic cooling haloes and calibrated in order to match the observed ultraviolet (UV) luminosity functions at z ∼4 −7 (the agreement holds up to z ∼13 as shown in Qin, Balu & Wyithe 2023 ) and the stellar mass functions at z ∼5 −8. Gi ven the lack of observ ations of Pop. III stars, the Pop. III parameters are largely unconstrained. The fiducial values adopted in this work are taken from V24 and their values are suggested from hydrodynamical simulations (e.g. Chon et al. 2021 ). 2.2 Reionization and radiati v e backgrounds Together with galaxy formation, MERAXES self-consistently computes the reionization and thermal evolution of the IGM using a modified version of the seminumerical code 21 cm FAST (Mesinger, Furlanetto & Cen 2011 ). In this work, we compute the backgrounds rele v ant for the computation of the 21 cm signal: the UV ionizing, Xrays, L ymanαand L yman–Werner (LW). The first is crucial to study the evolution of the reionization, while the X-ray and the Lymanα backgrounds determine the thermal state of the IGM. In particular the X-ray background is likely to be the dominant contribution to the heating of the IGM once the first galaxies form and the lymanα background is responsible for the coupling between the kinetic and the spin temperature of the neutral hydrogen. The LW background does not directly affect the IGM temperature, but determines whether or not mini-haloes have enough molecular hydrogen to cool the gas and form Pop. III stars. Hereafter, we briefly summarize the key quantities that determine the evolution of these backgrounds. For a more detailed explanation on the implementations of these backgrounds, we refer the reader to Balu et al. ( 2023a ) for the UV, Lymanαand X-ray and to V24 for the LW. The ionizing background is mostly dependent on the SFRD, the average number of ionizing photons per stellar baryon N γand the escape fraction of the UV photons f esc . The second quantity is mostly determined by the IMF: for Pop. II stars we adopt a Kroupa IMF which leads to N γ∼6000. Since the Pop. III IMF is a free parameter in our model, N γ, III is computed from the IMF adopted using the Pop. III stellar spectra from (Raiter, Schaerer & Fosbury 2010 ; N γ, III ∼ 20 000 −70 000). f esc is tuned to reproduce the EoR histories in agreement with observations. As per Balu et al. ( 2023a ), we adopt a redshift-dependent escape fraction defined as followed: f esc = f 0 esc 1 + z 6 αesc . (1) X-ray emission is mostly associated with high mass X-ray binaries (HMXB) and its contribution is proportional to the SFRD. In this work we use the widely adopted approximation for the comoving X-ray specific emissivity (erg s −1 Mpc −3 ) x ∝ L X / SFR ×SFRD. Finally, we need to account for the fact that only photons with an energy below 2 keV (soft X-rays) are able to heat the IGM. As a result the main free parameter that regulates X-ray emissivity is the soft X-ray luminosity per unit star formation L X < 2keV /SFR. For Pop. II stars, this quantity is estimated from theoretical studies of emission spectra of HMXBs in low-metallicity environments (e.g. Fragos et al. 2013 ; Das et al. 2017 ; Madau & Fragos 2017 ; Qin et al. 2020 ; Kaur et al. 2022 ). For Pop. III stars there are no observational constraints as this quantity depends on the unknown Pop. III IMF. Recently, Sartorio et al. ( 2023 ) estimated the L X /SFR for Pop. III stars and found that for more top-heavy IMFs this quantity can be up to two orders of magnitude higher than the Pop. II value. We highlight that in this work we separately compute the backgrounds from Pop. III and Pop. II stars due to the different spectra, properties, and star formation rate density of the two distinct populations. 2.3 21 cm physics Using the radiative backgrounds computed from the galaxy population in MERAXES , we can estimate the 21 cm signal. We encourage the reader to see Furlanetto, Peng Oh & Briggs ( 2006 ), Morales & Wyithe ( 2010 ), Pritchard & Loeb ( 2012 ), and Liu & Shaw ( 2020 ) for re vie ws on the topic. Hereafter, we only summarize the key equations used in this work (see also Balu et al. 2023a ). We start with the 21 cm brightness temperature field ( δT b ) which measures the deviation of the spin temperature of the neutral hydrogen ( T S ) from the cosmic background T γ(i.e. the CMB). This is given by (Furlanetto et al. 2006 ): δT b = T S −T γ 1 + z (1 −e −τν0 ) ≃ 27 x HI (1 + δnl ) H d v r / d r + H 1 −T γ T S  ×1 + z 10 0 . 15 m h 2 b h 2 0 . 023 mK, (2) where τν0 is the optical depth at the 21 cm transition frequency ν0 , x H is the neutral hydrogen fraction, 1 + δnl is the density contrast in the dark matter field, H ( z) is the Hubble parameter at the redshift z, and d v r / d r is the radial deri v ati ve of the line-of-sight component of the peculiar velocity. Once the cosmological model (Planck Collaboration XIII 2016 ) and the velocity and density field (from the N -body simulation) are fixed, δT b is determined by the ionization and the spin temperature fields. The latter quantifies the population ratio of the two H I hyperfine energy levels and is sensitive to the thermal state (i.e. the kinetic temperature T K ) of the gas as Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 486 E. M. Ventura et al. MNRAS 540, 483–497 (2025) Table 3. Simulation parameters. Label Box side (cMpc) Mass resolution (M ) Pixel side L10 10 h −1 4.71 ×10 5 0.2 h −1 L210 210 h −1 3.16 ×10 7 0.2 h −1 follows: T −1 S = T −1 γ+ x αT −1 α+ x c T −1 K x α+ x c + 1 , (3) where T αis the colour temperature which we take equal to T K while x αand x c are the Lymanαand collisional coupling coefficients, respecti vely. These coef ficients quantify the strength of the processes (resonant scattering of Lymanαphotons, Wouthuysen 1952 , and collisions with free electrons) that drive the spin temperature towards the kinetic temperature (when x α+ x c >> 1, T S ∼T K otherwise T S ∼T γ). T K is sensitive to the adiabatic cooling and to all the processes able to heat up (or cool) the IGM with the most dominant coming from the X-ray emission. Hence, in this work we will consider only the X-ray heating neglecting the other source of heating such as primordial magnetic fields (Minoda, Tashiro & Takahashi 2019 ; Bera, Datta & Samui 2020 ; Cruz et al. 2024b ), Lymanα (Ciardi, Salvaterra & Di Matteo 2010 ; Mittal & Kulkarni 2021 ; Reis, Fialkov & Barkana 2021 ), shocks (Furlanetto & Loeb 2004 ; Gnedin & Shaver 2004 ; Ma et al. 2021 ), cosmic rays (Bera, Samui & Datta 2023 ), early accreting black holes (Mebane et al. 2020 ; Ventura et al. 2023 ) and decaying or annihilating dark matter (Liu & Slatyer 2018 ; Sun et al. 2023 ; Facchinetti et al. 2024 ; Hou & Mack 2025 ). Ultimately, the evolution of T S during the Cosmic Dawn and the EoR is mostly determined by the Lymanα(for the coupling between T S and T K ) and X-ray flux. Using equation ( 2 ) we can estimate both the all-sky averaged global signal and its fluctuations (i.e. the power spectrum). In this work, we will often use the reduced power spectrum  2 21 ( k ) = k 3 /2 π2 P 21 (k) unless otherwise stated. 3 PUTTING POP. III GALAXIES AND MINI-HALOES IN A LARGE-SCALE SIMULATION In this section, we present a no v el approach that enables us to efficiently estimate the SFRD from mini-haloes in a large box (that does not directly resolve these objects) using results from a small, high-resolution box. In Table 3 , we summarize the key parameters for both the small (L10) and large (L210) box. 3.1 Calibrating scaling relations from the L10 box Our starting point is the small ( L = 10 h −1 cMpc) high-resolution (halo mass resolution of M ∼4 . 7 ×10 5 M ) simulation used in V24. When building a scaling relation between the star formation rate (SFR) and other physical quantity, the first obvious choice is the dark matter density field δ. For instance Mu ˜ noz ( 2023 ) showed, as a first-order approximation, SFR scales as e δR where δR is the density field smoothed o v er a certain radius R and this relationship works quite well for δ∼0 and large R ( ≥3 Mpc). To link our SFR in mini-haloes with the density field, we first compute density, δ( x , z), and SFR grids for both Pop. III, SFR MC , III ( x , z), and Pop. II, SFR MC , II ( x , z), in the L10 box using the same grid resolution used in Balu et al. ( 2023a ) to compute reionization ( L pixel ∼0 . 3 cMpc) and accounting for the SFR within mini-haloes. We highlight that our L pixel is quite small compared to the smoothing radius R adopted by Mu ˜ noz ( 2023 ), hence we expect a significant scatter in the abo v e relation. We also split the contribution between Pop. III and Pop. II stars (a chemically enriched mini-halo will form Pop. II stars). In the left panel of Fig. 1 we show our Pop. III SFR distribution as a function of the o v erdensity δand with the gre y line we highlight the SFR ∝ e δR relation as in Mu ˜ noz ( 2023 ). Despite the significant scatter for the reasons outlined abo v e ( σ∼0 . 65), the analytical approximation agrees with our results. The results shown hereafter are obtained from our fiducial simulation in V24. We start by investigating the distribution function of log 10 (SFR) at a fixed overdensity δand redshift log 10 (SFR( δ, z)) finding that it follows a Gaussian distribution (or lognormal in the linear space). We show results for log 10 (SFR III ) and selected values of δin the smaller panels in Fig. 1 . Hence, we can write:  ( log 10 (SFR MC | δ, z)) = A ( δ, z) e ( log 10 (SFR) −log 10 ( SFR ( δ,z))) 2 2 σ( δ,z) 2 , (4) where the normalization A , the mean log 10 ( SFR ) and the standard deviation σall depend on the o v erdensity and redshift. The normalization is defined as the ratio between the number of star forming pixels and the total number of pixels. We found the best-fitting parameters for each δ(grouped in bins of width = 0.1) and snapshot of the simulation. We tested whether the log 10 (SFR) distribution function is indeed Gaussian by conducting a K-S test. P -values are calculated for each cell with SFR > 0 and taking a predefined significance level of 0.05 below which the null hypothesis will be rejected. Results are shown in Fig. 2 for both SFR III and SFR II . P - values al w ays exceed the significance level for both Pop. III (left panel) and Pop. II (right panel) SFR suggesting that the Gaussian distribution reproduces  ( log 10 (SFR MC | δ, z)) both in the Pop. III and Pop. II cases. As expected, we see that there are far more Pop. III star forming pixels than Pop. II ones as mini-haloes are more likely to form Pop. III stars. The next step is to study how the mean, standard deviation, and normalization evolve with δand z. In Fig. 3 , we show the redshift evolution of these parameters for δ= 0.5 (black), 1.0 (grey), 1.5 (purple), 2.0 (red), 2.5 (green), and 3.0 (blue). SFR MC , III exhibits an almost constant trend in redshift and a correlation with δ(higher δresults in higher SFR MC , III ). This demonstrates that SFR MC , III is mostly determined by the number of Pop. III star forming haloes in a pixel, which is higher for more o v erdense re gions. Since Pop. III star formation episodes in mini-haloes are often the first episode of star formation experienced by a galaxy, it is not impacted by supernova feedback 2 so the Pop. III SFR is almost constant at all redshift. This also explains why σis constant for all z and δ( σMC , III ∼0 . 65). The parameter that is more sensitive to both δand z is the normalization. F or v ery o v erdense re gions ( δ≥2) it is almost one, meaning that almost all the o v erdense pix els host Pop. III star formation minihaloes. For lower δthere is also an evolution in z as, with cosmic time, lower density regions will host a larger number of mini-haloes abo v e the minimum mass for star formation. We repeated the same analysis for Pop. II star forming pixels (see Fig. 4 ). In this case, the evolution is more scattered as Pop. II star formation episodes 2 Even though supernova feedback can be neglected, this is not true for the Lyman–Werner background that halts star formation in mini-haloes. For this reason we add a further condition that if a pixel is irradiated by a LW flux J LW ≥J crit abo v e a critical threshold defined as M crit, MC = M ato (the atomic cooling threshold),  (log 10 (SFR MC )) = 0. M ato is the virial mass correspondent to a halo with a virial temperature T vir = 10 4 K. While M crit, MC is defined in equation ( A2 ). Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 21 cm PS with Pop. III star formation 487 MNRAS 540, 483–497 (2025) Figure 1. Left panel shows the density distribution of Pop. III star formation rate in mini-haloes (M yr −1 in logarithm scale) versus the dark matter o v erdensity δfor each pixel at z = 15. The thick grey line shows the analytical fit SFR ∝ e δsimilar to the one adopted by Mu ˜ noz ( 2023 ) together with the 1 σdeviation (thin lines). For different values of δ(highlighted with the black rectangles) we show the distribution of Pop. III SFR in mini-haloes (in logarithm scale) together with the best Gaussian fit. Figure 2. P -value distribution of K-S tests conducted on the Pop. III (left) and Pop. II (right) star forming pixels. Vertical line highlights the significant level of 0.05. are not the first star forming episodes within a galaxy and so will be affected by both mechanical and chemical feedback from the previous history of the galaxy. This is also demonstrated by the larger standard deviation ( σMC , II ∼0 . 8). The average value of SFR MC , II is ∼1 order of magnitude higher than for Pop. III. This reflects the higher Pop. II star formation efficiency . Finally , A II is al w ays smaller than A III showing that it is less likely for a mini-halo to form Pop. II stars. We test this parametrization on the small box by estimating the mini-halo contribution to the Pop. III and Pop. II SFRD from the matter density field. To do this, we read the density grid at each z and for each pixel assign a value of Pop. III SFR drawn randomly from  ( log 10 (SFR MC , III )) using δof the pixel. We repeat the same procedure for SFR MC , II adding the constraint that in order to have SFR MC , II > 0 that pixel needs to already have experienced a Pop. III star formation episode. This latter condition ensures a more realistic enrichment model (a pixel cannot have Pop. II star formation if it has not previously hosted Pop. III star formation). We ran twenty different realizations and for each realization estimated the SFRD MC , III and SFRD MC , II and compared with results of the simulations. Results are shown in Fig. 5 . In the upper panels, we show the Pop. III (left) and Pop. II (right) SFRD from MERAXES output (black line) and from each realization using the method outlined abo v e (c yan shaded lines). In the lower panels, we show the ratio between the average of the 20 realizations and the true SFRD from MERAXES . All the realizations are in reasonable agreement with the data and the ratio is al w ays lower than ∼10 per cent for both populations. This demonstrates the validity of the method outlined abo v e. The main advantage of this method is that it allows estimation of the SFRD from mini-haloes in a simulation where these are not directly resolv ed. P arametrizing SFRD with a Gaussian distribution enables us to account for stochastic star formation (different pixels with same δcan have different SFR), without losing the correlation with the matter density field ( SFR , σand A all depend on δ). This method can be applied as long as the density field from both the lowand the high-resolution simulation share the same properties (mean, standard deviation) and it can be calibrated for any choice of parameters. 3.2 Applying scaling relations to the L210 box We can now apply the methodology outlined in the previous section to the L210 box that can only resolve the atomic cooling haloes by reading the density field and applying  ( log 10 (SFR MC )) calibrated on different models. In Fig. 6 , the mini-halo contribution to the Pop. III SFRD shows a good agreement between the two different simulations. We highlight that our scaling relations do not explicitly depend on the Lyman–Werner background (except for the regions that are strongly irradiated by LW flux for which  (log 10 (SFR MC )) = 0). This implicitly assumes that both the L10 and L210 box have similar LW backgrounds. We verify this assumption by computing the average LW background and LW maps in both simulations (see Fig. 7 ). In the bottom right panel, we show the redshift evolution of the LW background (in units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 ) in the L10 (black line) and in the L210 box (red line). The two lines share similar trends showing that both simulations have a similar average LW background. The left and top right panels show the 2D projections of the LW field in the large and small box. These maps illustrate that the LW background is roughly uniform (as expected given that the mean free path of LW photons is ∼100 Mpc). These Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 488 E. M. Ventura et al. MNRAS 540, 483–497 (2025) Figure 3. Redshift evolution of log 10 ( SFR ) (top panel), σ(mid), and A (bottom) for Pop. III for δ= 0.5 (black), 1.0 (grey), 1.5 (purple), 2.0 (red), 2.5 (green), and 3.0 (blue). Figure 4. Same as Fig. 3 for Pop. II. plots demonstrate that the LW fields in the small and large box are indeed comparable showing that the estimation of the star formation in mini-haloes with the scaling relations accounts for the radiative feedback. 3 Differently from what has been done by Hazlett et al. ( 2024 ) who calibrated a semi-analytical model to the Reinassance simulation in order to account for Pop. III star formation by adding the previous star formation history to each atomic cooling galaxy resolved in the simulation, the methodology described in this section allows us 3 In this discussion, we neglected the UV photo-ionizing feedback. This is justified by the fact that this effect impacts more the low-mass atomic cooling haloes at z  10 rather than mini-haloes. Figure 5. (top) Pop. III (left) and Pop. II (right) SFRD versus z from MERAXES (black) and estimated from the density field (cyan). (bottom) ratio between the average of the 20 realizations for the SFRD estimated from the density field and the SFRD from MERAXES . Figure 6. SFRD MC , III versus z from the L10 (solid) and L210 (dotted) box for two different Pop. III star formation models (see more details in text and Table 4 ). to estimate only the total SFRD occurring in mini-haloes within a certain pixel of the simulation. In this work, we focus on three models of Pop. III star formation by varying three main parameters: the star formation efficiency, the IMF, and the specific Pop. III X-ray luminosity per unit star formation while all the other free parameters (e.g. the escape fraction) are fixed at the fiducial value (see Tables 1 and 2 ). We chose to focus only on these three parameters since these have the strongest impact on both the evolution of the Pop. III SFRD and the amount of UV and X-ray photons emitted. The values chosen for the specific Pop. III X-ray luminosity per unit star formation are similar to those found in Sartorio et al. ( 2023 ) for the different IMFs explored in their work. Hereafter, we analyse three different Pop. III models designed to have the minimum, intermediate, and maximum impact from Pop. III star formation in mini-haloes, each model is summarized in Table 4 . The IMFs considered in this model are a Salpeter between 1 and 500 solar masses and a lognormal IMF centred at 60 M (see V24 for more details). We highlight that in our extreme Pop. III model we enhance star formation efficiency, L X /SFR and the top-heaviness of the IMF at the same time. Each of these parameters has a different impact on the evolution of the 21 cm signal (see Appendix B for a more detailed discussion of how each Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 21 cm PS with Pop. III star formation 489 MNRAS 540, 483–497 (2025) Figure 7. Left panel shows the 2D projections of the LW background (units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 in the L210 box at z = 15. Top right panel shows the same map but in the L10 simulation. Bottom right panel shows the redshift evolution of the average LW background (same units as abo v e) in both the L10 (black line) and L210 (red line) simulation. Table 4. Pop. III model parameters. Label IMF type a αSF , III L X < 2keV , II / SFR Weak Pop. III Salpeter 0.008 3 ×10 40 Moderate Pop. III Salpeter 0.008 3 ×10 41 Extreme Pop. III logE 0.08 3 ×10 42 High SFE Salpeter 0.08 3 ×10 40 LogE logE 0.008 3 ×10 40 a See Table 2 in V24 for the details. Pop. III parameter changes the evolution of the 21 cm global signal and power spectrum). 4 IMPACT OF POP. III STAR FORMATION ON 21 CM PHYSICS We can now estimate how different Pop. III star formation models in MERAXES affect the 21 cm signal. We started by verifying that, after introducing the additional Pop. III contribution to the fiducial Pop. II only model (Balu et al. 2023a ) we still obtain reionization histories consistent with the observational constraints on the Thomson scattering optical depth τe (Fig. 8 ) and ¯ x HI (Fig. 9 ). The reionization histories from the models with Pop. III stars are only slightly modified and this negligible contribution comes from the secondary ionizations from X-rays. This is expected given that at z ≤15 the Pop. III SFRD is at least one order of magnitude lower than Pop. II and their main contribution is expected from the X-ray emission rather than the UV. We compute the sk y-av eraged 21 cm global signal (see Fig. 10 ) without (black line) and with (gre y, c yan, and red line for weak, moderate and extreme Pop. III respectively) Pop. III star formation. Figure 8. Integrated Thomson scattering optical depth τe computed for model with weak (grey), moderate (cyan), extreme (red) Pop. III, and Balu et al. ( 2023a ) (black). The green curve and shaded region show the measurement of τe from the Planck 2018 collaboration (Planck Collaboration VI 2020 ). As expected, introducing a Pop. III population with the same Xray properties as Pop. II ones (i.e. weak Pop. III), simply shifts the absorption through to earlier epochs in virtue of the stronger coupling at higherz (see also Hegde & Furlanetto 2023 ; Ventura et al. 2023 ). Ho we v er, if Pop. III stars hav e a stronger X-ray emission (i.e. moderate and extreme models) as suggested by Sartorio et al. ( 2023 ), the absorption signal is quickly suppressed turning into an emission signal as early as z ∼18 for the extreme Pop. III model and at z ∼13 for the moderate Pop. III one. We note that a similar result has been found by a contemporaneous work by Gessey-Jones et al. ( 2025 ) who found an analogous variation in the timing ( z ∼3) and Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 490 E. M. Ventura et al. MNRAS 540, 483–497 (2025) Figure 9. Constraints on the reionization history (neutral hydrogen fraction vesus z ) for model with weak (grey), moderate (cyan), extreme (red) Pop. III, and Balu et al. ( 2023a ) (black). The observational data are from analyses of dark pixels (McGreer et al. 2015 ; Jin et al. 2023 ), damping-wing absorption in quasar spectra (Ba ˜ nados et al. 2018 ; Davies et al. 2018 ; Wang et al. 2020 ; Greig et al. 2022 ; Spina et al. 2024 ) and equi v alent width measurements (Mesinger et al. 2014 ; Hoag et al. 2019 ; Mason et al. 2019 ; Jung et al. 2020 ; Whitler et al. 2020 ). Figure 10. Effect of Pop. III star formation on the 21 cm global signal ( δT b versus z). Pop. III models with small X-ray heating cause a stronger absorption at earlier times, while having a stronger Pop. III X-ray heating causes the signal to be seen in emission earlier. Colour coding as in the previous figures. Brown and yellow dashed lines are taken from Gessey-Jones et al. ( 2025 ) for a Salpeter ( Sal ) and flat ( Int-0 ) IMF. depth ( δT b ∼50 mK) of the absorption trough when considering a stronger X-ray contribution from Pop. III stars (in their model the L X /SFR is self-consistently modelled from the IMF so that the difference between their Int-0 and Sal model is of 2 orders of magnitudes.) As shown in Fig. 11 earlier coupling and heating from Pop. III impacts the 21 cm power spectrum both at the large and small scales. The models considered here produce 21 cm signals that are different not only during the coupling and heating epoch ( z ∼10 −20), but also at lower redshift when the reionization is in progress. The impact is stronger at smaller scales where models with a stronger heating exhibit a larger power spectrum at z ∼7 −10. This is of crucial importance as current and upcoming facilities are improving their observations at z ≤10 (see more discussion in Section 5 ). We compare our results with the ones in Mu ˜ noz et al. ( 2022 ) who studied ho w dif ferent Pop. III parameters impact the 21 cm signal using 21cm FAST . In the left panel of Fig. 11 , the yellow and brown dashed lines are obtained by changing the specific Pop. III X-ray luminosity by a factor of 9. The global trends are similar with our models with stronger Pop. III X-ray emission (dashed yellow line and red/cyan lines) showing an earlier and weaker first peak compared to low X-ray models (dashed brown and grey lines). In their model, different Pop. III X-ray properties strongly impact the position and amplitude also of the second peak. In our models instead, the position of the second peak is only slightly anticipated in the strong X-ray models. This difference is likely due to the fact that our second peak occurs at much lower redshift ( z ∼10) when most of the emission comes from Pop. II stars. Finally, differently from Mu ˜ noz et al. ( 2022 ) at z ∼10 simulations with a stronger Pop. III X-ray heating have a significantly larger power spectrum compared to weak X-ray models (a similar result has been found also by Gessey-Jones et al. 2025 ). This analysis shows that Pop. III star formation not only impacts the 21cm global signal during the cosmic dawn as previously assessed (e.g. Qin et al. 2020 ; Gessey-Jones et al. 2022 ; Hegde & Furlanetto 2023 ; Ventura et al. 2023 ; Cruz et al. 2024a ) but also that an early ( z ≥15) heating of the IGM provided by this population leaves a strong signature on the power spectrum during the EoR ( z ≤10). The impact on the power spectrum is stronger for models that have a large Pop. III X-ray emissivity (i.e. moderate and extreme Pop. III models) while models with a milder X-ray emission (i.e. weak Pop. III) have a stronger effect on the global signal. Ultimately, this tells us that the power spectrum during the EoR can be used to disentangle different heating models and potentially constrain the properties of Pop. III stars. We highlight that here we considered only the X-ray emission from stellar remnants. While there might be other sources that significantly heat the IGM at z ≥15 (see discussion at end of Section 2.3 ), the other effects are likely to be either subdominant, still related to Pop. III stars (e.g. cosmic rays) or dominant only at the dark-ages (e.g. dark matter annihilation). Finally, we note that in this work we did not include the effect of the X-ray feedback on Pop. III star formation. As noted by Ricotti ( 2016 ) and Park, Ricotti & Sugimura ( 2021 ), X-rays have both a positive and ne gativ e effect on Pop. III star formation as the y heat the gas and increase the electron fraction. The heating makes gas accretion more difficult (hence delaying star formation) and free electrons promote the formation of H 2 making the molecular cooling more efficient. In presence of a strong X-ray background this latter effect is dominant at 10  z  20 (see e.g. fig. 9 in Hegde & Furlanetto 2023 ). Hence, including the X-ray feedback would likely increase the Pop. III SFRD making the impact of Pop. III stars even stronger than what predicted here especially for the moderate and extreme models. 5 OBSERVABILITY WITH SKA In the previous section, we showed that an early heating of the IGM significantly affects the 21cm power spectrum during the EoR. Now we want to investigate how our models compare with the currently available upper limits and whether their differences will be observable with SKA. In this section, we consider the power spectrum at z ≤10 as the current and upcoming interferometers are significantly less sensitive at higher redshift. To better appreciate this, in Fig. 12 we show the power spectrum noise (mK 2 ) with SKA as a function of redshift at k ∼0 . 2 and 0.9 cMpc −1 assuming a 1000 h (solid lines) and 180 h (dashed lines) observations with SKA. At higher redshift the noise steadily increases and already at z ≥10 the noise is of the order of 10s to 100s mK 2 with 1000 h observation. First, we compare our predictions with current upper limits from a number of facilities including the Murchison Widefield Array (MWA), LOw-Frequency ARray (LOFAR), Giant Metrewave Radio Telescope (GMRT), Precision Array for Probing the Epoch of Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025 21 cm PS with Pop. III star formation 491 MNRAS 540, 483–497 (2025) Figure 11. Effect of Pop. III star formation on the 21 cm power spectrum (  21 versus z ) at large (k ∼0 . 1 Mpc −1 ) and small (k ∼0 . 9 Mpc −1 ) scales. Colour coding as in the pre vious figures. Bro wn and yello w dashed lines are taken from Mu ˜ noz et al. ( 2022 ) for k = 0 . 23 Mpc −1 assuming a weak and strong X-ray luminosity per unit of Pop. III SFR (bottom right panel of fig. 17). Figure 12. 21 cm power spectrum sensitivity as a function of redshift k ∼0 . 2 (black line) and 0.9 cMpc −1 (orange line) assuming a 1000 (solid lines) and 180 h (dashed lines) observation with SKA. Reionization (PAPER) and Hydrogen Epoch of Reionization Array (HERA) (we focus at z = 7 −10 where most of the measurements have been taken). In Fig. 13 , we show the 21 cm power spectrum at z = 10, 9, 8, and 7 for all the four models together with the available upper limits (see label). All our models are below these upper limits. Ho we v er, the moderate and e xtreme Pop. III models are closer to the current HERA constraints at z = 8 suggesting that soon these models can be either detected or ruled out. In our model the main effect of different Pop. III models is to change the timing and amplitude of the peaks in the 21 cm power spectrum rather than introducing a specific feature. We note that we did not account for spatial variations in the velocity acoustic oscillations (see Section 2.1 ) which would introduce wiggles in the  2 21 at scales k ∼0 . 1 Mpc −1 (Cruz et al. 2024a ). While this effect can be important when focusing on the 21 cm power spectrum during the coupling and heating epochs, these fluctuations are quickly washed out at z  13 (see Fig. 13 and section 6C in Cruz et al. 2024a ). We next consider the upcoming SKA. In order to estimate the observability of the four models analysed so far, we perform a similar analysis as in Balu, Greig & Wyithe ( 2023b ) that we briefly summarize hereafter. The sensitivity of a radio interferometer is mostly regulated by the thermal noise (  N ) and the cosmic variance with the former dominating the noise at small scales and the latter at large scales. The thermal noise is related to the bandwidths of the instrument, beam factor (see Parsons et al. 2014 ), the integration time of the mode k and the temperature of the system (given as the sum of the sky and receiver temperature) (Morales 2005 ; McQuinn et al. 2006 ; Parsons et al. 2012 ). We can hence write the total noise by summing these two components: 1 σ[  2 21 ( k)] 2 =  i 1  2 N +  2 21 2 . (5) By doing so, we are ef fecti vely assuming that the errors are Gaussian distributed, which is reasonable for the rele v ant scales in this work (Qin et al. 2021a ; Prelogovi ´ c & Mesinger 2023 ). Finally, a 21 cm detection is heavily limited by the ability of removing the foregrounds. We used the python package 21CMSENSE 4 (Pober et al. 2013 , 2014 ) which, given the specifics of an interferometer, a mock 21 cm power spectrum and an observational campaign, computes the interferometer sensitivity to the 21 cm power spectrum under different assumptions of foreground removals. We used the assumption ‘moderate’ foreground removals and we focused on the first phase of SKA (i.e. SKA1-low), in particular we included the stations in the ‘Central Area’ of the SKA1-low, 5 resulting in 296 stations of diameter 35 m distributed across a circular area with 1.7 km diameter. We assumed two different observational campaigns: 4 https:// github.com/ rasg-affiliates/ 21cmSense 5 See the official SKA1 System Baseline Design document in https://www. skao.int/en for further details. Downloaded from https://academic.oup.com/mnras/article/540/1/483/8123416 by guest on 08 July 2025