text:
# The Chandra Deep Field North Survey. X. X-ray Emission from Very Red Objects

## 1 Introduction

Understanding the properties and nature of Very Red Objects (hereafter VROs111In many studies these sources are also referred to as Extremely Red Objects (EROs).) remains a challenge to observational cosmology. This is principally because VROs are often optically faint (i.e., $`I24`$) and undetected at most wavelengths to faint limits. Although the definition of a VRO varies from study to study, it is generally accepted that a source with $`RK5.3`$ or $`IK4`$ is a VRO. These colors are consistent with those of a passively evolving elliptical galaxy at $`z\stackrel{>}{}`$1 (e.g., see Pozzetti & Mannucci 2000; Daddi et al. 2000). Optical and near-IR spectroscopy and morphology studies have indeed shown that $`z\stackrel{>}{}`$1 elliptical galaxies appear to comprise at least 40–50% of the VRO population (e.g., Cimatti et al. 1999; Moriondo et al. 2000; Mannucci et al. 2001). However, sub-millimeter and infrared observations have also demonstrated that the VRO population contains dust-enshrouded galaxies at $`z\stackrel{>}{}`$1 (e.g., Cimatti et al. 1998; Barvainis, Antonucci, & Helou 1999; Dey et al. 1999; Smail et al. 1999; Afonso et al. 2001; Elbaz et al. 2001; Pierre et al. 2001; Smith et al. 2001a). These dusty VROs may dominate the extreme red-end of the VRO population (i.e., $`IK5`$; Smail et al. 1999; Smith et al. 2001b), although the current statistics are limited. Due to the extreme faintness of VROs, the earliest studies of source densities were limited to small areas and susceptible to “cosmic variance” (e.g., see Table 1 in Scodeggio & Silva 2000). Indeed, recent wide-area VRO surveys have shown evidence for source clustering (e.g., Daddi et al. 2000, 2001).

Deep ROSAT observations first showed that VROs are detectable at X-ray energies (e.g., Newsam et al. 1997; Lehmann et al. 2000, 2001). Deeper X-ray surveys performed with the Chandra X-ray Observatory (hereafter Chandra; Weisskopf et al. 2000) have uncovered more examples of X-ray detected VROs (e.g., Alexander et al. 2001, hereafter Paper VI, Barger et al. 2001; Brandt et al. 2001b, hereafter Paper IV; Cowie et al. 2001; Crawford et al. 2001; Hornschemeier et al. 2001, hereafter Paper II; Stanford et al. 2001). In the majority of cases this X-ray emission is consistent with that expected from obscured AGN activity. In Paper II, a 221.7 ks observation of the Chandra Deep Field North (hereafter CDF-N) was used to investigate the X-ray emission from a sample of $`K_\mathrm{s}`$-band selected VROs. Four of the 33 sources were individually detected at X-ray energies although only in the hard band (2.0–8.0 keV), suggesting obscured AGN activity. Conversely, stacking analyses of the X-ray undetected VROs produced a significant average detection in only the soft band (0.5–2.0 keV). The properties of this X-ray emission were found to be consistent with those expected from normal elliptical galaxies at $`z`$ 1 (i.e., X-ray binary and thermal gas emission), suggesting a dichotomy in the X-ray properties of VROs. However, the definition of a VRO in Paper II was $`K_\mathrm{s}5.0`$, and none of the individually detected X-ray sources had $`K_\mathrm{s}5.3`$.

In this study we provide specific constraints on the X-ray properties of VROs using the stricter VRO definition of $`IK4`$ (equivalent to $`RK5.3`$) and the much deeper 1 Ms observation of the CDF-N (e.g., Brandt et al. 2001a, hereafter Paper V). We also exploit the deep and extensive multi-wavelength observations within the HDF-N and environs to provide constraints on the multi-wavelength properties of VROs. The Galactic column density along this line of sight is $`(1.6\pm 0.4)\times 10^{20}`$ cm-2 (Stark et al. 1992). We adopt the cosmology used in Paper VI (i.e., $`H_0=70`$ km s-1 Mpc-1 and $`q_0=0.1`$) throughout this paper. All coordinates in this paper are J2000.

## 2 VRO sample definitions and basic properties

In this study we have defined three VRO samples with different $`K`$-band (or equivalent) magnitude limits. The first sample includes all VROs with $`HK^{}20.4`$ detected in the $`8.4^{}\times 8.4^{}`$ area used in Paper VI (hereafter referred to as the reduced Hawaii Flanking-Field area) and is referred to as the moderate-depth VRO sample (see §2.1).222As the relationship between the $`K`$ band and $`HK^{}`$ band is $`HK^{}K=0.13+0.05(IK)`$ (Barger et al. 1999), this is equivalent to $`K<20.1`$. The second sample includes all VROs with $`K22`$ detected in the $``$ 5.3 arcmin2 region of the HDF-N itself and is referred to as the deep VRO sample (see §2.2). These two samples are complete VRO samples and allow us to estimate the fraction of X-ray emitting VROs within the VRO source population. The final VRO sample is an intermediate-depth X-ray detected VRO sample that covers the same $`8.4^{}\times 8.4^{}`$ area as the moderate-depth VRO sample but only includes X-ray detected VROs with $`HK^{}21.4`$ (see §2.3). We have defined this sample to improve the statistics of X-ray detected VROs and provide a more detailed understanding of the X-ray properties of VROs. We are unable to construct a reliable parent VRO sample at the $`2\sigma `$ magnitude limit of this sample without introducing a significant number of spurious sources. See Table 1 for the sample definitions. Thumbnail images of the 42 VROs used in this study are shown in Figure 1.

### 2.1 The moderate-depth VRO sample

As stated above, the area in which the moderate-depth VRO sample is defined is the reduced Hawaii Flanking-Field area used in Paper VI. However, to reduce the number of spurious VROs, we have further excluded sources within 17′′ of the positions of three bright $`HK^{}<14`$ stars where we could not reliably identify VROs. Although optical-to-near-IR catalogs are available for this region (Barger et al. 1999)333These images are publicly available at http://www.ifa.hawaii.edu/$``$cowie/hdflank/hdflank.html., to retain consistency with our previous studies we have used the $`I`$-band and $`HK^{}`$-band photometric catalogs (with $``$ 2$`\sigma `$ limits of $`I=25.3`$ and $`HK^{}=21.4`$, respectively) produced in Paper VI. As this region also has $`V`$-band coverage (Barger et al. 1999), we have produced a $`V`$-band catalog of sources down to $``$ 2$`\sigma `$ magnitudes of $`V=26.5`$ using the sextractor photometry tool (Bertin & Arnouts 1996) and assuming the “Best” magnitude criteria.

To ensure a low fraction of spurious VROs we use the 5$`\sigma `$ magnitude limit of $`HK^{}=20.4`$. We find 29 VROs with $`IHK^{}3.7`$ (equivalent to $`IK4.0`$) and $`HK^{}20.4`$ within this 70.3 arcmin2 region, corresponding to a source density of 1500$`{}_{250}{}^{}{}_{}{}^{+350}`$ deg-2 (see Tables 1 and 2).444All errors are taken from Tables 1 and 2 of Gehrels (1986) and correspond to the $`1\sigma `$ level; these were calculated assuming Poisson statistics. We find this VRO source density to be consistent with that found in other studies (e.g., see Table 1 of Scodeggio & Silva 2000), although the scatter from study to study is large.

### 2.2 The deep VRO sample

The deep VRO sample covers the region of the HDF-N itself. We have used the seven-band photometric catalog of Fernandez-Soto et al. (1999) to construct the deep VRO sample which has a 5$`\sigma `$ magnitude limit of $`K=22`$. The optical observations used in this catalog were taken with the Hubble Space Telescope (HST; Williams et al. 1996), and the near-IR observations were taken with ground-based telescopes (Dickinson et al. 2000). We find nine VROs with $`IK4`$ and $`K22`$ within this $``$ 5.3 arcmin2 region (see Table 3), corresponding to a source density of 6100$`{}_{2000}{}^{}{}_{}{}^{+2800}`$ deg-2 (see Table 1).555We converted the HST $`F814W`$-band photometry to the standard $`I`$ band using $`I=0.39+AB_{814}0.09(AB_{606}AB_{814})`$ (Cowie, Songaila, & Barger 1999). There are few constraints on the source density of VROs at these faint $`K`$-band magnitudes, and all studies performed to date have covered small areas (e.g., see Table 5 in McCracken et al. 2000). However, our source density is in reasonable agreement with previous studies given the probable large effects of “cosmic variance”.

Five sources within the deep VRO sample have $`K`$-band magnitudes and colors that fall within the moderate-depth VRO sample definition: two sources are included in both samples, and the other three sources fall just outside the VRO sample definition with the moderate-depth VRO sample photometry. These small differences in source magnitudes are to be expected given the faintness of these VROs. There are no sources within the moderate-depth VRO sample that lie within the HDF-N and are not included in the deep VRO sample; this is probably due to a statistical fluctuation as there does not appear to be any systematic differences in the photometry of these two samples.

### 2.3 The intermediate-depth X-ray detected VRO sample

The intermediate-depth X-ray detected VRO sample covers the same area as the moderate-depth VRO sample and includes all X-ray detected VROs down to the $``$ 2$`\sigma `$ magnitude limit of $`HK^{}=21.4`$. We are unable to define a parent VRO sample at this magnitude limit without introducing a significant number of spurious sources. However, due to the relatively small fraction of VROs with X-ray counterparts, we significantly reduce the probability of spurious $`HK^{}`$ sources by selecting those sources with X-ray detections.

The X-ray data used to produce this sample were taken from the 1 Ms CDF-N catalogs (see Paper V), which contain point sources down to on-axis flux limits of $`3\times 10^{17}`$ erg cm-2 s-1 and $`2\times 10^{16}`$ erg cm-2 s-1 in the soft and hard bands, respectively. Within the reduced Hawaii Flanking-Field area that defines this sample there are 141 X-ray sources (see Paper VI). We show the adaptively smoothed full-band (0.5–8.0 keV) Chandra image of this area in Figure 2. These X-ray sources are matched to $`I`$-band and $`HK^{}`$-band counterparts with a matching radius of 1′′ (in practice all sources had a matching distance $`<0.8`$′′) to be consistent with Paper VI; we would conservatively expect 0.2 spurious matches with this criterion. In total there are 10 X-ray detected VROs with $`IK4`$ and $`HK^{}21.4`$, corresponding to a source density of 500$`{}_{150}{}^{}{}_{}{}^{+200}`$ deg-2 (see Tables 1 and 4 and Figure 2). Six of these sources are exclusively in the intermediate-depth X-ray detected VRO sample (the other four sources are also within the moderate-depth VRO sample) and have magnitudes close to the $``$ 2$`\sigma `$ limit of $`HK^{}=21.4`$; we would expect 0.3 of these $`HK^{}`$ sources to be spurious.

Only one source is an extreme VRO (i.e., $`IK5`$), although three sources are undetected in the $`I`$ band and have lower limits on their $`IK`$ colors. The fraction of VROs within the X-ray source population at these X-ray depths is $`8_2^{+3}`$% and, with the exception of one source, all are optically faint with $`I24`$ (see Paper VI).

### 2.4 X-ray detected VROs in the moderate-depth and deep VRO samples

Four of the 29 VROs in the moderate-depth VRO sample are detected at X-ray energies, and none of the nine VROs in the deep sample is detected at X-ray energies (see Tables 1–4). All these X-ray sources were detected by wavdetect (Freeman et al. 2002) with a false-positive probability threshold of $`10^7`$ in at least one X-ray band (Paper V).

However, given the comparitively low source density of VROs, we are able to search for lower significance X-ray sources associated with VROs without introducing a significant number of spurious X-ray sources. We ran wavdetect with a false-positive probability threshold of 10-4 and found two further X-ray sources with positions within 1′′ of sources in the moderate-depth VRO sample as well as one X-ray source with a position within 1′′ of a source in the deep VRO sample (see Tables 1–4 and Figure 2). We would conservatively expect a total of 0.15 spurious X-ray sources at this threshold. The X-ray detected sources in the moderate-depth VRO sample had matching distances of $`<0.8`$′′ although the X-ray detected source in the deep VRO sample had a matching distance of 0.9′′. This larger matching distance is not unreasonable given the difficulty of source centroid determination when a low number of counts are measured. We inspected these low-significance sources to verify that they are not caused by “cosmic-ray afterglows” (Chandra X-ray Center 2000, private communication) or background fluctuations.

With these data we have detected X-ray emission from $`21_8^{+12}`$% of the $`HK^{}20.4`$ VRO population and $`11_9^{+26}`$% of the $`K22`$ VRO population (see Table 1). All but two of these X-ray detected sources are optically faint with $`I24`$. Three of the six extreme VROs (i.e., $`IK5`$) in the moderate-depth VRO sample are detected with X-ray emission. Two of these sources are extremely faint X-ray sources and have irregular optical morphologies, suggesting they may be dust-enshrouded galaxies (see §3.4.1).

### 2.5 Multi-wavelength counterparts and source redshifts

In addition to the extremely deep X-ray and moderately deep optical-to-near-IR coverage, the reduced Hawaii Flanking-Field area has the deepest feasible optical spectroscopy (e.g., Cohen et al. 2000; Dawson et al. 2001), deep 1.4 and 8.5 GHz radio observations (Richards et al. 1998; Richards 2000), and shallow sub-millimeter observations (Borys et al. 2001).666The reduced Hawaii Flanking-Field area also has additional optical coverage (e.g., Hogg et al. 2000); however, we found the Barger et al. (1999) images to provide superior resolution and depth. A fraction of this region also has deep 15 $`\mu `$m and sub-millimeter coverage (Aussel et al. 1999a; Hughes et al. 1998; Barger et al. 2000), and the central HDF-N region has some of the deepest multi-wavelength observations taken to date (see Ferguson, Dickinson, & Williams 2000). We cross-correlated our VROs to the sources in these catalogs using a matching radius of 1′′ with the exception of the ISOCAM catalog of Aussel et al. (1999a) and the SCUBA catalogs of Hughes et al. (1998) and Borys et al. (2001); see Tables 2–4. Given the larger uncertainty of the ISOCAM and SCUBA source positions, we used matching radii of 3′′ and 10′′ respectively; with these matching distances there is a non-negligible probability of mis-associations, and we would expect 0.5 spurious ISOCAM matches and 0.6 spurious SCUBA matches. As all the matches to ISOCAM sources were those included in the ISOCAM supplementary catalog and, with the exception of CXOHDFN J123647.9+621019, all had weak emission ($`f_{15\mu \mathrm{m}}<100`$ $`\mu `$Jy), there is also an additional possibility that some of the ISOCAM sources are spurious (see Aussel et al. 1999b). With the exception of the tentative (2.2$`\sigma `$) sub-millimeter counterpart of CXOHDFN J123701.6+621145 (see §3.4.1), we found no matches to SCUBA sources.

All of the sources in the deep VRO sample have photometric redshifts in the Fernandez-Soto et al. (1999) catalog of $`z=`$ 1.0–1.8 (see Table 3); this redshift range is consistent with that found in other VRO studies to date (e.g., Cimatti et al. 1999; Martini 2001). The comparitively blue near-IR colors of these sources ($`JK<2.3`$) are consistent with those expected for elliptical galaxies at $`z=`$ 1–2 (Pozzetti & Mannucci 2000), and the best fitting spectral energy distribution (SED) template for the majority of these sources is that of an elliptical galaxy.

We are unable to provide reliable photometric redshifts for the sources in the moderate-depth VRO sample and the intermediate-depth X-ray detected VRO sample due to a limited number of optical-to-near-IR bands.

## 3 Constraints on the nature of VROs

In this section we provide constraints on the nature of VROs from an analysis of their multi-wavelength properties, focusing on the X-ray emission characteristics of VROs. As most of our X-ray detected VROs do not have redshifts we have also defined an X-ray detected VRO literature sample which includes all X-ray detected VROs from the literature with spectroscopic or photometric redshifts (Fabian et al. 2000; Barger et al. 2001; Cowie et al. 2001; Crawford et al. 2001; Lehmann et al. 2001); this sample is clearly not complete and is only included to provide a comparison to the properties of our VROs.

### 3.1 Optical-to-near-IR constraints

Deep X-ray surveys have shown a correlation between the optical faintness and optical-to-near-IR color of X-ray sources (e.g., Hasinger et al. 1999; Giacconi et al. 2001; Lehmann et al. 2001). In Paper VI we showed that this correlation exists for the X-ray sources in the reduced Hawaii Flanking-Field area with $`>`$99.99% confidence. In Figure 3a we plot the mean $`IK`$ colors for X-ray sources with $`I23`$ showing this correlation. Clearly the majority of the X-ray detected VROs continue this trend, suggesting they are the extension of the X-ray source population to very red optical-to-near-IR colors. The only X-ray detected VROs that do not follow this trend are the two $`IK5.5`$ VROs which may be dust-enshrouded galaxies and different to the majority of the X-ray source population (see §3.4.1). In Figure 3b we show the mean redshifts for X-ray sources with $`I23`$ and compare them to those expected for $`M_I^{}`$ E, Sa and Sc galaxies; we determined $`M_I^{}`$ using the $`M_B^{}`$ magnitudes given in Folkes et al. (1999) and the $`BI`$ colors given in Mannucci et al. (2001). The redshifts of our X-ray undetected VROs are $`z`$ 1–2 (see Table 3), consistent with the few known spectroscopic and photometric redshifts of VROs to date (e.g., Cimatti et al. 1999; Martini 2001; Willott, Rawlings &

Blundell 2001); on average these sources appear to reside in approximately $`M_I^{}`$ galaxies. The redshifts of the VROs in the X-ray detected VRO literature sample cover $`z=`$ 1–3, similar to that found for X-ray undetected VROs.

We can also use the $`V`$-band observations to place further constraints on the properties of X-ray detected VROs. In Figure 3c we show the $`VI`$ versus $`IK`$ colors of all the VROs and compare them to those expected for different galaxy types.777We converted the HST $`F606W`$-band photometry to the standard $`V`$ band using $`V=0.16+AB_{606}+0.40(AB_{606}AB_{814})`$ (Cowie et al. 1999). All of the VROs in the deep VRO sample have $`V`$-band detections, and almost all lie within the region between the Sa and E galaxy tracks.888We have applied $`K`$-corrections to these galaxy tracks using the $`K`$-correction factors given in Poggianti (1997). However, we have not applied evolutionary corrections; hence these tracks should be considered extreme in the sense that a typical $`z=1`$ elliptical should have bluer colors than those shown here. Since our analysis is not sensitive to the accuracy of these tracks, our overall conclusions will not change even if the true evolutionary corrections are found to be large. The only X-ray detected source in the deep VRO sample appears to be an elliptical galaxy in the HST images and lies between the predicted tracks for $`z`$ 1 Sa and E galaxies (see §3.2.2 and Figures 3c and 3d). Due to the shallower $`V`$-band coverage of the moderate-depth VRO sample, only 15 of the 29 sources are detected, restricting the constraints we can place on the host galaxy types of these VROs. With the exception of the two $`IK>5.5`$ VROs, which are very faint in the optical band and may be dust-enshrouded galaxies (see §3.4.1), all of the X-ray detected VROs are also detected in the $`V`$-band; however, the average $`VI`$ colors of these sources are consistent with those found for the X-ray undetected VROs with $`V`$-band detections. Unfortunately the depth and resolution of these images are not sufficient to determine accurately the host-galaxy types of these sources.

### 3.2 X-ray constraints for X-ray detected VROs

Key distinctions between AGN and “normal” galaxies can be made at X-ray energies. AGN are generally distinguishable by having moderate-to-high X-ray luminosities, high X-ray-to-optical flux ratios, and, in the case of obscured AGN, flat X-ray spectral slopes.999The flat X-ray spectral slopes of obscured AGN are due to absorption and are not intrinsic. In Figure 4a we show the hard-band flux versus the X-ray band ratio, defined as the ratio of hard-band to soft-band counts, for all the X-ray detected VROs. The range of spectral slopes for the nine VROs detected in the hard band is $`0.3<\mathrm{\Gamma }<1.5`$, and the origin of the X-ray emission from these sources is almost certainly obscured AGN activity (see §3.2.1). There are, however, four VROs undetected in the hard band; we argue below that at least some of these sources are not obscured AGN (see §3.2.2).

#### 3.2.1 Hard-band detected VROs

The constraints on the X-ray spectral slope for the nine hard-band detected VROs are poor in most cases. However, we can derive an average X-ray spectral slope by stacking the individual sources using the stacking technique presented in §3.3 of Paper VI. We find a mean band ratio of $`0.94\pm 0.07`$, corresponding to a mean spectral slope of $`\mathrm{\Gamma }`$ 0.9. This is very similar to that found for the optically faint X-ray source population (Paper VI), suggesting that these X-ray detected VROs are simply the extension of the optically faint X-ray source population to very red optical-to-near-IR colors. Such a flat X-ray spectral slope is strongly indicative of obscured AGN activity and, assuming the observed X-ray spectral slope is due to obscuration of $`\mathrm{\Gamma }=2.0`$ power-law emission, corresponds to an intrinsic absorption column density at $`z=`$ 1–3 (the probable redshift range of these sources; see §3.1) of $`N_\mathrm{H}`$ (5–35)$`\times 10^{22}`$ cm-2. Three of the VROs are undetected in the soft band and have inverted spectra ($`\mathrm{\Gamma }<0.0`$); these characteristics can be indicative of Compton-thick absorption. We can derive an average spectral slope for these sources by stacking the individual sources in the same manner as performed above; this gives a mean band ratio of $`4.2_{1.0}^{+1.3}`$, corresponding to a mean spectral slope of $`\mathrm{\Gamma }`$ –0.4 and an average intrinsic column-density at $`z=`$ 1–3 of $`N_\mathrm{H}`$ (0.2–1.2)$`\times 10^{24}`$ cm-2. The upper-limit absorption estimate is just below that expected for Compton-thick absorption ($`1.5\times 10^{24}`$ cm-2). However, we note that this column density range should be considered a lower limit since additional components other than the intrinsic AGN emission (e.g., scattered AGN emission and star-formation) can raise the soft-band flux.

Another key signature of AGN activity is luminous X-ray emission. In determining source luminosities we have corrected for apparent intrinsic absorption in each VRO by assuming that the observed X-ray spectral slope is due to obscuration of $`\mathrm{\Gamma }=2.0`$ power-law emission; the equivalent $`z`$ 0 absorbing column densities for a given band ratio are given in Figure 4a. For those sources that do not have a measured redshift we have assumed $`z=1`$ and $`z=3`$. The estimated rest-frame unabsorbed 0.5–8.0 keV luminosities of all the X-ray detected VROs are given in Table 4. The range in X-ray luminosity for the hard-band detected VROs is $``$ 3$`\times 10^{42}`$ erg s-1 to $``$ 5$`\times 10^{44}`$ erg s-1, and two sources (CXOHDFN J123615.8+621515 and CXOHDFN J123714.0+620916) would be obscured QSOs (i.e., $`>3\times 10^{44}`$ erg s-1) at $`z`$ 3; see also Paper VI. As star-formation activity is rarely more luminous than $`10^{42}`$ erg s-1 in the X-ray band (e.g., Moran, Lehnert, & Helfand 1999 and references therein), it is probable that all the hard-band detected VROs are obscured AGN; these results are consistent with those found in Paper II. If we assume that we have detected all luminous AGN sources in the moderate-depth VRO sample then the fraction of such sources within the $`K20.1`$ VRO population is $`14_7^{+11}`$%.

Further evidence that these sources host AGN activity is found from their X-ray-to-optical flux ratios (see Figure 4b). With the exception of one source, all of the VROs detected in the hard band lie within the boundaries expected for typical AGN (i.e., $`1<\mathrm{log}(\frac{f_\mathrm{X}}{f_\mathrm{I}})<1`$; Akiyama et al. 2000; Paper VI). The only source found to lie outside these boundaries is CXOHDFN J123615.8+621515. This source is optically blank ($`I>25.3`$) and has $`\mathrm{log}(\frac{f_\mathrm{X}}{f_\mathrm{I}})>1`$, an inverted X-ray spectral slope ($`\mathrm{\Gamma }<`$ –0.3) and a $`z=3`$ rest-frame unabsorbed 0.5–8.0 keV luminosity of $``$ 4$`\times 10^{44}`$ erg s-1. This source may be an obscured QSO; see §6.2 in Paper VI for further discussion of this source.

#### 3.2.2 Hard-band undetected VROs

The four X-ray detected VROs that are undetected in the hard band are also faint in the soft band, limiting the X-ray spectral constraints we can provide. We can, however, derive an average X-ray spectral slope by stacking the individual sources using the technique given in §3.2.1. This still does not provide an average detection in the hard band, although it places a constraint on the average band ratio of $`<0.56`$, corresponding to $`\mathrm{\Gamma }>1.4`$. This X-ray spectral slope is considerably steeper than that found for the hard-band detected VROs, suggesting the X-ray emission mechanisms of these sources are different. Since we have redshift determinations for two of these VROs (CXOHDFN J123651.9+621225 and CXOHDFN J123701.6+621145) we can provide better constraints on their X-ray emission mechanism.

CXOHDFN J123651.9+621225 (optical identification 123652.0+621226) lies within the HDF-N itself and has a photometric redshift in the Fernández-Soto et al. (1999) catalog of $`z=1.0`$; three independent photometric redshift studies also place this galaxy at $`z=`$ 1.0–1.2 (see Hughes et al. 1998). We made our own determination using the publicly available photometric redshift code hyperz Version 1.1 (Bolzonella, Miralles, & Pelló 2000).101010The hyperz code is available at http://webast.ast.obs-mip.fr/hyperz/. In performing the photometric redshift fitting, we used only the Bruzual & Charlot (1993) spectral templates and allowed up to 1.6 mags of visual extinction. The best-fitting model was found to be a $`z=1.00_{0.06}^{+0.02}`$ elliptical galaxy with an age of 6.5 Gyr and 0.3 mags of visual extinction; the uncertainty in the redshift determination corresponds to the 90% confidence level. Indeed, based on optical morphology, this source is clearly an elliptical galaxy and has $`VI`$ versus $`IK`$ colors consistent with that of an elliptical at $`z`$ 1 (see Figures 3c and 3d). Assuming $`\mathrm{\Gamma }=2.0`$, the rest-frame 0.5–2.0 keV luminosity is $``$ 9$`\times 10^{40}`$ erg s-1, and the rest-frame 0.5–8.0 keV luminosity is $`<3\times 10^{41}`$ erg s-1. Furthermore, with a $`K`$-corrected $`B`$-band magnitude of $`M_B=20.7`$, the X-ray emission from this source is consistent with that found for local elliptical galaxies (e.g., O’Sullivan, Forbes, & Ponman 2001; S. Pellegrini 2001, private communication).

CXOHDFN J123701.6+621145 (optical identification 123701.5+621146) is an optically faint $`\mu `$Jy radio source (Richards et al. 1999) and has a spectroscopic redshift of $`z=0.884`$, based on a single emission line assumed to be \[OII\] $`\lambda `$3727 (Cohen et al. 2000).111111The SED of CXOHDFN J123701.6+621145 provides corroborating evidence for this redshift (see §3.4.1). Adopting $`\mathrm{\Gamma }=2.0`$, the rest-frame 0.5–2.0 keV luminosity is $``$ 1$`\times 10^{41}`$ erg s-1, and the rest-frame 0.5–8.0 keV luminosity is $``$ 4$`\times 10^{41}`$ erg s-1. As the irregular morphology of this source (van den Bergh et al. 2000) clearly indicates that it is not an elliptical galaxy, the X-ray emission is probably produced by either luminous star-formation or low-luminosity AGN activity; we provide further constraints in §3.4.1.

We do not have detailed constraints for the other two hard-band undetected VROs, CXOHDFN J123639.6+620936 and CXOHDFN J123642.8+620958. However, their irregular $`HK^{}`$-band morphologies and extremely faint optical magnitudes ($`I\stackrel{>}{}25`$) suggest they may be similar to CXOHDFN J123701.6+621145 and not normal elliptical galaxies. The extreme optical-to-near-IR color of CXOHDFN J123642.8+620958 ($`IK=6.2`$) and the lower limit optical-to-near-IR color of CXOHDFN J123639.6+620936 ($`IK>4.5`$) provide further support for this interpretation.

If we assume the hard-band undetected sources in the moderate-depth VRO sample with irregular morphologies are luminous star-forming galaxies, the fraction of such sources within the $`K20.1`$ VRO population is $`7_5^{+8}`$%. Clearly, this fraction should be used with caution since we have not proved that these sources are luminous star-forming galaxies, and a large number of such sources may lie below the detection limit of our X-ray observations.

### 3.3 X-ray constraints for X-ray undetected VROs

The majority of our VROs do not have detectable X-ray emission. We can, however, search for evidence of X-ray emission from these sources on average by stacking the individually undetected sources in the same manner as performed in previous studies (e.g., Paper II; Paper IV; Paper VI). We conservatively excluded from the stacking analyses five of the 23 X-ray undetected VROs in the moderate-depth VRO sample and one of the eight X-ray undetected VROs in the deep VRO sample that lay within 10′′ of a high-significance X-ray source (i.e., detected with a wavdetect false-positive probability threshold of $`10^7`$) or within 3′′ of a low-significance X-ray source (i.e., detected with a wavdetect threshold of $`10^4`$ but undetected with a wavdetect threshold of $`10^7`$); these could be affected by the low-level wings of the Chandra point-spread function (see Tables 2 and 3). In all the stacking analyses we used the 30-pixel aperture and restricted ACIS grade set defined in Paper IV. To determine the significance of the stacked X-ray emission we performed Monte-Carlo stacking simulations using the technique given in Appendix A of Paper IV. The Monte-Carlo simulations consisted of 10,000 trials where we stacked 18 random positions (seven random positions in the case of the deep VRO sample) within each VRO sample area using the same photometric aperture as for the stacking analyses. We verified the significances determined from these simulations using Poisson statistics; all significant detections (i.e., $`99`$% confidence level) were found to be so with both methods. The results of the stacking analyses are given in Table 5; we have assumed $`\mathrm{\Gamma }=2.0`$ when determining all fluxes and luminosities.

#### 3.3.1 Results from the moderate-depth VRO sample

The effective Chandra exposure time achieved from the stacking analyses of the moderate-depth VRO sample is 15.8 Ms. We find significant detections ($``$99% confidence level) in the soft and full bands corresponding to average fluxes of $`6.0\times 10^{18}`$ erg cm-2 s-1 and $`1.5\times 10^{17}`$ erg cm-2 s-1, respectively.

We may expect the optically brightest sources to contribute more to the stacked X-ray emission as for a given X-ray-to-optical flux ratio these sources would have brighter X-ray fluxes. To investigate this possibility we performed separate stacking analyses for the $`I24.1`$ and $`I>24.1`$ VROs; the $`I=24.1`$ threshold was chosen to give equal numbers of VROs in both sub-samples. The effective Chandra exposure times from these stacking analyses are similar (8.1 Ms for the $`I24.1`$ sub-sample and 7.7 Ms for the $`I>24.1`$ sub-sample) and allow for a relatively unbiased comparison. We find significant detections in the soft and full bands and a tentative (96.4% confidence level) detection in the hard band for the $`I24.1`$ VROs; there are no significant detections for the $`I>24.1`$ VROs in any X-ray band (see Table 5). This clearly shows that the $`I24.1`$ VROs are contributing the majority of the stacked X-ray emission, and deeper Chandra observations should preferentially detect the optically brightest VROs. Given the average soft-band flux ($`8.7\times 10^{18}`$ erg cm-2 s-1) and mean $`I`$-band magnitude ($`I=23.7\pm 0.3`$) of the $`I24.1`$ VROs, the mean X-ray-to-optical flux ratio of $`\mathrm{log}(\frac{f_\mathrm{X}}{f_\mathrm{I}})=1.9`$ is consistent with that expected from normal elliptical galaxies. Furthermore, given the mean $`I`$-band magnitude for the $`I>24.1`$ sources ($`I=24.8\pm 0.6`$) and assuming the same X-ray-to-optical flux ratio, the expected soft-band flux of $`3.3\times 10^{18}`$ erg cm-2 s-1 is consistent with our non-detection of the $`I>24.1`$ VROs.

The fainter $`I`$-band magnitudes of the $`I>24.1`$ VROs could be due to these sources lying at higher redshifts than the $`I24.1`$ VROs. Although we do not have redshifts for the majority of these VROs, assuming they are $`M_I^{}`$ elliptical galaxies the average redshifts of the $`I24.1`$ and $`I>24.1`$ sources are $`z=1.1`$ and $`z=1.4`$, respectively. At $`z=1.1`$, the average rest-frame 0.5–2.0 keV and 0.5–8.0 keV luminosities of the $`I24.1`$ VROs are $`5.0\times 10^{40}`$ erg s-1 and $`1.2\times 10^{41}`$ erg s-1, respectively. This X-ray emission is at the level expected from normal elliptical galaxies, although it could also be produced by star formation and low-luminosity AGN activity. At $`z=1.4`$, the average rest-frame 0.5–2.0 keV and 0.5–8.0 keV luminosity upper limits of the $`I>24.1`$ VROs are $`<6.2\times 10^{40}`$ erg s-1 and $`<2.1\times 10^{41}`$ erg s-1, respectively. Both of these X-ray luminosity upper limits are higher than those determined for the $`I24.1`$ VROs and therefore the non-detection of X-ray emission from the $`I>24.1`$ VROs is consistent with these sources lying at higher redshifts. Deep imaging in multiple optical-to-near-IR bands is required to determine source morphologies and photometric redshifts to test these hypotheses (see §4.2).

#### 3.3.2 Results from the deep VRO sample

The effective Chandra exposure time achieved from the stacking analyses of the deep VRO sample is 6.5 Ms. We find a significant detection ($``$99% confidence level) in the soft band corresponding to an average flux of $`6.5\times 10^{18}`$ erg cm-2 s-1. At the mean photometric redshift of these sources ($`z=1.5\pm 0.2`$; see Table 3), the rest-frame 0.5–2.0 keV luminosity is $`9\times 10^{40}`$ erg s-1, the same as that found for the one X-ray detected VRO in the deep VRO sample (CXOHDFN J123651.9+621225 at $`z=1.0`$; see §3.2.2). As CXOHDFN J123651.9+621225 is the lowest redshift source in the deep VRO sample, this suggests that the individual non-detections of the other VROs are mostly due to these sources lying at higher redshifts. This result is entirely consistent with that suggested in §3.3.1 based on the stacking analyses of X-ray undetected sources in the moderate-depth VRO sample, suggesting that deeper Chandra observations should detect more examples of X-ray emitting normal elliptical galaxies in the VRO population.

As only 7 VROs were used in these stacking analyses we are unable to investigate distinctions in average X-ray properties with optical brightness or redshift.

### 3.4 Multi-wavelength constraints

#### 3.4.1 Broad-band SED of CXOHDFN J123701.6+621145

Most of our VROs are only detected at optical-to-near-IR wavelengths; however, CXOHDFN J123701.6+621145 is remarkable in that it is detected at most of the wavelengths investigated here. This source is an extreme VRO (i.e., $`IK5`$), has an irregular morphology (van den Bergh et al. 2000), and has a spectroscopic redshift of $`z=0.884`$ (Cohen et al. 2000); see also §3.2.2. The extreme VRO colors and irregular morphology suggest this source may be a member of the dust-enshrouded VRO population.

In Figure 5 we compare the SED of CXOHDFN J123701.6+621145 to that of the archetypal nearby dusty starburst galaxy Arp 220. The fluxes of Arp 220 have been adjusted to $`z=0.884`$ to show the relative luminosities of both sources from radio to X-ray wavelengths; however, to produce the fit shown, we had to reduce the luminosity of the Arp 220 model by a factor of 2.2. Clearly CXOHDFN J123701.6+621145 has a similar SED to Arp 220. The bolometric luminosity of CXOHDFN J123701.6+621145 is lower than that of Arp 220, although the X-ray luminosity is $`2.5`$ times higher (see §3.2.2). The X-ray emission from Arp 220 is consistent with that produced by both a star-formation component and a hard X-ray power-law component (Iwasawa et al. 2001); this latter component could be due to AGN emission or a population of X-ray binaries. We do not have sufficient X-ray information to determine the origin of the X-ray emission from CXOHDFN J123701.6+621145; although we note that it is less luminous than that found for NGC 3256 (Moran et al. 1999; Lira et al. 2001), the most X-ray luminous local starburst.

The similarity of the $`IK`$ colors and X-ray properties of CXOHDFN J123639.6+620936 and CXOHDFN J123642.8+620958 (see §3.2.2) with those of CXOHDFN J123701.6+621145 suggests they may also be dust-enshrouded galaxies.

#### 3.4.2 Radio detections of VROs

Five VROs have radio counterparts, four of which are also detected at X-ray energies. This high matching fraction between radio and X-ray sources is consistent with that found in other studies (i.e., Paper II; Paper IV; Bauer et al. 2001a) suggesting an association between the production of X-ray and radio emission. Only three of the five VROs with radio counterparts have well-defined radio spectral slopes; two sources (123641.4+621143 and CXOHDFN J123707.2+621408) have flat radio spectra (i.e., $`\alpha _R0.3`$) while a third source (CXOHDFN J123701.6+621145) has a steep radio spectrum ($`\alpha _R=0.82_{0.10}^{+0.09}`$) and is described in §3.4.1. The origin of the radio emission in the sources with flat radio spectra can be either AGN or recent starburst activity (e.g., Condon 1992; Richards 2000). The flat X-ray spectrum of CXOHDFN J123707.2+621408 favors the former case for this object. The source 123641.4+621143 is undetected in all the X-ray bands; given a full-band upper limit of $`2\times 10^{16}`$ erg cm-2 s-1, its rest-frame 0.5–8.0 keV luminosity can be as high as $`3\times 10^{42}`$ erg s-1 at $`z=1.524`$ (Dawson et al. 2001). As this source has a disk-like morphology, it is most likely a star-forming galaxy or a low-luminosity AGN.

#### 3.4.3 ISOCAM detections of VROs

Four of the 14 VROs in the ISOCAM region of Aussel et al. (1999a) have 15 $`\mu `$m counterparts, two of which are also detected at X-ray energies (see Tables 2 and 4). This infrared radiation can be due to either dust or starlight emission. In Figure 5, we show the SEDs of both a dusty starburst (Arp 220) and a starlight-dominated galaxy (a giant elliptical galaxy) normalized to the $`HK^{}`$-band emission of CXOHDFN J123701.6+621145. The predicted 15 $`\mu `$m emission from starlight is an order of magnitude below that expected from the dusty starburst, and clearly the origin of the infrared radiation in CXOHDFN J123701.6+621145 is dust emission. As CXOHDFN J123701.6+621145 is the brightest ISOCAM-detected VRO in the $`HK^{}`$-band, this suggests that starlight emission at 15 $`\mu `$m for $`HK^{}>19.5`$ VROs at $`z`$ 1 is $`<1\mu `$Jy. This is 1–2 orders of magnitude fainter than the 15 $`\mu `$m flux limit, and the origin of the infrared emission in all the ISOCAM-detected VROs is likely to be dust emission.

### 3.5 VRO clustering

Wide-area VRO studies have shown that VROs are strongly clustered (e.g., Daddi et al. 2000, 2001; McCarthy et al. 2001). The solid angles covered by our VRO samples are too small to perform a detailed analysis of source clustering; the results of a two-point correlation analysis on the moderate-depth VRO sample with the estimator of Landy & Szalay (1993) did not show significant evidence for source clustering.

We also took an alternative approach and searched for VRO overdensities in the vicinities of three known moderate-redshift (i.e., $`z`$ 1) clusters and cluster candidates (see Figure 2). Three VROs lie within the $``$ 45′′ radius region of the cluster-candidate X-ray source CXOHDFN J123620.0+621554 (Bauer et al. 2001b, hereafter Paper IX) when we would expect 0.9 VROs in the case of a uniform distribution. This unusual X-ray source also contains an overdensity of optically faint X-ray sources and optically faint $`\mu `$Jy radio sources, and it is likely to be a poor cluster at $`z\stackrel{>}{}0.7`$; see Paper IX for further discussion. We found one VRO within 45′′ of the center of the $`z=0.85`$ cluster ClG 1236+6215 discovered by Dawson et al. (2001) and two VROs (excluding the VROs in the deep VRO sample) within 45′′ of the center of the $`z=1.01`$ cluster/group associated with the FR I radio galaxy VLA J123644.3+621133 (e.g., Richards et al. 1998; Paper IX). The evidence for VRO clustering within the vicinities of all these sources is not significant although we are limited by statistics. If we combine these results we find six VROs when we would expect 2.7 in the case of a uniform distribution; the evidence for VRO clustering is significant at the 94% confidence level.

## 4 Summary and Further work

### 4.1 Summary of main results

We have used the 1 Ms CDF-N observation of an $`8.4^{}\times 8.4^{}`$ region within the Hawaii Flanking-Field area to provide constraints on the X-ray emission from Very Red Objects (VROs; $`IK4`$). We defined three VRO samples: a moderate-depth VRO sample ($`HK^{}20.4`$) and an intermediate-depth X-ray detected VRO sample ($`HK^{}21.4`$) both covering the 70.3 arcmin2 reduced Hawaii Flanking-Field area, and a deep VRO sample ($`K22`$) covering the 5.3 arcmin2 region of the HDF-N itself. Our main results are the following:

(i) We detected X-ray emission from 13 VROs down to 0.5–2.0 keV and 2.0–8.0 keV point-source fluxes of $`3\times 10^{17}`$ erg cm-2 s-1 and $`2\times 10^{16}`$ erg cm-2 s-1, respectively. Six ($`21_8^{+12}`$%) of the 29 VROs within the moderate-depth VRO sample and one ($`11_9^{+26}`$%) of the nine VROs within the deep VRO sample were detected with X-ray emission. See §2.3 and §2.4.

(ii) The $`IK`$ color versus $`I`$-band magnitude of all the $`IK<5.5`$ X-ray detected VROs follows the trend found for optically brighter ($`I23`$) X-ray sources, suggesting the VROs are the $`z1`$ very red analogs of the optically bright X-ray source population. The two $`IK>5.5`$ X-ray detected VROs have redder $`IK`$ colors than those expected given their $`I`$-band magnitudes and may be dust-enshrouded VROs. See §3.1 and §3.4.1.

(iii) The hard-band detected VROs have flat X-ray spectral slopes ($`\mathrm{\Gamma }0.9`$) and X-ray properties consistent with those expected from obscured AGN; three of these VROs are undetected in the soft band, and the absorption to the X-ray nucleus in some of these sources may be Compton-thick. Assuming all the hard-band detected VROs lie at $`z=`$ 1–3, they have rest-frame 0.5–8.0 keV unabsorbed luminosities of $``$ $`3\times 10^{42}`$ to $`5\times 10^{44}`$ erg s-1; the two brightest sources would be obscured QSOs (i.e., $`>3\times 10^{44}`$ erg s-1) if they lie at $`z`$ 3. Assuming these hard-band detected VROs comprise all the luminous AGNs within the moderate-depth VRO sample, the fraction of such sources in the $`HK^{}20.4`$ VRO population is $`14_7^{+11}`$%. Unfortunately, the depth and resolution of the optical-to-near-IR images are not sufficient to determine the host-galaxy types of these sources. See §3.2.1.

(iv) The hard-band undetected VROs have steeper X-ray spectral slopes and are more likely to contain less energetic X-ray emission processes such as star formation, low-luminosity AGN, and normal elliptical galaxy emission (i.e., X-ray binary activity and thermal gas emission). One of these sources appears to be a normal elliptical galaxy at $`z`$ 1, while another source is a $`z=0.884`$ $`IK>5.5`$ VRO that has a SED very similar to that of the archetypal dusty starburst galaxy Arp 220. The X-ray emission from this source is likely to be due to either luminous star-formation or low-luminosity AGN activity. The irregular morphologies and optical-to-near-IR colors of the other hard-band undetected VROs are similar to those of the $`z=0.884`$ extreme VRO, and the X-ray emission from these sources is also likely to be produced by the same processes. Assuming the hard-band undetected sources with irregular morphologies in the moderate-depth VRO sample are luminous star-forming galaxies, the fraction of such sources within the $`K20.1`$ VRO population is $`7_5^{+8}`$% or larger. See §3.2.2 and §3.4.1.

(v) Stacking analyses of VROs not individually detected with X-ray emission produced significant detections ($`99`$% confidence level) in the 0.5–2.0 keV band for both the moderate-depth and deep VRO samples, and a significant detection in the 0.5–8.0 keV band for the moderate-depth VRO sample. These average X-ray detections are consistent with the emission expected from normal elliptical galaxies, star formation or low-luminosity AGN activity at $`z`$ 1. We found that the optically brightest VROs contribute the majority of this X-ray emission. The simplest explanation of these results is that we have detected the X-ray emission from non-active VROs with low X-ray-to-optical flux ratios \[$`\mathrm{log}(\frac{f_\mathrm{X}}{f_\mathrm{I}})`$ –2\]. This is consistent with expectations if the majority of these VROs are $`M_I^{}`$ elliptical galaxies. See §3.3.

(vi) Five VROs have radio counterparts, and four of these are detected with X-ray emission. This high matching rate between radio and X-ray sources is consistent with that found from other deep X-ray and radio studies. See §3.4.2.

(vii) Four of the 14 VROs within the ISOCAM HDF-N region have 15 $`\mu `$m counterparts, and two of these are detected with X-ray emission. The origin of the infrared radiation from all these ISOCAM sources is likely to be dust emission. See §3.4.3.

### 4.2 Further work

Given our results we can provide suggestions for further work to extend the understanding of X-ray emission from VROs. Arguably the most important issue to address is identifying the differences between X-ray detected and X-ray undetected VROs. There are three main approaches that may provide useful constraints in this regard: (i) galactic host types and environments, (ii) redshift distributions, and (iii) large-scale clustering properties.

The determination of galactic host types and environments will be most effectively pursued with high-resolution deep (i.e., $`I`$ 27) optical/near-IR imaging; the excellent angular resolution of HST is preferable in this regard. As the photometric redshift technique will probably be the only viable method of determining the redshifts for most VROs, deep imaging in multiple optical-to-near-IR bands will be necessary to determine redshift distributions. However, given that many of the X-ray detected VROs are probably AGN, deep optical spectroscopy may successfully detect high equivalent-width emission lines in some sources. It is becoming clear that VROs are strongly clustered although the clustering properties of X-ray detected VROs are almost completely unknown. The advantage of an investigation into the clustering properties of X-ray detected VROs is that it can be undertaken with comparatively moderate depth optical-to-near-IR imaging (i.e., $`I`$ 25). However, to place useful constraints and construct statistically significant samples, imaging observations over a large area (i.e., a few hundred arcmin2) would be required.

The 1 Ms Chandra observation used in the present study has shown that VROs without obvious AGN activity are probably only detected at faint X-ray fluxes (i.e., 0.5–2.0 keV fluxes $`<10^{16}`$ erg cm-2 s-1) and are therefore only detectable in extremely deep X-ray surveys. Our ultimate aim with the CDF-N survey is to reach a Chandra exposure time of $`5`$ Ms. Such an observation should reach 0.5–2.0 keV fluxes of $`6\times 10^{18}`$ erg cm-2 s-1 and will be deep enough to detect a large number of VROs without obvious AGN activity.

## Acknowledgments

This work would not have been possible without the support of the entire Chandra and ACIS teams; we particularly thank P. Broos and L. Townsley for data analysis software and CTI correction support. We thank L. Pozzetti for interesting discussions, L. Silva and G.L. Granato for providing model SEDs, and S. Dawson, M. Dickinson, S. Pellegrini, and P. Severgnini for providing useful information. We are grateful to A. Barger and L. Cowie for making their optical and near-IR images publicly available and to M. Bolzonella, J.-M. Miralles and R. Pelló for making hyperz available. We thank the referee for concise comments. We acknowledge the financial support of NASA grant NAS 8-38252 (GPG, PI), NSF CAREER award AST-9983783 (DMA, CV, FEB, WNB), NASA GSRP grant NGT5-50247 and the Pennsylvania Space Grant Consortium (AEH), and NSF grant AST-9900703 (DPS). C.V. also acknowledges partial support from the Italian Space Agency, under the contract ASI 00/IR/103/AS, and the Ministry for University and Research (MURST) under grant Cofin-00-02-36.


format: text

text:
# How Neutral is the Intergalactic Medium at 𝑧∼6?

## 1 Introduction

Recent spectroscopic observations of $`z>\mathrm{\hspace{0.33em}4.5}`$ quasars discovered by the Sloan Digital Sky Survey (SDSS) have opened up new windows into the study of the high redshift intergalactic medium (IGM) (Fan et al. 2000, Zheng et al. 2000, Schneider et al. 2001, Anderson et al. 2001 Fan et al. 2001a, Becker et al. 2001, Djorgovski et al. 2001). In particular, Becker et al. (2001) observed Gunn-Peterson troughs (Gunn & Peterson 1965) in the spectrum of a $`z=6.28`$ quasar, which were interpreted as suggesting that the universe was close to the reionization epoch at $`z6`$.

That the absorption increases quickly with redshift is not by itself surprising: ionization equilibrium tells us that the neutral hydrogen density is proportional to the gas density squared, which is proportional to $`(1+z)^6`$ at the cosmic mean. The evolution of the ionizing background and gas temperature will modify this redshift dependence, but the rapid evolution of absorption remains a robust outcome. What is interesting, as Becker et al. (2001) emphasized, is that the observed mean transmission at redshift $`z6`$ is lower than what one would expect based on an extrapolation of the column density distribution and its redshift evolution (number density of clouds scaling as $`(1+z)^{2.5}`$) from lower redshifts. On the other hand, the popular gravitational instability theory of structure formation provides detailed predictions for how the IGM should be clustered, and how this clustering evolves with redshift, which has been shown to be quite successful when compared with observations at $`z24`$ (See e.g. Cen et al. 1994, Zhang et al. 1995, Reisenegger & Miralda-Escudé 1995, Hernquist et al. 1996, Miralda-Escudé 1996, Muecket et al. 1996, Bi & Davidsen 1997, Bond & Wadsley 1997, Hui et al. 1997, Croft et al. 1998, Theuns et al. 1999, Bryan et al. 1999, McDonald et al. 2000a). These predictions allow us to directly infer the neutral fraction of the IGM from the observed absorption (the relation between the two depends on the nature of clustering of the IGM), and so can further inform our interpretations of the recent $`z6`$ results.

How neutral is the IGM at $`z6`$, and how different is the neutral fraction compared to lower redshifts? These are the questions we would like to address quantitatively, making use of the gravitational instability model of the IGM.

The paper is organized as follows. First, we start with a brief description of the gravitational instability model for the IGM and the simulation technique in §2. In §3.1, we derive the neutral hydrogen fraction $`X_{\mathrm{HI}}`$, and equivalently the level of ionizing flux $`J_{\mathrm{HI}}`$, at several different redshifts leading up to $`z6`$ from the observed mean Lyman-alpha (Ly$`\alpha `$) transmission. This exercise using the Ly$`\alpha `$ spectrum is similar to the one carried out in McDonald & Miralda-Escudé (2001), except for the addition of new high redshift data. 111This part of the calculation involving the matching of the mean Ly$`\alpha `$ transmission is also similar to a number of earlier papers where the primary quantity of interest is the baryon density (e.g. Rauch et al. 1997, Weinberg et al. 1997, Choudhury et al. 2000, Hui et al. 2001). Here, we fix the baryon density and study the ionizing background or the neutral fraction instead (see §2). We then examine in §3.2 the constraints on the same quantities $`X_{\mathrm{HI}}`$ and $`J_{\mathrm{HI}}`$ from the observed mean Lyman-beta (Ly$`\beta `$) transmission, Ly$`\beta `$ being particularly useful at high Ly$`\alpha `$ optical depth, because the Ly$`\beta `$ absorption cross-section is a factor of $`5`$ smaller than the Ly$`\alpha `$ cross-section. The goal here is to use Ly$`\beta `$ absorption to obtain constraints on $`X_{\mathrm{HI}}`$ and $`J_{\mathrm{HI}}`$ that are as stringent as possible. In §3.2, we also examine the sensitivity of our conclusions to the power spectrum normalization.

An intriguing question is: instead of focusing on the mean transmission, can one make use of the fact that the observed spectrum at $`z6`$ contains a continuous and long stretch ($`200300\AA `$) with little or no detected flux to obtain more stringent limits on the neutral fraction or $`J_{\mathrm{HI}}`$? The idea is that since the IGM gas density naturally fluctuates spatially, it seems a priori unlikely to have no significant upward fluctuation in transmission for many pixels in a row – unless of course the neutral fraction $`X_{\mathrm{HI}}`$ is indeed quite high. We will show in §3.3 this provides constraints that are slightly weaker to those obtained using the mean transmission.

In all the simulations discussed in this paper, the ionizing background is assumed uniform spatially, just as in the majority of high redshift IGM simulations. A natural worry is that as the universe becomes more neutral at higher redshifts, the ionizing background would be more non-uniform. One way to test this is to use several lines of sight, available at $`z5.5`$, and compare the observed line-of-sight scatter in mean transmission against the predicted scatter based on simulations with a uniform background. We discuss this in §4, estimate the level of ionizing background fluctuations, and make predictions for the scatter at $`z6`$. Here, we also introduce a technique to handle the problem of limited box-size.

Readers who are not interested in details can skip to §5 where we summarize the constraints obtained. We also discuss in §5 the issue of continuum placement, and how the associated uncertainties can be estimated. While the work described in this paper was being carried out, several papers appeared which investigate related issues (Barkana 2001, Razoumov et al. 2001, Cen & McDonald 2001, Gnedin 2001, Fan et al. 2001b). Where there is overlap, our results are in broad agreement with these papers. We present a comparision with other authors at the end of §3.2. Our approach here is most similar to that of Cen & McDonald (2001). In addition to obtaining constraints on the ionizing background from the Ly$`\alpha `$ and Ly$`\beta `$ transmission as was considered by Cen & McDonald, we consider the possible constraint from the Gunn-Peterson trough itself, examine the dependence on power spectrum normalization, and develop a method to predict the scatter in mean transmission by relating it to the power spectrum, which might be of wider interest. We also place a slightly stronger emphasis on the neutral fraction, which is more robustly determined compared to the ionizing background or photoionization rate.

## 2 The Gravitational Instability Model for the IGM

The Ly$`\alpha `$ optical depth is related to the IGM density, assuming ionization equilibrium , via

$$\tau _\alpha =A_\alpha (1+\delta )^{20.7(\gamma 1)}$$
(1)
where $`\delta `$ is the gas overdensity ($`\delta =(\rho \overline{\rho })/\overline{\rho }`$, where $`\rho `$ is the gas density and $`\overline{\rho }`$ its mean), $`\gamma `$ is the equation of state index for the IGM222The photoionized IGM at overdensity of a few or less is expected to follow a tight temperature-density relation of the form $`T=T_0(1+\delta )^{\gamma 1}`$, where $`T`$ is the gas temperature and $`T_0`$ is its value at the cosmic mean density (see Hui & Gnedin 1997). We caution that close to reionization, these quantities may not be a function of $`\delta `$ alone. The IGM may be heated inhomogenously, causing spatial fluctuations in $`T_0`$ and $`\gamma `$., and $`A_\alpha `$ is given by (see e.g. Hui et al. 2001 and references therein):

$`A_\alpha =51\left[{\displaystyle \frac{X_{\mathrm{HI}}}{1.6\times 10^4}}\right]\left[{\displaystyle \frac{\mathrm{\Omega }_bh^2}{0.02}}\right]\left[{\displaystyle \frac{0.65}{h}}\right]\left[{\displaystyle \frac{1+z}{7}}\right]^3\left[{\displaystyle \frac{11.7}{H(z)/H_0}}\right]`$ (2)
where $`X_{\mathrm{HI}}n_{\mathrm{HI}}/n_\mathrm{H}`$ ($`n_\mathrm{H}`$ is the total density of neutral and ionized hydrogen, and $`n_{\mathrm{HI}}`$ is the neutral hydrogen density) is the neutral hydrogen fraction at mean density ($`\delta =0`$). 333The neutral fraction at arbitrary $`\delta `$ is given by $`X_{\mathrm{HI}}`$ times $`(1+\delta )^{10.7(\gamma 1)}`$. Throughout this paper, whenever we quote values for $`X_{\mathrm{HI}}`$, we refer to the neutral hydrogen fraction at the cosmic mean density $`\delta =0`$.  Here $`H(z)`$ is the Hubble parameter at redshift $`z`$, $`H_0`$ is the Hubble parameter today, $`H_0=100h`$ km/s/Mpc, $`\mathrm{\Omega }_b`$ is the baryon density in units of the critical density. The value of $`11.7`$ for $`H(z)/H_0`$ above corresponds to that appropriate for a cosmology with $`\mathrm{\Omega }_m=0.4`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$ at $`z=6`$, where $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are the matter and vacuum densities in units of the critical density today.

The neutral fraction $`X_{\mathrm{HI}}`$ is related to the ionizing background by 444This equation assumes that hydrogen is highly ionized and that helium is largely doubly ionized. If helium is only singly ionized, the relation between $`J_{\mathrm{HI}}`$ and $`X_{\mathrm{HI}}`$ will be changed slightly: the right hand side of eq (3) will be multiplied by $`0.93`$.

$$X_{\mathrm{HI}}=1.6\times 10^4\left[\frac{\mathrm{\Omega }_bh^2}{0.02}\right]\left[\frac{2.55\times 10^2}{J_{\mathrm{HI}}}\right]\left[\frac{T_0}{2\times 10^4\mathrm{K}}\right]^{0.7}\left[\frac{1+z}{7}\right]^3$$
(3)
where the dimensionless quantity $`J_{\mathrm{HI}}`$ is related to the photoionization rate $`\mathrm{\Gamma }_{\mathrm{HI}}`$ by

$$\mathrm{\Gamma }_{\mathrm{HI}}=4.3\times 10^{12}J_{\mathrm{HI}}\mathrm{s}^1$$
(4)
The quantity $`J_{\mathrm{HI}}`$ provides a convenient way of describing the normalization of the ionizing background, without specifying the exact spectrum, in a way that is directly related to the physically relevant $`\mathrm{\Gamma }_{\mathrm{HI}}`$ (e.g. Miralda-Escudé et al. 1996). It is related to the specific intensity at $`912\AA `$ $`j_{\nu _{\mathrm{HI}}}`$ by $`J_{\mathrm{HI}}=j_{\nu _{\mathrm{HI}}}\times [3/(\beta +3)]`$, where $`\beta `$ is the slope of the specific intensity just blueward of $`912\AA `$ ($`j_\nu \nu ^\beta `$ where $`\nu `$ is frequency), and $`j_{\nu _{\mathrm{HI}}}`$ is measured in the customary units of $`10^{21}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{Hz}^1\mathrm{ster}^1`$ (for non-power law $`j_\nu `$, eq.  provides the exact definition for $`J_{\mathrm{HI}}`$; see e.g. Hui et al. 2001).

Two more ingredients should be mentioned to complete the specification of our model for the Ly$`\alpha `$ absorption (see e.g. Hui et al. 1997 for details). First, the optical depth as a function of velocity is computed by taking the right hand side of eq. (1) in velocity space (i.e. taking into account peculiar velocities) and smoothing it with a thermal broadening window. Second, the gas density and velocity fields are predicted by some Cold Dark Matter (CDM) cosmological model using numerical simulations.

There are obviously a number of free parameters in our model. Let us discuss each of them in turn.

Throughout this paper, we assume $`\mathrm{\Omega }_bh^2=0.02`$, as supported by recent cosmic microwave background measurements (Netterfield et al. 2001, Pryke et al. 2001) and the nucleosynthesis constraint from primordial deuterium abundance (Burles et al. 2001). We also assume throughout $`h=0.65`$, $`\mathrm{\Omega }_m=0.4`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$. Variations of these parameters within the current bounds do not contribute significantly to the uncertainties of the constraints obtained in this paper (see Hui et al. 2001).

The temperature $`T_0`$ and equation of state index $`\gamma `$ at the redshifts of interest in this paper are somewhat uncertain. There are no direct measurements of the thermal state of the IGM at our redshifts of interest, $`z>\mathrm{\hspace{0.33em}4}`$. Measurements at $`z<\mathrm{\hspace{0.33em}4}`$ yield values consistent with $`T_0=2\times 10^4`$ K and $`\gamma =1`$ (McDonald et al. 2000b, Ricotti et al. 2000, Zaldarriaga et al. 2001. Schaye et al. 2000, however, measure a slightly lower temperature). Given that the temperature right after reionization is expected to be about $`25000`$ K with $`\gamma =1`$ (with some dependence on the hardness of the ionizing spectrum; see e.g. Hui & Gnedin 1997), which is not too different from the measurements at $`z<\mathrm{\hspace{0.33em}4}`$, we will assume throughout this paper, when making use of eq. (3) to infer $`J_{\mathrm{HI}}`$, that $`T_0=2\times 10^4`$ K and $`\gamma =1`$. Note that while the theoretically allowed range for $`\gamma `$ is from $`1`$ to $`1.6`$ (Hui & Gnedin 1997), what matters for our purpose is $`20.7(\gamma 1)`$ (eq. ), which only ranges from $`1.58`$ to $`2`$, and does not significantly affect our results. It is also important to emphasize that the inference of $`X_{\mathrm{HI}}`$ from observations, unlike the case for $`J_{\mathrm{HI}}`$, is not subject directly to uncertainties in the temperature $`T_0`$. This is because observations constrain $`A_\alpha `$ from which we can obtain $`X_{\mathrm{HI}}`$ without knowing $`T_0`$ (see eq. ). 555 The above statement is subject to two small caveats. First, the optical depth given in eq. (1) has to be smoothed with a thermal broadening window whose width depends on $`T_0`$. We find that in practice, the exact width of the thermal broadening kernel does not affect very much quantities such as the mean transmission, which is what we will be interested in. Second, $`T_0`$ also affects the gas dynamics via the pressure term in the equation of motion. As we will discuss below, the effect of varying $`T_0`$ also appears to be small in this regard.

To generate realizations of the density and velocity fields for a given cosmology, we run Hydro-Particle-Mesh (HPM) simulations (Gnedin & Hui 1998). The HPM algorithm is essentially a Particle-Mesh code, modified to incorporate a force term due to gas pressure in the equation of motion. 666The temperature-density relation has to be specified as a function of redshift in the HPM code to compute the pressure term. We follow McDonald $`\&`$ Miralda-Escudé (2001) and linearly interpolate between $`T_0=19000`$ K and $`\gamma =1.2`$ at $`z=3.9`$ and $`T_0=25000`$ K and $`\gamma =1`$ at the redshift of reionization $`z_{\mathrm{reion}}`$. We found that assuming $`z_{\mathrm{reion}}=7`$ versus $`z_{\mathrm{reion}}=10`$ results in negligible difference in our results, in particular concerning the mean decrement and the probability distribution of transmission. All results in this paper are quoted from the $`z_{\mathrm{reion}}=7`$ HPM simulations. Note that in inferring $`X_{\mathrm{HI}}`$ and $`J_{\mathrm{HI}}`$ from eq. (2) and (3), we always use $`T_0=20000`$ K and $`\gamma =1`$ for simplicity, as mentioned before. For the initial power spectrum, we use a Cold Dark Matter (CDM) type transfer function, as parameterized by Ma (1996), which is very similar to the commonly used Bardeen et al. (1986) transfer function. For the primordial spectral slope, we adopt $`n=0.93`$ (Croft et al. 2000, McDonald et al. 2000a). For the linear power spectrum normalization, we employ the range suggested by measurements from the Ly$`\alpha `$ forest of Croft et al. (2000): $`\mathrm{\Delta }^2(k)4\pi k^3P(k)/(2\pi )^3=0.74_{0.16}^{+0.20}`$ at $`z=2.72`$ at a velocity scale of $`k=0.03(\mathrm{km}/\mathrm{s})^1`$.777This normalization corrects an error in an earlier draft of Croft et al. (2000). (R. Croft, private communication.) We, however, caution that the error-bar given is somewhat dependent on the assumed error of the mean transmission measurements, which is sensitive to the accuracy of the continuum-fitting procedure (see e.g. Zaldarriaga et al. 2001 for a slightly different assessment of the error-bar). The power spectrum in this model has a similar shape to that of favored cosmological models, but slightly lower amplitude (Croft et al. 2000). In §3.2 we demonstrate that our main conclusion, that the neutral fraction increases dramatically near $`z6`$, is insensitive to our assumptions about the amplitude of the power spectrum. In practice, we examine models with different normalizations by running a simulation with outputs at several different redshifts: each redshift then corresponds to a different power spectrum normalization, and linear interpolation is performed to reach any desired normalization. 888We do not vary the primordial spectral index $`n`$ here. Quantities such as the mean transmission which we are interested in here are generally sensitive to power on only a small range of scales. Varying $`n`$ is therefore largely degenerate with varying $`\mathrm{\Delta }^2`$.

Our simulations have a box size of $`8.9`$ Mpc/h, with a $`256^3`$ grid. McDonald & Miralda-Escude (2001) found this box size and resolution to be adequate for IGM studies up to $`z5`$. We have verified that the same is true up to $`z=6`$, in the sense that the transmission probability distribution has converged for our choice of simulation size and resolution.

Finally, we should say a few words about simulations of the Ly$`\beta `$ region. In regions of the quasar spectrum that are between $`973\AA (1+z_{\mathrm{em}.})`$ and $`1026\AA (1+z_{\mathrm{em}.})`$, where $`z_{\mathrm{em}.}`$ is the redshift of the quasar, two kinds of absorption can exist: one is Ly$`\beta `$ due to material at redshift $`0.948(1+z_{\mathrm{em}.})<1+z<1+z_{\mathrm{em}.}`$, the other is Ly$`\alpha `$ due to material at redshift $`0.800(1+z_{\mathrm{em}.})<1+z<0.844(1+z_{\mathrm{em}.})`$. In other words, in such a region, the observed optical depth would be given by $`\tau =\tau _\beta +\tau _\alpha `$ where $`\tau _\beta `$ and $`\tau _\alpha `$ arises at different redshifts. The Ly$`\alpha `$ optical depth can be computed as before. The Ly$`\beta `$ optical depth $`\tau _\beta `$ can be computed using eq. (1), except that $`A_\alpha `$ is replaced by $`A_\beta `$:

$$A_\beta =\frac{1}{5.27}A_\alpha $$
(5)
The factor of $`5.27`$ reflects the fact that the Ly$`\beta `$ transition has a cross-section that is $`5.27`$ times smaller than Ly$`\alpha `$.

## 3 Constraints on the Neutral Hydrogen Fraction and the Ionizing Background

### 3.1 Constraints from the Ly$`\alpha `$ Mean Transmission

Using eq. (1) and (2), we compute the $`X_{\mathrm{HI}}`$, which also fixes $`J_{\mathrm{HI}}`$ (eq. ), necessary to match the observed Ly$`\alpha `$ mean transmission $`e^{\tau _\alpha }`$ at $`z=4.56`$ (see Table 1 for a summary of the measurements). The results of our calculation are presented in Fig. 1. This plot also contains a point at $`z=6.05`$ which is the result of matching the mean transmission in the Ly$`\beta `$ forest, as we describe in §3.2.

Also shown in the figure is a dotted line which shows $`X_{\mathrm{HI}}(1+z)^3`$, which appears to be a good fit to the data from $`z=4.5`$ to $`z=5.7`$. From eq. (3), one can see that such a trend for the neutral fraction is equivalent to assuming constant $`J_{\mathrm{HI}}`$ (or more accurately, constant $`J_{\mathrm{HI}}T_0^{0.7}`$; see eq. ).

As one can see, ignoring for now the Ly$`\beta `$ point, the neutral fraction does appear to have a modest jump around $`z6`$: it increases by a factor of $`4.0`$ from $`z=5.7`$ to $`z=6.05`$, while it changes at most by $`1.9`$ from $`z=4.5`$ to $`z=5.7`$. A similar trend (but opposite in sign) can be seen in the ionizing flux $`J_{\mathrm{HI}}`$. The 1 $`\sigma `$ error-bar here takes into account the measurement error in mean transmission, and the range of power spectrum normalization stated in §2. As we have explained in §2, while $`X_{\mathrm{HI}}`$ is not sensitive to the assumed temperature of the IGM ($`T_0`$), our constraints on $`J_{\mathrm{HI}}`$ are directly influenced by it. As emphasized before, we assume $`T_0=2\times 10^4`$ K. In other words, our constraints on $`J_{\mathrm{HI}}`$ are really constraints on the quantity $`J_{\mathrm{HI}}(T_0/2\times 10^4\mathrm{K})^{0.7}`$ (see eq. ). It is therefore straightforward to rescale our constraints on $`J_{\mathrm{HI}}`$ if the temperature were a little bit different. 999The temperature also affects the thermal broadening window, but we find that in practice its effect on our constraints on $`A_\alpha `$ (eq. ) is small; see discussion in §2. It is an interesting question to ask whether the apparent jump in the ionizing flux can instead be attributed to a jump in the temperature. In general, the temperature is expected to evolve slowly with redshift after reionization (Hui & Gnedin 1997).

Regarding the measurement error, we should also emphasize that the Becker et al.’s 2 $`\sigma `$ error actually includes the possibility of having zero transmission at $`z=6.05`$. This means that at 2 $`\sigma `$, we would only have a lower limit on $`X_{\mathrm{HI}}`$, or an upper limit on $`J_{\mathrm{HI}}`$, for the highest redshift point in Fig. 1, allowing the possibility that the IGM is neutral at $`z6`$, $`X_{\mathrm{HI}}=1`$.

### 3.2 Constraints from the Ly$`\beta `$ Mean Transmission

In this section we consider the constraints placed by Becker et al.’s measurement of the mean transmission in the Ly$`\beta `$ region. Absorption in the Ly$`\beta `$ region has two components: $`\tau =\tau _\alpha +\tau _\beta `$, where the Ly$`\alpha `$ optical $`\tau _\alpha `$ and the Ly$`\beta `$ optical depth $`\tau _\beta `$ originate at different redshifts. Ly$`\beta `$ absorption due to material at $`z=6`$ coincides in wavelength with Ly$`\alpha `$ absorption due to material at $`z=(1+6)\times 1026/12161=4.9`$. Because the points of origin are so widely separated, they can be effectively treated as statistically independent i.e. $`e^\tau =e^{\tau _\alpha }e^{\tau _\beta }`$. Becker et al. measured $`e^{\tau _\beta }`$ at $`z6`$ by dividing the net mean transmission $`e^\tau `$ in the Ly$`\beta `$ region by the mean transmission in Ly$`\alpha `$ at $`z5`$. They obtained $`e^{\tau _\beta }=0.002\pm 0.020`$. Clearly, this measurement is consistent with a completely neutral IGM. However, the interesting question is: what kind of lower limit does it set on the neutral fraction, and does it improve upon the lower limit from the mean Ly$`\alpha `$ absorption ?

We carry out a calculation that is analogous to what is described in §3.1, except for the key difference that in computing $`\tau _\beta `$ from the simulated density and velocity fields, we employ $`A_\beta `$ which is a factor of $`5.27`$ smaller than $`A_\alpha `$ (see eq.  & ). The results of our calculation are shown in Fig. 1, as the highest redshift points in the plot, which have error bar arrows pointing towards a completely neutral IGM and a vanishing ionizing background. It can be seen that the ($`1\sigma `$) lower limit on $`X_{\mathrm{HI}}`$ is $`X_{\mathrm{HI}}>4.7\times 10^4`$. This is a factor of $`3`$ larger than the neutral fraction required to match the upper limit on the mean transmission in Ly$`\alpha `$ for our fiducial normalization, and a slightly stronger constraint than that obtained in section 3.1, including the uncertainty in power spectrum normalization. Similarly, the upper limit on $`J_{\mathrm{HI}}`$ is $`J_{\mathrm{HI}}<9.0\times 10^3`$. The moral here is that because the Ly$`\beta `$ absorption cross-section is a factor of $`5.27`$ smaller than the Ly$`\alpha `$ cross-section, Ly$`\beta `$ offers a more sensitive probe of the neutral fraction, especially when the Ly$`\alpha `$ optical depth is high.

The neutral fraction at $`z6`$ is thus a factor of $`10`$ higher than that at redshift $`z5.7`$, where it is $`X_{\mathrm{HI}}=4.9\times 10^5`$. This dramatic change in the neutral fraction is suggestive, probably indicating that the reionization epoch is nearby.

Furthermore, this conclusion is not sensitive to our assumptions about the amplitude of the power spectrum. Although the neutral fraction at redshift $`z=6.05`$ is itself sensitive to the amplitude of the power spectrum, we find that the factor by which the neutral fraction increases from $`z=5.7`$ to $`z=6.05`$ depends only weakly on the amplitude. In Fig. 2 we plot both the neutral fraction at $`z=6.05`$ and the jump in the neutral fraction for a range of different power spectrum normalizations. The jump is defined as the ratio $`X_{\mathrm{HI}}(z=6.05)/X_{\mathrm{HI}}(z=5.7)`$. Here $`X_{\mathrm{HI}}(z=6.05)`$ is the lower limit resulting from the 1 $`\sigma `$ error in the mean transmission in Ly$`\beta `$ at $`z=6.05`$ and the error bars in the jump arise from the 1 $`\sigma `$ error in the mean transmission in Ly$`\alpha `$ at $`z=5.7`$. As one can see in the plot, the lower limit on the neutral fraction at $`z=6.05`$ varies from $`X_{\mathrm{HI}}>3\times 10^4`$ to $`X_{\mathrm{HI}}>9\times 10^4`$ as $`\mathrm{\Delta }^2(k=0.03(\mathrm{km}/\mathrm{s})^1,z=2.72)`$ varies from $`0.5`$ to $`1.3`$. The neutral fraction itself varies significantly with power spectrum normalization, scaling approximately as $`X_{\mathrm{HI}}[\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)]^{1.1}`$, for this range of normalizations. The jump, however, changes only slightly over a large range of normalizations. As $`\mathrm{\Delta }^2(k=0.03(\mathrm{km}/\mathrm{s})^1,z=2.72)`$ varies from $`0.5`$ to $`1.3`$, the jump changes only from $`9.7`$ to $`11.1`$. Our conclusion that the neutral fraction of the IGM increases dramatically near $`z6`$ seems robust.

One can also consider the absorption in the Ly$`\gamma `$ region, or even the higher Lyman series. In practice, the accumulated amount of absorption from Ly$`\alpha `$ as well as Ly$`\beta `$ at different redshifts makes it harder to measure the Ly$`\gamma `$ transmission itself with good accuracy.

Our constraints on the neutral fraction and the intensity of the ionizing background are consistent with those found by other authors, given our different choices of power spectrum normalization. Fan et al. (2001b) found, from the mean Ly$`\beta `$ transmission, that $`\mathrm{\Gamma }_{12}<0.025`$, where $`\mathrm{\Gamma }_{12}`$ is the photoionization rate of eq (4) in units of $`10^{12}\mathrm{s}^1`$. Although this constraint is somewhat stronger than the constraint implied by our fiducial model, $`\mathrm{\Gamma }_{12}<0.039`$, we expect the difference due to our different power spectrum normalizations. Fan et al.’s (2001b) constraint comes from semi-analytic arguments, shown consistent with an LCDM simulation with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$,$`h=0.65`$, $`\mathrm{\Omega }_bh^2=0.02`$, and $`\sigma _8=0.9`$. This model has a substantially larger normalization, $`\mathrm{\Delta }^2(k=0.03(\mathrm{km}/\mathrm{s})^1,z=2.72)=1.25`$, than our fiducial model of $`\mathrm{\Delta }^2(k=0.03(\mathrm{km}/\mathrm{s})^1,z=2.72)=0.74`$. The difference in normalization reflects some tension between the normalization derived from the observed cluster abundance, which Fan uses, and that from the Lyman-$`\alpha `$ forest which our model is based on (Croft et al. 2000). Fan et al. (2001b) assume $`T_0=2.0\times 10^4`$ K in placing their constraint. Their limit, $`\mathrm{\Gamma }_{12}<0.025`$, includes only uncertainties in the mean transmission and not additional uncertainties from the power spectrum normalization. From figure (2), we infer that Fan et al.’s (2001b) normalization implies $`X_{\mathrm{HI}}>8.8\times 10^4`$ in our cosmology. Rescaling this result from our assumed $`\mathrm{\Omega }_m=0.4`$ to an $`\mathrm{\Omega }_m=0.3`$ cosmology, and using eqs. (3) and (4) we predict $`\mathrm{\Gamma }_{12}<0.024`$, or $`X_{\mathrm{HI}}>7.6\times 10^4`$, for Fan et al.’s (2001b) model. The constraint of Fan et al. (2001b) is thus consistent with our constraint given our different choice of normalization. Cen & McDonald (2001), using a model similar to that of Fan et al. (2001b), obtained the constraint $`\mathrm{\Gamma }_{12}<0.032`$, using the Ly$`\beta `$ mean transmission. The constraint is slightly weaker than that of Fan et al. (2001b) and our extrapolation to their normalization, because Cen & McDonald (2001) consider a larger upper limit to the observed mean transmission, including an estimate of the uncertainty due to sky subtraction. At slightly lower redshifts, we can also compare with the results of McDonald & Miralda-Escudé (2001) derived from matching the mean Ly$`\alpha `$ transmission. For example, at $`z=5.2`$, these authors found $`\mathrm{\Gamma }_{12}=0.16`$ to match the observed mean transmission of $`e^{\tau _\beta }=0.09`$. McDonald & Miralda-Escudé (2001) consider a model whose normalization we infer to be $`\mathrm{\Delta }^2(k=0.03(\mathrm{km}/\mathrm{s})^1,z=2.72)=0.98`$. To match the same mean transmission with this normalization we infer a somewhat higher photoionization rate, $`\mathrm{\Gamma }_{12}=0.19`$. Part of the difference may be that the $`\mathrm{\Gamma }_{12}`$ necessary to match a given mean transmission varies by $`5\%`$ between two different realizations of the density field. The remaining difference may come from the procedure of linearly interpolating between outputs or from some modeling difference. At any rate, our results are roughly consistent with those of other authors given our different power spectrum normalizations.

### 3.3 Constraints from the Gunn-Peterson Trough Itself – the Fluctuation Method

The fact that Becker et al. (2001) observed a Gunn-Peterson trough, where a long stretch of the spectrum contains little or no flux, can conceivably be used to further tighten the constraints obtained from the previous sections. Since the IGM is expected to have spatial fluctuations, the probability of having many pixels in a row turning up a very small transmission must be low, unless the neutral fraction is intrinsically quite high. The same reasoning can be applied to either the Ly$`\alpha `$ or Ly$`\beta `$ absorption. We will discuss our method for Ly$`\alpha `$ in detail. The method for Ly$`\beta `$ is a straightforward extension. For simplicity, we will call this method, the fluctuation method.

Becker et al. (2001) finds from the spectrum of SDSS 1030+0524, the $`z=6.28`$ quasar, the Ly$`\alpha `$ transmission is consistently below about 0.06 for a region that spans $`260\AA `$, between $`8450\AA `$ to $`8710\AA `$. The noise level per $`4\AA `$ pixel is $`\sqrt{n^2}0.02`$, where $`n`$ represents the photon noise fluctuation. 101010We estimate the noise per pixel from Becker et al.’s error-bar in the mean transmission, which is $`0.003`$. This is estimated from a chunk of the spectrum which is $`260\AA `$ long, and so the dispersion per $`4\AA `$ pixel should be approximately $`\sqrt{n^2}0.003\times \sqrt{65}0.02`$. Note that the actual dispersion varies across the spectrum, but this should suffice as a rough estimate. This estimate also agrees with an estimate of the error by comparing Fig. 1 and Fig. 3 of Becker et al. The observed transmission $`F`$ at a given pixel is $`F=e^\tau +n`$, where $`e^\tau `$ is the true transmission. The noise here should be dominated by Poisson fluctuations of the subtracted sky background (as well as perhaps read-out error). Let $`P(F_1,F_2,\mathrm{}F_N)dF_1\mathrm{}dF_N`$ be the probability that $`N`$ consecutive pixels have observed transmission fall into the range $`F_1\pm dF_1/2`$$`F_N\pm dF_N/2`$. In our case, $`N=65`$ for the pixel size of $`4\AA `$. The problem is then to find the probability $`_{<0.06}\mathrm{}_{<0.06}P(F_1\mathrm{}F_N)𝑑F_1\mathrm{}𝑑F_N`$ as a function of $`J_{\mathrm{HI}}`$, and ask what maximal $`J_{\mathrm{HI}}`$ (or equivalently, minimal $`X_{\mathrm{HI}}`$) would give an acceptable probability. By choosing the “acceptable probability” to be within $`68\%`$ of the maximum likelihood (maximum likelihood is achieved when the neutral fraction is unity), we obtain the 1 $`\sigma `$ upper limit on $`J_{\mathrm{HI}}`$ or 1 $`\sigma `$ lower limit on $`X_{\mathrm{HI}}`$.

Our simulation has a comoving box size of $`8.9`$ Mpc/h, corresponding to $`42\AA `$ for Ly$`\alpha `$ at $`z6`$, which falls short of the wavelength range we need for this problem, which is $`260\AA `$. In other words, the probability $`P(F_1,F_2,\mathrm{}F_N)dF_1\mathrm{}dF_N`$ can be estimated directly from our simulation only for $`N10`$. However, the mass correlation length scale at this redshift ($`<\mathrm{\hspace{0.33em}1}`$ Mpc/h) is actually a fraction of the box size, which means one can treat fluctuations on scales beyond the box size as roughly uncorrelated. Assuming so, we estimate $`_{<0.06}\mathrm{}_{<0.06}P(F_1,F_2,\mathrm{}F_{10})𝑑F_1\mathrm{}𝑑F_{10}`$ from the simulation, and then take its $`6`$-th power, which would give us the probability that $`60`$ consecutive pixels have transmission below $`0.06`$. This is slightly smaller than the number $`65`$ that we need, but at least will provide us conservative constraints on $`X_{\mathrm{HI}}`$ and $`J_{\mathrm{HI}}`$. We have also tested our approach by using fractions of the box-size as a unit, and find that our results do not change significantly (less than $`10\%`$).

Fig. 3a (dotted curve) shows our estimate of the probability $`_{<0.06}\mathrm{}_{<0.06}P(F_1,F_2,\mathrm{}F_N)𝑑F_1\mathrm{}𝑑F_N`$ for $`N=60`$ and pixel size $`4\AA `$, as a function of $`J_{\mathrm{HI}}`$. Our simulated spectra have been convolved with the observation resolution (full-width-at-half-maximum of $`1.8\AA `$), rebinned into pixels of $`4\AA `$ each and added Gaussian noise with a dispersion of $`0.02`$. From the dotted curve in Fig. 3a, applying a likelihood analysis, we obtain a 1 $`\sigma `$ upper limit on $`J_{\mathrm{HI}}`$ of $`J_{\mathrm{HI}}<0.014`$, and a corresponding lower limit on $`X_{\mathrm{HI}}`$ of $`X_{\mathrm{HI}}>2.95\times 10^4`$. This is for a model with a power spectrum normalization of $`\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)=0.74`$ (see §2). The mean Ly$`\alpha `$ transmission constraints for the same model are $`J_{\mathrm{HI}}<0.028`$ and $`X_{\mathrm{HI}}>1.5\times 10^4`$. 111111Do not confuse these constraints, which are for the particular power spectrum normalization mentioned above, to the constraints discussed in earlier sections, which include the uncertainty in the power spectrum normalization. We focus on a single model in this section for simplicity. This means that considering Ly$`\alpha `$ alone, the fluctuation method yields somewhat stronger constraints compared to using simply the mean transmission.

Fig. 3b (dotted curve) shows the same methodology applied to the Ly$`\beta `$ Gunn-Peterson trough. A new ingredient here is that one needs an additional simulation of the same model at redshift $`z=4.9`$ to produce the Ly$`\alpha `$ absorption that can be overlaid on top of the Ly$`\beta `$ absorption from $`z=6.05`$. This additional simulation should have different initial phases to mimic the fact that fluctuations at $`z=4.9`$ and those at $`z=6.05`$ should be uncorrelated. We obtain 1 $`\sigma `$ limits of $`J_{\mathrm{HI}}<0.012`$ and $`X_{\mathrm{HI}}>3.4\times 10^4`$. This is about $`40\%`$ weaker than the constraints we obtain from the Ly$`\beta `$ mean transmission. In other words, from Ly$`\beta `$ absorption, the fluctuation method yields slightly weaker constraints compared to using the mean transmission. It is also only slightly stronger than the constraint obtained from the fluctuation method applied to Ly$`\alpha `$.

It is an interesting question to ask how many sightlines one would need to improve the constraints by, say a factor of $`2`$. Our approach can be easily extended to multiple (uncorrelated) sightlines, and we find that about $`5`$ sightlines (each containing a Gunn-Peterson trough of the same length and same signal to noise) are necessary for such an improvement.

Part of the difficulty with obtaining stronger constraints, in addition to the small number of sightlines, is the dominance of noise. The lower panel of Fig. 4 shows the one-pixel ($`4\AA `$) probability distribution function (PDF) of the true transmission $`e^\tau `$ (i.e. no noise added), for three different values of $`J_{\mathrm{HI}}`$ (the power spectrum normalization is the same as that in Fig. 3). The upper panel shows the corresponding PDF’s of the observed transmission $`F`$ (i.e. after convolving $`P(e^\tau )`$ with a Gaussian of dispersion $`0.02`$). As expected, noisy data make the PDF’s more similar. Nonetheless, as we pointed out above, with sufficient number of sightlines, there might be a non-negligible chance of seeing pixels with high transmission that take place at the tail of the PDF’s, hence allowing us to distinguish between the different levels of the ionizing background. Alternatively, one can try improving the signal-to-noise per pixel. In Fig. 3b, we show with a dashed curve the corresponding probability if the noise per pixel is lowered by a factor of $`4`$. The constraints improve by a little more than a factor of $`2`$. We should emphasize, however, systematic errors are likely important here – we will discuss them in the next two sections.

## 4 The Variance of the Mean Transmission

If, as is suggested by our discussion in §3.2 (see Fig. 1), the IGM is close to the epoch of reionization at $`z6`$, one might expect large fluctuations in the ionizing background near that time. For instance, one line of sight might probe a region of the IGM where the ionized bubbles around galaxies or quasars have percolated, while another might probe the pre-percolation IGM. As mentioned before, the simulations we employ do not take into account fluctuations in the ionizing background. (For simulations incorporating radiative transfer see e.g., Gnedin & Abel 2001, Razoumov 2001). One useful check would then be to predict the sightline to sightline scatter in mean transmission from our simulations, and compare that against the observed scatter. At $`z5.5`$, 4 lines of sight are available for a measurement of the scatter. We will examine this, as well as make predictions for the scatter at $`z6`$, which more high redshift quasars in the future will allow us to measure.

Our estimate relies on simulation measurements of the transmission power spectrum. This is in contrast to an estimate of the same quantity made by Zuo (Zuo 1993) who makes a prediction based on extrapolations of the column density distribution and of the number of absorbing clouds per. unit redshift (Zuo & Phinney 1993). Zuo also assumes that the clouds are Poisson distributed, while our measurement incorporates the clustering in the IGM via our numerical simulation.

An immediate problem presents itself: sightlines from which the mean transmission is measured are typically longer than the usual simulation box. We tackle this problem by expressing the variance of mean transmission in terms of the transmission power spectrum, and making use of a reasonable assumption about the behavior of the power spectrum on large scales.

The mean transmission from one sightline is estimated using

$$\overline{F}=\frac{1}{N}\underset{i=1}{\overset{N}{}}F_i$$
(6)
where $`N`$ is the number of pixels, $`F_i`$ is the observed transmission at pixel $`i`$, $`F_i=f_i+n_i`$ where $`f=e^\tau `$ is the true transmission and $`n`$ is the noise fluctuation. We use the symbol $`\overline{F}`$ to represent the estimator, and $`\overline{f}`$ to denote the true mean transmission. The variance of the estimated mean transmission is then

$`\sigma _T^2\overline{F}^2\overline{F}^2={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j}{}}[F_iF_jF_iF_j]`$ (7)
$`={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j}{}}\xi _{ij}+{\displaystyle \frac{1}{N}}\sigma _n^2=2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{2\pi }}\left[{\displaystyle \frac{\mathrm{sin}(kL/2)}{kL/2}}\right]^2P_f(k)+{\displaystyle \frac{\sigma _n^2}{N}}`$
where $`\sigma _n^2n^2`$ is assumed roughly independent of position, and $`\xi _{ij}`$ is the un-normalized two-point correlation of the transmission i.e. $`\xi _{ij}f_if_j\overline{f}^2`$, and $`P_f(k)`$ is its one-dimensional Fourier transform. The symbol $`L`$ denotes the comoving length of the spectrum from which the mean transmission is measured, and $`k`$ is the comoving wavenumber.

To evaluate $`\sigma _T`$, we need to know the transmission power spectrum on scales generally larger than the size of the typical simulation box. It is expected that the transmission power spectrum takes the shape (not the normalization) of the linear mass power spectrum on large scales (i.e. essentially linear biasing; see Scherrer & Weinberg 1998, Croft et al. 1997, Hui 1999). We therefore use this to extrapolate the simulation $`P_f(k)`$ to large scales (small $`k`$’s). We find that $`P_f(k)`$ is well approximated by $`P_f(k)=B\mathrm{exp}(ak^2)_k^{\mathrm{}}(dk/2\pi )kP_{\mathrm{mass}}(k)`$ where $`P_{\mathrm{mass}}`$ is the three-dimensional linear mass power spectrum.

Becker et al. gave an estimate of $`\sigma _T0.03\pm 0.01`$, $`\overline{f}=0.1`$, at $`z=5.5`$ using 4 different sightlines, each spanning $`\mathrm{\Delta }z=0.2`$, which corresponds to $`L57`$ Mpc/h. 121212The error on $`\sigma _T`$ is estimated assuming Gaussian statistics and that the four lines of sight are independent. Then $`var(\sigma _T)=\sigma _T^2/2n`$ (see e.g. Kendall & Stuart 1958). An estimate of the noise term is provided by the error in the mean transmission, $`(\sigma _n^2/N)^{0.5}0.003`$. For the $`\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)=0.74`$ case, the fitting parameters are $`B=0.033`$ and $`a=0.013\mathrm{Mpc}^2/\mathrm{h}^2`$. Using eq. (5), we then find $`\sigma _T=0.030`$ for $`\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)=0.74`$, $`\sigma _T=0.031`$ for $`\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)=0.94`$, and $`\sigma _T=0.028`$ for $`\mathrm{\Delta }^2(k=0.03\mathrm{s}/\mathrm{km},z=2.72)=0.58`$. The variance is similar between the different normalizations because each normalization requires a different $`A_\alpha `$ in eq. (1) to match the mean transmission. This difference in $`A_\alpha `$ probably compensates for the effect of the different normalizations on $`\sigma _T`$. The predicted scatter of $`\sigma _T0.030`$ is consistent with the measured $`\sigma _T`$ of $`0.03\pm 0.01`$.

We apply the same methodology as the above to estimate $`\sigma _T`$ at $`z6`$. In Fig. 5, we show the results for a range of different $`J_{\mathrm{HI}}`$’s for each of our canonical power spectrum normalizations, ($`\mathrm{\Delta }^2(k=\mathrm{\hspace{0.17em}0.03}\mathrm{s}/\mathrm{km},z=2.72)=0.58,0.74`$ and $`0.94`$). 131313The comparison across normalizations is done here at fixed $`J_{\mathrm{HI}}`$ while at $`z=5.5`$ we compared the results of different normalizations at fixed mean transmission. We find that the dependence on normalization is larger at fixed $`J_{\mathrm{HI}}`$ than at fixed mean transmission. Photon noise is not included in the estimates of this figure. Even for relatively large $`J_{\mathrm{HI}}`$’s the scatter is small. For instance, for $`J_{\mathrm{HI}}=4.5\times 10^2`$, $`\sigma _T=1.1\times 10^2`$, assuming our fiducial normalization. This $`J_{\mathrm{HI}}`$ is large in that it already gives a mean transmission, $`\overline{f}=1.75\times 10^2`$, in excess of the observations. By $`J_{\mathrm{HI}}=1.4\times 10^2`$, the scatter is only $`\sigma _T=2.8\times 10^3`$ for our fiducial normalization. The scatter depends somewhat on normalization, as one can see in the figure. To measure the scatter well would require data that are less noisy than the one discussed here, which has photon noise of $`(\sigma _n^2/N)^{0.5}0.003`$, comparable to the predicted scatter.

On the other hand, the smallness of this scatter makes it a possibly interesting diagnostic. As we have emphasized before, this predicted scatter ignores fluctuations in the ionizing background. For sufficiently small $`J_{\mathrm{HI}}`$’s, the IGM should be close to the epoch of reionization, and one would expect large sightline by sightline variations. An observed scatter well in excess of what is predicted would be an interesting signature.

## 5 Discussion

Our findings are summarized as follows.

* The most stringent (1 $`\sigma `$) lower limit on the neutral hydrogen fraction $`X_{\mathrm{HI}}`$ (eq. ) or upper limit on the ionizing background $`J_{\mathrm{HI}}`$ (eq. ) at $`z6`$ is obtained from the observed mean Ly$`\beta `$ transmission: $`X_{\mathrm{HI}}>4.7\times 10^4`$. A comparison of this limit versus constraints at lower redshifts is presented in Fig. 1. The fact that the neutral fraction increases by a factor of $`10`$ from redshift of $`5.7`$ to $`6`$ even though it changes by no more than a factor of about $`2`$ from $`z=4.5`$ to $`z=5.7`$ suggests that $`z6`$ might be very close to the epoch of reionization. We emphasize that current constraints are still consistent with a highly ionized IGM at $`z6`$ – it is the steep rise in $`X_{\mathrm{HI}}`$ that is suggestive of dramatic changes around or just before that redshift. We should also mention that the constraints on $`X_{\mathrm{HI}}`$ are less subject to uncertainties in the IGM temperature compared to those on $`J_{\mathrm{HI}}`$ (see §2).
* The existence of a long Gunn-Peterson (Ly$`\alpha `$ or Ly$`\beta `$) trough at $`z6`$, where little or no flux is detected, can also be used to obtain constraints on $`X_{\mathrm{HI}}`$ or $`J_{\mathrm{HI}}`$. This we call the fluctuation method: the fact that a long stretch of the spectrum exhibits no large upward fluctuations in transmission provides interesting information on the neutral fraction or ionizing background. The constraints obtained this way turn out to be fairly similar to those obtained using the mean transmission. We estimate that a reduction in noise by a factor of 4, or an increase in number of sightlines to 5, would result in constraints that are 2 times stronger (§3.3).
* We develop a method to predict the dispersion in mean transmission measured from sightlines that are longer than the typical simulation box (eq.  and Fig. 5). Our predicted dispersion is consistent with that observed at $`z=5.5`$ (Becker et al. 2001). We also predict the scatter at redshift $`z=6`$, which can be measured when more sightlines become available. Assuming a spatially homogeneous ionizing background, we predict a small scatter at $`z=6`$, $`\sigma _T`$ a few $`\times 10^3`$, neglecting photon noise. The dispersion provides a useful diagnostic of fluctuations in the ionizing background – close to the epoch of reionization, one expects large fluctuations from one line of sight to another depending on whether it goes through regions of the IGM where percolation of HII regions has occurred.

There are at least three issues that will be worth exploring. First, with more quasars at $`z6`$ or higher discovered in the future, applying some of the ideas mentioned above would be extremely interesting, such as the measurement of the line of sight scatter in mean transmission, or the use of the Gunn-Peterson trough to obtain stronger constraints on the neutral fraction. Second, as we have commented on before, fluctuations in the ionizing background are expected to be important as we near the epoch of reionization. We have not discussed it here, but a calculation of the size of these fluctuations would be very interesting. Such a calculation will depend both on the mean free path of the ionizing photons as well as the spatial distribution of ionizing sources. The latter is probably quite uncertain, but useful estimates might be made (e.g. Razoumov et al. 2001). Lastly, a main source of systematic error which we have not discussed is the continuum placement. The mean transmissions at various redshifts given by Becker et al. are all obtained by extrapolating the continuum from the red side of Ly$`\alpha `$ by assuming a power law of $`\nu ^{0.5}`$. The continuum likely fluctuates from one quasar to another, and therefore, it would be very useful to apply exactly the same procedure to quasars at lower redshifts where the continuum on the blue side can be more reliably reconstructed. This will tell us how much scatter (and possibly systematic bias) the continuum placement procedure introduces to the measured mean transmission. This kind of error is especially important to quantify given the limited number of quasars available for high redshift measurements at the moment.

AL and LH are supported in part by the Outstanding Junior Investigator Award from the DOE, and AST-0098437 grant from the NSF. MZ is supported by NSF grants AST 0098606 and PHY 0116590 and by the David and Lucille Packard Foundation Fellowship for Science and Engineering. We thank Nick Gnedin for the use of an HPM code and Rupert Croft for providing us with a revised draft of his paper.