Signatures of Cosmic Ray Heating in 21-cm Observables

Cosmic rays generated by supernovae carry away a significant portion of the lifetime energy emission of their parent star, making them a plausible mechanism for heating the early universe intergalactic medium (IGM). Following a review of the existing literature on cosmic ray heating, we develop a flexible model of this heating mechanism for use in 3D semi-numerical 21-cm signal simulations and conduct the first investigations of the signatures it imprints on the 21-cm power spectrum and tomographic maps. We find that cosmic ray heating of the IGM is short-ranged, leading to heating clustered around star-forming sites, and a sharp contrast between heated regions of 21-cm emission and unheated regions of absorption. This contrast results in greater small-scale power for cosmic ray heated scenarios compared to what is found for X-ray heating, thus suggesting a way to test the nature of IGM heating with future 21-cm observations. Finally, we find an unexpectedly rich thermal history in models where cosmic rays can only escape efficiently from low-mass halos, such as in scenarios where these energetic particles originate from population III star supernovae remnants. The interplay of heating and the Lyman-Werner feedback in these models can produce a local peak in the IGM kinetic temperature and, for a limited parameter range, a flattened absorption trough in the global 21-cm signal.


INTRODUCTION
Observations of supernova remnants (SNRs) by the Fermi satellite (Ackermann et al. 2013) strongly support the theory that supernovae are a source of cosmic rays.These high-energy charged particles are estimated to carry away 10 to 50 per cent of the kinetic energy of the supernova shock (Drury et al. 1989;Berezinskii et al. 1990;Caprioli & Spitkovsky 2014;Sazonov & Sunyaev 2015) and so represent a substantial portion of the lifetime energy output of massive stars.
Due to the large amount of energy imparted into cosmic rays and the important role played by cosmic rays in determining the thermal and ionization state of the interstellar medium (ISM, Schlickeiser 2002), previous authors (Ginzburg & Ozernoi 1966;Nath & Biermann 1993;Sazonov & Sunyaev 2015;Leite et al. 2017;Jana & Nath 2018;Jana et al. 2019;Bera et al. 2023) have considered them as a potential mechanism for heating and ionizing the neutral intergalactic medium (IGM) at early times.For example, Sazonov & Sunyaev (2015) and Leite et al. (2017) demonstrate that cosmic rays produced by Population III (Pop III) and Population II (Pop II) star supernovae may be able to heat the IGM above the cosmic microwave background (CMB) temperature by redshift 15 and 10 respectively.
★ E-mail: tg400@cam.ac.ukHowever, neither study found cosmic rays contribute more than a few per cent to the ionization of the IGM.
Any efficient heating of the IGM before reionization is of particular interest in the context of the burgeoning field of 21-cm cosmology (Madau et al. 1997;Furlanetto et al. 2006;Pritchard & Loeb 2012;Barkana 2016;Mesinger 2019).21-cm cosmology aims to probe the dark age of the Universe, cosmic dawn, and the epoch of reionization via the degree of emission or absorption at the 21-cm spectral line of neutral hydrogen.The strength of this emission or absorption traces the evolution of the thermal and ionization states of the IGM providing a way to test the astrophysical properties of the first stars and galaxies (Yajima & Khochfar 2015;Cohen et al. 2016;Mirocha et al. 2018;Mebane et al. 2018;Tanaka et al. 2018;Schauer et al. 2019;Mirocha & Furlanetto 2019;Mebane et al. 2020;Tanaka & Hasegawa 2021;Muñoz et al. 2022;Gessey-Jones et al. 2022), and potentially the nature of dark matter (Barkana 2018;Muñoz et al. 2018;Fraser et al. 2018;Fialkov et al. 2018;Liu et al. 2019;Muñoz et al. 2020;Jones et al. 2021;Hibbard et al. 2022;Barkana et al. 2023).Optimistically, cosmic ray heating may therefore provide another mechanism by which early stars and galaxies impact the surrounding IGM, allowing additional information to be extracted about their properties from the 21-cm signal.Pessimistically, the impacts of cosmic ray heating on the 21-cm signal may be degenerate with those of other physical processes of interest thereby weakening possible constraints from a 21-cm signal detection.We leave investi-gations of such degeneracies to a future study, as this paper is focused on the detailed modelling of this heating contribution and providing some illustrative examples.
An additional motivation to explore the potential contribution of cosmic rays to the IGM heating is provided by the disputed global signal detected by the EDGES collaboration (Bowman et al. 2018;Hills et al. 2018;Singh et al. 2022).It was shown by Jana et al. (2019) that cosmic ray heating prevents the radio background from supernovae alone being an explanation for the anomalous depth of the reported signal, while Bera et al. (2023) argued models with cosmic ray heating alongside dark matter-baryon interactions could provide a potential explanation for some features of the best-fit absorption trough.It is, thus, evident that an understanding of the signatures of cosmic ray heating in the 21-cm signal may be necessary for the correct interpretation of the signal.
The aforementioned studies into cosmic ray heating found it to be dominated by sub-relativistic cosmic ray protons.Due to these particles moving at less than the speed of light, and potentially diffusing throughout the IGM by scattering from magnetic irregularities, heating from cosmic rays is anticipated to be more clustered around star-forming sources than alternatives such as X-ray heating (Fialkov et al. 2014b;Pacucci et al. 2014).Leite et al. (2017) and Yokoyama & Ohira (2023) both briefly discussed that this clustering of heating could result in specific features in the 21-cm power spectrum that would not be present in the case of an IGM heated by X-rays, potentially allowing a cosmic ray heated IGM to be distinguished from an X-ray heated one; however, no further study of these features was conducted.
In this paper, we fill this gap in the literature by performing the first 21-cm signal simulations to model the spatial distribution of cosmic ray heating.This is achieved by extending an existing seminumerical 21-cm signal simulation (e.g., Visbal et al. 2012;Fialkov et al. 2014a;Reis et al. 2020) with a parameterized model of cosmic ray heating.By using a semi-numerical simulation, rather than a globally averaged semi-analytic model like previous studies, we are able to simulate full 21-cm tomographic maps and thus calculate the 21-cm power spectrum.Hence, allowing for the first quantitative discussions on the impact of cosmic ray heating on these observables, and on whether or not cosmic ray heating may be distinguishable from other heating mechanisms of the early IGM.Furthermore, as the 21-cm signal is not linearly dependent on the IGM gas temperature, spatially resolved simulations may be necessary to produce accurate predictions of the 21-cm global signal, especially given the anticipated strongly clustered nature of cosmic ray heating.We are thus also able to assess the extent to which the globally uniform heating assumption has biased the conclusions reached by previous studies.Our approach additionally enables us to investigate the sensitivity of cosmic ray heating to other uncertain aspects of the high redshift universe such as the star formation history and the strength of the Lyman-Werner (LW) feedback (Haiman et al. 2000).We hence also consider how the imprints of cosmic ray heating vary with the uncertain astrophysical properties of the first stars and galaxies, and thus verify the robustness of our conclusions.
We begin this paper by reviewing the modelling of cosmic ray heating in previous studies in section 2, which a reader familiar with the literature may want to skip in the interest of time.In section 3 we then recap the theory of the 21-cm signal, outline our semi-numerical simulation code, and describe our parameterized model of cosmic ray heating.Using our extended simulations in section 4 we detail our findings regarding the signatures of cosmic ray heating in various 21cm observables and highlight unique features seen for specific types of cosmic ray heating scenarios.Finally, in section 5 we conclude with a summary of our results and discuss how they could be more generally relevant to the interpretation of the 21-cm signal.

COSMIC RAY HEATING OF THE IGM
A model of high-redshift cosmic ray heating requires several key ingredients: a production mechanism of cosmic rays in the early universe (subsection 2.1); a prescription for how these cosmic rays reach the IGM (subsection 2.2); channels through which cosmic rays lose their energy (subsection 2.3); a way to convert this energy loss into the IGM heat (subsection 2.4); and finally the heating distribution caused by the propagation of cosmic rays through the IGM (subsection 2.5).
Approaches taken to each one of these modelling steps have differed between previous studies (Nath & Biermann 1993;Stacy & Bromm 2007;Sazonov & Sunyaev 2015;Leite et al. 2017;Jana & Nath 2018;Samui et al. 2018;Jana et al. 2019;Bera et al. 2023), due in part to the uncertainties surrounding cosmic ray astrophysics in the early universe.Here we briefly review the prior approaches taken when investigating the impacts of cosmic ray heating on the 21-cm signal.

Production of cosmic rays
As a dying massive star undergoes a supernova it injects 10 51 to 10 53 erg of kinetic energy into ejected material (Woosley & Weaver 1986;Heger & Woosley 2002;Woosley & Heger 2007;Heger & Woosley 2010;Woosley 2010;Whalen et al. 2013a,b;Chen et al. 2014Chen et al. , 2017)).The exact amount of energy is dependent on the type of supernovae and hence the mass of the star, with core-collapse supernovae yielding the lower end of the stated energy range and pair-instability supernovae producing the highest amount of energy.The expelled high-velocity material travels outwards from the star, eventually reaching the boundary of the photoevaporated region surrounding the star.As the high-velocity material meets this denser medium a shock forms.This shock then continues to expand outwards as a supernovae remnant (SNR), gradually losing its energy until it finally dissipates.
A charged particle in the SNR can diffuse back and forth across the shock gaining energy each time, a process referred to as diffusive shock acceleration (DSA, Bell 1978a,b;Schlickeiser 2002;Schure et al. 2012).If a charged particle gains sufficient energy it can then escape the supernovae shock upstream (Ohira & Murase 2019).Lower energy charged particles remain trapped in the shock until the shock slows sufficiently for them to escape or for the shock to dissipate, at which point they are released into the surrounding medium.Except at the highest energies (> 10 7 MeV) cosmic rays released as the shock dissipates dominate the time-integrated cosmic ray spectrum of a SNR (Ohira et al. 2010;Caprioli et al. 2010).This cosmic ray spectrum is theoretically predicted (Bell 1978a;Longair 1994) to follow a power-law in particle kinetic energy  K , / K ∝  −2 K , with an upper-cutoff of ∼ 10 9 MeV (set by the SNR magnetic field strength and radius, Lagage & Cesarsky 1983;Schlickeiser 2002) and a lower-cutoff of ∼ 10 −3 MeV (set by the shock velocity, Sazonov & Sunyaev 2015).Simulations find this spectrum to be insensitive to the environment of the supernovae (Sazonov & Sunyaev 2015).Since cosmic rays with energies above the predicted upper limit of 10 9 MeV are observed, there is strong evidence for a second source of cosmic rays, potentially active galactic nuclei (AGN) or gamma-ray bursts (Alves Batista 2022).However, in this study, we limit our-selves to considering only SNRs as sources of cosmic rays, as AGN are expected to be rare at the high redshifts of interest.
Simulations of DSA and Milky Way observations suggest cosmic rays carry away between 10 and 50 per cent of the initial kinetic energy of the shock (Drury et al. 1989;Berezinskii et al. 1990;Caprioli & Spitkovsky 2014).The majority of this energy is in cosmic ray protons (Schlickeiser 2002;Leite et al. 2017), with per cent level portions of the energy in the form of alpha particles and electrons.Hence, we restrict our heating modelling and this review to considering cosmic ray protons only 1 .

Escape of cosmic rays into the IGM
Based on the findings of numerical simulations (Kitayama & Yoshida 2005;Greif et al. 2007;Whalen et al. 2008), Sazonov & Sunyaev (2015) argued that Pop III SNRs forming in low-mass mini-halos (≲ 10 7 M ⊙ ) would have escaped the virial radius due to their progenitor star having photoevaporated the gas out of the halo.The SNR then only forms a shock front and begins to suffer significant radiative energy losses once it caught up to the material previously ejected from the halo.As a result, the SNR would not dissipate until it has propagated out to the distance of several virial radii of the host minihalo which takes 1 − 10 Myr.Since the recombination time estimated for the gas photoionized by the progenitor star is less than a few million years, Sazonov & Sunyaev (2015) suggested that the cosmic rays produced by the SNR are directly released into the neutral IGM about the source mini-halo.For future ease of reference, we refer to this mechanism for cosmic rays reaching the IGM as direct injection.
In more massive halos, such as the present-day Milky Way, the stars are unable to photoevaporate the gas out to the virial radius.Hence, when supernovae occur in these more massive halos the SNR rapidly encounters the more dense gas causing the SNR to quickly lose its thermal energy radiatively and dissipate, releasing cosmic rays into the ISM.Observations in the Milky Way of the isotopic ratios of Beryllium (Schlickeiser 2002) suggest that cosmic rays diffuse out of the ISM into the surrounding IGM on a timescale of 10 Myr.Leite et al. (2017) argued that since star-forming halos in the early universe are anticipated to be smaller and have weaker magnetic fields than the present-day Milky Way, a large proportion of cosmic rays would be able to diffuse out into the IGM In contrast to the previously discussed direct injection mechanism where cosmic rays are only able to escape from the lowest mass halos, this mechanism which we refer to as diffusive escape does not have an upper halo mass limit, and so cosmic rays can escape into the IGM from all star-forming halos.
The third mechanism by which cosmic rays can escape into the IGM is galactic outflows, studied by Samui et al. (2018).The authors demonstrated that 75 per cent of the IGM could be metal enriched and heated at  = 8 by galactic outflows driven by the pressure of confined cosmic rays if star formation occurs in molecular cooling halos.Due to the growth of these outflows, a significant portion of the cosmic ray energy is lost to adiabatic expansion.As in diffusive escape, there is no halo mass limit from which cosmic rays can be emitted by galactic outflows, though the spectrum of cosmic rays reaching the IGM will differ between the two mechanisms.Due to these similarities between the last two emission mechanisms, we do not distinguish between them in this paper.All conclusions of our 1 Cosmic ray electrons are however important to consider when modelling SNR-produced radio backgrounds due to electrons producing a fraction (  /  )2 stronger synchrotron emission than protons (Jana et al. 2019).
study that apply to cosmic rays diffusing out of halos should be also taken to apply to cosmic ray heating from galactic outflows.

Energy-loss mechanisms of cosmic ray protons
Once in the IGM, cosmic ray protons can lose their energy through numerous mechanisms including excitation, ionization, Coulomb energy exchange with charged particles, and particle production in collisions (Schlickeiser 2002).However, when modelling cosmic ray heating, we only need to consider mechanisms contributing to efficient energy loss by cosmic ray protons in the energy range of interest.Hence, since in the early universe SNRs are believed to be the predominant cosmic ray sources, we limit our considerations to energies below 10 9 MeV, the highest energy SNRs can accelerate cosmic ray protons to.
To determine whether an energy-loss mechanism for cosmic ray protons is efficient we compare it to the energy-loss rate due to the expansion of the universe.As the universe expands the momentum of a cosmic ray proton decreases as the reciprocal of the scale factor.Consequently, the kinetic energy,  K (the energy potentially available for heating), of a proton decreases as where  p is the mass of the proton.In this equation, the first term entails the non-relativistic limit, and the second is the relativistic correction.We can thus define the Hubble cooling timescale to be and for any other energy-loss mechanism (with placeholder label i) we can similarly define its energy-loss timescale as The dominant energy-loss mechanism for a given value of  K at a given redshift  is thus the mechanism with the lowest energy-loss time scale.In the case where Hubble cooling is dominant, the energy of the proton is principally lost to the expansion of the universe and will not contribute to the heating of the IGM.Previous studies have shown point-like electromagnetic interactions (e.g.bremsstrahlung, inverse Compton scattering, and synchrotron) and collisions with CMB photons to be inefficient mechanisms for cosmic ray proton energy loss 2 (Schlickeiser 2002).Furthermore, pion production by collisions between neutral hydrogen and cosmic ray protons is only efficient at energies above the range of interest to us, while resistive heating is only dominant in the highly ionized regions around galaxies invisible to the 21-cm signal (Yokoyama & Ohira 2022, 2023).Alongside Hubble cooling this leaves three energy-loss mechanisms to consider3 : excitation and ionization of neutral hydrogen atoms (Sazonov & Sunyaev 2015;Jana & Nath 2018), Coulomb interactions with free electrons (Leite et al. 2017), and Alfvén wave emission (Samui et al. 2018;Bera et al. 2023).
If a cosmic ray proton collides with a neutral hydrogen atom it can directly ionize the atom.In each such interaction the cosmic ray proton losses ≈ 60 eV of energy to the liberated electron (Spitzer & Scott 1969;Sazonov & Sunyaev 2015), which can then ionize or excite surrounding hydrogen atoms in secondary collisions.The energyloss rate of cosmic ray protons to this process is described by the Bethe-Bloch equation, which for protons with  K < 8.4 × 10 5 MeV is well-approximated by (Schlickeiser 2002) In the above equation, Θ is the Heaviside step function,  HI the hydrogen neutral fraction,  H the physical number density of hydrogen,  = v/ the normalized velocity of the proton, and  0 = 0.01 the loss rate peak.Due to the small value of  0 , this process is anticipated to be the strongest for non-relativistic protons, and the strength of this mechanism will decrease as the IGM ionizes owing to its scaling with  HI .
Cosmic ray protons can also transfer energy directly to any free electrons in the IGM via Coulomb interactions, with this energy again being dissipated into heat in the IGM via the subsequent collisions of these energetic electrons.The equations to describe the rate at which cosmic ray protons lose energy to this mechanism were derived by Butler & Buckingham (1962), Sivukhin (1965), andGould (1972) and are well approximated this time by (Schlickeiser 2002) where   is the electron fraction,   the number density of baryons, and  e is the temperature of the free electrons which we take to be the kinetic temperature of the IGM,  K .Hence, unlike excitation and ionization interactions, the efficiency of this mechanism increases as the IGM becomes more and more ionized.Lastly, since cosmic rays are charged they gyrate along magnetic field lines B generating Alfvén waves (Samui et al. 2018) which carry away a portion of their energy.The dissipation of these Alfvén waves then transfers this energy to the IGM as heat at rate |B • ∇ cr |/(4) where  cr is the pressure from cosmic rays and  is the IGM physical density.Due to the dependence on  cr this is a bulk effect, unlike the energy-loss mechanisms mentioned previously, which can be considered in terms of individual protons.Consequently, we cannot define an energy-loss timescale for this mechanism in the same way we have done above.However, Bera et al. (2023) provides an estimate for when Alfvén wave dissipation is an efficient heating process by comparing the IGM heating timescale in this case ( A =  K /  K ) with the Hubble time ( H = 1/), note the latter differs by a factor of order unity from  H defined in equation (1).For the Planck 2018 best-fit ΛCDM cosmology (Planck Collaboration et al. 2020) this ratio becomes Here  0 is the magnitude of the primordial IGM magnetic field seen today,  cr the cosmic ray energy density, and  the typical length-scale of cosmic ray pressure gradients taken by Bera et al. (2023) to be the inter-halo distance.There is much uncertainty surrounding  0 with upper bounds of ∼ 1 nG from CMB observations at 1 cMpc scale (Planck Collaboration et al. 2016) and lower limits from second-order perturbation theory of ∼ 10 −11 nG at 10 ckpc (Ichiki et al. 2006).Recently there have also been attempts to measure the magnetic field in voids today with a disputed lower bound due to Blazar observations of ∼ 10 −6 nG (Tavecchio et al. 2011), and potential measurement of ∼ 10 −7 nG from gamma-ray bursts (Xia et al. 2022), though it is still unclear if the present-day void magnetic field is primordial in origin or has been enhanced by outflows from galaxies (Samui et al. 2018).Thus, the reference value of  0 adopted above is at the upper limits of the 11 orders of magnitude permitted range, with the majority of this range ( 0 ≲ 0.01 nG) leading to inefficient Alfvén wave heating of the IGM for the redshift range of interest.Therefore, we elect to ignore Alfvén wave losses from cosmic rays in our modelling and instead revisit this mechanism briefly in section 5 to discuss how our results would be changed by its inclusion due to a high  0 .Note as the magnetic fields essential to DSA (see Section 2.1) are produced in the SNR shock itself by the Weibel instability, a weak background magnetic field would not prevent cosmic ray acceleration by the first supernovae (Ohira & Murase 2019).
The relative importance of energy-loss mechanisms for cosmic ray protons can now be determined via a comparison of their energy-loss timescales   , with the lowest value corresponding to the dominant process.Such a comparison is shown in Fig. 1 for  = 8, 12, 20, and 30.In the figure, we normalize all energy-loss timescales via  H to allow for easier comparison between redshifts.Since reionization is anticipated to be patchy (Choudhury 2022),  e (and with it efficiencies of some of the aforementioned mechanisms) is expected to vary greatly between different regions in our simulations.To demonstrate the importance of this effect for the cosmic ray heating, we show Coulomb heating timescales for both a fully ionized region (the case referred to as Coulomb Ionized Bubbles in the figure) and for the average value of  e outside of ionized bubbles taken from the reference simulation presented in subsection 4.1.1(the case referred to as Coulomb Outside Bubbles).Note that inside the ionized bubbles no excitation or ionization interactions occur, so the timescale depicted is that for cosmic rays outside of ionized bubbles.As we see from the figure, Hubble cooling dominates above  K ∼ 30 MeV at all redshifts, suggesting higher energy relativistic cosmic rays do not contribute substantially to the IGM heating in agreement with the findings of Sazonov & Sunyaev (2015).At lower energies, outside ionized bubbles, excitation and ionization energy losses are dominant at all redshifts, while inside ionized bubbles Coulomb losses are dominant.Our findings at lower energies and low redshifts ( ∼ 8) seemingly contradict the conclusions of Leite et al. (2017), who showed Coulomb energy-losses dominate excitation and ionization losses at lower energies towards the end of reionization.However, this is not a discrepancy as Leite et al. (2017) used a globally averaged  e while we have distinguished between fully-ionized and non-fully-ionized regions of the IGM.We reach the same conclusion as Leite et al. (2017) when the globally averaged  e is used.
Since this paper is principally concerned with modelling the 21cm signal of neutral hydrogen, we are primarily interested in heating outside of fully ionized bubbles.Therefore, the mechanisms of interest to us are excitation and ionization losses, and Hubble cooling.However, due to the similarity between the energy-loss timescales of Coulomb interactions inside ionized bubbles and excitation and ionization interactions outside of them, the resulting energy evolution  or electron-positron pairs (brown), and radiative losses (purple).For Coulomb energy losses we distinguish between the IGM in an ionized bubble (dot-dashed) and outside of ionized bubbles (solid).Inside ionized bubbles we take  e = 1, and outside we use values of 0.0072, 0.0014, 0.0002, and 0.0002 in ascending redshift order from the output of the reference simulation (subsection 4.1.1). H is taken as the cosmological average value.At high energies,  K ≳ 10 MeV, we find that Hubble cooling is the dominant mechanism (the lowest energy-loss timescale) both inside and outside ionized bubbles.For lower energies, Coulomb dominates in ionized bubbles, and excitation and ionization dominates outside of ionized bubbles.Pion-producing collisions with matter, photon collisions, and radiative losses remain subdominant for the energy range considered here, with the latter three being too inefficient to appear on the scale shown.
equation for cosmic ray protons is anticipated to be a good model for cosmic rays propagating through ionized regions as well.

Cosmic ray heating and ionization
Not all of the energy lost by cosmic rays, discussed above, is converted to the IGM heat.None of the energy lost to Hubble cooling is transferred to the IGM, and of the energy lost to excitation and ionization only a portion  heat ultimately ends up as heat energy in the IGM with the rest lost to the ionization energy and the de-excitation spectral line emission of excited hydrogen atoms.Spitzer & Scott (1969) calculated the heat transfer in excitation and ionization interactions of cosmic ray protons for several different free electron fractions  e .We note the value of  heat is quite sensitive to small changes in  e , varying from 10% to 53% between  e values of 10 −4 and 0.5.Therefore, while recent studies found cosmic rays do not contribute significantly to reionization (Sazonov & Sunyaev 2015;Leite et al. 2017), ionization, including that from cosmic rays, may strongly impact the cosmic ray heating rate.Thus, it should be considered and modelled self-consistently.
Per primary ionization of a neutral hydrogen atom, a cosmic ray proton losses ≈ 60 eV (Sazonov & Sunyaev 2015), of which the fraction  heat contributes to heating the IGM.The primary ionization rate per baryon Λ primary cr can then be found from the cosmic ray heating rate per baryon,  cr , Primary electrons freed by cosmic ray interactions can then cause secondary ionizations in their subsequent collisions.The average number of secondary ionizations per primary ionization  for various values of  e were also calculated by Spitzer & Scott (1969).Thus the full cosmic ray ionization rate per baryon is given by

Propagation of cosmic rays through the IGM
Since the early universe is not homogeneous alongside the rate of cosmic ray heating in the IGM, we must also consider its spatial distribution.Most previous studies of cosmic ray heating (Sazonov & Sunyaev 2015;Leite et al. 2017;Jana et al. 2019;Bera et al. 2023) have assumed the heating to be spatially uniform.This assumption is partially motivated by arguments that cosmic ray heating should reach inter-halo distances (≳ 0.1cMpc at  = 20), based upon the Bohm diffusivity lower bound and Larmor radius of cosmic rays (Stacy & Bromm 2007;Sazonov & Sunyaev 2015).
A more detailed study of the distribution of cosmic ray proton and the heating pattern around the source halos was undertaken by Jana & Nath (2018), who modelled cosmic rays as diffusing out of their source halos with a constant diffusivity .This is physically motivated by cosmic rays scattering from magnetic fields leading to them undergoing a random walk and is similar to models of cosmic ray propagation through the Milky Way ISM such as GALPROP (Strong & Moskalenko 1998;Ptuskin 2012).From their study Jana & Nath (2018) conclude that cosmic ray heating around minihalos between  = 10 and 20 was uniform, while it would be inhomogeneous around high-mass galaxies with low star formation rates.However, importantly, they conclude their results are not definitive due to the large uncertainties in the value of the cosmic ray diffusivity in the early universe.These uncertainties come about in part due to the large, up to 11 orders of magnitude, uncertainty in the IGM magnetic field discussed in the previous subsection, and due to the scaling of cosmic ray diffusivity with cosmic ray energy and magnetic field being unknown.Even if cosmic ray heating is uniform on the typical distance scale between minihalos as Jana & Nath (2018) tentatively concluded, this does not necessarily mean cosmic ray heating is uniform throughout the universe since at early times star-forming minihalos are highly clustered around rare overdensities (Barkana & Loeb 2004;Barkana 2016).Instead, it would imply cosmic ray heating is uniform on the scales on which minihalo number density is uniform.
The uncertainty surrounding the propagation of cosmic rays in the early universe IGM leads us to consider multiple different modes of cosmic ray propagation in our study, both to bracket the range of physically plausible scenarios and also to allow us to investigate the potential of the 21-cm signal to distinguish between these possibilities.

SEMI-NUMERICAL 21-CM SIGNAL SIMULATIONS
Having outlined the five necessary steps to modelling cosmic ray heating in the previous section, we introduce our 21-cm signal simulation code and detail the changes we make to include the new heating mechanism.

Background theory and current simulation code
21-cm cosmology aims to observe the absorption or emission at the 21-cm spectral line of atomic hydrogen (Furlanetto et al. 2006).The strength of the 21-cm signal is dependent on the number density of atoms in the hyperfine states  1 , where the magnetic moments of proton and electron are aligned, compared to the state  0 , where the magnetic moments are anti-aligned.Additionally, the signal depends on the background flux of photons at the 21-cm line, commonly quantified in terms of the radiation temperature   .We can conveniently express the relative occupancy of the hyperfine states in terms of spin temperature,  s , via the use of the Boltzmann distribution (Scott & Rees 1990) where the factor of three comes from the three-fold degeneracy of the higher-energy states, and  21 = 5.87 eV (corresponding to  21 = 1420 MHz) is the small energy difference between the aligned and anti-aligned states.The relative value of  s to   determines whether the signal is to be seen in emission ( s >   ) or absorption ( s <   ).A more detailed treatment of the radiative transfer problem, accounting for the number density of hydrogen  H , hydrogen neutral fraction  HI , and the gradual redshifting of photons into and out of the spectral line allows for the calculation of the differential 21-cm brightness temperature that would be seen today formally defined as the difference in radiation temperature seen at  =  21 /(1 + ) due to the emission/absorption at the 21-cm line by the IGM at redshift .In the above equation  21 is the 21-cm optical depth given by where  10 = 2.85×10 −15 s −1 is the spontaneous emission rate of the 21-cm transition, and v ∥ / ∥ the proper velocity gradient along the line of sight including the Hubble flow.In summary, simulating the 21-cm signal requires determining   ,  s , and  HI .In this study, we assume the background radiation temperature   to be just the CMB radiation, 2.725(1 + ) K.We leave the simultaneous consideration of excess radio backgrounds (Feng & Holder 2018;Ewall-Wice et al. 2018;Fialkov & Barkana 2019;Reis et al. 2020) and cosmic ray heating to future investigations.The calculation of the spin temperature  s of the IGM is complicated due to it being acted upon by three competing influences: collisional coupling to the kinetic temperature of the IGM ( K ), radiative coupling to the background 21-cm radiation temperature (  ), and indirect coupling to  K via the scattering of Lyman line photons, the Wouthuysen-Field effect (WF, Wouthuysen 1952;Field 1958).The efficiency of each of these processes at any given redshift is encapsulated by the corresponding coupling coefficient,  c ,   , and   , which are in turn derived from atomic physics and the intensity of the Lyman line radiation fields (Furlanetto et al. 2006;Venumadhav et al. 2018).Once these coupling coefficients are found the spin temperature can be computed iteratively using where  c is the colour temperature of Ly  photons (Barkana 2016).
At cosmic dawn, the first generation of stars is anticipated to have filled the early universe with a background of Ly  photons causing the WF coupling to dominate the other processes (with   ≫   ,  c ) resulting in  s ≈  K .Hence before reionization when  HI ≈ 1, but after cosmic dawn we anticipate (from equation 10) that the 21-cm signal should trace variation in the gas temperature of the IGM, making it particularly sensitive to cosmic ray heating along with other heating/cooling mechanisms.However, since our study spans a wider range of redshifts than those for which the aforementioned approximation is valid we always utilize equation 12 when calculating  s .
We thus also need to determine  K and the Lyman line radiation fields to compute  s . K can be calculated by solving the ordinary differential equation describing the balance of heating and cooling influences on the IGM (Mesinger et al. 2011) The first term encompasses the heating/cooling rate per baryon,  p , due to various astrophysical processes, in our simulations, we include X-ray heating (Fialkov et al. 2014b), Compton heating (Madau et al. 1997) we extend this list to include cosmic ray heating as well.Next, the second term represents heating due to structure formation, and the third term takes into account cooling due to the increasing number of particles caused by the reionization of the IGM.Lastly, the fourth term represents the adiabatic cooling of the gas due to the expansion of the universe.
To compute the Lyman line radiation fields, and thus   , we follow the methodology of Reis et al. (2021), an extension of Barkana & Loeb (2005) and Fialkov et al. (2014a).Hence, we also assume the Wouthuysen-Field effect is dominated by Ly  photons due to such photons scattering orders of magnitude more times than higher Lyman line photons (Furlanetto et al. 2006).As a result, in this model, the WF coupling is taken to simply be directly proportional to the Ly  line radiation field   (Madau et al. 1997) where   = 0.4162 is the oscillator strength of the Ly  transition.Two separate contributions to   are modelled, photons that directly redshift into Ly  and photons produced in the cascading decays of neutral hydrogen excited by higher Lyman line photons.For the latter, it is assumed that photons emitted between Ly  and the Lyman limit travel on radial geodesics away from sources until they redshift into the Lyman line below their emission frequency and cascade into a Ly  photon with a line dependant recycling fraction.Whereas the paths of photons emitted between Ly  and Ly  that directly redshift into Ly , are treated as random walks due to their many scatterings in the tails of the Ly  line.In the simulation, both radiative transfer processes are implemented using isotropic window functions.These window functions are analytic spherical shells for the cascading photons and fits to Monto Carlo simulations for directly redshifting photons.By convolving the window functions with emissivity fields, which are in turn derived from the past star formation rate and stellar population,   and   can be calculated efficiently throughout the simulation box.Finally, we need to compute the neutral fraction  HI of the IGM, for which we adopt an approach similar to Mesinger et al. (2011).We use an excursion set-based formalism (Furlanetto et al. 2004) to identify regions that have been fully ionized ( HI = 0) by galactic UV emission, which is anticipated to be the dominant mechanism driving reionization (Yajima et al. 2011(Yajima et al. , 2014;;Wise et al. 2014;Ma et al. 2020).A region is considered fully ionized if there exists a spherical volume of radius  centred on that region in which the time-integrated ionizing UV emission exceeds the effective number of neutral atoms within said volume.Hence, for an effective galactic ionization efficiency per baryon , and collapse fraction of baryons into galaxies  coll (x, ) averaged over a volume of radius  centered at x, the point x is considered ionized if Where  mfp is some maximum radius UV photons can travel to, normally set to the mean free path of ionizing photons in the ionized IGM at the end of reionization (Furlanetto & Oh 2005), hence the choice of notation.However, such an approach on its own would ignore the residual ionization leftover from recombination and any ionization from other sources such as X-rays (Pacucci et al. 2014) or the cosmic rays we describe in this paper.Hence, we use a modified version of equation ( 15) introduced by Mesinger et al. (2013) to take into account the region already being partially ionized by these other mechanisms.With x considered ionized if where  e,oth (x, ) is the ionization fraction of the IGM when not including the UV contribution, averaged over the same sphere as  coll (x, ). e,oth (x) in turn is determined by solving the ionization differential equation with the first term encapsulating ionization at a rate of Λ ion per baryon from non-UV sources (e.g.X-rays and cosmic rays), and the second term modelling type-B recombination (Furlanetto et al. 2006) with recombination coefficient  B = 2.6×10 −13 ( K /104 K) −0.7 cm 3 s −1 .It is assumed in our simulations that the ionization fraction of hydrogen and helium is the same, as a result helium double ionization is not modelled.This reionization model is hence built on the assumption that ionization from UV is isolated to fully ionized bubbles, whereas other sources of ionization can travel into the neutral IGM and give rise to partial ionization.In practice, such bubbles may be smaller than the resolution of a simulation, which if not accounted for could lead to erroneous 21-cm signal predictions.To mitigate this issue regions that are not fully ionized are modelled as a two-phased medium, one fully ionized and one ionized to  e,oth (x), an approximation that was validated by Zahn et al. (2011) against radiative transfer simulations.The relative proportion of the region taken up by each phase being given by the region's overall neutral fraction Combined this multi-stage approach to ionization allows for the modelling of the propagation of fully ionized bubbles on both subresolution and resolved scales, while also accounting for the evolution of the ionization fraction of the mostly-neutral IGM outside these bubbles.
The above equations combined with a parameterized prescription of star formation rate and stellar/galactic emission form the foundations of our semi-numerical 21-cm signal simulation code (e.g., Visbal et al. 2012;Fialkov et al. 2014a;Cohen et al. 2016;Reis et al. 2021).The simulation starts from cosmological initial conditions for overdensity and baryon-dark matter relative velocity computed using CAMB (Lewis et al. 2000;Lewis & Bridle 2002;Lewis & Challinor 2011), and initial conditions for the gas temperature and residual ionization fractions computed using RECFAST (Seager et al. 2011).
For the purposes of this study, all these fields are created on a 128 3 grid of cells, each cell being a 3 cMpc sided cube.These initial conditions are then evolved forward in time, at each step an analytic prescription being used to compute the expected halo mass distribution within each simulation cell (Barkana & Loeb 2004;Fialkov et al. 2012).The code then employs the star formation prescription of Magg et al. (2022b) to translate these halo mass distributions into Pop III and Pop II star formation rates for each cell and in turn into an emission rate for various types of electromagnetic radiation.Our radiative transfer model then propagates the radio, Lyman band, and X-ray radiation fields via pre-computed window functions 4 .With the radiation fields in each cell, equations ( 10), ( 12), ( 13), ( 14), ( 16), ( 17) and ( 18) are then solved simultaneously to determine the expected Due to the uncertainty surrounding the universe between recombination and reionization, many aspects of our semi-numerical simulations are parameterized as exact values are not known.Except when stated otherwise in this paper we use the parameters listed in Table 1.
Our ionization model is described by two of these parameters  the effective galactic ionization efficiency per baryon, and  mfp the mean free path of ionizing photons in the ionized IGM at the end of reionization.While  mfp has a direct physical interpretation  is phenomenological, encapsulating the star formation efficiency of ionizing sources, the ionizing UV escape fraction, the stellar UV emissivity, and the reduction in effective ionization efficiency due to recombinations and clumping of the IGM (Furlanetto et al. 2004).In an upcoming paper, we refine this reionization model by modelling the individual physically-interpretable components of .Since we observe the impacts of cosmic ray heating on the 21-cm signal to be strongest prior to  = 10 we find varying these reionization parameters has no impact on our qualitative conclusions.Hence throughout this paper, we fix the values of these two parameters.With  set to 15, to recover optical depths to the CMB of  ≈ 0.06 consistent with the Planck 2018 measurements (Planck Collaboration et al. 2020), and  mfp fixed to 50 cMpc, motivated by the theoretical expectation of a mean free path of ionizing photons in the ionized IGM of 70 cMpc (Wyithe & Loeb 2004) at  = 6.
For all of our simulations, we enable baryon-dark matter relative velocities (Fialkov et al. 2012), Ly  heating (Reis et al. 2021), CMB heating (Venumadhav et al. 2018), Pop II star formation rate suppression in low-mass halos (Fialkov et al. 2013), Ly  multiple scattering (Reis et al. 2021), photoheating feedback (Cohen et al. 2016), and model the non-instantaneous emission of Pop III stars (Gessey- Jones et al. 2022).In this work, we do not enable Poisson fluctuations of star-forming halos (Reis et al. 2022).We derive the Pop III star Lyman band and Lyman-Werner band emissivities self consistently from the power-law IMF specified in Table 1.The ionizing UV emissivity and X-ray emissivity of Pop III stars are not currently derived from the IMF but will be in future studies (Gessey-Jones et al. in preparation; Liu et al. in preparation).
Unlike future experiments such as SKA2-LOW (Koopmans et al. 2015), which are expected to make full tomographic maps of  21 , experimental efforts so far have concentrated on measuring either the sky-average of the signal ⟨ 21 ()⟩ (Voytek et al. 2014;Sokolowski et al. 2015;Bowman et al. 2018;Price et al. 2018;Philip et Price et al. 2018;Mertens et al. 2020Mertens et al. , 2021;;Abdurashidova et al. 2022).Throughout this paper, we use the Δ 2 power spectrum convention, where k is the comoving wavevector and   the Dirac delta-function.
Given the experimental interest, we concentrate our study on cosmic ray heating signatures in the global 21-cm signal, power spectrum and tomographic maps.For our analysis, we include redshift space distortions when computing the 21-cm global signal and power spectrum but not when depicting 21-cm tomographic maps.

Updating Lyman-Werner feedback prescription
In addition to introducing cosmic ray heating as described in the subsequent subsection, we make a further modification to our seminumerical simulation code.Previous versions of our simulation code utilized the fitting formula from Fialkov et al. (2012) where  21 is the LW band intensity normalized to 10 −21 erg s −1 cm −2 sr 1 Hz −1 , and where v rms , is the root mean square velocity of the baryon-dark matter relative velocities at the relevant redshift.The critical halo mass for molecular cooling star formation is then given by where the redshift dependant term is motivated by the simulation findings of Kulkarni et al. (2021) and theoretical predictions of Tegmark et al. (1997).In this study we take  =20 = 5.8 × 10 5 M ⊙ .This value was chosen due to it being 1 below the mean mass at which molecular cooling halos start forming stars in the  21 = 0, v bc = 0 simulation of Schauer et al. (2021).We use 1 below the mean (instead of using the minimum or average mass) to account for the fact that we have a sharp cutoff with mass in halos that are Pop III star-forming (as opposed to the gradual transition found in the simulation).Therefore, using the minimum or average would lead us to over-predict or under-predict the star formation rate respectively.Additionally, in our simulations, we correct for the delayed response between molecular cooling star formation and the exponentially growing  21 by employing the methodology of Fialkov et al. (2013).Instead of using the  21 field of the current simulation time step in the feedback calculation, we use the field from when the universe was a fraction  LW = 0.75 of its age at that timestep. LW can be varied from this fiducial value to simulate stronger and weaker than expected LW feedback.The above feedback prescription allows us to model the minimum halo mass required for star formation due to molecular cooling.However, star formation can also occur via atomic cooling or, if a previous generation of stars has metal-enriched a halo, metal cooling.To account for atomic cooling we take the minimum halo mass for Pop III stars to form as the lower of the molecular cooling threshold  mol,crit and the atomic cooling threshold  atm,crit (Fialkov et al. 2013), which is not impacted by the LW feedback.Since we follow the star formation prescription of Magg et al. (2022b) we model metal-cooling Pop II star formation separately to Pop III stars.This prescription is based upon fits to A-SLOTH (Magg et al. 2022a) merger tree models, and assumes Pop II star formation occurs rapidly due to metal-cooling once a halo has had time to recover from the ejection of material by the supernova of the first generation of stars.External metal enrichment of halos is not included in the model as the authors of that study found including such a mechanism had little impact on the Pop II and Pop III star formation rates they derived.Since Pop III and Pop II star formation are triggered by distinct cooling mechanisms our simulations take separate star formation efficiencies  * ,III and  * ,II for these processes.

Our cosmic ray heating model
We now describe the model for cosmic ray heating we have integrated into our semi-numerical simulations.The model was intentionally developed to be versatile and include free parameters in order to encapsulate the range of different approaches of previous studies and allow for the uncertainties in the underlying physics discussed in section 2.

Cosmic ray sources
In previous studies, Sazonov & Sunyaev (2015) consider only Pop III stars in low-mass halos as cosmic ray sources, whereas Leite et al. (2017) consider all Pop II star-forming halos as sources, and Bera et al. ( 2023) consider both.We encompass all these possibilities in our model.Furthermore, to account for the fact that in some models cosmic rays can only efficiently escape from low-mass halos (≲ 10 7 M ⊙ ), we introduce a parameter  max cr , and restrict the source of cosmic rays in our simulation to a user-specified stellar population(s) in halos of mass less than  max cr .In our simulations, we thus define the cosmic ray emitting star formation rate density field, SFRD cr ( max cr ), as the star formation rate density (in each simulation cell) of objects that emit cosmic rays and are hosted in halos of mass less than  max cr .Due to the uncertainties in the fraction of stars that underwent supernovae, supernova kinetic energy yields, the conversion rate of this energy into cosmic rays, and the escape fraction into the IGM, we convert SFRD cr ( max cr ) into a cosmic ray energy injection rate via a combined efficiency factor  cr .In other words,  cr is formally defined to be the proportionality constant between star formation rate and the rate of cosmic ray energy injected into the IGM.Since massive stars that undergo supernovae have short lives on cosmological timescales, we ignore the time delay between star formation and supernovae in our modelling.Motivated by theoretical predictions, we take the spectrum of cosmic rays injected into the IGM as a power-law with exponent  cr , lower kinetic energy cutoff  K,min and upper kinetic energy cutoff  K,max .We keep  cr ,  K,min and  K,max as free parameters to account for the potential difference between the spectrum of cosmic rays produced at the shock front and that of the rays injected into the IGM.The parameterization also allows us to explore the sensitivity of cosmic ray heating to the cosmic ray spectrum, which was previously found to be significant for variations in  K,min (Sazonov & Sunyaev 2015) and  cr (Leite et al. 2017) 5 .
Thus we compute the injection rate of cosmic rays into the IGM per time interval  and volume  at cell x to be for  K,min ≤  K ≤  K,max , where the normalization factor is defined to be (25)

Energy-loss and propagation
With our source term specified, we now turn our attention to how cosmic ray protons behave inside the IGM.
By combining equations ( 1) and ( 4), accounting for the energy losses to Hubble cooling and excitation and ionization respectively, we can construct a differential equation that describes the evolution of a cosmic ray proton energy in the neutral IGM which we can solve numerically to find the kinetic energy of a cosmic ray emitted at redshift  0 with initial kinetic energy  K,0 .
The above equation describes the temporal evolution of our cosmic ray distribution but we also need the spatial distribution.Given the large uncertainty in how cosmic rays move through the primordial IGM, we consider several alternative models for the spatial distribution of cosmic rays about their source, (x, ;  K,0 ,  0 ).
Physically, the furthest a cosmic ray could have reached from its source is set by its comoving path length , with the cosmic ray free-streaming on a straight line path from its source.This can be calculated by integrating where  K is in turn found by solving equation ( 26).Such a freestreaming scenario could occur if the IGM magnetic field is very weak and so magnetic scattering of cosmic rays is rare.In this model (x, ;  K,0 ,  0 ) would become a shell window function at a comoving distance  from the origin.To gauge the length-scale of this free-streaming propagation, which gives us the maximum possible range of cosmic ray heating, we calculate the comoving path length of protons when they are absorbed  abs .If the protons are not absorbed by  = 6 we calculate their distance (at  = 6) to the emission source.The resulting values are depicted in Fig. 2a.We find lower energy cosmic rays with initial kinetic energy ≲ 2 MeV travel less than the length of one of our simulation cells,  pix = 3 cMpc.Above this energy threshold but below ∼ 200 MeV cosmic rays travel at least a simulation cell length and up to 1000 cMpc while still being absorbed before  = 6.Therefore, we find that in principle cosmic rays can travel far from the source distributing the energy as heat over large cosmological scales (of the order of a few hundred cMpc).
Cosmic rays with higher initial energies (≳ 200 MeV) are not absorbed by  = 6.From our comparison of energy-loss timescales in subsection 2.3, we know for cosmic rays with  K > 200 MeV their energy losses are dominated by Hubble cooling.Hence, while these highest-energy cosmic rays travel large (> 1 cGpc) cosmological distances they are not anticipated to contribute to cosmic ray heating due to their energy being lost to the expansion of the universe.
For cosmic ray spectrum with the theoretically predicted exponent of  cr = 2 the injected cosmic ray energy is split evenly over log initial kinetic energy.Hence, the low energy cosmic rays with  K ≲ 2 MeV are anticipated to contain a significant fraction (30 per cent, using the  K,min and  K,max of Sazonov & Sunyaev 2015) of the total injected cosmic ray energy.Since these cosmic rays can only travel short distances on cosmological scales we would thus anticipate that even in this limiting case of free-streaming a substantial portion of cosmic ray heating occurs within a few cMpc of the source halo, and thus cosmic ray heating is likely strongly clustered around starforming halos.This localized nature of cosmic ray heating would in turn suggest that the globally averaged models used in previous works are insufficient to fully model the impact of cosmic ray heating on the 21-cm signal.For comparison with other studies, we will also consider a spatially uniform case of cosmic ray heating with (x, ;  K,0 ,  0 ) being a constant.However, we note that this is not a physically plausible model of cosmic ray propagation due to its acausal nature.
The opposite logical extreme to free-streaming would be for cosmic rays to not escape their parent halos.This would correspond to no cosmic ray heating of the neutral IGM, and hence no impact on the 21-cm signal (easily modelled by setting  cr = 0).However, cosmic rays can escape the Milky Way (Schlickeiser 2002), a much larger and thus harder to escape galaxy than those expected at high redshifts.In addition, Sazonov & Sunyaev (2015) argued that the Bohm lower diffusivity bound showed cosmic ray heating should extend to distances greater than the average separation between minihalos at  ≤ 20.Similarly, detailed simulations by Jana & Nath (2018) of the diffusion of cosmic rays out of individual halos led them to conclude that cosmic ray heating from minihalos at  = 10 and  = 20 is homogeneous on inter-halo scales, even when the diffusivity of cosmic rays is assumed to be two orders of magnitude less than is seen in the Milky Way today.The physical expectation is thus that even in the most confined scenarios some heating cosmic rays will escape into the IGM.
Previous studies would thus indicate the minimum physically plausible range of cosmic ray heating is comparable to the inter-halo distance at  < 20.For our model of star formation, we find each cell of our simulation should contain several star-forming halos at  ≲ 20 (Reis et al. 2022).Hence for the epochs during which the 21-cm signal is driven by heating, the cosmic ray minimum range is below the resolution of our simulations.As a result, in our shortest plausible range model, hereby called the locally-confined heating model, cosmic ray heating should not travel between neighbouring simulation cells, e.g.(x, ;  K,0 ,  0 ) is 1 in the origin cell and 0 elsewhere.
In the locally-confined heating model the distribution of cosmic ray heating thus needs to be handled by a sub-grid model.Motivated by Jana & Nath (2018) we assume cosmic ray heating is uniform within each cell in our locally-confined heating model.Note this is not to say that the range of cosmic ray heating is the 3 cMpc size of the simulation cell, but instead that the overlapping cosmic ray heating regions around different star-forming halos within the cell leads to approximately uniform heating of the IGM within the cell.However, variations in heating rate between cells are still present due to the differences in matter overdensity and hence star-forming halo density from cell to cell.Furthermore, we assume the heating is partitioned between the ionized and non-ionized phases of each cell in proportion to their masses.In Section 4.4 we discuss the anticipated consequences of these assumptions breaking down, which simulations by Jana & Nath (2018) indicate is likely around large galaxies at redshifts lower than our simulations ( ∼ 4).
Short-ranged cosmic ray heating such as that we are attempting to describe via our locally-confined heating model would occur if the high-redshift magnetic field in the IGM were strong.These strong magnetic fields could either be primordial in nature or originating from weak seed fields amplified by structure formation (Sur et al. 2010) and then expelled into the IGM via SNRs or galactic outflows (Samui et al. 2018).Cosmic rays may even somewhat self-confine via the generation of magnetic fields through instabilities or the Biermann battery (Yokoyama & Ohira 2022), though the former process is potentially inefficient in the neutral IGM due to magnetosonic wave damping (Leite et al. 2017).
With the two propagation mechanisms (the locally-confined and free-streaming models) bracketing the range of physically plausible distributions for cosmic ray heating, we can now compute the cosmic ray number density in a given cell of our simulation.We combine the injection rate, equation ( 24), and the chosen propagation window functions where * is the spatial convolution, and  K,0 ( 0 ,  K ) is the initial kinetic energy computed for a cosmic ray emitted at redshift  0 which has kinetic energy  K at redshift .The function  K,0 ( 0 ,  K ) being defined implicitly as the solution to equation 26.

Heating and Ionization Rate
Finally, we need to convert  2  ( K , )/  K into cosmic ray heating and ionization rates.By using the heating fractions  heat defined in subsection 2.4 we can compute the cosmic ray heating rate per baryon as where −  K /| E&I is the rate of kinetic energy loss by a cosmic ray proton to excitation and ionization interactions.This heating rate can then simply be added to the sum of heating rates in equation ( 13) to model cosmic ray heating in our simulations, with the corresponding ionization rate from equation ( 8) used in equation ( 17).To properly account for the dependence of heating on spatial variation of   (rather than using a globally averaged or fixed value like previous works), we interpolate between the tabulated values of  heat for each cell of our simulation individually.
With our full model established, we can now infer the energy range of cosmic rays that dominate heating.For an individual cosmic ray proton, we consider the fraction of its initial kinetic energy that contributes to the IGM heating by  = 6 (shown in Fig. 2b) assuming the particles travel through the IGM of mean density and ionization fraction.As has been found in previous studies (Sazonov & Sunyaev 2015), the resulting fraction of energy that ends up as IGM heat is strongly dependent on the initial kinetic energy of the cosmic ray proton.Protons with  K,0 < 10 MeV deposit at least 10 per cent of their initial energy, with this contribution rising to above 30 per cent at  < 10 when the free electron fraction increases in the course of reionization.Above  K,0 = 10 MeV there is an initial gradual decrease in the energy fraction; however, once the cosmic rays have high enough  K,0 to not be fully absorbed by  = 6 the fraction rapidly falls below 1 per cent.At the highest energy  K,0 > 1000 MeV less than 0.1 per cent of a cosmic ray initial kinetic energy contributes to heat.If indeed the cosmic ray injected spectrum has an exponent around −2, as is predicted by theory, these results reaffirm the conclusions of Sazonov & Sunyaev (2015) and Leite et al. (2017) that cosmic ray heating is dominated by lower energy particles with < 30 MeV.However, we also find that cosmic rays with much higher energies up to 200 MeV have a non-negligible contribution, heating up the IGM on large scales out to ∼ 100 cMpc away from their sources.

Summary of our cosmic ray model
For clarity let us summarise our final model for cosmic ray heating.In our framework, cosmic rays are emitted from star-forming halos below a mass threshold  max cr at a rate proportional to either the Pop II, Pop III, or the total star formation rate in that halo, as specified by the user.Particles are injected into the IGM with an efficiency of  cr (energy injection per star formation rate), with the injected cosmic ray spectrum taking the form of a power-law in kinetic energy with exponent  cr and lower/upper cutoff of  K,min / K,max .While in the IGM, cosmic rays lose energy to Hubble cooling and excitation and ionization interactions, as described by equation ( 26).Through these processes they heat the IGM at a rate described via equation ( 29), propagating outward via one of three mechanisms: globally uniform, free-streaming, or locally-confined.This propagation is modelled via equation (28) using the relevant transfer function: a uniform distribution, a spherical shell (equation 27), and 1 in the origin cell with 0 elsewhere respectively.Cosmic ray heating rate is then integrated into our simulations by extending the sum in equation ( 13).The model is flexible by construction with free parameters describing the cosmic ray emitting population ( max cr ,  cr ,  cr ,  K,min ,  K,max ) and the propagation mechanisms (globally uniform, free-streaming, locally-confined).
As discussed in Section 3.3.1 the efficiency of cosmic ray injection into the IGM,  cr , is dependant on the fraction of stars that underwent supernovae, supernova kinetic energy yields, the conversion rate of this energy into cosmic rays, and the escape fraction of cosmic rays into the IGM.The first two of these are highly sensitive to the initial mass function of the stellar population of interest.Only higher mass stars undergo supernovae, and the type and hence kinetic energy yield of the supernovae is in turn mass dependent.In this study, we treat  cr as a free parameter and decouple it from our modelling of the Pop III initial mass function due to the large uncertainty in the other components of  cr .Future works that aim to constrain early universe models with 21-cm signal data, could model the link between IMF and  cr using appropriate conditional priors.
While a correspondence cannot always be made between our model parameters and those of other works when it can be it is instructive to consider the parameters used by previous studies.Here for ease of comparison, we attempt to convert the parameters used in various studies to the parameters used in our work, with the values listed in Table .2. Several of the previous studies consider one supernova of a given yield occurring per halo, which is then converted into cosmic ray energy injected into the IGM with some efficiency parameter or a combination thereof.To convert such values to  cr we have assumed that this supernova is occurring in a ∼ 10 6 M ⊙ halo, with a baryonic collapse fraction of 0.1 and a star formation efficiency  * = 0.01.Another difference in modelling is the usage of either kinetic energy or momentum power laws.When previous studies have used momentum power laws, we convert their spectral exponents to equivalent kinetic energy power law exponents by equating the proportion of cosmic ray energy that ends up in cosmic rays with  K < 30 MeV, since we anticipate these cosmic rays to dominate heating.The parameter values listed in Table .2 show that the largest differences between previous works are in the values of  cr which varies between 1×10 47 and 5×10 49 erg M −1 ⊙ .The discrepancy arises due to the uncertainty in the initial mass function of high redshift stellar populations (and thus type/frequency of supernovae), the efficiency of conversion of supernovae kinetic energy to cosmic rays, and the escape fraction of these cosmic rays into the IGM.
A further potential impact of cosmic rays at high redshift is positive feedback on star formation via enhancing the H 2 concentration (Stacy & Bromm 2007;Jasche et al. 2007), however, we do not include this effect in our simulations due to these results being refuted by Hummel et al. (2016).We also do not include the negative feedback from cosmic rays heating gas near the halos virial radius which might suppress the rate of gas accretion onto the halo (Lacki 2015;Jana & Nath 2018), as it is only efficient at low redshifts ( ≤ 4).
Finally and importantly, we do not include excess radio backgrounds from cosmic rays.Although Jana et al. ( 2019) demonstrated that cosmic ray electrons in Pop III SNRs can produce an excess radio background via synchrotron emission, their analysis assumed an ambient magnetic field of 0.32 G at  = 17 (equivalent to ∼ 1 nG comoving).As we discussed earlier, such a strong magnetic field might be too extreme as it is at the upper bound of the extremely broad (11 orders of magnitude) range of the experimentally allowed primordial magnetic field values.The excess radio background found by Jana et al. ( 2019) scales as the magnetic field to the 8/5th power, and so an order of magnitude lower magnetic field would reduce the extra radio background relative to the CMB from 300 per cent to only 8 per cent.If we assume a comoving magnetic field of ≲ 0.01 nG consistent with our treatment of Alfvén wave heating, the excess radio background from cosmic rays should be less than a per cent of the CMB temperature at the redshifts of interest, and thus will have negligible impact on the 21-cm signal (Sikder et 2023).

RESULTS
We now describe our findings regarding the impacts of cosmic ray heating on 21-cm observables, splitting our investigations into four main themes.First, in subsection 4.1 we explore the implications of the short-range nature of cosmic ray heating on the 21-cm signal and how this could be leveraged to distinguish a cosmic ray heated IGM from an IGM heated by other mechanisms.Secondly, we discuss the biases introduced by assuming spatially uniform cosmic ray heating in subsection 4.2.Thirdly in subsection 4.3, we compare the efficiency of cosmic ray heating to X-ray heating.Finally, we contrast the impact of cosmic ray heating with other high-redshift astrophysical processes in subsection 4.5.

Impacts of cosmic ray heating clustering
Previously (subsection 3.3) we arrived at the same conclusions as Sazonov & Sunyaev (2015); Leite et al. (2017); Jana & Nath (2018) that cosmic ray heating is dominated by sub-relativistic protons and, thus, is likely clustered around overdense regions.This is in contrast to heating by X-ray sources (Furlanetto et al. 2006;Fialkov et al. 2014b;Pacucci et al. 2014) which affects the gas temperature on large cosmological scales of several hundred cMpc.To isolate the characteristic patterns of cosmic ray heating in the 21-cm signal from the effects of other processes we compare models with similar IGM thermal histories.
Our main case (reference model) is similar to the model considered by Leite et al. (2017), with cosmic rays emitted from all star-forming halos regardless of population or mass (with parameters  cr = 10 48 erg M −1 ⊙ ,  K,min = 10 −2 MeV,  K,max = 10 9 MeV,  cr = −2) and assuming a free-streaming mode of cosmic ray propagation.We contrast this model with locally-confined cosmic ray heating, X-ray heating with a soft spectral energy distribution (SED), and X-ray heating with a hard SED.In our simulations, a soft X-ray SED is modelled as a truncated power-law with exponent −1.5 and lower cutoff of 0.2 keV corresponding to a mean free path of ∼ 1.3 cMpc at  = 15 (Furlanetto et al. 2006).The majority of X-ray photons in this case have mean free paths of several cMpc and, thus, inject energy relatively close to the sources (although are still longer-range compared to cosmic rays).The hard SED is modelled using a truncated power-law with exponent −1 and lower-cutoff of 3 keV, corresponding to a mean free path of ∼ 4300 cMpc at  = 15.In this case, most photons travel large cosmic distances, lose most of their energy to cosmic expansion and are hardly absorbed by the IGM.Therefore, this model is expected to yield relatively uniform and weak heating.
For ease of comparison, we calibrate heating efficiency parameters (either  cr or  X 6 ) to achieve similar thermal histories between the models.We achieve this by minimizing the root-meansquare difference between the reference global 21-cm signal 7 and the global signals of each other model over the redshift range  = 10 − 25.Through this procedure, we find signals with minimal differences as shown in Fig. 3.We find that the reference model with  cr = 10 48 erg M −1 ⊙ produces a similar global signal and a similar thermal history as the locally-confined cosmic ray heating model with  cr = 1.39 × 10 48 erg M −1 ⊙ , soft X-rays with  X = 0.332 and hard X-rays with  X = 60.0.For comparison, we note that   = 1 is the typical value calibrated to the present-day population of X-ray binaries and taking into account the more efficient X-ray emission in metal-poor environments of the high redshift universe (Fragos et al. 2013).
Closely comparing the matched global signals we see that the high-redshift parts regulated by the onset of star formation and WF coupling ( ∼ 25 − 35) are identical for all the calibrated models.Small differences in the shapes of the signals arise at redshifts affected by heating processes ( ∼ 10−20, we refer to this part of the signal as 6 All other astrophysical parameters (listed in Table 1) and cosmological parameters, including the cosmological initial conditions, are kept the same between the simulations. 7We also considered minimizing the difference between average kinetic temperatures and came to the same conclusions as we do here for matching global 21-cm signals.
Table 2. Parameters used by previous comic ray heating studies converted into our notation, in order: the emitting stellar population, the maximum halos mass cosmic rays can escape from  max cr , the efficiency of cosmic ray emission  cr , the exponent of the cosmic ray spectrum  cr , and the minimum/maximum cutoff of the cosmic ray spectrum  K,min / K,max .The propagation mechanism of cosmic rays is not included as all the previous studies assumed uniform heating, with the exception of Jana & Nath (2018) who considered individual halos.Values marked with a red asterisk are approximate equivalents.If there is no close equivalent value, or no value is given in the corresponding paper, we leave that table entry unpopulated (dash).Two separate sets of parameters are given for Bera et al. (2023) as they model Pop II and Pop III stars separately.For ease of comparison, we also list the ranges of values used in our subsequent investigation broken down by section.

Study
Emitting Pop.The efficiency of cosmic ray/X-ray emission was tuned for the models to minimize the root-mean-square error between their predicted global signals and that of the reference free-streaming cosmic ray heating model with  cr = 10 48 erg M −1 ⊙ .This procedure gave  cr = 1.39 × 10 48 erg M −1 ⊙ for the locallyconfined cosmic ray heating model,  X = 0.332 for the soft SED X-ray heating model, and  X = 60.0 for the hard SED X-ray heating model.All other astrophysical parameters listed in Table 1 and cosmological parameters, including the cosmological initial conditions, are kept the same between the simulations.The four models show good agreement at  > 22 as heating of the IGM has yet to become significant, as well as in the heating arm of the global signal around  ∼ 15.Differences are still seen between models at the global signal minimum and emission maximum.A vertical grey dashed line is shown at  = 16, illustrating the redshift at which we compare tomographic 21-cm maps, see Fig. 4. the heating arm) suggesting that differences in the spatial distribution of heating do not fully average out.In the following, to produce the clearest visual comparison of heating signatures, we contrast signals at  = 16 when the four models have approximately equal global signals ⟨ 21 ⟩ = 82-85 mK.

Imprints of cosmic ray heating in 21-cm tomography
The 21-cm signal is sensitive to gas temperature and, therefore, the character of the IGM heating is expected to be manifested in the spatial distribution of the brightness temperature.To visualize the differences in heating patterns between cosmic rays and X-rays we begin by considering 21-cm tomography.
Fig. 4 shows slices of the gas temperature cubes at  = 16 produced by the four simulations outlined above, alongside the corresponding slices of the 21-cm signal.As anticipated, the maps show that in cosmic ray heated models, the heating of the IGM is much more clustered compared to the cases with X-ray heating.We see a sharp contrast between the localized high-temperature regions and the vast regions of cooler gas, which is manifested in the 21-cm maps as regions of mild emission signal and deep absorption.The contrast is the sharpest for locally-confined cosmic ray heating as is expected given it has the shortest range among the considered scenarios.As we consider increasingly longer-ranged heating mechanisms (i.e.free-streaming cosmic rays, X-ray heating with soft SED, and, finally, hard X-rays), fluctuations in the gas temperature are reduced, and, consequently, the brightness temperature maps show less contrast.Ultimately, in the case of heating by hard X-rays (the rightmost column), the gas temperature is nearly uniform and there are practically no visible emission regions in the 21-cm maps.This trend is explained by the increasingly larger mean free path of the energy carriers (either cosmic ray protons or X-ray photons).For X-ray heating to appear as clustered as cosmic ray heating, the mean-free path of X-rays would need to be at most ∼3 cMpc which at  = 16 (and assuming neutral IGM) corresponds to X-ray energies of ∼ 280 eV or lower.However, high-redshift X-ray sources, such as X-ray binaries, are thought to have much harder SEDs peaking at a few keV (Fragos et al. 2013;Sartorio et al. 2023).Therefore, our simulations suggest that, if observed, the character of the heated regions in the 21-cm maps could be used to probe the nature of the dominant heating mechanism, e.g.distinguish between cosmic ray and X-ray heating.

Signatures of heating clustering in the 21-cm power spectrum
The 21-cm power spectrum is a signal targeted by radio interferometers including HERA (Abdurashidova et al. 2022), LOFAR (Mertens  et al. 2021), and the future SKA (Koopmans et al. 2015).Examples of the 21-cm power spectra for the models explored here are shown in Fig. 5.We start by considering the redshift evolution of the power spectra at a fixed comoving wavenumber  = 0.1 cMpc −1 (panel a).
In general, we expect to see two well-defined peaks in the power spectra vs redshift, one marking the dominance of fluctuations imprinted by the non-uniform WF coupling and the other reflecting the gas temperature fluctuations prevailing at lower redshifts.We find the power spectra produced by the two cosmic ray heating scenarios to be very similar while being dramatically different from the X-ray heating scenarios.As a result of their localized strong heating, the two cosmic ray models show clear narrow peaks due to the WF coupling at  ∼ 25 and heating at  ∼ 15.In the soft X-ray SED case heating is longer-range which results in a slight delay in energy injection, a broader WF coupling peak and a factor of 1.8 lower heating peak.Finally, in the hard X-ray SED case the Lyman-coupling peak is even broader and the heating peak is completely erased owing to the large mean free path of X-ray photons.
Next, we compare the shape of the signals as a function of wavenumber  at a fixed global signal value.For convenience, we choose the value of ⟨ 21 ⟩ = −84 mK8 (which for all the models happens at  ∼ 16) to show results for in panel b.The differences seen in the shapes of the power spectra for various heating mechanisms are much more striking than what we saw for the corresponding global signals and can be understood in terms of the typical length scale below which the structure is washed out.In the locally-confined cos-mic ray heating, this length scale is at the resolution limit of our simulation and so the spherical 21-cm power spectrum traces the matter power spectrum across our entire  range.For cosmic ray free-streaming the heating length scale is larger leading to a suppression of the power spectrum above  ∼ 0.2 cMpc −1 , which is evident from the figure.In the case of X-ray heating, the high energy particles (present even in the soft SED case) affect the shape of the signal at a broad range of scales leading to a suppression of power.The effect is visible even at the largest scales (smallest ) considered, though the suppression remains strongest at smaller scales (high ).Finally, since for the hard SED X-ray model the heating is nearly uniform, the strong suppression of the power spectrum is observed at all the considered scales.In summary, we find that the strongest differences in the signatures of heating mechanisms are produced at high , with the locally-confined cosmic ray heated model having a 60 per cent, 420 per cent, and 1080 per cent larger power spectrum than the free-streaming cosmic ray, soft X-ray SED, and hard X-ray SED models respectively at  = 0.5 cMpc −1 .
Our results suggest that a measurement of the 21-cm power spectrum at  > 0.2 cMpc −1 during the heating-driven stage ( ∼ 10−20) could be used to probe the dominant IGM heating mechanism.Stronger signals at high  values would indicate a more clustered heating, allowing us to distinguish X-ray from cosmic ray heating, and (if the environment was cosmic ray heated) the degree of confinement of cosmic ray particles.To illustrate the experimental feasibility of using the 21-cm power spectrum as such a diagnostic tool we have included 1000 hr SKA thermal noise sensitivities (Koopmans et al. 2015) in Fig. 5, which yield signal-to-noise-ratios > 9 at  = 0.5 cMpc −1 between each pair of models.In the heating-driven era of the 21-cm signal the suppression of the 21-cm power spectrum above a length-scale set by the dominant heating mechanism Comparison of the 21-cm power spectrum for the locally-confined cosmic rays, free-streaming cosmic rays, soft and hard X-rays (models are specified in subsection 4.1.1).Panel (a) shows the redshift evolution of the 21-cm power spectrum at a fixed wavenumber  = 0.1 cMpc −1 .Panel (b) shows the power spectra versus wavenumber  at the redshift when the corresponding global 21-cm signal equals to ⟨ 21 ⟩ = −84 mK (which occurs at  ∼ 16 for all the considered cases).We find that the small-scale (high ) 21-cm signal varies significantly between the models, suggesting that the variation of the 21-cm power spectrum with  could provide a diagnostic tool for the dominant heating mechanism.Also shown in both panels is the thermal-noise estimate for the SKA with 1000 hrs of observations (grey curve, Koopmans et al. 2015), illustrating the theoretical sensitivity of next-generation 21-cm power spectrum experiments.
is a generic effect, previously proposed by Fialkov et al. (2014b) as a way to probe the SED of X-ray sources.Thus the potential to distinguish the nature of high-redshift heating mechanisms based on the shape of the high  21-cm power spectrum should be robust to astrophysical uncertainties, though the exact signal-to-noise ratios will depend on the astrophysical scenario.Additionally, this analysis does not consider experimental systematics, or any potential degeneracies with new physics such as properties of dark matter (Sitwell et al. 2014;Barkana 2018;Muñoz et al. 2018;Fraser et al. 2018;Fialkov et al. 2018;Liu et al. 2019;Muñoz et al. 2020;Jones et al. 2021;Hibbard et al. 2022;Barkana et al. 2023), all of which could weaken the ability for the high  power spectrum to pin down the dominant IGM heating mechanism.Given these uncertainties, further work is required to reliably evaluate the statistical significance with which the nature of high-redshift heating mechanisms can be determined from the projected 21-cm power spectrum measurements of upcoming experiments.

Uniform heating assumption and the global signal
When attempting to match thermal histories between simulations (Fig. 3), we found that an exact correspondence between the produced global signals could not be achieved in models with locallyconfined and free-streaming cosmic ray heating.The inevitable discrepancy suggests that the spatial distribution of IGM heating is reflected in the global 21-cm signal due to ⟨ 21 ⟩ being a non-linear tracer of  K .Given that previous studies assumed uniform cosmic ray heating, rather than the theoretically expected highly clustered behaviour, such an assumption may have biased the global signal predictions in these works.
Here we investigate the degree to which the global signal is affected by the locality of heating.Unlike in the previous subsection, here we do not attempt to match thermal histories between models, instead, we keep the cosmic ray emission efficiency parameter  cr = 10 48 erg M −1 ⊙ constant so that the energy injected into cosmic rays that reach the IGM is the same between simulations.This allows us to isolate the effect of the heating distribution by only varying the cosmic ray propagation mode between locally-confined, freestreaming and globally uniform heating.All other parameters are kept the same as for the reference model used in the previous subsection (and listed in Table 1).The 21-cm global signals predicted for the three different cosmic ray propagation modes are shown in Fig. 6.
As expected, the differences in the predicted global signals with various heating mechanisms are manifested at redshifts affected by heating ( ≲ 20), while the signals are identical at higher redshifts dominated by the WF coupling.In the case of uniform heating, the global signal evolution happens faster, preceding the other two models by Δ ∼ 1.The uniform heating model also features a higher and earlier emission maximum of ⟨ 21 ⟩ = 17.9 mK, at  = 11 compared to ⟨ 21 ⟩ = 12.8 mK, at  = 10 and ⟨ 21 ⟩ = 8.7 mK, at  = 10 in the free-streaming and locally-confined cases respectively.While the effect is smaller we also observe similar differences between the free-streaming and locally-confined propagation modes.
Even though the energy of cosmic rays reaching the IGM is the same in all cases, our results suggest that clustered heating is less efficient at increasing the average 21-cm signal than uniform heating.This is contrary to what would be expected from the mean IGM kinetic temperature, which is found to be highest between  = 13 and 20 for the locally-confined model, due to the correlations between star formation and  e which increase the average efficiency with which cosmic ray energy is deposited into the IGM as heat (  heat ).The non-linear relationship between ⟨ 21 ⟩ and  K explains the fact that uniform heating is the most efficient at increasing ⟨ 21 ⟩.When the WF coupling is efficient (  ≫ 1) but before reionization becomes significant (i.e. HI ≈ 1), the 21-cm signal scales as  21 ∼ 1 −  / K .From this concave functional dependence, we see that  21 saturates at a positive value for  K ≫   but can take large negative values for  K ≪   .Consequently, when a volume average of  21 is taken the weighting favours low  K .As a result, owing to the concentration of heating into a small volume (that then saturates in 21-cm emission), models with clustered IGM heating are less efficient at raising the mean  21 signal than the models in which heating is distributed evenly across the entire IGM.This explains our findings that the Locally-confined Free-streaming Globally uniform Figure 6.Global 21-cm signal for different cosmic ray propagation models.All three models have the same cosmic ray emission efficiency and astrophysical parameters, and thus the energy in cosmic rays reaching the IGM is the same between the models.At redshifts affected by heating ( ∼ 12 − 20), the signal in the globally uniform cosmic ray heating model evolves faster and reaches a higher emission peak than the other two propagation models.Similar but smaller differences can also be seen between the free-streaming and locally-confined models.These differences suggest that clustered heating is less efficient than diffused heating at driving the 21-cm global signal into emission.
shorter range cosmic ray heating we consider the less effective the heating appears to be at raising ⟨ 21 ⟩.
The assumption of uniform cosmic ray heating has thus somewhat biased the predictions of previous studies, with cosmic ray heating appearing too effective at increasing the 21-cm global signal.It should be noted, however, that this effect is moderate with the largest discrepancy being 34.6 mK at  = 15, between locallyconfined and globally uniform heating mechanisms, comparable to the 25 mK EDGES residuals (Bowman et al. 2018), though at a potentially measurable level for future experiments (e.g. the 5 mK sensitivity level projected for 2500 hours of REACH observations, de Lera Acedo et al. 2022).Hence, this bias while present is not anticipated to greatly impact the conclusions of previous studies.

Comparison of heating efficiencies
So far we have considered simulations with either cosmic ray heating or X-ray heating active.However, in reality, both mechanisms will take place simultaneously.Here we contrast the two processes to gain insight as to under which circumstances each process would have a dominant contribution to thermal history.We consider three simulations modelling cosmic ray heating only, with inefficient ( cr = 10 47 erg M −1 ⊙ ), standard ( cr = 10 48 erg M −1 ⊙ ), or efficient ( cr = 10 49 erg M −1 ⊙ ) cosmic ray emission, and a further three simulations modelling X-ray heating only, with inefficient (  X = 0.1), standard (  X = 1), or efficient (  X = 10) X-ray emission which use the X-ray SED calculated by Fragos et al. (2013) for early universe X-ray binaries.One should keep in mind that the selected efficiencies are illustrative and both heating processes are subject to orders of magnitude uncertainties.For comparison, we also consider a seventh simulation with both standard efficiency cosmic ray heating and standard X-ray heating modelled.All other parameters of these simulations are the same as the reference simulation from subsection 4.1.1.The predicted evolution of the volume-averaged IGM temperature outside of ionized bubbles for these simulations is shown in Fig. 7. ⊙ , 10 48 erg M −1 ⊙ , and 10 49 erg M −1 ⊙ respectively; inefficient, standard, and efficient X-ray heating, corresponding to  X = 0.1, 1, and 10; as well as a model with both standard efficiency cosmic ray heating ( cr = 10 48 erg M −1 ⊙ ) and X-ray heating (  X = 1).We find that at standard efficiencies considered in the literature, cosmic ray heating heats up the neutral IGM more rapidly than X-rays.However, as both mechanisms at high redshifts are very uncertain, we advocate for the inclusion of both.
We find that for the range of typical efficiencies considered in the literature, cosmic rays are roughly as efficient as X-rays in heating up the neutral IGM at  < 20.Comparing the cosmic ray and X-ray standard efficiencies that we adopt here, cosmic rays are found to be slightly more efficient than X-rays at raising the temperature of the IGM.Furthermore, we see the temperature for the simulation with both cosmic ray heating and X-ray heating modelled is noticeably above that of the simulations with the standard efficiency heating mechanisms modelled individually.This fact is a manifestation of the comparable heating rates and shows that for some plausible efficiency parameter values, neither heating mechanism dominates and modelling of both is required to produce an accurate thermal history.This comparison indicates the potential importance of cosmic ray heating to the modelling of the 21-cm signal.

Sub-grid clustering of cosmic ray heating
In the above, motivated by the findings of Jana & Nath (2018), we have assumed that cosmic ray heating is uniform within simulation cells.However, this assumption is not expected to hold around massive galaxies at  ∼ 4 (Jana & Nath 2018) and may break down at  > 20 when there is a greater distance between star-forming halos.If the assumption were not valid, we would be artificially smoothing cosmic ray heating, and thus the 21-cm signal, on the scale of the resolution of our simulations.It is thus necessary to consider the impact on our conclusions of clustering of cosmic ray heating on scales smaller than our resolution.
There are two main impacts of sub-grid heating clustering to consider the inhomogeneity of the temperature of the neutral IGM within a cell, and the proportion of cosmic ray heat deposited into the ionized IGM.In section 4.2, we previously discussed how increased inhomogeneity in the temperature of the IGM reduces the effectiveness of a heating mechanism at raising the 21-cm global signal, thereby delaying the rise out of the global signal absorption trough (see Fig. 6).The theoretical explanation for this phenomenon ultimately stems from the non-linear relationship between  21 and  K , and so will also apply to these subgrid scales.Hence sub-grid clustering of heating would lead to a slower rise of the 21-cm global signal post absorption minimum compared to our current models with the same  cr .Consequently, a higher  cr would be needed to match cosmic ray heating scenarios to X-ray heating scenarios as we did in Section 4.1.1.Additionally, we would expect the trends seen when comparing 21-cm tomographic maps in Section 4.1.1 to continue, with more clustered heating leading to greater contrast in the 21-cm maps around radiative sources and larger voids with nearly uniform 21-cm signal.As a result, the wavenumber above which the 21-cm power spectrum is suppressed would now be set by the sub-grid length scale at which cosmic ray heating is clustered, leading the 21-cm power spectrum to trace the matter power spectrum to even higher wavenumbers than were seen in our locally-confined case.
If heating is clustered around sources such that the heating of the IGM between them is not approximately uniform, then we would also anticipate a greater proportion of the cosmic ray energy to be deposited in the ionized regions close to the sources.This is in contrast to our current assumption that heating in a cell is split between the ionized and non-ionized IGM in proportion to their masses.Since fully-ionized regions have zero 21-cm signal this would further diminish the efficiency of cosmic ray heating at raising the 21-cm signal from its minimum.However, providing that cosmic rays can still escape the ionized regions surrounding galaxies the length scale at which the 21-cm power spectrum begins to deviate from tracing the matter power spectrum will still be set by the heating length scale.As a result, as long as there is still a well-defined heatingdriven era of the 21-cm signal the 21-cm power spectrum should provide a diagnostic tool for probing the heating mechanism through its characteristic length scale.
Hence, overall we find that should cosmic ray heating not be uniform within cells we would expect a greater contrast in 21-cm signal tomographic maps near sources and larger unheated voids with nearly uniform 21-cm signal away from sources.In addition, we anticipate there would be a slower increase in the 21-cm global signal, and the 21-cm power spectrum would deviate from the matter power spectrum at higher wavenumbers.These three changes would serve to enhance our conclusions that X-ray heating and cosmic ray heating can be distinguished using the small-scale 21-cm power spectrum, and that assuming globally uniform cosmic ray heating makes cosmic ray heating erroneously efficient at raising the 21-cm signal.

Sensitivity to astrophysics
In the previous subsections, we considered models with a fixed set of astrophysical parameters listed in Table 1, and only varied the type of the IGM heating, cosmic ray propagation model, and emission efficiencies.However, both the cosmic ray heating mechanism and the astrophysical parameters of the early universe are highly uncertain.In this subsection, we consider how the strength and behaviour of cosmic ray heating vary when varying other astrophysical parameters of the model.Throughout this subsection, all simulations have cosmic ray heating, Ly  heating and CMB heating enabled, X-ray heating disabled, and use the free-streaming propagation model for cosmic rays.

Diffusive escape
We start by considering models of cosmic ray heating similar to that proposed by Leite et al. (2017), with cosmic rays emitted by all stellar populations from all star-forming halos assuming  cr = 10 48 erg M −1 ⊙ ,  K,min = 10 −2 MeV,  K,max = 10 9 MeV and  cr = −2.As all stellar populations are emitting cosmic rays in this model, we anticipate heating to be dominated by Pop II stars given their greater star formation rate at later times.We thus consider the variation of the 21-cm signal with the parameters that determine the efficiency and timing of Pop II star formation,  * ,II and  recov respectively.In addition, we seek to verify the finding of Leite et al. (2017) that cosmic ray heating is quite sensitive to the spectral exponent of the cosmic ray spectrum by varying  cr as well.To isolate the impacts of the variation in each one of the aforementioned parameters we show  K and the 21-cm signal as separate columns in Fig. 8 for different values of  * ,II ,  recov and  cr .
As expected, increasing the Pop II star formation efficiency increases the rate of cosmic ray heating in our simulations.Thus as  * ,II rises from 0.002 to 0.2 we find the redshift of  K − cmb equality moves to higher redshifts, from  = 9.9 to  = 15.4,and the temperature of the IGM at  = 8 increases from 40 to 1070 K.This faster heating of the IGM is reflected in the 21-cm signal with higher  * ,II leading to a shallower absorption trough shifted to higher redshifts and a stronger emission peak.Similarly, as  * ,II increases the power spectrum heating peak shifts to higher redshifts and is stronger.The effect of  recov is found to be more modest.Higher values of the recovery time delay the onset of IGM heating by Δ ∼ 1.As a result, for the range of  recov considered here, the IGM always reaches a similar temperature of 257 to 291 K at  = 8.These small shifts in the IGM thermal history lead to small changes in the 21-cm global signal and power spectrum with  recov = 100 Myr having a global signal minimum and power spectrum heating maximum at corresponding delays of Δ = 1 and Δ = 2 compared to the equivalent features for  recov = 10 and 30 Myr.The discrepancies in the 21-cm signals for the explored  recov values are similar to, if smaller in magnitude than, the differences seen by Magg et al. (2022b) when they considered the same parameter variation assuming X-ray heated models.The smaller differences observed in our work can be attributed to the cosmic ray heating being less efficient than the X-ray heating in Magg et al. (2022b).Additional effect comes from our updated LW feedback prescription discussed earlier (see subsection 3.2).
Our simulations are found to support the finding of Leite et al. (2017) that cosmic ray heating is quite sensitive to the spectral exponent  cr of the cosmic rays injected into the IGM.The aforementioned  K > 250 K at  = 8 is achieved for a spectrum flat in log( K ) with  cr = −2.As  cr increases a greater portion of the cosmic ray energy is in the form of higher energy cosmic rays ( K > 200 MeV), which as we found previously do not efficiently transfer their energy to the neutral IGM.Consequently, the efficiency of cosmic ray heating decreases quite sharply as  cr increases, with  K only exceeding  CMB at  = 7.4 in the case of an exponent with  cr = −1.7.The decrease in heating efficiency with  cr leads to deeper absorption troughs shifted to lower redshifts and the gradual elimination of the power spectrum heating peak.
The strong dependence of the efficiency of cosmic ray heating on the spectrum of the injected cosmic rays suggests that accurate modeling of the escape mechanism is required to understand the role of cosmic rays in thermal history.Such simulations are made more challenging by the uncertainty in the strength of magnetic fields in the early universe and the nature of the astrophysical population in the first star-forming halos (Jana & Nath 2018).Detection of cosmic ray heating through 21-cm observations would thus potentially provide insight into the magnetic field environment within the star-forming  (left column), recovery time between Pop III and Pop II star formation  recov (middle column), and exponent of the cosmic ray spectrum  cr (right column). * ,II is found to have a strong impact on the rate of the IGM heating and thus the 21-cm signal.We find an increase in the star formation efficiency from 0.002 to 0.2 results in the IGM temperature at  = 8 raising from 39.7 to 1070 K. Conversely,  recov is found to only have a small impact.Larger values of  recov delay Pop II star formation and thus the IGM heating, leading to Δ ∼ 1 shifts in the 21-cm global signal and power spectrum across the range of  recov considered.The depicted variations with  cr replicate the findings of Leite et al. (2017) that the cosmic ray heating rate is quite sensitive to the spectral shape of the cosmic rays injected into the IGM, with increasing  cr from −2 to −1.7 reducing the  = 8 IGM temperature from 284 to 19 K, and annihilating the power spectrum heating peak.
regions prior to reionization, as we would be able to test whether or not the low-energy cosmic rays were able to escape from their parent halos.

Direct injection
Until now we have assumed cosmic rays can escape from halos of any mass following the diffusive escape model discussed in subsection 2.2.However, Sazonov & Sunyaev (2015) proposed an intriguing alternative wherein cosmic rays are released directly into the IGM due to the energetic SNRs of Pop III stars escaping their fully photoevaporated host halos.Since in this mechanism the SNR has to escape the host halo for cosmic rays to be directly injected into the IGM, this sets an upper limit on the mass of star-forming halos that contribute to cosmic ray heating, which in their work Sazonov & Sunyaev (2015) take as 10 7 M ⊙ .
To allow for easier comparison between our findings and those of Sazonov & Sunyaev (2015) in this subsection we assume cosmic rays can only be injected into the IGM by Pop III stars, turn off X-ray heating, and adopt the same values for  K,min = 10 −3 MeV,  K,max = 10 8 MeV, and  cr = −2 as in their work.As reference values we take  max cr = 10 7 M ⊙ , to match the aforementioned study, and use  cr = 3 × 10 48 erg M −1 ⊙ , which is in the middle of the range considered by Sazonov & Sunyaev (2015).Since in this model cosmic rays are sourced from Pop III star supernovae, we consider the variation in cosmic ray heating with Pop III star formation efficiency  * ,III and the strength of the LW feedback (set by the delay parameter  LW , see Fialkov et al. 2013) due to the LW feedback impacting the minimum halo mass in which Pop III stars can form.Finally, we also consider the sensitivity of this model to its unique feature, the halo mass threshold below which cosmic rays can escape into the IGM  max cr .Similarly to the last section, we calculate the IGM kinetic temperature, global 21-cm signal, and power spectrum for variations of these three parameters and show the results in Fig. 9.
For many of the parameter values we consider, including our fiducial values, an unexpectedly complex redshift evolution of the IGM kinetic temperature occurs.Rather than the familiar picture in which the gas temperature is driven by monotonic heating preceded by adiabatic cooling (as we found in the previous subsection Fig. 8), the direct injection models develop a local peak in the gas temperature at  ∼ 10 − 20.In scenarios with the aforementioned complex  K evolution, we find adiabatic cooling initially dominates the IGM temperature evolution (as expected).Around the first minimum in gas temperature ( ≳ 20) cosmic ray heating becomes efficient causing the kinetic temperature to rise.However, soon after cosmic ray heating becomes efficient,  mol,crit rises above  max cr due to the LW feedback (which was not modelled in Sazonov & Sunyaev (2015), hence they did not observe this unusual behaviour).Consequently, cosmic rays can no longer escape from star-forming halos, and so cosmic ray heating of the IGM is cut off.Adiabatic cooling thus begins to dominate the IGM temperature evolution for the second time in cosmic history, the IGM temperature drops creating a local maximum in the thermal history.Finally, Ly  heating becomes efficient at later times causing the secondary minimum in  K as the IGM temperature rises.Ly  heating remains efficient for the rest of our simulation resulting in a further monotonically increasing gas temperature.
In Fig. 9 we explore how this peculiar behaviour changes with the efficiency of the LW feedback.The differences primarily come from the fact that  LW regulates the redshift at which  mol,crit exceeds the maximal halo mass for cosmic ray escape.As we increase  LW the LW feedback is stronger and so  mol,crit exceeds  max cr earlier.Consequently, the peak in gas temperature is shifted to earlier times and is lower due to the shorter phase of cosmic ray heating.Curiously, the presence of the peak in  K is manifested in the global 21-cm signal as an additional inflexion point (marked by black crosses in Fig. 9).We find that an earlier and weaker peak in  K in simulations with larger  LW values leads to the inflexion point moving to higher redshifts.In addition, models with stronger LW feedback (larger  LW ) produce deeper global 21-cm signals due to less efficient cosmic ray heating.This effect is also manifested in the power spectra as enhanced power at and around the heating peak.
Next, we explore how changing  max cr affects the shape of the signals.Decreasing the free parameter  max cr increases the redshift at which  mol,crit =  max cr thus leading to an earlier cut-off of cosmic ray heating.For the lowest explored value of  max cr ∼ 10 6 M ⊙ , cosmic ray heating is cut-off very early (at  > 20) and before its rate can exceed the adiabatic cooling rate.As a result, in this model, the IGM continues to cool down adiabatically until Ly  heating becomes efficient at  ∼ 13.The resulting thermal history is standard (without visible local maximum/minima).This picture qualitatively changes as  max cr rises above ∼ 10 7 M ⊙ .In such scenarios, cosmic ray heating is efficient for a sufficient amount of time to imprint changes in thermal history.At extremely high values of  mol,crit (at and above ∼ 10 8 M ⊙ ) cosmic ray heating becomes the dominant mechanism and is not cut off even by  = 8.In this simulation the resulting IGM temperature is monotonically rising after the onset of cosmic ray heating without developing any unexpected features.The resulting 21-cm global signal and power spectrum are also quite sensitive to the choice of  max cr with low values ( max cr = 10 6 M ⊙ ) resulting in a deep and late absorption trough (−158 mK at  = 17).
Owing to the absence of efficient cosmic ray heating in this case the heating peak in the power spectrum is very weak.For the intermediate values ( max cr ∼ 10 7 M ⊙ ) we find a more complex behaviour with additional inflexion points in the global 21-cm signal between  = 20 and 10 and an enhanced power spectrum peak.Finally, as expected from the thermal history, the highest explored values (10 8 M ⊙ ) with efficient heating result in an earlier and shallower absorption trough (−111 mK at  = 20) and a well-defined power spectrum heating peak.
The impact of  * ,III on the 21-cm signal is more complex than for the other two explored parameters since it affects both the rate of cosmic ray heating and the strength of the LW background.In our simulations, we find the local peak in  K occurring for models with  * ,III ≳ 0.02, while models with lower values of  * ,III feature a flattened extended  K minimum.This dependence illustrates that the stronger LW feedback, and thus the earlier cut-off of cosmic ray heating (e.g., from more vigorous Pop III star formation), is not sufficient to fully counteract the more efficient cosmic ray heating resulting from the higher Pop III star formation efficiency.As a result, we observe an earlier but higher  K peak as  * ,III increases (in the regime  * ,III ≳ 0.02).Consequently, higher  * ,III values lead to earlier and slightly shallower absorption troughs with inflexion points at progressively higher redshifts owing to the increasingly more efficient cosmic ray heating.On the other hand, we find that inefficient Pop III formation with the values  * ,III < 0.005 results in a deep global 21-cm signal with no additional inflexion point.The power spectrum trends are somewhat simpler, with higher  * ,III leading to earlier and stronger peaks from both the increased Ly  emission (peak at  > 20) and the more efficient heating (peak at 12 <  < 17).We observe no apparent new behaviour in the power spectra, unlike in  K and ⟨ 21 ⟩, around  * ,III ∼ 0.02.
We find that indeed direct injection of cosmic rays can heat up the IGM.However, we do not find this particular model of cosmic ray heating to be as efficient as in the original study by Sazonov & Sunyaev (2015).In part, this can be attributed to differences in star formation prescription, the relaxation of the global heating assumption that we introduced here, and the truncation of cosmic ray heating by the LW feedback.Consequently, whereas Sazonov & Sunyaev (2015) found cosmic rays to heat the IGM by 10 to 100 K at  = 15, in our simulations 10 K of heating at  = 15 is only achieved for high values of star formation rates (  * ,III ≥ 0.02) or high escape mass thresholds ( max cr ≥ 10 7 M ⊙ ) with the largest  = 15 temperature increase above the adiabatic cooling solution Δ K = 20.3K achieved for  max cr = 10 8 M ⊙ .While the amount of heat deposited into the IGM by cosmic rays is found to be smaller than in Sazonov & Sunyaev (2015), the 21-cm signal is still sensitive to this mechanism.Our results suggest that the 21-cm signal could provide a probe of the nature of Pop III star supernovae through the small amount of cosmic ray heating which they produce around  = 15 − 20.
To recap, we find that the direct injection model can lead to a rich IGM thermal history and unusual features in the 21-cm signal.Above we showed that the shape of the global signal undergoes a qualitative transition as the maximal halo mass changes from  max cr = 10 6 to 10 7 M ⊙ .A new inflexion point is developed as a manifestation of the complex heating history owing to the contribution of cosmic rays.Now we consider this transition in more detail by densely sampling the parameter  max cr .We run the code 50 times with equal logarithmic spacing in  max cr in the range 10 6 − 10 7 .The results for  K and ⟨ 21 ⟩ are shown in Fig. 10.
Examining the resulting thermal histories we find a clear peak in temperature for  max cr = 10 7 M ⊙ , the upper end of our considered range.As  max cr decreases, the peak shifts to higher redshifts and cr (right column).Some of the models are found to display a complex thermal history, with a maximum and a secondary minimum in  K .The height and redshift of this  K maximum are found to be quite sensitive to  LW and  max cr .This is due to the maximum being caused by  mol,crit (the critical halo mass for molecular cooling star formation) being increased above  max cr by the LW feedback.Due to the LW feedback cutting off cosmic ray heating, we do not find it to be as efficient as was found in the previous study (Sazonov & Sunyaev 2015).These complex IGM thermal histories in turn produce unusual inflexion points in the heating arm of the 21-cm signal (marked with black crosses), which we explore further in Fig. 10.lower temperature values.Eventually, the local maximum merges with the higher redshift minimum transforming into two inflexion points that can be seen in the highlighted signal (black line) corresponding to  max cr = 4.57 × 10 6 M ⊙ .For lower values of the critical mass, thermal evolution shows the classical behaviour driven by the adiabatic cooling and Ly  heating, with no extra features and no clear sign of cosmic ray heating.
Considering the effect on the global 21-cm signal, we see that as  max cr decreases the inflexion point shifts to higher redshifts.Eventually, the inflexion point merges with the minimum of the absorption trough producing a flattened profile (black line in the right panel of Fig. 10).As we reduce the value of  max cr even further, the absorption trough becomes deeper, narrower, and is no longer flattened.
Intriguingly, the flattened absorption feature that we find for  max cr = 4.57 × 10 6 M ⊙ has a similar shape to the disputed EDGES Low Band detection of the global signal (Bowman et al. 2018) which has a flattened minimum at a similar redshift.We note, however, that cosmic ray heating cannot provide a self-consistent explanation of the EDGES best-fit signal on its own as the resulting absorption trough does not have the required depth or the steep sides observed by EDGES.However, it does illustrate that such flattened global signals are not inherently pathological and can be achieved in models with two separate heating mechanisms, one active at higher redshifts and only mildly pre-heating the IGM and the other representing a more sustained heating that becomes efficient at later times (such as X-ray heating, Ly  or CMB heating).The high-redshift preheating could be achieved also in other models such as e.g., fractional annihilating dark matter with a cross-section that decreases with time (Liu & Slatyer 2018).

DISCUSSION AND CONCLUSIONS
In this paper, we have developed a model of cosmic ray heating for use in semi-numerical 21-cm signal simulations.Unlike previous globally-averaged models our approach includes realistic modelling of the spatial distribution of cosmic ray heating.We bracket the range of possible cosmic ray heating scenarios by considering freestreaming versus locally-confined models of particle propagation, exploring the dependence of heating on the cosmic ray energy spectrum, and considering different channels for escape from host halos.Our method thus allows us to perform the first investigations of the signatures of this heating mechanism in 21-cm tomographic maps and power spectra, as well as the previously studied global 21-cm signal.It also enables us to compare the contribution of cosmic rays to the thermal history of the IGM to other heating mechanisms, such as X-ray heating, and thus to propose novel methods to potentially distinguish an X-ray heated IGM from a cosmic-ray heated IGM using upcoming observations.We find the spatial range of cosmic ray heating to be much shorter than that of X-ray heating, even when cosmic rays are assumed to travel along straight lines.This short-range nature manifests in the IGM heating being clustered around regions of efficient star formation and results in a sharp contrast in the 21-cm signal between localized heated regions of emission and vast unheated regions of absorption.
Comparing our inhomogeneous cosmic ray heating model to a uniform model commonly used in the literature, we find that its localized nature results in a slightly delayed evolution of the global 21-cm signal.This demonstrates the existence of a small assumption-induced bias in the results of previous globally-averaged studies.Much larger differences are seen in the power spectra and tomographic maps, as with cosmic rays as a dominant heating source the 21-cm signal traces matter fluctuations to smaller scales than in the case of X-rays.Specifically, we find that the 21-cm power spectrum at higher  is increasingly suppressed for the longer-ranged heating mechanisms allowing us to constrain the propagation mechanism of cosmic ray particles and distinguish cosmic ray from X-ray heating.Any clustering of cosmic ray heating on scales smaller than our simulation resolution would enhance these effects and thus strengthen our conclusions.The potential of the high- end of the power spectrum to provide information about the spatial range of the dominant heating mechanism is anticipated to apply more generally and be robust to astrophysical uncertainties, making it a useful diagnostic tool.For example, it has previously been suggested (Fialkov & Barkana 2014;Fialkov et al. 2014b) that the high- end of the power spectrum could be used to probe the SED of early X-ray binaries.
Due to the uncertainties surrounding cosmic ray escape into the IGM, we explored several plausible mechanisms.In models where cosmic rays can escape from halos of any mass via diffusion or outflow advection, we found cosmic ray heating of the IGM can be more efficient than X-ray heating.However, similarly to Leite et al. (2017), we found that the efficiency of cosmic rays as a heating source is strongly dependent on the spectrum of the particles injected into the IGM.A change in the source kinetic energy spectral exponent from −2 to −1.7 decreases the IGM temperature from 284 to 19 K at  = 8.This is because only the lower energy cosmic ray protons are able to efficiently transfer their kinetic energy to the IGM as heat, as was previously demonstrated in Sazonov & Sunyaev (2015).
For cosmic ray heating models similar to that of Sazonov & Sunyaev (2015), wherein cosmic rays can only escape from lower mass Pop III star-forming halos, we found cosmic ray heating to be relatively inefficient.Due to Sazonov & Sunyaev (2015) not considering the LW feedback in their model, we were only able to replicate the significant increase in temperature they observed for unphysically weak LW feedback or if we raised the threshold halo mass of cosmic ray sources.Instead at theoretically motivated parameter values we find the interplay between this mass threshold and the LW feedback, which raises  mol,crit , leads to unexpectedly complex thermal histories.In such models, the IGM kinetic temperature was found to have a local peak due to cosmic ray heating being cut off by the LW feedback leading to a secondary episode of cosmic cooling before Ly  heating from Pop II stars becomes efficient and re-heats the IGM.This unusual thermal evolution results in an inflexion point in the heating arm of the global 21-cm signal.By fine-tuning the mass threshold for cosmic ray escape, this inflexion point can be merged with the minimum of the absorption trough resulting in a flattened profile.Thus, we demonstrate that flat-based global 21-cm signals (like the one detected by EDGES) can in principle be achieved in models where one heating mechanism becomes inefficient at a similar time to the gas becoming fully coupled via the WF effect and another heating mechanism turns on at a later time.We do not advocate for this model to be a self-consistent explanation of the EDGES best-fit signal as it does not explain the depth or the steep sides of the detected absorption trough (although a dual-heating scenario could form a part of a more complex explanation which also invokes a radio contribution to create the deep absorption).
There remains a major uncertainty in cosmic ray heating, namely the strength of the IGM primordial magnetic field which scatters cosmic rays.This uncertainty leads to a large uncertainty in the cosmic ray diffusivity within the IGM.Furthermore, in the case of primordial magnetic fields at the higher end of the experimentally allowed range (> 0.01 nG comoving field strength), the transfer of energy from cosmic rays to the IGM via Alfvén waves (Bera et al. 2023) and synchrotron radio emission from cosmic ray electrons (Jana et al. 2019) are anticipated to become efficient.Since these effects are only efficient for the largest experimentally allowed primordial magnetic fields, we do not include them in our modelling.However, should the primordial magnetic field take such a high value, these effects would result in stronger heating of the IGM from cosmic rays, clustered around star-forming halos, and enhanced radio background, the former acting to diminish the 21-cm signal and the latter enhancing it.This sensitivity of the 21-cm signal, in particular the 21-cm power spectrum, to cosmic ray heating of the IGM provides an indirect probe of the primordial IGM magnetic field.
Overall we found that the distinct nature of cosmic ray heating means it leaves characteristic signatures in 21-cm tomographic maps, the high- end of the power spectrum, and in some cases the 21-cm global signal.These features are of sufficient magnitude to potentially be probed by the 21-cm signal measurements from current and next-generation experiments.Thus, an understanding of cosmic ray heating is necessary for the correct interpretation of the 21-cm signal, and conversely, we may be able to constrain early universe astrophysics that impacts cosmic rays (for example, the primordial IGM magnetic field) through the 21-cm signal.

Figure 2 .
Figure2.Comoving path length (panel a) of a cosmic ray proton at  = 6 emitted at  0 with kinetic energy  K,0 , and the fraction of its initial kinetic energy that is converted into the IGM heat (panel b).Here we have assumed the cosmic rays are moving in ionized IGM, with cosmic mean density and  e taken from the reference simulation (subsection 4.1.1).The black line indicates the  K,0 below which a cosmic ray is fully absorbed before  = 6.In panel a, the side length of one of our simulation cells is indicated as a green dashed contour for a reference, showing that cosmic ray protons emitted with energies below ∼ 2 MeV should not be able to escape their simulation cell of origin.Free-streaming cosmic rays emitted above ∼ 2 MeV but below ∼ 200 MeV can travel between cells but are still absorbed before  = 6, while protons with  K,0 > 200 MeV are never absorbed.As shown in panel b, low energy cosmic rays ( K,0 < 10 MeV) deposit a significant portion (> 10 per cent) of their initial kinetic energy as heat, with this fraction increasing with redshift due to the increase in  heat with  e .Conversely, at  K,0 > 1000 MeV less than 0.1 per cent of a cosmic ray kinetic energy becomes heat in the IGM, at least by  = 6.

Figure 3 .
Figure 3. Matched 21-cm global signals with different heating mechanisms.The efficiency of cosmic ray/X-ray emission was tuned for the models to minimize the root-mean-square error between their predicted global signals and that of the reference free-streaming cosmic ray heating model with  cr = 10 48 erg M −1 ⊙ .This procedure gave  cr = 1.39 × 10 48 erg M −1 ⊙ for the locallyconfined cosmic ray heating model,  X = 0.332 for the soft SED X-ray heating model, and  X = 60.0 for the hard SED X-ray heating model.All other astrophysical parameters listed in Table1and cosmological parameters, including the cosmological initial conditions, are kept the same between the simulations.The four models show good agreement at  > 22 as heating of the IGM has yet to become significant, as well as in the heating arm of the global signal around  ∼ 15.Differences are still seen between models at the global signal minimum and emission maximum.A vertical grey dashed line is shown at  = 16, illustrating the redshift at which we compare tomographic 21-cm maps, see Fig.4.

Figure 4 .
Figure 4. Slices of tomographic maps of IGM gas temperatures (top) and the 21-cm signals (bottom) for different cosmic ray (CR) and X-ray heating mechanisms shown at  = 16 where all four simulations predict approximately the same global signal ⟨ 21 ⟩ ≈ −84 mK.From left to right the IGM heating becomes more diffused as the heating carriers become longer-ranged (locally-confined CRs, free-streaming CRs, soft X-rays and hard X-rays).In addition, the mean IGM gas temperature decreases from 20.5 K to 16.0 K, 13.4 K, and finally 12.4 K.The change in heating results in a reduced contrast between emitting and absorbing regions in the 21-cm maps.All four simulations used the same cosmological initial conditions and the same slice of the simulation box is shown in each case.
Figure5.Comparison of the 21-cm power spectrum for the locally-confined cosmic rays, free-streaming cosmic rays, soft and hard X-rays (models are specified in subsection 4.1.1).Panel (a) shows the redshift evolution of the 21-cm power spectrum at a fixed wavenumber  = 0.1 cMpc −1 .Panel (b) shows the power spectra versus wavenumber  at the redshift when the corresponding global 21-cm signal equals to ⟨ 21 ⟩ = −84 mK (which occurs at  ∼ 16 for all the considered cases).We find that the small-scale (high ) 21-cm signal varies significantly between the models, suggesting that the variation of the 21-cm power spectrum with  could provide a diagnostic tool for the dominant heating mechanism.Also shown in both panels is the thermal-noise estimate for the SKA with 1000 hrs of observations (grey curve,Koopmans et al. 2015), illustrating the theoretical sensitivity of next-generation 21-cm power spectrum experiments.

Figure 7 .
Figure7.Evolution of the volume-averaged IGM temperature outside of ionized bubbles for different heating mechanisms.Shown are the kinetic temperature for models with inefficient, standard, and efficient cosmic ray (CR) heating, corresponding to  cr = 10 47 erg M −1 ⊙ , 10 48 erg M −1 ⊙ , and 10 49 erg M −1 ⊙ respectively; inefficient, standard, and efficient X-ray heating, corresponding to  X = 0.1, 1, and 10; as well as a model with both standard efficiency cosmic ray heating ( cr = 10 48 erg M −1 ⊙ ) and X-ray heating (  X = 1).We find that at standard efficiencies considered in the literature, cosmic ray heating heats up the neutral IGM more rapidly than X-rays.However, as both mechanisms at high redshifts are very uncertain, we advocate for the inclusion of both.

Figure 8 .
Figure 8. Variation of the IGM temperature and the 21-cm signal with astrophysical and cosmic ray parameters for diffusive escape cosmic rays.The kinetic temperature of the IGM (top row), the global 21-cm signal (middle row), and the 21-cm power spectrum at  = 0.1 cMpc −1 (bottom row) are shown for varying Pop II star formation efficiency  * ,II (left column), recovery time between Pop III and Pop II star formation  recov (middle column), and exponent of the cosmic ray spectrum  cr (right column). * ,II is found to have a strong impact on the rate of the IGM heating and thus the 21-cm signal.We find an increase in the star formation efficiency from 0.002 to 0.2 results in the IGM temperature at  = 8 raising from 39.7 to 1070 K. Conversely,  recov is found to only have a small impact.Larger values of  recov delay Pop II star formation and thus the IGM heating, leading to Δ ∼ 1 shifts in the 21-cm global signal and power spectrum across the range of  recov considered.The depicted variations with  cr replicate the findings ofLeite et al. (2017) that the cosmic ray heating rate is quite sensitive to the spectral shape of the cosmic rays injected into the IGM, with increasing  cr from −2 to −1.7 reducing the  = 8 IGM temperature from 284 to 19 K, and annihilating the power spectrum heating peak.

Figure 9 .
Figure 9. Variation of the IGM temperature and the 21-cm signal with astrophysical and cosmic ray parameters for the direct injection cosmic rays model.The kinetic temperature of the IGM (top row), the global 21-cm signal (middle row), and the 21-cm power spectrum at  = 0.1 cMpc −1 (bottom row) are shown for varying Pop III star formation efficiency  * ,III (left column), LW feedback delay parameter  LW (middle column), and maximum halo mass cosmic rays can escape from  max cr (right column).Some of the models are found to display a complex thermal history, with a maximum and a secondary minimum in  K .The height and redshift of this  K maximum are found to be quite sensitive to  LW and  max cr .This is due to the maximum being caused by  mol,crit (the critical halo mass for molecular cooling star formation) being increased above  max cr by the LW feedback.Due to the LW feedback cutting off cosmic ray heating, we do not find it to be as efficient as was found in the previous study(Sazonov & Sunyaev 2015).These complex IGM thermal histories in turn produce unusual inflexion points in the heating arm of the 21-cm signal (marked with black crosses), which we explore further in Fig.10.

Figure 10 .
Figure10.Flattened global 21-cm signal minimum from the cosmic ray heating cut off.Variation of the IGM kinetic temperature and the global 21-cm signal is shown colour-coded with respect to the maximum halo mass that cosmic rays can escape from.At the highest end of the range considered,  max cr = 10 7 M ⊙ ,  K develops a local peak due to the contribution of cosmic rays to heating.In turn, an inflexion point appears in the global signal.As  max cr is decreased the  K maximum merges with the first minimum, and the inflexion point in the global signal coincides with the minimum of the absorption trough producing a flattened absorption feature at  max cr = 4.57 × 10 6 M ⊙ (highlighted in black).For lower values of the critical mass, the global signal minimum gets narrower and deeper and is no longer flattened as cosmic ray heating becomes insignificant.

Table 1 .
Default values of 21-cm signal simulation parameters used throughout this work.If a simulation uses parameter values other than those listed here those that differ and their values will be explicitly stated.Ultimately the simulations provide full tomographic maps of  21 at each redshift step from which other 21-cm observables can be computed.
Muñoz et al. (2022)n the mass threshold due to the LW feedback and baryon-dark matter relative velocities can be each accounted for by factors  LW and  v cb which combine multiplicatively.By comparison to numerical simulationsMuñoz et al. (2022)find a best-fit model of Kulkarni et al. (2021)23) mass required for star formation via molecular cooling, this formula included modelling of the increase of this mass threshold due to the Lyman-Werner feedback and baryon-dark matter relative velocities.Muñoz et al. (2022),Nebrin et al. (2023), andHegde & Furlanetto (2023)recently highlighted theFialkov et al. (2012)formula does not account for the mitigation of the LW feedback by H 2 self-shielding, and consequently over-predicts the strength of the LW feedback compared to what is seen in the numerical simulations ofSchauer et al. (2021)andKulkarni et al. (2021).Hence, we adopt the fitting formula ofMuñoz et al. (2022)which has been calibrated to the aforementioned recent numerical simulations.Muñoz et al. (2022), motivated byKulkarni et al. (2021), assume