Measuring the properties of reionised bubbles with resolved Lyman alpha spectra
MMNRAS , 1–11 (2020) Preprint 29 April 2020 Compiled using MNRAS L A TEX style file v3.0
Measuring the properties of reionised bubbleswith resolved Lyman alpha spectra
Charlotte A. Mason (cid:63) † and Max Gronke , † Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA, 02138, USA Department of Physics and Astronomy, University of California, Santa Barbara, 93106, USA Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Identifying and characterising reionised bubbles enables us to track both their sizedistribution, which depends on the primary ionising sources, and the relationship be-tween reionisation and galaxy evolution. We demonstrate that spectrally resolved z ∼ > α ) emission can constrain properties of reionised regions. Specifically,the distant from a source to a neutral region sets the minimum observable Ly α ve-locity offset from systemic. Detection of flux on the blue side of the Ly α resonanceimplies the source resides in a large, sufficiently ionised region that photons can es-cape without significant resonant absorption, and thus constrains both the sizes ofand the residual neutral fractions within ionised bubbles. We estimate the extent ofthe region around galaxies which is optically thin to blue Ly α photons, analogousto quasar proximity zones, as a function of the source’s ionising photon output andsurrounding gas density. This optically thin region is typically (cid:46) . (cid:38) −
250 km s − ), ∼ <
10% of the distance to the neutralregion. In a proof-of-concept, we demonstrate the z ≈ . α peak – likely resides in an ionised region > . < − . . To ionise its own proximity zone we infer COLA1 has a steep UVslope ( β < − . f esc > .
48) and moderateline-of-sight gas density ( < . × the cosmic mean). Key words: dark ages, reionisation, first stars – galaxies: high-redshift
Understanding the process of Hydrogen reionisation is oneof the frontiers of astronomy. It occurred neither homoge-neously nor instantaneously, as ionising photons propagat-ing from nascent galaxies reionised the most overdense re-gions first, carving out ionised ‘bubbles’ within the thenneutral Universe, gradually reionising the entire intergalac-tic medium (IGM). Measuring the timeline and topology ofreionisation, i.e., studying the redshift evolution and spatialdistribution of these ionised regions, is key to understand-ing how reionisation occurred (e.g., Furlanetto et al. 2004;McQuinn et al. 2007; Mesinger 2016).A key question is what drove reionisation, that is,where did the ionising photons originate from? Identifyingreionised or neutral regions of the IGM not only charac-terises the topology of reionisation but enables us to ad-dress this question by comparing the properties of observed (cid:63)
E-mail: [email protected] † Hubble Fellow galaxies in those regions to the local ionisation state (e.g.,Beardsley et al. 2015). Regions which reionise early are likelythe first overdensities where galaxy formation is accelerated,thus identifying those regions helps to identify the first gen-erations of galaxies.Mapping reionised bubbles and measuring their size dis-tribution is a goal of future 21 cm intensity experiments. Thisrequires spatial resolution capable of discerning ionised hy-drogen gas on scales of < i regions (e.g., SKA-low, Koopmans et al.2015), and detailed spectroscopic follow-up of the galax-ies within H ii regions to determine their ionising proper-ties. However, estimates of bubble sizes on small scales arecurrently feasible with Lyman-alpha (Ly α , rest wavelength1216 ˚A) spectroscopy of high redshift sources.Due to its high cross-section for absorption by neutralhydrogen, Ly α is a sensitive probe of neutral gas. Neutralhydrogen affects both the strength and lineshape of Ly α (see, e.g., Dijkstra 2014, for a review). With the advent of c (cid:13) a r X i v : . [ a s t r o - ph . GA ] A p r Mason and Gronke sensitive near-IR spectroscopy, Ly α emission from galaxiesand quasars at z > α emission from galaxies at z ∼ > α also encodes information aboutneutral hydrogen structures the photons encountered alongtheir path. Within or in close proximity to the emittinggalaxy the Ly α spectrum is shaped by resonant scatter-ing with H i that is not necessarily along the line-of-sight(Eide et al. 2018), and typically produces double-peakedemission line profiles due to the high optical depth at linecentre (e.g., Neufeld 1990). However, at larger distances,when the probability of scattering back into the line-of-sight becomes negligible, the impact of intervening H i canbe treated more simply as absorption. The smooth damp-ing wing, due to n hi ∼ > − cm gas that can be at largedistances, is commonly interpreted as a signature of reioni-sation ( e.g., Miralda-Escude 1998). In this case, the opticaldepth due to damping wing absorption is a function of thedistance to the nearest neutral patch and thus could be usedto recover the size of ionised bubbles (Malhotra & Rhoads2006).Due to the short recombination time at z ∼ > ∼ theage of the universe) even within ‘ionised’ bubbles there canbe significant residual neutral gas. The amount of neutralhydrogen depends on the local ionisation field and can leadto resonant absorption on the blue side of the Ly α resonance(e.g., Gunn & Peterson 1965; Zheng et al. 2010; Laursenet al. 2011) – which makes the detection of blue Ly α peakstowards higher redshift increasingly unlikely. In a numberof rare sightlines, however, blue Ly α flux has been observedat z ∼ > z > α emission line-shapes encode information about the ionised bubbles theirhost galaxies reside in. While previous works have shownthat Ly α can be visible early in the epoch of reionisationif galaxies reside in ionised bubbles (Haiman 2002; Masonet al. 2018b) we show here that spectroscopic measurementsof Ly α emitters enable us to calculate the minimum size ofthe ionised bubble such that Ly α at a given frequency off-set is visible to us. Furthermore, we demonstrate that blue-peaked Ly α lines observed at z ∼ > α emittedCOLA1 (Hu et al. 2016; Matthee et al. 2018).This paper is organised as follows: we describe ourmodel for the Ly α optical depth in Section 2 and present ourresults in Section 3. We discuss our results in Section 4 andpresent conclusions in Section 5. We use the Planck Collab-oration et al. (2015) cosmology: (Ω Λ , Ω m , Ω b , n, σ , H ) =(0 . , . , . , . , . ,
68 km s − Mpc − ). Magnitudes are in the AB system. Distances, volumes, and densities areproper unless otherwise stated. In this section, we describe the two components of ourmodel: the Ly α optical depth as a function of the distancefrom a galaxy ( § § The Ly α optical depth through hydrogen gas for photonobserved at λ obs = λ em (1 + z s ) to a source at redshift z s ,observed at z obs , is given by: τ α ( λ obs ) = (cid:90) z s z obs d z c d t d z x hi ( z ) n h ( z ) σ α (cid:18) λ obs z , T (cid:19) (1)where n h is the total number density of hydrogen and x hi is the fraction of hydrogen which is neutral. σ α ( λ, T ) is theLy α scattering cross-section through an ensemble of hydro-gen atoms with a Maxwell-Boltzmann velocity distribution,usually expressed as function of the dimensionless frequency x = ( ν − ν α ) / ∆ ν d : σ α ( x, T ) = σ × φ ( x ) , (2) σ = 1∆ ν d √ π f α πe m e c ≈ . × − (cid:18) T K (cid:19) − / cm where f α = 0 .
416 is the Ly α oscillator strength, m e and e are the mass and charge of an electron, and φ ( x ) is the Voigtfunction: φ ( x ) = a v π (cid:90) ∞−∞ d y e − y ( y − x ) + a v . (3)Here, ν α ≈ . × Hz is the resonant frequency of Ly α ,at wavelength λ α ≈ ν d = ν α (cid:112) k b T /m p c ≡ ν α v th /c is the thermally broadened frequency, and the Voigtparameter a v ≈ . × − ( T / K) − / . Equation (3) isnormalised such that (cid:82) d x φ ( x ) = 1. The cross-section istightly peaked around the core of the line, but has dampingwings which extend out to > − from the line cen-tre (e.g., Dijkstra 2014). We use the approximation for φ ( x )given by Tasitsiomi (2006).Approximating Equation (3) as a Dirac delta function,and assuming constant x hi , we obtain the Gunn & Peterson(1965) optical depth for blue photons emitted from a sourceat z s : τ gp ( z s ) = f α πe m e ν α x hi n h ( z s ) H ( z s ) ≈ . × x hi (cid:18) n h n h (cid:19) (cid:18) z s (cid:19) / . (4)We assume that the source galaxy resides inside an ionisedregion embedded in a neutral homogeneous intergalacticmedium at a distance R ion . This is representative of reioni-sation’s early pre-overlap phases when ionised bubbles growaround sources of ionising photons. However, the assump-tion of an isolated bubble breaks down as reionisation pro-gresses, meaning our method provides only a lower limit onthe bubble size.We construct the optical depth to a source galaxy by MNRAS , 1–11 (2020) roperties of reionised bubbles Rest wavelength [Å] L y t r an s m i ss i on , e x HI ∝ r x HI const. x HI ∝ r x HI const. Velocity offset, v [km/s] x HI (0.1 pMpc) = 10 R i on [ p M p c ] Rest wavelength [Å] L y t r an s m i ss i on , e Velocity offset, v [km/s] R ion = 1 pMpc 9876543210 l og x H I ( r = . M p c ) Figure 1. Ly α transmission for a source at z s = 7 using a x hi ∝ r − profile (solid lines). Left:
Varying the distance to first neutral patchwhile fixing x hi ( r = 0 . − ( n hi ∼ − cm − ). For decreasing bubble size, the transmission on the red side of the Ly α linecentre decreases, due to the increasing damping wing absorption. However, even in a fully neutral IGM, Ly α can be visible providing itis emitted at >
200 km s − . Thin dashed lines show the transmission if x hi = 10 − is constant inside the ionised region. The small blackarrow shows R ion = 1 pMpc to compare with the right panel. Right:
Changing the residual neutral fraction inside the ionised regionwhile holding the distance to the fully neutral patch, R ion = 1 pMpc, fixed. For increasing residual neutral fraction, the transmission onthe blue side decreases as the gas becomes optically thick to Ly α photons redshifting to the resonant frequency by x hi ∼ > − . For higherneutral fractions, the damping wing absorption due to residual neutral gas inside the ionised region also become significant, reducingtransmission on the red side of the line. modelling the two media separately, i.e. breaking the inte-gral into two components: from z s to z ion and z ion to z obs (following, Haiman 2002; Cen & Haiman 2000; Mesingeret al. 2004). In the ionised bubble we set T = 10 K (ap-propriate for photoionised gas at the mean density, e.g.,Hui & Gnedin 1997) and n h ( z ) = C hii n h ( z ). C hii is theclumping factor of ionised hydrogen, C hii ≡ (cid:104) n hii (cid:105) / ¯ n hii , afree parameter which enhances the rate of recombinations,and n h ( z ) is the cosmic mean hydrogen number density: n h ( z ) = ≈ . × − (1 + z ) cm − . In the neutral IGM,we set T = 1 K, assuming gas decouples from the CMB at z ∼
150 (Peebles 1993), and n h ( z ) = n h ( z ). Our results arenot strongly sensitive to these temperatures choices withina physically motivated range ( T < K).The left panel of Figure 1 shows the Ly α transmission, e − τ ( λ ) , as a function of wavelength – commonly expressedas velocity offset, ∆ v ≡ c ( λ em /λ α −
1) – and ionised bub-ble radius, assuming the residual neutral fraction inside theionised bubble is very low ( x hi = 10 − at 0.1 pMpc fromthe source, assuming x hi ∝ r – see § − of the red side can be very low. Ly α linesobserved with low velocity offsets from systemic must residein large ionised regions.Note that even in a fully neutral IGM ( R ion = 0) Ly α flux can still be transmitted on the red side: it is possibleto observe Ly α lines at very high redshifts, providing theyemit Ly α ∼ >
300 km s − from systemic (Dijkstra et al. 2011).Therefore, even at very high redshifts, merely detection ofLy α is not sufficient to identify a reionised bubble: theremust be flux <
300 km s − .The right panel of Figure 1 shows the transmissionthrough a bubble of fixed size (1 pMpc), but changing theresidual neutral fraction in the ionised region around thesource. The damping wing set by R ion acts as an enve-lope for the maximum possible transmission: for x hi ( r =0 . < − , the bubble is fully optically thin andthe maximum blue flux allowed given the damping wingshape can be transmitted. As the residual neutral frac-tion increases, more flux on the blue side of the line is ab-sorbed, with the ionised region becoming optically thick for x hi ( r = 0 . ∼ > − (corresponding to a neutral hydro-gen number density of n hi (cid:38) × − cm − ). For x hi > − the transmission displays a strong damping wing on the redside, though note the transmission in this case is lower thanthe R ion = 0 Mpc case in the left panel: in the right panel weassume gas in the bubble is heated to 10 K, which generatesa higher optical depth than the cold neutral IGM.
MNRAS , 1–11 (2020)
Mason and Gronke
The optical depth can be calculated for any values ofionised bubble size, R ion , and residual neutral fraction insidethe bubble, x hi , to estimate those parameters in a model-independent way. In the limiting case of a single ionisingsource (plus uniform ionising background) we can also esti-mate those quantities for a physical model. Assuming ionisation by a single source at redshift z s at thecentre of the ionised region, the proper radius of the regioncan be obtained by solving for the evolution of a ionisationfront (e.g, Shapiro & Giroux 1987; Cen & Haiman 2000;Yajima et al. 2018):d R d t = 3 f esc ˙ N ion πn h ( z s ) − C hii n h ( z s ) α b ( T ) R + 3 H ( z ) R (5)where the first term is due to ionisations from a source withionising photon output ˙ N ion (in units of s − ) and ionisingescape fraction f esc , the second term is due to recombina-tions – assuming Case B recombination in a clumpy mediumwith clumping factor C hii . For α b , we use the approximationfrom Hui & Gnedin (1997) for the hydrogen recombinationcoefficient as a function of temperature. The third term ofEq. (5) is the expansion of the region due to the Hubbleflow.For constant ˙ N ion and f esc , and simplified cosmology,Equation (5) can be solved analytically (Shapiro & Giroux1987). For instance, for a luminous source at z ∼ < R ion ≈ (cid:18) f esc ˙ N ion t age πn h ( z s ) (cid:19) / (6)where t age is the time since the ionising source has switchedon. In reality, due to the Hubble expansion the ionised radiusgrows more rapidly after ∼ years than Equation (6).Here, we solve Equation (5) numerically.The source emissivity ˙ N ion , in s − can be written as˙ N ion ( t ) = (cid:90) ∞ ν H d ν L ν ( t ) hν (7)where L ν is the ionising spectrum of the source inerg s − Hz − , which can be approximated as a power-lawfor hydrogen ionising photons: L ν ∝ ν − α for ν > ν h , where ν h ≈ . × Hz is the frequency of hydrogen photoion-isation, α is the spectral slope of the ionising continuum,typically 1 ∼ < α ∼ < ∼ N ion for galaxies we assume the spectrum isa double power law: L ν ∝ ν − α for ν > ν H , and L ν ∝ ν − β for ν < ν H . Thus we can estimate ˙ N ion from a UV magnitude(measured at 1500 ˚A) as:˙ N ion ≈ . × (cid:18) z (cid:19) − . M uv +20) α (1 + z ) β +2 (cid:18) (cid:19) β +2 s − Distance from source [pMpc], z s = 7 N eu t r a l h y d r ogen f r a c t i on , x H I optically thickincluding bg no bg RM UV = 16 M UV = 18 M UV = 20 M UV = 22 Figure 2.
Residual neutral hydrogen fraction inside an ionisedregion, with (solid lines) and without (dashed lines) includingionising background flux from Khaire & Srianand (2019), as afunction of the distance from the source and source UV mag-nitude. R ion is determined by solving Equation (5). Calculatedassuming f esc = 1 , C hii = 3 , β = − , α = 2. The dotted verticallines mark the radius of the optically thin region, R α , defined by x hi ∼ < . × − ( § (8)Galaxy spectra typically have a steeper drop-off beyond theHeII ionising limit (54.4 eV), which isn’t captured in oursimple power-law approximation. However, we note that thishas only a small impact on our estimation of ˙ N ion : for 1 <α < ν > ν HeII produces > .
75 ˙ N ion as defined in Equation (8). Due to the recombination of ionised hydrogen, inside theionised region there will be some residual neutral fraction.This can be computed by equating the recombination rateto the ionisation rate, assuming ionisation equilibrium. As-suming ionisations due to the central source and some diffuseionising background, the residual neutral fraction at a properradius r from the source is (e.g., Mesinger et al. 2004): x hi ( r ) = C hii n h ( z ) α b ( T ) (cid:18) Γ bg ( z ) + J s πr (cid:19) − . (9)Here, Γ bg is the hydrogen ionising rate due to the back-ground within the ionised region in s − , and J s is the hy-drogen ionising emissivity of the central source in cm s − .The gas reaches ionisation equilibrium with a characteris-tic timescale t − = Γ bg + J s / πr , so the bubble will be inionisation equilibrium within ∼ years assuming constantemissivity (e.g., Davies et al. 2020). Thus Equation (9) holdsfor sources with ionising populations > years, which is MNRAS , 1–11 (2020) roperties of reionised bubbles Wavelength [Å] L y t r an s m i ss i on , e M UV = 16 M UV = 18 M UV = 20 M UV = 221000 500 0 500 1000 Velocity offset [km/s]
Figure 3. Ly α transmission curves for the models shown in Fig-ure 2, calculating the optical depth as described in § reasonable for massive galaxies at z ∼ −
8, though maybreak down for galaxies with short bursts of star formation.We use the ionising background model by Khaire & Sri-anand (2019) but note that it does not significantly impactthe residual neutral fraction for bubbles at z > x hi ∝ r . Strictly, Γ bg accounts for other ionis-ing sources nearby and will therefore vary depending on thedensity of the environment. We expect Γ bg ∼ < (cid:104) Γ bg (cid:105) basedon fluctuations of density and mean free path (Mesinger& Dijkstra 2008; Davies & Furlanetto 2016). For reference (cid:104) Γ bg ( z ∼ (cid:105) ≈ . × − s − in the Khaire & Srianand(2019) model, and has been measured to be ∼ < . × − s − at z ∼ J s = f esc (cid:90) ∞ ν H d ν L ν hν σ ion ( ν ) (10)with the hydrogen photoionisation rate σ ion = σ ion , ( ν/ν h ) − , where σ ion , ≈ . × − cm (e.g.,Draine 2011). Assuming as above an ionising spectrum L ν ∝ ν − α yields J s = f esc ˙ N ion αα + 3 σ ion , (11)where ˙ N ion is given by Equation (8).Figure 2 shows some typical neutral fraction profilesinside a HII region. Here, and below, we assume that theneutral fraction is unity outside the HII region at the ra-dius determined by Equation 5. We see that more luminousgalaxies produce bubbles which are both larger (Equation 6)and more highly ionised at a fixed distance from the source.Figure 3 shows the Ly α for the same set of models. Only UVbright galaxies are capable of producing a sufficiently largeionised region to allow blue flux to be observed. Importantly, due to the high resonant cross-section of Ly α (Equation (2)) an ionised bubble can still be optically thickto Ly α . Thus blue Ly α flux can be suppressed by residualneutral gas within an ionised bubble. The proper radius atwhich the bubble becomes optically thick to Ly α is the ra-dius where the Gunn & Peterson (1965) optical depth (Equa-tion (4), using x hi ( r ) given by Equation 9) exceeds an opticaldepth threshold τ lim ∼ . ∼ R α = (cid:18) J s π (cid:19) / (cid:20) C hii n h ( z ) α b ( T ) H ( z ) τ lim f α πe m e ν α − Γ bg ( z ) (cid:21) − / (12) ≈ . (cid:18) τ lim ˙ N ion f esc . × . × s − (cid:19) / (cid:18) . αα + 3 (cid:19) / × (cid:18) C hii (cid:19) − (cid:18) T K (cid:19) . (cid:18) z (cid:19) − / Mpcwhere J s is given by Equation (11). For the latter equality weassumed Γ bg = 0, α b ( T ) ≈ . × − ( T / K) − . cm s − and used ˙ N ion for a M uv = −
20 galaxy (Equation 8).This radius corresponds to reaching a neutral hy-drogen number density of n hi ∼ > . × − ( τ lim / . z ) / / cm − – or x hi ∼ > . × − ( C hii / − ( τ lim / . z ) / − / – in the ionised region. At higher densities/neutralfractions the gas is optically thick to Ly α photons.This is analogous to the quasar near/proximity zonesdescribed by Bolton & Haehnelt (2007), except here weinclude the contribution of other, diffuse sources of ionis-ing photons. As discussed in Section 2.2.2 we assume thereionised region is in ionisation equilibrium, which is validfor for sources with ionising populations > years. SeeDavies et al. (2020) for discussion of the time evolution ofsuch optically thin regions around quasars.A lower limit on R α can thus be estimated from theminimum observable blue Ly α velocity offset ∆ v min α : R α > | ∆ v min α | H ( z s ) (13)Previous works, which assumed ionised bubbles are opticallythin to Ly α (e.g., Matthee et al. 2018; Hashimoto et al. 2018)estimated R ion > | ∆ v min α | /H ( z s ). From the above, we seethis is actually measuring R α and is an underestimate of R ion . We will show below in Section 3.1 that R α (cid:28) R ion . The Ly α optical depth (Equation 1) decreases with decreas-ing redshift, due to the increasing ionising output of sources,and the reducing density of neutral gas due to cosmic ex-pansion. Thus, we expect the optically thin regions aroundgalaxies in reionising bubbles to grow with decreasing red-shift, increasing the observable blue flux.Figure 4 shows the total ionised radius ( R ion ) and theoptically thin radius ( R α ) as a function of source redshift,fixing f esc = 1 , C hii = 3 , α = 2 , β = −
2. We comparethe sizes of the bubbles and proximity zones for sourceswith different UV luminosities ( M uv = − , − , − , − MNRAS000
2. We comparethe sizes of the bubbles and proximity zones for sourceswith different UV luminosities ( M uv = − , − , − , − MNRAS000 , 1–11 (2020)
Mason and Gronke
Source redshift, z s R / R i on t s = 10 yr t s = 10 yr t s = 10 yr t s = 10 yr t s = 10 yr t s = 10 yr 5 6 7 8 9 10 Source redshift, z s R [ p M p c ] COLA1MACS1149-JD1 5 6 7 8 9 10
Source redshift, z s v m i n = RH ( z ) [ k m / s ] M UV = 16 M UV = 18 M UV = 20 M UV = 22 M UV = 16 M UV = 18 M UV = 20 M UV = 22 Figure 4.
Radii of ionised bubbles and Ly α proximity zones around galaxies as a function of source redshift and age. Note theseshow optimistic upper limits on the sizes assuming f esc = 1 , C hii = 3 , α = 2 , β = − Left:
Ratio of the Ly α proximity zone ( R α ,Equation 12) to the radius of ionised sphere ( R ion , Equation 5) for galaxies with M uv = − , − , − , −
22 (black, blue, pink, orangelines respectively). We show the time evolution of R α /R ion for the brightest source. Early on, before the ionisation front can growsignificantly, R α ∼ . − . R ion , however, for sources >
100 Myr old, R α ∼ < . R ion . Center:
The proper size of the optically thin region R α in Mpc as a function of source redshift. All lines are the same as the left panel. Right:
Radius of optically thin region expressed asa Ly α velocity offset, i.e. the maximum Ly α blue peak velocity offset that would be observable from a Ly α emitter inside this ionisedregion. We also show the observed blue-peaked Ly α emitters COLA1 (Hu et al. 2016; Matthee et al. 2018), MACS1149-JD1 (Hashimotoet al. 2018). NEPLA4 (Songaila et al. 2018) is at the same redshift as COLA1 and has a similar blue peak velocity. and age (10 , , yrs). Except for very young sources R α ∼ < . R ion . R α does not change in size with age for con-stant emissivity once ionisation equilibrium is reached. Asnoted above, previous works, which assumed ionised bubblesare fully ionised when estimating blue peak transmission,underestimated the total extent of the ionised region whenusing the observed blue Ly α peak. By including recombina-tions we see the blue Ly α flux only probes the much smalleroptically thin region. In the next section, we show that thereare model-independent ways to estimate a lower bound onthe size of the full ionised region. α lines in a mostly neutralmedium Figure 5 shows the minimum Ly α velocity offset observableas a function of the distance to the first neutral region ( R ion )and residual neutral fraction in the bubble ( x hi ). We com-pute Ly α transmission e τ (∆ v ) on a grid of R ion and x hi , en-abling a model-independent estimate of the properties of anionised bubble based on Ly α transmission. We show the ve-locity offsets observable if >
10% of the flux emitted at thatvelocity offset is transmitted through the IGM. We choose10% as assuming an emitted Ly α EW of 200 ˚A, this trans-mitted flux should be observable with current facilities.For small bubbles with high x hi , it is only possible toobserve Ly α which is significantly redshifted. Conversely, itis only possible to observe blue Ly α flux if there is a signif-icant distance to the first neutral patch ( > . x hi < − ).The two panels in Figure 5 compare the minimum ob-servable Ly α velocity offsets in the case of a homogeneousresidual neutral fraction in the bubble ( x hi = constant), ap-proximating reionisation by a uniform ionising backgroundof ultra-faint sources) and in the limiting case of the Ly α emitter as the sole reionising source ( x hi ∝ r , Equation 9), with the value on the y -axis x hi ( r = 0 . α emitters COLA1 is a z ≈ . α peak with fluxup to −
250 km s − from systemic (Hu et al. 2016; Mattheeet al. 2018). From Figure 5 we see that this requires it to re-side in an ionised region at least 0 . x hi < − . Ourestimate of the extent of the ionised region is roughly doublethan that of Matthee et al. (2018), who estimated 0.3 pMpc(2.3 cMpc). This is due to their assumption that the entirebubble is optically thin. As discussed in Section 3.1, withjust the minimum observed blue flux velocity we can onlycalculate R α (Equation 13), but the total ionised region ismuch larger.Figure 6 shows the radius of the proximity zone asa function of source magnitude. We compare R α to theblue peak velocity offsets from Matthee et al. (2018) andHashimoto et al. (2018), and vary f esc = { , . } , C hii = { , } , α = { , } and β = {− , − . } . For COLA1 toionise its own proximity zone likely requires a high escapefraction, steep UV spectral slope β , and relatively low gasdensity, whilst the ionising spectral slope makes a negligibleimpact on the proximity zone size.To investigate in more detail the necessary conditionsfor COLA1’s blue peak to be observable we perform aBayesian inference to infer the parameters in Equation (12).We define the likelihood to observe blue Ly α flux at − ±
60 km s − (the minimum observed blue flux, using the in-strumental resolution R ∼ M uv = − . ± . z = 6 .
6. To observe a
MNRAS , 1–11 (2020) roperties of reionised bubbles Distance to neutral patch, R ion [pMpc] l og x H I - - - - - - - - z s = 7, T > 10%, constant x HI M i n i m u m ob s e r v ab l e v [ k m / s ] Distance to neutral patch, R ion [pMpc] l og x H I ( r = . M p c ) - - - - - - - - z s = 7, T > 10%, x HI ∝ r M i n i m u m ob s e r v ab l e v [ k m / s ] Figure 5.
Minimum observable Ly α velocity offset, ∆ v α as a function of bubble size and residual neutral fraction, assuming > Left for a constant residual neutral fraction, x hi , inside the ionised bubble. Right x hi ∝ r , withthe quoted value at 0.1 Mpc.
24 22 20 18
UV magnitude v m i n = RH ( z = ) [ k m / s ] f esc = 1, C HII = 3= 2, = 2 f esc = 0.2 C HII = 10= 1= 2.5COLA1MACS1149-JD1
Figure 6.
Radius of optically thin region produced by a singlesource of a given UV magnitude. We show the minimum observ-able blue-shifted velocity offset from systemic. The black solid lineshows a fiducial galaxy model ( f esc = 1 , C hii = 3 , α = 2 , β = − f esc = 0 . C hii = 10(orange dots), α = 2 (green dot-dash), β = − . blue peak at ∆ v , the proximity zone R α > | ∆ v | /H ( z ). Thelikelihood is: p ( R α ≥ | ∆ v | /H ( z ) | θ ) = 12 erfc (cid:18) | ∆ v | /H ( z ) − R α, mod ( θ ) √ σ R (cid:19) (14)where we have assumed the probability of a proximity zone having radius R α , p ( R α | θ ), is a normal distributionwith mean R α, mod ( θ ) (Equation 12) and variance σ R =[ σ ∆ v /H ( z )] . We use uniform priors on the parameters θ = [ f esc , C hii , α, β, Γ bg ]: 0 ≤ f esc ≤ , . < C hii < , ≤ α ≤ . , − < β < −
1, and 0 < Γ bg / − s − <
10. Werun two versions of the model: one where we fix the ionis-ing background Γ bg = 0 (i.e. assuming COLA1 ionises itsproximity zone alone), and one where Γ bg is a free param-eter. To estimate the posteriors and the evidence for eachmodel, Z = (cid:82) d θ p ( θ | R α ), we use Dynamic Nested Samplingimplemented in dynesty (Speagle 2019).Figure 7 shows the posteriors for these parameters andtheir median and 16 −
84% credible intervals or 68% up-per/lower limits, and the evidence Z . In both cases, gas den-sity is inferred to be relatively low ( ∼ < . n h , σ ), the UVslope β is inferred to be steep ( β < − .
27 1 σ for Γ bg = 0)and the spectral slope of the ionising continuum, α , is notwell constrained by the proximity zone, due to the smallerrange of possible α having a minimal impact on the sizeof the proximity zone (see Figure 6). In the model withoutan ionising background, high escape fractions are inferred( > . , σ ), while when we include an ionising backgroundthere is a degeneracy between high single source f esc andlow ionising background, or low f esc and high ionising back-ground. Neither model is strongly preferred, with a Bayesfactor Z bg /Z nobg = 1 . z > blue peaks Hashimoto et al. (2018) reported a 4 σ detection of a Ly α linein the z = 9 .
11 source MACS1149-JD1. The Ly α line is offsetby − ±
60 km s − from their detection of [OIII]88 µ m.Based on Figure 5, if the Ly α comes from the samesource as the [OIII], the environment of MACS1149-JD1must be extremely highly ionised ( > − ) and in a bubble ∼ > MNRAS000
60 km s − from their detection of [OIII]88 µ m.Based on Figure 5, if the Ly α comes from the samesource as the [OIII], the environment of MACS1149-JD1must be extremely highly ionised ( > − ) and in a bubble ∼ > MNRAS000 , 1–11 (2020)
Mason and Gronke f esc >0.48 f esc >0.43 No background (ln Z = 1.60)With ionizing background (ln Z = 1.39) C H II C HII <2.60 C HII <2.44 . . . . . =1.72 +0.520.49 =1.72 +0.530.49 . . . . . < 2.27< 2.10 . . . . . f esc bg [ s ] C HII . . . . . . . . . . bg [10 s ] bg [10 s ]>3.5 Figure 7.
Posterior distributions for f esc , C hii , α, β and Γ bg in-ferred from the observed maximum blue Ly α peak of COLA1( −
250 km s − ). We show 1 σ and 2 σ contours of the 2D posteriors,and histograms of marginalised 1D posteriors for the parameters.Blue lines show the model with Γ bg = 0, grey lines the model withΓ bg as a free parameter. The likelihood and priors are describedin Section 3.3.1. produced by galaxies of a given M uv . Given the observedfaintness of MACS1149-JD1 ( M uv = 18 . ± . f esc , steep α andlow gas density. We thus agree with a possible interpreta-tion by Hashimoto et al. (2018) that the Ly α emission comesfrom a different, slightly lower redshift, source compared tothe [OIII] emission in MACS1149-JD1.NEPLA4 (Songaila et al. 2018) is a narrow-band se-lected Ly α emitter at z = 6 . ∼−
250 km s − , similar to COLA1. As the UV continuum isnot known we cannot place it on Figure 6. α constraints on bubble properties Blue Ly α peaks can reveal conditions inside individualreionised bubbles, however we expect blue peaks to be rareat z >
6, due to the high IGM opacity (see § α lines which are red-shifted with respect to systemic may arise more often athigh redshift due to outflows, which aids the transmissionof photons through the IGM (Dijkstra et al. 2011). Figure 5demonstrates that the velocity offset of an observed red peakcan place lower limits on the size of an ionised region (seealso Malhotra & Rhoads 2006).Not only can a single source place constraints on bub-bles sizes, with a deep spectroscopic survey in a single field it could be possible to map an ionised bubble directly. Due tothe the ionisation gradient across bubbles the transmissionof Ly α will vary radially across a bubble (see Figure 8). Aswe demonstrate in Figure 5, the observable minimum veloc-ity offset from systemic varies as a function of the distancefrom the nearest neutral region. Likewise, the transmittedLy α flux will decrease for sources further from the centerof bubbles. If the faint-end of Ly α luminosity function issteep (e.g., α ∼ − . z ∼
7, using theluminosity function model by Gronke et al. (2015) (setting α = − .
8) we expect ∼
12 Ly α emitting galaxies with aluminosity L (cid:38) erg s − (flux ∼ > × − erg s − cm − )located within R α ∼ . α emitters live in over-dense regions (Ouchi et al. 2017) we expect higher numbercounts. Within this optically thin region we would expectthe fraction of galaxies with blue Ly α peaks to be compara-ble to those seen at lower redshifts, when the IGM is highlyionised.With a large near infra-red spectroscopic survey to mea-sure Ly α flux and high S/N resolved lineshapes, as well assystemic redshift from other emission lines (e.g. rest-frameoptical lines visible with JWST ) it could be possible to di-rectly map ionised bubbles and place constraints on the bub-ble size distribution during reionisation. These deep mea-surements will be feasible with 30 m telescope spectroscopy(e.g., E-ELT/MOSAIC is expected to reach 1 × − ergs − cm − in 40 hrs, Evans et al. 2015).To fully interpret such data requires a more realistictreatment of the neutral gas distribution than the uniformmodel used here. For example, Gronke et al. (in prep) ex-plores the prevalence of blue Ly α peaks in the cosmologicalradiative hydrodynamical simulation CoDaII (Ocvirk et al.2018), and finds large line-of-sight variation in Ly α trans-mission due to inhomogeneous gas distributions.While the visibility of blue peaks (or red peaks withsmall velocity offset) can put lower limits on R α , measuringthe size more accurately is difficult. However, given suffi-cient spectral quality, it might be possible to detect a sharp cutoff on the blue side of a blue peak – a signature of absorp-tion, since frequency redistribution yields smoother profilestowards the wings (Neufeld 1990). A sharp cutoff would beexpected if there is a sharp transition from an optically thinionised region to one that is optically thick, for example,from galaxies sustaining their own R α surrounded by ho-mogeneous 10 − ∼ < x hi ∼ < − reionised gas (which may betypical of the IGM at the end of reionisation e.g., Fan et al.2006). This would enable a direct measurement of R α andthus tighter constraints on f esc . In contrast, a sharp cutoffof the red peak towards line centre can be due to eitherthe IGM or radiative transfer effects and is commonly foundalso at low redshifts (e.g., Rivera-Thorsen et al. 2015; Yanget al. 2017), thus, this is harder to use as a measurementfor R α . However, one could use this signature in a statisticalsense: intragalactic radiative transfer sets a characteristiccorrelation between the width and velocity offset of red Ly α peaks (e.g., Neufeld 1990; Verhamme et al. 2018), and thistrend should be altered at high- z due to the IGM absorptionyielding a flatter slope in the width-offset relation. MNRAS , 1–11 (2020) roperties of reionised bubbles R α R ion A optically thick HIIoptically thin HII neutral IGM B C DA B C D
ISM emission IGM transmission observed line
Figure 8.
Illustration of IGM attenuation for galaxies at different radial positions from the center of an ionised bubble. We show relative(not to scale) positions and observed line profiles of galaxies: at the center of the bubble, inside the optically thin region (A), close to theedge of the optically thin region (B), and within the optically thick region (C), in the neutral IGM (D). The plots show the observed Ly α emission lines expected if all the galaxies had the same double-peaked emission line emerging from the ISM. Galaxy A can be observedwith significant blue flux, whereas galaxies C and D have no observed blue flux. A small fraction of Ly α is still visible from galaxy D,but only at ∆ v ∼ >
300 km s − . The increase in the size of the optically thin region meansthat more blue flux will be observed at lower redshifts, evenif the size of the ionised bubble remains the same. The evolv-ing flux distribution of Ly α emission from z > α flux distribution at z > z <
6, Fan et al. 2006), we can ask underwhat conditions could the increased optical depth due to anincreased residual neutral fraction in ionised bubbles causethe observed decline in Ly α emission.Assuming a double-peaked Gaussian Ly α lineprofile,with red:blue flux ratio R : 1, the resulting Ly α transmissionfraction through the IGM can be written as: T ( z = 7) T ( z = 6) = e − τ gp ( z =7) + Re − τ gp ( z =6) + R (15)Figure 9 shows this as a function of the relative increasein the average residual neutral fraction between z ∼ z ∼ x hi ( z = 7) /x hi ( z = 6). In order to produce a drop inthe observed Ly α fraction of T /T ∼ .
5, either the neutralfraction at z = 6 must be < − , i.e. not optically thick,or if x hi ( z = 6) ∼ − the blue peak flux must be ∼ > × stronger than the red peak. Both of these scenarios are un-likely: the z ∼ α forest is optically thick on average (e.g,Fan et al. 2006, find x hi ( z ∼ ∼ > − ) and the observedLy α line shapes of galaxies at all redshifts show dominantred peaks (e.g., Trainor et al. 2015; Rivera-Thorsen et al.2015; Yang et al. 2017; Steidel et al. 2018). Therefore, whilethe residual neutral fraction in ionised regions will increasebetween z ∼ z ∼ α on average, andthus has a small impact on reionisation inferences. We have used an analytic model to estimate the Ly α opticaldepth within reionised bubbles and investigate the impactof reionisation on Ly α lineshapes. Our conclusions are asfollows:(i) Both the size of, and residual neutral hydrogen frac-tion within, reionised bubbles affect the observed lineshapeof Ly α emission during reionisation. As such, measurementsof the Ly α velocity offset from systemic can provide lowerlimits on the source’s distance to the first large neutral patchalong the line of sight, and upper limits on the residual neu-tral fraction inside its HII bubble.(ii) Galaxies with Ly α lines observed with low velocityoffsets from systemic must reside in large reionised regions.Detecting blue Ly α peaks during reionisation requires thesource galaxy live in a R ∼ > . MNRAS000
5, either the neutralfraction at z = 6 must be < − , i.e. not optically thick,or if x hi ( z = 6) ∼ − the blue peak flux must be ∼ > × stronger than the red peak. Both of these scenarios are un-likely: the z ∼ α forest is optically thick on average (e.g,Fan et al. 2006, find x hi ( z ∼ ∼ > − ) and the observedLy α line shapes of galaxies at all redshifts show dominantred peaks (e.g., Trainor et al. 2015; Rivera-Thorsen et al.2015; Yang et al. 2017; Steidel et al. 2018). Therefore, whilethe residual neutral fraction in ionised regions will increasebetween z ∼ z ∼ α on average, andthus has a small impact on reionisation inferences. We have used an analytic model to estimate the Ly α opticaldepth within reionised bubbles and investigate the impactof reionisation on Ly α lineshapes. Our conclusions are asfollows:(i) Both the size of, and residual neutral hydrogen frac-tion within, reionised bubbles affect the observed lineshapeof Ly α emission during reionisation. As such, measurementsof the Ly α velocity offset from systemic can provide lowerlimits on the source’s distance to the first large neutral patchalong the line of sight, and upper limits on the residual neu-tral fraction inside its HII bubble.(ii) Galaxies with Ly α lines observed with low velocityoffsets from systemic must reside in large reionised regions.Detecting blue Ly α peaks during reionisation requires thesource galaxy live in a R ∼ > . MNRAS000 , 1–11 (2020) Mason and Gronke x HI ( z = 7)/ x HI ( z = 6) T / T intrinsicblue:red1 : 110 : 1100 : 1 x HI ( z = 6)10 Figure 9. Ly α transmission ratio z = 7 to z = 6 given by Equa-tion (15), assuming all of the optical depth to Ly α is due toresonant absorption of the blue peak (i.e. no damping wing), asa function of the relative increase in the residual neutral fraction x hi ( z = 7) /x hi ( z = 6). The coloured lines show the ratio for dif-ferent values of x hi ( z = 6) and the linestyles for different intrinsicblue:red peak ratios. For an increase in residual neutral fractionto explain the drop in Ly α transmission T /T ∼ . z ∼ x hi ∼ < − and the blue peak tobe significantly stronger than the red peak, contrary to observa-tions at z ∼ < ( x hi > − ) bubble. By contrast, Ly α can be detected evenfrom fully neutral regions, providing it was emitted at highvelocity offset ( ∼ >
300 km s − ).(iii) Around individual galaxies, we predict the regions oflow Ly α optical depth – proximity zones – to be typically < . α to be detected out to ∼ > −
250 km s − .(iv) The observed blue-peaked Ly α emitter COLA1, withblue flux up to −
250 km s − , must reside in a highly ionisedregion ( x hi < − . ) at least 0.7 pMpc from the nearestlarge scale neutral patch. For COLA1 to have generated itsown proximity zone requires it to have a high escape fraction, f esc > .
48, and steep UV spectral slope, β < − .
27, and forthe total gas density along the line of sight to be relativelylow < . n h . Including an ionising background alleviatesthe need for a high f esc and steep β , but still requires arelatively low gas density.Detailed measurements of Ly α lineshapes and velocityoffsets can be used to constrain the properties of reionisedregions, including to place lower limits on ionised bubblesizes. With rest-UV – optical spectroscopy with, e.g., theJames Webb Space Telescope ( JWST ) these methods en-able direct comparison between galaxy properties and theirionised regions.
ACKNOWLEDGEMENTS
The authors thank Zoltan Haiman, Steve Finkelstein, DanStark and Tommaso Treu for useful discussions. CAMacknowledges support by NASA Headquarters throughthe NASA Hubble Fellowship grant HST-HF2-51413.001-A awarded by the Space Telescope Science Institute, whichis operated by the Association of Universities for Researchin Astronomy, Inc., for NASA, under contract NAS5-26555.MG was supported by NASA through the NASA HubbleFellowship grant HST-HF2-51409 and acknowledges sup-port from HST grants HST-GO-15643.017-A, HST-AR-15039.003-A, and XSEDE grant TG-AST180036.
Software : IPython (P´erez & Granger 2007), matplotlib (Hunter 2007),
NumPy (Van Der Walt et al. 2011),
SciPy (Oliphant 2007),
Astropy (Robitaille et al. 2013), dynesty (Speagle 2019).
REFERENCES
Beardsley A. P., Morales M. F., Lidz A., Malloy M., Sutter P. M.,2015, ApJ, 800, 128Bolton J. S., Haehnelt M. G., 2007, MNRAS, 374, 493Bolton J. S., Haehnelt M. G., 2013, MNRAS, 429, 1695Bosman S. E. I., Kakiichi K., Meyer R. A., Gronke M., LaporteN., Ellis R. S., 2019, p. arXiv:1912.11486Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011,MNRAS, 412, 2543Cen R., Haiman Z., 2000, ApJ, 542, L75Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328Davies F. B., et al., 2018, ApJ, 864, 142Davies F. B., Hennawi J. F., Eilers A.-C., 2020, MNRAS, 493,1330Dijkstra M., 2014, Publ. Astron. Soc. Aust., 31, e040Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS, 414,2139Draine B. T., 2011, Physics of the Interstellar and IntergalacticMediumDrake A. B., et al., 2017, A&A, 608, A6Eide M. B., Gronke M., Dijkstra M., Hayes M., 2018, ApJ, 856,156Evans C., et al., 2015, p. arXiv:1501.04726Fan X., et al., 2006, AJ, 132, 117Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613,16Geil P. M., Mutch S. J., Poole G. B., Duffy A. R., Mesinger A.,Wyithe J. S. B., 2017, MNRAS, 472, 1324G¨otberg Y., de Mink S. E., McQuinn M., Zapartas E., Groh J. H.,Norman C., 2020, A&A, 634, A134Greig B., Mesinger A., Ba˜nados E., 2019, MNRAS, 484, 5094Gronke M., Dijkstra M., Trenti M., Wyithe S., 2015, MNRAS,449, 1284Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633Haiman Z., 2002, ApJ, 576, L1Hashimoto T., et al., 2018, Nature, 557, 392Hoag A., et al., 2018, ApJ, 854, 39Hoag A., et al., 2019, ApJ, 878, 12Hu E. M., Cowie L. L., Songaila A., Barger A. J., RosenwasserB., Wold I. G. B., 2016, ApJ, 825, L7Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27Hunter J. D., 2007, Comput. Sci. Eng., 9, 99Khaire V., Srianand R., 2019, MNRAS, 484, 4174Koopmans L., et al., 2015, in Advancing Astrophysicswith the Square Kilometre Array (AASKA14). p. 1( arXiv:1505.07568 ) MNRAS , 1–11 (2020) roperties of reionised bubbles Laursen P., Sommer-Larsen J., Razoumov A. O., 2011, ApJ, 728,52Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., VignaliC., Stern J., O’Meara J. M., 2015, MNRAS, 449, 4204Malhotra S., Rhoads J. E., 2006, ApJ, 647, L95Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pen-tericci L., de Barros S., Vanzella E., 2018a, ApJ, 856, 2Mason C. A., et al., 2018b, ApJ, 857, L11Mason C. A., et al., 2019, MNRAS, 485, 3947Matthee J., Sobral D., Gronke M., Paulino-Afonso A., StefanonM., R¨ottgering H., 2018, A&A, 619, A136McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldar-riaga M., 2007, MNRAS, 377, 1043Mesinger A., 2016, Understanding the Epoch of CosmicReionization: Challenges and Progress. Vol. 423, Springer,doi:10.1007/978-3-319-21957-8Mesinger A., Dijkstra M., 2008, MNRAS, 390, 1071Mesinger A., Haiman Z., Cen R., 2004, ApJ, 613, 23Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A.,Dijkstra M., 2015, MNRAS, 446, 566Miralda-Escude J., 1998, ApJ, 501, 15Neufeld D. A., 1990, ApJ, 350, 216Ocvirk P., et al., 2018, p. arXiv:1811.11192Oliphant T. E., 2007, Comput. Sci. Eng., 9, 10Ouchi M., et al., 2017, Publ. Astron. Soc. Japan, 00, 1Peebles P. J. E., 1993, Principles of Physical CosmologyP´erez F., Granger B. E., 2007, Comput. Sci. Eng., 9, 21Planck Collaboration et al., 2015, A&A, 594, A8Rivera-Thorsen T. E., et al., 2015, ApJ, 805, 14Robitaille T. P., et al., 2013, A&A, 558, A33Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014,ApJ, 795, 20Scott J. E., Kriss G. A., Brotherton M., Green R. F., HutchingsJ., Shull J. M., Zheng W., 2004, ApJ, 615, 135Shapiro P. R., Giroux M. L., 1987, ApJ, 321, L107Songaila A., Hu E. M., Barger A. J., Cowie L. L., Hasinger G.,Rosenwasser B., Waters C., 2018, ApJ, 859, 91Speagle J. S., 2019, p. arXiv:1904.02180Steidel C. C., et al., 2014, ApJ, 795, 165Steidel C. C., Bogosavljevic M., Shapley A. E., Reddy N. A.,Rudie G. C., Pettini M., Trainor R. F., Strom A. L., 2018,ApJ, 869, 123Stevans M. L., Shull J. M., Danforth C. W., Tilton E. M., 2014,ApJ, 794, 75Tasitsiomi A., 2006, ApJ, 645, 792Trainor R. F., Steidel C. C., Strom A. L., Rudie G. C., 2015, ApJ,809, 89Treu T., Schmidt K. B., Trenti M., Bradley L. D., Stiavelli M.,2013, ApJ, 775, L29Trotta R., 2008, Contemporary Physics, 49, 71Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Comput.Sci. Eng., 13, 22Verhamme A., et al., 2018, MNRAS, 478, L60Whitler L. R., Mason C. A., Ren K., Dijkstra M., Mesinger A.,Pentericci L., Trenti M., Treu T., 2019, p. arXiv:1911.03499Wyithe J. S. B., Hopkins A. M., Kistler M. D., Y¨uksel H., BeacomJ. F., 2010, MNRAS, 401, 2561Yajima H., Sugimura K., Hasegawa K., 2018, MNRAS, 477, 5406Yang H., et al., 2017, ApJ, 844, 171Zheng Z., Cen R., Trac H., Miralda-Escud´e J., 2010, ApJ, 716,574This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000