A Model-Insensitive Baryon Acoustic Oscillation Feature in the 21 cm Signal from Reionization
Christopher Cain, Anson D'Aloisio, Vid Iršič, Matthew McQuinn, Hy Trac
.. Draft version April 23, 2020
Typeset using L A TEX default style in AASTeX62
A Model-Insensitive Baryon Acoustic Oscillation Feature in the 21 cm Signal from Reionization
Christopher Cain, Anson D’Aloisio, Vid Irˇsiˇc,
2, 3
Matthew McQuinn, and Hy Trac Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA Kavli Institute for Cosmology, University of Cambridge Department of Astronomy, University of Washington, Seattle, WA 98195 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
ABSTRACTWe examine the impact of the baryon-dark matter relative velocity on intergalactic small-scalestructure and the 21 cm signal during the epoch of reionization (EoR). Streaming velocities reducedthe clumping of the intergalactic medium (IGM) on mass scales of ∼ − M (cid:12) . This effectproduced a distinct baryon acoustic oscillation (BAO) feature in the EoR 21 cm power spectrum atwave numbers k ∼ . ≈ ≈
10 % ionization). During this phase,the signal may appear at the 1 % (5 %) level at k ∼ . .
06) h/Mpcwith an amplitude that varies byless than a factor of 10 across a range of reionization histories. We also provide a simple model for thesignal originating from streaming velocity’s impact on ionizing sources, which can vary by 4 orders ofmagnitude depending on highly uncertain source properties. We find that the clumping-driven signalis likely to dominate the source-driven one unless population III star formation in halos of masses10 − M (cid:12) was efficient enough to drive the first 10% of reionization.1. IntroductionThe EoR was the last major phase transition in the Universe, during which the first sources of ionizing photonsre-ionized the IGM. In recent years, observational progress has been made towards constraining this epoch. The timingof Reionization has been constrained by Cosmic Microwave background (CMB) optical depth measurements (PlanckCollaboration et al. (2018)). Additional constraints are provided by observations of high-redshift quasars (e.g. Fanet al. 2006; Becker et al. 2015; McGreer et al. 2016; Davies et al. 2018; Bosman et al. 2018; Eilers et al. 2018; Beckeret al. 2019) and the population of high- z Lyman- α (Ly α ) emitters (e.g. Kashikawa et al. 2006; Schenker et al. 2012;Pentericci et al. 2014; Mesinger et al. 2015; Inoue et al. 2018; Weinberger et al. 2019). These observations have beeneffective at constraining timing of Reionization, but currently little is known in detail about the Reionization process.Forthcoming observations of the 21 cm spin-flip transition of neutral hydrogen promise a definitive window into theEoR (see e.g. Furlanetto et al. 2006, and references therein).Following a first detection of the EoR 21 cm signal, early efforts will focus on characterizing its brightness temper-ature fluctuations with the power spectrum. In this paper, we will investigate whether baryon-dark matter relativevelocities (or “streaming velocities”, or just v bc ; Tseliakhovich & Hirata 2010; henceforth T10), which were sourced Corresponding author: Christopher [email protected] a r X i v : . [ a s t r o - ph . C O ] A p r Cain et al. at recombination, were able to significantly impact the EoR 21 cm power spectrum. Though formally a second-ordereffect in perturbation theory, v bc was several times the baryon sound speed at decoupling, so its effect on baryonicstructure formation is important. Previous work has shown that v bc impacts a number of astrophysical processesin the early universe, including star formation in low mass halos (Maio et al. 2011; Greif et al. 2011; Schauer et al.2019), gas content of halos (Naoz et al. 2012), formation of direct-collapse black holes (Tanaka & Li 2014), the BAOfeature in the Ly α forest (Blazek et al. 2016; Hirata 2018), and possibly the formation of globular clusters (Naoz &Narayan 2014; Chiou et al. 2019). Of particular relevance for the current paper, v bc has been shown to modify the pre-reionization
21 cm signal at z ∼
20 through its impact on the properties of the first stars and galaxies (e.g. Dalalet al. 2010; Fialkov et al. 2012; McQuinn & O’Leary 2012; Ali-Haimoud et al. 2014; Cohen et al. 2016; Mu˜noz 2019).These papers have demonstrated that v bc can imprint distinct BAO features in the 21 cm power spectrum thatcould be detectable by future experiments.Streaming velocities impacted the universe at two scales that are particularly important for our investigation: (1)near the baryon jeans scale k J ∼ − h / Mpc, and (2) at the peak of the power spectrum of fluctuations in v , P v ( k ), which occurs at k ∼ − h / Mpc (T10, O’Leary & McQuinn (2012)). The former is roughly the minimumclumping scale of the pre-EoR gas. Recently, D’Aloisio et al. (2020) (henceforth D20) showed that gas clumpinesson this scale contributes significantly to the ionizing photon budget required to reionize the IGM. The suppressionof small-scale clumpiness caused by v bc , together with any impact v bc has on ionizing photon sources, will thereforetranslate into fluctuations in the neutral fraction that trace P v ( k ). Moreover, near 10 − h/ Mpc, P v ( k ) is a factor of ∼ larger than the EoR linear matter power spectrum. These facts suggest that v bc may have a pronounced effecton the EoR 21 cm power spectrum P ( k ) at large scales if it can be written in the form P ( k ) = b ,v P v ( k ) + matter terms (1)where b ,v is a linear bias factor coupling fluctuations in v bc to fluctuations in the 21 cm signal. This v bc -sourcedterm may be detectable in measurements of P even if b ,v < − . Since P v exhibits strong BAO features, itsappearance in measurements of P could serve as a “smoking gun” signature of Reionization. Such a signature wouldbe particularly helpful given the relatively featureless nature of the expected EoR 21-cm power spectrum, and theextreme difficulty of the measurement. Whereas previous studies have explored the coupling of P v with P throughthe reionization sources (Bittner & Loeb 2011; Cohen et al. 2016), the effect from gas clumpiness (or the “sinks”) hasnot been quantified in detail before.Employing a modified version of the code used in D20, we investigate with fully coupled radiative transfer andhydrodynamics simulations the impact of v bc on the sinks. We will use our simulation results to quantify their contri-bution to b ,v . We will also assess the potential impact of v bc on source properties during the EoR to see whetherthis effect may contribute at a similar level to sinks. As we will show, the contribution from the sinks is relativelyinsensitive to the details of reionization, being mostly fixed by the spectrum of cosmological density fluctuations andgas dynamics. This will allow us to more tightly bracket the coupling between P v and P from the sinks, as opposedto the highly uncertain effect from the sources.This paper is organized as follows. In §
2, we present the details of our simulation code and how we set up our initialconditions. In §
3, we present the results of our simulations and describe the impact of v bc during reionization. In § P analytically and estimate the magnitude of this term using our simulationresults. We also discuss potential contributions to b ,v coming from v bc ’s impact on source properties and assessthe detectability of the v bc -sourced signal. We summarize in §
5. Throughout this work, we assume the followingcosmological parameters: Ω m = 0 . Λ = 1 − Ω m (flat universe), Ω b = 0 . σ = 0 . n s = 0 . h = 0 . v bc . These are set up to track the response of a highly resolved patch of the IGM to ionizing radiation produced byexternal sources of constant intensity. Hence, we do not explicitly model galaxy formation in our simulations. Thisapproach allows us to assess the evolution of the sinks in a controlled manner, whereas it can be difficult to disentanglethe physical effects at play in a full simulation of reionization that also models the sources. Our approach also allowsus to achieve the required resolution for robustly modeling the clumpiness of the un-relaxed gas (see discussion inD20).This section describes the code (2.1), initial conditions (2.2), and simulations included in this work (2.3). We used a modified version of the Eulerian hydrodynamics code of Trac & Pen (2004) that includes the plane-parallelversion of the ray-tracing radiative transfer from Trac & Cen (2007). Our simulations were run on one large node ofthe Bridges supercomputer (Towns et al. 2014) and ran for 200 −
400 wall clock hours per simulation. We assume agas of primordial composition with H mass fraction X = 0 . Y = 1 − X . The number of hydrocells, RT cells, and DM mass elements are all equal to N (with N = 1024 in our fiducial runs), and our fiducial boxlength is L = 1 . /h , for a cell length of 1 kpc /h . The radiation is handled via plane-parallel ray tracing with anadaptive reduced speed-of-light approximation. Following D20, the radiation sources lie on the boundaries of cubical“sub-domains” of side length L dom = 32 kpc /h . The sub-domain structure of our boxes allows us to ionize all thegas at approximately the same redshift ( z re ) and to maintain a nearly constant photo-ionization rate throughout thebox, which simplifies interpretation of the gas evolution. We use a power law spectrum with intensity ∝ ν − . and 5frequency bins between 1 and 4 Ryd, roughly typical of the expected energy spectrum of reionization-era galaxies. Allthe important heating/cooling processes relevant to primordial gas are tracked by the code after the radiation turnson (see D’Aloisio et al. (2018)). In addition, we keep track of Compton scattering off the CMB at z > z re (whichis important for z (cid:38) . For a moredetailed description of the setup of our simulations, we refer the reader to § We generated Gaussian random field initial conditions at Recombination ( z = 1080), the time v bc was sourced, usingCAMB transfer functions (TFs). We did this to capture the cumulative effect of v bc self-consistently rather thanstarting from linear theory solutions at lower redshift (as was done in O’Leary & McQuinn (2012) and Ahn & Smith(2018)). Following those works, we used separate TFs for Baryons and DM to compute density and velocity growthfactors. The initial density and velocity fields were generated using the Zel’dovich approximation (see Padmanabhan1993, for a description). We modeled v bc by adding a constant velocity to the gas along the x direction at the initialredshift. This approximation is appropriate on scales (cid:46) several Mpc /h because v bc is coherent on those scales (T10).We tested the accuracy of our initial conditions prescription by showing that the matter power spectrum producedby our simulations agrees with the linear theory prediction at redshifts when it should (see Appendix B). We alsocompared our results to simulations initialized at a lower redshift to see whether starting from such a high redshiftproduced spurious shot noise. We found (as did Hirata (2018)) that this was not a significant effect. Our simulations are run with only hydrodynamics until z re . At this time, the box is rapidly filled with radiationand all the gas that cannot self-shield is ionized within a few times steps. The hydrogen photo-ionization rate Γ − (in units of 10 − s − ) at the boundaries of the sub-domains is a free parameter. We note that Γ − is nearly constantthroughout the box due to our sub-domain method. We ran simulations with Γ − = 0 . . z re = 6, 8, and12, v bc = 20, 41, and 65km / s. We used the simulations from D20 with these values of Γ − and z re to allow forcomparison to the no- v bc case. Note that D20 also considered box-scale density fluctuations by adding a constantbackground overdensity to some of their simulations. We do not do this here because it would make the parameterspace unmanageable given the computational cost of our simulations. In addition, simulating over-dense regions withnon-zero v bc requires a more complex treatment of the initial conditions (see Ahn (2016); Ahn & Smith (2018)).Simulations with v bc were run down to z stop = 5, 5 .
5, and 8 for z re = 6, 8, and 12, respectively, while the simulationstaken from D20 are all run to z = 5. Note that throughout this work, the quoted v bc values are those at z = 1080,after which v bc ( z ) ∝ (1 + z ). The values of v bc used here were chosen to facilitate evaluation of integrals of the form (cid:104) X (cid:105) v bc = (cid:90) ∞ dv bc X ( v bc ) P v bc (2) For z > z re , we set x e ( z ) = 0 . z − / .
6) + 9 . × − × θ (600 − z ) × z . + 0 . θ is the heaviside function) wasa good fit to the RECFAST free electron fraction. http://camb.info/ Cain et al. where X ( v bc ) is any quantity of interest and P v bc is the probability distribution of v bc in the universe, given by (Tseli-akhovich et al. 2011; Fialkov 2014). P v bc = (cid:18) πσ bc (cid:19) × πv e − v σ bc (3)where σ bc = 30 km/s is the RMS value. Equation 2 is the average of quantity X over the distribution of v bc in theuniverse. Assuming X ( v bc ) can be well-approximated by an order ≤ v bc , Equation 2 can be evaluatedexactly via Gaussian Quadrature with only the three v bc values used here.Our goal is to capture the impact of v bc on the formation of gas structures at high redshift and quantify howimportant this effect is once the gas becomes ionized. Our simulation setup is well-suited to achieve this goal. Oursimulations have resolution high enough to capture the impact of v bc on the gas at k ≥ Mpc /h , while being largeenough to include structures on mass scales of 10 − M (cid:12) , which should be relatively unaffected (Dalal et al. 2010).Hence, it is unlikely that our simulations will significantly under or over-estimate v bc ’s effect on the gas distribution(see Appendix C for some convergence tests). Second, our numerical setup allows us to isolate the effects of v bc on thesinks independently of its effect on sources, allowing for a straightforward interpretation of our results. Finally, ouruse of fully coupled hydro/RT will provide a realistic picture of how v bc ties into the reionization process. By modelingthe response of the sinks to reionization as in D20, we can make a physically realistic assessment of how important v bc is to their evolution. 3. Results We begin by visualizing the gas structure in runs with different v bc . Figure 1 shows 2D slices through the gas densityfield at redshifts of 7 .
9, 7 .
5, and 6 . v bc = 0, 41, and 65 km/s (top to bottom) for ( z re , Γ − ) = (8 , . z re , the gas ionizes quickly, reaching high temperatures. This rapidly increases thepressure in the high density gas filaments, which respond by expanding (“relaxing”) out of their DM potential wells,smoothing the gas density field considerably (see D20 for a detailed discussion). At redshift 7 .
9, most of the gas is stilltightly bound in these filaments, but by z = 6 . v bc is still visible, reflecting the integrated history of the un-relaxed gas. However, thedifferences largely disappear after the gas has relaxed. Even by z = 7 .
5, it is difficult to detect by eye any differencebetween the three runs, and in the relaxed limit at z = 6 . v bc affects the gas distribution the most on small scales, and it is precisely these scales that are smoothed bythe gas relaxation. Thus, the majority of the v bc effect does not survive the relaxation process. All this suggeststhat shortly after the gas is ionized, the recombination rate should be appreciably modulated by v bc because of itsimpact on small-scale structure. However, after some time passes the differences should largely disappear owing tothe smoothing effect of the relaxation process.For further clarity, in Figure 2 we show a 80 ×
80 kpc zoom-in of the gas density field at the same redshifts shownin Figure 1. This figure directly compares the initial and relaxed state of the gas. At z = 7 .
9, the structures aremuch more diffuse in the high- v bc runs, and the missing gas fills in some of the voids between structures. The effectis less prominent at z = 7 .
5, and almost absent by z = 6 .
5. Note also that the structures themselves are differentin the relaxed plot. This highlights the fact that the relaxation process effectively erases the initial conditions of theun-relaxed gas on small scales, including the v bc effect. This relaxation process makes it unlikely that any integratedhigh-redshift effect that affects only small scales will survive Reionization. Based on these results, we expect the recombination rate in a patch of the IGM to be substantially affected by v bc only relatively soon after z re . Here, we quantify the recombination rate by the clumping factor C R , defined to be the Figure 1.
2D slices through the gas density field for redshifts 7 .
9, 7 .
5, and 6 .
5, (left to right) and v bc of 0, 41, and 65 km/s(top to bottom). The relaxation of the gas is seen going from left to right. At z = 7 . v bc is visible, but after the gas relaxes the differences are too small to easily detect by eye. We show the results for z re = 8 andΓ − = 0 . Cain et al.
Figure 2.
Zoom-in on an 80 ×
80 kpc region in Figure 1. From left to right, the v bc values are 0, 41, and 65 km/s. The toprow shows a marked reduction in structure moving from smaller to higher v bc (left to right). At the lower redshift, there ismuch less of a difference ratio of the true recombination rate to that in a uniform-density IGM with constant temperature T ref , C R ≡ (cid:104) α B n e n HII (cid:105) α B ( T ref ) (cid:104) n e (cid:105)(cid:104) n HII (cid:105) (4)where α B is the case B recombination rate for hydrogen, n e is the free electron density, n HII is the HII number density,and T ref = 10 K. Since all of our simulations have mean densities equal to the global mean and the hydrogen isalmost completely ionized after z re , we approximate (cid:104) n HII (cid:105) ≈ n H ( z ) and (cid:104) n e (cid:105) ≈ n H ( z )(1 + n He ( z ) /n H ( z )), i.e. assum-ing singly ionized helium, where n H ( z ) and n He are the cosmological average number densities of H and He, respectively.We plot C R vs. cosmic time (∆ t ) since z re in the left panels of Figure 3 for z re = 12 (top), 8 (middle) and 6 (bottom)and Γ − = 0 . . C R ( v bc ) /C R ( v bc = 0) i.e. C R as a fractionof the no- v bc case. We see that the percentage difference between the different v bc values is largest ≈ −
10 Myr after z re , the time at which C R is also at a maximum. At this time, the ionizing radiation has penetrated deep into themost overdense regions, but the gas has not yet had time to dynamically relax. So, the recombination rate is set bythe clumpiness of the initial density field, which is significantly modulated between patches with different v bc . After ∼
200 Myr , the gas has had time to relax and the fluctuations in C R sourced by v bc have largely disappeared . Weemphasize that if our code did not capture the relaxation process, we would significantly over-estimate how much v bc reduces the recombination rate. Still, the effect on recently ionized gas is not insignificant, reaching ∼ −
20% for z re = 12 and ∼
10% for z re = 6 for v bc = 41 km/s. Because of the patchy nature of Reionization, at any time therewill always be some regions in the IGM that were ionized recently and haven’t had time to relax. In these regions,the recombination rate will depend non-negligibly on v bc , potentially leading to detectable fluctuations in the IGMneutral fraction (see the next section).It has been shown that X-ray heating prior to Reionization also reduces the clumpiness of the gas. D20 ran asimulation in which they set the pre-Reionization temperature to a uniform 1000 K to gauge the maximum effect ofX-ray pre-heating. They found that C R was suppressed in a fashion similar to what we find here due to v bc . Thisoccurs because X-ray pre-heating raises the pre-EoR Jeans mass, which eliminates structure on the smallest scales. Inthe event that pre-heating is significant, we expect the importance of v bc to be reduced somewhat as the two processesaffect structure at the same mass scales.Another important quantity during the EoR is the MFP of ionizing photons, which quantifies the typical distancean ionizing photon can travel before being absorbed. We calculated the MFP from our simulations using the approachof Emberson et al. (2013) (see D20 for details). Figure 4 shows the MFP from the same simulations as in Figure 3.We find that v bc modulates the MFP by roughly the same percentage that it does the clumping factor, but in theopposite direction. This result is consistent with Figure 3 because a less clumpy IGM should allow ionizing photonsto travel further on average before being absorbed. The behavior with time is also qualitatively the same as for C R ;early on, the MFP is modulated by 10 − v bc runs starts out slightly below the v bc = 0 case. This is likely because the first regions in the box to ionize arethe under-dense ones, which are slightly less dense in the v bc = 0 case because more of the gas is locked up in small,dense structures. Unlike for C R , the percentage difference in the MFP from v bc is small compared to the differencebetween the runs with high and low Γ − . Thus, spatial variations in MFP sourced by v bc should be subdominant tothose coming from fluctuations in the photo-ionization rate.We conclude from these results the IGM recombination rate and MFP are impacted by v bc at the 10 −
20% levelin patches of the universe that have reionized recently. After the gas has had time to relax, the effect we observe islargely erased. During Reionization, there will always be a small percentage of the IGM that was ionized recently,and will therefore display a moderate variation in C R and MFP due to v bc . These patches will consume fewer ionizingphotons due to recombinations than they would in the absence of v bc because of the reduced clumpiness, which willspeed up the local Reionization process. Patches with different values of v bc will be affected differently, leading to Note that for ∆ t >
100 Myr , the clumping factor is actually larger for the v bc = 41 case than for v bc = 0, particularly in the z re = 6case. We believe this offset is due to the difference in starting redshift between the v bc = 0 simulations and the others. We tested this byrunning a set of small box simulations starting at different redshifts, and found that starting at z = 300 produces a ∼
2% suppression in C R relative to starting at z = 1080. This difference is not large enough to impact our results. Cain et al.
Figure 3.
Left: Clumping factor vs. cosmic time since z re = 12 (top), 8 (middle), and 6 (bottom) for Γ − = 3 . . v bc . The difference in C R is largest shortly after the radiation turns on, but attenuates as the gasrelaxes. Right: Ratio of C R to C R ( v bc = 0) for both values of Γ − . The difference peaks 5 −
10 Myr after z re and steadilydeclines thereafter, reaching ∼ a few percent 200 Myr after z re . We note that the kinks in the right panels are due to thesparse time-stepping at small ∆ t . fluctuations in the local ionized fraction that trace fluctuations in v bc . Granted, these fluctuations should be small,but we also argued in § P v and P could result in a detectable BAO feature inthe latter. Exploring this possibility is the subject of the next section. Figure 4.
Same as Figure 3, but for the MFP of ionizing photons. The effect of v bc on the MFP is roughly equal and oppositeto its effect on the clumping factor. After ∼ −
10 Myr, the difference reaches 25% in the most extreme case, but has largelydisappeared after ∼
200 Myr. This is much less than the factor of ∼
10 difference between the runs with high and low Γ − . Cain et al.
4. Impact of v bc on the 21 cm SignalIn this section, we model the impact of v bc on the EoR 21 cm signal and estimate its detectability. We begin withsome preliminaries in § § v bc -dependent term. Using this model, we quantify how large the bias factor b ,v must be to produce a detectable imprint on P . In § § b ,v fromionizing photon sinks ( § § § The EoR 21 cm signal is produced by neutral patches of the IGM that have yet to ionize, so its spatial fluctuationsset by those of the neutral fraction and the density field (Furlanetto et al. 2006; Pritchard & Loeb 2012). For simplicity,we neglect redshift-space distortions (Jensen et al. 2013, 2016) and adopt the typical assumption that Ly α couplingand the first X-ray sources drove the 21 cm spin temperature to be T s >> T CMB by the time reionization largelybegins (Wouthuysen 1952; Madau et al. 1997; Pritchard & Furlanetto 2007). Under these assumptions, the 21 cmbrightness temperature T can be written T ( r , z ) = ˆ T ( z ) x HI ( r )(1 + δ ρ ( r )) (5)where δ ρ is the (nonlinear) matter overdensity, x HI is the neutral hydrogen fraction, and ˆ T depends only on cosmo-logical parameters and redshift. To first order in over-densities, fluctuations in T are proportional to δ HI + δ ρ , where δ HI is the overdensity in the neutral hydrogen fraction. Because the highest-density regions ionized first, δ HI and δ ρ will generally have opposite signs early in the EoR (e.g. Giri et al. 2019). For a given wavenumber k , the signal willreach a local minimum when ˜ δ ρ ( k ) = − ˜ δ HI ( k ) where the tildes denote the Fourier Transform (FT). At this time, thedominant density and ionization terms will cancel out and the signal will be sourced entirely by higher-order terms,one of which should be the v bc term in Equation 1. The signal will later reach a local maximum before disappearingentirely when there is no more neutral hydrogen. Forthcoming surveys will characterize the EoR brightness temperature fluctuations with the power spectrum, definedas P ( k ) ≡ (cid:104) ˜ δ ( k )˜ δ ( k ) (cid:105) where δ ≡ x HI (1 + δ ρ ). While modelling the signal is quite complicated and requiresnumerical simulations (see e.g. Furlanetto et al. 2019; Koopmans et al. 2019; Parsons et al. 2019, for a generaldiscussion), MQ18 showed that on large scales the power spectrum can be described surprisingly well with perturbationtheory. They modeled the 21 cm signal using a multi-parameter bias expansion, keeping the minimum number of termsthat produced a reasonable fit to the signal in numerical simulations of Reionization. At large scales and early times,they obtained a good fit using a model with only three parameters; their “minimal model” is given by˜ δ = b (cid:18) − R k (cid:19) ˜ δ ρ + b ˜ δ ρ (6)where ˜ δ ρ is the FT of the total matter over-density, and b , b and R eff are time-dependent but scale-independentbias factors. In what follows, we will approximate δ ρ ≈ δ in Equation 6, where δ is the linear matter over-density;this approximation is valid at the redshifts and scales considered here. R eff roughly characterizes the size of ionizedbubbles, which should be small compared to 1 /k at times and scales considered here, so we will drop it.We will build upon this model by adding a term proportional to the v “overdensity”, δ v ≡ [ v − σ bc ] /σ bc . First,we write x HI as x HI ≡ (cid:104) x HI (cid:105) (1 + δ HI ) (7)where the angle brackets denote an average over the whole IGM. Next, we assume that δ HI is a biased tracer of δ , δ ,and δ v and that δ ρ traces δ v . Then we have δ HI = b HI , δ + b HI , δ + b HI , v δ v δ ρ = δ + b ρ, v δ v (8) Note that redshift-space distortions have been shown to distort the 21 cm power spectrum at the scales we consider in what follows.Our conclusions are somewhat dependent on these effects being small. Figure 5.
Comparison of the power spectra entering Equation 11 at redshifts 5.8 (solid curves), 8.0 (dashed curves), and 10(dotted curves). P v (black dashed curve) is the same at all three redshifts. At scales k (cid:47) × − h/ Mpc, P v dominates overthe matter terms. where the coefficients are bias parameters. Combining Equations 7 and 8 with the definition of δ and dropping allterms 3rd order or higher yields δ = (cid:104) x HI (cid:105) (1 + [1 + b HI , ] δ + b HI , δ + [ b HI , v + b ρ,v ] δ v ) (9)Comparing this to Equation 6, we identify b = (cid:104) x HI (cid:105) (1 + b HI , ) and b = (cid:104) x HI (cid:105) b HI , in the case with no δ v term, sowe can substitute accordingly to get δ = (cid:104) x HI (cid:105) + b δ + b δ + b ,v δ v (10)where b ,v ≡ (cid:104) x HI (cid:105) [ b HI , v + b ρ,v ]. Taking the Fourier transform of both sides of Equation 10 and squaring gives ,assuming cross-terms are negligible, P ( k ) = b ,v P v ( k ) + b P ( k ) + b P ( k ) (11)where P and P ∝ P (cid:63) P are the first and second order total matter power spectra, respectively, and P v ( k ) is theFourier transform of (cid:104) δ v ( x ) δ v ( x + r ) (cid:105) . Ali-Haimoud et al. (2014) found that the linear and quadratic density fieldsare uncorrelated on all scales because δ ( δ ) have odd (even) dependence on v bc . By the same reasoning, δ and δ v should be uncorrelated as well. Thus our assumption of negligible cross-terms is exact for the terms involving δ . As long as the cross-term between the quadratic terms is sub-dominant to b P ( k ), we may safely ignore it whencomparing the v bc term to the contribution from the density terms, as we will do shortly . Figure 5 plots the dimen-sionless power spectra (∆ ≡ k P ( k ) / π ) for δ , δ , and v at redshifts 5 .
8, 8, and 10 (note that P v is independentof redshift). For all times shown here, P v > P > P for k (cid:47) × − h/ Mpc, with the differences growing largerfor with increasing redshift and decreasing k . Moreover, P v shows strong baryon acoustic oscillation (BAO) features, Note that the zeroth-order term becomes a delta function at k = 0 in Fourier space and thus does not contribute. The neglected cross term will either be featureless and can therefore be absorbed into the δ term or will have BAO features, in whichcase it may contribute to the signal we are studying. In either case, we can ignore it as long as it is small. Cain et al. suggesting that it’s appearance in P would be distinct even if it only contributes to the total signal at the ∼
10% level.Figure 7 of MQ18 shows how b and b evolve with time for three idealized models of Reionization. t. The signalreaches a maximum amplitude when b is largest (since the P term dominates). Depending on the model, themaximum value of b is between 0 . and 2 , at which time b ∼ − . As discussed in 4.1, there is a time earlyin Reionization when b = 0, and the amplitude of P is dominated by the second order term, with b ∼ − .In Figure 6, we show how the v bc term impacts P in two representative cases for several values of b ,v . Inthe top left panel, we plot the signal at z = 8 assuming b = b = 1, representative of the P maximum, with b ,v ∈ { , − , − , − , − } . The top right panel shows the same plot at z = 10, b = 0, and b = 1, represen-tative of the P minimum, with b ,v ∈ { , − , − , − , − } . The bottom panels show the range of b ,v P v as a fraction of the total signal without v bc for the bias parameters considered in the top panel. At the P maximum, b ,v (cid:38) − is required to produce a ∼
10% effect on the signal at 10 − Mpc/h, whereas at the P minimum thesame effect is achieved with b ,v (cid:38) − . The takeaway here is that v bc has its largest fractional effect at the P minimum, where the linear contributions from density and ionization cancel, and P is set by higher order terms inthe bias expansion.We emphasize that the results shown in Figure 6 are entirely agnostic about the cause of b ,v . In general, we canwrite b ,v = (cid:104) x HI (cid:105) (cid:2) b ρ,v + b sink HI,v + b source HI,v (cid:3) (12)where the second and third terms come from v bc ’s impact on the sink and source properties, respectively. In whatfollows, we will make the assumption that b ρ,v is essentially 0, since v bc does not affect the shape of the linear matterpower spectrum at scales near 10 − h/ Mpc (T10). In the next two subsections, we will derive a rigorous model for thecontribution to b ,v from ionizing photon sinks informed by our simulation results in §
3, and assess analytically thepotential contribution from sources. b ,v from Sinks In this section, we present our model for the contribution to b ,v from the sinks, which we can evaluate using thesimulation results presented in §
3. Our approach is to relate fluctuations in v to fluctuations in the ionized fraction x i via the effect on the clumping factors that we measure from our simulations. From this relationship, we will derivean expression for the sinks bias b sinkHI , v and relate this to the corresponding 21 cm sinks bias b sink21 ,v .We begin with the reionization “accounting equation”, given by (Madau et al. 1999),˙ x i = (cid:15)n H − (cid:104) C R (cid:105) α B n e x i (13)where (cid:104) C R (cid:105) is given by (cid:104) C R (cid:105) ( t ) = (cid:90) z ( t ) z dz re P z re ( x i ( t )) (cid:90) ∞ dv bc P v bc C R ( z re , v bc , t ) (14)Here, (cid:15) is the emissivity of ionizing photons, x i is the ionized fraction, and n e , n HII , and n H are the mean free electrondensity, HII number density, and H number density, respectively. The integral runs from the start of Reionization at z to redshift z ( t ). The clumping factor, C R ( z re , v bc , t ), which depends on the local values of v bc and z re , is extracteddirectly from our simulations. Equation 3 gives P v bc , and P z re ( x i ( t )) is the probability distribution of z re , P z re ( x i ( t )) = dx i /dz re x i ( t ) − x i ( z ) (15)Note that (cid:104) C R (cid:105) depends on the ionization history x i ( t ) through P z re ( x i ( t )); this is how our model accounts for the“patchiness” of Reionization. However, there are two important dependencies missing from Equation 14: Γ − and thelocal over-density. Because Reionization starts is “inside-out” i.e. moves from higher to lower density regions (Ciardiet al. 2003; Furlanetto et al. 2004; Mesinger et al. 2011) over-dense regions are more likely to ionized at higher z .Because these regions have a higher density of sources, they will also have higher-than average Γ − , and the impactof v bc will be different from what we measure here at mean density. We leave a thorough investigation of these3 Figure 6.
Top Left: Dimensionless 21 cm power spectrum at z = 8 for b = b = ( − . for b ,v = 0 (solid black), 10 − (red), 10 − (cyan), 10 − (yellow), and 10 − (green). This plot shows a conservative estimate of P near its maximum. TopRight: The same at z = 10 for b = 0, b = 1, and b ,v = 0 (solid black), 10 − (magenta), 10 − (blue), 10 − (red), and 10 − (cyan). This plot shows an estimate of P at its minimum. The bottom panels show the v bc -sourced component as a fractionof the signal without v bc for each of the bias factors. correlations to future work.We perturb Equation 13 by assuming it holds for a spherical patch of the IGM of radius r with mean ionized fraction x ri ≡ (cid:104) x i (cid:105) (1 + δ rx i ) and clumping factor C r R = (cid:104) C R (cid:105) (1 + δ rC R ). In doing so, we take Equation 13 to be locally true withinthe patch. This is only strictly true if the MFP of ionizing photons is << r . Since we are primarily interested in therange 10 − h/ Mpc < k < − h/ Mpc, we will take 60Mpc /h (cid:46) r (cid:46) /h for the perturbation scale. Figure 11of D20 shows how the MFP in ionized regions evolves with time for our assumed ionization history (their solid redcurve). For Γ − = 0 .
3, MFP << /h at z = 10 and at z = 8 it is still an order of magnitude smaller. In theΓ − = 3 . ∼ z = 10 and is comparable to 60Mpc /h at z = 8. However,since this is the MFP for ionized regions only, the MFP with neutral regions included will be considerably smaller,especially early in Reionization. Moreover, typical values for Γ − extracted from the Ly α forest are in the range0 . − . Cain et al. during the majority of the EoR.Since Equation 13 is also satisfied by the IGM mean values (cid:104) x i (cid:105) and (cid:104) C R (cid:105) , we can solve for the perturbation δ rx i ( t )(see Appendix A for details). Assuming P z re ( x i ( t )) is roughly scale-independent, we may write δ rx i ( t ) = b x i ,v ( t ) δ rv (16)where δ rv is time-independent and b x i ,v ( t ) is scale-independent. An expression for b x i ,v ( t ) can be obtained byTaylor-expanding δ rC R to first order in δ rv . Since b HI ,v = −(cid:104) x i (cid:105) / (cid:104) x HI (cid:105) b x i ,v and b ,v = (cid:104) x HI (cid:105) b HI ,v , we have b sink21 ,v = −(cid:104) x i (cid:105) b x i ,v (17)Hence, we obtain a model for the v bc term in Equation 11. The assumption of a scale-independent P z re ( x i ( t )) is validprovided that spatial fluctuations in P z re ( x i ( t )) on the perturbation scale are small compared to the global mean (givenby plugging (cid:104) x i (cid:105) into Equation 15) at each redshift. This is not immediately obvious because δ rx i implicitly containsnot only the v bc perturbation term, but also matter terms analogous to those in Equation 10. So, it is important tocheck that P z re ( x i ( t )) is roughly homogeneous on the scales considered here. Figure 3 of Nasir & D’Aloisio (2019)plots the distribution of z re for three different models of Reionization. In these plots, the distribution of z re appearsto be roughly homogeneous at scales r (cid:38) /h , justifying our approximation of a scale-independent bias factor.To compute b sink21 ,v from our simulations, we must first solve Equation 13 for (cid:104) x i (cid:105) and (cid:104) C R (cid:105) by plugging in oursimulation results for C R ( z re , v bc , t ) in Equation 14. The integral over v bc in this equation can be done via GaussianQuadrature as discussed in 2.3. Integrating over z re requires interpolating in two dimensions between the z re = 12,8, and 6 C R data as was done in D20 (see their Figure 11). To solve Equation 13, we assume the uniform emissivityfunction from Robertson et al. (2015) for (cid:15) and that n e and n H assume their cosmological mean values at each redshift.Once we have the global history, we can compute the bias (see Appendix A). Figure 7 shows the results of this exercisefor several Reionization histories. The left panel plots | b sink21 ,v | vs. x i for each history and the right panel plots x i vs. redshift. Our fiducial history (red solid curve) starts Reionization at redshift z = 12 and uses the C R ( z re , v bc , t )from our simulations with Γ − = 0 .
3. We also include histories using C R ( z re , v bc , t ) from our Γ − = 3 . − ”, solid cyan), z = 10 (Late Start, magenta dashed), an emissivity that is 50% higher (“High Emissivity”,blue dashed). Lastly, we include a “Best Case” model (green dashed) in which the first ∼
10% of the IGM is ionizedalmost instantly, and we take C R ( z re , v bc , t ) from our Γ − = 3 . b sink21 ,v at an ionized fraction of 10 − x i for which P reaches a minimum, where v bc has its largest fractional effect. This is because the “flash-ionized” patches reachthe time at which the v bc effect is largest coherently, so they contribute maximally to b sink21 ,v all at once.The bias squared varies in the range 10 − − − depending on ionization history near x i = 0 .
1, but approaches afew times 10 − in all histories by x i = 0 .
5. To get a 1% level effect in P ( z = 8) at k = 10 − h/Mpc would require | b sink21 ,v | ∼ − (see Figure 6), so for any of these histories the effect would be sub-percent level at the epoch ofmaximum P . Even at k ∼ − Mpc/h, where the difference between P v and the linear terms is much larger, theeffect would still only be a few percent for the bias factors measured here. However, at the epoch of minimum P the results are more promising. The Best Case model gives | b sink21 ,v | ≈ − at this time, which is enough to changethe signal by a few percent at 10 − h/Mpc and by ∼ − h/Mpc. The other histories (which are physicallyrealistic) give a ∼
1% effect at 10 − h/Mpc and tens of percent at 10 − h/Mpc. Note that the curves in the left panelof Figure 7 are very similar (even at low ionized fraction) despite the significant differences in the ionization historiesin the right panel. This suggests that the sinks bias is constrained to be | b ,v | ∼ − − − regardless of thedetails of Reionization, e.g. the nature of the source population. As we will see in the next section, this is not true ofthe v bc term coming from the sources themselves. b ,v from Sources Here, we discuss possible contributions to the source bias, b source21 ,v . Since our simulations do not model the sources, wewill proceed analytically and use results from the literature where appropriate. Previous work has demonstrated that v bc has an important effect on the star formation rate in “minihalos,” with masses 10 − M (cid:12) , within which it is5 Figure 7.
Left: | b sink21 ,v | vs. ionized fraction for each of the EoR histories discussed in the text. At low ionized fractions,the bias varies in the range 10 − − − between the different models, but settles down to a few × − for all the models latein Reionization. None of these models produce bias factors large enough to produce a detectable signal at the epoch of 21 cmmaximum (see Figure 6). At the 21 cm minimum near x i = 0 .
13, the Best-Case model gives | b sink21 ,v | ≈ − , which wouldalter the minimum signal by a few percent at k = 10 − h/ Mpc and by a factor of 2 at k = 10 − h/ Mpc. Among the physicallyrealistic histories, the one with high Γ − gives the largest bias. Right: The ionized fraction as a function of redshift for eachhistory. Note that Γ − changes the ionization history very little early on, but raises | b sink21 ,v | by a factor of ∼ − believed the first stars (Pop III) formed (Dalal et al. (2010), McQuinn & O’Leary (2012), Fialkov et al. (2012), Mu˜noz(2019), to name a few). Primarily, v bc raises the minimum halo mass able to form stars, thereby modulating theradiation backgrounds produced. As a result, v bc could modulate the 21 cm spin temperature T s near the start ofReionization. This could change the signal significantly if the assumption T s >> T CMB is not valid near the startof Reionization, which would introduce an extra factor of 1 − T CMB /T s in Equation 5. Additionally, the suppressionof ionizing photon production by v bc would work opposite the direction of the sinks by slowing down Reionizationlocally, thereby increasing the 21 cm signal in patches with higher v bc . Although it is widely believed that halos moremassive than 10 M (cid:12) (which are less affected by v bc ) drove Reionization, the degree to which Pop III star formationin minihalos contributed to its early phases is highly uncertain. In this section we continue to assume T s (cid:29) T CMB andwe attempt to quantify the coupling of P v with P through the source bias, b source21 ,v .To estimate the impact of v bc on halos, we begin with a general expression for the star formation rate density(SFRD) during Reionization (Sun & Furlanetto (2016), Irˇsiˇc et al. (2019)).˙ ρ SFRD = (cid:90) ∞ M c ( v bc ) dM n ( M ) f (cid:63) ( M ) ˙ M Ω b Ω m (18)where n ( M ) is the halo mass function, ˙ M is the halo accretion rate, and f (cid:63) ( M ) is the mass-dependent star formationefficiency. This integral contains contributions from both Population II and III stars, with the mass cutoff at M atom ≡ M vir ( T vir = 10 K), is the atomic cooling threshold, given by inverting Eq. 26 of Barkana & Loeb (2001). For the lowerlimit M c ( v bc ), we assume M c ( v bc , z ) = M vir (500K , z ) × (cid:18) V cool ( z, v bc , J LW = 0) V cool ( z, v bc = 0 , J LW = 0) (cid:19) × [1 + B (4 πJ LW ) β ] (19)where J LW is the specific intensity of the Lyman-Werner (LW) background in units of 10 − erg / s / cm / Hz / sr, ( B, β ) =(7 , .
47) for the regular feedback model in Mu˜noz (2019), and V cool ( z, v bc , J LW = 0) is the minimum circular velocity6 Cain et al. for star formation in the absence of LW feedback, derived from simulations by Fialkov et al. (2012) (their Eq. 2). Toobtain J LW , we combine Eq. 8 of Mebane et al. (2018) with Eq. 6 of McQuinn & O’Leary (2012) and include the LWopacity correction from Irˇsiˇc et al. (2019) to obtain J LW = 7 . π × (1 + z ) H ( z ) e − τ LW ( N IILW ˙ ρ IISFR + N IIILW ˙ ρ IIISFR ) (20)where we take N IILW = 9690, N IIILW = 10 , and e − τ LW = 0 . H ( z ) andSFRD are km / s / Mpc and M (cid:12) / yr / Mpc , respectively. For f (cid:63) ( M ), we used the form in Fialkov et al. (2014) for PopIII stars (their Eq. 2) and the form in Furlanetto et al. (2017) for Pop II stars (their Eq. 10), where we have tunedthe parameters of the latter to give Pop II SFRDs that agree well with the results of Visbal et al. (2020). The PopIII star formation efficiency at M atom , f (cid:63) , is a free parameter in our model. To evaluate Equation 18, we use theSheth-Torman mass function and the halo accretion rate given by Trac et al. (2015) (their Eq. 11) which is calibratedfrom high-redshift simulations. Since J LW and ˙ ρ SFR are interdependent, we use an iterative scheme to simultaneouslysolve for them given a value of f (cid:63) .Following the same formalism as Dalal et al. (2010) (see their section 2 for details), it can be shown that˙ ρ SFR = (cid:104) ˙ ρ SFR (cid:105) (1 + b SFR , v δ v ) (21)where b SFR , v = − (cid:104) v ˙ ρ SFR (cid:105) σ bc (cid:104) ˙ ρ SFR (cid:105) (22)and the averages are over P v bc . Depending on the magnitude of J LW , M c may be larger or smaller than this cutoff, sowe will include contributions from both populations of stars to the bias. From here, we can work out an expressionfor b source21 ,v using the same strategy as in § (cid:104) (cid:15) (cid:105) n H δ (cid:15) = A He N IIIion f IIIesc (cid:104) ρ − m ˙ ρ IIISFR (cid:105) b IIISFR , v δ v ≡ b III (cid:15),v δ v (23)where A He = 1 . N ion is the number of ionizing photons produced per stellar baryon, f esc is the escape fraction, and ρ m = Ω m ρ crit ( z = 0) is the present-day matter density of the universe. For Pop III stars we assume f esc = 0 . N ion = 40000. We do not include a contribution from Pop II stars to the bias because these stars are expected to formin atomic cooling halos that are unaffected by LW feedback (Fialkov et al. 2012). We therefore do not expect v bc byitself to raise the minimum circular velocity for star formation in these halos above the threshold set by the atomiccooling limit. If we ignore the recombination term in Equation 13, the 21 cm bias b source21 ,v is simply the time integralover b III (cid:15),v . We start this integral at z = 30, which is early enough to account for the full cumulative impact of Pop IIIstars on the signal (see Mu˜noz (2019)). Note that this bias factor has the opposite sign of the sinks bias because v bc reduces the number of ionizing photons being produced.The left panel of Figure 8 shows the result of this exercise for f (cid:63) = 10 − (cyan-dashed) and a range of f (cid:63) going from10 − to 10 − (blue shaded band). This range brackets the values generally considered in the literature (e.g. Trenti& Stiavelli 2009; Visbal et al. 2018; Irˇsiˇc et al. 2019; Visbal et al. 2020) as well as other sources of uncertainty (seebelow). At z = 10 for f (cid:63) = 10 − , the source bias is comparable to the sinks bias ( ∼ a few times 10 − ) and spans arange of about 2 orders of magnitude above and below this (since | b source21 ,v | ∝ | f (cid:63) | ).The right panel shows the average SFRD for both populations of stars, with the range given for Pop III corre-sponding to the range of bias factors in the left panel. The Pop III SFRDs in our model agree reasonably well withthose in Visbal et al. (2020) for the same Pop III star formation efficiencies . Note that the bias shown in Figure 8is formally an upper limit because we neglected the recombination term in Equation 13. It is therefore likely thatthe source bias is less important than the sinks bias early in the EoR . However, the result above depends strongly That work assumed a constant Pop III f (cid:63) , but because ours depends rather weakly on mass, their f III corresponds closely with our f (cid:63) Note that Mu˜noz (2019) finds much larger effective bias factors at z ∼
20 than we show in Figure 8. However, in that work, the mainsource of coupling between v bc and P was the coupling between spin temperature and gas temperature rather than between emissivityand ionized fraction. We briefly address the possibility of spin temperature fluctuations below. Figure 8.
Left: source bias squared for f (cid:63) = 10 − (cyan), and the range 10 − < f (cid:63) < − (red band). The bias grows withincreasing f (cid:63) , and spans 4 orders of magnitude over the range we consider here for that parameter. Right: SFRD for pop IIIand II stars in our model for the same range of f (cid:63) shown in the left panel. The SFRD for both populations agrees well withFig. 2 Visbal et al. (2020) for f (cid:63) = 10 − . For that choice of f (cid:63) , our star formation efficiency parameters for both populationsare similar to the fiducial values used in that work. on the assumed values of f (cid:63) , as well as on how the critical mass M c is modeled, the assumed accretion rate ˙ M , andthe exact relationship between star formation rate and J LW . For example, weaker LW feedback will result in a largerbias, and a smaller minimum Pop III star formation mass (in the absence of v bc or J LW ) will increase the bias asPop III stars will play a larger role in Reionization. In light of these large uncertainties, we cannot draw definitiveconclusions about the magnitude of the EoR source bias. Note that these uncertainties highlight the relative precisionof our model for the sinks bias.The case in which T s (cid:29) T CMB does not hold is even more difficult to assess. Conventional wisdom would suggest(O’Leary & McQuinn (2012)) that by the time Reionization starts, T s has been coupled to the gas temperature T K such that T s >> T CMB , as assumed in Equation 5. However, in scenarios where reionization starts very early, or Ly α scattering is inefficient at coupling T s to T K , this assumption may not be true at the times of interest. In this scenario,to first order Equation 5 would be multiplied by a factor 1 − T CMB (cid:104) T s (cid:105) + T CMB (cid:104) T s (cid:105) δ T s . The fluctuation δ T s would occur becausehigher v bc reduces the production of Ly α photons in 10 − M (cid:12) halos, making it harder for those photons to couple T s to T K . This bias would work in the same direction as the sinks bias because it would reduce the amplitude of the21 cm signal in regions with higher v bc . For reionization models in which T CMB = (cid:104) T s (cid:105) occurs early in reionization,the signal at that time would be dominated to linear order by the δ T s term, which may trace δ v . This may offeranother window of time during Reionization at which v bc could be detectable. Estimating the magnitude of this effectis beyond the scope of this paper, so we leave it to future research. v bc Here we briefly discuss the detectability of the expected signal near the EoR 21 cm minimum (where the fractionaleffect of v bc is likely to be largest). Cosmological 21 cm experiments like the Square Kilometer Array (SKA, Koopmanset al. (2015)) and the Hydrogen Epoch Reionization Array (HERA, DeBoer et al. (2017)) are expected to be able toprobe the scales discussed here with much higher sensitivity than current experiments, which are struggling to detectthe EoR signal at its maximum (EDGES, LOFAR, etc.). Trott & Pober (2019) show, in their Figure 1, the levels ofnoise in the dimensionless 21 cm brightness power spectrum ∆ b expected for several 21 cm experiments at z = 8 . We note that previous work (e.g. McQuinn & O’Leary 2012; Mu˜noz 2019) has studied this effect at 15 < z <
30, but not duringReionization. Cain et al. assuming 1000 hours of integration, including SKA and HERA. In the best-case scenario of thermal noise only, theyfind uncertainties of ≈ × − (1 × − ) mK for SKA (HERA) at the smallest wavenumber that both experimentscan detect, k = 0 . h/ Mpc. Assuming the thermal noise power spectrum scales as ν − α where α = 2 .
55 (Greig et al.2020), the thermal noise will be larger by a factor of ∼ . z = 10 than at z = 8 .
5. The resulting uncertaintiesare a factor of ∼ k = 0 . h/ Mpc. This suggests that the signal near the 21 cm minimum is marginally below the currentdetection limit of SKA and HERA, assuming these modes are not inaccessible due to foregrounds (see Lanmanet al. (2019) for a discussion of foreground contamination). We therefore suggest that future versions of these ex-periments may be able to detect the EoR 21 cm minimum signal at these wave-numbers if foregrounds can be removed.Figure 9 illustrates the possible contribution of the v bc sourced signal to the total at the EoR 21 cm minimum.The black-dashed line denotes the second-order density term in Equation 11 with b = 1 at z = 10. The blue (red)shaded regions denote the range of signal contributions from the sinks (sources) that we compute in § k = 0 .
06 h/Mpc and the dotted magenta (green) lines denote the thermal noise limits of SKA(HERA) at that wavenumber. At this wavenumber, the sinks term contributes 3 −
17% of the signal for an ionizedfraction of 13% depending on the Reionization history and assumed value of Γ − . The lower end of this range comesfrom physically realistic histories with Γ − = 0 .
3, and the high end comes from our “Best Case” scenario and shouldbe treated as an upper limit. A physically realistic history with Γ − = 3 . − and ionized more quickly than auniform-emissivity model would suggest (see Figure 3 of Nasir & D’Aloisio (2019)). On the other hand, it is unclearhow the effect of v bc on gas structures would be different in over-dense regions. A more detailed study would simulateclumping in over-dense regions and would account for a physically motivated relationship between z re and over-density.Despite these uncertainties, we suggest that a ∼
5% contribution to the signal from the sinks term is not unrealisticfor k = 0 .
06 h/Mpc, provided the linear order term in Equation 10 is close to 0. At k = 0 . ∼ .
06 h/Mpc, so we expect a ∼
1% contribution at thiswavenumber. Note that the range of percentages we find for the sinks term varies by only a factor of a few, whereasthe source term varies by 4 orders of magnitude for the range of f (cid:63) we consider. We therefore interpret the sinks termas a lower bound on the v bc -sourced signal, except in the very unlikely case that the source and sinks terms happento exactly cancel each other. 5. SummaryIn this work, we studied the impact of baryon-dark matter relative velocities on the small-scale clumpiness of theIGM during reionization, and how this effect impacts the EoR 21 cm signal. Although the streaming velocities weresmall ( ∼ . v bc was previously large relative to the sound speed. To quantify these effects, we usedhigh-resolution radiation hydrodynamics simulations that tracked the hydrodynamic response of the IGM to reioniza-tion. We found that the peak suppression of the clumpiness occurs within the first 5-10 Myr after the gas becomesmostly ionized, before the small-scale structure is erased by Jeans pressure smoothing of the gas. The clumpingfactor of ionized gas shows a peak suppression of 5 −
10% in regions that had streaming velocities of 30 km/s atrecombination (approximately the RMS value). Differences between regions with and without v bc fall to the percentlevel by ∆ t = 300 Myr, after the gas has had sufficient time to relax in response to the photo-heating from reionization.To quantify the impact of v bc on the EoR 21 cm power spectrum, we constructed a model for the signal that includesa term coupling P to fluctuations in v bc through a corresponding bias parameter. We modelled contributions to thisparameter from ionizing photon sinks and sources. Using our simulation results for the former, we found that thecontribution from sinks is relatively insensitive to the details of Reionization, as it is set mainly by the spectrum ofprimordial density fluctuations and pressure smoothing of the gas. We found that the characteristic BAO featureimprinted on P through coupling with the sinks is likely to appear at only the sub-percent level when P is at itsmaximum, roughly halfway through reionization. The feature is most pronounced at ≈
10% ionization, when P is ata minimum. At this time, the near cancellation of fluctuations in density and ionization allows power from higher-orderterms (i.e. from v bc ) to contribute more significantly. At the epoch of minimum P , we expect the BAO feature toappear at the 1% (5%) level at k ∼ . .
06) h/Mpc due to modulation of the sinks. The signal due to sources may9
Figure 9.
Components of the 21 cm brightness power spectrum at z = 10 when the signal is at its minimum during the EoR,using the same density bias parameters as in the upper right panel of Figure 6. We show the contribution from the second-orderdensity term in Equation 11 with b = 1 at z = 10 (black dashed). The blue (red) shaded region denotes the range of possiblecontributions from the sinks (sources) term discussed above at z = 10. The blue vertical line denotes k = 6 × − h/Mpc, andthe magenta (green) dashed lines denote the thermal noise limits of SKA (HERA) at that wavenumber. be larger than this, but it is subject to a large uncertainty because it depends on poorly-constrained source propertieslike the star formation efficiency. At these wave numbers, the minimum P that we estimate is close to the thermalnoise sensitivity limits of 21 cm experiments like SKA and HERA, so the prospect of detecting the signal in the nearfuture seems low. However, it may well be within the capability of the next generation of 21 cm instruments.This research was supported by HST award HST-AR15013.005-A. This work used the Extreme Science and En-gineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. We used the Bridges Large supercomputer at the Pittsburgh Supercomputing Center through allocationTG-AST120066. We also thank Jonathan Pober and Hyunbae Park for helpful comments on the draft version of thiswork. REFERENCES Ahn, K. 2016, ApJ, 830, 68Ahn, K., & Smith, B. D. 2018, ApJ, 869, 76Ali-Haimoud, Y., Meerburg, P. D., & Yuan, S. 2014,PhRvD, 89, 083506Barkana, R., & Loeb, A. 2001, PhR, 349, 125Becker, G., D’Aloisio, A., Davies, F. B., Hennawi, J. F., &Simcoe, R. A. 2019, BAAS, 51, 440 Becker, G. D., Bolton, J. S., Madau, P., et al. 2015,MNRAS, 447, 3402Bittner, J. M., & Loeb, A. 2011, arXiv e-prints,arXiv:1110.4659Blazek, J. A., McEwen, J. E., & Hirata, C. M. 2016,PhRvL, 116, 121303 Cain et al.
Bosman, S. E. I., Fan, X., Jiang, L., et al. 2018, MNRAS,479, 1055Chiou, Y. S., Naoz, S., Burkhart, B., Marinacci, F., &Vogelsberger, M. 2019, arXiv e-prints, arXiv:1904.08941Ciardi, B., Stoehr, F., & White, S. D. M. 2003, MonthlyNotices of the Royal Astronomical Society, 343, 1101.https://doi.org/10.1046/j.1365-8711.2003.06797.xCohen, A., Fialkov, A., & Barkana, R. 2016, MonthlyNotices of the Royal Astronomical Society: Letters, 459,L90. https://doi.org/10.1093/mnrasl/slw047Dalal, N., Pen, U.-L., & Seljak, U. 2010, JCAP, 2010, 007D’Aloisio, A., McQuinn, M., Davies, F. B., & Furlanetto,S. R. 2018, MNRAS, 473, 560D’Aloisio, A., McQuinn, M., Trac, H., Cain, C., &Mesinger, A. 2020, arXiv e-prints, arXiv:2002.02467Davies, F. B., Hennawi, J. F., Ba˜nados, E., et al. 2018, TheAstrophysical Journal, 864, 142.https://doi.org/10.3847%2F1538-4357%2Faad6dcDeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017,PASP, 129, 045001Eilers, A.-C., Davies, F. B., & Hennawi, J. F. 2018, ApJ,864, 53Emberson, J. D., Thomas, R. M., & Alvarez, M. A. 2013,The Astrophysical Journal, 763, 146. https://doi.org/10.1088%2F0004-637x%2F763%2F2%2F146Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132,117Fialkov, A. 2014, International Journal of Modern PhysicsD, 23, 1430017Fialkov, A., Barkana, R., Pinhas, A., & Visbal, E. 2014,MNRAS, 437, L36Fialkov, A., Barkana, R., Tseliakhovich, D., & Hirata,C. M. 2012, MNRAS, 424, 1335Furlanetto, S., Carilli, C. L., Mirocha, J., et al. 2019, arXive-prints, arXiv:1903.06204Furlanetto, S. R., Mirocha, J., Mebane, R. H., & Sun, G.2017, MNRAS, 472, 1576Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, PhR,433, 181Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004,The Astrophysical Journal, 613, 1Giri, S. K., D’Aloisio, A., Mellema, G., et al. 2019, Journalof Cosmology and Astro-Particle Physics, 2019, 058Greif, T. H., White, S. D. M., Klessen, R. S., & Springel,V. 2011, Apj, 736, 147Greig, B., Mesinger, A., & Koopmans, L. V. E. 2020,MNRAS, 491, 1398Hirata, C. M. 2018, MNRAS, 474, 2173Inoue, A. K., Hasegawa, K., Ishiyama, T., et al. 2018,PASJ, 70, 55 Irˇsiˇc, V., Xiao, H., & McQuinn, M. 2019, arXiv e-prints,arXiv:1911.11150Jensen, H., Majumdar, S., Mellema, G., et al. 2016,MNRAS, 456, 66Jensen, H., Datta, K. K., Mellema, G., et al. 2013,MNRAS, 435, 460Kashikawa, N., Shimasaku, K., Malkan, M. A., et al. 2006,ApJ, 648, 7Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, inAdvancing Astrophysics with the Square KilometreArray (AASKA14), 1Koopmans, L., Barkana, R., Bentum, M., et al. 2019, arXive-prints, arXiv:1908.04296Lanman, A. E., Pober, J. C., Kern, N. S., et al. 2019, arXive-prints, arXiv:1910.10573Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429Maio, U., Koopmans, L. V. E., & Ciardi, B. 2011, MonthlyNotices of the Royal Astronomical Society: Letters, 412,L40.https://dx.doi.org/10.1111/j.1745-3933.2010.01001.xMcGreer, I. D., Eftekharzadeh, S., Myers, A. D., & Fan, X.2016, AJ, 151, 61McQuinn, M., & D’Aloisio, A. 2018, JCAP, 2018, 016McQuinn, M., & O’Leary, R. M. 2012, The AstrophysicalJournal, 760, 3.https://doi.org/10.1088%2F0004-637x%2F760%2F1%2F3Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2018,MNRAS, 479, 4544Mesinger, A., Aykutalp, A., Vanzella, E., et al. 2015,MNRAS, 446, 566Mesinger, A., & Furlanetto, S. 2009, Monthly Notices of theRoyal Astronomical Society, 400, 1461.https://doi.org/10.1111/j.1365-2966.2009.15547.xMesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411,955Mu˜noz, J. B. 2019, PhRvD, 100, 063538Naoz, S., & Narayan, R. 2014, The Astrophysical Journal,791, L8. https://doi.org/10.1088%2F2041-8205%2F791%2F1%2Fl8Naoz, S., Yoshida, N., & Gnedin, N. Y. 2012, ApJ, 747, 128Nasir, F., & D’Aloisio, A. 2019, arXiv e-prints,arXiv:1910.03570O’Leary, R. M., & McQuinn, M. 2012, The AstrophysicalJournal, 760, 4.https://doi.org/10.1088%2F0004-637x%2F760%2F1%2F4Padmanabhan, T. 1993, Structure Formation in theUniverse, 499Parsons, A., Aguirre, J. E., Beardsley, A. P., et al. 2019, inBAAS, Vol. 51, 241 Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, TheAstrophysical Journal, 793, 113. https://doi.org/10.1088%2F0004-637x%2F793%2F2%2F113Planck Collaboration, Aghanim, N., Akrami, Y., et al.2018, arXiv e-prints, arXiv:1807.06209Pritchard, J. R., & Furlanetto, S. R. 2007, Monthly Noticesof the Royal Astronomical Society, 376, 1680.https://doi.org/10.1111/j.1365-2966.2007.11519.xPritchard, J. R., & Loeb, A. 2012, Reports on Progress inPhysics, 75, 086901. https://doi.org/10.1088%2F0034-4885%2F75%2F8%2F086901Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop,J. S. 2015, ApJL, 802, L19Schauer, A. T. P., Glover, S. C. O., Klessen, R. S., &Ceverino, D. 2019, MNRAS, 484, 3510Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ,744, 179Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417Tanaka, T. L., & Li, M. 2014, MNRAS, 439, 1092Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computingin Science & Engineering, 16, 62 Trac, H., & Cen, R. 2007, The Astrophysical Journal, 671,1. https://doi.org/10.1086%2F522566Trac, H., Cen, R., & Mansfield, P. 2015, ApJ, 813, 54Trac, H., & Pen, U.-L. 2004, NewA, 9, 443Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879Trott, C. M., & Pober, J. C. 2019, arXiv e-prints,arXiv:1909.12491Tseliakhovich, D., Barkana, R., & Hirata, C. M. 2011,MNRAS, 418, 906Tseliakhovich, D., & Hirata, C. 2010, Phys. Rev. D, 82,083520.https://link.aps.org/doi/10.1103/PhysRevD.82.083520Visbal, E., Bryan, G. L., & Haiman, Z. 2020, arXive-prints, arXiv:2001.11118Visbal, E., Haiman, Z., & Bryan, G. L. 2018, MNRAS, 475,5246Weinberger, L. H., Haehnelt, M. G., & Kulkarni, G. 2019,Monthly Notices of the Royal Astronomical Society, 485,1350. https://doi.org/10.1093/mnras/stz481Wouthuysen, S. A. 1952, AJ, 57, 31Wu, X., McQuinn, M., Kannan, R., et al. 2019, MonthlyNotices of the Royal Astronomical Society, 490, 3177.https://doi.org/10.1093/mnras/stz2807 Cain et al.
APPENDIXA. Bias Factor DerivationWe provide here the derivation for the sinks contribution to b x i ,v in Equation 16. Assuming x ri and C r R satisfyEquation 13, we can write a differential equation for δ rx i with the simple form˙ δ rx i + A ( t ) δ rx i = B ( t ) δ rC R (A1)where A ( t ) ≡ ˙ (cid:104) x i (cid:105)(cid:104) x i (cid:105) + α B n e (cid:104) C R (cid:105) B ( t ) ≡ − α B n e (cid:104) C R (cid:105) (A2)are functions only of the IGM mean values of x i and C R for which the solution to Equation 13 is known already.The fluctuation δ rC R is given by the average of local fluctuations over the distribution of z re and v bc within the patch,analogous to Equation 14. Assuming the right-hand-side of Equation A1 is a function of time only and not of δ rx i (which will be justified momentarily), the solution is δ rx i ( t ) = D ( t ) (cid:90) tt ( z ) dt (cid:48) F ( t (cid:48) ) δ rC R ( t (cid:48) ) (A3)where D ( t ) ≡ e − (cid:82) tt ( z dt (cid:48) A ( t (cid:48) ) F ( t ) ≡ B ( t (cid:48) ) e (cid:82) tt ( z dt (cid:48) A ( t (cid:48) ) (A4)where z is the redshift at which Reionization starts. Finally, we may write to first order in δ rv ≡ ( v − σ bc ) /σ bc , δ rC R = σ bc (cid:104) C R (cid:105) (cid:28) ∂C R ∂v (cid:12)(cid:12)(cid:12) v bc = σ bc (cid:29) δ rv + matter terms (A5)where the partial derivative is averaged as in Equation 14. Equation A5 is the statement that v bc is a biased tracer of C R . Combining Equations A3 and A5 yields δ rx i = b rx i ,v δ rv (A6)where b rx i ,v is the scale-dependent ionized fraction bias factor. Provided r is large enough that spatial fluctuations in P z re ( x i ( t )) are unimportant, we may write the scale-independent bias factor (Equation 16) as b x i ,v ( t ) = lim r →∞ b rx i ,v ( t ) = σ bc D ( t ) (cid:104) C R (cid:105) (cid:90) tt ( z ) dt (cid:48) F ( t (cid:48) ) S ( t (cid:48) ) (A7)where S ( t ) ≡ (cid:90) z ( t ) z dz re P z re ( x i ( t )) ∂C R ∂v (cid:12)(cid:12)(cid:12) v bc = σ bc ( z re , t ) (A8)Note that P v bc was absorbed in the definition of δ rv .B. Test of Initial ConditionsWe tested the initial conditions prescription used in this work by comparing the simulated matter power spectrum atvery high redshifts to the expectation from linear theory (LT). We did this primarily to verify that v bc is implementedcorrectly in our simulations, but also to confirm that starting from z = 1080 produces correct results. To do this, weran a set of hydro-only test simulations down to z = 30, which are listed in Table 1. These simulations were initializedat z start , which is either 1080 (as in our production runs) or at z = 300 (as in D20).The matter power spectrum for the tests starting from z = 1080 are shown in Figure A1 at redshifts 270, 145, 68,and 45. The top (bottom) set of curves show the DM (baryon) power spectrum. The solid blue (black) curves arethe LT predictions evolved from redshift 1080 CAMB TFs using the LT approximation from O’Leary & McQuinn3 Simulation z init N L (Mpc /h ) v bc (km / s)Fiducial 1080 256 .
256 0Fiducial + v bc .
256 30High Res 1080 512 .
256 0High Res + v bc .
256 30Low z 300 256 .
256 0
Table 1.
List of simulations run to test the accuracy of the initial conditions setup used in this work. Here, z init is theinitialization redshift of the simulation. (2012) and employed in their initial conditions code CICsASS. Runs with and without v bc have indistinguishable DMpower spectra, while the baryon power spectrum is suppressed significantly in the cases with v bc . In all cases, thesimulations agree well with the LT expectation until z = 45 when non-linear effects begin to become important. Thehigher resolution runs do a better job at small scales, as expected. Importantly, the simulations with v bc reproducethe CICsASS prediction very well at scales that are captured by the simulations. These results demonstrate that v bc is implemented correctly in our simulations.We also checked how our results are affected by using different starting redshifts. In Figure A2, we plot the Fiducial(cyan dashed), High res (magenta dotted), Low z(red dashed) power spectra. We compare these to the CICsASSLT expectation without v bc (black solid curve). We see that the simulations initialized at z = 1080 (Fiducial andHigh Res) agree well the LT approximation. The one started from the z = 300 CAMB TF deviates slightly fromthe other two initially. However, after some time has passed, the relative difference decreases, indicating that theslight difference in initial conditions does not affect the results significantly at much later times. This is impor-tant for us because it indicates that we can reasonably compare our results to the simulations in D20 (which wereinitialized in the same was as the Low z run). It also demonstrates that initializing simulations at z = 1080 doesnot introduce significant shot noise, as has been suggested by many previous authors (e.g. O’Leary & McQuinn (2012)).C. Effects of Resolution and Box SizeWe also assessed the sensitivity of our clumping factor results to numerical resolution and box size. This is importantbecause v bc impacts small-scale gas structures appreciably but leaves the larger structures unaffected. This suggeststhat too-small boxes would fail to capture the large-scale structures, producing an over-estimate of the v bc effect.Conversely, large boxes with poor resolution would fail to resolve the structures that are most affected, leading toan under-estimate. We quantified these differences by running a set of simulations with N = 512 and L = 256 kpc,which gives 8 times the resolution and 1 / v bc = 0 , − = 0 . , . z re = 6. We compared these results to our full box-size run with v bc = 65 km/s, z re = 6,and Γ − = 0 .
3. We see a maximum suppression in C R relative to the no- v bc case with the same parameters of ∼ − = 0 . . z re = 12, but this time varying the resolution and box size one at a time. We found that increasingbox size at fixed resolution reduces the relative v bc effect after about ∆ t ∼ t (cid:47)
10 Myr but not much after this. These results are consistent with the picture thatsmall structures that are affected by v bc dominate the recombination rate early, but after relaxation is complete therecombination rate is set by larger structures that are not appreciably affected by v bc .In Figure A3, we plot the number of hydrogen recombinations per hydrogen atom since z re for the convergence testsat z re = 6 alongside our production runs (the fiducial case) with ( z re , Γ − ) = (6 , . v bc = 0 (65 km/s). Thedifference between the runs with and without v bc increases for smaller box size/higher resolution and increasing Γ − ,suggesting that the systems that are resolved in those simulations are more strongly impacted by v bc . In addition,the number of recombinations is higher at later times in fiducial case, suggesting that large structures not captured inthe smaller simulations contribute a large fraction of the recombinations. This result confirms our suspicion that boxsizes that are too small to capture a representative sample of absorbing systems will over-estimate the importance of4 Cain et al.
Figure A1.
Power spectra of baryons (bottom curves) and DM (top curves) for Fiducial (cyan dashed), Fiducial + v bc (yellowdashed), High Res (magenta dotted), High Res + v bc (greed dotted) at redshifts 270, 145, 68, and 45 compared to the LTexpectation from CICsASS. The blue (black) solid curves are the CICsASS LT approximation with(out) v bc . All the DM curvesare indistinguishable, and the simulations with and without v bc agree well with their respective LT predictions, especially whenthe resolution is increased. v bc . However, it may be that some of the difference comes from the additional resolution these boxes, in which caseour fiducial runs may slightly under-estimate v bc ’s importance in patches that have been recently ionized.5 Figure A2.
Baryon and CDM power spectra showing how our results vary with different initialization schemes. The simulationsshown are Fiducial (cyan dashed), High Res (magenta dotted), and Low z + CAMB (red dashed). The Low z run deviatesslightly from the others initially, but all three converge at lower redshifts. This implies that starting from z = 1080 and 300give very similar results especially at low redshifts. Cain et al.
ΔtΔ[Myr] R e c o m b i n a i o n s Δ p e r Δ H Δ a o m v bc = 0ΔFiducialv bc = 65ΔFiducialv bc = 0ΔΓ −12 = 0.3v bc = 65ΔΓ −12 = 0.3v bc = 0ΔΓ −12 = 9.2v bc = 65ΔΓ −12 = 9.2 Figure A3.
Number of recombinations per hydrogen atom for our convergence runs (dotted/dashed curves) and our fiducial z re = 6, Γ − = 0 . v bc is much larger in the smaller boxes, especially the high Γ − case. This is likely due to a combination of the lack of large systems and better resolution of small ones, both of which enhancethe importance of v bcbc