ApJ,in press PreprinttypesetusingLATEXstyleemulateapjv.6/22/04 CONSTRAINTS ON REIONIZATION AND SOURCE PROPERTIESFROM THE ABSORPTION SPECTRA OF Z >6.2 QUASARS Andrei Mesinger1,2 & Zolta´n Haiman2 ApJ, in press ABSTRACT We make use of hydrodynamical simulations of the intergalactic medium (IGM) to create model quasar absorptionspectra. We compare these model spectra with the observedKeck spectra of three z > 6.2 quasars with full Gunn-Peterson troughs: SDSS J1148+5251 (z = 6.42), SDSS J1030+0524 7 (z = 6.28), and SDSS J1623+3112 (z = 6.22). We fit the probability density distributions (PDFs) 0 of the observed Lyman α optical depths (τα) with those generated from the simulation, by adopting 0 a template for the quasar’s intrinsic spectral shape, and exploring a range of values for the size of 2 the quasar’ssurrounding HII region, R , the volume-weightedmean neutral hydrogenfraction in the S n ambient IGM, x¯HI, and the quasar’sionizing photon emissivity, N˙Q. In order to avoidaveragingover a possibly large sightline–to–sightline fluctuations in IGM properties, we analyze each observed quasar J independently. We find the following results for J1148+5251, J1030+0524, and J1623+3112: The 0 best–fit sizes RS are 40, 41, and 29 (comoving) Mpc, respectively. For the later two quasars, the 3 valueissignificantlylargerthantheradiuscorrespondingtothewavelengthatwhichthequasar’sflux vanishes. These constraints are tight, with only 10% uncertainties, comparable to those caused by 2 redshift–determinationerrors. Thebest–fitvalue∼sofN˙ are2.1,1.3,and0.9 1057 s−1,respectively, Q v with a factor of 2 uncertainty in each case. These values are a factor of ×5 lower than expected 8 fromthetemplat∼espectrumofTelferetal. (2002). Finally,thebest–fitvalues∼ofx¯ are0.16,1.0,and HI 5 1.0,respectively. The uncertainty in the case of J1148+5251is large,and x¯ is not well constrained. 2 HI However,for both J1030+0524and J1623+3112 we find a significantlower limit of x¯ >0.033. Our 0 HI method is different from previous analyses of the GP absorption spectra of these quas∼ars, and our 1 results strengthen the evidence that the rapid end–stage of reionization is occurring near z 6. 6 ∼ 0 Subject headings: cosmology: theory–earlyUniverse–galaxies: formation–high-redshift–evolution / – quasars: spectrum h p - o 1. INTRODUCTION ization. TheLyαopticaldepthsassociatedwithsuchGP r troughs, and more importantly their accelerated evolu- The epoch of reionization, when the radiation from st early generations of astrophysicalobjects ionized the in- tionatz >6,suggeststhatthereionizationepochisend- a tergalacticmedium(IGM),offersawealthofinformation ing at z ∼6, with a lower limit on the volume-weighted v: aboutcosmologicalstructureformationandphysicalpro- IGM neu∼tral fraction, x¯HI(z 6)>10−3, along at least ∼ i cessesintheearlyuniverse. Onlyrecentlyhavewebegun thelineofsight(LOS)totwoquasa∼rs(White et al.2003; X to gather clues concerning this epoch. The Thomson Fan et al. 2006a). The sharp decline in flux at the edge r scattering optical depth, τ = 0.09 0.03 (Page et al. of the GP trough of quasar J1030+0524 can be used to a e 2006; Spergel et al. 2006b), measured±from the polariza- infer the stronger constraint of x¯HI(z 6) > 0.2, along ∼ tion anisotropies in the cosmic microwave background thatLOS(Mesinger & Haiman2004). Theo∼bservedsize (CMB) by the Wilkenson Microwave Anisotropy Probe of the quasars’ surrounding HII regions can be used to (WMAP), suggests that cosmic reionization began at infer a similar, independent lower limit (Wyithe & Loeb redshift z > 10. However, the current uncertainty of 2004; Mesinger & Haiman 2004), although this requires this value, a∼nd the fact that the optical depth measures additionalassumptionsconcerningthequasaranditsen- only an integrated electron column density, means that vironment, and the inferred neutral fraction scales di- many possible competing reionization scenarios can still rectly with these assumptions (Bolton & Haehnelt 2006; be accommodated. Maselli et al. 2006). The rapid evolution in the sizes of In order to chart the progress of reionization, several the HII regions for different quasars between z 5.7 ∼ attempts have been made to measure the neutral frac- and z 6.4 suggests that the average neutral frac- ∼ tion of the IGM at high redshifts (z > 6). The detec- tion increased by a factor of 10 during this interval ∼ tion of dark spectral regions, so-calle∼d Gunn-Peterson (Fan et al.2006a),althoughtherearealternativewaysto (GP) troughs, in the Lyman line absorption spectra of accountfor the steep size–evolution(Bolton & Haehnelt high-redshiftquasarsdiscoveredinthe SloanDigitalSky 2006). The distribution of spectral sizes of dark gaps Survey (SDSS) offers a probe of the late stages of reion- should help distinguish these scenarios in the future (Gallerani et al. 2006; Fan et al. 2006a). Upper limits 1 YaleCenterforAstronomyandAstrophysics,YaleUniversity, ontheneutralfractionatz 6havealsobeensuggested 260WhitneyAvenue, NewHaven, CT06520 by other means. The lack o∼f significant evolution in the 2 Department of Astronomy, Columbia University, 550 West luminosity function (LF) of Lyman α emitters between 120thStreet,NewYork,NY10027 z 5.7 and z 6.5 was used to infer x¯ (z 6) < 0.3 HI ≈ ≈ ∼ ∼ 2 (Malhotra & Rhoads 2004; Haiman & Cen 2005), if the the sightline-to-sightline fluctuations in these quantities galaxies are assumed to lie in isolation, or x¯ (z 6)< (Mesinger & Haiman 2004). HI ∼ 0.5(Furlanetto et al.2006),ifananalyticprescriptiono∼f The purpose of this paper is to derive these param- the clustering ofHII bubbles (Furlanetto et al.2004b)is eters, i.e. the hydrogen neutral fraction, the size of taken into account. However, a recent, more accurate surrounding HII region, and the ionizing emissivity, by measurementofthe LFatz 6.5revealedevidence ofa separately modeling the spectra of each of the z > 6.2 decrease,weakeningtheseup≈perlimits(Kashikawa et al. quasars. Specifically,weanalyzeJ1148+5251(z =6.42), 2006). Finally, the lack of noticeable absorption due J1030+0524 (z = 6.28) and J1623+3112 (z = 6.22). At to the Lyα damping wing in the afterglow spectrum of present, a 4th quasar is known with a full GP trough. the gamma-ray burst 050904 at z = 6.3 can constrain This source, however, is a broad absorption line (BAL) x¯ (z 6)<0.2 along that LOS (Totani et al. 2006). quasarJ1048+4637(z =6.2),withcomplexspectralfea- HI The∼prep∼onderance of constraints mentioned above tures (Fan, X., private communication) that preclude a suggeststhatthe combinationofallexistingdatais con- simple interpretation in terms of IGM properties. We sistent with x¯ (z 6) 0.1. While this is plausible, have therefore omitted this quasar from our study. HI the constraints liste∼d abo∼ve allrely onvarious model as- Our analysis differs from our previous work sumptions. For example, the conversion of an observed (Mesinger & Haiman 2004) where we modeled the average flux decrement to a neutral fraction (Fan et al. gross transmission features in Lyα and Lyβ close to 2002; White et al. 2003; Fan et al. 2006a) is difficult, the edge of the HII region surrounding J1030+0524. In and requires precise ab–initio knowledge of the IGM the present study, we perform a much more detailed density distribution at high-redshifts (Oh & Furlanetto analysis, by fitting the quasar’s intrinsic emission, and 2005; Fan et al. 2006a). All of the other constraints modeling the distribution of optical depths, pixel–by– mentioned above rely on the IGM density distribution, pixel, blueward of the Lyα line center, following the as well, albeit less sensitively. Techniques involving the procedure proposed in Mesinger et al. (2004). Unlike in sizes of HII regions (Wyithe & Loeb 2004; Maselli et al. Mesinger & Haiman (2004), here we do not make use 2006; Fan et al. 2006a) are subject to degeneracies be- of the Lyβ absorption spectrum, except to set a lower tweenquasarlifetime,ionizingflux,andtheIGMneutral limit in the allowed size of the HII region surrounding fraction. Further complicating matters for the HII bub- J1030+0524, the only one among the three quasars ble size techniques is the fact that the onset of the GP that has a noticeable offset between the Lyα and Lyβ troughneednotcorrespondtotheedgeoftheHIIregion. GP troughs. Our reason for ignoring the Lyβ region is The GP through just corresponds to the region where that the majority of the pixels we analyze (see Figure 2 the flux drops below the detection threshold, and thus below) have flux detected in the Lyα absorption region, merely provides a lower limit on the size of the HII re- which can directly be converted to a Lyα optical depth. gion (Mesinger & Haiman 2004; Maselli et al. 2006; see This avoids the need for a statistical modeling of the also Bolton & Haehnelt 2006 for a recent, comprehen- foregroundLyα absorption, necessary to convert Lyβ to sive review of these issues). Furthermore, determining a Lyα optical depth. thenumberdensityandevolutionofLyαemittinggalax- This paper is organized as follows. In 2, we describe § ies is complicated and depends on knowledge of galaxy our analysis technique, including generating the simu- clustering(Haiman & Cen2005;Furlanetto et al.2004a, lated absorption ( 2.1), fitting for the quasars’intrinsic § 2006),absorberclustering,aswellasthelarge-scaleden- emission spectrum ( 2.2), and the statistical method § sity (e.g. Barkana & Loeb 2004b,a) and velocity fields used to compare the observed and simulated spectra (e.g. Cen et al. 2005), which are all unconstrained at ( 2.3). In 3, we present our results for each of the § § the relevant high redshifts. quasars we analyze. In 4, we discuss various aspects § More generally, the drawback of many of these meth- of our results, such as the origin of the constraints we ods is that they rely on the combination of data from obtain, and the robustness of our results in the face of multiple LOSs to obtain statistical significance. How- uncertaintiesinthequasartemplatespectrum. In 5,we § ever, assuming constancy in the IGM properties along summarizeour key findings andpresentour conclusions. disparateLOSsis dangerousathigh-redshifts,especially approachingtheendofreionization,whenlargesightline- 2. ANALYSIS to-sightline fluctuations are expected in the density field In this section, we describe the three components of andthe ionizingbackground(see,e.g., recentreviewsby our analysis: (i) a hydrodynamical simulation to model Djorgovskiet al. 2006and Lidz et al. 2006). Indeed, the theabsorption( 2.1);(ii)modelingthequasars’intrinsic current sample of high-redshift quasars already reveals emission spectru§m ( 2.2); and (iii) the statistical com- large fluctuations in the IGM properties along different parisonofobserveda§ndsimulatedspectra( 2.3). Unless LOSs (Fan et al. 2006a). This would suggest that the stated otherwise, we adopt the backgroun§d cosmologi- most rewarding technique would be able to analyze and cal parameters (Ω , Ω , Ω , n, σ , H ) = (0.76, 0.24, Λ M b 8 0 obtain constraints from each source independently from 0.0407, 1, 0.76, 72 km s−1 Mpc−1), consistent with the other sources. This is statistically less powerful (e.g. three–yearresults by the WMAP satellite (Spergel et al. Mesinger et al. 2004), but the rewards can be rich, po- 2006a). Weuseredshift,z,andtheobservedwavelength, tentially providing additional constraints on the proper- λ , interchangeably throughout this paper as a mea- obs tiesofindividualquasarsandtheirenvironments,suchas sureofthedistancealongthe lineofsight. Unlessstated thequasar’semissivityofionizingphotons,sizeofitssur- otherwise, we quote all quantities in comoving units. rounding HII region, sharpness of the HII regionbound- ary, the quasars’X-ray emission, as well as a measure of 2.1. Simulating Absorption in the IGM 3 In order to simulate the absorption spectra, we make Since high-redshift quasars are surrounded by their useofahydrodynamicalsimulationboxinaΛCDMuni- own highly ionized HII regions, it is useful to express verse at redshift z = 6. The details of the simulation τ as a sum of contributionsfrom inside (τ ) andoutside R are described in Cen et al. (2003), and we only briefly (τ ) the HII region, τ = τ +τ . Defined as such, the D R D discuss the relevant parameters here. The box is 11 h−1 resonant optical depth, τ , would be given by Eq. (1), R Mpc on each side, with each pixel being about 25.5 h−1 with the lower limit of integration changed to the red- kpc. This scale resolves the Jeans length in the smooth, shiftcorrespondingto the edge ofthe HII region,z(R ). S partially-ionized3 IGM by more then a factor of 10. The The damping wing optical depth τ , would be given by D box also has a slightly different set of cosmological pa- Eq. (1), with the upper limit of integration changed to rameters: (Ω ,Ω ,Ω ,n,σ ,H )=(0.71,0.29,0.047,1, z(R ) 5. Λ M b 8 0 S 0.854 , 70 km s−1 Mpc−1). Since cosmic hydrogen is highly ionized at low red- Density and velocity information was extracted from shifts, one can replace the lower limit of integration in the simulation box along randomly chosen lines of sight. Eq. (1) with z , denoting the redshift below which HI end The step size along each LOS was taken to be 5.1 absorption becomes insignificant along the LOS to the h−1 kpc, which resolves the Lyman α Doppler width source. In our analysis, we chose the value of z sep- end in the partially-ionized IGM by more than a factor of arately for each quasar to correspond approximately to 40. The exact value of the step size was chosen some- the blue edge of its GP trough. This roughly translates what arbitrarily, and does not influence the results as toz 0.2–0.3forourthreequasars. We findthatthe end ∼ long as it adequately resolves the Doppler width. At exactchoiceofz isnotimportant,becausemostofthe end each step, the density and velocity values were aver- contributiontoτ comesfromneutralhydrogencloseto D aged for the neighboring pixels, and weighted by the z(R ). Inside the HII region, we assume an IGM tem- S distance to the center of the pixels. We extended each perature of T = 2 104 K. Outside the HII region, for × LOS by the common practice of randomly choosing a thecaseofamostlyionizeduniverse,weassumeanIGM LOS through the box, and stacking several segments temperature of T = 104 K, while the temperature in a together (Cen 1994; Mesinger et al. 2004; Kohler et al. neutraluniversewastakentobeT =2.73 151(1+z)2K, × 151 2005; Bolton & Haehnelt 2006). The r.m.s. mass fluc- valid for z < 150 (Peebles 1993). Our results are fairly tuation on the scale the box is σ(z = 6.3) 0.15. insensitive to the exact temperature used. The Doppler ∼ Thisisasmallenoughvaluethatneglectingwavelengths width of the Lyman-α absorption cross section scales as longer than the box size is justified in our analysis. ν T1/2, but the totalintegratedarea under the cross D Note however, that this is not the case for modeling sect∝ion is independent of temperature. highly non-linear processes, such as reionization (e.g. Assumingionizationequilibrium,wecomputetheneu- Barkana & Loeb2004b). WhenusingsuchaLOSinsim- tralhydrogenfractionateachpointalongthe LOSwith: ulatingagivenquasarspectrum,densityvaluesfromthe simulation were scaled as (1+z)3. nx (Γ +Γ )=n2(1 x )2α , (2) HI Q BG HI B The total optical depth due to Lyman α absorption, − τ, between an observerat z = 0 and a source at z = zQ, where ΓQ and ΓBG are the number of ionizations per at an observedwavelengthof λ =λ (1+z ), is given second per hydrogen atom due to the quasar’s and the obs 0 Q by: background’s ionizing fluxes, respectively, and αB is the case-B hydrogen recombination coefficient. The LHS of Eq. (2) accounts for the number of hydrogenionizations zQ dt λ per cm3, while the RHS accounts for the number of hy- obs τ(λobs)=Z0 dzcdz n(z)xHI(z)σ"(1+z)(1− v(cz))# donroegoebntareincos:mbinations per cm3. Solving Eq. (2) for xHI (1) where n(z), v(z), xHI(z) are the total hydrogen number b √b2 4 density,thecomponentofthepeculiarvelocityalongthe xHI = − − − , (3) 2 LOS, and the local hydrogen neutral fraction, respec- tively, at redshift z. The Lyα absorption cross section, where b 2 (Γ + Γ )/(nα ). We further Q BG B ≡ − − to first order in v(z)/c, is given by σ[λ /(1+z)/(1 parametrize the quasar’s ionization rate with: obs − v(z)/c)]. Γ f Γe (4) 3ThetemperatureofaneutralIGMbeforeanyionizationwould Q ≡ Γ Q beseveral orders of magnitude lower than inthe partially-ionized where Γe = (4πr2)−1 ∞1.55 1031 (ν/ν )−1.8 [(1+ case,andtheJeanslengthintheIGMwouldbeunder-resolvedby Q νH × H af4acTtohreofth∼r6e.e–year WMAP results yielded a smaller value of zb)e/in(1g+thzeQi)o]−n0iz.8a(tσio/nhνfr)edRqνuse−n1c,ywoifthhyνdHro=ge3n.2a9n×d1r0b15eiHngz σ8=0.76±0.05(Spergeletal.2006b). Asmallerσ8 wouldresult the luminosity distance between the source and redshift in a narrower distribution of optical depths than those obtained z. Γe results from redshifting a power–law spectrum from the box we use. Although this does not exactly mimic the Q effect of increased damping wing absorption from a more neutral with a slope of ν−1.8, normalizedsuch that the emission universe,whichshiftstheopticaldepthdistributiontohigherval- ues, it is likely that if confirmed, a smaller σ8 would weaken our 5 Notethattheterms“resonant”and“dampingwing”areonly lowerlimitsontheIGMneutral fractionbelow. Estimatesforthe appropriate in describing the contributions of τR and τD inside sisenqsuitainvtitityatoifveoluyrdreiffisuclutsltowniothuoructhoruicnenoinfgcoasmsuoiltoegiocfalmpoadrealms.etWeres the HII region (i.e. at λobs >∼ λα[1+z(RS)]). Outside of the note however that analysis combining the three–year WMAP re- HIIregion(λobs<∼λα[1+z(RS)]),τD actuallyintegratesoverthe sultswithdatafromtheLyman–αforestsuggestsahighervalueof resonantpartoftheabsorptioncrosssectionandτRintegratesover σ8=0.86±0.03(Lewis2006),whichisconsistentwithourchoice. thedampingwing(seeFig. 1). 4 rate of ionizing photons per second is 1.3 1057, match- 1. R , the radius of the quasar’s HII region, S × ing the emissivity one obtains (Haiman & Cen 2002) for SDSS J1030+0524 using a standard template spectrum 2. xIGM, the neutralhydrogenfractionin the IGM at HI (Elvis et al. 1994). Note that the normalization inferred mean density, from the more recent template of Telfer et al. (2002) is 3. f , the quasar’s ionizing luminosity, in units of a factor of of 5 higher. This roughly translates to an Γ ∼ 1.3 1057 s−1. emission rate of ionizing photons of × To obtainanintuition forthe effect of varyingthese free N˙ =f 1.3 1057 s−1 . (5) Q Γ parameters, in Figure 1 we plot τ (solid curve) and × × R τ (dashed curve) for a typical simulated LOS, for pa- We ignore radiative transfer along the LOS (and also D theionizationofhelium). Instead,outsideofthequasar’s rameter choices (RS, xIHGIM, fΓ) = (40 Mpc, 1, 1). As HII region, we compute x assuming the quasar con- stated above, the total Lyman α optical depth is the HI tributes nothing to the ionizing flux, i.e. ΓQ(R>RS)= sum of these two contributions, τ = τR + τD. Roughly 0,while inside the HII region,weassume the gasis opti- speaking,changing RS movesthe dashed (τD) curve left callythininthe ionizingcontinuum. Thisassumptionof and right,while changingxIGM moves this curve up and HI a sharp HII region boundary is an approximation, valid down in Figure 1. Changing fΓ approximately corre- in the regime where either the quasar’s spectrum is soft, spondstomovingthesolid(τR)curveupanddown. The or where a strong soft ionizing background flux domi- dotted line demarcates roughly the maximum total op- nates over that of the quasar already at radii smaller tical depth, τ 5.5, at which a signal can be detected ∼ that the mean free path of the typical ionizing quasar withaconfidencegreaterthan1-σ,inthespectralregion photon (as is the case of the proximity effect at lower just bluewardof the Lyα emissionline (see discussion in redshift; Bajtlik et al. 1988) If the quasar has a hard 2.3). Finally, after the optical depth and the associ- spectrum, and the background intensity is low, then the §ated transmission (e−τ) is calculated according to the edge of the HII region will be “blurred”. A similar blur- procedure outlined above, the transmission is averaged ring would arise if additional ionizing sources inside the over 3.3 ˚A wavelength bins to match the Keck ESI ≈ HII region, with an extended spatial distribution, would spectra with which the simulated spectra are compared contributesignificantflux (Wyithe & Loeb2006),a pos- (note that the actual smoothing length ∆λobs differs by sibility we do not address here. However,already strong a few percent for each source, since they are at different upper limits can be placed on the extent of such blur- redshifts). ring from the offset between the wavelengths of the Lyα ThefreeparameterRS wasvariedinlinearincrements and Lyβ GP troughs. Since Lyβ absorption is less sen- of∆RS =2Mpc,withthelowerlimitsetbytherededge sitive than Lyα, the amount of offset between the GP of the Lyα GP through in the observed spectra. The troughs gives a measure of the spatial (i.e. wavelength) otherparameterswerevariedintheranges: 3.2 10−4 evolution of τ close to the edge of the HII region, with xIHGIM ≤1, in multiples of 5, and 0.1≤fΓ ≤10,×in linea≤r a large offset indicating slow τ evolution, and a small increments of ∆fΓ = 0.3 between 0.1 fΓ 2 and ≤ ≤ offset indication rapid τ evolution. Mesinger & Haiman ∆fΓ =2 between 2<fΓ 10. ≤ (2004) studied this forthe case ofJ1030+0524,but even stronger constraints could be placed on the spectra of 2.2. Modeling the Intrinsic Emission Spectrum J1148+5251 and J1623+3112, since they do not show As observational input in our analysis, we make any offsets between the GP troughs (Kramer et al., in use of Keck ESI spectra of the three highest redshift preparation). Furhtermore, neglecting possible obscura- quasars discovered to date: J1148+5251 (z = 6.42), Q tioneffectsarisingfromopticallythickgasanddustclose J1030+0524 (z = 6.28), and J1623+3112 (z = 6.22). Q Q to the quasarcan bias our determination of the quasar’s These quasars are the only quasars discovered to date luminosity to lower values. This means that our deter- with z > 6.2, as well as the only quasars (aside from Q mination of the quasar’s luminosity should be viewed as the BAL quasar J1048+4637) exhibiting complete GP theobservedluminosityalongthe LOS,andequaltothe troughs, with no detectable flux over a wide wavelength intrinsic luminosity in the low opacity limit. range (Fan et al. 2006b). Readers interested in the de- Finally, we define the neutral hydrogenfraction in the tailsofthecorrespondingobservationsareencouragedto IGM, xIHGIM, with Eq. (3) using the mean density of the consult Fan et al. (2006b) and White et al. (2003). universeattheedgeoftheHIIregion,n=n0[1+z(RS)]3, Inordertocomparethequasars’observedspectrawith and ΓQ = 0. Although ΓBG is the more physically rel- our simulated spectra, we must know the intrinsic emis- evant quantity and we use it in our analysis, henceforth sionspectrumfromeachofthe observedsources. Obser- we willuse xIHGIM as the proxyfor ΓBG in orderto match vations of lower redshift quasars (e.g. Telfer et al. 2002) convention. As defined above, xIGM represents the neu- suggestthattheirintrinsiccontinuumemissioniswellfit HI tral fraction at the mean density. The more common by a broken power–law, with a mean spectral slope of representation of the neutral fraction in the literature is α 0.7(f να)redwardoftheLyαline,andaslope ν ∼− ∝ inthe formofthevolume(ormass)weightedmean,x¯ . ofα 1.8bluewardof the Lyα line. The high-redshift HI ∼− This definition, of course, depends on the assumed den- SDSS quasars do not show any evidence for deviating sity distribution. Hence, for comparison purposes, we from this trend (e.g. White et al. 2003). The precise express our results both in terms of xIGM and x¯ , with location of the wavelength of the break in the power– HI HI the latterobtainedby averagingoveroursimulationbox law is difficult to determine for any given source, due densities and assuming ionization equilibrium. to a strong broad Lyα emission line superimposed on To summarize, our analysis has three free parameters: the continuum spectrum. However,for our purposes, we 5 the fits to the Lyα line of J1148+5251and J1623+3112. However,the Lyα line ofJ1030+0524requirestwogaus- sians for an accurate fit, as it contains a broad and nar- row component, similar to the shape of a Voigt profile (Ho et al. 1997; Greene & Ho 2004). The redshifted line centersofthegaussianswereheldfixedinthefittingpro- cedure, set to (1+z )λ , where λ =1215.67˚A for Lyα Q 0 0 andλ =1240.81˚AforNV.Theheightandwidthofthe 0 gaussianswerethenfoundbyfitting the fluxinpixelson the red side of the line center. Care was taken, in addi- tion, to excise the spectral region immediately redward of the Lyα line center, as this region could be suscepti- ble to non-negligibleabsorptionby dense infalling gas in the foreground, close to the quasar (e.g. Mesinger et al. 2004; Barkana & Loeb 2004a). Indeed, the quasar spec- tra do exhibit sharp drops in the observed flux as one approaches the line center from the red side (see Fig. 2; see also Barkana & Loeb 2003). This region was excised by eye from the fitting procedure. The amplitude of the continuumwasfoundseparately,byfittingthefluxinre- gionsofthespectrumfurthertothered,thatarea-priori knownnottosuffermetallineabsorption,corresponding to roughly 1300 – 1350 ˚A in the rest frame. We explic- itly avoid spectral sub–regions with obvious emission or absorption lines. More specifically, the spectral regions Fig. 1.— Model from a hydrodynamical simulation for the in the observed frame used in the fit for the continuum opticaldepthcontributionsfromwithin(τR)andfromoutside(τD) the local ionized region for a typical lineof sight towards a zQ = were (demarcated with horizontal lines in Figure 2): 6.42quasar. Thisfigure was created withparameter choices (RS, xanIHGdIMth,efΓs)ol=id(c4u0rvMepcco,rr1e,s1p)o.nTdshetodaτRsh.eTdhceurtvoetacloLrryemspaonndαsotpotiτcDal, • J9510104˚A8+,59275010:–10090100˚A0–;9125˚A, 9128–9190˚A, 9220– depth is the sum of these two contributions, τ = τR + τD. The dotted linedemarcates roughlythemaximumtotal optical depth, J1030+0524: 8880–8990˚A, 9020–9100˚A, 9840– tτh∼an51.5-,σa,tinwthhicehsapescigtrnaallrceagniobnejduesttebctluedewwairtdhoafctohnefiLdeynαceemgriesasitoenr • 9930˚A, 9960–10020˚A; line(seediscussionin§2.3). Inouranalysis,thetransmission(e−τ) isaveragedover∼3–4˚Awavelengthbins(1/10oftheresolution J1623+3112: 8830–8895˚A, 8980–9040˚A, 9790– inthefigureabove),thusdecreasingthefluctuationintheeffective • 9863˚A, 9880–9907˚A; optical depth. For reference, the redshifted Lyα wavelength is at 9020˚A,offtheplot. InFigure2,weshowthedatafortheobservedspectra, F(λ ), with 1-σ error bars as well as the continuum + obs areinterestedinthe opticaldepthatwavelengthswithin Lyα + N V fit to the intrinsic emission, F0(λobs), gen- just a few tens of rest-frame Angstroms blueward of the erated by our procedure (solid curve). The Lyα optical Lyαlinecenter,wheretheLyαemissionlineisstilldom- depth can then be inferred in each pixel separately as: inant over the continuum. Therefore, we do not model τobs(λobs)= ln[F(λobs)/F0(λobs)]. − thebreak,anduseasinglepower–lawcontinuumcompo- nent instead in our emission template. Throughout our 2.3. Comparing Simulated and Observed Spectra analysis, we keep the logarithmic slope α = 0.7 fixed After the observed optical depths, τ (λ ), are ob- − obs obs (aswewilldemonstratebelow,ourresultsareinsensitive tained as discussed in 2.2, and the simulated optical to this choice). depths,τ (λ ),areo§btainedforeveryLOSandpoint sim obs We deriveourconstraintsbelow entirelyfromthe blue in our parameter space as discussed in 2.1, we wish side of the Lyα line in our analysis, where the effects of to statistically compare the simulated an§d observed op- theIGMcanbefeltmoststrongly(Mesinger et al.2004). tical depths for each quasar. More precisely, for each To this end, we make use of “clean” spectral regions on quasarandeverypoint in our three–dimensionalparam- the redsideofthe Lyαline,wherethe resonantattenua- eter space, (R , xIGM, f ), we want to answer the ques- S HI Γ utiaotnioins niseglelisgsibthlean(τR10∼%0()e,−aτnDd>th0e.9d)aemvpeninignwthinegcaatsteeno-f utieosn:anwdhathteislisthteofliτkelih(oλod t)hvaatlutheeslwisetreofdτroabwsn(λforbosm) vtahle- sim obs afullyneutralIGMandasmallRS. WefittheNVemis- same underlying probability distribution? Note that we sionlineat(1+zQ)1240.81˚Awithasinglegaussian,and have a fixed list of about 30 τobs(λobs) values for each thecontinuumwithasinglepower–law. Additionally,we quasar,butwegenerateadifferentlistof3000τ (λ ) sim obs fit the Lyα line with a single gaussian for J1148+5251 values, combining 100 LOSs, for each choice of (R , S ((wwiitthh aa fifitttteeddwwiiddtthhooff∆∆λλoobbss==8111˚A3˚A) )anadndwiJt1h62a3s+u3m11o2f xIHGWIMe, fuΓs)e. the Kolmogorov-Smirnov (K-S) test (e.g. two independent gaussians for J1030+0524 (with fitted Press et al.1992)to comparethe cumulativeprobability widths of ∆λ = 45˚A and ∆λ = 113˚A). We have distributions (CPDFs) – i.e., histograms – of τ (λ ) obs obs obs obs found that adding a second gaussian does not improve and τ (λ ). The usefulness of the K-S test is that sim obs 6 4 6 10 3 8 4 2 6 1 2 4 2 0 0 0 8900 8950 9000 8750 8800 8700 8720 8740 Fig. 2.—Observedspectra(shownaspointswith1-σerrorbars)andthecontinuum+Lyα+NVfitoftheintrinsicemissionspectrum generated byourprocedure(solid curve). Thespectralregions(ontheredsideoftheLyαlinecenter)that weusedinthefitaremarked with horizontal lines. The inlaid panels zoom in on the spectral region (on the blue side of the Lyα line center) that were then used to compare the absorption optical depth τ with hydrodynamical simulations. The vertical dashed lines bracket the wavelength ranges used toconstructτ histograms(see§2.3). ThespectrashownarethoseofJ1148+5251, J1030+0524, J1623+3112 (left to right). it provides a figure of merit which can be directly in- formation from the wavelength dependence is fully pre- terpreted as the likelihood that two distributions were served; however,only a single observedpixel is available drawn from the same underlying distribution, without to compare with a corresponding pixel from each LOS, theneedforthead-hocbinningofdata,andtheresulting rendering the use of the K-S test unreliable. The con- loss of information, required by such methods as the χ2 verse is true in the other extreme, where the wavelength test. ThedrawbackofusingtheK-Stestisthatcomput- binsizeapproachestheentirespectralregionweanalyze: ing absolute confidence levels around our best fit values the shape of the CPDF is accurately captured, but in- would require prohibitively expensive Monte-Carlo sim- formationfromthewavelength–dependenceisdestroyed. ulations. We caution the reader that, as a result, the Here we (somewhat arbitrarily) divide the spectral re- relative likelihoods we quote below do not representtra- gion of analysis into wavelength bins of width 25–30 ditional “error–bars”. ˚A, so that each bin contains between 6–8 pixe∼ls. The Our statistical analysis is quite similar to that em- two blue–most bins were chosen bracket the edge of the ployed by Rollinde et al. (2005) to model density struc- GP trough in the observed spectrum. Furthermore, we turesurroundingmoderateredshift(z 2)quasars. The remove the region 20–30 ˚A blueward of the line center ∼ main difference in technique is that while our present from our analysis, corresponding to the expected mean approach uses a parameterized model, Rollinde et al. radius of the large–scale overdense region surrounding (2005) fit for the intrinsic spectrum directly from the such quasars (Barkana & Loeb 2004a). This step was low–τ regions of the spectrum. Using simulated spec- necessary,as our simulation box size is too smallto con- tra, they show that by applying the K-S test on optical tain a statistical sample (or even a single occurrence) of depth CPDFs, one can statistically discriminate among suchlargeoverdensities,requiredfortheiraccuratemod- different models of the distance–evolution of the optical eling. The wavelength bins we use in our analysis are depth. In Mesinger et al. (2004), we use a similar statis- shownbythe verticaldashedlinesinthe inlaidpanelsof ticaltesttoshowthatthemoderatedegeneracybetween Figure2. We verifythatourresultsarefairlyinsensitive R and xIGM can be broken with a single spectrum, if to our choice of bin size, i.e. the relative K-S likelihoods S HI R ismoderatelyconstrained. Thisprioranalysisdidnot between models change little with decreasing bin size. S take into account the additional constraint arising from In order to incorporate an uncertainty into our emis- the evolution of the CPDF with distance, as our current sion template, we add Gaussian-distributed, uncorre- analysis does. Nevertheless, we caution the reader that lated flux variations to each pixel in the template. We thesetestsonlyconfirmtheabilitytostatisticallyrecover do this by regarding our fitted value of the template, input parameters using our statistical analysis. They do F (λ ), asthe meanvalue ofthe true flux distribution, 0 obs not test the validity of our model. andthe“true”emissionatλ isdrawnfromaGaussian obs The optical depth CPDF should have a strong depen- distributionaroundthatmean,withastandarddeviation dence on wavelength (since the gas is more highly ion- σ (λ ). Hence, we replace every τ (λ ) value by flux obs sim obs ized near the quasar). In principle, the information con- a set of 1000 values, defined by tainedinthisdependencecouldbeextractedandusedto further constrain model parameters. For example, one τsampled(λ ) ln[e−τsim(λobs) A] , (6) could compare the CPDFs of τ (λ ) and τ (λ ) sim obs ≡− × obs obs sim obs in several wavelength bins for each quasar. Unfortu- where A represents our emission uncertainty, and is nately, in practice, one also needs to have many points drawn from a Gaussian distribution with a mean of 1, in each individual bin, to be able to sample the under- and a standard deviation of σ (λ )/F (λ ). We flux obs 0 obs lying distribution in that bin. In one extreme, as the set σ (λ )/F (λ ) = 0.3, which corresponds ap- flux obs 0 obs wavelength bin size approaches the pixel width, the in- proximately to the typical quasar-to-quasar scatter in 7 the continuum level immediately redward of Lyman α bottom–right with increasing f . Both of these trends Γ (Vanden Berk et al. 2001). However, we have verified areevidentinFigure3. Forexample,thecoloredsquares that our results are not sensitive to this choice. For ex- (correspondingtoK-Slikelihoodsthatarewithinafactor ample, in the case of J1623+3112, we find that in the of27ofthepeaklikelihood)inthemiddlerowcovermost range 0.1 < σ (λ )/F (λ ) < 0.5, the peak K-S of the displayed parameter space in the left–most panel, flux obs 0 obs probability∼remains unchanged, and∼the confidence lim- though there are no pixels with high likelihoods. In- its quoted below change at most by a single pixel in our creasingthequasar’sionizingflux(i.e. goingthroughthe parameter space. panels from left to right), shifts the colored squares in- Finally,incomparingspectra,oneneedstodefinewhat creasinglytowardsthe bottom–rightcornerofthe panel. oneconsiderstobe“zeroflux”,assimulationshavenear- infinite precision and observations do not. To this end, 3.1. J1148+5251 wedefine amaximumallowedvalueoftheopticaldepth, Likelihood estimates for J1148+5251 (z = 6.42) can Q MAX[τ] 5.5. This approximatelycorrespondstoa1-σ beseeninthetoprowofthreepanelsinFig. 3. Theleft, ≡ detectionin ourobservedspectra atwavelengthsaround center,andrightpanels correspondto valuesoff = 0.7, Γ the edge of the GP trough, where the majority of pixels 1.6, and 1.9, respectively. The peak likelihood of 3.0% with fluxes less than this value are located (see Fig. 2). occursat(R ,xIGM,f )=(40Mpc,0.2,1.6). Thisisthe S HI Γ In our comparison, all values of τ > 5.5 are treated the smallest peak likelihood we obtain (c.f. peak likelihoods same, regardless of their actual value. of 30% for the two quasars below), and most likely To summarize, for each quasar we divide the spectral mea∼ns that our model isn’t as accurate in describing the region roughly blueward of (1+zQ)λα–25˚A and red- physical spectra of this particular quasar. The volume– ∼ ward of ∼ [1+z(RS)]λα–25˚A into 3–5 bins of approx- averagedneutralfractionat the peak is x¯HI =0.16. The imately equal wavelength width, containing 6–8 pixels rangeofparametervalueswhoselikelihoodestimatesare each. In each bin we compare the CPDFs of τobs(λobs) withinafactorof3ofthepeaklikelihoodisencompassed andτsample(λ )(thelatterconcatenatedover100simu- by: sim obs latedLOSs). Foreverypointinour3Dparameterspace, 40 Mpc <R < 42 Mpc, these two distributions are compared via the K-S test, • S ∼ ∼ which produces a probability to be interpreted as the xIGM < 1.0 or x¯ < 1.0, likelihood that the two distributions were drawn from • HI HI ∼ ∼ the same underlying distribution. Finally, the K-Sprob- 0.7 <f < 1.9. abilities for each wavelength bin are multiplied to give a • Γ ∼ ∼ likelihood estimate for that point in the 3D parameter Note that here as well as below, the end points of the space. quoted ranges correspond to the minimum or maximum value of the parameterfor whichthere exists a combina- 3. RESULTS tionoftheothertwoparametersgivingaK-Sprobability Our results for the likelihoods of various parameter– within a factor of three from the peak probability. Note combinations are shown in Figure 3 for J1148+5251 that these ranges are similar in spirit to “marginalized (top), J1030+0524 (middle), and J1623+3112 (bottom). errors”,inthattheotherparametersareallowedtovary, Thecrossineachcasemarksthepointinthe3Dparam- rather than fixed (but they are not actual marginalized eter space where the likelihood peaks. In order from errors;see caveat mentioned in the footnote above). lighter to darker, the yellow, green, and blue squares correspond to points in our parameter space with like- 3.2. J1030+0524 lihoods that are a factor of 27, 9, and 3 below the peak Likelihood estimates for J1030+0524 (z = 6.28) are Q likelihood, respectively.6 Parameter combinations with shown in the middle row of panels in Figure 3. The left, likelihoodssmallerthan1/27thofthe peak areshownas center,andrightpanels correspondto valuesoff = 0.4, Γ empty squares. 1.0, and 1.3, respectively. The peak likelihood of 34% Below, we will discuss each quasar in turn. However, occurs at (R , xIGM, f ) = (41 Mpc, 1, 1.0). The range S HI Γ it is encouraging to note first that the isocontours in ofparametervalueswhoselikelihoodestimatesarewithin Figure 3 behave as expected, following the moderate de- a factor of 3 of the peak likelihood is encompassed by: generacies in the parameters. Namely, since a stronger damping wing (τ ) contribution can be obtained either 37 Mpc <R < 45 Mpc, D S • by increasing the neutral fraction or by decreasing the ∼ ∼ sizeoftheHII region,theisocontoursshouldbe oriented xIGM < 1.0 or x¯ < 1.0, • HI HI frombottom-left to upper-rightineachpanel. Similarly, ∼ ∼ sinceasmallerionizingflux(i.e. greaterτR)canroughly • 0.7 <fΓ < 1.9. be compensated for in the total optical depth, τ +τ , ∼ ∼ R D by increasing the size of the HII region or decreasing As mentioned earlier, J1030+0524 has a unique (at the neutral fraction (both decrease τ ), the isocontours least among our sample of three quasars) feature. D should shift monotonically from the upper–left to the Namely,theLyβGPthroughisnoticeablyoffsetfromthe LyαGPthrough. Themerepresenceofflux(irrespective 6 In order to obtain absolute confidence limits, we would need of its properties) in the Lyβ region of the spectrum can to performa set of computationally expensive Monte-Carlosimu- be used to set a minimum size of the surrounding HII lations. For simplicity, we avoid such computations, and instead region,R >40 Mpc. The regionexcluded by this prior contendourselvesbyquotingresultsintermsoftherelativeoffset S fromthepeaklikelihood. isshadedov∼erinthemiddlerowinFigure3. Takingthis 8 Fig. 3.— Likelihoods for J1148+5251 (top), J1030+0524 (middle), and J1623+3112 (bottom) in the 3–dimensional parameter space of (RS, xIHGIM, fΓ). Each panel shows contours ina2D sliceof the 3D likelihood surface. Thecross marks the parameter combination with the peak likelihood value (see § 2.3 for how the likelihoods were derived). In order from lighter to darker, the yellow, green, and blue squares correspond to points inour parameter space withlikelihoods that are afactor of 27, 9, and 3below the peak value, respectively. Parametercombinationswithlikelihoodsbelow1/27thofthepeakareshownasemptysquares. J1148+5251: Thepeaklikelihoodof3.0% occurs at (RS, xIHGIM, fΓ) = (40 Mpc, 0.2, 1.6). The left, center, and right panels correspond to values of fΓ= 0.7, 1.6, 1.9, respectively. J1030+0524: Thepeaklikelihoodvalueof34%occursat(RS,xIHGIM,fΓ)=(41Mpc,1,1.0). Theleft,center,andrightpanelscorrespond tovaluesoffΓ=0.4,1.0,1.3,respectively. TheshadedregionshowsvaluesruledoutbythepresenceoffluxintheLymanβ region,which setstheadditionalconstraintRS >∼40Mpc. J1623+3112: Thepeaklikelihoodvalueof39%isat(RS,xIHGIM,fΓ)=(29Mpc,1,0.7). The left,center,andrightpanelscorrespondtovaluesoffΓ=0.4,0.7,1.0,respectively. lowerlimitonR intoaccount,therangeofvalueswhose 3.3. J1623+3112 S likelihood estimates are within a factor of 3 of the peak Likelihood estimates for J1623+3112 (z = 6.22) are Q likelihood shrinks to: shown in the bottom row of panels in Figure 3. The left,center,andrightpanelscorrespondtovaluesoff = 41 Mpc <R < 45 Mpc, Γ • S 0.4, 0.7, and 1.0, respectively. The peak likelihood of ∼ ∼ 39% occurs at (R , xIGM, f ) = (29 Mpc, 1, 0.7). The 0.04 <xIGM or 0.033 <x¯ , S HI Γ • HI HI rangeofparametervalueswhoselikelihoodestimatesare ∼ ∼ withinafactorof3ofthepeaklikelihoodisencompassed 0.7 <fΓ < 1.3. by: • ∼ ∼ 9 solidblue,long-dashedred,short-dashedgreen,anddot- tedpurple curvescorrespondtoparametercombinations of (R , xIGM, f ) = (29 Mpc, 1.0, 0.7), (29 Mpc, 0.04, S HI Γ 0.7), (33 Mpc, 1.0, 0.7), and (29 Mpc, 1.0, 0.1), respec- tively. The crosses are located at the set of observed optical depth values. The solid blue curve corresponds to the parameter combination with the peak likelihood value, while the parameter combinations represented by the other curves have likelihood values below 1/27th of the peak. The sharp jump corresponding at the largest τ binandtheoverlayingofsevencrossesinthetoppanel areduetothefinitedetectionthreshold,wherebyallval- ues of τ > 5.5 are indistinguishable from one another, and hence fall in the same bin in this figure. The fig- ure clearly demonstrates that our lower limit on xIGM HI arises from the low–τ tail predicted in models with low neutralfractionsorthosecorrespondingtoregionsofpa- rameter space which could mimic a low neutral fraction (seethediscussiononisocontoursin 3). Whilethesolid § bluecurvesareconsistentwiththedistributionofpoints, these other curves all have tails that are “too long”: no points corresponding to these tails are seen in the data. Note also that while some models might perform well in onewavelengthbin,onlythesolidbluecurveisconsistent with the observed data in all wavelength bins. Fig. 4.— Distributions of simulated optical depths for the three wavelength bins used in our analysis of J1623+3112. The 4.1. Sanity Checks solidblue,long-dashedred,short-dashedgreen,anddottedpurple curvescorrespondtoparametercombinationsof(RS,xIHGIM,fΓ)= As our analysis technique is as yet untested, it is use- (29Mpc,1.0,0.7),(29Mpc,0.04,0.7),(33Mpc,1.0,0.7),and(29 ful to take a step back, and verify that our results make Mpc, 1.0, 0.1), respectively. The crosses are located at the set of sense. Asalreadymentioned,thelikelihoodiso–contours observedopticaldepthvalues. Thesolidbluecurvecorrespondsto behave as expected. Here we discuss several other en- the parameter combination with the peak likelihood value, while the parameter combinations represented by the other curves have couragingchecksonthevalidityofouranalysis. Wecau- likelihoodvaluesbelow1/27thofthepeak. Thesharpjumpatthe tion the reader, however, that these are not meant to largestτbinandtheoverlayingofsevencrossesinthetoppanelare be a proof of the accuracy of our analysis, and instead duetothefinitedetectionthreshold,wherebyallvaluesofτ >5.5 should be regarded merely as consistency checks. areindistinguishablefromoneanother,andhencefallinthesame bin in this figure. The figure demonstrates that our lower limit First, the procedure here recovers fairly well the re- on xIGM comes from the low–τ tail predicted in models with low sult of our previous analysis in the case of J1030+0524 HI neutral fractions or those corresponding to regions of parameter (Mesinger & Haiman 2004). This need not be the case, spacewhichcouldmimicalowneutralfraction(seethediscussion since the nature of the two analysis procedures are dif- on isocontours in § 3). No points corresponding to these tails are seeninthedata. Notealsothatwhilesomemodelsmightperform ferent; namely, our previous work focused only on gross, wellinone wavelength bin, only the solidbluecurve isconsistent approximateLyαandLyβ transmissionfeaturesnearthe withtheobserveddatainallwavelengthbins. edge of the HII region. Nevertheless, in both cases, the peak likelihood occurs at the same value of xIGM and HI 25 Mpc <RS < 29 Mpc, RS, with only the value of fΓ being somewhat larger • ∼ ∼ here (fΓ=1.0 instead of fΓ=0.4 in Mesinger & Haiman • 0.04 <∼xIHGIM or 0.033 <∼x¯HI, 2cl0u0s4io).nWofewaatvtreilbeuntgeththbisinsshfiuftrtihnefrΓinpwriamrdarfirlyomtotthheeHinII- 0.4 <f < 1.3. region edge, where the effects of R and xIGM are felt • Γ S HI ∼ ∼ less strongly and f wields the most statistical weight. Γ In particular,the likelihoods in the two blue-most wave- 4. DISCUSSION length bins are comparable for f = 1.0 and f = 0.4 Γ Γ In this section, we discuss a few important aspects of (given xIGM =1.0 and R =41Mpc); however, the like- HI S the above results. First, it is encouraging that the sim- lihoods in the two red-most wavelength bins are 7–8 ∼ ple methodology we outlined can deliver statistically in- times greater for f =1.0 than f =0.4. Γ Γ teresting constraints, even from the few dozen available Second, our results on the neutral fraction agree in spectral pixels. Clearly, it will be interesting to apply relative terms with estimates obtained from transmis- a similar analysis to a larger sample of sources in the sion in GP troughs (Fan et al. 2006a). J1030+0524and future. J1623+3112 both have very dark GP troughs in Lyα Perhaps our most interesting result is the lower limit and Lyβ, and thus are only able to provide lower limits we obtain, for two of the quasars, on the neutral frac- onthe opticaldepth andcorrespondingneutralfraction; tion. We would like to thereforeunderstand where these however, J1148+5251 has transmission gaps in the GP constraints actually come from. To this end, in Figure 4 troughs,thus favoringasmalleropticaldepthandcorre- we show the simulated optical depth PDFs for the three spondingneutralfraction. Thisisqualitativelysimilarto wavelengthbinsusedinouranalysisofJ1623+3112. The our results: we were unable to constrain xIGM from the HI 10 spectra of J1148+5251 obtaining only a shallow peak in the Lyα line width. This does not, of course, preclude the likelihoodatxIGM =0.2;however,bothJ1030+0524 the possibility that the intrinsic Lyα lines of the z 6 HI ≈ and J1623+3112 had peak likelihoods at xIGM = 1, and quasars possess some deviations from a gaussian shape. HI strongly prefer xIGM >0.04. For example, if there is an asymmetry, in the sense of HI a steeper drop of flux on the blue side of the line than Third, the maximu∼m likelihoods of f (i.e. the rela- Γ on the red side, our technique would have mistakenly tive strength of the quasars’ ionizing luminosity) match attributed this steeper drop to absorption. This implies therelativestrengthsofthequasars’continuaredwardof that we might have overestimated the neutral fractions. Lyα. The values of f corresponding to the peak likeli- Γ At present, there is, however, no evidence of such an hoods are 1.6, 1.0, and 0.7 for J1148+5251 J1030+0524 asymmetry in the emission lines of low redshift quasars and J1623+3112 respectively. Looking at Figure 2, (where both sides of the lines can be seen unabsorbed), one can note that this is also the order of the relative nor is there evidence, from the red side of the lines, that strengthsofthequasars’continuaredwardofLyα,asone thehigh–zquasars’lineshaveshapesdifferentfromthose wouldexpectiftheUVpower-lawindicesandamountof of their lower–redshiftcounterparts. obscuration didn’t vary greatly among the quasars. Finally, we also note that our spectral fits, describing Finally, the quasar with the smallest HII region, the intrinsic emission line on the unabsorbed red side J1623+3112 with a peak likelihood at R = 29 Mpc, S of the line, leave typical flux residuals of only 1%. also has the highest favored value of xIHGIM = 1.0 and Similarresidualsarepresumablypresentontheblu∼eside the lowest favored value of f . This makes sense, since Γ of the line, as well. However, the absorption exp( τ ) a fainter quasar, embedded in a more neutral medium − D causedbytheIGMdampingwingonthebluesideofthe shouldhaveasmallerHIIregion,givenasimilarlifetime. line (see Fig. 1), has an amplitude of (0.1 1)x , HI exceeding the residuals for x > 0.01.∼Furthe−rm−ore, in HI 4.2. Uncertainties in the Intrinsic Emission Spectra orderfor the residualsto mimic∼the IGMdamping wing, Althoughwehaveaddedgaussiannoisetoourinferred theywouldhavetobecoherent(monotonicallyincreasing template spectrum, as described above, we wish to get with wavelength). a sense of how correlated errors, modifying the broader shapeofthetemplatespectrum,wouldaffectourresults. 4.3. Biases in Local Density and Velocity Fields The general concern is that mis–estimates of the flux, correlatedovermany pixels, couldpotentially mimic the As mentioned above, we remove the region 20–30 ˚A effect of shifts in our parameter space. To get a sense of bluewardofthelinecenterfromouranalysis,correspond- the magnitude of this effect, we chose to vary the width ing to the expected mean radius of the large–scale over- ofthefittedLyαgaussianforquasarJ1623+3112andsee denseregionsurroundingsuchquasars(Barkana & Loeb how this modifies our results. Naively,one wouldexpect 2004a). This is a necessary step as present-day simu- a narrower Lyα emission line and an accompanying de- lations cannot statistically model such rare, large–scale crease in effective τ to cause a shift in the parameter overdensities. However, if the density bias extended obs space to regions favoring smaller optical depths (larger into our region of analysis, our inferred xIHGIM lower lim- R , smaller xIGM, larger f ), with the opposite trends its would be conservative and our inferred fΓ values S HI Γ would most likely be underestimated. If the source is for a wider Lyα emission line. hosted by a large–scale overdensity, then density would We find that such degeneracies are very weak. For a decrease with decreasing λ (increasing distance from 10% narrower line, we find that the peak likelihood re- obs the source). If such a decreasing density extended into mains at the same parameter–combinations. The range ourregionofanalysis,itwouldtranslateintoashallower ofparametervalueswhoselikelihoodestimatesarewithin rise in τ n2/Γ. In matching the total optical depth afactorof3ofthepeaklikelihoodchangesbyatmostone R ∝ pixelinour3Dparametergrid: 25Mpc<R <29Mpc; τLyα = τR +τD, our analysis, which assumes constant S mean density, would translate this shallower rise in τ 0.008 < xIGM < 1.0; 0.4 < f < 1.6. ∼This i∼s a rather R HI Γ to a shallower rise in τ , i.e. a lower value of xIGM than small c∼hange, a∼nd we emp∼hasiz∼e that a 10% narrower appropriate if the densDity bias was taken into aHcIcount7. Lyα line would actually be an unacceptably poor fit to Note that our results on xIGM are driven by the wave- the(redsideofthe)observedspectrum. AwiderLyαline HI length dependence of τ; a constant offset in the density ismorephysicallymotivated, sincestrongdampingwing field can be absorbed by a shift in Γ (i.e. f ), which Q Γ absorptionmightbe able tosomewhatattenuateflux far dominates over Γ throughout most of our region of BG into the red side of the line. Hence, we vary the width interest and is predominantly set by the amount of ab- of the Lyα line, well in excess of the absorbing potential sorptioninbinsfurtherfromtheHIIregionedgewhereit of a strong damping wing. Again, we find little change has the most statistical weight. Furthermore, it is possi- in our results. Even for a 50% wider line, the param- ble that an infalling IGM could introduce a similar bias, eter combination at peak likelihood remains the same, buttheeffectislikelytobesmallforabrightquasarand while the range of parameter values whose likelihood es- far away from the line center on the blue side. More- timates are within a factor of 3 of the peak likelihood over,the effect from the infall would have the same sign changes by at most one pixel to: 25 Mpc < R < 31 S as the density bias above, making the flux fall off on Mpc; 0.2<xIGM <1.0;0.4<f <1.3. Ase∼xpected∼,an HI Γ intrinsical∼ly wider∼Lyα line s∼hows∼a slight preference for 7 Note that a shallower rise in τD could also be achieved if we a more neutral universe. were overestimating RS, instead of underestimating xIHGIM. How- ever,our inferredRS values arealreadyatthelowestlimitsetby Overall, the above exercise suggests that our results the onset of GP troughs for two of our quasars and thus can not are robust and conservative,at least to mis–estimates of bedecreasedfurther.