Systematic uncertainties in models of the cosmic dawn
MMon. Not. R. Astron. Soc. , ?? – ?? (2020) Printed 15 December 2020 (MN L A TEX style file v2.2)
Systematic uncertainties in models of the cosmic dawn
Jordan Mirocha, (cid:63) † Henri Lamarre, Adrian Liu ‡ Department of Physics & McGill Space Institute, McGill University, 3600 Rue University, Montréal, QC, H3A 2T8
15 December 2020
ABSTRACT
Models of the reionization and reheating of the intergalactic medium (IGM) at red-shifts z (cid:38) continue to grow more sophisticated in anticipation of near-future datasetsfrom 21-cm, cosmic microwave background, and galaxy survey measurements. Thougha relatively simple picture of high- z galaxy evolution has emerged in recent years thatshould be well constrained by upcoming observations, there are many potential sourcesof systematic uncertainty in models that could bias and/or degrade these constraintsif left unaccounted for. In this work, we examine three commonly-ignored sources ofuncertainty in models for the mean reionization and thermal histories of the IGM:the underlying cosmology, the halo mass function (HMF), and the choice of stellarpopulation synthesis (SPS) model. We find that cosmological uncertainties are small,affecting the Thomson scattering optical depth at the few percent level and the ampli-tude of the global 21-cm signal at the ∼ few mK level. The differences brought aboutby choice of HMF and SPS models are more dramatic, comparable to the σ error-baron τ e and a ∼ mK effect on the global 21-cm signal amplitude. We compare modelswith comparable empirical predictions by jointly fitting galaxy luminosity functionsand global 21-cm signals, and find that (i) doing so requires additional free param-eters to compensate for modeling systematics and (ii) the spread in constraints onparameters of interest for different HMF and SPS choices, assuming mK noise in theglobal signal, is comparable to those obtained when adopting the “true” HMF and SPSwith (cid:38) mK errors. Our work highlights the need for dedicated efforts to reduce theuncertainties in common modeling ingredients in order to enable precision inferencewith future datasets. Key words: galaxies: high-redshift – intergalactic medium – galaxies: luminosityfunction, mass function – dark ages, reionization, first stars – diffuse radiation.
Many ongoing and near-future experiments are designed inlarge-part to provide astrophysical and cosmological con-straints on the “cosmic dawn,” when the first stars and galax-ies formed. For example, galaxy surveys searching for highredshift z (cid:38) galaxies will continue piecing together theluminosity function of galaxies over cosmic time (Williamset al., 2018; Behroozi et al., 2020), while cosmic microwavebackground (CMB) measurements constrain the rough tim-ing and duration of reionization from the Thomson scatter-ing optical depth, τ e , and kinetic Sunyaev-Z’eldovich effect(Planck Collaboration et al., 2018; Reichardt et al., 2020).21-cm experiments target neutral hydrogen itself, and arethus direct probes of the mean ionization and thermal his- (cid:63) [email protected] † CITA National Fellow ‡ [email protected] tories, as well as their topology (Furlanetto, 2006; Pritchard& Loeb, 2012).Of course, the native form of measurements providedby the aforementioned probes, e.g., luminosity functions, τ e ,and 21-cm spectra, are not straightforward to interpret. As aresult, forward models of galaxies and reionization are usedto predict these quantities, and, once coupled to samplingalgorithms like Markov Chain Monte Carlo (MCMC), canbe used to infer the underlying parameters of the modelfrom measurements. For example, UV luminosity functions(UVLFs) of high- z galaxies can constrain the star formationefficiency (SFE; f ∗ ) of galaxies (Mason et al., 2015; Sun &Furlanetto, 2016), while the reionization and thermal histo-ries are sensitive to, e.g., the product of f ∗ , the escape frac-tion of ionizing photons f esc , and assumptions about galax-ies beyond detection limits (Robertson et al., 2015; Bouwenset al., 2015a). Meanwhile, 21-cm probes are sensitive also tothe minimum mass of halos capable of forming stars, M min ,and the efficiency of X-ray and non-ionizing UV photon pro- © 2020 RAS a r X i v : . [ a s t r o - ph . C O ] D ec Mirocha, Lamarre, & Liu duction in high- z galaxies (Barkana & Loeb, 2005; Furlan-etto, 2006; Pritchard & Furlanetto, 2007; Mesinger et al.,2013). Additionally, 21-cm experiments can in principle beused to constrain fundamental physics, given that warm andinteracting dark matter models affect the formation of struc-ture on small scales (Sitwell et al., 2014; Schneider, 2018;Muñoz et al., 2020) and modulate the IGM thermal historyrelative to the predictions of Λ cold dark matter (CDM)models (Tashiro et al., 2014; Barkana, 2018; Muñoz & Loeb,2018; Boddy et al., 2018; Fialkov et al., 2018). Indeed, manyexperiments are underway, and have provided interestinglimits on both the reionization and thermal histories (e.g.,Bowman et al., 2018; Monsalve et al., 2017; Singh et al.,2018).Several parameter inference pipelines have emerged inrecent years to interpret the impending deluge of constraintson the cosmic dawn. Semi-empirical or semi-analytic modelsof galaxy formation are now commonly employed to connectconstraints on UVLFs, the mean reionization history, andionizing background after reionization (Finkelstein et al.,2019; Tacchella et al., 2018; Behroozi et al., 2019; Yunget al., 2019), as well as global 21-cm signals (Mirocha et al.,2017). Semi-numerical models of reionization have emergedto efficiently model the 21-cm background in 3-D (Mesinger& Furlanetto, 2007; Mesinger et al., 2011; Santos et al., 2010;Fialkov et al., 2013; Hutter et al., 2020), and have recentlybeen coupled to MCMC samplers (Greig & Mesinger, 2015,2017) and emulators (Kern et al., 2017; Cohen et al., 2020)to enable joint galaxy survey–21-cm inference (Park et al.,2019; Qin et al., 2020).Recent forecasts largely find that precision constraintson parameters of interest are well within the reach of cur-rent facilities. For example, 21-cm measurements can con-strain the population-averaged values of f ∗ , f esc , above M min (Mirocha et al., 2015; Greig & Mesinger, 2017), warm darkmatter (Muñoz et al., 2020), and even cosmological param-eters (Liu et al., 2016; Kern et al., 2017). If limited onlyby thermal noise, current experiments are expected to reachthe cosmic variance floor in relatively short order (Muñoz& Cyr-Racine, 2020). Of course, such forecasts are based onideal instruments, foregrounds, and noise—luxuries yet tobe achieved in real experiments, for which systematic un-certainties are the main impediment to progress (see, e.g.,Liu & Shaw 2020 for a description of some of the main is-sues).Systematic uncertainties are not limited to the realmof experiment. Theoretical models are also uncertain, inpart due to numerical approximations made to acceleratecomputations, but also because of assumptions made in themodeling. Some of these assumptions are not necessarily pa-rameterizable, and thus may be difficult to marginalize overwhen computing parameter constraints. For example, thesemi-numerical approach to reionization is an approxima-tion employed to avoid running expensive radiative transfersimulations, and appears to work at the level of ∼ − %(Zahn et al., 2011; Majumdar et al., 2014; Hutter, 2018). Itis not clear that the level of uncertainty is invariant withrespect to model parameters or implementations, and so itis not clear how to generically inflate uncertainties on modelparameters of interest. Other aspects of cosmic dawn model-ing may be even more difficult to pin down as they are unre- lated to numerical approximations, e.g., choosing the “right”stellar population synthesis model or halo mass function.There is no doubt that observational issues are still themost pressing. However, given that the sensitivity of the cur-rent generation of 21-cm experiments is sufficient for detec-tion of the global 21-cm signal and power spectrum, a break-through in the reduction of systematics could lead to rapidtransition from crude upper limits to proper constraints onastrophysical and cosmological parameters. As a result, itis fruitful to consider the limitations of current models, sothat particularly weak areas of models can be identified andtargeted for priority work. We also hope that doing so mayalso inspire the community to critically examine the degreeof precision needed for different science questions.While there are many potential sources of system-atic uncertainty in models, in this work we focus on threemain effects associated with seemingly innocuous model-ing decisions: (i) the neglect of uncertainties in cosmolog-ical parameters, (ii) the halo mass function (HMF), and(iii) stellar population synthesis (SPS) codes. While cos-mological parameter variations can in principle be includedself-consistently and marginalized over, they rarely are,especially in semi-numeric codes for which re-generatingthe cosmology-dependent density and velocity fields canbe the most expensive part of the model. Unlike the cos-mological parameters, the choice of HMF and SPS modeldoes not lie along a continuum of possibilities. In otherwords, these “parameters” are switches rather than knobs.Other likely causes of uncertainty include numerical tech-niques (e.g., two-zone vs. semi-numeric vs. fully-numerical),or choice of parameterization for galaxy properties (e.g.,mass-dependent vs. mass-independent parameters, inclusionof multiple sources populations). These areas are perhapsequally deserving of attention, though we defer them to fu-ture work.The structure of the paper is as follows. In §2 we reviewour modeling approach, giving special attention to the com-ponents most likely to be uncertain. We present our mainresults in §3 and conclude in §4. Our approach to modeling the reionization and re-heatingof the IGM closely follows early two-zone models for theIGM (Furlanetto, 2006; Pritchard & Loeb, 2010) though weuse updated models for the X-ray background and galax-ies following Mirocha (2014) and Mirocha et al. (2017), re-spectively. We provide a brief summary of the model be-low, drawing attention to the components sucseptible to themodeling uncertainties we investigate. Our results can bere-created using the publicly-available ares code .To begin, we assume that galaxies inhabit dark matterhalos in a 1:1 fashion and that star formation is fueld by theinflow of pristine gas from the IGM. This means that thestar formation rate ˙ M ∗ (SFR) in galaxies is directly relatedto the mass accretion rate ˙ M h (MAR) of dark matter halos, ˙ M ∗ = f ∗ f b ˙ M h (1) https://ares.readthedocs.io/en/latest/ © 2020 RAS, MNRAS , ?? ––
Many ongoing and near-future experiments are designed inlarge-part to provide astrophysical and cosmological con-straints on the “cosmic dawn,” when the first stars and galax-ies formed. For example, galaxy surveys searching for highredshift z (cid:38) galaxies will continue piecing together theluminosity function of galaxies over cosmic time (Williamset al., 2018; Behroozi et al., 2020), while cosmic microwavebackground (CMB) measurements constrain the rough tim-ing and duration of reionization from the Thomson scatter-ing optical depth, τ e , and kinetic Sunyaev-Z’eldovich effect(Planck Collaboration et al., 2018; Reichardt et al., 2020).21-cm experiments target neutral hydrogen itself, and arethus direct probes of the mean ionization and thermal his- (cid:63) [email protected] † CITA National Fellow ‡ [email protected] tories, as well as their topology (Furlanetto, 2006; Pritchard& Loeb, 2012).Of course, the native form of measurements providedby the aforementioned probes, e.g., luminosity functions, τ e ,and 21-cm spectra, are not straightforward to interpret. As aresult, forward models of galaxies and reionization are usedto predict these quantities, and, once coupled to samplingalgorithms like Markov Chain Monte Carlo (MCMC), canbe used to infer the underlying parameters of the modelfrom measurements. For example, UV luminosity functions(UVLFs) of high- z galaxies can constrain the star formationefficiency (SFE; f ∗ ) of galaxies (Mason et al., 2015; Sun &Furlanetto, 2016), while the reionization and thermal histo-ries are sensitive to, e.g., the product of f ∗ , the escape frac-tion of ionizing photons f esc , and assumptions about galax-ies beyond detection limits (Robertson et al., 2015; Bouwenset al., 2015a). Meanwhile, 21-cm probes are sensitive also tothe minimum mass of halos capable of forming stars, M min ,and the efficiency of X-ray and non-ionizing UV photon pro- © 2020 RAS a r X i v : . [ a s t r o - ph . C O ] D ec Mirocha, Lamarre, & Liu duction in high- z galaxies (Barkana & Loeb, 2005; Furlan-etto, 2006; Pritchard & Furlanetto, 2007; Mesinger et al.,2013). Additionally, 21-cm experiments can in principle beused to constrain fundamental physics, given that warm andinteracting dark matter models affect the formation of struc-ture on small scales (Sitwell et al., 2014; Schneider, 2018;Muñoz et al., 2020) and modulate the IGM thermal historyrelative to the predictions of Λ cold dark matter (CDM)models (Tashiro et al., 2014; Barkana, 2018; Muñoz & Loeb,2018; Boddy et al., 2018; Fialkov et al., 2018). Indeed, manyexperiments are underway, and have provided interestinglimits on both the reionization and thermal histories (e.g.,Bowman et al., 2018; Monsalve et al., 2017; Singh et al.,2018).Several parameter inference pipelines have emerged inrecent years to interpret the impending deluge of constraintson the cosmic dawn. Semi-empirical or semi-analytic modelsof galaxy formation are now commonly employed to connectconstraints on UVLFs, the mean reionization history, andionizing background after reionization (Finkelstein et al.,2019; Tacchella et al., 2018; Behroozi et al., 2019; Yunget al., 2019), as well as global 21-cm signals (Mirocha et al.,2017). Semi-numerical models of reionization have emergedto efficiently model the 21-cm background in 3-D (Mesinger& Furlanetto, 2007; Mesinger et al., 2011; Santos et al., 2010;Fialkov et al., 2013; Hutter et al., 2020), and have recentlybeen coupled to MCMC samplers (Greig & Mesinger, 2015,2017) and emulators (Kern et al., 2017; Cohen et al., 2020)to enable joint galaxy survey–21-cm inference (Park et al.,2019; Qin et al., 2020).Recent forecasts largely find that precision constraintson parameters of interest are well within the reach of cur-rent facilities. For example, 21-cm measurements can con-strain the population-averaged values of f ∗ , f esc , above M min (Mirocha et al., 2015; Greig & Mesinger, 2017), warm darkmatter (Muñoz et al., 2020), and even cosmological param-eters (Liu et al., 2016; Kern et al., 2017). If limited onlyby thermal noise, current experiments are expected to reachthe cosmic variance floor in relatively short order (Muñoz& Cyr-Racine, 2020). Of course, such forecasts are based onideal instruments, foregrounds, and noise—luxuries yet tobe achieved in real experiments, for which systematic un-certainties are the main impediment to progress (see, e.g.,Liu & Shaw 2020 for a description of some of the main is-sues).Systematic uncertainties are not limited to the realmof experiment. Theoretical models are also uncertain, inpart due to numerical approximations made to acceleratecomputations, but also because of assumptions made in themodeling. Some of these assumptions are not necessarily pa-rameterizable, and thus may be difficult to marginalize overwhen computing parameter constraints. For example, thesemi-numerical approach to reionization is an approxima-tion employed to avoid running expensive radiative transfersimulations, and appears to work at the level of ∼ − %(Zahn et al., 2011; Majumdar et al., 2014; Hutter, 2018). Itis not clear that the level of uncertainty is invariant withrespect to model parameters or implementations, and so itis not clear how to generically inflate uncertainties on modelparameters of interest. Other aspects of cosmic dawn model-ing may be even more difficult to pin down as they are unre- lated to numerical approximations, e.g., choosing the “right”stellar population synthesis model or halo mass function.There is no doubt that observational issues are still themost pressing. However, given that the sensitivity of the cur-rent generation of 21-cm experiments is sufficient for detec-tion of the global 21-cm signal and power spectrum, a break-through in the reduction of systematics could lead to rapidtransition from crude upper limits to proper constraints onastrophysical and cosmological parameters. As a result, itis fruitful to consider the limitations of current models, sothat particularly weak areas of models can be identified andtargeted for priority work. We also hope that doing so mayalso inspire the community to critically examine the degreeof precision needed for different science questions.While there are many potential sources of system-atic uncertainty in models, in this work we focus on threemain effects associated with seemingly innocuous model-ing decisions: (i) the neglect of uncertainties in cosmolog-ical parameters, (ii) the halo mass function (HMF), and(iii) stellar population synthesis (SPS) codes. While cos-mological parameter variations can in principle be includedself-consistently and marginalized over, they rarely are,especially in semi-numeric codes for which re-generatingthe cosmology-dependent density and velocity fields canbe the most expensive part of the model. Unlike the cos-mological parameters, the choice of HMF and SPS modeldoes not lie along a continuum of possibilities. In otherwords, these “parameters” are switches rather than knobs.Other likely causes of uncertainty include numerical tech-niques (e.g., two-zone vs. semi-numeric vs. fully-numerical),or choice of parameterization for galaxy properties (e.g.,mass-dependent vs. mass-independent parameters, inclusionof multiple sources populations). These areas are perhapsequally deserving of attention, though we defer them to fu-ture work.The structure of the paper is as follows. In §2 we reviewour modeling approach, giving special attention to the com-ponents most likely to be uncertain. We present our mainresults in §3 and conclude in §4. Our approach to modeling the reionization and re-heatingof the IGM closely follows early two-zone models for theIGM (Furlanetto, 2006; Pritchard & Loeb, 2010) though weuse updated models for the X-ray background and galax-ies following Mirocha (2014) and Mirocha et al. (2017), re-spectively. We provide a brief summary of the model be-low, drawing attention to the components sucseptible to themodeling uncertainties we investigate. Our results can bere-created using the publicly-available ares code .To begin, we assume that galaxies inhabit dark matterhalos in a 1:1 fashion and that star formation is fueld by theinflow of pristine gas from the IGM. This means that thestar formation rate ˙ M ∗ (SFR) in galaxies is directly relatedto the mass accretion rate ˙ M h (MAR) of dark matter halos, ˙ M ∗ = f ∗ f b ˙ M h (1) https://ares.readthedocs.io/en/latest/ © 2020 RAS, MNRAS , ?? –– ?? osmic dawn modeling systematics were f ∗ is the star formation efficiency (SFE) and f b isthe cosmic baryon fraction. We assume that the SFE as adouble-power law (DPL) in halo mass, f ∗ ( M h ) = f ∗ , C (cid:16) M h M p (cid:17) − α lo + (cid:16) M h M p (cid:17) − α hi (2)where f ∗ , is the SFE at M (cid:12) , M p is the mass at which f ∗ peaks, and α hi and α lo describe the power-law index atmasses above and below the peak, respectively. The addi-tional constant C ≡ (10 /M p ) − α lo + (10 /M p ) − α hi isintroduced to re-normalize the standard DPL formula to M (cid:12) , rather than the peak mass. We model the MARfollowing Furlanetto et al. (2017), and assume halos growat fixed number density, in which case ˙ M h can be inferredfrom the evolution in the HMF. Though idealized, it hasthe advantage of preserving self-consistency, and does ap-pear to agree well with the results of N-body simulations(Trac et al., 2015), at least for relatively massive halos at z (cid:46) . The MAR is yet another potential source of mod-eling uncertainty that may be comparable to the effects weexplore here (Schneider et al., 2020), though we defer inves-tigation of this possibility to future work.In this work, we only consider the rest-ultrioviolet con-tinuum emission of galaxies ( Å), which is probed byrecent observations with
Hubble (Bouwens et al., 2015b;Finkelstein et al., 2015) targeting galaxies at z (cid:38) . Be-cause the rest-UV continuum probes young stars, it is agood tracer of recent star formation and thus the SFR ofgalaxies. Under the assumption of a constant SFR, the UVluminosity of a galaxy will asymptote to a constant valueon timescales of order ∼ Myr. In this limit, the UVluminosity can be modeled simply as L UV = l UV ˙ M ∗ (3)where l UV carries units of erg s − Hz − ( M (cid:12) yr − ) − andis often written instead via its reciprocal κ UV ≡ l − .Finally, with these assumptions, the UV luminosityfunction is simply dφ ( L UV ) = dn h dM h dM h dL UV dL UV . (4)Here, dn h /dM h is the halo mass function, which we com-pute using the hmf code (Murray et al., 2013a), which re-lies on camb for the matter power spectrum (Lewis et al.,2000). For an assumed value of l UV , one can calibrate f ∗ empirically by fitting UVLF constraints from any numberof studies (Bouwens et al., 2015b; Finkelstein et al., 2015)using Eqs. (1)-(3) to compute dM h /dL UV . We note thatseveral aspects of this model have since been generalized,e.g., self-consistent dust reddening, “bursty star” formationhistories (Mirocha et al., 2020; Mirocha, 2020), though weneglect these complications for the duration of this paper, asthey will not obviously exacerbate (or mitigate) the sourcesof uncertainty we explore here. For all fits presented in theremainder of this paper, we assume a standard dust cor-rection using the Meurer et al. (1999) IRX- β relation and M UV - β relations from Bouwens et al. (2014), and use em-cee (Foreman-Mackey et al., 2013) to explore the multi-dimensional parameter space.With a model for sources in hand, we can proceed tomodeling the mean reionization history and global 21-cm signal. The differential brightness temperature δT b is givenby δT b (cid:39) − x H ii ) (cid:18) Ω b , h . (cid:19) (cid:18) . m , h z (cid:19) / (cid:18) − T γ T S (cid:19) , (5)where x H ii = Q HII + (1 − Q HII ) x e is the volume-averagedionized fraction, and T − S ≈ T − γ + x c T − K + x α T − α x c + x α . (6)is the spin temperature. The spin temperature quantifies thelevel populations in the hyperfine singlet and triplet states,and is set by a combination of collisional coupling (quan-tified by x c ; Zygelman, 2005), radiative coupling throughthe Wouthuysen-Field effect (quantified by x α ; Wouthuy-sen, 1952; Field, 1958) and microwave background temper-ature, T γ . Equation (5) assumes the cosmic mean densityand neglects relative velocities along the line of sight, as isappropriate for the global 21-cm signal.Our model partitions the IGM into two distinct phases:an ionized phase, characterized by its volume-filling factor Q HII , and a “bulk IGM” phase, characterized by its tem-perature and Ly- α intensity. To evolve the gas temperaturein the bulk IGM phase, we solve the uniform backgroundproblem as outlined in Mirocha (2014), essentially a high- z implementation of the Haardt & Madau (1996) algorithm.Assuming a multi-colour disk spectrum for X-ray emission(Mitsuda et al., 1984) with no host-galaxy absorption, wesolve for the mean intensity of the cosmic X-ray backgroundas a function of redshift and photon energy, J ν ( z ) , and in-tegrate over J ν (weighted by bound-free absorption cross-sections) to obtain the photo-heating and ionization ratesas a function of redshift. We assume a fully-neutral IGM,in which case the bound-free optical depth of the IGM toX-ray photons can be tabulated ahead of time for an ef-ficiency boost. The intensity of the Ly- α background, J α ,which sets the strength of Wouthuysen-Field coupling, iscomputed in an analogous fashion. Note that this part ofthe model is cosmology-dependent, as J ν depends on theHubble parameter—we self-consistently include the cosmol-ogy dependence in this computation. See §2.2-2.3 in Mirocha(2014) for more details.Because the mean-free path of ultraviolet photons isshort, we relate the growth rate of the fractional ionizedvolume, Q HII , directly to the ionizing photon productionrate, i.e., dQ HII dt = f esc N ion ˙ ρ ∗ n H − α H ii Cn e Q HII (7)where f esc is the escape fraction of ionizing photons, N ion isthe number of ionizing photons emitted per unit stellar mass,and ˙ ρ ∗ ≡ (cid:82) dM h ˙ M ∗ dn h /dM h is the cosmic star formationrate density. The recombination rate is the product of theclumping factor C , which we set to unity for simplicity, theelectron number density n e , and the case A recombinationcoefficient, α H ii . Having established our basic underlying model of UVLFsand the global cm signal, we now highlight several poten-tial sources of theoretical uncertainty. © 2020 RAS, MNRAS , ?? – ?? Mirocha, Lamarre, & Liu
Latent in virtually every expression in this sectionare assumptions about cosmological parameters. The 21-cmbackground itself depends directly on Ω b, and the Hubbleparameter, as it involves an integral along the line of sight,but also indirectly, as the opacity of the IGM to X-rays sim-ilarly depends on Ω b, , the primordial helium abundance,and the Hubble parameter.The HMF, while possible to compute analytically underidealized circumstances, is likely more complicated in real-ity. Recent numerical simulations (e.g., Tinker et al., 2010)suggest that the HMF lies somewhere between the classic an-alytic models (Press & Schechter, 1974; Sheth et al., 2001)often used in reionization models. The results of N-bodysimulations can be reproduced via fitting functions, thoughthese vary from model to model due to differences in gravitysolver and/or halo finder used in the analysis (Knebe et al.,2013a,b), and themselves may be insufficient for precisionwork (Bhattacharya et al., 2011). Most numerical simula-tions run to date are not optimized for high- z (though see,e.g., Reed et al., 2007), nor has there been a systematicexploration of the agreement between different fitting func-tions at z (cid:38) . The fact that variations from one HMFfitting function to the next are larger than those caused bycosmological parameter variations for a single HMF at low- z (Murray et al., 2013b) suggests this is a problem worthinvestigating at high- z as well. In this work, we investigatethe Press & Schechter (1974), Sheth et al. (2001), and theTinker et al. (2010) form of the mass function, the latter ofwhich is a representative example that lies between. Otherrecent results are comparable to Tinker et al. (2010) (War-ren et al., 2006; Reed et al., 2007; Bhattacharya et al., 2011),at least at low redshift where comparisons are most common(Murray et al., 2013b).Next, the conversion factor from SFR to luminosity, κ UV (see Eq. 3), relies on assumptions about stellar evolu-tion, atmospheres, and initial mass function (IMF). Like theHMF, pieces of κ UV are in principle “knobs” just like the cos-mological parameters, e.g., the stellar initial mass function(IMF), which is generally taken to be a power-law or brokenpower-law. However, in practice, it is common to choose anSPS code and implicitly adopt the IMF used within it. Inthis work we use the bpass version 1.0 models (Eldridge &Stanway, 2009) as well as the original starburst99 models(Leitherer et al., 1999).In Figure 1 we show the HMFs adopted in this work.We show both the cumulative HMF (top) and fraction ofmatter in collapsed halos, f coll , above the atomic coolingthreshold [ M h ( T vir = 10 K) ≡ m ; bottom]. At z = 6 ,the disagreement between models is most clear in massivehalos, M h (cid:38) M (cid:12) , though upon close inspection differ-ences are present in low-mass halos, and the sense of thedisparity between fitting functions reverses. At z (cid:38) , thedisagreement is noticeable at all masses. The bottom panelis perhaps the clearer metric for reionization models, wherethe production rate of UV and X-ray photons is often linkeddirectly to the rate at which mass collapses into halos. At z ∼ , the spread is nearly an order of magnitude, whileat z ∼ it is still a factor of ∼ . In all cases, variations inthe cosmological parameters (shaded regions) has a smallereffect than the choice of fitting function (different colours).In Figure 2, we compare the SPS models used in thiswork. First, in the top panel, we show an example spectral M h / M n ( > m ) PSSTTinker10 z = 6, 10, 156 8 10 12 14 16 18 20 z f c o ll ( m > m ) Figure 1. Difference in cumulative HMF (top) and col-lapse fraction (bottom) for different fitting functions.
Width of each line is due to variations in cosmological param-eters ( and σ contours drawn for each), while each colour indi-cates a different fitting function: Press & Schechter (1974) (PS;cyan), Sheth et al. (2001) (ST; blue), and Tinker et al. (2010)(Tinker10; magenta). In the top panel, the three sets of curvesindicate redshifts z = 6 , , and 15 from top to bottom. energy distribution (SED) for a constant SFR and metallic-ity of Z = 0 . using both bpass (black) and starburst99 (blue). There is a clear offset in normalization between thetwo codes, with starburst99 producing more overall emis-sion (per M (cid:12) yr − ) at all wavelengths. To investigate the λ -dependent nature of the offset, we scale the starburst99 spectrum so that it matches bpass at Å (blue dashed).It is now clear that, in addition to producing a stronger non-ionizing continuum, starburst99 also has a higher ratio ofnon-ionizing to ionizing emission than bpass .The SED effects are explored further in the bottompanel of Fig. 2. For a series of metallicities, we computethe ratio of ionizing (black) and Lyman-Werner band (blue)luminosities to the
Å luminosity for each code, and thenplot the ratio of values for each code, i.e., R ≡ ( N × κ UV ) bpass ( N × κ UV ) s99 (8) © 2020 RAS, MNRAS , ?? – ?? osmic dawn modeling systematics / Å L [ e r g s Å ] = 0.004 bpass-sins99s99 (scaled)0.005 0.010 0.015 0.0201.01.21.41.61.82.02.22.42.6 R bp a ss / s N ion UV N LW UV bpass-sinbpass-bin
Figure 2. Difference in rest-UV SEDs of galaxies as com-puted with bpass (black) vs. starburst99 . Top panelshows the SED of a metallicity Z = 0 . stellar populationforming stars continuously for the last 100 Myr, while the bottompanel compares the ratio of ionizing ( N ion ) and Lyman-Werner( N LW ) luminosities to the Å continuum for each code at aseries of stellar metallicities. where N is either N ion , the number of ionizing photons emit-ted per stellar baryon, or N LW , the equivalent quantity inthe Lyman-Werner band.The rationale here is that the rest-ultraviolet ∼ Åemission is the only part of the spectrum that is (generally)observable—at fixed κ UV , UVLFs constrain f ∗ , so any dif-ferences in mean reionization history and global 21-cm sig-nal will largely be determined by the product N ion κ UV and N LW κ UV , respectively. Indeed, such differences can be manytens of percent, perhaps a factor of ∼ in extreme cases,particularly at low metallicity Z (cid:46) . and when adopt-ing the bpass binaries models (“bpass-bin”; open symbols).We focus only on single star models throughout this workin order to make a more fair comparison between bpass and starburst99 .
24 22 20 18 16 14 M UV ( M U V )[ m a g c M p c ] z z Figure 3. Difference in UVLFs caused by cosmology andHMF variations.
Linestyle and colour conventions are the sameas Figures 4 and 6. Data at z ∼ are from Bouwens et al. (2015b)(circles) and Finkelstein et al. (2015) (squares), while z ∼ measurements are from Oesch et al. (2018). Note that there arethree bands shown for each redshift. We now present our main results, consisting of mod-els run with elements from the Planck chain( planck_TTTEEE_lowl_lowE ), with three different HMFs, us-ing both starburst99 and bpass
SPS models. For eachcosmological model, we re-generate (i) the recombinationhistory using cosmorec (Chluba & Thomas, 2011), (ii) theoptical depth of the IGM used to compute the X-ray back-ground, and (iii) the HMF using hmf (Murray et al., 2013a).The cosmological model number is subsequently used as afree parameter, which indicates the appropriate lookup ta-ble to use for each of these three quantities, as well as thecosmological parameter values used elsewhere in the model.We begin with a pure forward modeling approach to ex-amine the extent to which these different assumptions affectpredictions for the mean reionization history, global 21-cmsignal, and UVLFs. Then, we turn our attention to the infer-ence problem, and quantify the degree to which parameterconstraints are worsened when including these “new” sourcesof uncertainty, while holding roughly fixed the UVLFs andglobal 21-cm signal via joint fitting. In the former case, weadopt the Mirocha et al. (2017) model calibration, whichfit the Bouwens et al. (2015b) UVLFs with the model de-scribed in §2 under the assumption of no dust reddening (forillustrative purposes). For the inference problem, we employthe standard IRX- β -based dust correction to more closelymimic the procedure likely to be applied on upcoming mea-surements.In Figure 3, we show the UVLF predictions subjectedto cosmological and HMF uncertainties. At z ∼ , the dif-ferences between Sheth et al. (2001) and Tinker et al. (2010)mass functions is comparable to the difference betweenUVLF measurements reported in Bouwens et al. (2015b) andFinkelstein et al. (2015), though we caution the reader thatthese two sets of UVLF measurements are obtained via dif- © 2020 RAS, MNRAS , ?? – ?? Mirocha, Lamarre, & Liu
50 100 150 200 /MHz T b / m K PSSTTinker10 68101215203080 z Figure 4. Global 21-cm signal predictions subject to cos-mological and HMF uncertainties.
The width of each bandis set by uncertainties in the cosmological parameters, while eachcolour indicates a different HMF fitting function. Open contoursuse bpass models, while filled contours use starburst99 . ferent means and, e.g., make slightly different assumptionsabout the wavelength corresponding to M UV .In Figure 4, we show predictions for the global 21-cmsignal subject to variations in the cosmological parametersand HMF fitting functions. Qualitatively, the curves are verysimilar—all exhibiting a deep absorption trough at frequen-cies of ν (cid:39) MHz, in line with the predictions of Mirochaet al. (2017). However, the level of variation is not insignif-icant, as we show more quantiatively in Figure 5. Here, weshow only the position of the absorption minimum ( ν min , δT b, min ) , and the Thomson scattering optical depth, τ e . Theposition in frequency varies by ∼ − MHz between thedifferent fitting functions, while cosmological uncertaintiesconstitute ∼ − MHz of these differences for each HMF.Similarly, the amplitude of the absorption trough is uncer-tain at the level of ∼ − mK, about ∼ mK of which isdue solely to uncertainties in cosmological parameters. Fi-nally, τ e uncertainties are ∼ . .In Figure 6, we show the full reionization history foreach model. Once again, differences are apparent by eye.For example, at the midpoint of reionization, uncertaintiesin the cosmological parameters contribute an uncertainty of ∼ % in the mean neutral fraction, while the position of thereionization midpoint itself varies by ∆ z rei (cid:39) . from Shethet al. (2001) to Tinker et al. (2010) mass functions. At firstglance, the latter seems to produce a reionization history in-compatible with current constraints. However, we emphasizethat the escape fraction is held fixed in all of these models.Increasing f esc from 0.2 to ∼ . − . is enough to bring theTinker et al. (2010) models into agreement with the Press &Schechter (1974) and Sheth et al. (2001) models. Note thatour neglect of dust for this exercise is in part responsiblefor needing such high escape fractions, since dustier galax-ies require more star formation to preserve agreement withUVLFs. We include dust following standard recipes (see §2)in all that follows.
100 105 110 115 120 125 max /MHz . . . . . e - - - - - - - T b ,max /mK .
04 0 .
045 0 .
05 0 .
055 0 . e - - - - - - - T b , m a x / m K P D F PSSTTinker10
Figure 5. Posterior PDFs of the position of the global21-cm absorption trough ( ν min , δT b, min ) and Thomsonoptical depth. Different HMF fitting functions are indicatedby colour, while open (filled) contours indicate the use of bpass ( starburst99 ). The width of the 68% and 95% confidence inter-vals for each individual HMF/SPS pairing are set by uncertaintiesin the cosmological parameters. z N e u t r a l F r a c t i o n , x H I PSSTTinker100.02 0.04 0.06 0.08 e Figure 6. Mean reionization history predictions subjectto cosmological and HMF uncertainties.
Conventions arethe same as Fig. 4. The inset shows the Thomson scattering opti-cal depth shown for each model, with gray regions indicating therange of values excluded at σ according to the Planck constraint τ e = 0 . ± . . © 2020 RAS, MNRAS , ?? – ?? osmic dawn modeling systematics We now turn our attention to the effect these uncer-tainties have on our ability to constrain the main parame-ters of our model. To proceed, we create a global 21-cm sig-nal mock generated assuming the Sheth et al. (2001) massfunction, bpass
SPS models, and the best-fit Planck cos-mology, with f ∗ parameters calibrated to the Bouwens et al.(2015b) UVLFs assuming the Meurer et al. (1999) IRX- β and Bouwens et al. (2014) M UV - β relations to correct fordust. We then add spectrally-flat Gaussian random noisewith standard deviation σ = 5 mK to the global 21-cm mock,which is comparable to the level of uncertainty caused bycosmological parameter variations. We adopt the Bouwenset al. (2015b) UVLFs throughout as our fiducial set of obser-vations to be fit, though similar conclusions would be drawnif we instead adopted the Finkelstein et al. (2015) measure-ments. We do not include the z ∼ UVLFs of Oesch et al.(2018) in the analysis.There is no guarantee that the global 21-cm signal and,e.g., the z ∼ − UVLFs from Bouwens et al. (2015a)can both be well-fit simultaneously if we adopt an HMF orSPS model different from that assumed when generating the21-cm mock. In order to accommodate this possibility, weemploy a simple extension to the standard double power-law f ∗ model, multiplying Eq. (2) by a modulation factor (as inSchneider et al., 2020), i.e., f ∗ → f ∗ (cid:20) (cid:18) M h a (cid:19) b (cid:21) c (9)which allows f ∗ to depart from a power-law at the low-massend. We further allow f ∗ , to evolve with redshift as apower-law with index γ z , resulting in a total of four newparameters. In this work, we treat these as nuisance param-eters, though of course departures from power-law and/orredshift-independent f ∗ likely correspond to very interest-ing scenarios worth constraining. Finally, we allow the fourusual SFE parameters to vary freely, as well as the normal-ization of the X-ray luminosity–SFR relation ( L X / SFR ), es-cape fraction of ionizing radiation ( f esc ), and minimum virialtemperature of star-forming halos, T min . These last three pa-rameters are only constrained by the global 21-cm signal inour framework, as we do not include measurements of thereionization history in our likelihood.The results of this exercise are shown in Figures 7 and8. First, in each panel of Fig. 7, we see six posterior dis-tributions: one for each combination of HMF (colours) andSPS model (open vs. filled contours). It is immediately clearthat the differences from posterior to posterior are gener-ally larger than the width of any individual posterior dis-tribution. In other words, these systematic modeling uncer-tainties are larger than the uncertainty caused by the er-ror budget of each individual fit (5 mK Gaussian randomnoise in global 21-cm signal and errors in UVLF reportedby Bouwens et al. 2015b) as well as uncertainty due tomarginalizing over cosmological parameters. For reference,we show also the 95% contours obtained under the assump-tion of 20 mK level noise in the global 21-cm signal as dottedcontours in each interior panel of Fig. 7. In some dimensions,the dotted contours enclose a subset of the individual PDFs,indicating that systematic modeling uncertainties are com-parable to those caused by a 20 mK noise level. However,generally speaking, the biases causd by modeling systemat- ics are not mimicked well by inflated uncertainties in theglobal 21-cm signal.Though Figure 7 paints a somewhat grim picture inwhich many parameters are constrained much more poorlythan idealized forecasts would suggest, a few panels givereasons to be optimistic. For example, the low-mass slope ofthe SFE (third column), α lo , while affected by systematics ina relative sense differs little from case to case in an absolutesense—the overall spread is only ∼ . . For reference, thedifferences in α lo for different stellar feedback scenarios is / (see, e.g., Dayal et al., 2013; Furlanetto et al., 2017),while purely empirical models for the SFE vary by δα lo ∼ . depending on assumptions about, e.g., dust (Mirochaet al., 2017; Tacchella et al., 2018; Mirocha et al., 2020).As a result, it seems that the part of the SFE most easilyrelated to physical arguments is imminently constrainable.Similarly, the L X -SFR normalization (sixth column) isalso quite impervious to the systematic uncertainties we ex-plored here—though bimodal with respect to SPS model,the difference is order ∼ . dex, which is hardly a causefor concern in differentiating models. More worrisome is theminimum virial temperature, T min , which varies by an orderof magnitude from end to end, and the high-mass behaviorof the SFE (set by M peak and α hi ). However, if we can con-fidently discard the Press & Schechter (1974) mass function(cyan contours), the range of viable possibilities in each slicethrough parameter space is substantially reduced.Figure 8 focuses on the ‘nuisance’ parameters intro-duced to alleviate tension caused by modeling systematics.In reality, these parameters are of interest, as they quan-tify the degree to which the efficiency of star formation de-parts from a double power-law in mass and whether it mayevolve with redshift. Both effects are expected to some de-gree, as stellar feedback models predict redshift evolution in f ∗ , while reionization feedback and/or efficient Pop III starformation are likely to cause a departure from the doublepower-law form at low-mass. In this work, we investigated three potential sources of sys-tematic uncertainty in models for high- z galaxies and themean ionization and thermal history of the IGM.First, by drawing samples from the Planck
MCMCchains, we showed that cosmological uncertainties alone re-sult in ∼ % errors on the mean ionized fraction of the IGM.The global 21-cm signal—which encodes the IGM thermaland ionization histories—is affected by cosmological uncer-tainties at the level of ∼ mK in the amplitude of its ab-sorption trough, with ∼ − MHz uncertainties in posi-tion. These cosmological uncertainties are the smallest ofthe three effects we explored.Next, we adopted three commonly-used halo mass func-tion models, Press & Schechter (1974), Sheth et al. (2001),and Tinker et al. (2010), as well as two commonly-used stel-lar population synthesis models, starburst99 (Leithereret al., 1999) and bpass (Eldridge & Stanway, 2009), re-sulting in a total of 6 possible pairings. Differences due toone’s choice of HMF and SPS model are larger than thosecaused by cosmological uncertainties by virtually any metric.For example, variations in the Thomson scattering optical © 2020 RAS, MNRAS , ?? – ?? Mirocha, Lamarre, & Liu - . - . - . - . - . log f , 10 . . . . l o g T m i n / K . . . log M peak / M . . . lo - . - . - . . hi . . f esc . . . . log L X /SFR . . . . . log T min /K .
23 9 .
43 9 .
63 9 . l o g L X / S F R . . f e s c - . - . - . . h i . . . l o .
31 1 .
71 2 . l o g M p e a k / M P D F PSSTTinker10
Figure 7. Posterior distributions for six different fits to the same mock global 21-cm signal with mK Gaussian noiseand the z ∼ UVLF from Bouwens et al. (2015b).
Colours indicate the different HMFs, while filled (open) contours use of starburst99 ( bpass ) SPS models. Marginalized 1-D distributions lie along the diagonal, where dashed lines correspond to the opencontours of the interior panels. Dotted contours in the interior panels indicate 95% contours obtained when assuming 20 mK, rather than5 mK, Gaussian random noise in the global 21-cm signal mock. depth, τ e , fill roughly half of the σ error bar reported by Planck , τ e = 0 . ± . (Fig. 6). The strongest feature ofthe global signal shifts by up to ∼ MHz in frequency and ∼ mK in amplitude (Fig. 4).To explore how these uncertainties propagate to con-straints on model parameters, we performed a joint MCMCfit of a mock global 21-cm signal with currenty high- z UVLFconstraints. Global 21-cm forecasts based on thermal noisegenerally predict tight, (cid:46) % level constraints on modelparameters of interest (see, e.g., Mirocha et al., 2015); how- ever, in many dimensions of parameter space, we find thatsystematic modeling uncertainties result in factor of ∼ to uncertainties in best-fit parameter values. Fortunately, twohighly sought-after parameters— α lo and L X /SFR—are rel-atively impervious to systematic uncertainties, leaving theexpectations from idealized forecasts largely unchanged.There are many as-yet-unexplored sources of systematicmodeling uncertainty in the reionization context. For exam-ple, we have completely neglected any numerical shortcom-ings of our model, e.g., the accuracy of the two-zone model of © 2020 RAS, MNRAS , ?? – ?? osmic dawn modeling systematics - . . z - . - . . . . c . . . log a - . - . . . b - . - . . . . c - . - . . . b . . . l o g a P D F PSSTTinker10
Figure 8. Posterior distributions for four nuisance param-eters introduced to relieve tension caused by modelingsystematics.
For the mock global 21-cm signal being fit, f ∗ doesnot evolve with redshift ( γ z = 0 ) and the parameters a , b , and c (see Eq. 9) are unused. the IGM, correlations between the fields composing the 21-cm background, or the idealized approach to high- z galaxyevolution. Of course, more detailed 3-dimensional reioniza-tion calculations sourced with galaxies modeled via semi-analytic prescriptions are likely to be more realistic in suchregards. However, more sophisticated models have system-atic uncertainties as well, caused by, e.g., approximate radia-tive transfer techniques and sub-grid models for sources andsinks. Similarly, extending models to redshifts z (cid:28) canin principle tighten constraints on parameters of interest,but may also introduce yet more systematic uncertainties(Behroozi et al., 2010).Given the ever-improving constraints on reionization,it is important to continue efforts to quantify the theoret-ical modeling error budget. Doing so will help to properlyreport uncertainties on parameters of interest, and identifythe biggest problem areas and thus guide dedicated effortsto mitigate them. For some applications, eliminating sys-tematics may be unnecessary, since some parameters of in-terest cease to be meaningful at arbitrarily high levels ofprecision. For example, most parameters in 21-cm modelsare averages over entire galaxy populations, with referencepoints anchored to uncertain empirical scaling relationships( L X /SFR) for which even factor of ∼ few level constraintsare interesting. However, tests of the cosmological model,e.g., non-cold dark matter, is an intrinsically more quantita-tive enterprise, in which all sources of systematic uncertaintymust be understood or eliminated. In these contexts, thereis much work still to be done before limits (or constraints)on various exotic models can be trusted to high precision.In this paper we have conducted a first study of howuncertainties in theoretical models can affect interpretationsof current and upcoming probes of cosmic dawn. Although the conclusions may seem pessimistic at first glance, theyin fact demonstrate the excitement that will accompany theexpected observational progress of the next few years: ratherthan a scenario where experiments will merely fine tune ahandful of model parameters, we will be in a situation wherethere exists a truly symbiotic relationship between theoryand experiment, with each side informing the other in fu-ture improvements. Together, they will lead the communitytoward a new understanding of our cosmic dawn. ACKNOWLEDGMENTS
The authors thank Saurabh Singh, Brad Greig, and Pe-ter Sims for helpful conversations that improved thismanuscript. J.M. acknowledges support from a CITA Na-tional Fellowship. H.L. was supported by a Science Un-dergraduate Research Award from the McGill UniversityFaculty of Science. A.L. acknowledges support from theNew Frontiers in Research Fund Exploration grant pro-gram, a Natural Sciences and Engineering Research Coun-cil of Canada (NSERC) Discovery Grant and a DiscoveryLaunch Supplement, a Fonds de recherche Nature et ech-nologies Quebec New Academics grant, the Sloan ResearchFellowship, the William Dawson Scholarship at McGill, aswell as the Canadian Institute for Advanced Research (CI-FAR) Azrieli Global Scholars program. Computations weremade on the supercomputer Cedar at Simon Fraser Univer-sity managed by Compute Canada. The operation of thissupercomputer is funded by the Canada Foundation for In-novation (CFI).
Software: numpy (Van Der Walt et al., 2011), scipy(Virtanen et al., 2020), matplotlib (Hunter, 2007), h5py ,and mpi4py (Dalcín et al., 2005). Data Availability:
The data underlying this article isavailable upon request, but can also be re-generated fromscratch using the publicly available ares code.
References
Barkana R., 2018,
Nat , 555, 71Barkana R., Loeb A., 2005,
ApJ , 626, 1Behroozi P. et al., 2020,
MNRAS , 499, 5702Behroozi P. et al., 2019,
MNRAS , 488, 3143Behroozi P. S., Conroy C., Wechsler R. H., 2010,
ApJ , 717,379Bhattacharya S. et al., 2011,
ApJ , 732, 122Boddy K. K. et al., 2018,
Phys Rev D , 98, 123506Bouwens R. J. et al., 2015a,
ApJ , 811, 140Bouwens R. J. et al., 2014,
ApJ , 793, 115Bouwens R. J. et al., 2015b,
ApJ , 803, 34Bowman J. D. et al., 2018,
Nat , 555, 67Chluba J., Thomas R. M., 2011,
MNRAS , 412, 748Cohen A. et al., 2020,
MNRAS , 495, 4845Dalcín L., Paz R., Storti M., 2005, Journal of Parallel andDistributed Computing, 65, 1108Dayal P. et al., 2013,
MNRAS , 434, 1486Eldridge J. J., Stanway E. R., 2009,
MNRAS , 400, 1019Fialkov A., Barkana R., Cohen A., 2018,
PRL , 121, 011101 © 2020 RAS, MNRAS , ?? – ?? Mirocha, Lamarre, & Liu
Fialkov A. et al., 2013,
MNRAS , 432, 2909Field G. B., 1958, Proceedings of the IRE, 46, 240Finkelstein S. L. et al., 2019,
ApJ , 879, 36Finkelstein S. L. et al., 2015,
ApJ , 810, 71Foreman-Mackey D. et al., 2013,
PASP , 125, 306Furlanetto S. R., 2006,
MNRAS , 371, 867Furlanetto S. R. et al., 2017,
MNRAS , 472, 1576Greig B., Mesinger A., 2015,
MNRAS , 449, 4246Greig B., Mesinger A., 2017,
MNRAS , 472, 2651Haardt F., Madau P., 1996,
ApJ , 461, 20Hunter J. D., 2007, Computing in Science & Engineering,9, 90Hutter A., 2018,
MNRAS , 477, 1549Hutter A. et al., 2020, arXiv e-prints, arXiv:2004.08401Kern N. S. et al., 2017,
ApJ , 848, 23Knebe A. et al., 2013a,
MNRAS , 428, 2039Knebe A. et al., 2013b,
MNRAS , 435, 1618Leitherer C. et al., 1999,
ApJS , 123, 3Lewis A., Challinor A., Lasenby A., 2000,
ApJ , 538, 473Liu A. et al., 2016,
Phys Rev D , 93, 043013Liu A., Shaw J. R., 2020,
PASP , 132, 062001Majumdar S. et al., 2014,
MNRAS , 443, 2843Mason C. A., Trenti M., Treu T., 2015,
ApJ , 813, 21Mesinger A., Ferrara A., Spiegel D. S., 2013,
MNRAS , 431,621Mesinger A., Furlanetto S., 2007,
ApJ , 669, 663Mesinger A., Furlanetto S., Cen R., 2011,
MNRAS
Meurer G. R., Heckman T. M., Calzetti D., 1999, The As-trophysical Journal, 521, 64Mirocha J., 2014,
MNRAS , 443, 1211Mirocha J., 2020,
MNRAS , 499, 4534Mirocha J., Furlanetto S. R., Sun G., 2017,
MNRAS , 464,1365Mirocha J., Harker G. J. A., Burns J. O., 2015,
ApJ , 813,11Mirocha J., Mason C., Stark D. P., 2020,
MNRAS , 498,2645Mitsuda K. et al., 1984,
PASJ , 36, 741Monsalve R. A. et al., 2017,
ApJ , 847, 64Muñoz J. B., Cyr-Racine F.-Y., 2020, arXiv e-prints,arXiv:2005.03664Muñoz J. B., Dvorkin C., Cyr-Racine F.-Y., 2020,
PhysRev D , 101, 063526Muñoz J. B., Loeb A., 2018,
Nat , 557, 684Murray S. G., Power C., Robotham A. S. G., 2013a, As-tronomy and Computing, 3, 23Murray S. G., Power C., Robotham A. S. G., 2013b,
MN-RAS , 434, L61Oesch P. A. et al., 2018,
ApJ , 855, 105Park J. et al., 2019,
MNRAS , 484, 933Planck Collaboration et al., 2018, arXiv e-prints,arXiv:1807.06209Press W. H., Schechter P., 1974,
ApJ , 187, 425Pritchard J. R., Furlanetto S. R., 2007,
MNRAS , 376, 1680Pritchard J. R., Loeb A., 2010, Physical Review D, 82,23006Pritchard J. R., Loeb A., 2012, Reports on Progress inPhysics, 75, 086901Qin Y. et al., 2020,
MNRAS , 495, 123Reed D. S. et al., 2007,
MNRAS , 374, 2Reichardt C. L. et al., 2020, arXiv e-prints,arXiv:2002.06197 Robertson B. E. et al., 2015,
ApJL , 802, L19Santos M. G. et al., 2010,
MNRAS , 406, 2421Schneider A., 2018,
Phys Rev D , 98, 063021Schneider A., Giri S., Mirocha J., 2020, arXiv e-prints,arXiv:2011.12308Sheth R. K., Mo H. J., Tormen G., 2001,
MNRAS , 323, 1Singh S. et al., 2018,
ApJ , 858, 54Sitwell M. et al., 2014,
MNRAS , 438, 2664Sun G., Furlanetto S. R., 2016,
MNRAS , 460, 417Tacchella S. et al., 2018,
ApJ , 868, 92Tashiro H., Kadota K., Silk J., 2014,
Phys Rev D , 90,083522Tinker J. L. et al., 2010,
ApJ , 724, 878Trac H., Cen R., Mansfield P., 2015,
ApJ , 813, 54Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Com-puting in Science & Engineering, 13, 22Virtanen P. et al., 2020, Nature Methods, 17, 261Warren M. S. et al., 2006,
ApJ , 646, 881Williams C. C. et al., 2018, The Astrophysical Journal Sup-plement Series, 236, 33Wouthuysen S. A., 1952, AJ , 57, 31Yung L. Y. A. et al., 2019, MNRAS , 483, 2983Zahn O. et al., 2011,
MNRAS , 414, 727Zygelman B., 2005,
ApJ , 622, 1356 © 2020 RAS, MNRAS , ?? ––