text:
# Anisotropies in the redshift-space correlations of galaxy groups and clusters I: Simulated catalogues

## 1 Introduction

Rich clusters of galaxies form at the highest peaks of the mass distribution through physical processes which are simpler than those involved in galaxy formation. As a consequence, models show that the correlation function and power spectrum of clusters of galaxies have the same shape as that of the mass, with an amplitude depending on a small number of parameters (see for example Colberg et al. 2000, Sheth, Mo & Tormen 2001, Padilla & Baugh, 2002). In fact, the space density of clusters of galaxies is usually the only parameter needed to obtain a direct estimate of the underlying mass distribution (Padilla & Baugh, 2002). This shows the clear advantage of studying clusters of galaxies, which has made them popular tracers of the large scale structure of the Universe in the literature.

Cluster distances are derived from member galaxy redshifts, which are affected by galaxy peculiar velocities. Even if all the member galaxy redshifts are used in the determination of the cluster distance, this will be still affected by the cluster peculiar motion, and will not be a true distance measurement. The distribution function of cluster peculiar motions depends upon the value of the mass density parameter $`\mathrm{\Omega }`$, for models with similar mass fluctuation amplitudes (Croft & Efstathiou 1994a; Bahcall, Cen & Gramann 1994). These motions, either for galaxies or clusters of galaxies, produce an apparent distortion of the clustering pattern as measured by the two-point correlation function in redshift space, $`\xi (\sigma ,\pi )`$, where $`\sigma `$ and $`\pi `$ are the separations perpendicular and parallel to the line of sight, respectively. At small separations, non-linear, nearly virialized regions produce elongations along the line of sight which allow measurements of the one-dimensional pairwise rms velocity dispersion, $`<w_{12}^2>^{1/2}`$ (Davis & Peebles 1983). Larger scales are dominated by the infall onto overdense regions in the form of bulk motions which result in a compression of the $`\xi `$ contours along the direction of the line of sight.

Evidence of strong anisotropies in the correlation function of Abell clusters has been given by several authors (see for instance Postman, Huchra & Geller 1992 and references therein). Systematic effects in the Abell catalogue originate in the superposition of clusters along the line of sight and generate a bias in the observed correlation function of the catalogue (Sutherland 1988; Sutherland & Efstathiou 1991; see also Lucey 1983). Using mock catalogues from numerical simulations, van Haarlem et al. (1997) showed that a large fraction of Abell clusters of richness class $`R1`$ would not be physically bound systems and would suffer from contamination by galaxies and groups along the line of sight. Therefore it is likely that the observed large anisotropies could be mainly produced by systematics in cluster detection algorithms from angular galaxy catalogues.

The construction of objectively defined cluster catalogues drawn from machine-scanned survey plates with better calibrated photometry brought a new generation of cluster catalogues which are believed to be less affected by identification biases (APM: Dalton et al. 1992, 1994, 1997; Cosmos: Lumsden et al. 1992). The typical radius used to define clusters in the machine based catalogues is significantly smaller than that used by Abell, reducing the enhancement of cluster richness by projection effects. The clustering signal found in these more recent cluster redshift surveys does not display large enhancements along the line of sight; furthermore, the trend of increasing correlation amplitude with decreasing space density of clusters is weaker than that found for Abell clusters (Croft et al. 1997) and is similar to the trend expected in current models of structure formation. A similarly weak dependence of correlation length on cluster space density is found in X-ray selected cluster catalogues, such as the X-ray Bright Abell Cluster Sample (XBACS) and the REFLEX sample (Collins et al. 2000), which are less susceptible to line of sight projection effects than the optically selected Abell catalogue (Abadi, Lambas & Muriel 1998; Borgani et al. 1999). However, Miller et al. (1999) argue that the early redshift surveys of Abell clusters contain large fractions of low richness clusters (Abell richness class $`R=0`$), which were not intended to form complete samples suitable for statistical analyses. Miller et al. (1999) present the clustering analysis of a new redshift survey of Abell clusters with richness $`R1`$, and with the majority of cluster positions determined using several galaxy redshifts. The clustering signal along the line of sight is greatly reduced in the new redshift surveys compared with the Bahcall & Soneira (1983) results, and is comparable to the amount of distortion of the clustering pattern found for APM clusters (see Fig 5 of Miller et al. 1999). The anisotropy is further reduced after the orientation of two superclusters that are elongated along the line of sight is changed. Peacock & West (1992) also found that restricting the attention to higher richness Abell clusters removed the strong radial anisotropy seen in the clustering measured in the earlier surveys.

Another possible source of systematics could rely on the fact that several cluster distances are determined by a single galaxy redshift, usually the brightest cluster member. Therefore it is important to explore the effects of cluster distance errors on the analysis of correlation function of rich clusters.

The gravitational amplification of small primordial fluctuations has been analysed through numerical simulations to explore the spatial distribution of clusters (e.g. White et al. 1987; Bahcall & Cen 1992; Croft & Efstathiou 1994b; Watanabe, Matsubara & Suto 1994; Eke et al. 1996). These early studies do not reach a consensus on the predicted clustering of clusters in cold dark matter cosmologies. Part of the reason for this discrepancy is due to differences in the way in which clusters are identified in the simulations (Eke et al. 1996). More recent studies have made use of much larger volumes than in these earlier studies, with sufficient resolution to allow the reliable extraction of massive dark matter haloes that can be identified as rich clusters (Governato et al. 1999; Colberg et al. 2000, Jenkins et al. 2001).

In this paper, we analyse the redshift space clustering of massive dark matter haloes in mock catalogues extracted from the $`\tau `$CDM Hubble Volume simulation. In particular, we use the mock catalogues made available by the Virgo Consortium, which contain galaxy angular positions and distances, with a redshift distribution set by a selection function. In this work we will make use of these mock catalogues to analyse in detail the effects of projection biases on catalogues of clusters identified in two dimensions. We also predict the results of the new generations of cluster catalogues with a high degree of spectroscopic completeness such as the 2dF Galaxy Redshift Survey (Norberg et al. 2002).

The outline of this paper is as follows: section 2 describes the equations which allow us to obtain a statistical description of the redshift-space correlation function as a function of coordinates parallel and perpendicular to the line of sight, through the pairwise velocity dispersions, redshift-space correlation length and bias parameters; section 3 presents the correlation function measurement technique. Section 4 presents the different mock cluster samples constructed using several cluster identification algorithms, and studies the results arising from using different number of redshifts to derive cluster distances and also different values of the search radius in cluster identification. In this section, we also study the dependence of the relative pairwise velocities, redshift-space correlation length and bias factors as a function of the number of redshifts used in the determination of the cluster distances. In section 5 we discuss our results and present the main conclusions drawn from this work.

## 2 Redshift-space correlation function anisotropy

The redshift-space correlation function can be calculated as a function of the pair separation parallel and perpendicular to the line of sight. This approach has led to the quantification of characteristics of the redshift-space distribution of galaxy and clusters of galaxies, such as the ”fingers of God”, which are elongated structures seen in redshift surveys, originated from the random motions of galaxies inside clusters.

The strength of the “fingers of god” effect depends on the pairwise velocities of galaxies, which so far has been measured using two different methods. The first method corresponds to that adopted by Loveday et al. (1996) and Ratcliffe et al. (1998), which compares a theoretical expression for the correlation function with the measured values, $`\xi ^m(\sigma ,\pi )`$ (the index $`m`$ indicates measured quantities), in a grid of values of distances parallel and perpendicular to the line of sight. The other approach presented by Padilla et al. (2001), compared the isopleths (curves of equal correlation function amplitude) in the measured and predicted correlation functions. The $`\xi (\sigma ,\pi )`$ contour levels are approximated by the functions $`r^m(\theta )`$ and $`r^p(\theta )`$ for the measured and predicted correlations respectively. Here $`\theta `$ is the polar angle measured from the direction perpendicular to the line of sight (Padilla et al. 2001), such that $`\xi ^x(\sigma _r,\pi _r)\xi ^{fix}=0`$, where $`\sigma _r=r_{\xi ^{fix}}^x(\theta )\mathrm{cos}(\theta )`$, $`\pi _r=r_{\xi ^{fix}}^x(\theta )\mathrm{sin}(\theta )`$, and the index $`x`$ indicates either measured ($`m`$) or predicted quantities ($`p`$).

In either method, the measured correlation function $`\xi ^m(\sigma ,\pi )`$ is compared with the convolution of the real-space correlation function, $`\xi (r)`$ with the pairwise velocity distribution function, $`f(w)`$, following Bean et al. (1983). We calculate

$$1+\xi ^p(\sigma ,\pi )=_{\mathrm{}}^{\mathrm{}}[1+\xi (r)]f[w^{}w^s(r,r^{})]dw^{},$$
(1)
where $`r^2=r^2+\sigma ^2`$, and $`H_0`$ is the Hubble constant, $`r^{}=\pi w^{}/H_0`$ (the prime denotes the line-of-sight component of a vector quantity) and

$$w^s(r,r^{})H_0\beta \xi (r)(1+\xi (r))^1r^{},$$
(2)
is the mean streaming velocity of galaxies at separation $`r`$.

Following usual procedures we calculate the best-fit rms peculiar velocity, $`<w^2>^{1/2}`$, for an exponential distribution,

$$f(w)=\frac{1}{\sqrt{2}<w^2>^{1/2}}\mathrm{exp}\left(\sqrt{2}\frac{|w|}{<w^2>^{1/2}}\right).$$
(3)
We adopted this pairwise velocity distribution as it has shown to be the most accurate fit to the results from numerical simulations (Ratcliffe et al. 1998).

Two estimates of the real-space correlation function are usually adopted in equation 1: the inversion of the projected correlation function (Baugh 1996, Ratcliffe et al. 1998), or a simple power-law fit obtained from the angular correlation function (Padilla et al. 2001). In this work we use the theoretical correlation function of the simulated mass distribution and adopt a linear bias parameter $`b`$ that relates clusters and mass by Fourier transforming the model power spectrum

$$\xi ^{CDM}(r)=\frac{1}{2\pi ^2}b^2_0^{\mathrm{}}P(k)\frac{\mathrm{sin}(kr)}{kr}k^2dk,$$
(4)
We use this to evaluate the theoretical prediction for $`\xi ^p(\sigma ,\pi )`$ using equation (1). We search for the optimum values of the scale independent bias parameter $`b`$ and $`<w^2>^{1/2}`$ by minimising the quantity $`\chi ^2`$,

$$\chi ^2=\underset{i}{}[r_l^m(\theta _i)r_l^p(\theta _i)]^2,$$
(5)
where we have chosen to compare a set of discrete levels, $`l=0.6,0.8,1.0,1.2`$ and $`1.4`$, of the redshift-space correlation function $`\xi (\sigma ,\pi )`$ amplitude instead of comparing values of the correlation function on a grid of $`\sigma `$ and $`\pi `$ distances. Our choice is based on the fact that more reliable and stable results are obtained using this technique (Padilla et al. 2001).

## 3 Measuring $`\xi (\sigma ,\pi )`$

We have computed $`\xi ^m(\sigma ,\pi )`$ as a function of the separation perpendicular ($`\sigma `$) and parallel ($`\pi `$) to the line of sight. To compute $`\xi ^m(\sigma ,\pi )`$, we generate a random catalogue with the same angular limits and radial selection function as the sample of objects. We cross correlate data-data and random-random pairs ($`N_{dd}`$ and $`N_{rr}`$ respectively) binning them as a function of separation in the two variables $`\sigma `$ and $`\pi `$. Our estimate of $`\xi ^o(\sigma ,\pi )`$ is (Davis & Peebles 1983):

$$\xi ^m(\sigma ,\pi )=\frac{N_{dd}n_R^2}{N_{rr}n_D^2}1$$
(6)
where $`n_D`$ and $`n_R`$ are the number of data and random points respectively.

This estimator is affected by uncertainties in the mean density, in particular on large scales where $`\xi (r)`$ is small (eg. Hamilton 1993). However, since our analysis is confined to small separations where the correlation function has a large amplitude, our results are insensitive to the choice of estimator.

Figure 1 shows the contours of equal amplitude of the 2-point correlation function of groups identified in a numerical simulation (left panel) and simulation particles (right panel) in the coordinates $`\sigma `$ and $`\pi `$ (mock galaxy and group samples from Padilla et al. 2001). The compression of the iso-correlation curves for groups in the $`\pi `$ direction can be appreciated, indicating a lack of high pairwise velocities in this sample. From the comparison with results derived from clusters identified in two-dimensional surveys (see for instance Bahcall, Soneira & Burgett, 1986) the existence of either large pairwise velocities or large projection biases is clear, as also discussed by Sutherland (1988).

This elongation could originate in the systematic presence of groups along the line of sight, in the fields of clusters identified from angular data. From the theoretical point of view, such strong elongations along the line of sight are not expected in a hierarchical scenario of structure formation. This is confirmed by the compression observed in left hand side panel of figure 1 and provides a clear evidence for the infall of groups onto larger structures (see also Peacock et al. 2001, for galaxies). As the contours are well defined for all the correlation function levels shown in this figure ($`\xi =0.6`$, $`0.8`$, $`1.0`$, $`1.2`$ and $`1.6`$), the analysis developed in the next sections will be based on the results at these contours. In the right hand side panel, the large distortion seen in the iso-correlation curves in the $`\pi `$ direction is due to the larger peculiar velocities of the simulation particles and their smaller correlation amplitude, associated with the non-linear evolution of the density field on small scales.

## 4 Results from mock catalogues extracted from N-body simulations

In this section, we analyse the statistical properties of cluster samples drawn from mock catalogues. These mock catalogues have been obtained from one of the mock galaxy catalogues extracted from the $`\tau `$CDM Hubble Volume simulation by the Durham Extragalactic Group (Evrard et al., 2002) following a procedure similar to that used in Cole et al. (1998). This simulation follows the evolution of cold dark matter (CDM) density fluctuations in a $`\tau `$CDM cosmology, with parameters $`\mathrm{\Omega }_0=1.0`$, a cosmological constant $`\mathrm{\Lambda }_0c^2/(3H_0^2)=0`$, and a power spectrum described by a shape parameter of $`\mathrm{\Gamma }=0.21`$ and $`\sigma _8=0.6`$. The huge volume of the simulation ($`8h^3\mathrm{Gpc}^3`$) and the large number of particles employed ($`10^9`$) allow cluster statistics to be studied with unprecedented accuracy (Colberg et al. 2000; Jenkins et al. 2001).

The particular mock catalogue we chose to apply the different identification methods to, is that of the APM survey, with a selection function consistent with a limiting magnitude $`b_J=20.0`$.

This particular mock catalogue has the advantage that the derived galaxy apparent magnitudes have $`k+e`$ corrections which reproduce the results of the 2dFGRS (see Norberg et al. 2002, for details). Also we notice that the observer is not chosen at random but instead constrained to lie in a region with similar properties to those of the Local Group.

### 4.1 Mock Groups identification

We use three main algorithms to search for groups in the numerical simulations, and construct four samples of clusters from them:

* Sample 1, is constructed using a friends-of-friends algorithm applied to the 3-dimensional distribution of particles in the simulation cube. We used a linking length $`b=0.2`$, expressed in terms of the mean interparticle separation, adequate for the $`\tau `$CDM cosmology (Helly et al., 2003). This sample is almost completely free of spurious clusters since we have chosen to study only those groups with at least $`60`$ member particles. The fraction of unbound associations found using FOF drops quickly as the number of member particles rises over $`10`$.

Given the lack of a significant number of spurious clusters, the measurements of anisotropies in the correlation function should clearly show the expected infall pattern. As we do not have computational access to the full Hubble Volume simulations, we applied the FOF identification algorithm to a numerical simulation with box size $`600`$h-1Mpc with exactly the same cosmological parameters as the $`\tau `$CDM simulation. This ensures that this particular sample of clusters is statistically comparable to those obtained from the mock galaxy catalogue produced by the Durham group.
* Sample 2 comprises the results from the identification of groups in a mock catalogue that includes a selection function, and where the position of particles is determined from redshifts, that is including particle peculiar velocities. We use a friends-of-friends algorithm with different linking lengths in the directions parallel and perpendicular to the line of sight. This procedure emulates the identification process used in the construction of the group sample obtained from the Updated Zwicky Catalogue, which is based on an algorithm described by Huchra & Geller (1982). In this case we use different linking lengths in the directions perpendicular and parallel to the line of sight, which vary linearly with distance, $`d_0=V/V_{fid}0.229`$h-1Mpc and $`v_0=V/V_{fid}350`$km/s, where $`V`$ is the distance, and $`V_{fid}=7000`$km/s. These values correspond to a density contrast $`\delta \rho /\rho 80`$.
* Sample 3 is obtained using the algorithm presented by Lumsden et al.(1992), which is used in identifying clusters from the COSMOS catalogue. Following the prescription described by Lumsden et al. (1992), we use an angular search radius $`r_c=1.0`$h-1Mpc, and apply this procedure to the same mock catalogues as the second algorithm. We tested this identification algorithm by applying it to the COSMOS galaxies, and comparing the outcome of the identification algorithm to the positions of the EDCC clusters (Lumsden et al. 1992). We found out that $`75\%`$ of the EDCC clusters were re-obtained by our identification procedure. We also applied this to a numerical simulation, and found that only $`70\%`$ of the clusters identified in the simulation cube with a FOF algorithm were found by our angular identifier. More importantly, only $`60\%`$ of the clusters identified are real clusters. These numbers illustrate the degree of projection effects sample 3 is subject to, and its study will serve as a test of how these effects translate into the cluster correlation function in redshift-space.
* Sample 4: We extract from sample 3 the subset of clusters whose angular positions are coincident (within the cluster identification radius) to those in sample 2. This sample is derived from a 2 dimensional catalogue and is confirmed to have a redshift space identification at the same angular position. Cluster distances are determined by using the galaxy members identified from the angular position data. This sample should be less affected by projection effects than that obtained from the angular data alone. However, cluster distances may include particles that are not physically related to any particular cluster at all. We find that only $`50\%`$ of the clusters from sample 3 are present in this sample indicating that this procedure provides a significant restriction. We also notice that this fraction is slightly different than that resulting from comparing the sample 3 clusters with those obtained from using the FOF algorithm on a full simulation box.

All the samples defined in this section include groups with a minimum mass threshold. This threshold is set so as to make the space density of all the samples as close to that of sample 1 as possible.

In figure 2. we show the different results from the identification of clusters of galaxies in the mock catalogue using angular and 3-D galaxy positions. This figure shows the density of galaxies in the mock catalogue in pixels, and clusters of galaxies in circles. The gray scale of the pixels indicate number of galaxies; the darker the pixel the higher is the number of galaxies in it. As can be seen, the circles enclose regions of high density of galaxies. The different line styles with which the circles are drawn correspond to the different algorithms used in the identification of clusters. The double-line circles correspond to clusters from sample 3, that is, clusters identified using angular data. In this case, the radius of the circle corresponds to the radius out to which the counting of galaxies was done when measuring the cluster richness. The thick line circles show clusters from sample 2, identified from 3-dimensional positions affected by the peculiar velocities of the clusters. In this case, the radius of the circle is proportional to the richness of the clusters. By inspection to this figure, we can appreciate the large fraction of clusters from sample 3 whose angular positions are not coincident with any 3-dimensional cluster within one angular search radius.

The aim of constructing the mock cluster samples presented here is to understand and reproduce the results obtained from observational cluster samples. In particular, we intend to establish the reason for the large elongations along the line of sight present in the measured correlation functions. The different behaviours of the correlation functions found for the four mock samples are shown in figure 3, where we used clusters of galaxies with distances measured using $`10`$ member galaxies. The results from samples 1 and 2 show the expected infall pattern. This simply indicates that the FOF algorithm still works well when applied to a catalogue which incorporates a selection function and includes the effects of galaxy peculiar velocities. The infall is not so evident when inspecting the remaining panels, which simply show the different degrees of projection effects the different samples are subject to. It can be seen that samples 3 and 4 show somewhat large elongations along the line of sight. Taking into account that all panels show the results for the same value of $`n_z=10`$, this is a clear signature of projection effects. The fact that sample 4 shows elongated contours along the line of sight is indicative of substantial projection contamination. Given that we have constructed this sample by requiring an angular coincidence between 2-dim and 3-dim identified clusters, such anisotropies reflect the fact that there is a strong contribution by other structures along the line of sight in producing the observed distortion pattern.

In the following section we provide detailed descriptions of the characteristics of each sample, analysing the differences in elongations arising from using different $`n_z`$.

### 4.2 Effects of observed number of redshifts

In order to make a thorough comparison between the mock cluster samples and the observational correlation functions, we study in detail the effects of using different $`n_z`$ in each of our mock samples.

We first show the results from sample 1 which comprises clusters of galaxies identified from the 3-dimensional distribution of particles in the full numerical simulation. Figure 4 shows the correlation function contours when the distances to the clusters in this sample have been calculated using $`1`$ (left) and $`20`$ (right panel) member redshifts. As can be seen, the differences between the panels are quite small. The left panel only shows slight evidence of an elongation along the line of sight. This is an effect produced by the errors in the distance to these clusters arising from the use of only one redshift to determine their distances. By comparison between the left upper panel in figure 3 and the right panel in figure 4 (that is, Sample 1 with $`n_z=10`$ and $`20`$), it can be seen that the use of a small $`n_z=10`$ is sufficient to erase spurious elongations along the line of sight, and that considering more redshifts per cluster makes no significant difference. This result once more reflects the reliability of the FOF algorithm in finding clusters of galaxies in numerical simulations.

We now investigate the effects of changing $`n_z`$ in sample 2, which consists of clusters identified from a mock galaxy catalogue which incorporates a selection function. We show in figure 5 the correlation function contours for samples with distances obtained using $`1`$ and $`20`$ redshifts in the left and right hand side panels respectively. As can be seen in this figure, the effects of using a small $`n_z`$ are not critical. In the case of using $`1`$ galaxy, the infall pattern can still be seen, specially in the correlation contours corresponding to $`\xi (\sigma ,\pi )<1`$. The inclusion of further redshifts only makes a small change, although it can be seen that in this case, this signature is also visible for levels of higher correlation function. When comparing with figure 4 one can notice a neat infall pattern in the clusters from sample 1 for high levels of correlation even when $`n_z=20`$ (a very clear signature for $`\xi =1`$ contour levels, and mildly visible for $`\xi =1.4`$). The case of sample 2 is less ideal, since this infall pattern can only be marginally seen for a correlation function contour level corresponding to $`\xi =1`$, which simply reflects the higher degree of difficulty in identifying clusters from a galaxy catalogue which incorporates a selection function.

It should be noted that sample 2 is representative of cluster or group catalogues identified in galaxy redshift surveys, as is the case with UZC groups. Therefore, the infall of structures which is a signature of the hierarchical clustering obtained in the correlation function contours measured using the UZC groups sample (Padilla et al. 2001), can be assessed by the results of sample 2. In other words, the study arising from sample 2 results in a correlation function which is not badly affected by projection effects, even for small $`n_z`$, which also indicates that the use of a small $`n_z`$ does not induce important elongations along the line of sight.

We next consider a mock cluster catalogue drawn from an angular galaxy catalogue (sample 3). Here the distances to each cluster were obtained using the redshifts of galaxies identified as cluster members using angular data. The correlation function of these catalogues are shown in figure 6, where we show the results when using $`1`$, $`5`$, $`10`$, and $`20`$ galaxy redshifts in determining cluster distances in the upper left, upper right, lower left and lower right panels respectively.

The results in this figure show significant differences in the elongations for the different samples. In the case where only $`1`$ galaxy redshift is used to assess the cluster distance, the elongation along the line of sight is severe, and not entirely different to what was found for Abell clusters with distances estimated using less than $`10`$ galaxy redshifts (see for instance Bahcall, Soneira & Burgett, 1986). The use of $`5`$ galaxies already makes a marked difference but would still not be sufficient for an infall pattern to be found. For comparison, the number of Abell clusters with $`n_z5`$ is small, about 20 % of the total sample. The first, faint traces of an infall can be seen when considering $`10`$ galaxy redshift measurements, and only for small values of $`\xi (\sigma ,\pi )<1`$, which are difficult to obtain from an observational sample without large contributions from noise.

The main sources of elongation cannot be considered to be the error in the cluster position arising from using a small $`n_z`$. This is so since, for instance, we have already shown in figure 5 that the use of a very small number of galaxies is still not enough to erase the infall pattern. It is more likely that the elongation is produced by the inclusion of non member galaxies in the determination of cluster distances, the worst scenario being that in which only a fraction of the identified clusters correspond to physically bound galaxies, and the remainder being constructions of projection effects that cause the spreading of galaxies along the line of sight to appear as coherent structures, as only angular data is available.

We now consider what are the expected effects of the number of redshift measurements in samples of clusters identified in two dimensions with angular positions coincident with 3-dimensional identification.

Figure 7 shows the correlation function contours for clusters in mock sample 4. The different panels correspond to different $`n_z`$, in the same order as the previous figure.

As can be seen in figure 7, the projection effects in sample 3 are still present, even after the clusters have coincident angular positions with a cluster sample identified in 3D. The elongations in $`\xi (\sigma ,\pi )`$ are quite severe, and reflect the contamination of foreground and background structures, which introduces biases in the cluster distances.

### 4.3 Effects of projected radius in the identification algorithm

The algorithms used to identify clusters from angular data depend quite sensitively on the search radius $`r_c`$ used to find member candidates. In the case of Abell clusters this corresponds to $`r_c=1.5`$h-1Mpc.

A possible cause of projection effects becomes clear when one considers that the density of galaxies diminishes as a function of cluster radius. This means that in the outskirts of the cluster, the angular population of galaxies will be composed of ever larger fractions of foreground and background galaxies. Therefore it is increasingly probable to consider as member galaxies, objects that are not physically bound to the cluster. We investigate this effect by using different cluster search radius in our 2-dimensional cluster finding algorithms. In figure 8, we show the correlation function contours for clusters identified using $`r_c=0.5`$, $`1.0`$ and $`1.5`$h-1Mpc in the left, middle and right panels respectively. In all these cases, $`n_z=10`$, which ensures us that any elongations observed in the $`\xi (\sigma ,\pi )`$ isopleths will be mainly produced by projection effects arising from the choice of radius.

As can be seen, the use of a small search radius is not entirely convenient as it produces spurious clusters, probably from confusing satellite haloes with proper clusters. This is in agreement with the elongated iso-correlation contours in the left panel, and also with the small amplitude of the correlation length $`\sigma _0`$, that is $`\xi (\sigma _0,\pi =0)=1`$, and is supported by figure 9, which shows the clusters of galaxies identified using $`r_c=1.0`$h-1Mpc (circles), and those identified using $`r_c=0.5`$h-1Mpc (dots); the gray scale of the pixels indicate number of galaxies, and the radii of the circles correspond to the radius out to which the counting of galaxies was done when finding the cluster richness. As can be seen in three cases in figure 9, the identification using a small search radius produces pairs of clusters whereas the use of a larger search radius produces a single cluster. The satellite cluster identified with a small search radius corresponds to a lower peak in the density of projected galaxies in all cases; also, the richness of a satellite is smaller than the main cluster in general, confirming the hypothesis of low mass satellites stated earlier in this paragraph.

The other values considered, $`r_c=1.0`$ and $`1.5`$h-1Mpc, produce similar results, both largely free from projection effects, and with a larger correlation length $`\sigma _0`$ (figure 8). There is a slight difference in the noise level in these samples, the correlation levels in the $`r_c=1.5`$h-1Mpc sample being less smooth than the smaller $`r_c`$ case. This is not a severe difference, but is enough to support our choice of $`r_c=1.0`$h-1Mpc used throughout this paper.

### 4.4 Relative velocities and correlation lengths in mock cluster samples

In this section we present determinations of relative velocities, redshift-space correlation lengths and bias factor from the $`\xi (\sigma ,\pi )`$ anisotropies. Specifically, we study the dependence of the relative pairwise peculiar velocities obtained using eq. 5 for the different mock samples as a function of $`n_z`$. We recall that the cosmological parameters of the simulation from which the mock catalogues were extracted ($`\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.0`$, $`\mathrm{\Gamma }=0.21`$ and $`\sigma _8=0.6`$) are those used in the calculations.

Figure 10 shows the pairwise velocities found in the different samples (sample 1 in upper left panel, sample 2 in the upper right, sample 3 in the lower left, and sample 4 in the lower right panel). The different thin lines correspond to results using different values of $`\xi (\sigma ,\pi )`$. The thick solid line shows the mean pairwise velocity dispersion from the clusters identified in three dimensions in the simulation box.

This plot shows the great range of results which can be obtained from the contours in the correlation function if the mean $`n_z`$ used for the clusters in the sample is not known. It is remarkable the fact that the relative velocities found using clusters from sample 2 which were obtained in redshift space are in good agreement with the actual value measured from the full simulation.

The disagreement between the value of true relative cluster velocities with the measurements from the redshift-space correlations is most severe when the sample has been identified using angular data. It is also interesting that the velocities obtained from sample 4 (clusters identified from 2 dimensions with angular confirmation from inspection over the 3-dimensionally identified sample) are in some cases still very high even when their distances are obtained using a large $`n_z`$. It should be noted though, that the velocities favoured by low noise contour levels, $`\xi 1`$, suggest small cluster relative velocities. The problem of a wide range of predictions for $`w_{12}`$ is not present in the case of clusters in sample 3, in part due to the larger number (almost a factor of 2) of clusters in this sample. The results obtained from different correlation levels are more compatible, showing values differing by less than $`1000`$ km/s among the different correlation levels. The observed behaviour of the relative velocities does not favour sample 4 over sample 3, indicating that the confirmation by angular coincidence with real clusters can not remedy the systematic problems of samples of clusters identified in two dimensions.

This problem can be dug into a little further by studying the cluster redshift-space correlation lengths, which also show great variation with $`n_z`$. We show these quantities in figure 11, where marked differences between samples of clusters identified from angular and 3-D data can be seen. Both samples obtained using 3-D information show a steady value of $`s_0`$, without any significant variation with $`n_z`$. The first important result from this figure is that the results from clusters in sample 2 are in excellent agreement with the true values. We notice that the results for sample 1 are indistinguishable from those of sample 2 and are not shown for clarity.

An interesting feature in figure 11, is the difference between confirmed and unconfirmed 2-D clusters, which also show controversial results in figure 10. The observed correlation length of clusters in sample 4 is higher than that of redshift space identified clusters in sample 2. Furthermore, this correlation length is also higher than that of the clusters in sample 3, which is expected to be more contaminated by foreground/background galaxies. This differences could be explained if several spurious clusters in sample 3, were associations of a small number of galaxies along the line of sight. The higher correlation amplitude derived for sample 4 clusters with respect to those in sample 2, could indicate that the identification of clusters using angular data is biased towards massive systems. This follows from the larger values of $`s_0`$ found from the sample of confirmed clusters. The inclusion of spurious small associations of galaxies in sample 3, would make the value of $`s_0`$ in this case smaller than those obtained for clusters in sample 2. We have also confirmed this by inspection to the mean cluster mass of samples 2 and 4.

These findings are in agreement with the results of Miller et al. (1999), who showed that by restricting a sample of Abell clusters to high richnesses the degree of projection effects is significantly reduced. The difference in correlation lengths also explains the marginally larger cluster relative velocities found in sample 4 compared with sample 3, at least when considering high levels of correlation, where the noise is not playing an important role in the results from sample 4. Equation 2 shows that the streaming velocities of clusters are proportional to their correlation functions. This means that the clusters in sample 4 have larger streaming velocities than those in sample 3. Taking into account that the elongations along the line of sight seen in figure 7 are similar for both samples, it is clear that the results from equation 5 will yield larger values of cluster relative velocities for the sample 4 than for sample 3.

The results for the values of the bias parameter obtained from minimising equation 5 as a function of $`n_z`$ show a behaviour similar to the results from the redshift-space correlation length. Figure 12 shows the resulting bias parameter for the different mock samples (the symbols are indicated in the figure). The dotted lines show the acceptable range of effective bias parameters from the simulation, obtained using the Sheth, Mo & Tormen (2001) mass function, and the number density of clusters in the mock samples along with its uncertainty.

As it can be seen, the bias found for the clusters in mock sample 2 is within the acceptable values for the CDM effective bias. This result holds for any value of $`n_z`$. As in the previous figure, the results for sample 1 are consistent with the true values, and are not shown in order to preserve clarity. The result from sample 3 is in agreement for roughly $`n_z>10`$. Sample 4 shows significantly larger values of bias and correlation length, suggesting that this sample comprises a high mass cluster population in agreement with the results found earlier in this section.

## 5 Conclusions

We have analysed the consequences of different effects on the distortion of the correlation function of clusters in redshift space for different mock samples of galaxy clusters. We take into account the pairwise velocity distribution of galaxy systems, coherent bulk motions, as well as redshift errors, and different cluster identification systematics. The mock cluster samples were extracted from numerical simulations, following algorithms that closely match the procedures used to identify the clusters in observational samples. We find that the correlation function of galaxy systems is influenced by wrongly assigned cluster distances due to a small number of galaxy redshift measurements in the fields of clusters. However, the most important effect is that associated with the cluster identification procedure from two dimensional surveys. Due to these effects, the estimated mean pairwise velocity dispersion can be an order of magnitude larger than the actual value.

This is consistent with the observed anisotropy of the correlation function contours of UZC groups, which shows the flattening produced by this infall motion (Padilla et al. 2001). For samples of clusters identified in mock galaxy redshift surveys we find an infalling pattern even in the case where the distances to groups or clusters of galaxies are obtained using a small number of member redshifts.

In order to provide predictions for new generation cluster catalogues we have also analysed a mock sample of groups resembling those obtained by Merchán & Zandivarez (2002) from the 2dF100K release galaxy redshift survey. This mock catalogue is subject to the same complicated angular mask as the observations and serves as a further test of the reliability of the information that can be obtained from observational cluster and group catalogues. The results shown in figure 13 are similar to those of Sample 2, and indicate that the angular mask does not affect severely the results so that consistent values of $`w_{12}`$ and $`s_0`$ can be derived.

In order to explore the projection effects afflicting cluster samples identified in two dimensions we analysed different mock cluster samples identified from mock angular catalogues using different cluster search radius. Our findings favour the use of smaller radius than that used in the identification of Abell clusters. The use of a too small value produces a sample increasingly affected by projection effects and the inclusion of smaller groups of galaxies which are probably cluster satellite haloes.

We notice that 2-dim identified clusters whose angular positions coincide with a true 3-dim cluster show a strongly distorted clustering pattern. This fact shows the large influence of groups and clusters along the line of sight in the observed elongation. Contamination is therefore mainly arising from structures present in the filamentary large scale distribution. These results indicate the difficulties of obtaining unbiased cluster samples even when using accurate photometric data such as present in new and forthcoming surveys. Redshift information is crucial for constructing cluster samples which lack significant projection effects, since galaxy colours alone can remove foreground and background contamination but not distinguish membership to clusters at relatively smaller separations along the line of sight.

## Acknowledgments

This work was supported in part by CONICET, Argentina, and the PPARC rolling grant at the University of Durham. DGL acknowledges support from the John Simon Guggenheim Memorial Foundation. We thank the Referee for invaluable comments and advice which greatly improved the previous version of the paper. We have benefited from helpful discussions with Carlton Baugh. We acknowledge the Durham Extragalactic Astronomy Group and the Virgo Consortium for making the Hubble Volume simulation mock catalogues available.


format: text

text:
# The Thermal Memory of Reionization History

## 1. Introduction

The detection by the Wilkinson Microwave Anisotropy Probe (WMAP) of a large optical depth $`\tau _e`$ to electron scattering has opened a new window in studies of the ultra–high redshift ($`z15`$) universe (Kogut et al. 2003, Spergel et al. 2003). Taking the value $`\tau _e=0.17\pm 0.04`$, inferred from the polarization and temperature anisotropies of the cosmic microwave background (CMB), at face value, implies the reionization of the universe must begin very early.

The optical depth to electron scattering is given by (e.g. Dodelson 2003)

$`\tau _e={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{(1+z)H(z)}}\sigma _Tn_e(z)`$ (1)
where $`H(z)`$ is the Hubble parameter at redshift $`z`$, $`\sigma _T`$ is the Thompson cross-section, and $`n_e(z)`$ is the proper free electron density. This can be rewritten as

$`\tau _e=0.0691\times \mathrm{\Omega }_bh{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{(1+z)^2dz}{\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}}}`$ (2)
$`[(1Y_P)X_{\mathrm{HII}}+{\displaystyle \frac{1}{4}}Y_P(X_{\mathrm{HeII}}+2X_{\mathrm{HeIII}})]`$
where $`\mathrm{\Omega }_b`$, and $`\mathrm{\Omega }_m`$ are the baryon and matter densities in fraction of the critical, $`h`$ is the Hubble constant in units of $`100`$ km/s/Mpc, and $`Y_P=0.244\pm 0.002`$ is the helium mass fraction (Burles, Nollett & Turner 2001). The fractions of ionized hydrogen $`X_{\mathrm{HII}}`$, singly ionized helium $`X_{\mathrm{HeII}}`$, and doubly ionized helium $`X_{\mathrm{HeIII}}`$, are functions of $`z`$.

Ignoring helium, the observed $`\tau _e=0.17\pm 0.04`$ is consistent with a universe in which $`X_{\mathrm{HII}}`$ changes from essentially zero to close to unity at $`z=17\pm 3`$, and $`X_{\mathrm{HII}}1`$ since (Kogut et al. 2003). If helium is fully ionized together with hydrogen, the reionization redshift changes slightly to $`z=15.3\pm 2.7`$. 111Throughout this paper, we adopt $`\mathrm{\Omega }_bh^2=0.024`$, $`\mathrm{\Omega }_mh^2=0.14`$, and $`h=0.72`$ the central best–fit values measured by WMAP (Spergel et al. 2003).

Two features are noteworthy. First, helium reionization has a sub-dominant effect on $`\tau _e`$ compared to hydrogen. Second, since the electron scattering optical depth is controlled by the free electron density, it is insensitive to the neutral fractions of hydrogen ($`X_{\mathrm{HI}}=1X_{\mathrm{HII}}`$) and helium ($`X_{\mathrm{HeI}}=1X_{\mathrm{HeII}}X_{\mathrm{HeIII}}`$), as long as they are small.

In contrast, the (hydrogen) Lyman-alpha (Ly$`\alpha `$) optical depth, inferred from the spectra of distant quasars, is extremely sensitive to small amounts of neutral hydrogen. The Ly$`\alpha `$ optical depth at mean density is (Gunn & Peterson 1965)

$`\tau _\alpha =41.8\left[{\displaystyle \frac{X_{\mathrm{HI}}}{10^4}}\right]\left[{\displaystyle \frac{1+z}{7}}\right]^3\left[{\displaystyle \frac{H(z=6)}{H(z)}}\right],`$ (3)
Whether $`X_{\mathrm{HI}}`$ is $`10^4`$ or $`10^5`$, for instance, makes a big difference to $`\tau _\alpha `$. Using a model for large scale structure in the intergalactic medium (IGM), the observed mean Ly$`\alpha `$ transmission at $`z6`$ (or the lack thereof i.e. the Gunn-Peterson trough; Becker et al. 2001) implies a ($`1\sigma `$) lower limit on the hydrogen neutral fraction for a fluid element at mean density: $`X_{\mathrm{HI}}>10^4`$. The analog for Ly$`\beta `$ provides a stronger constraint, due to the smaller absorption cross-section: $`X_{\mathrm{HI}}>5\times 10^4`$ (taken from Lidz et al. 2002; see also Cen & McDonald 2002, Fan et al. 2002).222 Note that some authors quote volume– or mass–weighted neutral fractions. We here use the neutral fraction for a fluid element that happens to be at the cosmic mean density, motivated by the fact that we will discuss the temperature evolution of fluid elements with the same property in this paper.

Taken at face value, these two separate observations are therefore consistent with the simplest “vanilla” model in which the universe is reionized in a single step: the neutral fraction experiences a drop of order unity at at $`z\text{ }>15`$, and it remains $`1`$ thereafter. Our goal in this paper is to confront this model with other existing observations, and to see if additional constraints can be put on the reionization history of our universe.

Several authors have pointed out that the evolution of the Ly$`\alpha `$ optical depth suggests reionization might take place not much earlier than $`z6`$ (Becker et al. 2001, Djorgovski 2001, Cen & McDonald 2002, Gnedin 2001, Razoumov et al. 2002; but see also Barkana 2001, Songaila & Cowie 2002). This is based on an extrapolation of the mean transmission measurements from lower redshifts. The small number of lines of sight used (a single quasar was employed for the measurement at $`z6`$; but see Fan et al. 2003 for three new $`z>6`$ sources with Gunn-Peterson troughs), the challenge of sky subtraction and continuum extrapolation in Gunn-Peterson trough measurements (see discussion in Becker et al. 2001; see also Hui et al. 2001), as well as our still maturing understanding of radiative transfer during reionization, motivates us to look for other clues for a late period of reionization.

The main idea is quite simple. Reionization typically heats up the IGM to tens of thousands of degrees, and the gas subsequently cools due to the expansion of the universe as well as due to other processes. If the universe was reionized early, and has stayed highly ionized thereafter, photo–ionization heating of the gas cannot overcome the overall cooling, and the IGM might reach too low a temperature at low redshifts compared to observations. This idea is not new (e.g. Miralda-Escude & Rees 1994, Hui & Gnedin 1997, Haehnelt & Steinmetz 1998, Hui 2000, Theuns et al. 2002). Theuns et al. (2002), in particular, deduced a limit on the reionization redshift ($`z<9`$) based on temperature measurements by Schaye et al. (2000), assuming a quasar-like ionizing spectrum, and that HeII reionization occurs at $`z3`$. Our objective here is to seek a formulation of this argument that is as clean and robust as possible, that reveals clearly the underlying assumptions, and to check the consistency with the IGM temperatures of specific models that produce the high value of $`\tau _e`$ measured by WMAP.

The Ly$`\alpha `$ forest temperature measurements we will use are taken from Zaldarriaga, Hui & Tegmark (2001; ZHT01 thereafter): $`T_0=2.1\pm 0.9\times 10^4`$ K at $`z=2.4`$, $`2.3\pm 0.7\times 10^4`$ K at $`z=3`$, and $`2.2\pm 0.4\times 10^4`$ K at $`z=3.9`$. Note that 2 $`\sigma `$ errorbars are quoted, and the temperature $`T_0`$ was derived for fluid elements at the mean density, consistent with our modeling of the IGM temperature in the rest of this paper. There have been several other measurements in the past (Ricotti, Gnedin & Shull 2000, Schaye et al. 2000 \[ST00\], Bryan & Machacek 2000, McDonald et al. \[MM01\], Meiksin, Bryan & Machacek 2001). The ones that are easiest to compare, because they are based on very similar datasets, are ST00, MM01 and ZHT01. The former two are based on line width measurements, while the last one makes use of the small scale transmission power spectrum. A virtue of the last method is that the temperature constraints come from marginalizing over a wide array of parameters, including the slope and amplitude of the primordial power spectrum (e.g. Hui & Rutledge 1999), and the equation of state index. It is reassuring that MM01 and ZHT01, using very different methods and employing different simulations (MM01 using hydrodynamic simulations, and ZHT01 using N-body simulations with a marginalized smoothing to mimic Jeans smoothing), agree well with each other. The results of ST00 are somewhat discrepant from these two works – the reader is referred to ZHT01 for further discussions.

An important issue in determining the temperature of the IGM is the second ionization of helium (HeII$``$HeIII). As we will see, this can be an important source of heating at low redshifts $`z34`$. There are several lines of evidence that suggest HeII might be reionized at $`z3`$, including observations of HeII patches that do not seem to correlate with HI absorption (Reimers et al. 1997, Anderson et al. 1999, Heap et al. 2000, Kriss et al. 2001, Jakobsen et al. 2003), increase in IGM temperature (Schaye et al. 2000; Theuns et al. 2002), the evolution of the hardness of the ionizing background spectrum (Songaila 1998), and evolution of the mean transmission (Bernardi et al. 2003). On the other hand, it is unclear if the fluctuations in HeII absorption observed in a few lines of sight might not be due to a naturally fluctuating IGM (Miralda-Escude, Haehnelt & Rees 2000); the IGM temperature measurements by MM00 and ZHT01 are consistent with no feature at $`z3`$; power spectrum evolution also seems to argue against HeII reionization at $`z3`$ (McDonald & Seljak, private communication). In this paper, when we consider constraints from the temperature measurement, primarily at $`z=3.9`$, we therefore leave two options open: HeII can be ionized or not ionized by $`z=3.9`$. With this explained, we can now state our vanilla model in more concrete terms. It has two variants: 1. both hydrogen and helium are fully reionized at $`z\text{ }>15`$, and they remain highly ionized thereafter; 2. the same as 1. except that helium is only singly ionized, and remains so until at least past $`z=3.9`$.

The paper is organized as follows. In §2, we explain the idea of a thermal asymptote for the IGM, and use it to derive a constraint on the hardness of the ionizing spectrum if the universe were to reionize before $`z=10`$, and remains highly ionized thereafter. In §3, we discuss limits on the hardness of the ionizing spectra, and show that these spectra fall short of making the thermal asymptote sufficiently hot to match observations. We then work out illustrative examples in §4 of how order unity changes in the neutral fractions at $`z\text{ }<10`$ can reproduce the temperature measurements, while being consistent with the WMAP data. We go on to offer a physically motivated model in §5, and discuss the implications of the fact that different fluid elements are reionized at different times (stochastic reionization). We conclude in §6.

## 2. Thermal Asymptotics

Several thermal processes are at work in a photoionized IGM. They are described in detail in Hui & Gnedin (1997). Here is a brief qualitative summary:

* Adiabatic heating/cooling. Gas elements can heat up or cool simply due to adiabatic contraction or expansion. The overall expansion of the universe drives a temperature fall-off as $`(1+z)^2`$ as $`z`$ decreases.
* Photoionization heating. Photons inject energy into the gas by ionizing hydrogen or helium.
* Recombination cooling. Protons and electrons (or ionized helium and electrons) can cool by recombining and radiating energy away.
* Compton cooling. At sufficiently high redshifts, $`z\text{ }>10`$, Compton scattering of free electrons with the lower temperature CMB can be an important source of cooling.

Typically, photoionization heating provides the dominant source of heating for the tenuous IGM. This occurs in primarily two forms. One is during what we call reionization ’events’ – these are (possibly extended) periods in which a given gas element experiences order unity changes in the fractions of HI/HeI/HeII. The other is photo-heating for an already highly ionized plasma – this occurs through photoionization of small amounts of HI/HeI/HeII, whose abundances are determined by photoionization equilibrium. The former provides a big boost to the temperature, while the latter determines the asymptotic thermal state of the IGM. Typical temperature evolutions are illustrated in several toy models in Fig. 1.

Fig 1 shows the evolution in temperature $`T_0`$ for a fluid element at mean density, with a variety of reionization histories. These histories have different reionization redshifts and initial temperature jumps, selected here only for illustration. The code for computing the thermal and chemical evolution is described in Hui & Gnedin (1997). For each thermal curve, one can see a rise to high temperatures after a reionization event. Thereafter, the varieties of reionization histories and initial reheat temperatures all result in a rather similar late time thermal asymptote (black dotted line). In fact, if reionization occurs before $`z10`$ (and the IGM stays highly ionized since, or to use our earlier terminology: no reionization events take place thereafter), the temperature is always very close to the asymptote for $`z<4`$, which is where temperature measurements exist (points with error-bars). This is true even if reionization is not a single step process, as illustrated by the dashed (magenta) line, as long as complexities in the reionization process occur before $`z10`$. The term ’complexities’ here has a rigorous meaning: it refers to two or more episodes during which the neutral fractions of hydrogen and/or helium undergo changes of order unity. In all thermal curves shown in Fig. 1, with one exception, such changes take place before about $`z10`$. The exception, i.e. the leftmost curve where reionization takes place at $`z7`$, is also the only one where the asymptote is not reached by $`z=4`$ – this is because there is not sufficient time for the IGM to cool after a recent episode of reionization.

Fig. 1 therefore illustrates two very important general facts:

1. The thermal state at $`z4`$ does not remember the part of reionization history prior to $`z10`$.
2. The IGM does, however, retain short-term memory of reionization events in its recent past.

What determines the late time thermal asymptote? It is the combination of photoionization heating, recombination cooling, and adiabatic cooling. An approximate analytic expression can be derived for the thermal evolution of the IGM (Hui & Gnedin 1997). We find that the following simple formula provides a more accurate ($`5\%`$) fit to the numerically computed thermal asymptote (for $`T_0`$ at $`z=24`$), for a variety of spectral shapes we have tested:

$`T_0=[B(1+z)^{0.9}]^{\frac{1}{1.7}}`$ (4)
where

$`B18.8\mathrm{K}^{0.7}\left[\sqrt{{\displaystyle \frac{0.14}{\mathrm{\Omega }_mh^2}}}{\displaystyle \frac{\mathrm{\Omega }_bh^2}{0.024}}\right]`$ (5)
$`\times \left(T_{\mathrm{HI}}+{\displaystyle \frac{X_{\mathrm{HeII}}}{16}}T_{\mathrm{HeI}}+{\displaystyle \frac{5.6X_{\mathrm{HeIII}}}{16}}T_{\mathrm{HeII}}\right)`$
$`T_ik_B^1{\displaystyle \frac{_{\nu _i}^{\mathrm{}}J_\nu \sigma _i(\mathrm{h}\nu \mathrm{h}\nu _i)𝑑\nu /(\mathrm{h}\nu )}{_{\nu _i}^{\mathrm{}}J_\nu \sigma _i𝑑\nu /(\mathrm{h}\nu )}}`$
Here, $`k_B`$ is the Boltzmann constant, the symbol $`J_\nu `$ denotes the ionizing intensity as a function of frequency (it has unit of energy per frequency per time per area per ster-radian), $`\sigma _i`$ denotes the ionization cross-section for the respective species ($`i=\mathrm{HI},\mathrm{HeI}`$ or $`\mathrm{HeII}`$), $`\nu _i`$ is the ionizing threshold frequency, and $`\mathrm{h}`$ is the Planck constant (in distinction from the Hubble constant $`h`$). It is useful to remember $`\mathrm{h}\nu _{\mathrm{HI}}/k_B=1.57807\times 10^5`$ K, $`\nu _{\mathrm{HeI}}=1.808\nu _{\mathrm{HI}}`$, and $`\nu _{\mathrm{HeII}}=4\nu _{\mathrm{HI}}`$.

The asymptote given in eq. (4) assumes that at late times, hydrogen is highly ionized. That is why the term due to photo-heating of HI, $`T_{\mathrm{HI}}`$, has the form shown in eq. (5). The denominator of $`T_{\mathrm{HI}}`$, when multiplied by $`4\pi `$, is the photoionization rate of neutral hydrogen – its inverse is hence proportional to the (small) hydrogen neutral fraction under ionization equilibrium. Its numerator gives the photo-heating rate per neutral atom. The combination of factors in $`T_{\mathrm{HI}}`$ therefore gives a temperature scale whose amplitude is proportional to the net amount of HI photo-heating. Similarly, the terms $`X_{\mathrm{HeII}}T_{\mathrm{HeI}}/16`$, and $`5.6X_{\mathrm{HeIII}}T_{\mathrm{HeII}}/16`$ in eq. (5) quantify the importance of photo-heating of HeI and HeII respectively. The factors of $`X_{\mathrm{HeII}}`$ and $`X_{\mathrm{HeIII}}`$ arise due to photoionization equilibrium. If, asymptotically, helium is doubly (singly) ionized, then $`X_{\mathrm{HeII}}`$ ($`X_{\mathrm{HeIII}}`$) can be set to zero.

It is important to emphasize that the thermal asymptote given in eq. (4) is determined completely by the shape (or ’hardness’) of the ionizing spectrum $`J_\nu `$, but not its amplitude. For a power-law $`J_\nu \nu ^\alpha `$, the thermal asymptote can be recast simply as:

$`T_0=2.49\times 10^4\mathrm{K}\times (2+\alpha )^{\frac{1}{1.7}}\left({\displaystyle \frac{1+z}{4.9}}\right)^{0.53}`$ (6)
assuming, asymptotically, helium is doubly ionized.

Given the above, we can ask the following question: how hard does the ionizing spectrum have to be for the thermal asymptote to match the observed temperatures at low redshifts? It suffices to use the measurement at the highest redshift, $`z=3.9`$: $`T_0=2.2\pm 0.4\times 10^4`$ K (2 $`\sigma `$ error-bar; ZHT01). More concretely, what constraint can be placed on $`\alpha `$ if the thermal asymptote were to reach $`T_0>1.6\times 10^4`$ K (3 $`\sigma `$ lower limit), by $`z=3.9`$? It is straightforward to show that this requires $`\alpha <0.12`$, from eq. (6).

One can also derive a similar limit if it is assumed helium is only singly ionized by $`z=3.9`$ (i.e. $`X_{\mathrm{HeI}},X_{\mathrm{HeIII}}1`$, $`X_{\mathrm{HeII}}1`$). Assuming again a power-law $`J_\nu \nu ^\alpha `$, but with a cut-off for $`\nu >\nu _{\mathrm{HeII}}`$, the requirement is $`\alpha <2.2`$. In other words, the spectrum needs to be even harder with no HeII photo-heating, and a spectrum as hard as $`\alpha <2.2`$ is unrealistic in comparison with stellar and quasar spectra.

To summarize: if reionization takes place at $`z>10`$, and the fractions of HI/HeI/HeII experience no significant (order unity) change after $`z=10`$, a sufficiently hard ionizing spectrum is necessary to keep the IGM temperature high enough to match observations at $`z=3.9`$. Parameterizing the ionizing background by a power-law $`J_\nu \nu ^\alpha `$, this requires $`\alpha <0.12`$ if HeII is reionized by $`z=3.9`$, or $`\alpha <2.2`$ if HeII is not reionized by then (the latter assumes $`J_\nu `$ is cut off for $`\nu >\nu _{\mathrm{HeII}}`$). One should keep in mind that the only relevant slope of the spectrum is the slope just above each ionization threshold (unless the spectrum is very hard), because the ionization cross section $`\sigma _i\nu ^3`$. In other words, $`\alpha `$ can deviate greatly from the values given above as long as the deviation takes place at frequencies far away from the ionization thresholds. Note also that the spectrum described here refers to the asymptotic spectrum at $`z=3.9`$. If the spectrum changes significantly after $`z10`$, one can view the above limits as applicable to the hardest spectrum between $`z=3.910`$.

The power index we find, $`\alpha <0.12`$, or $`\alpha <2.2`$, represents a very hard spectrum. We next turn to the question: how hard can a realistic ionizing spectrum be?

## 3. The Ionizing Spectrum

Two kinds of ionizing spectrum are generally discussed in the literature. One is quasar-like and the other is star-like.

Zheng et al. (1998) finds a quasar spectral shape of $`\nu ^{1.8}`$ for the relevant ionizing frequencies in high resolution HST spectra. In this paper, we will follow Haardt & Madau (1996) and consider, conservatively, a quasar spectrum that goes as $`\nu ^{1.5}`$. One should keep in mind that, at the relevant redshifts $`\text{ }>4`$, the known populations of quasars probably cannot contribute significantly to the (hydrogen) ionizing background. Here, we take the conservative view that there might be a population of dim quasars that still contribute significantly to a hard spectrum (Haiman & Loeb 1997).

Stellar spectra are generally softer than a typical quasar spectrum. An exception is the spectrum of metal-free stars (Tumlinson & Shull 2000, Bromm, Kudritzki & Loeb 2001, Schaerer 2002, Venkatesan, Tumlinson & Shull 2003). In their theoretical models, Bromm et al. (2001) find a spectrum that has the following form: $`J_\nu \nu `$ for $`\nu `$ just above $`\nu _{\mathrm{HI}}`$, $`J_\nu \nu ^0`$ at $`\nu \text{ }>\nu _{\mathrm{HeI}}`$, and $`J_\nu \nu ^{4.5}`$ for $`\nu \text{ }>\nu _{\mathrm{HeII}}`$. Such a spectrum is harder than the quasar spectrum at frequencies below the HeII threshold. Note that metal free stars probably cause reionization early on, but it is unlikely they contribute significantly to the asymptotic ionizing spectrum at low redshifts (Haiman & Holder 2003). We consider a metal free stellar spectrum here for the sake of being conservative i.e. assume a spectrum that is as hard as possible.

The actual ionizing spectrum seen by a fluid element is different from the above, due to processing by the IGM. Haardt & Madau (1996) have done a careful calculation of such effects. Absorption generally hardens the spectrum333Note the apparent paradox: absorption, by taking away ionizing photons, increases the photo-heating rate. This is a result of photo-ionization equilibrium, which makes the amplitude of the ionizing background irrelevant (eq. ). Rather it is the spectral shape that is important. However, diffuse recombination radiation from the absorbing medium tends to compensate for this hardening. Haardt & Madau (1996) found that the spectrum above $`\nu _{\mathrm{HeII}}`$ hardens by $`1.5`$ (to be precise, $`\alpha \alpha 1.5`$) (see also Zuo & Phinney 1993). The spectrum just above $`\nu _{\mathrm{HI}}`$ and $`\nu _{\mathrm{HeI}}`$ maintains roughly the same slope as the source.

We are therefore led to consider the following IGM modified spectra. Let us use the symbols $`\alpha _{\mathrm{HI}}`$, $`\alpha _{\mathrm{HeI}}`$ and $`\alpha _{\mathrm{HeII}}`$ to denote the spectral slopes (more precisely, its negative), above the three relevant ionizing thresholds. A quasar-like processed spectrum has $`\alpha _{\mathrm{HI}}=\alpha _{\mathrm{HeI}}=1.5`$, and $`\alpha _{\mathrm{HeII}}=0.0`$. A metal-free-stellar spectrum gives rise to $`\alpha _{\mathrm{HI}}=1`$, $`\alpha _{\mathrm{HeI}}=0`$, and $`\alpha _{\mathrm{HeII}}=3`$. At the moment(s) of reionization, the relevant spectrum can, however, be even harder (Abel & Haehnelt 1999). In principle, the spectrum can be hardened relative to the source by as much as a power-law index of $`3`$, the $`3`$ coming from the scaling $`\sigma _i\nu ^3`$ (see arguments in Abel & Haehnelt (1999), and also Zuo & Phinney 1993). We therefore allow the quasar-spectrum to have $`\alpha _{\mathrm{HI}}=\alpha _{\mathrm{HeI}}=\alpha _{\mathrm{HeII}}=1.5`$, and the stellar-spectrum to have $`\alpha _{\mathrm{HI}}=4,\alpha _{\mathrm{HeI}}=3,\alpha _{\mathrm{HeII}}=1.5`$, at the initial moment(s) of reionization. This is not relevant for the thermal asymptote, but is relevant for the magnitude of temperature boosts during reionization events.

Comparing the above spectral slopes against the limits obtained in the last section suggests neither quasars, nor metal–free stars can match the observed temperature at $`z=3.9`$, if reionization occurs at $`z>10`$, and no reionization event takes place afterwards. However, those limits were based on a strict power-law spectrum. The thermal asymptotes (eq. ) for our more realistic spectra for quasars and metal–free stars are shown in Fig. 2 and Fig. 3 (black dotted lines), respectively. Indeed, in both cases, the thermal asymptotes fall short of the observed (3 $`\sigma `$) lower limit of $`T_0>1.6\times 10^4`$ K at $`z=3.9`$, regardless of whether or not helium is doubly ionized.

Hence, if the universe is reionized before $`z10`$, and remains highly ionized thereafter, neither reasonably hard spectra can reproduce the IGM temperatures inferred from the Ly$`\alpha `$ forest.

## 4. Illustrative Examples

The conclusion from the last section, together with the large electron scattering optical depth measured by WMAP, therefore implies that one or more of the fractions $`X_{\mathrm{HI}},X_{\mathrm{HeI}},X_{\mathrm{HeII}}`$ must change by order unity at $`z<10`$, to give the IGM temperature boosts above the thermal asymptote. We therefore rule out the vanilla reionization models laid out in §1. One should keep in mind, however, this is predicated upon (reasonable) assumptions about the hardness of the ionizing spectrum. In this section, we work out some illustrative examples of what it takes to match the IGM temperature constraints, using the spectra from the last section.

Fig. 2 shows the thermal evolution for a fluid element at mean density subject to a quasar-like ionizing spectrum, and experiencing a variety of ionization histories, for both the case of having helium doubly ionized (upper panel), and the case of only having singly ionized helium (lower panel; except for model E, see below). Models A and C both describe early reionization before $`z=10`$, and no significant (order unity) changes in $`X_{\mathrm{HI}}`$, $`X_{\mathrm{HeI}}`$, $`X_{\mathrm{HeII}}`$ thereafter. As explained before, both converge to their respective thermal asymptotes by $`z4`$, which fall short of the observed temperatures, especially for model C, which has no HeII reionization, and therefore lower temperatures. Model B demonstrates how early reionization can be made consistent with the forest temperature measurements. It has a late second episode of (both hydrogen and helium I and II) reionization at $`z7`$, which boosts the temperature significantly above the thermal asymptote. Model E is similar, except that in the second episode, only HeII is reionized (HI and HeI were already highly ionized before then). It shows that HeII reionization alone can provide a significant boost to the IGM temperature.

Model D is intriguing, as it shows a case where a second episode of hydrogen (and HeI, but not HeII) reionization occurs at the smallest redshift ($`z=6`$) allowed by SDSS observations (Fan et al. 2002). Even with such a late reionization epoch, the temperature can barely match the observations (the temperature at $`z=3.9`$ is $`1.5\times 10^4`$, which is just below the 3 $`\sigma `$ limit). This is due to the lack of HeII reionization in this model. Therefore, with a quasar spectrum like the one assumed here, some amount of HeII reionization is necessary prior to $`z4`$, to heat the IGM sufficiently. If the evidence that suggests that HeII reionization is occurring at $`z3`$ holds up (see §2), then this implies either that a spectrum harder than assumed here for quasars exists prior to $`z4`$, or else HeII reionization must last for a an extended period, from $`z3`$ back to at least $`z>4`$.

Fig. 3 shows a similar exercise for a metal free stellar spectrum (as defined in §3). We emphasize again, however, that it is unlikely that metal free stellar spectra remain a dominant contribution to the ionizing background down to low redshifts (Haiman & Holder 2003). This figure should only be viewed as an illustration of possibilities. The curves are exact analogs of those in Fig. 2. Two features are different from the previous figure. Model E, where a late period of HeII reionization occurs around $`z6`$, can no longer heat up the IGM sufficiently to be consistent with observations. This is because of the softness of a stellar spectrum above the HeII threshold. On the other hand, model D, which has a late hydrogen (and HeI, but no HeII) reionization at $`z6`$, has no problem matching the observed temperatures, unlike its quasar analog shown in Fig. 2. This is because the stellar spectrum we adopted is, in fact, harder than the quasar spectrum for frequencies just above $`\nu _{\mathrm{HI}}`$ and $`\nu _{\mathrm{HeI}}`$.

Finally, before we move on to the next section, it is interesting to explore a situation that is intermediate between the two extremes we have considered above i.e. HeII is partially ionized early on, and stays partially ionized thereafter. Because the recombination rate of HeIII to HeII is high, one might wonder if the continuous photoionization and recombination could give rise to a sufficiently high temperature.444We thank Andrea Ferrara for posing the question. We show in Fig. 4 the thermal evolution (upper panel) for a fluid element which experiences reionization at $`z>10`$, with the HeII fraction (lower panel) maintaining at the $`80\%`$ level thereafter. The thermal evolution assumes a quasar spectrum identical to that used in Fig. 2 (in particular, $`\alpha _{\mathrm{HeII}}=0`$ for the processed, post-reionization spectrum; see §3), except that a break in the ionizing flux just above the HeII threshold is tuned to give rise to $`X_{\mathrm{HeII}}0.8`$. Notice how the analytic asymptote (dotted line; from eq. , with $`X_{\mathrm{HeII}}=0.8`$, and $`X_{\mathrm{HeIII}}=0.2`$) remains an excellent approximation to the late time thermal evolution. One can see that the late time temperature is intermediate between the cases of full HeII reionization and of no HeII reionization, depicted in Fig. 2. In other words, partial HeII reionization is nicely bracketed by the cases we have considered already. A caveat one should keep in mind, however, is that when $`X_{\mathrm{HeII}}`$ is close to unity, the processed ionizing spectrum is expected to be still harder above the HeII threshold compared to what we have already assumed. Trying $`\alpha _{\mathrm{HeII}}=1.5`$ (which is the hardest possible; see §3) in place of $`\alpha _{\mathrm{HeII}}=0`$ in Fig. 4, results in a temperature at $`z4`$ of $`1.3\times 10^4`$ K, which still falls short of the observed temperature. For smaller $`X_{\mathrm{HeII}}`$, one expects the spectrum to be softer.

## 5. Stochastic Reionization History – Results from a Physically Motivated Model

In the above sections we have used toy models of reionization to illuminate the key issues that determine the temperature evolution of the IGM. It is interesting to consider this evolution in a physically motivated model of the reionization history that appears to fit all the relevant observations (including the electron scattering optical depth $`\tau _e=0.17`$ measured by WMAP). Based on assumptions about the nature and efficiency of the ionizing sources, the reionization history can be predicted from “first principles” in numerical simulations (e.g .Gnedin & Ostriker 1997; Nakamoto, Umemura, & Susa 2001; Gnedin 2001; Razoumov et al. 2002), and in semi–analytical models (e.g. Shapiro, Giroux & Babul 1994; Tegmark et al. 1994; Haiman & Loeb 1997, 1998; Valageas & Silk 1999; Wyithe & Loeb 2003; Cen 2003). Here we consider a semi–analytical model adopted from Haiman & Holder (2003), but modified to include the ionization of HeII$``$HeIII. For technical details, the reader is referred to that paper. Inclusion of HeIII is conceptually straightforward.

In this model, we follow the volume filling fractions $`x_{\mathrm{HII}}`$ and $`x_{\mathrm{HeIII}}`$ of HII and HeIII, assuming that discrete ionized Strömgren spheres are being driven into the IGM by ionizing sources located in dark matter halos. As can be seen from eq. (5), photoionization of HeI (to HeII) plays a sub-dominant role in heating the IGM. For simplicity, we assume that HeI is reionized at the same time as HI, so there is no need to separately keep track of the filling fraction of HeI. Note that before the Strömgren spheres percolate, the radiation background is extremely inhomogeneous: fluid elements inside ionized regions see the flux of a single (or a cluster of a few) sources, whereas fluid elements in the still neutral regions see zero flux.555These statements would no longer hold if the early ionizing sources had a hard spectrum extending to $`\text{ }>`$1 keV energies; an interesting possibility (e.g. Oh 2001, Venkatesan, Giroux & Shull 2001) that we do not consider in this paper. In this picture, each fluid element is engulfed by an ionization front at a different time. In effect, each fluid element therefore has a different reionization history. Rather than considering a single temperature for a fluid element at the mean density, it is more appropriate to consider an ensemble of fluid elements at the mean density, each with a different reionization history and temperature evolution. Note that this stochasticity is in addition to the IGM having a distribution of temperatures due to density variations (see §6 for a discussion of the latter).

The evolutions of $`x_{\mathrm{HII}}`$ and $`x_{\mathrm{HeIII}}`$ in our model are shown by the solid and dashed curves in the top panel of Figure 5 (ignore the QSO curve to the left for the moment). Reionization has an interesting history that reflects contributions from three distinct types of ionizing sources, and two different feedback effects (all of which have physical motivations as described in detail in Haiman & Holder 2003). In short, ionizing sources (assumed to be massive metal–free stars) first appear inside gas that cools via $`\mathrm{H}_2`$ lines, and collects in the earliest non–linear halos with virial temperature of $`100\mathrm{K}\text{ }<T\text{ }<10^4`$K. These sources ionize $`50\%`$ of the volume in hydrogen, and they have sufficiently hard spectra (see discussion above) that they reionize $`10\%`$ of the helium. However, at this stage (redshift $`z17`$) the entire population of these first generation sources effectively shuts off due to global $`\mathrm{H}_2`$–photodissociation by the cosmic soft UV background they had built up. Soon more massive halos, with virial temperatures of $`10^4\mathrm{K}\text{ }<T\text{ }<2\times 10^5`$K form, which do not rely on $`\mathrm{H}_2`$ to cool their gas (they cool via neutral H excitations), and new ionizing sources turn on in these halos. These are assumed to be “normal” stars, since the gas had already been enriched by heavy elements from the first generation. These sources continue ionizing hydrogen, but since they produce little flux above the HeII edge, helium starts recombining. This second generation population is also self–limiting: gas infall to these relative shallow potential wells is prohibited inside regions that had already been ionized and photo–heated to $`10^4`$K. As a result, hydrogen reionization starts slowing down around $`z10`$ (see the solid curve in the bottom panel of Fig. 5). However, at this stage, still larger halos with virial temperatures of $`T\text{ }>2\times 10^5`$K start forming. These relatively massive, third generation halos are impervious to photoionization feedback, and complete the reionization of hydrogen.

As far as the late time thermal state is concerned, the relevant ionizing spectrum is that of the sources that turn on after $`z17`$. These are Population II stars – a reasonable spectrum is $`\alpha _{\mathrm{HI}}=1`$, and $`\alpha _{\mathrm{HeI}}=4`$, and heavily truncated beyond $`\nu _{\mathrm{HeII}}`$ (Leitherer et al. 1999). As we have illustrated with examples in the last section, a spectrum as soft as this, even with some amount of late (HI, HeI but not HeII) reionization after $`z10`$, has difficulty heating up the IGM to high enough temperatures (even after accounting for processing of spectrum by the IGM, see §3).

We are therefore led to consider the effect of quasars. Sokasian, Abel & Hernquist (2002) showed, based on a 3D radiative transfer simulation, that HeII reionization should occur around $`z4`$, driven by the population of known quasars. 666Wyithe & Loeb (2002) considered HeII reionization by the known quasar population in a semi-analytical model and obtained a similar result: HeII reionization extends beyond $`z4`$, but is complete by around then. According to this work, HeII$``$HeIII reionization completes within a relatively short redshift range (see the left dashed curve in the top panel of Fig. 5). We approximate their probability distribution of (late) HeII reionization redshifts as a Gaussian centered at $`4.1`$ with a full-width-at-half-maximum of $`0.5`$ (dashed curve in the bottom panel of Fig. 5). The quasar spectrum is as described in §3.

We generate reionization histories for an ensemble of fluid elements (at mean density) in a stochastic way: determining the redshifts of HI (and HeI) reionization, and HeII reionization by drawing from the probability distributions described in the bottom panel of Fig. 5. The resulting temperatures at $`z=2.4`$, $`3.0`$, and $`3.9`$ have a scatter because of the stochastic history, and their probability distributions are shown in Fig. 6.

The results shown in this figure are intriguing. At $`z=3.9`$, the distribution (red solid histogram) peaks around $`T_0=2.2\times 10^4`$ K, but it has a significant bump below $`10^4`$ K as well. This bimodal distribution is a result of the stochastic reionization history: at $`z=3.9`$, there is a minority of fluid elements that have not undergone HeII reionization, and so they are significantly colder, by more than a factor of $`2`$. There have been attempts in the past to look for temperature fluctuations (beyond what one expects from variations with density) in the Ly$`\alpha `$ forest (Zaldarriaga 2002, Theuns et al. 2002b). So far there has been no detection. Our results indicate that the stochasticity of reionization can detectably increase the width of the temperature distribution. It will be very interesting to confirm our results by a more detailed analysis that folds in a realistic density distribution, and to apply the observational search techniques to larger datasets and to higher redshifts ($`z3`$ so far).

The formal mean and rms scatter of $`T_0`$ at $`z=3.9`$ is $`1.8\times 10^4\pm 7\times 10^3`$ K. It can be seen from Fig. 6 that the scatter gets progressively smaller as one goes to lower redshifts: at $`z=3.0`$, $`T_0=1.7\times 10^4\pm 8\times 10^2`$ K; at $`z=2.4`$, $`T_0=1.5\times 10^4\pm 5\times 10^2`$ K. The scatter at these lower redshifts are no larger than the expected scatter from shock-heating as well as dynamics (e.g. Croft et al. 1997, Hui & Gnedin 1997, Dave & Tripp 2001).

## 6. Discussion

We find that the temperature of the IGM, as inferred from Ly$`\alpha `$ forest spectra at redshifts $`z24`$, leads to several interesting conclusions. Especially interesting are the conclusions we obtain when the Ly$`\alpha `$ forest temperature is considered together with the WMAP results. Our main conclusions are summarized as follows:

* A vanilla reionization model, where $`X_{\mathrm{HI}},X_{\mathrm{HeI}}`$, and/or $`X_{\mathrm{HeII}}`$ undergo order unity changes at $`z>10`$, but suffer no such changes at $`z<10`$, is ruled out by temperature measurements of ZHT01, especially at $`z=3.9`$. This relies on assumptions about the ionizing spectrum, which are discussed in §3. For a power-law spectrum, rigorous requirements can be put on the hardness of the ionizing background to evade the the above argument: $`\alpha <0.12`$ if helium is doubly ionized, or $`\alpha <2.2`$ if helium is singly ionized, for $`J\nu ^\alpha `$ (with a cutoff at $`\nu _{\mathrm{HeII}}`$ if helium is only singly ionized). The idea of an asymptotic evolution for the thermal state of the IGM is quite useful in formalizing the above argument. The asymptote (at $`2z4`$) is described quite accurately ($`5\%`$) by eq. (4). It is reached as long as any order unity changes in $`X_{\mathrm{HI}},X_{\mathrm{HeI}}`$, and/or $`X_{\mathrm{HeII}}`$ occur prior to $`z=10`$.
* Conversely, the requirement by WMAP that the universe reionizes early, at $`z>10`$, implies the reionization history is complex. To fulfill the temperature constraints, there must be additional periods of order unity changes in one or more of the fractions $`X_{\mathrm{HI}},X_{\mathrm{HeI}}`$, $`X_{\mathrm{HeII}}`$ at redshift below 10. This is an argument separate from the argument based on the evolution of mean transmission (Haiman & Holder 2003; based on the Gunn-Peterson troughs observed by Becker et al. 2001 and Fan et al. 2003 at $`z6`$, and comparison against an extrapolation from lower redshifts).
* Exactly what kind of changes at $`z<10`$ are necessary to match the temperature constraint at $`z=3.9`$ depends critically on the ionizing spectrum. If the spectrum is harder than $`\nu ^{1.5}`$ above the HI and HeI ionization thresholds (i.e. harder than the ’quasar’ spectrum discussed in §3), then order unity changes in $`X_{\mathrm{HI}}`$ and $`X_{\mathrm{HeI}}`$ in the period $`6\text{ }<z\text{ }<10`$ (the lower limit of $`6`$ being set by the SDSS mean transmission measurements; Fan et al. 2002), even without accompanying changes in $`X_{\mathrm{HeII}}`$ (i.e. no HeII reionization), are sufficient to heat up the IGM to within the constraint at $`z=3.9`$. Conversely, if the spectrum is not hard enough (for instance, a population II stellar spectrum described in §5), some amount of HeII reionization-heating is necessary prior to $`z=3.9`$.
* The IGM temperature can have a broad distribution (beyond what one expects based on variation with density) close to reionization events, but the scatter diminishes with time. The reionization history of the universe is almost certainly stochastic, in the sense that different fluid elements reionize at different times, depending on when a given element becomes engulfed in expanding HII (or HeIII) regions around ionizing sources. Based on simple semi-analytic models, we predict the probability distribution of ionization redshifts. The resulting temperature scatter is particularly large close to reionization events, but diminishes quickly with time. This is illustrated in Fig. 6. It would be very interesting to search for the broad initial temperature scatter that is predicted by physically motivated reionization models.
* The known quasar population at $`z<5`$ may heat the IGM to sufficiently high temperatures by reionizing HeII. In §5, we describe a physically motivated reionization model with a stochastic history. Hydrogen (and HeI) reionization is accomplished early on by metal-free stars ($`17\text{ }<z\text{ }<30`$), and later on by population II stars ($`z\text{ }<17`$), and completes by $`z6.5`$. This fulfills the dual requirements of a large CMB optical depth, and the Gunn-Peterson trough seen by SDSS. Because the population II stars have a rather soft spectrum, photo-heating of HI and HeI alone is not enough to bring the temperature up to the observed level, and this model fails to predict sufficiently high IGM temperatures – despite the percolation occurring at the relatively low redshift of $`z7`$. We therefore make use of the predictions of Sokasian et al. (2002) for a brief period of HeII reionization around $`z4`$, based on the known population of quasars. We find that such a model can satisfy the temperature constraints. It is important to note, however, that this is not the only way to achieve the observed temperatures. For instance, a population of dim quasars can turn on at $`z>6`$ (Haiman & Loeb 1997), which has a sufficiently hard spectrum to either heat the IGM via photo-heating of HI and HeI alone ( but their spectrum has to be harder than $`\nu ^{1.5}`$ if HeII is not reionized), or boost the temperature via HeII reionization as well. If HeII reionization takes place via these mini-quasars, the turn-on of the known population of quasars at $`z\text{ }<4`$ then has a much weaker influence on the thermal state at low redshifts.
* HeII reionization makes little difference to the electron scattering optical depth, particularly if it occurs late. For instance, in the model shown in Fig. 5, extra electrons from HeII reionization changes $`\tau _e`$ at a $`2\%`$ level fractionally.

An important caveat in our discussion so far is the possible existence of additional heating mechanisms. Two such possibilities are galactic outflows and Compton heating by a hard X-ray background. Adelberger et al. (2002) recently observed signatures of galactic winds into the IGM. Such outflows can in principle heat up the IGM. However, observations by Rauch et al. (2001) of close pairs of lines of sight in lens systems suggest that the IGM is not turbulent on small-scales, arguing against significant stir-up of the IGM by winds. Moreover, galactic outflows can heat up the IGM to a variety of temperatures – the fact that the observed temperatures are in the range of expectations for a photo-ionized, and photo-heated, gas suggests photo-heating is the simplest explanation.

Another important question is whether Compton heating by a hard X-ray background (XRB) could be more important than photoelectric heating. This question was considered by Madau & Efstathiou (1999), who assumed that the redshift evolution of the sources of the hard XRB parallels the flat distribution that had been determined for the soft X–ray AGN luminosity function beyond $`z2`$. Under this assumption, they found the hard XRB to be an important source of heating relative to the UV background at redshifts $`z>2`$, raising the IGM temperature to $`1.5\times 10^4`$K at $`z4`$. Two new developments make it unlikely for the hard XRB to be an important source of heating. First, as pointed out by Abel & Haehnelt (1999), the photoelectric heating rate can be significantly increased in optically thick gas (when hydrogen and helium first gets ionized); we here adopt these increased rates. Second, the sources of the hard XRB have been resolved by the Chandra satellite, and, rather than paralleling the soft X–ray luminosity function, the hard X–ray sources exhibit a steep decline towards high redshift beyond $`z2`$ (Cowie et al. 2003).

Our investigation raises a number of interesting issues. Is HeII reionization prior to redshift $`4`$ necessary to heat up the IGM to the right level? This depends critically on what kind of sources and spectra are available. The prediction for a large scatter in temperature (factor of about $`2`$) close to the epoch of HeII reionization is something one could look for, if indeed HeII reionization occurs late (Zaldarriaga 2002, Theuns et al. 2002b). In this paper, we have focused entirely on the thermal state of fluid elements at the mean density. The formalism of Hui & Gnedin (1997) can be used to compute the same for elements at $`\delta \rho /\rho \text{ }<10`$. It is interesting to explore how the variation of temperature with density ($`T\rho ^{\gamma 1}`$, an effective equation of state) might place additional constraints on the reionization history. At present, measurements of $`\gamma `$ are quite noisy (e.g. ZHT01). New approaches to constrain it better will therefore be very useful (Dijkstra, Lidz & Hui 2003). Finally, there is the issue of exotic ionization mechanisms (e.g. Berezhiani & Khlopov 1990) which is particularly interesting in the face of early reionization.

We thank Renyue Cen, Andrea Ferrara, Wayne Hu, Avi Loeb, Jerry Ostriker, Joop Schaye, and especially Adam Lidz, for useful discussions. LH is supported in part by an Outstanding Junior Investigator Award from the DOE, an AST-0098437 grant from the NSF, and by the DOE at Fermilab, and NASA grant NAG5-10842.