The redshifted 21cm signal from the Cosmic Dawn is expected to provide unprecedented insights into early Universe astrophysics and cosmology. Here we explore how dark matter can heat the intergalactic medium before the first galaxies leaving a distinctive imprint in the 21cm power spectrum. We provide the first dedicated Fisher matrix forecasts on the sensitivity of the Hydrogen Epoch of Reionization Array HERA telescope to dark matter decays. We show that with 1000 hours of observation HERA has the potential to improve current cosmological constraints on the dark matter decay lifetime by up to three orders of magnitude. Even in extreme scenarios with strong Xray emission from earlyforming metalfree galaxies the bounds on the decay lifetime would be improved by up to two orders of magnitude. Overall HERA shall improve on existing limits for dark matter masses below 2 GeVc SUP2SUP for decays into eSUPSUPeSUPSUP and below few MeVc SUP2SUP for decays into photons.
['cosmology of theories beyond the SM', 'dark matter theory', 'particle physics - cosmology connection', 'physics of the early universe', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies', 'High Energy Physics - Phenomenology']
21cm signal sensitivity to dark matter decay
{'No Header': 'Prepared for submission to JCAP  \nULB-TH/23-09', 'G. Facchinetti, a L. Lopez-Honorez, a,b Y. Qin, c,d and A. Mesinger e': "a Service de Physique Th'eorique, C.P. 225, Universit'e Libre de Bruxelles, Boulevard du Triomphe, B-1050 Brussels, Belgium  \nb Theoretische Natuurkunde & The International Solvay Institutes,  \nVrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium  \nc School of Physics, University of Melbourne, Parkville, VIC 3010, Australia  \nd ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)  \ne Scuola Normale Superiore, 56126 Pisa, PI, Italy  \nE-mail: gaetan.facchinetti@ulb.be, laura.lopez.honorez@ulb.be, yuxiang.l.qin@gmail.com, andrei.mesinger@sns.it  \nAbstract. The redshifted 21cm signal from the Cosmic Dawn is expected to provide unprecedented insights into early Universe astrophysics and cosmology. Here we explore how dark matter can heat the intergalactic medium before the first galaxies, leaving a distinctive imprint in the 21cm power spectrum. We provide the first dedicated Fisher matrix forecasts on the sensitivity of the Hydrogen Epoch of Reionization Array (HERA) telescope to dark matter decays. We show that with 1000 hours of observation, HERA has the potential to improve current cosmological constraints on the dark matter decay lifetime by up to three orders of magnitude. Even in extreme scenarios with strong X-ray emission from early-forming, metal-free galaxies, the bounds on the decay lifetime would be improved by up to two orders of magnitude. Overall, HERA shall improve on existing limits for dark matter masses below 2 GeV /c 2 for decays into e + e -and below few MeV /c 2 for decays into photons.  \nKeywords: dark matter theory, particle physics - cosmology connection, physics of the early universe, cosmology of theories beyond the SM  \nArXiv ePrint:  \n2308.16656", '1 Introduction': 'Dark matter (DM) is a crucial component of our Universe that shapes its evolution. Cosmological probes are amongst the powerful tools at our disposal to shed light on the nature of DM. In particular, the Cosmic Microwave Background (CMB) [1] currently provides the strongest constraints on DM properties. From CMB data, the DM relic density is determined at the percent level. The latter surpasses the baryonic matter density by more than a factor of five. Cosmological observations can also efficiently probe weak couplings between DM and standard model (SM) particles that would leave specific imprints. In that regards, they are complementary to astro-particle experiments [1, 2].  \nThe CMB is especially useful at constraining DM annihilation or decay into SM particles, see e.g., [3-14]. Both scenarios result in an exotic injection of energy into the intergalactic  \nmedium (IGM). For most of the SM final states, this exotic injected energy partially ionizes hydrogen and helium, increasing the residual free electron fraction post recombination. This has an observational imprint in the temperature and polarization anisotropy power spectra of the CMB 1 , constraining DM decay lifetimes up to τ ∼ 10 -24 s [16] and dark matter annihilation efficiencies up to p ann ∼ 3 × 10 28 cm 3 /s/GeV [1]. 2  \nLate time probes (sensitive at z ≪ 1000), including the Lymanα forest in Quasi-Stellar Objects (QSO) spectra and the 21cm signal of the hyperfine HI transition, are expected to be particularly efficient in testing late time energy injection. They are more sensitive to the IGM temperature T k and, consequently, to exotic energy injections transferred to the IGM in the form of heating [12, 14, 17-28]. In particular, the authors of Ref. [2] have used the Lymanα forest sensitivity to T k at redshifts z ∼ 4 -6 to derive constraints on decaying DM, see also [2, 10, 20, 29]. They have improved on the CMB bounds for decaying DM lifetimes, disfavoring to τ ∼ 10 -25 s for DM decays into electron-positron pairs for m DM < MeV /c 2 . On the other hand, the 21cm signal will be sensitive to T k at even earlier times, during the so-called Cosmic Dawn (CD) of galaxies ( z ∼ 10-20) and the epoch of reionisation (EoR, z ∼ 5-10). Interestingly, the relative dearth of galaxies during the CD should make it easier to isolate an additional heating contribution from DM decay or annihilation.  \nHere we revisit the imprint of DM decays on the cosmological 21cm signal and obtain the first forecasts of DM lifetime constraints that will be enabled by 21cm power spectrum measurements. In particular, we focus on the Hydrogen Epoch of Reionisation Array (HERA) telescope which was designed to measure the 21cm power spectrum at a high signal to noise ratio (S/N). HERA has completed deployment [30] and is currently analysing data from an extended observational campaign. An initial observational result performed with 71 antennas (out of the total 331) and only 94 nights of measurement has already provided the most constraining upper bounds on the 21 cm power spectrum at redshifts z = 8 and 10 [31]. Combined with complementary observations of galaxy UV luminosity functions and the timing of the EoR, these upper limits imply significant early IGM heating [31]. If this heating was provided by high mass X-ray binary stars in CD galaxies, expected to dominate the X-ray background at high redshifts (e.g., [32]), they would need to be considerably brighter than those observed today (perhaps because they were born in extremely metal poor environments; see e.g., [33]). This clearly illustrates the potential of 21-cm measurements to constrain IGM heating during the Cosmic Dawn. Consequently, it also indicates that 21cm cosmology will soon be mature enough to quantitatively probe exotic heating scenarios like decaying DM models.  \nOur numerical approach is built on the 21cmFAST [34, 35] and the DarkHistory packages [36, 37], which we combine into a hybrid code, exo21cmFAST 3 . We organize this paper as follows. In section 2, we briefly review 21cm cosmology and how to model and simulate the 21cm signal with exo21cmFAST . In section 3, we show the imprint of DM decay on the IGM and the corresponding 21cm signal. In section 4, we present Fisher forecasts of constraints available  \nwith a 1000h observation using HERA. We forecast joint constraints on decaying DM properties and two different models of CD galaxies. Finally, in section 5, we discuss these results and conclude. In this work, we use a ΛCDM cosmology with parameters (Ω m , Ω b , Ω Λ , h, σ 8 , n s ) = (0.31, 0.049, 0.69, 0.68, 0.81, 0.97) following the TT,TE,EE+lowE+lensing+BAO result of Planck 2018 [38]).', '2 A glimpse into 21cm cosmology': 'We start with a brief review of 21-cm cosmology in section 2.1 before outlining in section 2.2 our model for the astrophysics of first galaxies, which we then expand to include imprints of DM decays.', '2.1 The 21cm signal and its power spectrum': "The redshifted cosmic 21cm signal, arising from the hyperfine spin-flip transition of neutral hydrogen, can be seen in emission or in absorption compared to the radio background. Here we fix the latter to the CMB whose temperature is denoted by T CMB ( z ). The differential brightness temperature of the 21cm signal can be expressed as [39]  \nδT b ≈ 20mK ( 1 -T CMB T S ) x HI (1 + δ b ) ( 1 + 1 H d v r d r ) -1 √ 1 + z 10 0 . 15 Ω m h 2 Ω b h 2 0 . 023 , (2.1)  \nwhere T S is the spin temperature of neutral hydrogen in the IGM, x HI is the neutral fraction, δ b is the relative density perturbation in the baryon number density, Ω m , b are the matter and baryon energy densities relative to the critical energy density today and h is the Hubble parameter in units of 100 km/s/Mpc. The relative motion of the neutral gas with regard to the Hubble flow, with H denoting the Hubble rate, is taken into account with the d v r d r term that is the gradient of the proper velocity along the line of sight. 4 The spin temperature, T S , quantifies the relative occupancy of the two hyperfine levels of the ground state of neutral hydrogen. It is obtained from the equilibrium balance of (i) absorption/emission of 21cm photons from/to the CMB background at a temperature T CMB ; (ii) collisions with atoms and electrons in the IGM, with the gas kinetic temperature denoted by T k ; and (iii) resonant scattering of Lymanα photons coupling the spin temperature to T k . When the spin temperature is coupled to IGM gas kinetic temperature through (ii) or (iii), T S can differ from T CMB . The signal can appear in absorption (if T S < T CMB ), in emission (if T S > T CMB ) or can be zero (if T S = T CMB and/or x HI = 0).  \nSpatial variation of IGM properties leads to fluctuations in the 21cm signal. In what follows, we refer to the 21cm global signal, δT b , as the sky averaged brightness temperature while the 21cm power spectrum refers to the dimensional quantity, δT 2 b ∆ 2 21 , obtained from:  \nδT 2 b ∆ 2 21 ( k, z ) = δT 2 b ( z ) × k 3 2 π 2 P 21 ( k, z ) (2.2)  \nwhere P 21 is defined as:  \n⟨ ˜ δ 21 ( k , z ) ˜ δ 21 ( k ' , z ) ⟩ = (2 π ) 3 δ D ( k -k ' ) P 21 ( k, z ) (2.3)  \nwith ⟨⟩ the ensemble average, k the comoving wave vector, and ˜ δ 21 ( k , z ) the Fourier transform of δ 21 ( x , z ) = δT b ( x , z ) /δT b ( z ) -1. Notice that x denotes the position vector.", '2.2 Modelling the 21cm signal': 'As seen above, the 21cm signal depends on the IGM temperature and ionization fraction. These quantities, in turn, are influenced by the radiation emitted by stars and sources of exotic energy in the late Universe. Here we review how astrophysical sources are modeled in 21cmFAST , before moving onto the inclusion of exotic heating in our modified version, exo21cmFAST , in the following section.', '2.2.1 Ionization, excitation and heating of the IGM': 'The EoR is an inhomogeneous process, with ionizing photons from galaxies carving out cosmic HII regions that grow and eventually overlap. We denote the volume filling factor of HII regions as Q HII ( z ) (thus filling factor of the mostly neutral regions is 1 -Q HII ). The 21cmFAST codes employ an excursion-set algorithm 5 to compute inhomogeneous reionization [44]. A cell is considered as ionized when  \nn ion ≥ (1 + n rec )(1 -x e ) , (2.4)  \nwhere n ion denotes the cumulative number of ionizing photons per baryon (see section 2.2.2), n rec gives the cumulative number of recombinations per baryon within spheres of decreasing radii and x e accounts for secondary ionization in mostly-neutral IGM. The redshift evolution of latter quantity and of the gas kinetic temperature ( T k ) is described by  \ndx e dz = dt dz ( Λ X ion +Λ DM ion -α A C x 2 e n b f H ) , (2.5)  \ndT k dz = 2 3 k B (1 + x e ) dt dz ∑ β ϵ β heat + 2 T k 3 n b dn b dz -T k 1 + x e dx e dz , (2.6)  \nwhere k B is the Boltzmann constant, ϵ β heat ( x , z ) is the heating rate per baryon from different sources β including Compton scattering (effective at high redshifts, z ≳ 300) as well as heating by X-rays ( ϵ X heat ) and DM ( ϵ DM heat ). Λ X , DM ion are the secondary ionization rates per baryon accounting for the contribution from either X-rays or DM, α A is the case-A recombination coefficient, C ≡ ⟨ n 2 e ⟩ / ⟨ n e ⟩ 2 is the clumping factor, set to its default value C = 2 in our analysis, with n e the electron number density, and f H = n H /n b is the hydrogen number fraction. The implementation of Λ DM ion and ϵ DM heat accounting for DM energy injection is specific to exo21cmFAST which is detailed in section 3. When considering ionization in both the mostly-ionized and mostly-neutral IGM, the total mean ionization fraction becomes:  \nx i ≈ Q HII + x e (1 -Q HII ) . (2.7)  \nThe strength of the Wouthuysen-Field effect that couples the spin temperature of gas to its kinetic temperature is determined by the total Ly α background, which has the following components:  \nJ α = J X α + J ⋆ α + J DM α , (2.8)  \nwhere J X α arises from X-ray sources, J ⋆ α encapsulates the contribution from stellar photons with energies between Lymanα and the Lyman limit, while J DM α accounts for the DM contribution introduced in section 3. We then evaluate the Lymanα coupling efficiency according  \nto  \nx α = 1 . 7 × 10 11 1 + z ( J α s -1 Hz -1 cm -2 sr -1 ) S α , (2.9)  \nwhere S α acts as a quantum mechanical correction of the order of unity [45] and further compute the spin temperature as  \nT -1 S = T -1 CMB +( x α + x c ) T -1 k 1 + x α + x c , (2.10)  \nwith x c denoting the spatially-varying collisional coupling efficiency in the IGM.  \nIn the next subsections, we discuss our parameterization of galaxy properties, which allows us to determine the budget of photons in various wavelength ranges and evaluate the aforementioned radiative backgrounds that are essential to the 21cm signal.', '2.2.2 Star formation and galaxy evolution': "As we wish to consider the 21cm signal during the initial stages of the CD, our models will include the very first galaxies, hosted by so-called minihalos (with a virial mass M vir ≲ 10 8 M ⊙ ) for which the dominant cooling channel is provided by roto-vibrational transitions of H 2 . Due to their metal free gas and shallower potential wells, these molecular-cooling galaxies (MCGs) might have different properties compared to the atomic-cooling galaxies (ACGs) that we have observed at later times. They are expected to predominately host Population-III stars 6 and likely possess different initial mass functions (IMFs), as well as potentially different star formation efficiencies (see e.g., [46-49] and references therein).  \nThe star formation rate (SFR) of an ACG/MCG inside a halo of virial mass M vir is parameterised by [43]  \n˙ M II / III ⋆ ≡ f II / III ⋆ M vir t ⋆ H ( z ) -1 (2.11)  \nwhere f II ⋆ and f III ⋆ represent different stellar-to-halo mass ratios for ACGs and MCGs, respectively. Here, t ⋆ H -1 is a characteristic star-formation time-scale, with t ⋆ a dimensionless parameter taking values between zero and one (for our analysis we chose a fiducial value t ⋆ = 0 . 5, see table 1 for all fiducial values). Through t ⋆ , the model encompasses scenarios where all stars are formed in an instantaneous burst event or as a gradual buildup over the age of the Universe. References [50, 51] have further inferred a nearly non-evolving, powerlaw relation between f II / III ⋆ and M vir for galaxies that are expected to dominate reionization. Therefore, we consider  \nf II ⋆ = Ω b Ω m min { 1 , f II ⋆, 10 ( M vir 10 10 M ⊙ ) α II ⋆ } f III ⋆ = Ω b Ω m min { 1 , f III ⋆, 7 ( M vir 10 7 M ⊙ ) α III ⋆ } , (2.12)  \nwhere f II ⋆, 10 , f III ⋆, 7 , α II ⋆ and α III ⋆ are free parameters. Similarly, the UV ionizing escape fractions of ACGs and MCGs are also parameterized as power-law relations with the halo mass  \nf II esc = min { 1 , f II esc , 10 ( M vir 10 10 M ⊙ ) α esc } f III esc = min { 1 , f III esc , 7 ( M vir 10 7 M ⊙ ) α esc } (2.13)  \nwhere f II esc , 10 , f III esc , 7 and α esc are additional free parameters (for simplicity, the same scaling index is considered in this work for f II / III esc ). The efficiency parameters, f II ⋆, 10 , f III ⋆, 7 , f II esc , 10 , f III esc , 7 can vary between zero and one. 7 In our Fisher Matrix analysis, we set log 10 f II ⋆, 10 = -1 . 3 and log 10 f II esc , 10 = -1 . 0 while we take log 10 f III ⋆, 7 = -2 . 0 and log 10 f III esc , 7 = -2 . 0 following [52] and [43]. On the other hand, values of the scaling index parameters are mostly guided by observational results or detailed simulations. For instance, fitting their modelled stellar mass function to observations, the authors of Ref. [52, 53] find that α II ⋆ is close to 0.5. This result is consistent with a supernova-regulated galaxy growth history [54]. In addition, α esc is likely negative as we expect more low-column density channels to be created by supernovae in low-mass galaxies and to allow UV ionizing photons to escape into the IGM. In our Fisher Matrix analysis, we set α II ⋆ = α III ⋆ = 0 . 5 and α esc = -0 . 5.  \nThe number density of ACGs and MCGs per unit halo mass can be described by the halo mass function (HMF) 8 (d n/ d M vir ) weighted by different duty cycles ( f II / III duty ). The duty cycles account for inefficient star formation in low mass ACGs and MCGs. This mass threshold for ACGs is given by either the atomic cooling limit ( M atom ) or by a characteristic mass ( M RE crit ) below which photoheating during reionization is able to significantly suppress gas content [56]. The duty cycle of ACGs can be described by [43, 52]  \nf II duty ≡ exp ( -M II turn M vir ) . (2.14)  \nWhen considering both ACGs and MCGs in our analysis, we self-consistently evaluate M II turn as max ( M atom , M RE crit ) where M atom = 5 × 10 7 M ⊙ ((1+ z ) / 10) -1 . 5 [57]. Furthermore, the duty cycle of MCGs can be written as [43, 58]  \nf III duty ≡ exp ( -M III turn M vir ) exp ( -M vir M atom ) . (2.15)  \nCompared to eq. (2.14), the additional factor of exp ( -M vir /M atom ) results from the assumption that there is a smooth transition from MCGs (below M atom ) hosting PopIII stars to ACGs (above that mass) hosting PopII stars. The lower mass cutoff M III turn , for MCGs star formation, is set by three types of feedback mechanisms. Similarly, MCGs' gas reservoir also becomes inadequate when their host halos are too small (i.e., < M RE crit ) to withstand the UV ionizing background. In addition to this, as molecular hydrogen can be easily dissociated by photons within the Lyman-Werner (LW) energy range (11 . 2 -13 . 6 eV), star formation is also quenched in regions that have a LW background above a critical value M LW crit . In our  \nanalysis, we have thus M III turn = max ( M RE crit , M LW crit ) and evaluate these critical masses using the excursion algorithm [56]. 9  \nNote that in this work, we first consider a simplified model in which ACGs and MCGs follow the same scaling relations. In this case, it is effectively a single-population model and therefore f II ⋆, 10 , α II ⋆ , f II esc , 10 , α esc are sufficient to describe their stellar mass/UV ionizing escape fractions. In addition, the threshold for star forming is set by f II duty as in eq. (2.14) with M II turn = M turn taken as a free parameter for more flexible interpretation of the feedback strength. In models without MCGs, M turn is set above 10 8 M ⊙ due to the atomic cooling limit at these redshifts. In addition, currently observed galaxy UV luminosity functions [6062] indicate M turn is likely to be lower than 3 × 10 9 M ⊙ [52]. In our Fisher Matrix analysis, we will use the fiducial value of M turn = 10 8 . 7 M ⊙ .  \nWith these definitions, we can evaluate the UV ionizing photon budget from eq. (2.4):  \nn ion = 1 ρ b ∑ i ∈{ II , III } ∫ d M vir d n d M vir f i duty M vir f i ⋆ f i esc n i γ , (2.16)  \nwhere ρ b is the baryon density, n II / III γ refers to the number of UV ionizing photons emitted per stellar baryon and, in this work, is chosen to be 5000 and 50000 for ACGs and MCGs, respectively. We can further estimate photoheating feedback and evaluate M RE crit following [56]. Similar approaches are also adopted to calculate the Lymanα background from stellar components (i.e., J ∗ α in equation 2.8) as well as the LW radiation and its feedback strength ( M LW crit ). We refer readers who are interested in these details to Ref. [43].", '2.3 X-ray heating and its impact on the 21cm signal': 'To estimate the heating due to X-rays as well as their contribution to ionizing the mostlyneutral IGM and coupling its spin temperature to the kinetic temperature, we further assume that the specific X-ray luminosity per unit star formation rate (denoted as L II / III X in units of s -1 per M ⊙ yr -1 ) follows power-law relations [63] with an energy index of α X , i.e., L II / III X ∝ E -α X (in this work α X = 1). It is normalized with the integrated soft-band ( E < 2 keV) luminosity per SFR (in units of erg yr s -1 M -1 ⊙ ), that we denote L II / III X :  \nL II / III X ≡ ∫ 2keV E 0 d E L II / III X , (2.17)  \nwhere E 0 is the X-ray energy threshold above which photons can not escape the host galaxy. L II X and L III X are considered as free parameters in our analysis while L III X is set equal to L II X in the single-population model. These allow us to evaluate the X-ray emissivity at the emitting redshift:  \nϵ X = ∑ i ∈{ II , III } ∫ d M vir d n d M vir f i duty ˙ M i ⋆ L i X . (2.18)  \nFigure 1 : 21cm signal for 3 illustrative scenarios. The fiducial model considered in this analysis is shown with a black line (with log 10 L II X = 40) and 2 other scenarios with higher values of the X-ray luminosity amplitude are shown with gray lines (dark gray and light gray: log 10 L II X = 40 . 5 and 41 respectively). The blue region around the fiducial model corresponds to the 2σ thermal noise for a 1000h observation with HERA, while the green region additionally includes cosmic variance and modeling uncertainty (see section 4.1 for details). Left panel: 21cm power spectrum as a function of the redshift for a fixed scale, k = 0 . 13/Mpc, relatively free of foregrounds (top), and the sky averaged differential brightness temperature (bottom). Right panel : 21cm power spectrum as function of the scale at redshift z = 7 . 5 (top) and z = 10 . 9 (bottom).  \n<!-- image -->  \nNote that, the X-ray threshold E 0 depends on the interstellar medium as well as the environment of X-ray sources (likely to be high mass X-ray binaries). Motivated by hydrodynamic simulations of the first galaxies [63], E 0 is assumed to range between 0.2 and 1.5 keV, which correspond to a column density of roughly 2 × 10 19 and 10 23 cm -2 , respectively. In our Fisher Matrix analysis, we take E 0 = 0 . 5 keV as a fiducial value. On the other hand, L II / III X for high-redshift galaxies is less understood. Local, star-forming galaxies typically have L X ∼ 10 39 . 5 [64], though CD galaxies due to their lower metallicity ISM are expected to have higher values of L X (e.g., [65-67]). We set log 10 L II / III X = 40 as the fiducial model in this work, motivated by recent analyses of the HERA observation [31, 68]. The consequences of a larger value log 10 L II X = 41 are discussed in appendix C.  \nFor illustration, figure 1 shows the 21cm signal for three values of L II X considering a single population of galaxies. The black line corresponds to a fiducial value of log 10 L II X = 40. The blue region around the fiducial model corresponds to the 2 σ thermal noise for a 1000h observation with HERA, while the green region additionally includes cosmic variance and modeling uncertainty (see section 4.1 for details). In the left panel, we show the power spectrum (top) at a fixed scale k = 0 . 13 Mpc -1 (that is expected to be relatively free from foreground contamination) and the global signal (bottom) as a function of redshift.  \nUsing the fiducial model as benchmark, we see that the global signal (bottom plot) displays an absorption trough between z ∼ 16 and 9, and in emission below z ∼ 8. Very roughly,  \nthe minimum of absorption separates the so-called epoch of Lymanα or Wouthuysen-Field (WF) coupling ( z ≳ 11) and the subsequent epoch of heating (EoH, for z ≲ 11). Between z ∼ 8 and z ∼ 6, the ionized fraction of IGM increases, i.e., x HI decreases, correspoding to the EoR. Around z ∼ 6 the IGM is fully ionized and δT b becomes zero. In parallel, we see that the power spectrum on these large scales, i.e., k ≲ 0 . 13 / Mpc, (upper plot) displays multiple peaks as a function of redshift. As discussed in previous works (e.g., [69-71]), three peaks in the evolution of the large-scale power correspond to the three astrophysical epochs: WF coupling, EoH and EoR, when the 21cm power spectrum is dominated by the auto-power spectrum of x α , T K and Q HII , respectively. Between these epochs, the corresponding cross power spectra are negative and have a significant contribution, causing the large scale power to drop. 10 For example, the first regions to ionize are those in the proximity of galaxies that are also exposed to a stronger X-ray background and are hotter than average. Thus, reionization effectively zeros the peaks of the 21cm signal during the EoH. This causes the large scale power to decrease, before increasing again when fluctuations in Q HII begin to dominate during the middle stages of the EoR.  \nIncreasing log 10 L II X shifts the EoH to earlier times. This decreases the overlap between the EoR and EoH, diminishing the importance of the cross terms and thus increasing the peak power spectrum amplitude during the EoR. However, it increases the overlap between the EoH and WF coupling, boosting cross terms and decreasing the associated peak power during the EoH. Similar effects (increase of the EoR peak and suppression of power in EoH) can be observed when considering two populations of galaxies and increasing L III X , the normalisation of X-ray flux from MCGs. This is illustrated in appendix A.2. 11', '3 The impact of Dark matter decay on the 21cm signal': 'Dark matter annihilations and decays inject energy into the IGM. This induces extra heating, ionization and excitation of the IGM that leave a footprint in the 21cm signal. In this paper, we focus on dark matter decays for two reasons. First, there is a large room for improvement on current constraints as DM decay mainly impacts cosmological observables at relatively late times (at z ≪ 1000) while DM annihilation is already be strongly constrained by CMB observations (see e.g., the discussion in [1, 10, 20, 38]). In addition, the DM annihilation signature in the 21cm signal is expected to strongly depend on the late time boost arising from structure formation (see e.g., [18, 19]), contrary to CMB bounds [9, 11]. We therefore leave the study of DM annihilation for future work.', '3.1 Dark matter energy injection and deposition': 'A priori, DM can decay into a plethora of SM final states. Yet, after subsequent decays, hadronization processes, etc., it is only the leftover photons, electrons, and positrons that will efficiently deposit their energy into the IGM (see e.g., [6, 8, 74, 75]). For this reason, we will focus on decays into electron-positron pairs and photon pairs.  \nHigh energy SM final state particles injected into the IGM from DM decays do not deposit their energy instantaneously. In addition, multiple channels (denoted with a subscript a ) of energy deposition have to be considered. These channels are IGM heating (with a = heat), Hydrogen ionization ( a = HII), Helium single or double ionization ( a = HeII or HeIII), and neutral atom excitation ( a = exc). The energy deposition rate (per average baryon number) can be expressed in terms of the energy injection rate as [7]  \nϵ DM a ≡ 1 n b ( z ) ( dE a ( x e , z ) dt dV ) deposited = f a ( x e , z ) 1 n b ( z ) ( dE ( z ) dt dV ) injected , (3.1)  \nwhere n b ( z ) denotes the baryon number density. If all of the dark matter decays with a lifetime, τ , longer than the age of the universe, the injected energy per unit time and volume takes the form:  \n( dE dV dt ) injected = (1 + z ) 3 ρ DM , 0 c 2 τ , (3.2)  \nwhere c denotes the speed of light. Below we also use the DM decay rate, Γ ≡ 1 /τ , to characterize DM energy injection.  \nThe coefficients, f a ( x e , z ), in eq. (3.1) are the DM energy deposition efficiencies. They account for all of the details associated with the delay in energy deposition and separation into different channels a at a given redshift z and free-electron fraction x e . In this analysis, we compute energy deposition efficiencies using the public DarkHistory package [36]. 12 We have implemented, in exo21cmFAST , a slightly adapted version of release 1.1 (with tables of the deposition fractions that have been upgraded from release 1.0) [37]. To account for early exotic energy injection into the IGM, DarkHistory not only evaluates the deposition efficiencies but also solves for the homogeneous IGM evolution between redshift 3000 and 35. Afterwards, 21cmFAST takes over to simulate the full IGM evolution including late-time astrophysical processes from cosmic dawn to reionisation.  \nThe energy deposition rate per channel per baryon ϵ DM a , defined in eq. (3.1), enters directly in the evolution of the ionized fraction, the IGM temperature and the Lymanα flux [17-19], which drive the 21cm signal (see section 2). Indeed, the energy injection from DM decays contributes to eq. (2.5), eq. (2.6), and eq. (2.8) via the heating rate ϵ DM heat , the ionization rate Λ DM ion , and the Ly α flux J DM α . While ϵ DM heat is directly computed from eq. (3.1) by DarkHistory , the DM-induced ionization rate and Ly α background flux are respectively derived as  \nΛ DM ion = f H ϵ DM HII E HI th + f He ϵ DM HeII E HeI th , (3.3)  \nJ DM α = c n b 4 πH ( z ) ν α ϵ DM Ly α hν α . (3.4)  \nDarkHistory accounts for the variation of the deposition efficiencies with the ionization fraction, x e . The DM source terms in eq. (2.5) and eq. (2.6) are thus not constants but depend on the state of the IGM. The authors of [2] refer to this as the backreaction of the IGM on the deposition efficiencies. In DarkHistory we can switch the backreaction off to  \nFigure 2 : Heating rates per baryon from PopII-dominated ACGs (blue), PopIII-dominated MCGs (red) and DM (green). For this figure, we have considered either log 10 ( L II X ) = log 10 ( L III X ) = 40 (left panel) or log 10 ( L II X ) = log 10 ( L III X ) = 41 (right panel), with the remaining parameters fixed to their fiducial values given in table 1. The DM contribution is computed assuming a 100 MeV /c 2 DM particle species decaying into electron positron pairs with a lifetime τ = 10 27 s (left panel) and τ = 10 26 s (right panel).  \n<!-- image -->  \ngauge its significance by considering the IGM evolution that would result from a scenario without exotic energy injection, ˜ x e and forcing f a ( x e = ˜ x e , z ). For the models considered in our Fisher Matrix analysis we have noticed that the backreaction has a negligible impact on the ionization and temperature evolution.  \nAs heating will play an important role in this analysis, we illustrate different contributions to the heating rate per baryon, ϵ β heat in eq. (2.6) in figure 2. The blue and red curves show the X-ray heating rate from ACGs and MCGs while the green curve is due to heating from DM decays discussed in section 3. We assume that DM particles have mass m χ = 100 MeV /c 2 and decay into electrons and positrons with a lifetime τ = 10 27 s in the left panel and τ = 10 26 s in the right panel. We further show the impact of the value of log 10 ( L II / III X ) fixing these parameters to 40 in the left panel or 41 in the right panel. The rest of the astrophysical parameters are set to the fiducial values listed in table 1. By definition, MCGs start to efficiently heat the medium before ACGs. However, due to rapid growth of the halo mass function, ACGs soon become dominant (see also, for example, [43, 58]). The DM contribution shows a flat dependence in redshift. Indeed, in the case of DM decays, the heating rate per baryons essentially scales as (1 + z ) 3 /n b ( z ), which is redshift independent, see eqs. (3.1) and (3.2). The prefactor f heat ( z, x e ) induces an extra redshift dependence. The latter is rather mild for a DM candidate decaying into electron-positron pairs between z = 5 and 25. Similar redshift dependence is expected for low mass DM decaying into a pair of photons. DM heating appears thus to be dominating at the earliest times.  \nThe imprint of DM decay on the 21cm power spectrum is strongest at early times, z ≳ 10 -15 for these examples.', '3.2 Imprint in the 21cm signal': 'As previously pointed out, dark matter energy injection mainly affects the 21cm signal as a new global heating source (e.g., refs. [17, 18]). For the fiducial model considered here,  \nFigure 3 : Same as figure 1 with the two gray lines illustrating the case of a 100 MeV /c 2 dark matter candidate decaying to electron positron pair with a lifetime of 10 28 s (dark gray) and 10 27 . 5 s (light gray).  \n<!-- image -->  \nk  \n[Mpc  \n-  \n1  \n]  \nthe ionizing photons can not compete with the ones from astrophysical source. 13 Figure 3 illustrates the effect of this new source of heating on the global signal and the power spectrum when considering a 100 MeV /c 2 dark matter candidate decaying at 100% into a electron positron pair with a lifetime of 10 27 . 5 (light gray) and 10 28 s (dark gray). As visible in eq. (3.2) a shorter lifetime implies a stronger heating of the IGM and hence induces a shallower absorption trough in the global signal.  \nNotice that figure 3 illustrates the case of a single population of galaxies (AGCs only). This is the scenario in which the DM heating can be more easily disentangled from the astrophysics sources as AGCs heating rate becomes comparable to the DM one at rather late times ( z ≲ 10 -15), see figure 2. In addition this scenario is better constrained by data with a minimal number of sources of heating and the least number of astrophysics parameters that could be degenerate with DM heating, see section 4.2. The case of a more complete, yet more complex, astrophysics model is illustrated in appendix A.2 where we consider both ACGs and MCGs sources of heating. The latter scenario involves more sources of uncertainties as MGCs properties are yet to be determined. We will also consider this scenario in our analysis. MGCs can heat the medium earlier than AGCs and the DM imprint becomes less easy to untangle from astrophysics. This is expected to mitigate the constraints on exotic sources of heating, see section 4.2 for a quantitative result.  \nIn the model illustrated in figure 3, DM energy injection dominates IGM heating at early times when z > 10-15. As a result, DM decays give rise to a more uniformly heated IGM at early times, which can decrease the large-scale 21cm power during the EoH compared with galaxy-only heating (see also [18, 19]. On the other hand, the impact of DM  \n]  \n2  \n[mK  \n221  \n∆  \n2  \nb  \nT  \nδ  \n]  \n2  \n[mK  \n221  \n∆  \n2  \nb  \nT  \nδ  \n10  \n2  \n10  \n1  \n10  \n10  \n3  \n2  \n10  \n1  \n0  \n.  \n1  \n0  \n.  \n2  \n0  \n.  \n3  \n0  \n.  \n4  \n0  \n.  \n5  \n0  \n.  \n6  \n0  \n.  \n7  \n0  \n.  \n8  \n0  \n.  \n9  \nz  \n= 7  \n.  \n5  \nχ  \nχ  \nz  \n= 10  \n.  \n9  \n→  \n→  \ne  \ne  \n+  \n+  \ne  \n-  \n,  \nm  \ne  \n-  \n,  \nm  \nχ  \nχ  \n= 100 MeV  \n/c  \n= 100 MeV  \n/c  \n2  \n2  \nτ  \n=  \nτ  \n= 10  \nτ  \n= 10  \n∞  \n(fid  \n.  \n)  \n28  \n.  \n0  \n27  \n.  \n5  \ns  \ns  \nFigure 4 : Lightcone slices of the differential brightness temperature in our (250 Mpc) 3 large simulation box. We show results for the fiducial model with log 10 ( L II X ) = 40 and τ = ∞ (upper panel), log 10 ( L II X ) = 41 (middle panel), and for a 100 MeV /c 2 decaying DM with τ = 10 27 . 5 s (bottom panel).  \n<!-- image -->  \ndecay on the late-time 21cm power ( z < 9) is far more modest. As mentioned earlier, the contribution to reionisation from DM decay is sub-dominant to that from stars, even for even the shortest lifetimes considered here. Our shortest lifetime model (gray curve) does show a slight enhancement of 21-cm power in the early stages of the EoR ( z ∼ 7-8), when the signal is in emission. The additional heating from DM decay increases the (1 -T CMB /T S ) temperature term from eq. (2.1). However, this term quickly saturates to unity (due to the X-rays from stars) and so the different DM lifetimes considered here do not impact the 21-cm power below z ≲ 7.  \nThe discussion above is further illustrated with figure 4 that shows a 2D slice through the brightness temperature lightcone. We show the fiducial model with L II X = 10 40 and no decaying DM (i.e., τ = ∞ ) (top panel), a model with larger L II X (central panel) and a model with a shorter lifetime τ (bottom panel). The figure illustrates how the different features of the power spectrum (in Fourier space) relate to the real space signal - e.g., a lower power spectrum is related to a weaker signal and/or low contrast.  \nWe conclude that heating from DM decay has a qualitatively similar imprint on the 21cm signal as heating from galaxies, but their redshift and scale dependencies can be very different. Given that experiments such as HERA will be able to probe a large range of redshifts and scales (see section 2), we expect them to be able to disentangle these two sources of heating. From the sensitivity estimates shown in figure 3, we expect HERA could probe lifetimes up to 10 27 -10 28 s, surpassing the sensitivity from CMB and Lymanα probes. We quantify this further below.', '4 Method and results': "Our goal is to forecast the sensitivity of the full HERA array with 331 antennas to DM decay. For that purpose multiple approaches are possible. A cost-efficient solution in terms of computational resources is the Fisher Matrix. For a given set of cosmological parameters θ the Fisher information matrix is defined as  \nF ij ≡ -〈 ∂ ln L ∂θ i ∂θ j 〉 (4.1)  \nTable 1 : Values of the astrophysical parameters in the fiducial models. Parameters specific to PopII-dominated ACGs and PopIII-dominated MCGs are listed in the first and second row, respectively. The other parameters are in the last row. Note that M turn is only used in the simplified one-population model (see more details in section 2.2.2).  \nwhere L denotes the likelihood, the expected distribution of the data given a certain model. The Cr'amer-Rao theorem [76-78] states that the marginalized error, σ θ i , on a given parameter θ i follows σ θ i ≥ √ ( F -1 ) ii , implying the Fisher matrix approach always gives a local optimistic estimate of the errors. In the following, we estimate the expected variance from the diagonal elements of the inverse Fisher Matrix components,  \nσ 2 θ i = ( F -1 ) ii . (4.2)  \nMore details on the Fisher matrix method are given in appendix B.1.  \nIn the following, we describe our treatment of the Fisher Matrix forecasts within a new automatized python package called 21cmCAST 14 - based on 21cmfish [79]. We then discuss the results of our analysis, which demonstrates that 21cm cosmology could provide key constraints on dark matter energy injection.", '4.1 Fisher matrix analysis': 'We evaluate the 21cm power spectrum δT b 2 ∆ 2 21 on a fixed grid of modes k and redshifts z . We assume a total bin number of N k × N z and consider uncorrelated z and k bins as in [79]. Assuming that the likelihood L can be described by a Gaussian distribution for each bin, the Fisher matrix elements reduce to the sum  \nF ij = N k ∑ i k N z ∑ i z 1 σ 2 ∆ ( k i k , z i z ) ∂δT b 2 ( z i z | θ )∆ 2 21 ( k i k , z i z | θ ) ∂θ i ∂δT b 2 ( z i z | θ )∆ 2 21 ( k i k , z i z | θ ) ∂θ j . (4.3)  \nIn accordance with the common practice in the field [52, 79, 80], we account for three contributions to the measurement error in the power spectrum:  \nσ 2 ∆ ≡ σ 2 exp + σ 2 Poisson + [ 0 . 2 δT b 2 ∆ 2 21 ] 2 . (4.4)  \nThe first contribution, σ exp denotes the experimental error (i.e., thermal noise), which is evaluated using the public code 21cmSense [81, 82] 15 . In that code we fix the HERA experiment as an hexagonal array of 331 antennas (11 on each side) with separation and dish size of 14 m, located at a latitude of ∼ 30 . 8 · . Moreover, we fix the bandwidth to B = 8 MHz and  \nthe spectral resolution to δν ∼ 100 kHz. We assume a total operating time of 1000 hrs (6 hrs per day during 166.7 days). We adopt the moderate foregrounds setting, with the default system temperature for the 21cm line, implying a system temperature with the following redshift dependence: 16  \nT sys ( z ) = 100 K + 260 K ( ν ( z ) 150 MHz ) -2 . 6 ∼ 100 K + 300 K ( 1 + z 10 ) 2 . 6 , (4.5)  \nwith σ exp ( k i k , z i z ) ∝ T 2 sys ( z iz )[81]. For illustration, this thermal noise is shown in blue in Figs. 1 and 3. The second contribution to the measurement error in eq. (4.4) is the cosmic variance, σ Poisson , arising from the evaluation of the power spectrum in a finite size simulation box. Finally, as a third contribution to σ ∆ , we consider a modeling uncertainty of 20% [80]. For illustration, the total measurement error is shown in green in figures 1 and 3.  \nAs shown in eq. (4.5), the thermal noise increases with the redshift. We have set the maximum redshift considered in this analysis to z max ∼ 20, after having checked that higher redshifts have a negligible signal-to-noise ratio, and do not impact our results. In addition, because at low redshifts reionization sharply suppresses the signal, we set the minimum redshift to z min = 6. The bin size in redshift is set by the bandwidth B . Because foregrounds increasingly contaminate large scale modes (e.g., [81]), we set the minimum k value to k min = 0 . 1 Mpc -1 and the bin size to ∆ k ∥ ( z min , B ) = 0 . 053 Mpc -1 . In addition, since the thermal noise drastically increases at large k , we further restrict the bin range to k max = 1 Mpc -1 . We check that our results are not sensitive to this choice. For the exact values used in our analysis, see tables 2.  \nAdopting a similar approach as 21cmfish [79], 21cmCAST generates a series of configuration files in which all parameters of the model are varied around their fiducial values, θ fid , by a few percent so as to numerically calculate the derivatives of the power spectrum. The configuration files are then given as input for 21cmFAST which produces lightcones of the corresponding simulated universe. Here we work with simulation boxes of size 250 Mpc and divided into 128 3 cells. Eventually, 21cmCAST gathers all the lightcones , bins the power spectra, evaluates the experimental noise using 21cmSense , computes the derivatives numerically and outputs the Fisher matrix. The main differences between 21cmCAST used here and 21cmfish are: (i) more flexibility with the binning choice and (ii) the complete integration of 21cmSense to evaluate the noise directly from the chosen fiducial model.  \nThrough this Fisher matrix analysis, we evaluate the best lower limits that are expected to be set on the dark matter lifetime τ by the HERA experiment within two possible astrophysics scenarios and fixed dark matter mass m χ . The first scenario considers only radiation from PopII-dominated ACGs. Unlike MCGs, ACGs have been observed and their stellar to halo mass relation is well constrained by UV luminosity functions. In this case, the minimal set of astrophysical (with one single population of galaxies) and dark matter parameters considered in our Fisher Matrix analysis is:  \nθ II = { log 10 ( f II ⋆, 10 ) , α II ⋆ , t ⋆ , log 10 ( f II esc , 10 ) , α esc , log 10 ( M turn /M ⊙ ) , log 10 ( L II X ) , E 0 , Γ } , (4.6)  \ni.e., nine parameters with only the last one, the DM decay rate Γ, that sets the DM decay contribution with lifetime τ = 1 / Γ.  \nWe also consider an extended model that includes an additional putative contribution from unseen MCGs. This model is parametrized by:  \nθ II+III = { log 10 ( f II ⋆, 10 ) , α II ⋆ , t ⋆ , log 10 ( f II esc , 10 ) , α esc , log 10 ( L II X ) , E 0 , Γ , log 10 ( f III ⋆, 7 ) , α III ⋆ , log 10 ( f III esc , 7 ) , log 10 ( L III X ) } . (4.7)  \ni.e., 12 parameters characterizing PopII-dominated ACGs and PopIII-dominated MCGs 17 . Because MCGs appear before ACGs, their contribution to the 21-cm PS would be more degenerate with that of DM decay (see figure 2). The bounds on the DM decay properties obtained from θ II+III should then be more conservative than from θ II . The fiducial values considered in this work are tabulated in table 1.  \nNotice that the constraint on the DM lifetime τ , or equivalently its decay rate Γ, is derived assuming a fiducial model with no exotic energy injection. We thus have a fiducial value of the decay rate set to Γ fid = 0, i.e., τ = ∞ . More details on the treatment of this parameter in our Fisher analysis, which comes with some technical difficulties, is provided in appendix B.1.', '4.2 Results': "Triangle plots resulting from our Fisher matrix analysis for the set of parameters θ II and θ II+III are shown in the left and right panel of figure 5 respectively. More precisely we show in both cases the marginalised posterior distributions of the parameters that are most degenerate with the DM decay rate Γ (assuming that they are Gaussian). See appendix B.4 for the full triangle plots with all parameters in θ II and θ II+III . Here we consider decay into electron positron pairs and a DM mass fixed to m χ = 100 MeV /c 2 . Similar plots can be obtained for different DM masses or decays into photons. The dark (light) blue contours represent the 1- (2-)sigma confidence intervals of the two-dimensional marginalized posterior probability distributions while, in the right-most plots, the blue lines are the one dimensional marginalized posteriors for the full set of cosmological parameters θ II and θ II+III . Above the latter plots we provide the fiducial values of the parameters (in black), θ i, fid , and the corresponding one-sigma error (in blue), σ θ i .  \nTo assess the validity our Fisher-matrix analysis, we also provide the posteriors and the estimated 1-sigma errors (in red) that we obtain when the decay parameter, Γ, is not considered in the analysis. The corresponding curves can easily be compared to the results of the Fisher Matrix analysis of [79]. Our estimated errors for the common set of parameters are in very good agreement.  \nThe marginalized error on log 10 ( L II X ) increases by a factor of ∼ 3 when introducing Γ, for our ACG only model. This is expected from the discussion in Secs. 2.3 and 3.2 as both galaxies and DM decay play an important role in the evolution of the IGM temperature in the redshift range of interest. Understandably, we see that Γ and log 10 ( L II X ) are anti-correlated. Γ is also degenerate, although less strongly, with α II ⋆ and E 0 . Indeed, in those cases the marginalized 1-D posteriors are much more weakly affected and the 1-sigma error changes by at most 30%. As visible in eq. (2.17), increasing E 0 makes the X-ray photon spectrum heating the IGM harder. Harder X-rays have longer mean free path (MFP). If the MFP is sufficiently long and heating happens in a highly homogeneous manner, then lowering E 0 leads to more  \n2  \n.  \n0  \nFigure 5 : Marginalized posterior from our Fisher matrix analysis for reduced sets of parameters. The blue contours corresponds to forecasts for a model with 100 MeV /c 2 DM decaying into electron positron pairs. For comparison, we show in red the results of our Fisher forecasts without DM decays. The top right plots show 1D marginalized posteriors and the lower triangle plots show the 2D marginalized posteriors. In the 2D plots, we show with dark (light) colors the 1 (2) σ confidence intervals for the posteriors. Above the top right plots, fiducial values of the parameters are given in black while the 1 sigma errors with (without) decaying dark matter quoted in blue (red). The left (right) panel shows our results for θ II (for θ II+III ). For completeness, the full triangle plots are provided in appendix B.4.  \n<!-- image -->  \nΓ  \nefficient heating. Indeed more photons have MFP shorter than the Hubble length and can be absorbed in the IGM during the EoH. In the latter case, we need less efficient heating from DM, i.e., E 0 and Γ are positively correlated as in the left panel of figure 5. On the other hand, from eq. (2.12), we see that increasing α ∗ penalizes the star formation rate (SFR) of small mass halos that would contribute more to an early heating of the IGM than large mass halos. It is thus positively correlated with Γ that induces an early heating.  \nThe parameters with the strongest degeneracies with Γ in the extended parameter space θ II+III are shown in the right panel of figure 5. This time it is the posterior of the X-ray normalisation factor associated to the PopIII-dominated MCGs (log 10 ( L III X )), instead of that of the PopII-dominated ACGs (log 10 ( L II X )), that is more strongly affected by the introduction of decaying DM. This is to be expected as, by definition, PopIII-dominated MCGs form earlier than PopII-dominated ACGs. Thus, MCGs contribute earlier to IGM heating, similar to decaying DM. Because our log 10 ( L III X ) normalization is defined as the X-ray luminosity above E 0 , lowering E 0 at a fixed value of log 10 ( L III X ) decreases the luminosity at the highest X-ray energies that heat the IGM homogeneously. This decrease in the galactic homogeneous IGM heating can be compensated by increasing the (homogeneous) DM decay heating, i.e., increasing Γ. E 0 and Γ are thus anticorrelated as seen in the right panel of figure 5. On  \n+  \nFigure 6 : Compilation of constraints on the dark matter lifetime (at 95% level) for decay into an electron/positron pair (left panel) and photons (right panel). We superpose the forecasts for the HERA telescope assuming PopII-dominated ACGs only (light blue area, round markers) or PopIIdominated ACGs + PopIII-dominated MCGs (dark blue, cross markers) with existing constraints. Green and red areas show the cosmological bounds set by Lymanα forest [2, 13] and CMB [13, 16] data. Gray areas show astrophysical constraints. In the left panel we have reported the bound from the Voyager I observation of cosmic rays [83, 84] (light gray) and the constraint from X-/ γ -ray experiments [85-88] (dark gray). In the right panel we show the X-ray limit from [89] in light gray colors, as well as bounds from INTEGRAL/SPI [90], COMPTEL, EGRET [85] and Fermi [91] with different shades of gray. The purple areas show the conservative constraints from the Leo T dwarf galaxy [92]. The dashed grey line on the left panel shows the latest XMM-Newton constraint with an improved treatment of cosmic ray propagation [93].  \n<!-- image -->  \nthe other hand, higher f III ∗ , 7 implies stronger feedback which could eventually suppress the low-mass ACGs at lower redshifts. This again suppresses homogeneous heating from stars and can be compensated by larger DM heating. 18  \nWhen considering the extended parameter space, the 1-sigma upper bound on Γ is weaker by approximately one order of magnitude when fixing the decaying DM mass to 100 MeV /c 2 . Indeed, when considering the minimum set θ II , the heating effect of ACGs like galaxies only, which appear relatively late compared to the heating of DM, and the wide redshift range probed by HERA allows to disentangle the DM effect from the astrophysics and sets strong constraints on the DM decay rate. In the case of MCGs+ACGs, with MCGs heating the IGM at earlier times than ACGs, the effect of DM is more easily drawn by the effect of the astrophysics parameters and the limits on the decay rate are much weaker.  \nOur main results, presented in figure 6, showcase the lower bound at a 95% confidence level (CL) on the lifetime of dark matter (DM) derived from our Fisher matrix forecasts based on HERA specifications. Black lines with bullets are obtained for one single population of galaxies while black lines with crosses assume ACGs+MCGs. The blue area below these curves are excluded at 95% CL. As expected, when considering the ACGs+MCGs scenario, heating from galaxies competes with DM heating earlier and the DM heating parameters  \nbecome more difficult to constrain. The lower bound on the DM life time becomes thus less stringent in ACGs+MCGs scenario (crosses) than in the AGCs only case (bullets).  \nMore precisely, our bounds are actually obtained conducting a Fisher matrix analysis at each value of the DM mass marked with a cross or a bullet. We show the results both for decays into electron positron pairs (left panel) and photons (right panel). The 95% CL bound just simply corresponds to τ ( m χ ) < (2 × σ Γ ( m χ )) -1 , where σ Γ is obtained saturating the Cr'amer-Rao inequality as in (4.2). It is thus an optimistic limit. 19 For decays into e + e -, a few 100 MeV /c 2 DM gets the most stringent 21cm bounds on Γ while, for decays into photons, it is the case for the lowest DM masses (with m χ < MeV /c 2 ). The mass dependency of our 21cm forecasts actually traces back to the dependency of energy deposition efficiencies, f a ( x e , z ), on the injected energy (see e.g., [8, 74, 75]). 20  \nIn figure 6, we also show the lower bounds arising from cosmology probes such as Lymanα forest [2] (green) or CMB [13, 16] (red). 21-cm measurements with HERA could improve by up to 3 orders of magnitude the current limits on the DM lifetime, when considering PopII-dominated ACGs only (black line with bullets). Note that mass dependency seen in the Lymanα constraints (green area) is similar to the one of our 21cm bound forecasts (bullet or crosses). This is because Lymanα data also probe the DM heating and is thus sensitive to f heat energy injection dependency. Yet, these data probe the IGM temperature at lower z , when T k is orders of magnitude larger, resulting into weaker constraints compared to our 21cm forecast (for f a as a function of redshift and energy injection see also e.g. [8]). Even when PopIII-dominated MCGs (dark blue area, line with crosses) are included in our analysis, the 21cm constraints still improve by a factor of 10 to 100 compared to existing cosmology constraints in the lower DM mass range.  \nIn figure 6, we additionally compare our forecast against constraints from the Voyager I observation of cosmic rays [83, 84] and the result of X- or γ -ray experiments [85-87, 89] including those coming from INTEGRAL/SPI [90], COMPTEL, EGRET [85] and Fermi [91]. We find that existing constraints for dark matter heavier than 1 GeV /c 2 (or 100keV /c 2 ) for decays into e + e -(or photons) remain competitive and our 21cm forecast for 1000 hours of HERA observation is unlikely to improve these limits in the higher DM mass range.", '5 Conclusion and future perspectives': 'Determining the nature of DM is one of the major goals of particle physics and cosmology today. It has brought us to explore all possible ways to shed new light on its properties. Cosmology probes are well known to set stringent bounds on DM energy injection into the IGM. In this work, we have carefully studied the prospects for 21cm cosmology to probe dark matter decays into electron-positron pairs and photons. For that purpose, we have developed a new version of the public code 21cmFAST , that interfaced with the public package DarkHistory , accounts for DM energy deposition into heat, ionization and excitation of the medium and its effect on the 21cm signal. In particular, DM energy injection is expected to  \nmainly affect the 21cm global signal and its power spectrum as an exotic source of heating, already efficient before PopII-dominated ACGs light on.  \nAs is well known, DM decays give rise to a relatively late time energy deposition into the medium. This makes late time probes, such as Lymanα forest or 21cm cosmology, very interesting targets to detect the DM imprint. In this work, we focus on the effect of DM on the 21cm signal power spectrum and prospects for constraints on the DM lifetime by the HERA interferometer. This telescope will enable us to explore a vast range of redshifts, stretching from the Epoch of Reionization to Cosmic Dawn, with exceptional precision. This capability is of paramount importance because DM is not the sole contributor to the heating process, and it is crucial to distinguish its distinct signature from that generated by X-rays emitted from the first galaxies. In our work, we argue that their different imprint in the probed redshift range is the key to obtain competitive constraints with respect to existing probes, both from cosmology and from astro-particle physics experiments.  \nIn order to provide quantitative forecasts, we have performed for the first time a dedicated Fisher Matrix analysis considering the HERA telescope and two different astrophysics scenarios including DM decays. For the Fisher matrix forecast we have developed our own tool, 21cmCAST , that interfaces with 21cmSense to evaluate the anticipated experimental errors from 1000 hours of observation using all 331 antennas from HERA, extracting the expected marginalized error on a set of astrophysics parameter of our choice.  \nOur results, summarized in figure 6, are very promising. When considering the minimal astrophysics scenario, HERA is expected to improve on existing cosmology constraints (from CMB and Lymanα probes) on the DM lifetime by up to 3 orders of magnitude. We also compare these prospects to the case where the astrophysics model includes both PopIIdominated ACGs and PopIII-dominated MCGs. Similarly to DM, MCGs give rise to a new source of IGM heating before PopII-dominated ACGs light on, partially drowning the DM signal. Nevertheless, even in the latter case, HERA can improve on existing cosmology constraints by a factor of 10 to 100. Finally, compared to existing γ -ray and cosmic-ray limits, HERA is expected to be a key player in constraining DM candidates decaying to e + e -in the mass range m χ ≲ 2 GeV /c 2 . For decays into photons, HERA improves on other searches in the low mass range for m χ ≲ few MeV /c 2 .  \nIn the light of these very good prospects, it would be necessary to refine this study making use of a more advanced statistical analysis such as e.g., Bayesian inference. The latter is much more time consuming/computationally expensive but would give rise to more realistic estimates of the prospect for constraining DM decays. We will also vary the complexity of the astrophysical models, comparing the Bayesian evidences (e.g., [43]). In addition, DM could also heat the IGM through annihilations. In this case, one should take special care in the treatment of the late time boost arising for structure formation. This induces yet another modeling uncertainty, but is worth exploring given the sensitivity of upcoming 21cm measurements to IGM heating during the CD.', 'Acknowledgments': "We thank S. Junius for collaboration at early stages of this work as well as T. Slatyer and Q. Decant for useful discussions on 21cm cosmology and DM imprint. GF acknowledges support of the ARC program of the Federation Wallonie-Bruxelles and of the Excellence of Science (EoS) project No. 30820817 - be.h 'The H boson gateway to physics beyond the Standard Model'. LLH is supported by the Fonds de la Recherche Scientifique F.R.S.-FNRS  \nthrough a research associate position and acknowledges support of the FNRS research grant number F.4520.19, the ARC program of the Federation Wallonie-Bruxelles and the IISN convention No. 4.4503.15. Computational resources have been provided by the Consortium des Equipements de Calcul Intensif (CECI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region of Belgium. YQ is supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project #CE170100013. Part of this work was performed on the OzSTAR and Gadi national computational facilities in Australia.", 'A Extra information on IGM heating sources': 'Here we provide further details on X-ray energy injection in appendix A.1 and we discuss the impact of MCGs and the DM heating when considering both ACGs and MCGs in appendix A.2.', 'A.1 X-ray heating contributions': "The X-ray emissivity ϵ X introduced in eq. (2.18) contributes to secondary ionisations, heating and the Lyman alpha background appearing in eqs. (2.5) and (2.6). We first introduce the total integrated radiation intensity  \nI X = (1 + z ) 3 4 π ∫ ∞ z d z ' c d t d z ' ϵ X exp( -τ X ) , (A.1)  \nwhere τ X stands for the X-rays optical depth in the IGM (see more in Ref [34]). The latter depends on the specific X-ray luminosity per SFR whose amplitude and softness/hardness is set by the parameters L I / II X and E 0 , see eq. (2.17). The ionization and heating rates from X-ray photons (i.e., Λ X ion in eq. (2.5) and ϵ X heat in eq. (2.6)) as well as their contribution to the Lymanα radiation (i.e., J X α in equation 2.8) follows  \nΛ X ion = ∫ ∞ E 0 d E 4 πI X E ∑ j x j σ j f j [ ( E -E j th ) ∑ k f k ion E k th +1 ] , (A.2)  \nϵ X heat = ∫ ∞ E 0 d E 4 πI X E ∑ j x j σ j f j ( E -E j th ) f heat , (A.3)  \nJ X α = cn b 4 πH ( z ) ν α ∫ ∞ E 0 d E 4 πI X E ∑ j x j σ j f j ( E -E j th ) f α hν α . (A.4)  \nand  \nHere, f j refers to the number fraction of each species, j , in HI, HeI,and HeII, with σ j and E j th being their corresponding cross-section and energy for ionization; f heat , f k ion and f α represent the fraction of the electron energy after ionization, E -E j th , that contributes to heating, secondary ionization or emitting Lymanα photons ( ν α ≡ 2 . 47 × 10 15 Hz) of each species [94]; and x j ≡ 1 -x e for HI and HeI or x e for HeII represents the secondary ionization fractions (See more in [43]).  \nFigure 7 : Similar as figure 1 including MCGs and varying L III X . We show the 1 σ error (instead of 2 σ ) for a more readable figure.  \n<!-- image -->  \nFigure 8 : Same as figure 3 including MCGs and showing the 1 σ error (instead of 2 σ for a more readable figure).  \n<!-- image -->", 'A.2 MCGs vs DM decays imprint on 21cm signal': 'As shown in figure 2, the X-ray heating rate from PopIII-dominated MCGs decreases more smoothly in redshift than that from ACGs. This is similar to the case of the DM heating rate that is roughly constant. It is thus expected that the heating parameters associated to MCGs, including L III X , shall be more strongly degenerate with Γ, which normalises the DM decay rate, than the ones associated to ACGs. Figures 7 and 8, illustrating the dependence of the 21cm signal on L III X and Γ in an ACGs+MCGs scenario, can be compared to figures 1 and 3 for ACGs only. Including MCGs, we see that the fiducial scenario (black curve) present a shallower and earlier absorption in the global signal. Furthermore, the power  \nFigure 9 : Time slices of the differential brightness temperature in our (250 Mpc) 3 large simulation box . We show results for the ACGs + MCGs fiducial model - with log 10 ( L III X ) = 40 and τ = ∞ - (upper panel), log 10 ( L III X ) = 41 (middle panel), and for a 100 MeV /c 2 DM with τ = 10 27 . 5 s (bottom panel).  \n<!-- image -->  \nspectrum is suppressed at larger redshifts (between redshifts 12 and 16) and the EoR peak is wider (between redshifts 6 and 10). Figure 7 illustrates the effect of varying L III X , instead of L II X in figure 1, on the 21cm signal while figure 8 shows variations in the signal due to the DM life time. In the MCGs+ACGs case, it is more difficult to disentangle DM imprint from astrophysics. Indeed, the features in the power spectrum that can help discriminate the DM heating from first galaxies X-ray heating are pushed to larger redshifts (where the instrument is less sensitive) and are more difficult to detect due to a smaller signal-to-noise ratio. In figures 7 and 8, we show the noise at 1 σ level with the blue (experimental noise only) and green (total noise) bars - while it was shown at 2 σ in figures 1 and 3. Nonetheless, we observe that modifying the fiducial values of L III X and τ to 10 40 . 5 and 10 27 . 5 s are already sufficient to shift the signal outside the 1 σ uncertainty area for some of the probed redshifts and k values. We also provide time slices of the differential brightness temperature in the MCGs+ACGs scenario in figure 9 that can be compared to figure 4.', 'B Fisher matrix analysis': 'In this appendix we give more details on the mathematics behind the Fisher analysis, the binning choice we choose and the treatment of the decay parameter. Furthermore, we provide full triangle plots for all parameters considered in the analysis.', 'B.1 General method': "Here we detail the formalism of the Fisher matrix inference method. Our analysis is inspired and similar to that introduced in ref. [79]. For the following discussion we assume that the power spectrum is binned into a grid of N k k -modes over N z redshifts. We define the vector of length N ≡ N k × N z ,  \nX ≡ ( δT b 2 ( z 1 )∆ 2 21 ( k 1 , z 1 ) , δT b 2 ( z 1 )∆ 2 21 ( k 1 , z 1 ) , . . . , δT b 2 ( z N z )∆ 2 21 ( k N k , z N z )) (B.1)  \ncapturing all the information that can experimentally be obtained about the 21cm powerspectrum. From measured X , one could reconstruct the posterior distribution on the parameters using an inference algorithm. However, because we aim at forecasting sensitivities  \nTable 2 : Left panel . Redshift bin edges ˆ z n and associated centers z n . In parenthesis we put the associated frequencies of the redshifted 21cm line. Right panel . Mode bin edges ˆ k n and associated centers k n .  \nof forthcoming experimental runs, one can only rely on the theoretical predictions obtained from 21cmFAST and that depend on the choice of model to describe the Universe. A Fisher analysis is thus adequate to obtain reliable estimates in a relatively short time (without requiring extensive computational resources).  \nWe assume that X is randomly distributed according to the likelihood L : ( X , θ ) ↦→ L ( X | θ ) with θ a vector of parameters. The Fisher information matrix F associated to L and evaluated for the parameters θ is given by  \nF ij = -〈 ∂ 2 ln L ∂θ i ∂θ j 〉 ≡ -E X [ ∂ 2 ∂θ i ∂θ j ln L ( X | θ ) ∣ ∣ ∣ θ ] (B.2)  \n∣ where E X is the expectation value associated to L . The Cram'er-Rao limit [76-78] ensures that the covariance matrix associated to the posterior distribution on the set of parameters θ satisfies the inequality  \n( C θ ) ij ≥ ( F -1 ) ij ∀ ( i, j ) . (B.3)  \nWe assume the optimistic case and define the covariance of the posterior as the inverse of the Fisher matrix. We consider that X follows a multivariate Gaussian distribution with covariance matrix C X (independent of θ ) 21 and mean µ X ( θ ): X ∼ N ( µ X ( θ ) , C X ). Thus, the likelihood and the Fisher matrix elements take the form,  \nL ( x | θ ) = e -1 2 [ x -µ X ( θ )] T C -1 X [ x -µ X ( θ )] √ (2 π ) N | det C X | and F ij = ∂ µ T X ( θ ) ∂θ i C -1 X ∂ µ X ( θ ) ∂θ j , (B.4)  \nwhere T indicates vector transposition. Therefore, we can estimate the covariance matrix of the parameters θ from a matrix-sum on the derivative of the mean value of the data. Having no real data set at our disposal, we can only assume that the true Universe is similar to that simulated in 21cmFAST for a chosen fiducial model θ = θ fid . Hence, µ X is directly given by the binned output obtained from that fiducial model and one computes derivatives by slightly varying the parameters around the fiducial values.", 'B.2 Details on the binning choice': 'Notwithstanding that 21cmCAST can work for different binning choice, by default, we bin the mock data according to the frequency bandwidth B of the instrument. The redshift bin width is straightforwardly related to the frequency bandwidth through the redshifted 21cm line frequency ν ( z ) = ν 21 / (1+ z ). We write the minimum redshift considered for the analysis as z min = ˆ z 0 and define frequency bins with edges { ˆ ν n } n ≥ 0 and centers { ν n } n ≥ 1 by  \n{ ˆ ν n = ν (ˆ z 0 ) -nB ν n = ν (ˆ z 0 ) -( n -1 2 ) B. (B.5)  \nSubsequently, z -bins are defined with edges ˆ z n = ν 21 / ˆ ν n -1 and centers z n = ν 21 /ν n -1. The sequence runs up to n = N z , until ˆ z N z +1 becomes greater than a fixed threshold z max . A mode k is given by the quadratic sum of its two components, k ∥ and k ⊥ respectively parallel and perpendicular to the line of sight. The minimal accessible value of | k ⊥ | is fixed by the minimal baseline distance between the antennas. However, that of k ∥ happens to be orders of magnitude larger and is fixed by the bandwidth [30]. One can associate a frequency range ∆ ν to a range of the parallel projection of modes using  \n∆ k ∥ ( z, ∆ ν ) ≡ 2 π ν 21 ∆ ν H ( z ) (1 + z ) 2 , (B.6)  \ndecreasing function of the redshift. Thus, we define the k -bin width from the maximal possible value of ∆ k ∥ at ˆ z 0 for ∆ ν = B . On a range [ k min , k max ], the bin edges are  \nˆ k n ≡ max { k min , ∆ k ∥ (ˆ z 0 , B ) } + n ∆ k ∥ (ˆ z 0 , B ) , (B.7)  \nup to n = N k where ˆ k n becomes larger than k max . The instrument being limited by its spectral resolution δν , k max can not be arbitrary and should at least satisfy k max < ∆ k ∥ ( z max , δν ) / 2.  \nWe put in table 2 the fixed grid in redshifts and modes that we use to bin the output of 21cmFAST and perform the Fisher matrix analysis. Note that it appears we could also have simply defined the center of the redshift bins z n as  \n˜ z n = ˆ z n + ˆ z n -1 2 since ˜ z n z n = (ˆ ν n +1 / 2) 2 ˆ ν n (ˆ ν n +1) (B.8)  \ngoes to 1 when ˆ ν n ≫ 1. Indeed, in practice ˆ ν n > 65 for the redshift range we consider, thus the two definitions are completely equivalent.  \nFigure 10 : Triangle plot result of the Fisher analysis with the contribution of a single polulation of PopII dominated ACGs. In blue we consider a 100 MeV /c 2 DM decaying into electron positron pairs. For comparison, in red we show the result of the Fisher forecast without DM decay. We can conclude that DM decay mostly impacts on the reconstruction of L II X .  \n<!-- image -->', 'B.3 Treatment of dark matter decay in the Fisher forecast': "The constraint on the DM lifetime τ , or equivalently its decay rate Γ, is derived assuming a fiducial model with no exotic energy injection. We thus have a fiducial value of the decay rate set to Γ fid = 0 s -1 . That is to say, we want to quantify by how much Γ can be larger than zero before it cannot be neglected in the parameter reconstruction. This gives rise to technical difficulties in the numerical evaluation of the derivative of the power spectrum with  \nFigure 11 : Same as figure 10 but adding the contribution of PopIII-dominated MCGs. Here, DM decay impacts more strongly the reconstruction of L III X than the one of L II X .  \n<!-- image -->  \nrespect to the decay rate - see eq. (4.3). Firstly, because the fiducial value for Γ equal zero, 22 one cannot consider its logarithm even-though Γ has to vary on many orders of magnitude. For the that reason, varying the parameters by a few percent to compute the derivative also does not make sense here. Secondly, it is not possible to evaluate the power spectrum for Γ < 0. As a result, we use a one side derivative scheme (instead of a centered scheme for the other parameters),  \n∂f ∂ Γ ∣ ∣ ∣ Γ=0 = f ( ϵ Γ ) -f (0) ϵ Γ , (B.9)  \n+  \nFigure 12 : Constraints on the dark matter lifetime (at 95% level) for decay into an electron/positron pair (left panel) and photons (right panel). We superpose the forecasts for the HERA telescope assuming PopII-dominated ACGs only and log 10 ( L II X ) = 41 (dashed blue line with round markers), PopII-dominated ACGs only and log 10 ( L II X ) = 40 (solid dark line with round markers), or PopIIdominated ACGs + PopIII-dominated MCGs with log 10 ( L II , III X ) = 40 (solid dark line with crosses) with existing cosmological constraints. We refer to figure 6 for more details.  \n<!-- image -->  \nto evaluate the derivative of a function f with respect to Γ, with ϵ Γ , a small parameter. In practice here we have f = δT b 2 ∆ 2 21 .  \nThe choice of 'small value' for ϵ Γ used in the numerical analysis is quite critical as Γ vary over orders of magnitude and one needs ∂f/∂ Γ to be independent of the choice of ϵ Γ . Indeed, for choices of ϵ Γ too small, f ( ϵ Γ ) is expected to be essentially identical to f (0) up to the numerical noise and one has then ∂f/∂ Γ ∝ 1 /ϵ Γ . On the other hand, when ϵ Γ is too large, the expression (B.9) is no longer valid for a derivative. In order to make a sensible choice, we have computed the marginalized error on Γ, σ Γ = √ ( F -1 ) i Γ i Γ , for various values of ϵ Γ . We have then selected a value of ϵ Γ at which σ Γ is independent ϵ Γ , i.e., for which d σ Γ / d ϵ Γ ∼ 0. In the analysis performed in this work, we find that the correct choice for ϵ Γ is a few percent of the error σ Γ , i .e, ϵ Γ = 10 -27 -10 -29 s -1 .", 'B.4 Full triangle plots': 'The full triangle plots resulting from our Fisher Matrix forecasts are shown in figure 10, completing the left plot of figure 5 for the analysis involving PopII-dominated ACGs. Furthermore, figure 11, completes the right plot of figure 10 for the analysis involving both PopII-dominated ACGs and PopIII-dominated MCGs.', 'C Impact of the X-ray normalisation on the forecast': 'In this section we show the results of our Fisher analysis in the ACGs-only scenario but considering log 10 ( L II X ) = 41 (instead of 40) as the fiducial value. In figure 12 the 95% CL bound for the new fiducial value of L II X is shown (dashed blue line with round markers) together with the other cosmological constraints from CMB and Lymanα . For comparison,  \nthe bounds obtained in the ACGs-only scenario with log 10 ( L II X ) = 40 (solid dark line with round markers) and in the ACGs+MCGs scenario with log 10 ( L II / III X ) = 40 (solid dark line with crosses) are also plotted. With a larger X-ray injection the contribution of dark matter decay is harder to disentangle from the astrophysical sources (see figure 2). Consequently, there is a drop of almost an order of magnitudes between the fiducial log 10 ( L II X ) = 40 and the fiducial log 10 ( L II X ) = 41. Nonetheless, the constraints with ACGs and log 10 ( L II X ) = 41 remains slightly stronger than considering MCGs and than the CMB and Lymanα bounds.', 'References': "- [1] Planck collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters , Astron. Astrophys. 641 (2020) A6, [ 1807.06209 ].\n- [2] H. Liu, W. Qin, G. W. Ridgway and T. R. Slatyer, Lymanα constraints on cosmic heating from dark matter annihilation and decay , Phys. Rev. D 104 (2021) 043514, [ 2008.01084 ].\n- [3] J. M. Shull and M. E. van Steenberg, X-ray secondary heating and ionization in quasar emission-line clouds , Astrophys. J. 298 (1985) 268-274.\n- [4] J. A. Adams, S. Sarkar and D. W. Sciama, CMB anisotropy in the decaying neutrino cosmology , Mon. Not. Roy. Astron. Soc. 301 (1998) 210-214, [ astro-ph/9805108 ].\n- [5] X.-L. Chen and M. Kamionkowski, Particle decays during the cosmic dark ages , Phys. Rev. D70 (2004) 043502, [ astro-ph/0310473 ].\n- [6] N. Padmanabhan and D. P. Finkbeiner, Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects , Phys. Rev. D 72 (2005) 023508, [ astro-ph/0503486 ].\n- [7] T. R. Slatyer, N. Padmanabhan and D. P. Finkbeiner, CMB Constraints on WIMP Annihilation: Energy Absorption During the Recombination Epoch , Phys. Rev. D 80 (2009) 043526, [ 0906.1197 ].\n- [8] T. R. Slatyer, Indirect Dark Matter Signatures in the Cosmic Dark Ages II. Ionization, Heating and Photon Production from Arbitrary Energy Injections , Phys. Rev. D 93 (2016) 023521, [ 1506.03812 ].\n- [9] L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz and A. C. Vincent, Constraints on dark matter annihilation from CMB observationsbefore Planck , JCAP 07 (2013) 046, [ 1303.5094 ].\n- [10] R. Diamanti, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz and A. C. Vincent, Constraining Dark Matter Late-Time Energy Injection: Decays and P-Wave Annihilations , JCAP 02 (2014) 017, [ 1308.2578 ].\n- [11] V. Poulin, P. D. Serpico and J. Lesgourgues, Dark Matter annihilations in halos and high-redshift sources of reionization of the universe , JCAP 12 (2015) 041, [ 1508.01370 ].\n- [12] B. Bolliet, J. Chluba and R. Battye, Spectral distortion constraints on photon injection from low-mass decaying particles , Mon. Not. Roy. Astron. Soc. 507 (2021) 3148-3178, [ 2012.07292 ].\n- [13] F. Capozzi, R. Z. Ferreira, L. Lopez-Honorez and O. Mena, CMB and Lymanα constraints on dark matter decays to photons , 2303.07426 .\n- [14] H. Liu, W. Qin, G. W. Ridgway and T. R. Slatyer, Exotic energy injection in the early universe II: CMB spectral distortions and constraints on light dark matter , 2303.07370 .\n- [15] M. Lucca, N. Schoneberg, D. C. Hooper, J. Lesgourgues and J. Chluba, The synergy between CMB spectral distortions and anisotropies , JCAP 02 (2020) 026, [ 1910.04619 ].\n- [16] T. R. Slatyer and C.-L. Wu, General Constraints on Dark Matter Decay from the Cosmic Microwave Background , Phys. Rev. D 95 (2017) 023010, [ 1610.06933 ].  \n- [17] M. Vald'es, C. Evoli and A. Ferrara, Particle energy cascade in the Intergalactic Medium , Mon. Not. Roy. Astron. Soc. 404 (2010) 1569-1582, [ 0911.1125 ].\n- [18] C. Evoli, A. Mesinger and A. Ferrara, Unveiling the nature of dark matter with high redshift 21 cm line experiments , JCAP 11 (2014) 024, [ 1408.1109 ].\n- [19] L. Lopez-Honorez, O. Mena, A. Molin'e, S. Palomares-Ruiz and A. C. Vincent, The 21 cm signal and the interplay between dark matter annihilations and astrophysical processes , JCAP 08 (2016) 004, [ 1603.06795 ].\n- [20] H. Liu, T. R. Slatyer and J. Zavala, Contributions to cosmic reionization from dark matter annihilation and decay , Phys. Rev. D 94 (2016) 063507, [ 1604.02457 ].\n- [21] H. Liu and T. R. Slatyer, Implications of a 21-cm signal for dark matter annihilation and decay , Phys. Rev. D 98 (2018) 023501, [ 1803.09739 ].\n- [22] G. D'Amico, P. Panci and A. Strumia, Bounds on Dark Matter annihilations from 21 cm data , Phys. Rev. Lett. 121 (2018) 011103, [ 1803.03629 ].\n- [23] K. Cheung, J.-L. Kuo, K.-W. Ng and Y.-L. S. Tsai, The impact of EDGES 21-cm data on dark matter interactions , Phys. Lett. B 789 (2019) 137-144, [ 1803.09398 ].\n- [24] A. Mitridate and A. Podo, Bounds on Dark Matter decay from 21 cm line , JCAP 05 (2018) 069, [ 1803.11169 ].\n- [25] S. Clark, B. Dutta, Y. Gao, Y.-Z. Ma and L. E. Strigari, 21 cm limits on decaying dark matter and primordial black holes , Phys. Rev. D 98 (2018) 043006, [ 1803.09390 ].\n- [26] S. K. Acharya, B. Cyr and J. Chluba, The role of soft photon injection and heating in 21 cm cosmology , Mon. Not. Roy. Astron. Soc. 523 (2023) 1908-1918, [ 2303.17311 ].\n- [27] N. Hiroshima, K. Kohri, T. Sekiguchi and R. Takahashi, Impacts of new small-scale N-body simulations on dark matter annihilations constrained from cosmological 21-cm line observations , Phys. Rev. D 104 (2021) 083547, [ 2103.14810 ].\n- [28] S. Mittal, A. Ray, G. Kulkarni and B. Dasgupta, Constraining primordial black holes as dark matter using the global 21-cm signal with X-ray heating and excess radio background , JCAP 03 (2022) 030, [ 2107.02190 ].\n- [29] M. Cirelli, F. Iocco and P. Panci, Constraints on Dark Matter annihilations from reionization and heating of the intergalactic gas , JCAP 10 (2009) 009, [ 0907.0719 ].\n- [30] HERA collaboration, Z. Abdurashidova et al., First Results from HERA Phase I: Upper Limits on the Epoch of Reionization 21 cm Power Spectrum , Astrophys. J. 925 (2022) 221, [ 2108.02263 ].\n- [31] HERA collaboration, Z. Abdurashidova et al., Improved Constraints on the 21 cm EoR Power Spectrum and the X-Ray Heating of the IGM with HERA Phase I Observations , Astrophys. J. 945 (2023) 124, [ 2210.04912 ].\n- [32] T. Fragos, B. D. Lehmer, S. Naoz, A. Zezas and A. Basu-Zych, Energy Feedback from X-Ray Binaries in the Early Universe , ApJ 776 (Oct, 2013) L31, [ 1306.1405 ].\n- [33] H. D. Kaur, Y. Qin, A. Mesinger, A. Pallottini, T. Fragos and A. Basu-Zych, The 21-cm signal from the cosmic dawn: metallicity dependence of high-mass X-ray binaries , MNRAS 513 (July, 2022) 5097-5108, [ 2203.10851 ].\n- [34] A. Mesinger, S. Furlanetto and R. Cen, 21cmFAST: A Fast, Semi-Numerical Simulation of the High-Redshift 21-cm Signal , Mon. Not. Roy. Astron. Soc. 411 (2011) 955, [ 1003.3878 ].\n- [35] S. G. Murray, B. Greig, A. Mesinger, J. B. Mu˜noz, Y. Qin, J. Park et al., 21cmFAST v3: A Python-integrated C code for generating 3D realizations of the cosmic 21cm signal , J. Open Source Softw. 5 (2020) 2582, [ 2010.15121 ].  \n- [36] H. Liu, G. W. Ridgway and T. R. Slatyer, Code package for calculating modified cosmic ionization and thermal histories with dark matter and other exotic energy injections , Phys. Rev. D 101 (2020) 023530, [ 1904.09296 ].\n- [37] Y. Sun and T. R. Slatyer, Modeling early-universe energy injection with dense neural networks , Phys. Rev. D 107 (2023) 063541, [ 2207.06425 ].\n- [38] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., Planck 2018 results. VI. Cosmological parameters , A&A 641 (Sept., 2020) A6, [ 1807.06209 ].\n- [39] S. R. Furlanetto, S. P. Oh and F. H. Briggs, Cosmology at low frequencies: The 21 cm transition and the high-redshift Universe , Phys. Rep. 433 (Oct., 2006) 181-301, [ astro-ph/0608032 ].\n- [40] B. Greig and A. Mesinger, 21CMMC: an MCMC analysis tool enabling astrophysical parameter studies of the cosmic 21 cm signal , Mon. Not. Roy. Astron. Soc. 449 (2015) 4246-4263, [ 1501.06576 ].\n- [41] Y. Mao, P. R. Shapiro, G. Mellema, I. T. Iliev, J. Koda and K. Ahn, Redshift-space distortion of the 21-cm background from the epoch of reionization - I. Methodology re-examined , MNRAS 422 (May, 2012) 926-954, [ 1104.2094 ].\n- [42] H. Jensen, K. K. Datta, G. Mellema, E. Chapman, F. B. Abdalla, I. T. Iliev et al., Probing reionization with LOFAR using 21-cm redshift space distortions , MNRAS 435 (Oct., 2013) 460-474, [ 1303.5627 ].\n- [43] Y. Qin, A. Mesinger, J. Park, B. Greig and J. B. Mu˜noz, A tale of two sites - I. Inferring the properties of minihalo-hosted galaxies from current observations , MNRAS 495 (June, 2020) 123-140, [ 2003.04442 ].\n- [44] S. R. Furlanetto, M. Zaldarriaga and L. Hernquist, The Growth of H II Regions During Reionization , ApJ 613 (Sept., 2004) 1-15, [ astro-ph/0403697 ].\n- [45] C. M. Hirata, Wouthuysen-Field coupling strength and application to high-redshift 21-cm radiation , MNRAS 367 (Mar., 2006) 259-274, [ astro-ph/0507102 ].\n- [46] V. Bromm and R. B. Larson, The first stars , Annual Review of Astronomy and Astrophysics 42 (2004) 79-118, [ https://doi.org/10.1146/annurev.astro.42.053102.134034 ].\n- [47] H. Xu, M. L. Norman, B. W. O'Shea and J. H. Wise, Late Pop III Star Formation During the Epoch of Reionization: Results from the Renaissance Simulations , ApJ 823 (June, 2016) 140, [ 1604.03586 ].\n- [48] R. Schneider, A. Ferrara, P. Natarajan and K. Omukai, First stars, very massive black holes, and metals , The Astrophysical Journal 571 (may, 2002) 30.\n- [49] R. H. Mebane, J. Mirocha and S. R. Furlanetto, The Persistence of Population III Star Formation , Monthly Notices of the Royal Astronomical Society 479 (07, 2018) 4544-4559, [ https://academic.oup.com/mnras/article-pdf/479/4/4544/25180578/sty1833.pdf ].\n- [50] M. Stefanon, R. J. Bouwens, I. Labb'e, G. D. Illingworth, V. Gonzalez and P. A. Oesch, Galaxy Stellar Mass Functions from z ∼ 10 to z ∼ 6 using the Deepest Spitzer/IRAC Data: No Significant Evolution in the Stellar-to-Halo Mass Ratio of Galaxies in the First Gyr of Cosmic Time , 2103.16571 .\n- [51] M. Shuntov, H. J. McCracken, R. Gavazzi, C. Laigle, J. R. Weaver, I. Davidzon et al., COSMOS2020: Cosmic evolution of the stellar-to-halo mass relation for central and satellite galaxies up to z ∼ 5 , A&A 664 (Aug., 2022) A61, [ 2203.10895 ].\n- [52] J. Park, A. Mesinger, B. Greig and N. Gillet, Inferring the astrophysics of reionization and cosmic dawn from galaxy luminosity functions and the 21-cm signal , Mon. Not. Roy. Astron. Soc. 484 (2019) 933-949, [ 1809.08995 ].  \n- [53] S. J. Mutch, P. M. Geil, G. B. Poole, P. W. Angel, A. R. Duffy, A. Mesinger et al., Dark-ages reionization and galaxy formation simulation - III. Modelling galaxy formation and the epoch of reionization , MNRAS 462 (Oct., 2016) 250-276, [ 1512.00562 ].\n- [54] J. S. B. Wyithe and A. Loeb, A suppressed contribution of low-mass galaxies to reionization due to supernova feedback , MNRAS 428 (Jan., 2013) 2741-2754, [ 1209.2215 ].\n- [55] R. K. Sheth, H. J. Mo and G. Tormen, Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes , MNRAS 323 (May, 2001) 1-12, [ astro-ph/9907024 ].\n- [56] E. Sobacchi and A. Mesinger, How does radiative feedback from an ultraviolet background impact reionization? , MNRAS 432 (July, 2013) 3340-3348, [ 1301.6781 ].\n- [57] R. Barkana and A. Loeb, In the beginning: the first sources of light and the reionization of the universe , Physics Reports 349 (2001) 125-238.\n- [58] J. B. Mu˜noz, Y. Qin, A. Mesinger, S. G. Murray, B. Greig and C. Mason, The impact of the first galaxies on cosmic dawn and reionization , MNRAS 511 (Apr., 2022) 3657-3681, [ 2110.13919 ].\n- [59] W. Qin, J. B. Munoz, H. Liu and T. R. Slatyer, Birth of the first stars amidst decaying and annihilating dark matter , 2308.12992 .\n- [60] R. J. Bouwens, G. D. Illingworth, P. A. Oesch, M. Trenti, I. Labb'e, L. Bradley et al., UV Luminosity Functions at Redshifts z ∼ 4 to z ∼ 10: 10,000 alaxies from HST Legacy Fields , ApJ 803 (Apr., 2015) 34, [ 1403.4295 ].\n- [61] R. J. Bouwens, P. A. Oesch, I. Labb'e, G. D. Illingworth, G. G. Fazio, D. Coe et al., The Bright End of the z ∼ 9 and z ∼ 10 UV Luminosity Functions Using All Five CANDELS Fields* , ApJ 830 (Oct., 2016) 67, [ 1506.01035 ].\n- [62] P. A. Oesch, R. J. Bouwens, G. D. Illingworth, I. Labb'e and M. Stefanon, The Dearth of z ∼ 10 Galaxies in All HST Legacy Fields-The Rapid Evolution of the Galaxy Population in the First 500 Myr , ApJ 855 (Mar, 2018) 105, [ 1710.11131 ].\n- [63] A. Das, A. Mesinger, A. Pallottini, A. Ferrara and J. H. Wise, High-mass X-ray binaries and the cosmic 21-cm signal: impact of host galaxy absorption , MNRAS 469 (July, 2017) 1166-1174, [ 1702.00409 ].\n- [64] S. Mineo, M. Gilfanov and R. Sunyaev, X-ray emission from star-forming galaxies - I. High-mass X-ray binaries , MNRAS 419 (Jan., 2012) 2095-2115, [ 1105.4610 ].\n- [65] B. D. Lehmer, R. T. Eufrasio, A. Basu-Zych, K. Doore, T. Fragos, K. Garofali et al., The Metallicity Dependence of the High-mass X-Ray Binary Luminosity Function , ApJ 907 (Jan., 2021) 17, [ 2011.09476 ].\n- [66] T. Fragos, B. Lehmer, M. Tremmel, P. Tzanavaris, A. Basu-Zych, K. Belczynski et al., X-Ray Binary Evolution Across Cosmic Time , ApJ 764 (Feb., 2013) 41, [ 1206.2395 ].\n- [67] M. Brorby, P. Kaaret, A. Prestwich and I. F. Mirabel, Enhanced X-ray emission from Lyman break analogues and a possible LX-SFR-metallicity plane , Monthly Notices of the Royal Astronomical Society 457 (02, 2016) 4081-4088, [ https://academic.oup.com/mnras/article-pdf/457/4/4081/18516845/stw284.pdf ].\n- [68] H. Lazare, D. Sarkar and E. D. Kovetz, HERA bound on X-ray luminosity weakens when accounting for Population III stars , arXiv e-prints (July, 2023) arXiv:2307.15577, [ 2307.15577 ].\n- [69] A. Lidz, O. Zahn, M. McQuinn, M. Zaldarriaga and S. Dutta, Higher Order Contributions to the 21 cm Power Spectrum , Astrophys. J. 659 (2007) 865-876, [ astro-ph/0610054 ].  \n- [70] J. R. Pritchard and A. Loeb, 21-cm cosmology , Rept. Prog. Phys. 75 (2012) 086901, [ 1109.6012 ].\n- [71] A. Mesinger, A. Ewall-Wice and J. Hewitt, Reionization and beyond: detecting the peaks of the cosmological 21 cm signal , MNRAS 439 (Apr., 2014) 3262-3274, [ 1310.0465 ].\n- [72] A. Mesinger, A. Ewall-Wice and J. Hewitt, Reionization and beyond: detecting the peaks of the cosmological 21 cm signal , Mon. Not. Roy. Astron. Soc. 439 (2014) 3262-3274, [ 1310.0465 ].\n- [73] HERA collaboration, Z. Abdurashidova et al., HERA Phase I Limits on the Cosmic 21 cm Signal: Constraints on Astrophysics and Cosmology during the Epoch of Reionization , Astrophys. J. 924 (2022) 51, [ 2108.07282 ].\n- [74] T. R. Slatyer, Energy Injection And Absorption In The Cosmic Dark Ages , Phys. Rev. D87 (2013) 123513, [ 1211.0283 ].\n- [75] T. R. Slatyer, Indirect dark matter signatures in the cosmic dark ages. I. Generalizing the bound on s-wave dark matter annihilation from Planck results , Phys. Rev. D 93 (2016) 023527, [ 1506.03811 ].\n- [76] M. Fr'echet, Sur l'extension de certaines 'evaluations statistiques au cas de petits 'echantillons , Revue de l'Institut International de Statistique (1943) 182-205.\n- [77] G. Darmois, Sur les limites de la dispersion de certaines estimations , Revue de l'Institut International de Statistique (1945) 9-15.\n- [78] A. Aitken and H. Silverstone, Xv.-on the estimation of statistical parameters , Proceedings of the Royal Society of Edinburgh Section A: Mathematics 61 (1942) 186-194.\n- [79] C. A. Mason, J. B. Mu˜noz, B. Greig, A. Mesinger and J. Park, 21cmfish: Fisher-matrix framework for fast parameter forecasts from the cosmic 21-cm signal , 2212.09797 .\n- [80] O. Zahn, A. Mesinger, M. McQuinn, H. Trac, R. Cen and L. E. Hernquist, Comparison Of Reionization Models: Radiative Transfer Simulations And Approximate, Semi-Numeric Models , Mon. Not. Roy. Astron. Soc. 414 (2011) 727, [ 1003.3455 ].\n- [81] J. C. Pober, A. R. Parsons, D. R. DeBoer, P. McDonald, M. McQuinn, J. E. Aguirre et al., The Baryon Acoustic Oscillation Broadband and Broad-beam Array: Design Overview and Sensitivity Forecasts , AJ 145 (Mar., 2013) 65, [ 1210.2413 ].\n- [82] J. C. Pober, A. Liu, J. S. Dillon, J. E. Aguirre, J. D. Bowman, R. F. Bradley et al., What Next-generation 21 cm Power Spectrum Measurements can Teach us About the Epoch of Reionization , ApJ 782 (Feb., 2014) 66, [ 1310.7031 ].\n- [83] M. Boudaud, J. Lavalle and P. Salati, Novel cosmic-ray electron and positron constraints on MeV dark matter particles , Phys. Rev. Lett. 119 (2017) 021103, [ 1612.07698 ].\n- [84] M. Boudaud, T. Lacroix, M. Stref and J. Lavalle, Robust cosmic-ray constraints on p -wave annihilating MeV dark matter , Phys. Rev. D 99 (2019) 061302, [ 1810.01680 ].\n- [85] R. Essig, E. Kuflik, S. D. McDermott, T. Volansky and K. M. Zurek, Constraining Light Dark Matter with Diffuse X-Ray and Gamma-Ray Observations , JHEP 11 (2013) 193, [ 1309.4091 ].\n- [86] T. Cohen, K. Murase, N. L. Rodd, B. R. Safdi and Y. Soreq, γ -ray Constraints on Decaying Dark Matter and Implications for IceCube , Phys. Rev. Lett. 119 (2017) 021102, [ 1612.05638 ].\n- [87] A. Massari, E. Izaguirre, R. Essig, A. Albert, E. Bloom and G. A. G'omez-Vargas, Strong Optimized Conservative Fermi -LAT Constraints on Dark Matter Models from the Inclusive Photon Spectrum , Phys. Rev. D 91 (2015) 083539, [ 1503.07169 ].\n- [88] M. Cirelli, N. Fornengo, J. Koechler, E. Pinetti and B. M. Roach, Putting all the X in one basket: Updated X-ray constraints on sub-GeV Dark Matter , JCAP 07 (2023) 026, [ 2303.08854 ].  \n- [89] D. Cadamuro and J. Redondo, Cosmological bounds on pseudo Nambu-Goldstone bosons , JCAP 02 (2012) 032, [ 1110.2895 ].\n- [90] F. Calore, A. Dekker, P. D. Serpico and T. Siegert, Constraints on light decaying dark matter candidates from 16 yr of INTEGRAL/SPI observations , Mon. Not. Roy. Astron. Soc. 520 (2023) 4167-4172, [ 2209.06299 ].\n- [91] J. W. Foster, Y. Park, B. R. Safdi, Y. Soreq and W. L. Xu, Search for dark matter lines at the Galactic Center with 14 years of Fermi data , Phys. Rev. D 107 (2023) 103047, [ 2212.07435 ].\n- [92] D. Wadekar and Z. Wang, Strong constraints on decay and annihilation of dark matter from heating of gas-rich dwarf galaxies , Phys. Rev. D 106 (2022) 075007, [ 2111.08025 ].\n- [93] P. De la Torre Luque, S. Balaji and J. Koechler, Importance of Cosmic-Ray Propagation on Sub-GeV Dark Matter Constraints , Astrophys. J. 968 (2024) 46, [ 2311.04979 ].\n- [94] S. R. Furlanetto and S. J. Stoever, Secondary ionization and heating by fast electrons , MNRAS 404 (June, 2010) 1869-1878, [ 0910.4410 ]."}