Stacking Redshifted 21cm Images of HII Regions Around High Redshift Galaxies as a Probe of Early Reionization
James E Davies, Rupert A. Croft, Tiziana Di-Matteo, Bradley Greig, Yu Feng, Stuart Wyithe
MMNRAS , 1–11 (2019) Preprint 28 July 2020 Compiled using MNRAS L A TEX style file v3.0
Stacking Redshifted 21cm Images of HII Regions AroundHigh Redshift Galaxies as a Probe of Early Reionization
James E. Davies, , (cid:63) Rupert A. Croft, , , Tiziana Di-Matteo, , , Bradley Greig, , Yu Feng, and J. Stuart B. Wyithe , School of Physics, The University of Melbourne, Parkville, Victoria 3010, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) McWilliams Center for Cosmology, Carnegie Mellon University, Pittsburgh PA, 15213 Berkeley Center for Cosmological Physics, University of California at Berkeley, Berkeley, CA, 94720, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
A number of current and future experiments aim to detect the reionization of neu-tral hydrogen by the first stars and galaxies in the Universe via the redshifted 21cmline. Using the
BlueTides simulation, we investigate the measurement of an aver-age ionised region towards the beginning of reionization by stacking redshifted 21cmimages around optically identified bright galaxies using mock observations. We findthat with an SKA 1000 hour observation, assuming perfect foreground subtraction, a σ detection of a stacked HII region can be made with 30 images around some of thebrightest galaxies in bluetides (brighter than M UV < − . ) at z = (correspond-ing to a neutral fraction of 90.1 % in our model). If the HII regions are detected withhigh certainty, we can recover simulated relationships between the UV magnitude ofgalaxies and the sizes of the ionised regions they reside in. These mock observationscan also distinguish between scenarios where the IGM is in net emission or absorptionof 21cm photons. Once 21cm foreground contamination is included, we find that evenwith up to 200 images around these rare, bright galaxies, only a tentative > σ detec-tion will be possible. However, partial foreground subtraction substantially improvessignal-to-noise. For example, we predict that reducing the area of Fourier space dom-inated by foregrounds by 50 (80) percent will allow > σ ( > σ ) detections of ionisedregions at z = . Key words: dark ages, reionization, first stars – intergalactic medium – early universe
The reionization of cosmic hydrogen is believed to have orig-inated in the densest regions in the Universe, where brightstar-forming galaxies were more highly clustered. These firstionised regions are predicted to grow and overlap, filling theIGM by z ∼ (e.g: Greig & Mesinger (2017)).Radio observations of the redshifted 21cm line from thehyperfine splitting of atomic hydrogen currently present themost promising future probe of the epoch of reionization.Radio interferometers aim to detect the EoR in one of twoways. Most first generation low frequency radio telescopessuch as the Murchison Widefield Array (MWA) and theLow Frequency Array (LOFAR) , along with the more sen-sitive second generation telescope Hydrogen Epoch of Reion- (cid:63) E-mail: [email protected] ization Array (HERA) aim to detect the epoch of reioniza-tion statistically by measuring the 21cm power spectrum.Another second generation telescope, the Square KilometreArray (SKA) , also aims to image reionization directly in21cm emission (Mellema et al. 2013). In this paper we fo-cus on images of ionised regions around bright galaxies athigh redshift. A number of works have considered identifyingHII regions around individual sources using matched filtersand imaging (Datta et al. 2007, 2008; Majumdar et al. 2012;Ghara et al. 2017; Ghara & Choudhury 2019; Hassan et al.2019), and similar methods using high-redshift quasars wereproposed in Wyithe et al. (2005); Kohler et al. (2005); Geil &Wyithe (2008); Ma et al. (2020). Since these methods workwith single images, they require relatively large bubbles, andnearly perfect subtraction of radio foregrounds in order to https://reionization.org/ © a r X i v : . [ a s t r o - ph . C O ] J u l James E. Davies detect the EoR. Thus, imaging an individual HII region to-ward the beginning of reionization will be difficult. How-ever, since galaxies are much more common than quasars,we can stack 21cm images around selected bright galaxies, inanalogy to high-z HI detections (Blyth et al. 2015), takingadvantage of synergies between 21cm and optical-infraredobservations in order to average out noise and boost thesignal, facilitating observations of smaller HII regions.Simulations suggest that the brightest galaxies at highredshift will likely sit in the first ionised regions. Not only dothese galaxies produce a large number of ionising photons,but they are also more likely to be inside the dense, clusteredfilaments of the cosmic web, with many neighboring galax-ies and quasars which also emit ionising photons (Geil et al.2017). 21cm images around these objects would show anionised region approximately centered on the bright galaxy.If many of these images were stacked, we would see an av-erage ionised region around the brightest galaxies. Thesestacks will provide an alternative detection strategy whichcould increase the sensitivity of observations, allowing us toprobe further into the EoR and understand the nature ofthe galaxies that reionised the Universe.In this work, we investigate stacking of 21cm intensitymaps around bright galaxies using the reionization historyand galaxy population from the
Bluetides cosmologicalsimulation. We draw connections between the appearance ofthese stacked profiles, the properties of the sources they arecentred on, and the distribution of ionised hydrogen in theearly Universe. Using these relations from the simulation,we predict the signal for upcoming galaxy samples from the
Wide-Field Infrared Survey Telescope (WFIRST)
High Lati-tude Survey (HLS) and observations of the redshifted cosmo-logical 21cm signal from the low-frequency Square KilometerArray (SKA1-low).This strategy was proposed in Geil et al. (2017), whereanalytic spherical HII regions were constructed from sim-ulated luminosity functions and simulated relationships be-tween HII region size and UV luminosity. We expand on thisidea in this work by directly stacking 21cm intensity maps in
BlueTides , thereby including the effects of asymmetric HIIregions, galaxies laying off-center relative to the HII regionsthey reside in, as well as the effect of radio foregrounds.This paper is laid out as follows: Section 2 covers thedetails of the
BlueTides simulation, as well as the method-ology for calculating the 21cm brightness temperature pro-files and stacking them. Section 3 covers the details of ourmock observations, including number of available galaxiesand SKA1-low thermal noise. Section 4 shows our results,including how average HII regions are detected, inferenceof the mean bubble radius from stacked images, and corre-lations with galaxy luminosity and spin temperature. Sec-tion 5 explores the effects of observational limitations onthis method by including the effect of radio foregrounds. Weconclude in Section 6.
The
Bluetides simulation (Feng et al. 2016) is a cos-mological hydrodynamic simulation within a cube of sidelength h − comoving Mpc. It contains × darkmatter and baryonic particles of mass . × h − M (cid:12) and − − − − − − − UV Magnitude − − − − Φ ( h c M p c − m ag − ) z = 11.0z = 10.0z = 9.0z = 8.0 Figure 1.
UV Magnitude function predictions from
Bluetides .Luminosity functions are shown at redshifts z =
11, 10, 9 and 8. . × h − M (cid:12) , respectively. Particle snapshots are avail-able from z = down to z = (Marshall et al. 2019; Niet al. 2019). Bluetides uses a pressure-entropy formulation ofSmoothed Particle Hydrodynamics (Read et al. 2010; Hop-kins 2013). This includes a multi-phase star formation model(Springel & Hernquist 2003; Vogelsberger et al. 2013), cool-ing via radiative processes (Katz et al. 1996), metal cool-ing (Vogelsberger et al. 2014a), the effects of molecularhydrogen on star formation (Vogelsberger et al. 2014b),a type II supernovae feedback model (Okamoto et al.2010), and a super-massive black hole model (Di Mat-teo et al. 2005).
Bluetides is the largest high-redshifthydrodynamic simulation run to date. The large volumeallows us to have statistically significant samples of ob-servable galaxies at the beginning of reionization. TheWMAP9 cosmology (Hinshaw et al. 2013) is used in blue-tides and throughout this work: { Ω m , Ω b , Ω Λ , h , σ , n s } = { . , . , . , . , . , . } .In this work, we investigate how the properties of thebrightest galaxies correlate with the structure of the 21cmsignal. The UV magnitude distributions for these galaxies in BlueTides at redshifts 11, 10, 9, and 8 are shown in Figure1. These have been shown by Feng et al. (2016), Waters et al.(2016), and Wilkins et al. (2017) to be consistent with cur-rent observations. These distributions are used to estimatethe number of galaxies per SKA field which are availablefor stacking. The neutral fraction x HI is set by the paramet-ric reionization model used in BlueTides (Battaglia et al.2013), which takes as input the density field at the chosenmedian redshift of reionisation ¯ z = (consistent with con-straints from Greig & Mesinger 2017). A bias function isapplied to the density field, which accounts for the correla-tion between the redshift of reionisaiton and density fieldson different scales, calibrated to hydrodynamic and radiativetransfer simulations. This process results in reionization red-shift ( z reion ) field in a ( − cMpc resolution) grid for thesimulation. The neutral fraction in each cell ( x HI ) at a par-ticular redshift z is set to unity in cells where z reion < z , and MNRAS000
Bluetides is the largest high-redshifthydrodynamic simulation run to date. The large volumeallows us to have statistically significant samples of ob-servable galaxies at the beginning of reionization. TheWMAP9 cosmology (Hinshaw et al. 2013) is used in blue-tides and throughout this work: { Ω m , Ω b , Ω Λ , h , σ , n s } = { . , . , . , . , . , . } .In this work, we investigate how the properties of thebrightest galaxies correlate with the structure of the 21cmsignal. The UV magnitude distributions for these galaxies in BlueTides at redshifts 11, 10, 9, and 8 are shown in Figure1. These have been shown by Feng et al. (2016), Waters et al.(2016), and Wilkins et al. (2017) to be consistent with cur-rent observations. These distributions are used to estimatethe number of galaxies per SKA field which are availablefor stacking. The neutral fraction x HI is set by the paramet-ric reionization model used in BlueTides (Battaglia et al.2013), which takes as input the density field at the chosenmedian redshift of reionisation ¯ z = (consistent with con-straints from Greig & Mesinger 2017). A bias function isapplied to the density field, which accounts for the correla-tion between the redshift of reionisaiton and density fieldson different scales, calibrated to hydrodynamic and radiativetransfer simulations. This process results in reionization red-shift ( z reion ) field in a ( − cMpc resolution) grid for thesimulation. The neutral fraction in each cell ( x HI ) at a par-ticular redshift z is set to unity in cells where z reion < z , and MNRAS000 , 1–11 (2019) tacking 21cm Images Around High Redshift Galaxies zero otherwise. The global neutral fraction and topology ofreionization in BlueTides are shown in Figure 2.This method does not directly take into account sourceproperties, which may contribute to the positions of thebrightest galaxies within HII regions, and the asymmetryof those regions. It is possible that a model that computesionisation structure directly from emitted ionising photonswould result in early HII regions being more spherical, andcentred closer to the brightest galaxies within them (e.g:Geil et al. 2017), making them easier to detect in stackedimages.This reionization grid is applied in post-processing,meaning inhomogeneous feedback from reionization is notincluded. However, since we are studying the earliest phasesof reionization, around the brightest galaxies in the simula-tion which are well above the Jeans mass, we do not expectreionization feedback to have a large effect. Estimated bub-ble sizes for the brightest 15 galaxies at redshift 9, where theglobal volume-weighted neutral fraction is 90.1%, are shownin Figure 3. Although these are the brightest galaxies in thesimulation, they do not necessarily reside in the centre of thelargest HII regions in the simulation. Many of these brightgalaxies are off-centre in these regions, and many have largerHII regions nearby that the galaxy itself is not a part of.We calculate an estimate of the true size of HII regionsin the simulation via a ray-tracing method, taking the av-erage distance from each galaxy to the nearest ionisationboundary (where X HI < . ) in randomly selected di-rections (e.g: Geil et al. 2017). The distribution of HII regionsizes, calculated via this ray-tracing method, plotted againstthe UV magnitude of the brightest galaxy within them isshown in Figure 4.We calculate the 21cm brightness temperature in the Bluetides volume using (Furlanetto et al. 2006): δ T b , = K (cid:18) Ω b h . (cid:19) (cid:18) . Ω m h + z (cid:19) / (cid:18) − T CMB T s (cid:19) ( + δ ) x HI (1)where Ω b and Ω m are the abundances of baryons and mat-ter respectively, T s is the 21cm spin temperature, T CMB = . ( + z ) is the temperature of the cosmic microwave back-ground, and x HI is the local neutral fraction.Unless otherwise stated, we assume the IGM is ‘satu-rated’, where the spin temperature term in the above equa-tion T C M B T s ≈ , and the IGM is always seen in net emission of21cm radiation. In section 4.4, we examine an ‘unsaturated’case which follows the global spin temperature evolution ofthe ‘bright galaxies’ scenario from the Evolution of 21cmStructure (EOS) simulations (Mesinger et al. 2016), whichhas a similar neutral fraction evolution to our reionizationmodel. Using this spin temperature evolution, at redshifts11,10, 9 and 8 the brightness temperature of neutral hy-drogen at mean density is -98.6, -41.6, -1.06 and 18.7 mKrespectively. We centre mock δ T b observations around galaxies in Blue-Tides . The centre point of each image is slightly displaced from the galaxy along the line-of-sight to model redshift un-certainties (see section 3.1). Each grid cell is then assigneda redshift depending on its distance from this point alongthe line-of-sight direction, where the central point is alwaysat the redshift of the simulation snapshot. 21cm grids arethen stacked in × ×
80 h − Mpc sub-cubes, centred onthe closest cell to the central point, so that grid cells arealigned.While the reionization model produces asymmetric HIIregions, even towards the beginning of reionization, we ex-pect a signal that is roughly cylindrically symmetric afterstacking, since the amount and direction of asymmetry willbe uncorrelated between each HII region. The right panel offigure 3 shows the HII profile of a region stacked around thebrightest 15 galaxies in
BlueTides . Along with the asym-metry of the ionised regions, the mis-centreing of brightgalaxies within HII regions is the largest source of errorwhen inferring statistics of the bubble size distribution fromstacked profiles. Both of these effects will smooth out thestacked profiles depending on the average amount of asym-metry or the average distance of each galaxy from the centreof its HII region. However, the stacked profile itself will stillbe approximately symmetrical.
The signal for stacked 21cm images will be sensitive tothe number of visible galaxies in the SKA field. The num-ber of galaxies that we would expect to find in a surveycan be calculated from the bluetides luminosity functionand the SKA survey volume, determined by central red-shift z , depth ∆ z (such that the redshift boundaries areat z ± ∆ z / ), and the field of view Ω SK A = π (cid:16) λ D (cid:17) . Here λ is the observed wavelength of rest-frame 21cm radiationat the redshift of observation, D is the station diameterof 35m. At redshifts z = { , , , } , the field of view is Ω SK A ≈ { . , . , . , . } square degrees respectively.The predicted number of galaxies brighter than a lim-iting UV magnitude M cut at redshift z in the SKA fieldof view with depth ∆ z is shown in Figure 5. At redshift9 and with a depth of ∆ z = , we would expect to have1112 galaxies brighter than M cut = − . in the vol-ume. Magnitude cuts are chosen based on the WFIRSTHLS σ limit at the redshift upper bound in the chosenvolume, and those cuts at redshifts z = { , , , } are M cut = {− . , − . , − . , − . } respectively. At red-shifts z = { , , , } the number of galaxies above the mag-nitude cut in the bluetides volume that we sample fromis 132, 633, 3241, and 12837, and the number of sampledgalaxies within a SKA field of view is 54, 239, 1112, and3936 respectively.The instrument used to observe these galaxies is also im-portant in providing accurate locations, particularly alongthe line-of-sight. Uncertainty in the positions of brightgalaxies smooths out brightness temperature profiles whenstacked, resulting in shallower signals that are harder to de-tect. WFIRST grism spectroscopy will have an associatedredshift error of σ z = . ( + z ) (Spergel et al. 2015). Atredshift 9, this results in a line-of-sight position uncertainty MNRAS , 1–11 (2019)
James E. Davies x(h − cMpc) y ( h − c M p c ) Redshift . . . . . . X H I Figure 2.
Redshift of reionization slice (left) and global neutral fraction (right) in the
Bluetides reionization scenario we use. The EoRhas a median redshift of z = , and a duration in redshift of roughly ∆ z = . Figure 3.
The brightest 15 galaxies at redshift 9, x H I = . , the surrounding ionisation fields, and estimated individual bubble sizesfrom ray-tracing (shown as circles centred on each galaxy). The dimensions of each slice are 40 x 40 x 5 h − Mpc . The color representsthe average volume-weighted neutral fraction along the 5 h − Mpc depth of each slice. We also show the stacked profile of these 15 ionisedregions in the rightmost panel, illustrating how stacking can smooth over the asymmetry of individual HII regions of σ d = .
56 h − cMpc . Angular position uncertainty is as-sumed to be negligible. In order to create mock observations, we use the package (Pober et al. 2013, 2014) to create realisationsof thermal noise for the SKA1-low. takes theantenna configuration of SKA1-low , along with the system https://github.com/jpober/21cmSense The antenna layout used is the SKA1-low Configuration V4a temperature and dish size, to calculate sensitivity to line-of-sight and perpendicular modes at a certain frequency.We use the 21cm rest-frame frequency at the snapshot red-shift as our central frequency, and generate 2D sensitivitypower spectra using . Assuming cylindrical sym-metry, we interpolate the 2D power spectra to the samemodes sampled by our brightness temperature grids, gener-ating fluctuation amplitudes for each point in Fourier spacefor a 1080 hour observation (more detailed parameters usedin the package can be found in Table 1). This is then added,with a random phase for each mode (assuming uncorrelatednoise between modes), to the Fourier transform of the bright-
MNRAS000
MNRAS000 , 1–11 (2019) tacking 21cm Images Around High Redshift Galaxies − . − . − . − . − . − . − . − . − . UV Magnitude B ubb l e R a d i u s ( h − M p c ) Figure 4.
HII region radius distributions, binned by the UVmagnitude of the brightest galaxy within them, at redshift 9.The solid black lines show mean bubble radius in each bin, andthe dotted lines show the 16th and 84th percentiles, such thatthe central region contains 67% of HII regions in each bin. Wesee that dimmer galaxies generally host smaller ionised regions.There is also an increasing population of galaxies outside a HIIregion (with zero bubble radius) as UV magnitude increases. − − − − − M c u t z = 11 z = 10 . . . . ∆z − − − − − M c u t z = 9 . . . . ∆zz = 8 Figure 5.
Number of expected galaxies brighter than a UV mag-nitude M cut within a volume of a survey with the SKA field ofview and depth ∆ z . Contour lines are drawn at 1, 10, 30, 100,200, 500, 1000, and 5000 galaxies. ness temperature grids from BlueTides to simulate thermalnoise from SKA1-low.Regions of Fourier space in the simulation that are notprobed by the instrument are discarded in both the imageand noise cubes, including the k = modes. We make ad-ditional cuts in Fourier space above k ⊥ = .
5h cMpc − and k (cid:107) = .
5h cMpc − . This acts as a low-pass filter which re-moves the higher small scale noise while retaining most ofour signal, which will usually be several h cMpc − wide due Parameter ValueIntegration time
Observing days 180Hours per day 6Station diameter 35mReciever temperature 50KNumber of Antennas 513Longest baseline used 10km
Table 1.
Observing parameters to calculate the thermal noise forSKA1-low with . . . . . . . k ⊥ (h cMpc − ) . . . . . k || ( h c M p c − ) mK Figure 6. − h . Produced by the package at redshift9. to the size of the HII regions. This filter smooths over bothperpendicular and line-of-sight modes, since we expect thesignal to be roughly symmetrical, however in practice takingdifferent k (cid:107) cuts makes little difference to our results. Resultsin Section 4 assume perfect subtraction of radio foregrounds.In Section 5 we show results including model foregrounds byexcluding areas of Fourier space below a certain line (e.g:Hassan et al. 2019). Figure 6 shows the 2D sensitivity powerspectrum under the parameters laid out above, before cutsto Fourier space are made. We begin by considering the stacks of brightness tempera-ture fields around up to 128 galaxies in the
Bluetides sim-ulation. We build a sample of galaxies by first estimating thenumber of galaxies brighter than the WFIRST HLS σ UVmagnitude limit ( M cut = − . ) , in a volume determined bya central redshift z = with depth ∆ z = , and the field ofview for the Square Kilometer Array (9.28 square degrees at z = ). This results in 1112 galaxies. This number of galaxiesare sampled from the Bluetides galaxies at the central red-shift and above the same magnitude cut with replacement,
MNRAS , 1–11 (2019)
James E. Davies
Brightest Brightest 10 Brightest 30 − − − −
10 0 10 20 30 40 x k ( h − cMpc) − − δ T b ( m K ) Brightest Brightest 10 Brightest 30 − − m K Figure 7.
Top: × × − × − column. representing a sample of available galaxies in an SKA1-lowfield. The regions around up to the brightest 128 of thesegalaxies are then selected as our sample for stacking. We takea random rotation and reflection along two axes, and keepthe line-of-sight axis constant for the purposes of adding red-shift uncertainty (this gives us 8 possible configurations pergalaxy in Bluetides ). The top panel of Figure 7 shows acomparison of the brightness temperature around the bright-est galaxy in one sample compared to a stack of the 10 and 30brightest galaxies in the same sample. Given noise, redshifterrors, asymmetry and mis-centreing of the ionised regionas shown in figure 3, distinguishing a single ionised regionfrom the background could be difficult for some galaxy sam-ples. However stacking many regions diminishes backgroundfluctuations, asymmetries and noise, resulting in an averagebubble profile which is roughly Gaussian. As we increase thenumber of galaxies stacked, background variation decreasesand the average profile converges. The bottom panel of Fig-ure 7 shows the decrease in background noise and fluctua-tions as we stack greater numbers of images along the lineof sight spectra averaged within a 5 h − cMpc wide aperture. The average ionised region profile is fit to a 3D Gaussianwith separate widths along the line-of-sight and sky plane,using a Monte-Carlo Markov Chain fitting method of theform, f ( ρ, z ) = Ae − (cid:18) ρ σ ρ + z σ z (cid:19) , (2) − − − A m p li t ud e ( m K ) . . . . σ r ( c M p c h − ) Number of Images . . . . S N R Figure 8.
Convergence of the Gaussian amplitude (top) andwidth (middle) versus number of stacked images, the amplitudebecomes inconsistent with random pointings, where amplitude isfit closely to zero (see Appendix A). The SNR (bottom) reaches5 at ∼ stacked images. Errorbars are 1 σ . where ρ is the transverse distance from the centre of theimage, and z is the line-of-sight distance.We choose A as the parameter describing the depth (orheight) of the profile, and σ ρ as the parameter describing itswidth. We require separate widths along the two dimensionsin order to fit the profiles because redshift uncertainty willelongate the profile along the line-of-sight. We note that σ z isnot truly a free parameter, as it results from the convolutionof a spherically symmetric Gaussian with the redshift uncer-tainty. However, we take the conservative approach of allow-ing it to take any value in our fits. It also becomes necessaryto fit separately to σ z when modelling foreground contami-nation, where the small-scale sky-plane structure and large-scale line-of sight structure is removed, changing the shapeof the profiles. Since σ z contains no new information aboutthe underlying HII structure however, it is not used for infer-ence in this work. We do not fit to a background brightnesstemperature, since in our treatment of instrument noise, weremove the zero mode in Fourier space, so that the imageeffectively has zero mean.We can look at the effect of stacking more quantita-tively, in terms of the fit parameters ( A , σ ρ ) and bootstrapsampling of the image stacks. We take 128 bootstrap samplesof the same size, giving us different realisations of availablegalaxies, HII regions, and instrumental noise. We choose thedepth of the stacked absorption to measure the strength ofthe signal. Our Monte-Carlo fitting process gives us a mean µ A and standard deviation σ A in depth. For comparion, weperform the same analysis on stacked images of random lo-cations in the BlueTides volume. We find that amplitudesare distributed closely around zero, and widths are widelydistributed around the scale of density fluctuations and in-
MNRAS000
MNRAS000 , 1–11 (2019) tacking 21cm Images Around High Redshift Galaxies strument noise (see appendix A for distributions of param-eter for random pointings). Since this wide range includesthe width of many HII regions, we only use the amplitudeas a metric of detection.We calculate the signal to noise ratio (SNR) by com-paring the distributions of fitted depths between the stackedimages of galaxies and random pointings. We find that bothdistributions are roughly Gaussian, so we characterise themby their means µ A , µ R and standard deviations σ A , σ R wherethe subscript A refers to the galaxy case, and R to the ran-dom pointings.The signal to noise radio in this case is the mean of thedistribution of differences in fit amplitude between galaxyand random pointings, divided by the standard deviation ofthis difference: SNR = | µ diff | σ diff ≈ | µ R − µ A | (cid:113) σ + σ . (3)Figure 8 shows the convergence of our Gaussian fit param-eters as more images are stacked. The stacked images reacha SNR greater than at images.We find that the SNR increases more slowly as the num-ber of stacked images increases. This is because the smallerHII regions around dimmer galaxies are being added to thestack, resulting in a slightly shallower and narrower pro-file. Eventually, the SNR gain by reducing background fluc-tuations will be overcome by the decrease in fit amplitudecaused by smaller ionised bubbles around dimmer galaxiesbeing added to the stack. However, with <
100 galaxies, wedo not reach the turning point in SNR.
We next demonstrate that we can recover the
BlueTides
Bubble size to UV magnitude correlations by performingstacks of 30 galaxies using the same bootstrap samplingmethod as section 4.1. We then compare the properties ofthe resulting Gaussian fits with the true bubble sizes inthe simulation. We include samples from redshifts z = ( X HI = . , M cut = . ) , z = ( X HI = . , M cut = − . ) z = ( X HI = . , M cut = − . ) and z = ( X HI = , M cut = − . ) to test different stages of theEoR. Each stack of 30 images has a single fit depth andwidth, as well as a true mean bubble radius calculated fromthe ray tracing method described in section 2. The correla-tions between these quantities are plotted in Figure 9.As shown by the z = results, this method is bestapplied before the overlap phase of the EoR, when higheramounts of asymmetry smooth out the profiles, and the re-lationship between the width and depth of the stacks andbubble size distribution flattens. Similarly, if we stack at z = where bright galaxies do not necessarily reside in largeionised regions, the stacked image is dominated by the denseneutral regions surrounding the galaxies rather than ionisedbubbles, resulting in a fit approximately ∼ cMpch − widewith a positive height of ∼ + mK.At redshifts 9 and 10, where the fit width is linearlyprorportional to bubble radius, fitting a simple linear modelto this relationship results in a slope of σ ρ | R | ≈ . . Theslope within each redshift group is flattened slightly by therandom sampling, since each galaxy sample contains HII re- F i t W i d t h ( h − c M p c ) z = 11 z = 10 z = 9 z = 8 Mean Bubble Radius (h − cMpc) − − − F i t A m p li t ud e ( m K ) Figure 9.
Gaussian fit parameters versus mean HII region radiusfor various galaxy samples at z = { , , , } within BlueTides .Top: Gaussian fit amplitude. Bottom: σ Gaussian width. Thereare clear relationships between the width and depth of the fitsand the bubble sizes. The best fit linear realationship betweenbubble radius and fit width at redshifts z = and z = is shownas a dotted line. gions from a wider range of radii that are averaged withineach stack before fitting. While the depth of the fits is corre-lated with bubble radius, it is explicitly dependent on IGMdensity, spin temperature, and redshift (see equation 1). Be-cause of this, we only present the relation between fit widthand mean bubble radius. This section details the correlations between galaxy magni-tudes, and the resulting stacked 21cm brightness temper-ature profiles. We bin galaxies based on UV magnitude,and create stacked 21cm observations within each bin us-ing the method discussed in section 2.1. The resulting radiiand depth of the profiles at redshift 9 are compared withgalaxy UV magnitude in Figure 10.We find that fitted profiles have lower depths and widthsfor galaxy samples with lower luminosity recovering, thetrends in Figures 4 and 9. On average, galaxies of abso-lute UV magnitude -24 (-21) are hosted in ionised regionsof radius 5.8 (4.1) Mpc which, when stacked with simulatedSKA1-low thermal noise for a 1080h observation, are fit toGaussian profiles of depth -12 (-9) mK and ( σ ) width 4.7(4.0) h − cMpc respectively.The correlations between galaxy properties and HII re-gions imply that stacked 21cm observations could be used toprobe the contribution of different galaxies to the reioniza-tion of the Universe. However, distinguishing between reion-ization models using stacked 21cm images will require suf- MNRAS , 1–11 (2019)
James E. Davies
Muv [-25.] Muv [-23.] Muv [-21.] − . − . − . − . . . . . . m K − − − − − r a d i u s ( σ )( c M p c ) − − − − − UV Magnitude − − d e p t h ( m K ) Figure 10.
Top: 21cm brightness temperature slices of 30 stackedgalaxies, where the samples were chosen based on UV magni-tude. Dimmer galaxies reside in visibly smaller HII regions. Bot-tom: Gaussian fit parameters for binned galaxy samples. Brightergalaxies result in significantly deeper and slightly wider stackedprofiles. Vertical errorbars are σ variation in the fit parameter,horizontal errorbars cover the magnitude bins, with points at theaverage UV magnitude in each bin. ficiently high SNR detections, so that the properties of theimages correlate with the statistics of the EoR, and not withnoise in the image. The spin temperature evolution is determined by the his-tory of the first sources of Lyman alpha, UV and X-rayphotons. In this section we show that stacked 21cm imagescould distinguish between positive and negative brightnesstemperature cases. We test two models of spin temperature.In the first model, the IGM spin temperature is saturated( T s >> T CMB in equation (1)), and the IGM is emitting atits maximum value. The second model follows the ’BrightGalaxies’ case presented in Mesinger et al. (2016), wherethe IGM is seen in absorption for z ≥ .In Figure 11, we show the slices and fit parameters ineach spin temperature case, for up to 100 stacked imagesat redshift 10, using SKA1-low thermal noise described inSection 3.2. In the unsaturated case, the IGM is in net ab-sorption of 21cm photons, and the brightness temperature ofneutral hydrogen at mean density is -41.6 mK. Each stack isdetected at over σ with stacked images, indicating thatsuch observations could distinguish between these two cases.The depth of the stacked profiles is linearly proportionalto the spin temperature. As a result, distinguishing betweenIGM in emission or absorption will be possible with thesemeasurements. However, a more precise measurement of thespin temperature will be more difficult using this method, saturated unsaturated − A m p li t ud e ( m K ) saturated unsaturated10 N galaxies w i d t h ( c M p c h − ) − − m K Figure 11.
Top: Stacked brightness temperature slices in a sat-urated (left) and unsaturated (right) IGM at z = , where thecentral feature appears as a trough or peak in brightness temper-ature respectively. Bottom: Gaussian fit parameters for stacks ofup to 30 images, comparing the saturated and unsaturated cases.Distinguishing between these cases can be done with very fewgalaxies. The width of the stacked HII profiles are distributedalmost identically in each case, however the profile depths haveopposite sign. due to degeneracies with the distribution of sizes and shapesof early HII regions. For example, a decrease in fit amplitudemay be due to smaller HII regions, more asymmetrical re-gions, an increase in instrumental noise or foreground levels,or a lower spin temperature. Radio interferometers measure the angular component of the21cm signal in Fourier space, and it is therefore most naturalto present observations of a 2D power spectrum, represent-ing fluctuations on the sky, and along the line-of-sight, onseparate axes. Since our stacked profiles are approximatelyGaussian in real-space, so they are also approximately Gaus-sian in the 2D power spectrum, centred on zero.The EoR must be observed through radio foregroundswhich outshine the EoR signal by several orders of magni-tude (Liu et al. 2014). Fortunately foregrounds tend to besmooth in frequency, and occupy the low k | | region of fourierspace, although the inherent chromaticity of the telescopespush the foregrounds into higher k | | regions, creating a fore-ground ”wedge”. One strategy is to avoid the foreground MNRAS000
Top: Stacked brightness temperature slices in a sat-urated (left) and unsaturated (right) IGM at z = , where thecentral feature appears as a trough or peak in brightness temper-ature respectively. Bottom: Gaussian fit parameters for stacks ofup to 30 images, comparing the saturated and unsaturated cases.Distinguishing between these cases can be done with very fewgalaxies. The width of the stacked HII profiles are distributedalmost identically in each case, however the profile depths haveopposite sign. due to degeneracies with the distribution of sizes and shapesof early HII regions. For example, a decrease in fit amplitudemay be due to smaller HII regions, more asymmetrical re-gions, an increase in instrumental noise or foreground levels,or a lower spin temperature. Radio interferometers measure the angular component of the21cm signal in Fourier space, and it is therefore most naturalto present observations of a 2D power spectrum, represent-ing fluctuations on the sky, and along the line-of-sight, onseparate axes. Since our stacked profiles are approximatelyGaussian in real-space, so they are also approximately Gaus-sian in the 2D power spectrum, centred on zero.The EoR must be observed through radio foregroundswhich outshine the EoR signal by several orders of magni-tude (Liu et al. 2014). Fortunately foregrounds tend to besmooth in frequency, and occupy the low k | | region of fourierspace, although the inherent chromaticity of the telescopespush the foregrounds into higher k | | regions, creating a fore-ground ”wedge”. One strategy is to avoid the foreground MNRAS000 , 1–11 (2019) tacking 21cm Images Around High Redshift Galaxies wedge when making 21cm observations (Dillon et al. 2014).However removal of the foreground wedge will result in theloss of any signal inside of it.To simulate foreground contamination in the wedge, weremove portions of Fourier space in our image below a linegiven by k (cid:107) ≤ mk ⊥ (4)where m = DH E ( z ) sin ( θ ) c ( + z ) . (5)Here D is the comoving distance to the observation, H isthe Hubble constant, E ( z ) ≡ (cid:112) Ω m ( + z ) + Ω Λ , and θ is thebeam angle. We take the horizon θ = π / as a conservativelimit. This is the maximum slope of the wedge for spectrallysmooth foregrounds (Liu et al. 2014), which at redshifts z = { , , , } is m = { . , . , . , . } respectively. Fouriermode bins that lie on this limit are down-weighted by afactor equal to the proportion of their volume above thisline. We also consider models where the slope of the line ismultiplied by a factor < q < , to represent scenarios whereforeground subtraction or other techniques reduce the areaof Fourier space dominated by foregrounds (for examplesof foreground mitigation techniques, see Chapman & Jeli´c2019). We show stacked profiles of 30 images at redshift 9, for q = and q = . in Figure 12. While a clear central trough inbrightness temperature exists in the more optimistic case, itwill be more difficult to detect ionised regions with maximalforeground contamination, where the signal to noise ratio isreduced to approximately 1.Since most of the signal power for an ionised region isat large scales in both perpendicular and parallel modes,the amount of signal rendered unusable by foreground con-tamination is roughly proportional to the slope of the wedge.Comparing Gaussian fits for foreground models with ≤ q ≤ in Figure 13, we can see that a significant amount of fore-ground mitigation will be required to detect early ionisedregions with a stack of 30 redshifted 21cm images aroundbright galaxies.Figure 14 shows the signal to noise ratio (Equation 3)obtained when stacking larger numbers of galaxies, with dif-fering levels of foreground contamination at redshifts 10, 9and 8. Increasing the number of galaxies improves the SNRfor low numbers of galaxies. However, as in the case of per-fect foreground subtraction, beyond ∼ galaxies there arediminishing returns, because you will be stacking dimmergalaxies which add less to the signal.Tentative detections SNR > could be achieved forearly ionised regions at redshift z = with 100 galaxiesand full foreground contamination. However, for a detec-tion with SNR greater than 2, partial foreground subtractionwill be required. An SNR of ( , , , ) will require a frac-tion of foreground contaminated Fourier space of at most q = ( . , . , . , . ) for less than 200 galaxies at redshift z = . A similar parameterisation of foreground contamination wasused in Hassan et al. (2019) k ⊥ (cMpc − h) . . . . k || ( c M p c − h ) q = 0.5 k ⊥ (cMpc − h)q = 1.0 x (cMpch − ) y ( c M p c h − ) x (cMpch − ) − . − . − . − . . . . . . m K − m K Figure 12.
Two foreground models at redshift 9, showing thesignal in cylindrically averaged Fourier space (top), and an x-zslice through the centre of the stacked profile (bottom). We showresults for maximum foreground contamination, q = , and thecase where foreground subtraction methods decrease the horizonslope by 50 % ( q = . ) . − − − A m p li t ud e ( m K ) w i d t h ( c M p c h − ) . . . . . . foreground gradient coefficient (q) S N R Figure 13.
The effect of radio foregrounds contaminating differ-ent regions of the 2D power spectrum. Each panel plots a Gaus-sian fit parameter versus the slope coefficient of the foregroundwedge, where q = , represents the case of no contamination, and q = represents the case of maximal foreground contamination forspectrally smooth foregrounds. Errorbars represent the standarddeviation of the fit parameters. Top: Amplitude of the Gaussianfit. Center: σ width of the Gaussian fit, Bottom: Average signalto noise for stacks at each foreground coefficientMNRAS , 1–11 (2019) James E. Davies . . . . . q ga l a x i e s z = 10 . . . . . qz = 9 . . . . . qz = 8 SNR
Figure 14.
Signal to noise ratio (SNR) versus foreground level and number of galaxies stacked at redshifts 10,9 and 8. Contours fromright to left show
SNR = ( , , , , ) We have used the
Bluetides hydrodynamic simulation tocreate mock observations of 21cm brightness temperaturestacked around the brightest galaxies in the simulation. Wehave shown that stacking images of 21cm brightness tem-perature centred around bright galaxies provides a promis-ing way of detecting the first large cosmological HII regionsduring reionization. Our main findings are as follows: • The detection of HII regions can be achieved with im-age stacks around relatively low numbers of bright galaxies.For example, in an idealised case with no foreground con-tamination, a σ detection could result from 30 galaxies atredshift 9. • Stacked profiles will provide an estimate of the mean ofthe bubble size distribution, which could be recovered fromthe width of stacked profiles. • By binning sources according to their brightness, an es-timate of the relationship between UV magnitude and bub-ble size could be obtained. • The sign of the brightness temperature of the IGM af-fects the sign of the stacked profiles, which could be used todistinguish between IGM that is in net absorption or emis-sion of 21cm photons. • With enough stacked images, HII regions could be de-tected in the presence of foreground contamination. For ex-ample, if 100 21cm images are stacked around the brightestavailable galaxies in an SKA1-low field, a tentative . σ de-tection could be made with full foreground contaminationat redshift 9. However, if foreground subtraction reduces thecontaminated region of Fourier space by 50 (80) percent, a 3(5) σ detection of early ionised regions can still be achievedwith the same number of galaxies. ACKNOWLEDGEMENTS
JD is supported by the University of Melbourne under theMelbourne Research Scholarship (MRS). This research was also supported by the Australian Research Council Cen-tre of Excellence for All Sky Astrophysics in 3 Dimensions(ASTRO 3D), through project number CE170100013. Thiswork was performed on the OzSTAR national facility atSwinburne University of Technology. OzSTAR is funded bySwinburne University of Technology and the National Col-laborative Research Infrastructure Strategy (NCRIS). TheBlueTides simulation was run on the BlueWaters facility atthe National Center for Supercomputing Applications. TDMand RACC are grateful for the hospitality of the Physics De-partment of the University of Melbourne, and for the sup-port of a Shimmins Fellowship (TDM) and Lyle Fellowships(TDM and RACC). TDM and RACC acknowledge supportfrom NASA NNX17AK56G, NASA ATP 80NSSC18K101and NSF AST-1614853. RACC was also supported by NSFAST-1615940 and NSF AST-1909193.
REFERENCES
Battaglia N., Trac H., Cen R., Loeb A., 2013, The AstrophysicalJournal, 776, 81Blyth S., et al., 2015, in Advancing Astrophysics with the SquareKilometre Array (AASKA14). p. 128 ( arXiv:1501.01295 )Chapman E., Jeli´c V., 2019, arXiv e-prints, p. arXiv:1909.12369Datta K. K., Bharadwaj S., Choudhury T. R., 2007, MNRAS,382, 809Datta K. K., Majumdar S., Bharadwaj S., Choudhury T. R., 2008,MNRAS, 391, 1900Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604Dillon J. S., et al., 2014, Phys. Rev. D, 89, 023002Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., WilkinsS., 2016, MNRAS, 455, 2778Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433,181Geil P. M., Wyithe J. S. B., 2008, MNRAS, 386, 1683Geil P. M., Mutch S. J., Poole G. B., Duffy A. R., Mesinger A.,Wyithe J. S. B., 2017, Monthly Notices of the Royal Astro-nomical Society, 472, 1324Ghara R., Choudhury T. R., 2019, arXiv e-prints, p.arXiv:1909.12317 MNRAS000
Battaglia N., Trac H., Cen R., Loeb A., 2013, The AstrophysicalJournal, 776, 81Blyth S., et al., 2015, in Advancing Astrophysics with the SquareKilometre Array (AASKA14). p. 128 ( arXiv:1501.01295 )Chapman E., Jeli´c V., 2019, arXiv e-prints, p. arXiv:1909.12369Datta K. K., Bharadwaj S., Choudhury T. R., 2007, MNRAS,382, 809Datta K. K., Majumdar S., Bharadwaj S., Choudhury T. R., 2008,MNRAS, 391, 1900Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604Dillon J. S., et al., 2014, Phys. Rev. D, 89, 023002Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., WilkinsS., 2016, MNRAS, 455, 2778Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433,181Geil P. M., Wyithe J. S. B., 2008, MNRAS, 386, 1683Geil P. M., Mutch S. J., Poole G. B., Duffy A. R., Mesinger A.,Wyithe J. S. B., 2017, Monthly Notices of the Royal Astro-nomical Society, 472, 1324Ghara R., Choudhury T. R., 2019, arXiv e-prints, p.arXiv:1909.12317 MNRAS000 , 1–11 (2019) tacking 21cm Images Around High Redshift Galaxies Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2017,MNRAS, 464, 2234Greig B., Mesinger A., 2017, MNRAS, 465, 4838Hassan S., Liu A., Kohn S., La Plante P., 2019, MNRAS, 483,2524Hinshaw G., et al., 2013, ApJS, 208, 19Hopkins P. F., 2013, MNRAS, 428, 2840Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19Kohler K., Gnedin N. Y., Miralda-Escud´e J., Shaver P. A., 2005,ApJ, 633, 552Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023018Ma Q.-B., Ciardi B., Kakiichi K., Zaroubi S., Zhi Q.-J., Busch P.,2020, ApJ, 888, 112Majumdar S., Bharadwaj S., Choudhury T. R., 2012, MNRAS,426, 3178Marshall M. A., Ni Y., Di Matteo T., Wyithe J. S. B., WilkinsS., Croft R. A. C., 2019, arXiv e-prints, p. arXiv:1912.03428Mellema G., et al., 2013, Experimental Astronomy, 36, 235Mesinger A., Greig B., Sobacchi E., 2016, MNRAS, 459, 2342Ni Y., Di Matteo T., Gilli R., Croft R. A. C., Feng Y., NormanC., 2019, arXiv e-prints, p. arXiv:1912.03780Okamoto T., Frenk C. S., Jenkins A., Theuns T., 2010, MNRAS,406, 208Pober J. C., et al., 2013, The Astronomical Journal, 145, 65Pober J. C., et al., 2014, The Astrophysical Journal, 782, 66Read J. I., Hayfield T., Agertz O., 2010, MNRAS, 405, 1513Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757Springel V., Hernquist L., 2003, MNRAS, 339, 289Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V.,Hernquist L., 2013, MNRAS, 436, 3031Vogelsberger M., et al., 2014a, MNRAS, 444, 1518Vogelsberger M., et al., 2014b, MNRAS, 444, 1518Waters D., Wilkins S. M., Di Matteo T., Feng Y., Croft R., NagaiD., 2016, MNRAS, 461, L51Wilkins S. M., Feng Y., Di Matteo T., Croft R., Lovell C. C.,Waters D., 2017, MNRAS, 469, 2517Wyithe J. S. B., Loeb A., Barnes D. G., 2005, ApJ, 634, 715
APPENDIX A: FIT SYSTEMATICS
In this appendix we first test the systematics of our fittingprocess, and to justify our usage of the Gaussian fit ampli-tude as our signal-to-noise metric. We apply the same sam-pling and fitting method to random positions in the
Blue-tides volume, to examine the case of non-detection. FirgureA1 shows that random pointings are fit very closely to zeroamplitude, and shows a very high variance in width, sincedensity fluctuations and noise can appear on different scalesnear the center of the image. Since these widths include thoseof the real HII regions, we only use the height of the profilesto measure signal-to-noise.We have also tested the convergence of the model fordifferent numbers of bootstrap samples of 30 galaxies. FigureA2 shows that as more galaxy samples are taken, variationin both the parameters and signal to noise ratio decreases.However, both the parameters and SNR vary by less than10% for 100 - 1000 galaxy samples.The size of the sub-cubes that we stack has implica-tions for the Fourier-space resolution when we add thermalnoise and foregrounds. To test how sub-cube size affects thegaussian parameters, we compare stacks of sub-cubes of sidelength 40, 80, 120 and 160 h − cMpc . We show in Figure A3that above 80 h − cMpc , differences in parameters and SNRremain under 10 percent. Galaxies Random − − − − A m p li t ud e ( m K ) Galaxies Random10 N galaxies . . . . w i d t h ( c M p c h − ) Figure A1.
Fits around stacks of 30 bright galaxies and stacksof 30 random pointings compared. Top panels: × × − cMpc slices of stacked galaxy (left) and random (right) pointings. Mid-dle: Gaussian fit amplitude for both galaxy and random pointings.Bottom: Gaussian fit widths for both cases. Random pointings arefit closely to zero amplitude, and show a wide range of widths,consistent with noise.This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS , 1–11 (2019) James E. Davies × − − − A m p li t ud e ( m K ) gal rand0 . . . . . . × w i d t h ( c M p c h − ) galrand Number of Samples S N R Figure A2.
Fit parameters for different numbers of galaxy sam-ples. Top: Gaussian fit amplitudes. Middle: Gaussian fit widths.Bottom: signal to noise ratios, where the shaded area representsthe 95th percentile limits of SNR when the process is repeatedmany times. As more galaxy samples are taken, variation in pa-rameters and SNR decreases. However, variation in SNR remainsunder 10% for 50 - 1000 galaxy samples
40 60 80 100 120 140 160 − − − − A m p li t ud e ( m K )
40 60 80 100 120 140 1600246810 w i d t h ( c M p c h − )
40 60 80 100 120 140 160
Sub-cube side length (cMpc h − ) S N R Figure A3.
Fit parameters from different size sub-cubes. Top:Gaussian fit amplitudes. Middle: Gaussian fit widths. Bottom:Signal to noise ratio. Results for sub-cubes larger than cellsshows little change in both fit parameters and SNR MNRAS000