Impact of the Primordial Stellar Initial Mass Function on the 21-cm Signal

Properties of the first generation of stars (Pop III), such as their initial mass function (IMF), are poorly constrained by observations and have yet to converge between simulations. The cosmological 21-cm signal of neutral hydrogen is predicted to be sensitive to Lyman-band photons produced by these stars, thus providing a unique way to probe the first stellar population. In this paper, we investigate the impacts of the Pop III IMF on the cosmic dawn 21-cm signal via the Wouthuysen-Field effect, Lyman-Werner feedback, Ly-alpha heating, and CMB heating. We calculate the emission spectra of star-forming halos for different IMFs by integrating over individual metal-free stellar spectra, computed from a set of stellar evolution histories and stellar atmospheres, and taking into account variability of the spectra with stellar age. Through this study, we therefore relax two common assumptions: that the zero-age main sequence emission rate of a Pop III star is representative of its lifetime mean emission rate, and that Pop III emission can be treated as instantaneous. Exploring a bottom-heavy, a top-heavy, and intermediate IMFs, we show that variations in the 21-cm signal are driven by stars lighter than 20 solar masses. For the explored models we find maximum relative differences of 59% in the cosmic dawn global 21-cm signal, and 131% between power spectra. Although this impact is modest, precise modelling of the first stars and their evolution is necessary for accurate prediction and interpretation of the 21-cm signal.


INTRODUCTION
In 2018 the EDGES collaboration reported a detection of a skyaveraged (global) radio signal at 78 MHz which could be interpreted as the cosmological 21-cm signal of neutral hydrogen from ∼ 17 (Bowman et al. 2018). If this signal truly is of cosmological origin, it provides evidence of star formation at cosmic dawn (e.g. Schauer et al. 2019;Mirocha & Furlanetto 2019;Fialkov & Barkana 2019;Reis et al. 2020;Mebane et al. 2020). Debate remains, however, over the true nature of the detected signal (Hills et al. 2018;Singh & Subrahmanyan 2019;Bradley et al. 2019;Sims & Pober 2020;Singh et al. 2022). A multitude of low-frequency radio telescopes (Singh et al. 2018;Philip et al. 2019;Mertens et al. 2020;Garsden et al. 2021;Abdurashidova et al. 2022;Mertens et al. 2021;Rahimi et al. 2021;Singh et al. 2022) are attempting to detect or bound the skyaveraged 21-cm signal or the power spectrum of its variations in the hope of supporting or refuting the EDGES measurement.
The numerous efforts to detect the 21-cm signal are reflective of the ★ E-mail: tg400@cam.ac.uk fact that the signal is expected to be a rich source of information about the Universe between recombination and reionization (e.g. Madau et al. 1997;Furlanetto et al. 2006;Pritchard & Loeb 2012;Barkana 2016;Mesinger 2019). The insights into the evolution and variation in the properties of the intergalactic medium (IGM) provided by the 21-cm signal will, in turn, probe physical processes that impact the IGM during this period, such as the Wouthuysen-Field coupling (WF, Wouthuysen 1952;Field 1958), heating by X-ray sources (e.g. Fialkov et al. 2014;Pacucci et al. 2014), ionization (e.g. Park et al. 2020), and coupling between dark matter and gas (e.g. Barkana 2018;Muñoz et al. 2018;Liu et al. 2019), as well as test the existence of any radio background in addition to the Cosmic Microwave Background (CMB) present at these high redshifts (e.g. Feng & Holder 2018;Ewall-Wice et al. 2018;Fialkov & Barkana 2019;Reis et al. 2020).
The first generation of stars, commonly referred to as the population III (Pop III) stars (Bond 1981;Bromm & Larson 2004;Klessen 2019), are expected to have a strong effect on the high-redshift IGM, and, thus, the 21-cm signal from cosmic dawn (Cohen et al. 2016; Hibbard et al. 2022). It has been shown that the efficiency and rate of star formation for the Pop III stars strongly impact the depth of the 21-cm global signal absorption trough (Yajima & Khochfar 2015;Cohen et al. 2018;Mebane et al. 2018;Mirocha & Furlanetto 2019;Schauer et al. 2019;Chatterjee et al. 2020;Muñoz et al. 2022;Reis et al. 2022). Furthermore, recent studies of the transition between metal-free Pop III star formation and the later generation of metal-poor Pop II stars have predicted characteristic signatures of this transition in the 21-cm signal Magg et al. 2022). The sensitivity of the 21-cm signal to the amount of, and the duration of, Pop III star formation comes about due to the several mechanisms through which Pop III stars affect the surrounding IGM.
The Lyman-band emission, between Ly and the Lyman limit, of these first stars gives rise to the WF effect: Lyman line photons scattering off of atomic hydrogen couple the 21-cm spin temperature of hydrogen and the kinetic temperature of the IGM. This effect is dominated by the Ly line (e.g. Pritchard & Furlanetto 2006). Such Ly photons are both produced by redshifting of the continuum photons emitted between Ly and Ly , and as injected photons created in the decay of excited hydrogen atoms. Because Lymanband photons emitted between Lyman lines contribute to the WF effect once they have redshifted into a Lyman line, the emission of the first stars enables the WF effect in a finite region around each star creating Ly coupled bubbles (Reis et al. 2021). Inside these bubbles the spin temperature tends to the IGM temperature rapidly, resulting in stronger absorption in the 21-cm signal within the bubble. Consequently, a deeper global absorption signal and a corresponding peak in the power spectrum are expected to be seen around cosmic dawn as these first Ly coupled bubbles form.
In addition to the WF effect, scattering of the Lyman line photons transfers energy from these photons, and from the CMB, to the IGM, thus raising its kinetic temperature in processes called Ly heating (Madau et al. 1997;Chuzhoy & Shapiro 2006;Reis et al. 2021) and CMB heating (Venumadhav et al. 2018) respectively. Both of these significantly impact the IGM temperature and, thus, the 21-cm signal in scenarios where X-ray heating is weak (e.g. Chuzhoy & Shapiro 2006;Ciardi et al. 2010;Venumadhav et al. 2018;Mittal & Kulkarni 2021;Reis et al. 2021).
Furthermore, higher energy Lyman-band emission of the Pop III stars contributes to the Lyman-Werner (LW) feedback (Haiman et al. 2000). LW photons have energies between 11.2 eV and 13.6 eV and can disassociate molecular hydrogen in the Solomon process (Stecher & Williams 1967;Glover & Brand 2001) Since molecular hydrogen is the primary coolant of gas in the first star-forming halos (Haiman et al. 1996), a strong LW background has the potential to suppress this cooling, thus delaying further star formation and reducing the impact of the Pop III stars on the 21-cm signal. Higher-energy UV photons emitted by the first stars will also ionize the hydrogen gas in a finite region around the stars, eliminating the 21-cm signal from this region (Tanaka & Hasegawa 2021). However, at cosmic dawn, the impact of ionization on the large-scale 21-cm signal is expected to be subdominant. Furthermore, Pop III stars keep affecting the IGM even after the end of their lives via X-rays (Ricotti & Ostriker 2004;Ricotti et al. 2005;Yajima & Khochfar 2015) produced in supernovae explosions (Jana et al. 2019) and by Pop III remnant X-ray binaries (XRBs). The X-ray luminosity of both Pop III (Sartorio et al. in preparation) and Pop II (Fragos et al. 2013) XRBs is anticipated to be significant, potentially providing the dominant source of heating that drives the 21-cm signal into emission (Fialkov et al. 2014;Pacucci et al. 2014).
Due to the multitude of processes through which the Pop III stars affect the IGM and, thus, the 21-cm signal, modelling the highredshift stellar populations is key to making realistic predictions for the 21-cm signal from cosmic dawn. However, much about these stars remains unknown with numerical simulations (Hirano et al. 2014;Greif 2015;Stacy et al. 2016;Hirano & Bromm 2017;Susa 2019;Sugimura et al. 2020;Wollenberg et al. 2020;Sharda et al. 2020Sharda et al. , 2021 and observational constraints (de Bennassuti et al. 2014;Frebel & Norris 2015;Hartwig et al. 2015;Fraser et al. 2017;de Bennassuti et al. 2017;Ishigaki et al. 2018;Magg et al. 2019;Tarumi et al. 2020) having yet to converge on the mass distribution of the first stars, i.e. their initial mass function (IMF). For now, the typical mass scale of the first generation of stars remains uncertain.
In this paper, we aim to investigate how the uncertainties in the Pop III stellar IMF propagate into the uncertainties in the predicted 21-cm signal, taking into account the dominant effects at cosmic dawn, namely the WF effect, Ly heating, CMB heating, and Lyman-Werner feedback. As we are focusing on cosmic dawn, we ignore the impact of X-ray and ionizing photons, which are expected to have a delayed effect on the 21-cm signal (e.g. Cohen et al. 2018). In addition, we limit our investigation to ≥ 12, as at lower redshifts emission from Pop II stars is likely to dominate the Lyman radiation fields and effects of X-ray and UV photons are expected to be significant. Throughout this work, we model the Pop III to Pop II transition in our simulations using the methodology of Magg et al. (2022). Importantly, through this investigation, we relax two assumptions previously made in semi-numerical 21-cm simulations by our and other groups, namely: that (i) the zero-age main sequence (ZAMS) emission rate of a Pop III star is approximately equal to its lifetime mean emission rate, and that (ii) the lifetime of the Pop III stars are short compared to cosmological time scales and so their emission can be modelled as instantaneous. Instead, we calculate the emission spectra of star-forming halos for different IMFs by integrating over individual metal-free stellar spectra, computed from a set of stellar evolution histories and stellar atmospheres. We also implement finite stellar lifetimes in our modelling and take into account the variation of stellar spectra with age.
We begin in Sec. 2 with a brief outline of both the theory of the 21-cm signal and our semi-numerical simulation code. Then in Sec. 3 we discuss the assumptions previously made in 21-cm signal studies and detail our methodology for an improved model of Pop III stellar emission spectra. In Sec. 4, we present our calculation of the Ly and LW radiation fields, which are later used in our simulations to calculate the strength of the WF effect, LW feedback, Ly heating, and CMB heating. In addition, we describe our new developments that extend these calculations to include the effect of arbitrary IMFs and the impact of finite stellar lifetimes. Using our extended radiation field calculation, in Sec. 5 we investigate how important it is to properly model Pop III star emission rates and lifetimes, before assessing the variation of the 21-cm signal with the Pop III IMF. We conclude in Sec. 6.

THE 21-CM SIGNAL
The 21-cm signal is produced by neutral hydrogen gas in the IGM. At the hyperfine level, the = 1 ground state of atomic hydrogen splits into a lower energy singlet state and a higher energy triplet of states. An atom undergoing the forbidden transition between these levels emits or absorbs a photon with a wavelength of 21 cm. Whether a collection of hydrogen atoms will on net emit or absorb 21-cm photons is determined by the relative occupancy of the hyperfine levels typically quantified by the statistical spin temperature S (Scott & Rees 1990): If S exceeds the background radiation temperature (usually assumed to be the CMB temperature 1 ), then atomic hydrogen is a net emitter, otherwise hydrogen will on net absorb photons from the background radiation.
Before reionization, the IGM is primarily composed of neutral atomic hydrogen and so 21-cm absorption/emission occurs throughout the Universe. Furthermore, due to the expansion of the Universe any photons emitted at the 21-cm line are over time redshifted out of the frequency range of the transition preventing re-absorption. Hence, the emission or absorption of 21-cm photons by the IGM at redshift can be seen today at a frequency 21, obs = 21 /(1 + ), where 21 = 1420 MHz is the intrinsic frequency of the 21-cm transition. Commonly, the excess or deficit of photons at frequency 21, obs is quantified in terms of the differential brightness temperature, b , which measures the difference between the observed radiation temperature at 21, obs and the expected temperature of the background radiation. By solving the equations of radiative transfer in an expanding universe, b can be related to S and at the corresponding redshift where 21 is the 21-cm radiation optical depth at that redshift 21 = 3 32 with HI being the fraction of hydrogen in its atomic state, H the number density of hydrogen, 10 = 2.85 × 10 −15 s −1 the spontaneous emission rate of the 21-cm transition, and / the proper velocity gradient along the line of sight including the Hubble flow.
The spin temperature at a given redshift is determined by the balance between the competing influences of the background radiation temperature and gas kinetic temperature, K . The former couples to the spin temperature via scattering of background photons, and the latter via collisional coupling and the WF effect. From the relative coupling strengths of these three processes ( CMB , c , and ) the spin temperature can be determined Venumadhav et al. 2018) where C is the radiation temperature of the Ly photons which mediate the WF effect. Current generation experiments trying to detect the 21-cm signal typically fall into two categories, global signal experiments and power spectrum experiments. Global signal experiments attempt to detect the sky-averaged b , whereas power spectrum experiments aim to detect the spatial variations of the 21-cm signal quantified in the 21-cm power spectrum 21 ( , ) with (k−k ) being the 3D Dirac delta-function and the comoving wavevector. Throughout this paper we shall use the alternative form 1 However, models with excess radio background have been extensively dis- of the power spectrum In this work we utilize our own semi-numerical code to compute the expected 21-cm signal (e.g. Visbal et al. 2012;Fialkov et al. 2014;Cohen et al. 2016;Reis et al. 2021, and other related works). Starting from initial conditions for density, velocity, and temperature fields set over large cosmic volumes of 384 3 comoving Mpc 3 , the simulation evolves radiative backgrounds and computes cubes of b at specified redshifts. From these cubes, the global signal and power spectrum can then be evaluated. The simulation volume is split into 128 3 cubic cells with side lengths of 3 comoving Mpc, containing on average 1.65 × 10 11 M of baryonic mass and 9.02 × 10 11 M of dark matter. Within each pixel halo formation is analytically modelled using the approach of Barkana & Loeb (2004), a hybrid between the Press & Schechter (1974) and Sheth & Tormen (1999) prescriptions. The halo mass function within each cell depends on the local largescale density and velocity fields Visbal et al. 2012). Using a sub-grid prescription for star formation, the pixellevel halo distributions are then converted into star formation rates and ultimately radiation fields, which we describe in more detail in Sec. 3.6 and Sec. 4.
As inputs the 21-cm simulation accepts values of several parameters which parameterize fundamental astrophysical processes driving the evolution of the 21-cm signal (see Cohen et al. 2020;Reis et al. 2021, for recent summary). Unless otherwise stated, in all of our simulation runs presented here we use a minimum virial circular velocity of halos for star formation of c = 4.2 km s −1 (corresponding to the molecular cooling threshold), ionizing efficiency of galactic UV emission of = 15, strong LW feedback (see Fialkov et al. 2013, for details), and a mean free path of ionizing photons of mfp = 50 cMpc (comoving Megaparsec). In addition, we assume no X-ray heating (by setting X-ray emission efficiency = 0); however, all the simulations include both Ly heating and CMB heating (following the implementation of Reis et al. 2021;Fialkov & Barkana 2019, respectively). Furthermore, in our simulations, we model baryon dark matter relative velocities (Tseliakhovich & Hirata 2010;Fialkov et al. 2012;Visbal et al. 2012), assume that the Pop II star formation efficiency in molecular cooling halos (with 4.2 < c < 16.5 km s −1 ) is logarithmically suppressed (Fialkov et al. 2013), model Ly multiple scattering and photoheating (Reis et al. 2021;Cohen et al. 2016), do not take Poisson fluctuations (as introduced in Reis et al. 2022) into account, and assume the CMB is the background radiation. Our assumptions for the stellar IMF are explained in the next section.

MODELLING LYMAN EMISSION OF POP III STARS
In this paper, we focus on the impact of the Pop III star IMF on the 21-cm signal at the onset of star formation. As we are concentrating on cosmic dawn, we only consider effects introduced by Lymanband photons and ignore the impact of ionizing and X-ray radiative backgrounds which we plan to model in future work. The processes triggered by the Lyman-band photons explored here (including WF effect, LW feedback, as well as Ly and CMB heating) depend on the details of the Pop III emission spectra in the Lyman band. Therefore, it is essential to model the impact of the IMF on the stellar spectra. Such a capability was lacking in our 21-cm code and we introduce it here. The calculation is outlined in Fig. 1 and more details are provided in the rest of this section. In summary, we first simulated Model evolution of 100 individual metal-free stars using (Sec. 3.3).
Compute stellar atmospheres and spectra for each metal-free star using (Sec. 3.4).
Combine evolution histories and spectra to find lifetime photon-number emission spectra of each metal-free star ( ; ) (Sec. 3.5).
Store lifetime spectra as lookup-  the stellar evolution history for an ensemble of individual metalfree stars of different masses using the stellar evolution code version 12115 (Paxton et al. 2019, and references therein) following the method outlined in Mirouh et al. (submitted). We, then, computed the emission spectra of those individual stars throughout their lives using the stellar atmosphere code version 205 (Hubeny 1988;Hubeny & Lanz 2017a,b,c;Hubeny et al. 2021). From these spectra and stellar evolution histories we synthesised a dataset of lifetime Lyman-band emission spectra for individual metal-free stars. We then integrate over these individual emission spectra weighted by the assumed IMF to compute the Lyman-band emission spectra per stellar baryon for the desired IMF. The resultant spectra are used in the 21-cm signal simulation to calculate the global signal and the large-scale power spectrum.

Initial mass function
The mass distribution of stars is commonly quantified in terms of the IMF, i.e. the number distribution of stellar masses ( / ) at ZAMS. Since we aim to investigate how 21-cm signal predictions vary with the IMF of the primordial stars, we upgraded our 21-cm code to admit any arbitrary IMF of Pop III stars, and demonstrate the results for four example IMFs each of which is a truncated power-law where is the number of stars with a given mass , and is the exponent of the power-law. As their names suggest, min and max bracket the range of possible stellar masses. The model parameters ( , min , and max ) for each one of the four example IMFs are given in Table 1, along with the characteristic mass scale, 50% , defined so that 50 per cent of stellar mass is in stars of lower masses. Also shown for each IMF is its mean stellar mass¯, mass-weighted mean stellar lifetime¯l ife , and the redshift¯d eath at which a star born at = 20 and living for¯l ife would die. The four IMFs in Table 1 have very different characteristics and are intended to span the large uncertainty in the properties of Pop III stars and the resulting 21-cm signals. The Salpeter (Salpeter 1955) IMF represents an extreme case in which the majority of stars are lowmass, with half of its stellar mass being in stars lighter than 4.1 M . Such stars are very long-lived even by cosmological standards, with the mass-weighted mean stellar lifetime of the Salpeter IMF being 1320 Myr, and the lowest-mass stars it contains having lifetimes of up to 14 billion years (Marigo et al. 2001). In other words, if a population of such stars was formed at = 20 the majority of the stellar mass would be in stars that survive until at least ≈ 4.1. The Salpeter IMF is similar to the one observed for the present-day Pop I stars in the Milky Way.
The Top-Heavy IMF represents the opposite extreme case with half of its stellar mass in stars heavier than 354 M . Such massive stars are short-lived with a typical lifetime of just 2.2 Myr, corresponding to Δ ≈ 0.17 at = 20. We limit our investigation to stars of at most 500 M , as present evolution models suggest that more massive Pop III stars may experience an early onset of photo-disintegration in their core (Heger et al. 2003;Woosley & Heger 2007;Farmer et al. 2019) and, thus, will not contribute significantly to the Lyman output.
The real Pop III IMF is unknown and will likely lie between these two intentionally extreme cases. Therefore, we also explore the Log-Flat and the Intermediate IMF cases, which represent plausible Pop III IMFs suggested by simulations (Greif et al. 2011;Dopcke et al. 2013) and stellar archaeology studies of metal mixing (Tarumi et al. 2020) respectively. These two IMFs are composed of relatively massive stars with 50% = 75 M and 50% = 113 M , respectively, with rather short typical lifetimes of¯l ife = 36 Myr (corresponding to Δ ≈ 2.4 at = 20) and¯l ife = 4.3 Myr (Δ ≈ 0.33 at = 20).
In this work, we model a gradual transition between Pop III and Pop II stars, see Sec. 3.6 for further details. When adding the contribution of Pop II stars to the Lyman-band emission we follow our previous works and assume a Scalo IMF (Scalo 1998) using the stellar spectra for stars of 1/20 solar metallicity computed in Leitherer et al. (1999).

Assumptions in historic modelling of Pop III stars
Many previous 21-cm semi-numerical simulations (Mesinger et al. 2011;Visbal et al. 2012;Fialkov et al. 2014;Reis et al. 2020) have utilized the methods and results of Barkana & Loeb (2005) for modelling the Lyman-band emission of the Pop III stars, which are in turn based upon the stellar spectra calculated in Bromm et al. (2001).
We refer to the combination of these methods and results as our Legacy model of Pop III star emission. Built into this model are two assumptions that we need to reconsider when implementing a more general model for Pop III IMFs: Firstly, it is usually assumed that the ZAMS emission rate of a Pop III star is representative of its average lifetime emission rate (the ZAMS assumption). The second assumption is that the lives of the stars are short compared to the cosmological timescales of interest, and, thus, that the emission of the Pop III stars can be taken as instantaneous (the instantaneous emission assumption).
As we show below, to relax the first assumption we evolve stellar spectra throughout their lives. Our motivation is rooted in an early work by Schaerer (2002) who showed that the ZAMS assumption could introduce errors into the computation of the ionizing fluxes of Pop III stars. This claim suggests that similar errors may be introduced into the computation of Pop III Lyman-band emission. If so, any under/overestimate would in turn propagate through our calculations of the Ly and LW radiation fields into the magnitude of the various Lyman band mediated effects, and, thus, affect the 21-cm signal.
We relax the instantaneous emission assumption in Sec. 4.3. In the derivation of our Legacy model, it was assumed that all Pop III stars were 100 solar masses. Such stars are expected to have lifetimes of ≈ 2.7 Myr (Schaerer 2002) which indeed is a small fraction of the age of the Universe at = 20, 180 Myr (Wright 2006). In this work, however, we are interested in a broad range of Pop III stellar masses (between 0.8 M and 500 M ). As is shown in Table 1, the typical lifetime can be very long for the IMFs dominated by low-mass stars, and so the instantaneous emission assumption is expected to be only valid for IMFs dominated by very massive stars.

Evolution of Pop III stars
To calculate the Lyman-band emission from a population of primordial stars with a given IMF, we first need to calculate the Lyman-band emission spectra for individual metal-free stars with masses in the range between min and max . Since stellar spectrum evolves over the lifetime of a star (Schaerer 2002), we need to trace the evolutionary histories of individual metal-free stars of different masses and then sum up the spectra weighted by IMF.
We use evolutionary histories for 100 metal-free stars with fixed, logarithmically-spaced masses ranging from 0.5 M to 500 M each simulated using the Modules for Experiments in Stellar Astrophysics stellar evolution code version 12115 (Paxton et al. 2019, and references therein). The evolution tracks of individual initially metal-free stars are illustrated in Fig. 2 where we show the effective temperature eff versus the logarithm of surface gravity of the star. Technical details on the internal physics are provided in Mirouh et al. (submitted).
The stars were initialized with Big Bang nucleosynthesis (BBN) proportions of hydrogen and helium and then relaxed onto the ZAMS. The stars were then evolved assuming they were non-rotating, in isolation, and had no mass loss. The absence of metals results in these stars having peculiar nuclear reactions which make them hotter and brighter (for more details, we refer the reader to Marigo et al. 2001;Murphy et al. 2021a). Pop III stars thus exhibit shorter mainsequence lifetimes with respect to their more metallic counterparts.
From ZAMS stars with masses ≤ 310 M were simulated until hydrogen depletion in their core, after which the stars are expected to begin helium burning and expand, becoming red or blue giants. Truncating these stellar evolution tracks at hydrogen depletion is not anticipated to greatly impact the lifetime Lyman-band emission . Evolutionary tracks of individual initially metal-free stars projected into effective temperature eff vs logarithm of surface gravity space. Each line is colour-coded with the stellar mass (see colour-bar on the right). For stars of mass less than 310 M these evolution tracks stretch from ZAMS to central hydrogen depletion, while for more massive stars the evolution tracks end when the star begins to photoevaporate. Also illustrated as black dashed lines are contours indicating 10 per cent and 90 per cent of the time along the evolutionary track. These serve to highlight the evolutionary direction of the stars and the portion of their evolutionary tracks that form the majority of their lives. We thus observe that the evolutionary tracks show that lighter stars initially increase in temperature before cooling toward the end of their lives, whereas for more massive stars temperature monotonically decreases.
calculated for the stars, due to the giant branch representing a comparatively small portion of the life of a star immediately before its death. Our most massive stars ( > 310 M ) were instead evolved until the star began to photoevaporate. The assumption that the emission of massive stars while they are photoevaporating is negligible is not initially well motivated. Instead, this assumption is retroactively supported by our findings in Sec. 3.5 that for such stars the Lymanband emission rapidly decreases before they begin photoevaporating. This is attributed to the effective temperature of these stars sharply dropping as their outer layers expand, for a more detailed discussion see appendix A.

Stellar atmospheres and fluxes
From the stellar evolutionary histories shown in Fig. 2, we need to calculate the spectral flux of the Pop III stars at each point in their lives. Historically, two approaches have been used to address similar problems: either the spectra of the stars were assumed to have a black-body shape with effective temperature eff (Windhorst et al. 2018), or a stellar atmosphere code was used to calculate a detailed stellar spectrum (Bromm et al. 2001;Schaerer 2002). A comparison of these two approaches is shown in Fig. 3, from which we see that the black-body approximation does not match either the shape or the magnitude of the detailed spectra computed with a stellar atmosphere code. We, hence, adopted the latter approach since an accurate magnitude and shape of the Lyman-band emission is required for determining the strength and spatial distribution of the Ly mediated effects of interest. For our calculations, we employed the (Hubeny 1988  chosen due to it previously being successfully employed to model high-mass metal-free stars (Schaerer 2002;Lanz & Hubeny 2003. To aid convergence of the stellar atmospheres we adopted the iterative approach of Auer & Mihalas (1969) calculating first an atmosphere with local thermodynamic equilibrium (LTE), then a non-local thermodynamic equilibrium (NLTE) atmosphere without lines, and finally an NLTE atmosphere with lines, with each step adopting as its initial atmosphere model the result of the previous step. Each atmosphere computation requires three physical parameters: eff , log( ), and chemical abundances. For the latter we adopt the BBN proportions of hydrogen and helium in all of our stellar atmospheres computations, because our simulations predict negligible surface metallicity for all considered stars. This negligible metal mixing to the stellar surface is due to the structure of metal-free stars being primarily composed of radiative layers. All other settings, that determine the numeric routines used in the computations, were fixed to the fiducial values from Lanz & Hubeny (2003) for spectra at eff > 35, 000 K and those of Lanz & Hubeny (2007) for spectra at eff < 35, 000 K.
We computed the stellar flux for atmospheres on a 400 by 240 regular rectangular grid covering eff from 4,000 K to 110,000 K and log( ) from -0.5 to 5.5. These values were chosen to fully cover the stellar evolutionary histories previously shown in Fig. 2. Stellar fluxes at intermediate eff and log( ) were then computed using linear interpolation over the grid of fluxes (in log space). We verified by computing intermediate spectra using that this interpolation results in a root-mean-square error of ≤ 0.5 per cent in the final stellar spectra.
We note that not all the stellar atmospheres in the grid converged, with the failed atmospheres having been either super-Eddington or low-temperature ( eff < 15, 000 K). However, out of the relevant parameter space in Fig. 2, only stars below 2.2 M were affected. As we show below, such low-mass stars have negligible Lyman-band emission despite their long lifetime (shown in Table 1 and Fig. 5).
Therefore, we neglect the Lyman-band emission from stars with masses less than 2.2 M and the lack of convergence at low masses does not affect our results.

Individual stellar spectra
Finally, by combining our stellar evolution tracks and grid of stellar fluxes we calculated the lifetime Lyman-band emission spectra of each individual star of mass Hence the emission rate of a star across all frequencies scales linearly with its area. Similarly it was found increasing eff also increases the total emission rate, but this increase is chromatic with higher eff resulting in harder spectra. Conversely increasing produces softer Lyman-band spectra as well as broadens spectral lines.
Two key results are apparent from our computed individual stellar spectra. Firstly, we find the Lyman-band emission rate of a Pop III star increases over its life due to the area of the star increasing faster than its Lyman-band flux decreases. Thus taking the ZAMS emission as representative would result in an underestimation of the total stellar output. For example, for a Pop III star of mass 10 M , whose spectral emission rate is illustrated in Fig. 4, we find that the emission rate increases by a factor of 3.28 at a frequency of 3 × 10 15 Hz over the stellar lifetime. Second, we find that the total Lyman-band emission per stellar baryon, Ly b ( ), of stars with masses below ∼ 5 M decreases sharply with mass, while stellar lifetimes greatly increase, as is depicted in Fig. 5. Consequently, we confirm that lowmass stars (≤ 3 ) contribute little to the cosmic dawn era Lyman radiation field, as their low emission, caused by their low effective temperatures, is spread out over a long time period (≥ 215 Myr). This finding justifies our assumption made at the end of Sec. 3.4 that low-mass stars ( < 2.2 M ) do not significantly contribute to the Lyman-band emission.
It is worth noting here that we also find the peak in Ly b ( ) occurs for stars of mass 4 -7 M . At higher masses, Ly b decreases as approaches 10 M , then remains roughly constant out to 500 M with the only notable feature being a kink around 300 M , potentially due to the onset of photoevaporation (see discussion in Appendix A). The plateau in Ly b ( ) at 10 -500 M suggests that the total Lyman-band emission from a population of such heavy stars will not strongly depend on the IMF. Based on this finding alone, we anticipate a small difference in the resulting 21-cm signal for the Log-Flat, Intermediate, and Top-Heavy IMFs when only the effects of Lymanband photons are taken into account. However, the differences are expected to grow once the IMF-dependent impact of X-ray heating from Pop III X-ray binaries is taken into account.
We discuss several additional trends in how metal-free stellar Lyman-band emission varies with stellar mass in appendix A.

Star formation rate
Within our semi-numerical simulations we model the spatially varying Pop II and Pop III star formation rate densities (SFRD) via the methodology of Magg et al. (2022). In this star formation prescription each halo forms one generation of Pop III stars when it first crosses the critical virial velocity for star formation, with efficiency * ,III . After this initial burst of star formation, it is assumed haloes take a time recov to recover from the ejection of material from Pop III supernovae. We assume recov = 100 Myr which was shown by Magg et al. (2022) to produce a Pop III to Pop II transition redshift qualitatively consistent with the FiBY simulations (Johnson et al. 2013). Once this time has passed, the halo begins to form Pop II stars with an efficiency * ,II .
Unless otherwise stated throughout this paper we take * ,II = 0.05 and * ,III = 0.002 so that our Pop III SFRDs (shown in Fig. 6) are comparable within a factor of 2.5 to the simulations of Jaacks et al. (2019) at < 20. In this same redshift range the Jaacks et al. (2019) predictions are in general agreement with other numerical studies of Pop III star formation (Yoshida et al. 2004;Greif & Bromm 2006;Tornatore et al. 2007;Wise et al. 2012;Johnson et al. 2013;Pallottini et al. 2014;Xu et al. 2016;Sarmento et al. 2018). At higher redshifts our model predicts more Pop III star formation than that of Jaacks et al. (2019). This is to be expected due to our larger simulation box sizes, side-length of 384 cMpc versus 5.9 cMpc in Jaacks et al. (2019), meaning that our simulations contain more rare overdensity peaks that dominate Pop III star formation at the highest redshifts (Bromm & Larson 2004;Bromm 2013).
An important caveat to our star-formation prescription is that it does not include external metal enrichment of halos. Detailed in an appendix of Magg et al. (2022), the authors found that introducing external enrichment had little impact on the Pop II and Pop III fractions found in their semi-numerical models. Similarly, Visbal et al. (2018) found external metal enrichment to have a very small effect on the global star formation rates in their semi-analytic study of > 20 star . We see that below = 20 the two star formation rates are comparable, differing by at most a factor of 2.5. Due to our larger simulation boxes, which will include more rare overdensity peaks, our model predicts a greater rate of Pop III star formation at higher redshifts. The depicted SFRD is that of the Intermediate IMF model in Fig. 11. formation. However, simulations that fully model the chemical and thermodynamic evolution of the star-forming gas find a significant degree of external metal enrichment (Smith et al. 2015;Chen et al. 2017;Hicks et al. 2021). With externally enriched halos potentially being the site of formation for the first Pop II stars. It is thus probable that Pop III star formation in our simulations is overly clustered.

LYMAN-BAND RADIATION FIELDS
Having calculated the Lyman-band emission of individual Pop III stars we integrated these spectra into our semi-numerical simulations as a lookup table. Utilizing this lookup table we now extend the calculation of the spatially varying Ly and LW radiation fields in our simulation to be IMF dependent.

Ly radiation field
The previous model employed for the Ly radiation field in our simulation code was introduced in Barkana & Loeb (2005). The core idea of this model is that Lyman-band emitted photons can only contribute to the Ly radiation field once redshifted into the Lyman line immediately below their emission frequency. Hence, emission at em only affects the Ly radiation field at abs if it is emitted at a frequency for some , where Ly is the frequency of the th Lyman line. With the further condition that cannot exceed Ly +1 as otherwise those photons would have already cascaded on redshifting into that higher line. As a result, Lyman-band photons emitted above the th Lyman line have a maximum redshift from which they could have an impact on the Ly radiation field at a given absorption redshift of Historically, it was assumed that between their emission and absorption any emitted photons travel in a straight line for a cosmologydependent comoving distance ( abs , ems ). However, this assumption was recently revised in our code by Reis et al. (2021) to incorporate multiple scattering of Ly photons. When absorbed at a Lyman line, a photon cascades into a Ly photon with a recycling efficiency recycle ( ). For Ly continuum photons that redshift directly into the Ly line, recycle (1) is taken as 1. Under these assumptions, plus isotropic emission, it was shown by Barkana & Loeb (2005) that the Lyman radiation at the absorption redshift is where ( , em ) ( abs , ems ),x is the photon emissivity at the emission redshift averaged over a spherical shell of comoving distance about the point the field is being evaluated at x. In our 21-cm simulation code this spherical shell average was replaced by Reis et al. (2021) with a broader window function to take into account multiple scattering of Ly continuum photons.
Under the earlier discussed assumption that star formation is rapid and stellar lifetimes are short compared to the cosmological timescales of interest it was also argued by Barkana & Loeb (2005) that emissivity can be expressed as where¯0 b is the mean baryon density, * the star formation efficiency of the stellar population being modelled, coll / the rate of change of gas collapse fraction into halos at that location, and b ( ) the photon-number emissivity per stellar baryon.
For our purposes of generalising the calculation to an arbitrary IMF  we need to revise the calculation of b ( ). In Barkana & Loeb (2005) and subsequent studies, b ( ) is taken as a piece-wise power-law computed from the stellar spectra provided by Bromm et al. (2001) under the assumption that all Pop III stars were 100 solar masses (Legacy model). In our extended simulations we instead calculate b ( ) from an arbitrary IMF and our pre-calculated individual star lifetime photon-number emissivities ( ; ) as where is the mass of the proton. During implementation, it was found that sampling at 400 logarithmically spaced stellar masses was sufficient to ensure the integral in equation 14 converged to an accuracy greater than 1 per cent, for all IMFs considered in this paper.
The resulting b ( ) for our four example IMFs are shown in Fig. 7. As expected from the flat shape of Ly b ( ) at high masses (Fig. 5), we find the emissivity of the Log-Flat, Intermediate, and Top-Heavy IMFs to be very similar. Hence, we expect these three IMFs to produce similar 21-cm signals. Conversely, the Salpeter IMF Ly b ( ) is markedly different from the other three example IMFs, having lower emissivity at all frequencies, and a softer spectrum. Due to its lower b ( ), we can anticipate that a Salpeter IMF would result in a delayed WF effect compared to the other three example IMFs, and hence a different 21-cm signal.
In addition to the overall intensity, the differences in the spectral shape of b ( ) between IMFs can also affect the 21-cm signal. For example, the lower proportion of high-frequency emission seen in the Salpeter IMF leads to less injected Lyman cooling and thus stronger Lyman heating (Chuzhoy & Shapiro 2006;Reis et al. 2021;Mittal & Kulkarni 2021), as well as reduced LW feedback. Furthermore, the spectral shape of b ( ) impacts the spatial distribution of the Lyradiation field and thus the WF coupling, potentially altering the 21-cm power spectra (Santos et al. 2011). Hence, we anticipate that the 21-cm signal will vary with the Pop III IMF, though it remains to be demonstrated whether these variations are measurable.
For comparison we also plot b ( ) for our Legacy model in Fig. 7. We find that the Lyman-band photon emissivity for our Legacy model is much lower than those for our example IMFs, with the mean b ( ) being 0.26 to 0.34 times that of our example IMFs. By comparing our Legacy case to a model computed using our fully updated framework with the same delta-function Pop III IMF composed purely of 100 M stars, we find that most of this difference is due to the relaxation of the ZAMS assumption and the changes in stellar modelling. This comparison shows the mean b ( ) in our Legacy model is 72 per cent lower than the model using our methodology. Recalculating b ( ) using our stellar evolution histories, but now assuming that the ZAMS emission rate of a star is the same as its average lifetime emission, we find a 54 per cent drop in the emissivity. The remaining 18 per cent is thus due to differences in stellar modelling between our model and that of Bromm et al. (2001). We, hence, see the ZAMS assumption leads to a significant underestimate of b ( ), and, therefore, an underestimation of the magnitude of the Lyman mediated impacts on the IGM. In Sec. 5.1.1 we go on to quantify the resulting impact on the 21-cm signal due to this underestimation and further discuss how it may have biased the 21-cm predictions of our Legacy model.

Lyman-Werner feedback
The LW radiation field, LW (Holzbauer & Furlanetto 2012;Fialkov et al. 2013), is calculated in a similar manner to (equation 12), with a few key differences. These differences are driven by LW being the LW band mean spectral intensity, rather than the photon radiance at the Ly line, as is the case in the calculation. As the Ly line is not in the LW band, cascading photons now leave the band rather than inject into the line of interest. Hence, the sum over cascading lines in equation 12 is replaced with a multiplicative modification factor mod ( abs , em ) inside the integral. mod ( abs , em ) accounts for both the loss of LW photons to cascading but also the relative effectiveness of the different LW lines at triggering H 2 disassociation. More formally, mod ( abs , em ) is defined as the effectiveness of the surviving LW band emission at abs to disassociate molecular hydrogen relative to the emission at the source. In our calculations we use mod from Fialkov et al. (2013).
Because LW is an averaged energetic radiometric quantity, the photon emissivity ( , em ) is replaced in the calculation of the radiative background by the mean LW emissivity,¯( em ). Hence, LW is given by The quantity¯( em ) can be computed in an analogous manner to equation 13 in which b ( ) is replaced by the IMF-dependent LW band averaged emissivity per stellar baryon,¯b. In turn,¯b is computed by taking a photon energy-weighted integral of b ( ) given by equation 14 over the LW frequency range. The¯b calculated for each of our four example IMFs and our Legacy model are given in  Due to the large-scale overdensity and baryon-dark matter relative velocity differing between pixels in our simulations, the star formation rate and emissivity also spatially vary. Once convolved with the relevant window functions and integrated over time, as described by equations (12) and (15), this in turn leads to spatial variations in the computed radiation fields. Inclusion of such variations is essential for accurate calculation of the 21-cm power spectrum with the resulting variations in the strength of the WF effect leading to the cosmic dawn power spectrum peak. Furthermore, the correlation between the strength of the LW radiation field and the matter overdensity in the early Universe is required to properly account for the global suppression of star formation via the LW feedback due to the clustering of the first radiative sources (Ahn et al. 2009). Numerical simulations (e.g. Ahn et al. 2009;Johnson et al. 2013) show variations of the local LW radiation field on sub cMpc scales, an effect which we do not take into account in our work where the LW fluctuations are modelled on the scale of the simulation cell (3 cMpc). These unresolved fluctuations in the LW background will contribute to the spatial variation of star formation rate, and, thus, could boost the 21-cm power spectrum on corresponding scales. However, Johnson et al. (2013) found in their simulations that star formation rate density is mostly regulated by the background LW radiation field due to the large mean free path of LW photons. It is only in rare regions that the LW feedback is driven by the local small-scale star formation. Although the omission of these small-scale variation in the LW field and thus star formation rate is a limitation of our simulation, it is not expected to greatly affect our predictions for the large-scale observable signals targeted by the existing interferometers HERA (Abdurashidova et al. 2022

Non-instantaneous emission
We now relax the assumption of instantaneous stellar emission. For non-instantaneous stellar emission the cosmology, radiative transfer, and atomic physics remain the same and so equations 12 and 15 still hold. The differences in the non-instantaneous emission case instead enter into the calculation of the emissivity because the emission now depends not just on the stars forming at em but those that formed prior and are still emitting at em . Consequently, while equation 13 holds in the case of instantaneous emission it must be generalized to consider the star-formation history and the rate of stellar emission in the case of non-instantaneous emission.
Using our results for the emissivity rate of individual stars (equation 9) we derive a population photon emissivity rate per stellar baryon where ( , ; ) is 0 for larger than the lifetime of a Pop III star of mass . The results are shown in Fig. 8. Considering first the Salpeter IMF we find an initial increase in the emission rate before a plateau stretching from a few Myr to over 300 Myr, and finally a sharp decrease in emission. This plateau shows that a significant proportion of the emission from the low-mass dominated Salpeter IMF is spread out over hundreds of millions of years. Because, in the Planck best-fit cosmological concordance model (Planck Collaboration et al. 2020), redshift = 20 corresponds to 180 Myr after the Big Bang, the emission from stars in the Salpeter case is certainly non-instantaneous. This extended emission is expected to slow down the build-up of the cosmological Ly and LW radiation fields. As a result, the associated impacts on the IGM of these stars are anticipated to be delayed, leading to the characteristic feature of cosmic dawn in the 21-cm signal appearing at later times.
For our other example IMFs emission is concentrated during the first 2 to 3 Myr after star formation, with little emission beyond 10 Myr. The characteristic time-scales of emission for the Log-Flat, Intermediate, and Top-Heavy IMFs are, hence, much shorter than that for the Salpeter IMF. As a result, we expect the instantaneous emission assumption to be a better approximation for these IMFs.
The photon emissivity at a location ì in the simulation box can hence be calculated from b ( , ) by integrating over the previous rate of star formation at that position Here is the time between the onset of star formation and the moment at which the emission rate is sampled, and em is the time corresponding to the emission redshift em . In the limiting case of b ( , ) tending to a delta-function, the instantaneous emission limit is recovered (equation 13). The non-instantaneous stellar emission version of¯( em ) ì used for the LW radiation field calculation can be computed in an analogous manner to ( , em ) ì . Alternativelȳ ( em ) ì can be found by averaging ℎ ( , em ) ì over the LW band. Ultimately we are interested in the 21-cm signal so we quantify the significance of non-instantaneous emission in more detail in the next section using the results of our extended semi-numerical simulations.

PREDICTED 21-CM SIGNAL
In this section, we present our results for our full calculation of the 21-cm signal. We start by assessing the impacts of the highlighted assumptions on the 21-cm signal (Sec. 5.1). We then explore the sensitivity of the 21-cm signal to the IMF of Pop III stars. Firstly, we show the 21-cm signal for a range of models in which all stars are of equal masses (Sec. 5.2). This allows us to determine which stellar masses will create distinct signatures in the observable signal. Second, we compare the 21-cm signal for our four example IMFs (Sec. 5.3) under a range of star formation models. Finally, we briefly address potential parameter degeneracy in Sec. 5.4.

Impacts of the ZAMS assumption
As discussed in Sec. 4.1, the assumption that stellar emission at the ZAMS well represents stellar lifetime emission leads to an erroneously low photon-number emissivity with a 54 per cent reduction in mean b ( ) for a Pop III IMF composed of only 100 stars. The ZAMS assumption thus reduces the WF effect, Ly heating, CMB heating, and LW feedback. To quantify the implications of this assumption on the 21-cm signal we run our semi-numerical 21-cm simulation adopting an IMF composed of only 100 Pop III stars with the ZAMS assumption on and off. The results are shown in the top row of Fig. 9.
The global 21-cm signal predicted for both cases is shown in the top left panel (a) of Fig. 9. As expected the suppressed Lyman-band emission under the ZAMS assumption shifts the onset of the global signal absorption trough to higher frequencies (lower redshift), by Δ ≈ 2. With a resulting maximum relative difference (calculated as the difference divided by the mean of the two quantities) between the fixed redshift global 21-cm signals of 56 per cent, at = 23. The dominant mechanism causing this shift is the decreased WF effect, which is partially compensated for by the weakened LW feedback.
Furthermore, in Fig. 9 we depict the 21-cm power spectra (panel b) and the resulting differences between the power spectra for the two cases (panel c). Similarly to the global signal, the ZAMS assumption leads to a shift of the power spectrum cosmic dawn peak to lower redshifts by Δ ≈ 1. This shift leads to large relative differences between the power spectra computed with and without the ZAMS assumption, the maximal value of which is 103 per cent at = 28 for = 0.022 cMpc −1 . Owing to this impact, we relax the ZAMS assumption in all our subsequent calculations. As our Legacy model employs the ZAMS assumption, we anticipate that previous 21-cm signal predictions utilizing this model have thus been biased. While the exact magnitude of the effect depends on the specific astrophysical parameters, we generically expect a deficit in Lyman radiation which, in the cases explored here, leads to a shift of the high-redshift part of the 21-cm signal to lower redshifts by Δ of a few.

Instantaneous versus non-instantaneous emission
We now explore the implications of the instantaneous stellar emission assumption for the cosmic dawn 21-cm signal showing the results in the bottom row of Fig. 9. We compare the signal calculated with this assumption against the case in which this assumption is relaxed for each one of our four example IMFs. As in the case of non-instantaneous stellar emission a star's radiation output is spread over its lifetime, and stellar Lyman-band emission increases over said lifetime, ab initio we would expect the resulting Lyman-band radiation fields to take longer to reach the same magnitude when compared to instantaneous stellar emission. Consequently, we anticipate a delay in all Lyman-band mediated effects, and so a later global signal absorption trough and cosmic dawn power spectrum peak in the non-instantaneous emission case. Indeed, a mild shift of the signals to lower redshifts is observed in our results ( Fig. 9 panels d-f and Fig. B1).
We find that the signal shift is modest for all our explored IMFs, with the largest impact being on the Salpeter IMF ( Fig. 9 panels d-f) where the proportion of low-mass long-lived stars is greatest. In this case, both the global signal absorption trough and the cosmic dawn peak in the power spectrum shift to lower redshift by Δ ≈ 1 when non-instantaneous stellar emission is modelled. This shift results in a difference between the fixed redshift global 21-cm signal of at most 53 per cent (at = 24), and maximum relative difference between the power spectrum of 119 per cent measured at = 31.
With our other example IMFs, dominated by more massive and thus short-lived stars, the effect is much weaker (Fig. B1) with the maximal relative differences in the global 21-cm signals reaching 11 per cent at = 27, 6.6 per cent at = 28, and 4.1 per cent at = 29 for the Log-Flat, Intermediate, and Top-Heavy IMFs respectively. The Our results suggest that for the global signal instantaneous emission is a reasonable approximation for IMFs dominated by high-mass stars ( 70 M , Table 1), while the IMFs dominated by lower-mass stars (< 5 M ) require non-instantaneous stellar emission to be modelled. However, the effect on the power spectra is much stronger (≥ 22 per cent) for all our considered IMFs. Hence, including noninstantaneous stellar emission in our semi-numerical simulations is necessary if we hope to accurately predict the 21-cm power spectrum during the early phases of cosmic dawn. These results, therefore, motivate us to utilize the fully-modelled non-instantaneous stellar emission in all of our subsequent simulations.
In all our previous papers we used our Legacy model, which is based upon an IMF composed of 100 stars and thus is an example of a high-mass dominated IMF. Therefore we do not anticipate the instantaneous stellar emission assumption to significantly alter our previous conclusions 2 . Hence, while the effects of the ZAMS and instantaneous emission assumptions on the 21-cm signal are in opposite directions, they do not cancel out for our Legacy IMF as the ZAMS assumption impact is much greater. Therefore, previous predictions 2 At high redshifts, an order 10 per cent effect in the power spectra is expected, while the impact on the global signal is smaller. using our Legacy model were biased, under-predicting the Lyman radiation fields in the early Universe and producing lower-redshift cosmic dawn 21-cm signal features due to the ZAMS assumption.

Single-mass Pop III IMFs
We now consider the case of delta-function IMFs, each composed of stars of a single mass . Such simplified IMFs allow us to explore the implications of different stellar masses on the 21-cm signal and so aid in the interpretation of the effects of our more complicated four example IMFs which we show in the next subsection.
The 21-cm global signal and power spectrum at fixed = 0.1 cMpc −1 are shown in Fig. 10 (panels a and b respectively) for logarithmically spaced values of with masses in the range from 2.5 to 500 M . A visual inspection of these results (as well as power spectra at a wide range of scales 0.05 -1.0 cMpc −1 , not shown) suggests that high-mass dominated IMFs produce very similar 21-cm signals with little difference between the models with > 20 M . As decreases, divergence from the generic high-mass behaviour is first seen in the power spectrum, and later in the global signals as well. In both the power spectrum and global signals these deviations manifest themselves primarily at high redshifts > 20. However, the lower , the stronger the impact on the late-time signals.
For IMFs with lower the global signal absorption trough shifts to lower redshifts and the onset of the cosmic dawn power spectrum peak is seen at later times compared to the generic high-mass 21-cm . Testing the ZAMS assumption (top) and instantaneous emission assumption (bottom). We show the global 21-cm signals (left, panels a and d), power spectra at = 0.1 cMpc −1 (middle, panels b and e), and the relative differences in the power spectra (calculated as the difference divided by the mean of the two power spectra, right, panels c and f). Top: The 21-cm signals predicted for an IMF of purely 100 stars, with (dashed) and without (solid) the ZAMS approximation. With the ZAMS assumption, the 21-cm signals are shifted by Δ ≈ 1 to lower redshifts, compared to the model without. These shifts produce relative differences of up 56 per cent (at = 23) in the global 21-cm signals, and relative differences between the power spectra reaching 103 per cent at = 28 for = 0.022 cMpc −1 . Bottom: The 21-cm signals predicted for the Salpeter IMF under the instantaneous (solid) and non-instantaneous stellar emission (dashed) models. Relaxing the instantaneous emission assumption is found to shift the signals to lower redshifts by Δ ≈ 1, resulting in the maximal change of up to 53 per cent at = 24 in the global signal, and a maximal relative difference of 119 per cent between the power spectra (at = 31). A similar comparison for our other IMFs is shown in Fig. B1. signal. We also observe a slight increase in the magnitude of the power spectrum peak. The delay in the signals as we move to lower is attributed to the weaker WF effect, which stems from lower values of the Lyman-band radiation fields. As we have seen above, the dominant mechanism for the latter is the sharp increase in the lifetime of the Pop III stars at lower which leads to the emission being spread over a longer time period. Another mechanism that leads to the reduction in the Lyman-band fields is the inefficiency (over the lifetime) of stars lower in mass than 4 M at production of Lyman-band photons. However, this is a secondary effect, as even at the masses which yield the highest total number of photons per stellar baryon (4 -7 M , Fig. 5) the 21-cm signals are still delayed compared to the cases with less efficient higher mass stars.
By taking the derivative (at a fixed redshift) of both the 21-cm global signal and log power spectrum with respect to log( ), we can identify the threshold mass at which the signals start to deviate from the generic high-mass behaviour. Such derivatives evaluated for the global signal at = 24 and power spectrum at = 30 are shown in panels (c) and (d), respectively, of Fig. 10. Both derivatives fluctuate around zero at the high end, confirming that the signal does not vary with for the high-mass dominated cases. At lower masses, however, large peaks in the magnitudes of both derivatives can be seen, indicating the mass range at which the 21-cm signal becomes sensitive to the mass scale of the Pop III IMF. For the global signal at = 24 we find the derivative to be large (in absolute value) at 10 M , and for the power spectrum (at = 30) the derivative grows at 20 M .
In summary, through the suite of the single-mass test cases, we showed that the 21-cm signal is not sensitive to the shape of the Pop III IMF in the high-mass regime, 20 M . Therefore, if the majority of the Pop III stellar mass was in stars more massive than 20 M as some simulations suggest (Hirano et al. 2014;Hirano & Bromm 2017), the 21-cm signal would not be able to probe the shape of the Pop III IMF in detail. However, we should still be able to test whether or not the first stellar population was composed of predominantly massive stars, by checking if the signal traces the generic high-mass solution. Conversely, if the majority of the Pop III stellar mass is in stars below the 20 M threshold, the 21-cm signal from cosmic dawn can be used to probe the Pop III IMF in more detail. However, degeneracies with other parameters, such as the Pop III star formation efficiency * ,III which we consider in Sec. 5.4, as well as the instrumental noise and systematics, may weaken the constraints or measurements in practice.

21-cm signal with realistic Pop III IMFs
We now consider the variations in the 21-cm signal between our example Pop III IMFs (Table 1), showing the results in Fig. 11. As opposed to our test cases in the previous subsection, each example IMF represents a distribution of masses and, thus, the effect on the 21-cm signal is more complicated. However, owing to the typical masses of stars in the case of the Intermediate, Log-Flat, and Top-Heavy IMFs being very high (with 50% = 75, 113, and 354 M , respectively), we expect the 21-cm signals of these three IMFs to closely follow the generic high-mass behaviour. This is indeed what we find from our simulations (Fig. 11). On the other hand, the Salpeter IMF is dominated by low-mass stars (with 50% = 4.1 M ). Therefore, we expect it to have a delayed 21-cm signal compared to the other cases. This is exactly what we see in Fig. 11, where both the power spectrum and the global signal corresponding to the Salpeter IMF are shifted to lower redshifts by Δ ≈ 1 compared to the more massive-star dominated cases.
Comparing the two IMFs with the greatest difference in cosmic dawn 21-cm signals, namely Salpeter and Top-Heavy, we find power spectrum relative differences of up to 131 per cent at = 30, shown in the panel (c) of Fig. 11. With relative difference > 80% down to = 21. These large differences in the power spectrum predicted between the two IMFs suggest that the 21-cm power spectrum from = 20 -35 is potentially a useful observable for constraining the Pop III IMF if observations are sufficiently precise. Putting this into an experimental context, the redshift = 21 power spectrum differences of 10.2 mK 2 seen between these two IMFs exceeds the thermal-noise sensitivity predicted for 1000 hrs of SKA1-Low observation (Koopmans et al. 2015) by a factor of 2.9.
Furthermore, the shift in redshift of Δ ≈ 1 seen in the 21-cm global signal absorption trough between the two most extreme IMFs leads to a difference of up to 18.1 mK in the global signal at = 19. The EDGES collaboration reported root-mean-square residuals of 25 mK (Bowman et al. 2018) for their existing experiment and the upcoming global signal experiment REACH is forecast to be able to reach a sensitivity of 5 mK using ∼ 2500 hrs of observations (Anstey et al. 2022;de Lera Acedo et al. 2022). Hence, the differences seen between the predicted global signals are comparable to the sensitivity of current global signal experiments and above the sensitivity forecast for near-future experiments. The above suggests that at least in some scenarios the differences between 21-cm signals caused by variations in the Pop III IMF are measurable. However, so far all of the above has assumed one set of astrophysical simulation parameters that lead to the Pop III star formation rate density (SFRD) depicted in Fig. 6. While this SFRD is consistent with some of the existing numerical simulations, we are ignorant of the real star formation history at cosmic dawn due to the lack of observations. Therefore, it is reasonable to consider a series of alternative star formation models. The star formation model considered so far assumed Pop III stars form in molecular cooling halos (critical virial velocity for star formation c = 4.2 km s −1 ), subject to a strong LW feedback (Fialkov et al. 2013), with a star formation efficiency of * ,III = 0.002. Let us now consider alternative star formation histories with a lower star formation efficiency * ,III = 0.0002. We also consider Pop III star formation in atomic cooling halos ( c = 16.5 km s −1 ) as well as the case with negligible LW feedback in molecular cooling halos 3 (potentially due to halo self-shielding). The resulting 21-cm power spectra for the Salpeter and Top-Heavy IMFs at a comoving wavevector of = 0.1 cMpc −1 are shown in Fig. 12. In all of the considered cases we find differences in the power spectrum greater than a factor of 2 above ≈ 20, and so a confident detection of the 21-cm power spectrum in the range = 20 to = 25 would allow the two IMFs to be distinguished. Comparing to the expected sensitivity for 1000 hrs of SKA1-LOW observation (Koopmans et al. 2015;Mertens et al. 2021), we see that the differences are detectable in all the considered cases with * ,III = 0.002. In contrast, for * ,III = 0.0002 the 21-cm signals are considerably weaker and the IMF signature is below the sensitivity of the SKA1-LOW. Additionally, we find that self-shielding of molecular hydrogen in minihalos could result in enhanced differences between the 21-cm power spectrum of the two 3 Recent work by Latif & Khochfar (2019) suggested that the neglection of H 2 self-shielding in Machacek et al. (2001) has led to an overestimate of the strength of the LW feedback (including in our own simulations Fialkov et al. 2013).
IMFs. This is attributable to the enhanced Pop III star formation rate owing to the negligible LW feedback.

* ,III degeneracy
Our simulations predict potentially observable differences between the 21-cm global signals and power spectra of high and low-mass dominated Pop III IMFs. If the 21-cm signal alone is to constrain the Pop III IMF, an important question that needs to be addressed is whether the signature in the 21-cm signal of the IMF is degenerate with other astrophysical processes.
We found earlier that the differences in signal between IMFs arise because of the lower intensity of the cosmic-dawn Lyman-band radiation field, with the Salpeter IMF resulting in a lower intensity compared to the IMFs dominated by more massive stars. Although details of the IMF signature in the 21-cm signal could be unique, a general suppression in Lyman-band radiation could also be caused by a lower Pop III star formation efficiency * ,III or a larger c 4 .
Here, to demonstrate the degree of degeneracy, we concentrate on an illustrative example leaving a full analysis of the parameter space of Pop III star formation models to future work. We take the case of the Top-Heavy IMF with * ,III = 0.002 as a reference. Next, we calculate a series of 21-cm signals for the Salpeter IMF varying * ,III values and keeping all the other model parameters fixed to their reference case values. We then identify the value of * ,III which minimizes the maximum relative difference in the power spectra between the Salpeter case and the reference case. Through this process we find that the minimal maximum relative difference is obtained for * ,III = 0.00528 and is 28 per cent, down from 131 per cent when the models with the different IMFs have the same * ,III . A visual comparison of the 21-cm signal for the Salpeter IMF . We see in all cases the power spectra differ by around a factor of 2 between = 20 and = 25, with differences seen beyond = 30 in some cases. A confident 21-cm power spectrum signal detection in this redshift range would, hence, be able to distinguish between the two IMF models.
with * ,III = 0.00528 and the Top-Heavy IMF with * ,III = 0.002 is shown in Fig. 13.
We find a small but non-negligible difference between the global signals and the power spectra suggesting that the effects of * ,III and the shape of IMF are not fully degenerate. The larger * ,III in the Salpeter case required to match the 21-cm signal for the Top-Heavy IMF is expected because in the case of the Salpeter IMF stellar emission is spread over a longer time period and, thus, is effectively less efficient. With this boosted * ,III , the Salpeter case develops a strong global 21-cm signal faster, being 7.5 mK deeper than the Top-Heavy IMF case at = 20. Variations are also seen between the two power spectra, the maximum relative difference of 28 per cent is found at = 23, with differences of a similar magnitude (17 per cent) at the power spectrum peak at = 17. This 13.3 mK 2 difference seen between the power spectra of the two IMFs at = 17 = 0.1 cMpc −1 is 17.6 times the projected SKA1-Low thermal-noise sensitivity, demonstrating that, in this case, the partial * degeneracy does not eliminate the possibility of constraining the Pop III IMF using the SKA.

CONCLUSIONS
In this work we explore the impact of the IMF of Pop III stars on the cosmic dawn 21-cm signal focusing on the effects mediated by photons in the Lyman band. To bracket the uncertainty in Pop III IMF we consider four different scenarios which include extreme cases as well as IMFs motivated by simulations of star formation. As we focus on the early stages of cosmic dawn, our simulations of the 21-cm signal include WF coupling and LW feedback, as well as several subtle effects such as Ly and CMB heating, multiple scattering of Ly photons, and a gradual transition from Pop III to Pop II star formation. Throughout this work we ignore X-ray heating and ionization which are expected to drive the 21-cm signal at later times.
To rigorously model stellar emissivity in the Lyman band, we simulate stellar evolution histories for a range of initially metal-free stars using the stellar evolution code . We compute the emission spectra of these stars throughout their lives using a stellar atmosphere code . Finally, to compute the Lyman band emission spectra for a population of stars with each considered IMF, we integrate over the individual stellar-emission spectra weighted by the stellarmass distribution. The resultant spectra are used in our 21-cm signal   Figure 13. 21-cm signal star formation parameters degeneracy. We show the 21-cm global signals (left, a), power spectrum at = 0.1 cMpc −1 (middle, b), and power spectrum relative difference (right, c) predicted for the Salpeter IMF with * ,III = 0.00528 and Top-Heavy IMF with * ,III = 0.002. The projected 1000 hour SKA1-Low sensitivity is shown in panel (b) as a grey gashed line. While the positions of the global signal minima coincide, the shapes of the global signals are slightly different, with a broader absorption trough for the Salpeter IMF. These small variations in shape leave residual differences of 7.5 mK between the = 20 global signals. Also the shapes of the power spectra differ with a potentially measurable 13.3 mK 2 difference seen between the two IMFs at = 17 = 0.1 cMpc −1 . Together these differences indicate that the effects of * ,III and the IMF shape are not fully degenerate.
simulations to calculate the global signals and the large-scale power spectra.
Through this work, we also relax two assumptions usually made in the literature, namely that (1) the mean emission rate of a Pop III star is represented by its ZAMS emission rate and (2) Pop III stellar emission can be treated as instantaneous. These assumptions were found to introduce errors of up to 56 per cent and 53 per cent in the global signals and 103 per cent and 119 per cent in the power spectra, respectively. Because such errors are compatible with the expected precision of the upcoming experiments, our results show that accurately modelling the evolution of Pop III stars and their lifetimes is required for precision analysis and interpretation of the cosmic-dawn 21-cm signal. Our models rely on the following assumptions: all stars are assumed to form in isolation, be non-rotating, and have no mass loss on the main sequence. These three assumptions are not fully supported by simulations (Marigo et al. 2003;Stacy et al. 2010;Stacy 2013;Sugimura et al. 2020;Murphy et al. 2021b) and, thus, could be further relaxed in future work.
Using the methods developed in this work, we find that, in the absence of X-ray heating and ionization, the 21-cm signal is sensitive to the Pop III IMF and can probe the typical mass of stars if the IMF is dominated by stars lighter than 20 M . In contrast, the 21-cm signal is not sensitive to the details of the IMF if it is dominated by heavier stars. By considering two extreme examples of a Salpeter IMF dominated by low-mass stars (50 per cent of mass is in stars lighter than 4.1 M ) and a Top-Heavy IMF dominated by stars with > 354 M , we find relative differences of up to 131 per cent in the power spectrum and 59 per cent in the global signal. These factor of 2 differences in the power spectra were shown to be robust against changes in the star formation model parameters used in our simulations. By considering the degeneracy between the effects of the IMF and those of other parameters on the 21-cm signal such as Pop III star formation efficiency, we find that only a partial degeneracy is present in the signal. However, a more quantitative exploration of degeneracy is required and is left for future work. Adding the dependence of X-ray heating and ionization on stellar IMF will increase the discrepancies between different scenarios and help remove the degeneracy between the effects of the IMF and other model parameters.
The signatures introduced by the IMF in the 21-cm signal are found to be at the noise level of existing 21-cm global instruments. Future generations of both interferometers and radiometers are expected to have lower noise, thus allowing us to probe properties of the first stellar population via the 21-cm signal.

APPENDIX A: VARIATION OF LYMAN-BAND EMISSION WITH STELLAR MASS
In the main body of this paper, we present our key findings concerning the Lyman-band emission of individual Pop III stars. Here we continue this discussion in more detail, and show how the emission rate of these stars changes with time and exploring the variation in Lyman-band spectra between stars of different masses.
In Fig. A1 we show the rate of emission of Lyman-band photons for five Pop III stars of different masses. We see that the emission rate sharply rises with stellar age. In the explored scenarios, the increase in the emission rate is driven by the increasing area of  Figure A1. Lyman-band photon emission rate per stellar baryon for five Pop III stars of different masses (as indicated in the legend). As the area of the star increases with stellar age, so does the emission rate, as can be seen in the figure. The increase in emission between zero-age main sequence and hydrogen depletion of the star is larger than a factor of 5 in some cases. The 500 M star shows a peak and a sharp drop in the emission rate in the final period of its life. This sharp drop is believed to be due to the rapid decrease in the effective temperatures of such massive stars as they approach photoevaporation. The low temperature in the stellar atmosphere exponentially suppresses the Lyman flux of the star, resulting in the observed drop in its overall Lyman emission rate.
the stars as they expand when approaching hydrogen depletion or photoevaporation. For some stellar masses the increase is quite large, with the maximum emission rate exceeding 5 times the emission rate at ZAMS (i.e., the emission rate at the beginning of stellar life). The simultaneous decrease in effective temperature (and, hence, reduced Lyman-band flux) seen in these stars as they expand is subdominant to the impacts of the increasing stellar area. Thus, the overall trend of increasing Lyman-band emission is observed. The one exception to the overall observed increase in the emission rate with stellar age is the sharp drop seen for the most massive stars (≥ 310 M ) at the end of their lives. For the 500 M star, illustrated in the figure, the Lyman-band emission rate decreases from its peak to 52 per cent of that value in the final 53,000 years before the star begins to photoevaporate. As the surface area of this massive star is still increasing during this emission drop, the observed drop must be caused by a strong decrease in the Lyman flux of the star. This occurs due to the effective temperature of the star becoming low enough that all layers of the stellar atmosphere that are contributing to its emergent spectra are significantly cooler than the characteristic temperature of Ly photons. For any such layer of the stellar atmosphere the thermal Lyman-band flux is suppressed exponentially with temperature; whereas if the effective temperature is above this threshold, the flux varies linearly with eff . Therefore, as the effective temperature of such a massive star rapidly decreases at the end of its life (as illustrated in Fig. 2) the system enters the exponential regime, and, therefore, the stellar Lyman flux is greatly suppressed. As this suppression is exponential, it dominates the increasing stellar surface area leading to the overall sharp decline in stellar Lyman-band emission. Owing to this rapid decrease in the Lyman-band emissivity of massive stars with age, we can somewhat justify our earlier assumption that truncating stellar evolution tracks at photoevaporation will not greatly impact stellar lifetime Lyman emissivity. An improved model would incorporate this evaporation by considering stars of non-fixed mass and would thus verify whether this truncation is indeed a good approximation. So far our discussions of Lyman-band emission focused mostly on the integrated Lyman-band emission but we also find interesting variations within the Lyman-band spectra. Fig. A2 shows ( ; ) against Pop III stellar mass. The spectral shape of the Lyman-band emissivity is found to vary over the mass range. The highest mass stars ( > 50 M ) have fairly flat emission spectra, whereas the lower-mass stars (< 10 M ) have very little emission at the highest frequencies. Higher energy Lyman-band photons play a key role in Ly heating, as the injected photons from their cascades can cool the IGM (Chuzhoy & Shapiro 2006;Reis et al. 2021;Mittal & Kulkarni 2021) reducing the overall heating effect. Hence, this suggests that the Ly heating may be more efficient for lower-mass Pop III stars due to their low emission of such higher energy photons, this effect is automatically taken into account in our 21-cm signal simulations.
In the spectra shown in Fig. A2, a series of horizontal lines of low emission are apparent. These are absorption lines from the outer layers of the stellar atmospheres. Two separate sets of emission lines are visible, the H Lyman lines and the He Balmer lines. The He Balmer lines are only seen between stars of 5 M to 100 M and are typically fainter than the H Lyman lines. The absence of the He lines in the spectra of the lowest-mass stars is likely due to the low temperatures of these stars meaning they contain negligible He in their outer atmospheres. For the H Lyman lines on the other hand the absorption troughs are most prominent in the lowest-mass stars being very broad and deep, with the lines getting fainter in more massive stars as the rising stellar temperatures ionize atomic hydrogen and so reduces the prevalence of H in the outer atmospheres of the stars.
These variations between the lifetime emission spectra of various masses of Pop III stars are what leads to the IMF Lyman-band photon emissivities per baryon having different spectral shapes, as was previously seen in Fig. 7. The Salpeter IMF has prominent Lyman absorption lines and weak high-frequency emission due to the majority of its emission coming from low-mass stars, whereas the other example IMFs have relatively flat spectra with weak absorption lines as we see is typical of high-mass Pop III stars. Fig. 9 depicted the 21-cm signal computed for the Salpeter IMF with and without the assumption of instantaneous stellar emission. Here we provide an extended version of that figure, Fig. B1, which includes a similar comparison for the Log-Flat, Intermediate, and Top-Heavy IMFs alongside the previously depicted results for the Salpeter IMF. This paper has been typeset from a T E X/L A T E X file prepared by the author. The largest difference between the 21-cm signals computed with and without the instantaneous emission approximation are seen for the Salpeter IMF, with the approximation having a much smaller impact on the 21-cm signals of the other IMFs. For all IMFs the relative differences in the power spectrum (bottom) are most prominent at cosmic dawn around redshift 32 reaching 119 per cent at = 31, 41 per cent at z = 33, 31 per cent at z = 33, and 22 per cent at z = 33, for the Salpeter, Log-Flat, Intermediate, and Top-Heavy IMFs respectively. These modest differences suggest that simulating non-instantaneous stellar emission is necessary for accurate cosmic dawn 21-cm power spectrum predictions.