# The thermal history of the intergalactic medium The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

## 1 INTRODUCTION

According to the standard big bang model, the primordial hydrogen and helium comprising the intergalactic medium (IGM) was hot and highly ionized at early times. As the universe expanded, the hot plasma cooled adiabatically, becoming almost completely neutral at a redshift of $`z10^3`$. The IGM remained neutral until the first stars and quasars began to produce ionizing photons. Eventually, the ionizing radiation became intense enough to reionize hydrogen and later, because of its higher ionization potential, to fully reionize helium. Since the thermal evolution of the IGM depends strongly on its reionization history, it can be used as a probe of the end of the ‘dark ages’ of cosmic history, when the first stars and quasars were formed \[Miralda-Escudé & Rees 1994, Hui & Gnedin 1997, Haehnelt & Steinmetz 1998\].

The absence of Gunn-Peterson absorption \[Gunn & Peterson 1965\] in quasar spectra, i.e. the complete absorption of quasar light blueward of the H i and He ii Ly$`\alpha `$ wavelengths requires that hydrogen must have been highly ionized by $`z5`$ \[Schneider, Schmidt & Gunn 1991, Songaila et al. 1999\] and helium by $`z2.5`$ (Davidsen, Kriss & Zheng 1996). Measurements of the He ii Ly$`\alpha `$ opacity suggest that helium may have reionized around $`z3`$ \[Heap et al. 2000, Reimers et al. 1997, Jakobsen et al. 1994, Davidsen et al. 1996, Anderson et al. 1999\]. This would fit in with evidence for a hardening of the UV background around this time, as derived from the ratio of Si iv/C iv in high redshift quasar absorption lines \[Songaila & Cowie 1996, Songaila 1998\], although both the observational result and its interpretation are still controversial \[Boksenberg, Sargent & Rauch 1998, Giroux & Shull 1997\].

The resonant Ly$`\alpha `$ absorption by residual low levels of neutral hydrogen along the line of sight to a quasar produces a forest of absorption lines. Although many of the basic observational facts about the Ly$`\alpha `$ forest at high redshift ($`z2`$–5) had been established before the 10 m telescope era, the advent of the Keck telescope has lead to much larger data samples at much higher signal-to-noise ratio than hitherto available \[e.g. Hu et al. 1995, Lu et al. 1996, Kirkman & Tytler 1997\]. The observational progress has been matched on the theoretical side by semi-analytic models \[e.g. Bi & Davidsen 1997\] and cosmological hydro-simulations \[e.g. Cen et al. 1994, Zhang, Anninos & Norman 1995, Petitjean, Mücket & Kates 1995, Hernquist et al. 1996, Miralda-Escudé et al. 1996\] which together with the new data are now beginning to yield significant quantitative cosmological constraints (see Rauch 1998 for a recent review).

These simulations show that the low column density ($`N10^{14.5}\mathrm{cm}^2`$) absorption lines arise in a smoothly varying IGM of low density contrast ($`\delta 10`$), which contains most of the baryons in the universe. Since the overdensity is only mildly non-linear, the physical processes governing this medium are well understood and relatively easy to model. On large scales the dynamics are determined by gravity, while on small scales gas pressure is important. Since shock heating is unimportant for the low-density gas, the interplay between photoionization heating and adiabatic cooling due to the expansion of the universe results in a tight temperature-density relation, which is well described by a power-law for densities around the cosmic mean, $`T=T_0(\rho /\overline{\rho })^{\gamma 1}`$ \[Hui & Gnedin 1997\]. This relation is generally referred to as the ‘equation of state’ (even though the true equation of state is that of an ideal gas).

For models with abrupt reionization, the IGM becomes nearly isothermal ($`\gamma 1`$) at the redshift of reionization. After reionization, the temperature at the mean density ($`T_0`$) decreases while the slope ($`\gamma 1`$) increases because higher density regions undergo increased photoheating and expand less rapidly. Eventually, the imprints of the reionization history are washed out and the equation of state approaches an asymptotic state, $`\gamma =1.62`$, $`T_0\left[\mathrm{\Omega }_bh^2/\sqrt{\mathrm{\Omega }_mh^2}\right]^{1/1.7}`$ \[Miralda-Escudé & Rees 1994, Hui & Gnedin 1997, Theuns et al. 1998\]. However, the timescale for recombination cooling in the low density IGM is never small compared to the age of the universe for $`z20`$ and inverse Compton cooling of free electrons off the cosmic microwave background is only efficient for $`z5`$. Consequently, unless both hydrogen and helium were fully reionized at redshifts considerably higher than this, the gas will have retained some memory of when and how it was reionized.

A standard way of analyzing Ly$`\alpha `$ forest spectra is to decompose them into a set of distinct absorption lines, assumed to have Voigt profiles (e.g. Carswell et al. 1987). Various broadening mechanisms, such as Hubble broadening (the differential Hubble flow across the absorber), peculiar and thermal velocities contribute to the line widths (Meiksin 1994; Hui & Rutledge 1999; Theuns, Schaye & Haehnelt 2000). However, there exists a lower limit to the line-width, set by the temperature of the gas. Because the physical density of the IGM correlates strongly with the column density of the absorption lines, this results in a cut-off in the distribution of line widths ($`b`$-parameters) as a function of column density, which traces the equation of state of the gas (Schaye et al. 1999, hereafter STLE; Ricotti, Gnedin & Shull 2000; Bryan & Machacek 2000). Hence we can infer the equation of state of the IGM by measuring the minimum Ly$`\alpha `$ line width as a function of column density.

Here, we measure the $`b(N)`$ cut-off in nine high resolution, high S/N quasar spectra, spanning the redshift range 2.0–4.5. We use hydrodynamic simulations to calibrate the relations between the parameters of the $`b(N)`$ cut-off and the equation of state. By applying these relations to the observations, we are able to measure the evolution of the equation of state over the observed redshift range. We find that the thermal evolution of the IGM is drastically different from that predicted by current models. The temperature peaks at $`z3`$, which, together with supporting evidence from measurements of the He ii opacity and the Si iv/C iv ratios, we interpret as evidence for the second reionization of helium (He ii $``$ He iii). Ricotti et al. \[Ricotti et al. 2000\] recently applied a similar technique to published lists of Voigt profile fits. A comparison with the method and results of Ricotti et al. is given in section 8.

This paper is organized as follows. In sections 2 and 3 we describe the observations and the simulations respectively. We discuss the difference between evolution of the $`b`$-distribution and evolution of the temperature in section 4. In section 5 we briefly describe our method for measuring the equation of state, before we present our results in section 6. Systematic errors are discussed in section 7. Finally, we discuss and summarize the main results in section 8.

## 2 OBSERVATIONS

We analyzed a sample of nine quasar spectra, spanning the redshift range $`z_{\mathrm{em}}=2.14`$–4.55 (Table 1). The spectra of Q1100$``$264 and APM 08279+5255 were kindly provided by R. Carswell and S. Ellison respectively. All spectra were taken with the high-resolution spectrograph (HIRES, Vogt et al. 1994) on the Keck telescope, except the spectrum of Q1100$``$264, which was taken with the UCL echelle spectrograph of the Anglo Australian Telescope. Details on the data and reduction procedures, as well as the continuum fitting, can be found in Carswell et al. \[Carswell et al. 1991\] for Q1100$``$264, Ellison et al. \[Ellison et al. 1999\] for APM 08279+5255 and Barlow & Sargent \[Barlow & Sargent 1997\] and Rauch et al. \[Rauch et al. 1997\] for the others. The nominal velocity resolution (FWHM) was $`8\mathrm{km}\mathrm{s}^1`$ for Q1100$``$264 and $`6.6\mathrm{km}\mathrm{s}^1`$ for the others and the data were rebinned onto $`0.04\mathrm{\AA }`$ pixels on a linear wavelength scale. The signal to noise ratio per pixel is typically about 50, except for Q1100$``$264 for which it is about 20.

In order to avoid confusion with the Ly$`\beta `$ forest, only the region of a spectrum between the quasar’s Ly$`\beta `$ and Ly$`\alpha `$ emission lines was considered. In addition, spectral regions close to the quasar (typically 8–10$`h^1\mathrm{Mpc}`$, but $`21h^1\mathrm{Mpc}`$ for APM 08279+5255 and $`32h^1\mathrm{Mpc}`$ for Q1100$``$264) were omitted to avoid proximity effects. Regions thought to be contaminated by metals and damped Ly$`\alpha `$ lines were removed (metal line regions were identified by correlating with metal lines redwards of the quasar’s Ly-alpha emission line and with strong H i lines). The absorption features in the remaining spectral regions were fitted with Voigt profiles using the same automated version of VPFIT \[Webb 1987, Carswell et al. 1987\] as was used for the simulated spectra. Using a fully automatic fitting program invariably results in a few ‘bad fits’. However, given that there is no unique way of decomposing intrinsically non-Voigt absorption lines into a set of discrete Voigt profiles, it is essential to apply the same algorithm to simulated and observed spectra. Since ‘bad fits’ will also occur in the synthetic spectra, we have made no attempt to correct them.

The Ly$`\alpha `$ forest of a single quasar spans a considerable redshift range ($`\mathrm{\Delta }z0.5`$). In order to minimize the effects of redshift evolution and S/N variation across a single spectrum, we divided each Ly$`\alpha `$ forest spectrum into two parts of equal length. STLE showed that their algorithm for measuring the cut-off of the $`b(N)`$ distribution is relatively insensitive to the number of absorption lines, the statistical variance is almost the same for e.g. 150 and 300 lines. Hence little information is lost by analyzing narrow redshift bins if the absorption line density is high. The two halves of the spectra were analyzed separately and each was compared with its own set of simulated spectra (see section 3). For the two lowest redshift quasars the number of absorption lines is too small to split the data in half. Hence each quasar, except for Q1100$``$264 and Q2343+123, provides two nearly independent data sets.

The complete set of absorption line samples is listed in Table 2. The median, minimum and maximum redshifts of the absorption lines used to determine the $`b(N)`$ cut-off are listed in columns 2–4. The minimum column density considered was set to $`10^{12.5}\mathrm{cm}^2`$ for all samples, since blends dominate at lower column densities. The maximum column densities are listed in column 5, they were determined by the following considerations: (a) the cut-off can be measured more accurately if the column density interval is larger; (b) we only want to measure the cut-off for column densities that correspond to the density range for which the gas follows a power-law temperature-density relation. The total number of absorption lines in each sample is listed in column 6 (only lines for which VPFIT gives relative errors in both the $`b`$-parameter and column density less than 0.25 are considered). Finally, column 6 lists the pivot column density of the power-law cut-off, $`b=b_{N_0}(N/N_0)^{\mathrm{\Gamma }1}`$, that was fit to the data (c.f. section 5).

Scatter plots of the $`b(N)`$-distribution for all observed samples (Table 2) are shown in Fig. 1. The solid lines are the measured cut-offs, the vertical dashed lines indicate the maximum column density used for fitting the cut-off. There are clear differences between the samples. We will show in section 4 that even a non-evolving $`b`$-distribution would imply a strong thermal evolution. Several samples contain a few lines that fall far below the cut-off. These lines, which have no significant effect on the measured cut-off, are most likely blends or unidentified metal lines.

Fig. 2 shows the effective optical depth, $`\tau _{\mathrm{eff}}\mathrm{ln}(F)`$, of the observed samples as a function of (decreasing) redshift. The ionizing background in the simulations was rescaled to match these effective optical depths. The scatter is small, considering that most data points represent just half of a Ly$`\alpha `$ forest spectrum. Note that Rauch et al. \[Rauch et al. 1997\] studied the opacity of the forest using seven out of nine of the quasars from this sample. They found a slightly less rapid increase with redshift, because they rebinned the data into three redshift bins, centered on $`z=2`$, 3 and 4.

## 3 SIMULATIONS

In order to calibrate the relation between the parameters of the $`b(N)`$ cut-off and the effective equation of state, we have simulated eight variants of the currently favoured flat, scale invariant, cosmological constant dominated cold dark matter model, which vary only in their heating rates (Table 3). The calibration was repeated for each of the observed samples of absorption lines listed in Table 2. Synthetic spectra were computed along 1200 random lines of sight through the simulation box at the nearest redshift output ($`\mathrm{\Delta }z=0.25`$). The background flux was rescaled such that the mean effective optical depth in the simulated spectra matches that of the observed sample. Each spectrum was convolved with a Gaussian with full width at half maximum (FWHM) identical to that of the observations and resampled onto pixels of the same size. The noise properties of the observed spectrum were computed as a function of flux and imposed on the simulated spectra. The resulting spectra were continuum fitted as described in Theuns et al. \[Theuns et al. 1998\]. Finally, Voigt profiles were fitted using the same automated version of VPFIT as was used for the observed spectra. We will refer to the sample of lines drawn from the synthetic spectra of simulation X, designed to mimic the observed spectrum Y as as model X-Y, e.g. model L1-1442a.

All models have a total matter density $`\mathrm{\Omega }_m=0.3`$, vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, baryon density $`\mathrm{\Omega }_bh^2=0.019`$, present day Hubble constant $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and the amplitude of the initial power spectrum is normalized to $`\sigma _8=0.9`$. The IGM is assumed to be of primordial composition with a helium abundance of 0.24 by mass and is photoionized and photoheated by the UV-background from quasars. The exact cosmology and UV-background is unimportant for this analysis, as is the normalization (see section 7). Although these parameters may affect the equation of state of the IGM in the simulations, they do not change its relation to the $`b(N)`$ cut-off.

The numerical simulations used in this paper follow the evolution of a periodic, cubic region of the universe and are performed with a modified version of HYDRA \[Couchman, Thomas & Pearce 1995\], which uses smooth particle hydrodynamics \[Lucy 1970, Gingold & Monaghan 1977\]. The simulations employ $`64^3`$ gas particles and $`64^3`$ cold dark matter particles in a box of of comoving size 3.85 Mpc, so the particle mass is $`1.14\times 10^6\mathrm{M}_{}`$ for the gas and $`6.51\times 10^6\mathrm{M}_{}`$ for the dark matter.

Our reference model, L1, is photoionized and photoheated by the UV-background from quasars as computed by Haardt & Madau (1996, hereafter HM), using the optically thin limit. Models L0.3, L2 and L3 are identical, except that we have multiplied the helium photoheating rates (column $`ϵ_{\mathrm{He}}`$ in Table 3) by factors of 1/3, 2 and 3 respectively (keeping the ionization rates constant). The effective helium photoheating rate may be higher than computed in the optically thin limit because of radiative transfer effects (e.g. Abel & Haehnelt 1999). Model Lx is identical to model L1, except that we have included Compton heating by the hard X-ray background as computed by Madau & Efstathiou \[Madau & Efstathiou 1999\] (column $`_\mathrm{X}`$ in Table 3). For a highly ionized plasma, the energy input per particle from Compton scattering of free electrons is independent of the density. Hence, Compton heating tends to flatten the effective equation of state. We have used this fact to artificially construct models with low values of $`\gamma `$, by multiplying the X-ray heating rates by (unrealistic) factors of 2.5 and 5 for models Lx2.5 and Lx5 respectively. Finally, model L1e is identical to model L1, except that we have set the ionization and heating rates for H i and He i for redshifts between 6 and 10 equal to those at $`z=6`$. In this model H i and He i ionize early (at $`z=10`$), which drives $`\gamma `$ to larger values.

In addition to the models listed in Table 3, we have performed some simulations to investigate possible systematic effects. We simulated model L1 twice with lower resolutions ($`2\times 54^3`$ and $`2\times 44^3`$ particles) and model L3 in a larger box, but with the same resolution ($`5h^1\mathrm{Mpc}`$ and $`2\times 128^3`$ particles instead of $`2.5h^1\mathrm{Mpc}`$ and $`2\times 64^3`$ particles). Model L1 was simulated twice more with lower normalizations of the initial power spectrum ($`\sigma _8=0.65`$ and 0.4 instead of 0.9).

## 4 Evolution of the $`b`$-distribution vs thermal evolution

Williger et al. \[Williger et al. 1994\] and Lu et al. \[Lu et al. 1996\] found that the $`b`$-parameters at $`z4`$ are smaller than at $`z=2`$–3. Kim et al. \[Kim et al. 1997\] showed that the increase in the line widths with decreasing redshift continues over the range $`z=3.5`$ to 2.1. It is tempting to interpret these results as evidence for an increase in the temperature $`T_0`$ with decreasing redshift. However, we will show in this section that the $`b`$-values are smaller at higher redshift even for models in which $`T_0`$ is higher, as is the case for models in which the universe is fully reionized by $`z=4`$.

As pointed out by STLE, any statistic that is sensitive to the temperature of the absorbing gas, will in general depend on both the amplitude, $`T_0`$, and the slope, $`\gamma `$, of the equation of state. This is because temperature is a function of density and the absorbing gas is, in general, not all at the mean density of the universe. After reionization, $`T_0`$ will decrease and $`\gamma `$ will increase with time. Consequently, the evolution of the temperature at a given overdensity can be very different from the evolution of $`T_0`$. This is illustrated in Fig. 3, where the temperature $`T_\delta `$ at a density contrast $`\delta \rho /\overline{\rho }1`$ is plotted as a function of redshift for model L1. Even though the temperature at the mean density decreases with time (solid line), the temperature at a density contrast as little as 2 remains almost constant.

The general expansion of the universe ensures that the column density corresponding to a fixed overdensity is a strongly increasing function of redshift. In fact, most of the evolution of the Ly$`\alpha `$ forest can be understood in terms of the resulting scaling of the optical depth (e.g. Hernquist et al. 1996; Davé et al. 1999; Machacek et al. 1999). When interpreting the evolution of the $`b`$-distribution, one therefore has to keep in mind that: (a) at fixed column density, absorption lines at higher redshift will correspond to absorbers of smaller overdensities; (b) the evolution of the temperature at a fixed overdensity depends on the evolution of both $`T_0`$ and $`\gamma `$. Together these effects can conspire to make the $`b`$-parameters smaller at higher redshift, even when the temperature $`T_0`$ is higher. Fig. 4 shows that this will happen for models in which the IGM is fully reionized at the observed redshifts ($`z4`$).

The discussion in this section shows that to derive the evolution of $`T_0`$ using a statistic that is sensitive to the temperature of the absorbing gas, one needs to determine the evolution of: (1) the temperature of the gas; (2) the overdensity of the gas and (3) the slope of the equation of state. We will see that the uncertainty in $`\gamma `$ is the limiting factor.

## 5 Method

STLE demonstrated that the observed cut-off in the distribution of $`b`$-parameters as a function of column density can be used to measure the equation of state of the IGM. In particular, they showed that the $`b(N)`$ cut-off can be fitted by a power-law, $`b=b_{N_0}(N/N_0)^{\mathrm{\Gamma }1}`$, whose parameters $`\mathrm{log}b_{N_0}`$ and $`\mathrm{\Gamma }1`$ are proportional to $`\mathrm{log}T_{\delta (N_0)}`$ and $`\gamma 1`$ respectively. For each observed sample, we first calibrate these relations using the simulations and then use them to convert the observed cut-offs into measurements of the equation of state.

For each observed sample of absorption lines (Table 2) we go through the following procedure. First mock spectra are generated from the 8 simulations listed in Table 3. The synthetic spectra are processed to give them the same characteristics (mean absorption, resolution, pixel size and noise properties) as the corresponding observed spectra. These are then fitted with Voigt profiles using the same automated fitting package that was used for the observations. For each of the eight simulated sets of absorption lines, the $`b(N)`$ cut-off is fitted using the iterative procedure developed by STLE. We then use these 8 simulations to calibrate the relations between the parameters of the $`b(N)`$ cut-off, $`(b_{N_0},\mathrm{\Gamma })`$ and the parameters of the equation of state, $`(T_{\delta (N_0)},\gamma )`$. We measure the density contrast corresponding to the pivot column density, $`N_0`$, by using the fact $`\mathrm{log}T_{\delta (N_0)}\mathrm{log}b_{N_0}`$ (STLE), essentially because both thermal broadening and Jeans smoothing scale as the square root of the temperature.

Fig. 5 illustrates our method for measuring $`T_\delta `$ and $`\delta (N_0)`$. In the left panel the intercept of the cut-off, measured at $`N_0=10^{14.0}\mathrm{cm}^2`$ in this example, is plotted as a function of $`\mathrm{log}T_0`$ for the simulations of sample 1422a (filled circles). As expected, the relation between the intercept, $`\mathrm{log}b_{14.0}`$, and $`\mathrm{log}T_0`$ is not very tight. The large scatter arises because $`\mathrm{log}b_{14.0}`$ is not proportional to $`\mathrm{log}T_0`$, but to $`\mathrm{log}T_{\delta (N_0)}`$, where $`\delta (N_0)`$ is the density contrast corresponding to the pivot column density $`N_0=10^{14.0}\mathrm{cm}^2`$. Indeed, if we plot $`\mathrm{log}b_{14.0}`$ as a function of $`\mathrm{log}T_{\delta =1.6}`$ (right panel), the scatter becomes very small, implying that $`\delta (N_0)1.6`$.

Changing the value of $`\delta `$ shifts the data points in the plot horizontally, the amount depending on the value of $`\gamma `$ in the simulation corresponding to the data point. Because we do not know a priori what the density contrast corresponding to $`N_0`$ is, we vary $`\delta `$ and see for which value the scatter is minimal. The inset in the right panel shows the total $`\chi ^2`$ of the linear least squares fit to the data points as a function of the density contrast. For this example the scatter is minimal for $`\delta =1.6`$, which is the density contrast used in the right panel. Using the optically depth weighted density - column density relation, introduced by STLE, we find that the column density $`N=10^{14.0}\mathrm{cm}^2`$ does indeed correspond to a density contrast of about 1.6.

The dot-dashed lines in Fig. 5 show the $`b`$-value expected for pure thermal broadening, $`b=\sqrt{2k_BT/m_p}`$, where $`k_B`$ is the Boltzmann constant and $`m_p`$ is the proton mass. The relation between $`\mathrm{log}b_{14.0}`$ and $`\mathrm{log}T_{\delta =1.6}`$ lies slightly above the pure thermal broadening line, but has the same slope. This implies that the widths at the cut-off are dominated by thermal broadening, with an additional component that also scales as $`T^{1/2}`$. We identify this last component as the differential Hubble flow across the absorber, whose size is set by the Jeans smoothing scale and does indeed depend on the square root of the temperature. We consistently found that minimizing the scatter in the $`\mathrm{log}b_{N_0}`$-$`\mathrm{log}T_\delta `$ relation by varying $`\delta `$ results in a relation that has a slope of about 0.5 and that the value of $`\delta `$ found agrees well with direct measurements of the density contrast corresponding to the column density $`N_0`$. We therefore use this procedure to estimate $`\delta `$ and $`\mathrm{log}T_\delta `$ and conservatively estimate the error in the density contrast to be $`\sigma (\mathrm{log}(1+\delta ))=0.15`$. Although for this work the error from the determination of the density - column density relation does not contribute significantly to the total error in $`T_0`$, it may well become the limiting factor when a larger sample of quasars is used.

Having measured $`T_\delta `$, $`\delta `$ and $`\gamma `$ we can compute $`T_0`$,

$$\mathrm{log}T_0=\mathrm{log}T_\delta (\gamma 1)\mathrm{log}(1+\delta ).$$
(1)   
Each measurement of $`\mathrm{log}b_{N_0}`$ and $`\mathrm{\Gamma }`$ comes with associated errors, which are determined from the bootstrap distribution (c.f. STLE). These errors can be converted directly into errors in $`\mathrm{log}T_\delta `$ and $`\gamma `$ respectively, using the linear relations between the cut-off and the equation of state, determined from the simulations. To these errors we add in quadrature the residual scatter of the data points around the linear fit. The error in $`\mathrm{log}T_0`$ is then given by,

$`\mathrm{\Delta }^2(\mathrm{log}T_0)`$ $`=`$ $`\mathrm{\Delta }^2(\mathrm{log}T_\delta )+\left[\mathrm{\Delta }(\gamma )\mathrm{log}(1+\delta )\right]^2+`$ (2)   
$`\left[(\gamma 1)\mathrm{\Delta }(\mathrm{log}(1+\delta ))\right]^2.`$   
After having measured the thermal evolution of the IGM, we ran a simulation designed to match the observations (dashed lines in Fig. 6). The open square in Fig. 5 corresponds to this simulation. The difference between the evolution in the calibrating simulations and the true evolution is particularly large at $`z3`$ (c.f. Fig. 6a), the redshift of model 1422a. The fact that the open square follows the same relation as the other data points, confirms that a model whose equation of state matches the one determined from the observations using the methods described in this section, does indeed have the observed $`b(N)`$ cut-off.

## 6 RESULTS

The measured evolution of the temperature at the mean density and the slope of the effective equation of state are plotted in Fig. 6. From $`z4`$ to $`z3`$, $`T_0`$ increases and the gas becomes close to isothermal ($`\gamma 1.0`$). This behavior differs drastically from that predicted by models in which helium is fully reionized at higher redshift. For example, the solid curves correspond to our reference model, L1, which uses a uniform metagalactic UV-background from quasars as computed by HM and which assumes the gas to be optically thin. In this simulation, both hydrogen and helium are fully reionized by $`z4.5`$ and the temperature of the IGM declines slowly as the universe expands. Such a model can clearly not account for the peak in the temperature at $`z3`$ (reduced $`\chi ^2`$ for the solid curves are 6.8 for $`T_0`$ and 3.6 for $`\gamma `$). Instead, we associate the peak in $`T_0`$ and the low value of $`\gamma `$ with reheating due to the second reionization of helium (He ii $``$ He iii).

If reionization of He ii happens locally on a timescale that is short compared to the recombination timescale, which for He iii is of the order of the age of the universe at $`z3`$, then the energy density injected by photoionization will be proportional to the gas density. Consequently, the temperature increase will be independent of the density and the equation of state of the IGM will become more isothermal. The change in the slope of the equation of state at $`z3`$ is thus physically consistent with our interpretation of the peak in the temperature at the same redshift.

The dashed lines in Fig. 6 are for a model that was constructed to fit the data (reduced $`\chi ^2`$ is 0.22 for $`T_0`$ and 1.38 for $`\gamma `$). This model, for which stellar sources ionize H i and He i by $`z5`$ and quasars ionize He ii at $`z3.2`$, has a much softer spectrum at high redshift. Before reionization, when the gas is optically thick to ionizing photons, the mean energy per photoionization is much higher than in the optically thin limit \[Abel & Haehnelt 1999\]. We have approximated this effect in this simulation by enhancing the photoheating rates during reionization, so raising the temperature of the IGM.

Since the simulation assumes a uniform ionizing background, the temperature has to increase abruptly (i.e. much faster than the gas can recombine) in order to make $`\gamma `$ as small as observed. In reality, the low-density gas may be reionized by harder photons, which will be the first ionizing photons to escape from the dense regions surrounding the sources. This would lead to a larger temperature increase in the more dilute, cooler regions, resulting in a decrease of $`\gamma `$ even for a more gradual reionization. Furthermore, although reionization may proceed fast locally (as in our small simulation box), it may be patchy and take some time to complete. Hence the steep temperature jump indicated by the dashed line, although compatible with the data, should be regarded as illustrative only. The globally averaged $`T_0`$ could well increase more gradually which would also be consistent with the data.

Note that if the reionization is patchy, it would give rise to large spatial fluctuations in the temperature. Because we measure the temperature-density relation from the lower cut-off of the $`b(N)`$-distribution, our results should be regarded as lower limits to the average temperature. Absorption lines arising in a local, hot ionization bubble would not necessarily raise the observed cut-off in the $`b(N)`$-distribution.

The errors in Fig. 6 can be directly traced back to the corresponding $`b(N)`$-distributions (Fig. 1). Take for example 0827a, which has an extremely low value of $`\gamma `$, with a very large error. This is clearly due to the gap in the $`b`$-distribution at $`\mathrm{log}N12.5`$–13.0. The lack of lines in that region could be a statistical fluctuation, or it could be an indication of large variations in the temperature of the IGM.

## 7 SYSTEMATICS

In this section we will investigate whether there are any systematic effects that could affect our results.

### 7.1 Numerical resolution

The $`b`$-distribution has been shown to be very sensitive to numerical resolution \[Theuns et al. 1998, Bryan et al. 1999\]. It is therefore important to check that the lower limit to the line widths in our simulations is not set by the numerical resolution. We have resimulated our reference model, L1, which is the second coldest model, twice more at a lower resolution. The resolution was decreased by decreasing the number of particles from $`2\times 64^3`$ to $`2\times 54^3`$ and $`2\times 44^3`$ respectively, while keeping the size of the simulation box constant. The resulting probability distributions for the intercept and the slope of the $`b(N)`$ cut-off are plotted in Fig. 7 for our lowest and highest redshift samples. Only for the lowest resolution simulation do the differences become noticeable. The intercept increases slightly and the slope becomes slightly shallower, indicating that the lines at the low column density end are not resolved. We conclude that the simulations used for this work have sufficient resolution to provide an accurate determination of the effective equation of state of the IGM.

### 7.2 Simulation box size

In order to investigate the effect of the simulation box size, we resimulated model L3 using a larger box, but with the same resolution ($`5h^1\mathrm{Mpc}`$ and $`2\times 128^3`$ particles instead of $`2.5h^1\mathrm{Mpc}`$ and $`2\times 64^3`$ particles). The effect of increasing the size of the simulation box is more difficult to determine than the effect of numerical resolution. When the box size is increased, the gas becomes slightly hotter ($`T_0`$ increases by a few per cent), presumably because shock heating is more effective due to the larger infall velocities. This effect is small and since it does not affect the relation between the $`b(N)`$ cut-off and the equation of state, it is unimportant for this work. We find that using the larger simulation box results in higher values of $`T_0`$. The effect is negligible at $`z4`$ and increases to $`\mathrm{\Delta }\mathrm{log}T_00.05`$ ($`\mathrm{\Delta }T_0/T_00.12`$) at $`z2`$. The change in the derived value of $`\gamma `$ is very small ($`0.05`$). Although $`T_0`$ and $`\gamma `$ may change a bit more if we increase the box size further, the small difference between the two box sizes indicates that the effect of the box size is insignificant compared to the statistical errors.

### 7.3 Cosmology

STLE showed that the relation between the cut-off in the $`b(N)`$ distribution and the equation of state is independent of the assumed cosmology. However, the initial power spectra of the models investigated by STLE were all normalized to match the observed abundance of galaxy clusters at $`z=0`$. Bryan & Machacek \[Bryan & Machacek 2000\] claimed that the $`b`$-distribution depends strongly on the amplitude of the power spectrum, as predicted by the model of Hui & Rutledge \[Hui & Rutledge 1999\]. However, Theuns et al. \[Theuns et al. 2000\] found the dependence to be very weak, provided that the line fitting is done using an algorithm, like VPFIT, that attempts to deblend absorption lines into a set of thermally broadened components. Although the absorption features do become broader for models with less small-scale power, the curvature in the line centers does not change much. Consequently, the total number of Voigt profile components used by VPFIT will generally increase, but the fits to the line centers will change very little.

To quantify the effect of decreasing the amount of small-scale power, we resimulated our reference model L1 twice using a lower normalization of the initial power spectrum. We find that normalizing to $`\sigma _8=0.65`$ instead of $`\sigma _8=0.9`$, changes the derived values of $`T_0`$ by less than 3 per cent for all redshifts. The change in $`\gamma `$ is never greater than 0.1. For the extreme case of $`\sigma _8=0.4`$, the derived values of $`T_0`$ are about 15 per cent lower and $`\gamma `$ differs by about 0.15. We conclude that the effect of the uncertainty in the normalization of the primordial power spectrum is small.

### 7.4 Continuum fitting and the mean absorption

Errors in the continuum fit of the observed spectra, will lead to errors in the effective optical depth. Underestimating the observed continuum will decrease the measured effective optical depth. Decreasing the mean absorption in the simulations will increase the density corresponding to a given column density. Hence the slope of the cut-off will remain unchanged, but the intercept will increase, although the effect is small (STLE). Increasing the intercepts of the calibrating simulations will decrease the derived temperature $`T_\delta `$. The higher derived density contrast will work in the same direction (provided that $`\gamma >1`$), resulting in a lower $`T_0`$.

The measured effective optical depth for sample 2343 seems to be relatively high compared to the other samples (Fig. 2). We tried scaling the synthetic spectra to the effective optical depth corresponding to the dashed line in Fig. 2 at the redshift of 2343, which is about 30 per cent lower than the measured value. This resulted in an increase of $`\mathrm{log}T_0`$ by 0.04 ($`\mathrm{\Delta }T_0/T_00.09`$), while leaving $`\gamma `$ unchanged.

In addition to errors in the continuum fit of the observed spectra, the continuum fitting of the synthetic spectra could also lead to systematic errors. We checked this by repeating the analysis of the simulations corresponding to our highest redshift sample, 2237b, but this time without continuum fitting the synthetic spectra. In this case the derived value of $`T_0`$ would be 9 per cent lower, while the value of $`\gamma `$ would be higher by 0.07. Hence systematic effects in the continuum fitting of the observed and the synthetic spectra are unlikely to be important.

## 8 SUMMARY AND DISCUSSION

We have measured the cut-off in the distribution of line widths ($`b`$) as a function of column density ($`N`$) in a set of nine high-quality Ly$`\alpha `$ forest spectra, spanning the redshift range 2.0–4.5. We emphasized that the evolution of the temperature of the intergalactic medium (IGM) cannot be derived directly from the evolution of the $`b(N)`$ distribution, the decrease of the overdensity corresponding to a fixed column density with redshift has to be taken into account. We therefore used hydrodynamic simulations to calibrate the relations between the $`b(N)`$ cut-off and the temperature-density relation of the low density gas. The calibration was done separately for each observed spectrum, using synthetic spectra that were processed to give them identical characteristics as the observed spectrum. Crucially, Voigt profiles were fitted to the real and simulated spectra using the same automated fitting package (a modified version of VPFIT). We have checked possible systematic errors arising from the finite numerical resolution, the finite size of the simulation box, the amplitude of the initial power spectrum and continuum fitting errors in both synthetic and observed spectra. In all cases, the effects were small compared to the statistical errors.

The measured thermal evolution differs drastically from the scenario predicted by current models of the ionizing background from quasars, in which helium is fully reionized by $`z4.5`$. The temperature at the mean density, $`T_0`$, increases from $`z4`$ to $`z3`$, after which it decreases again (Fig. 6). The slope of the equation of state reaches a minimum at $`z3`$, where it becomes close to isothermal. More data at $`z3`$ is needed to determine whether the rise in $`T_0`$ is sharp or gradual. These results suggest that the low density IGM was reheated from $`z4`$–3, which we interpret as reheating associated with the reionization of He ii.

These results are in qualitative agreement with those reported recently by Ricotti et al. \[Ricotti et al. 2000\]. They used a method which relies on the assumption that only thermal broadening contributes to the line widths of the absorption lines at the peak of the $`b`$-distribution and used approximate simulation techniques to determine the density - column density relation. By applying their method to published lists of Voigt profile fits they found that at $`z3`$, $`\gamma `$ is smaller than would be expected if the reionization of helium had been completed at high redshift. Ricotti et al. have only three data points in the range $`z=2`$–4 which all overlap at the 0.5 sigma level, so we cannot compare the shape of the temperature evolution. However, it is interesting that they measure a temperature at $`z=2`$ and $`z=4`$ that is about seventy per cent higher than is reported here, although their error bars are sufficiently large to agree with our results at the 1$`\sigma `$ level. Since Jeans smoothing contributes to the line widths, as Fig. 5 indicates (see also Theuns et al. 2000), their method should lead to an overestimate of the temperature (by about sixty per cent for the example in Fig. 5). The contribution of Jeans smoothing to the line width may be larger for the lines near the peak of the $`b`$-distribution than indicated in Fig. 5, which is for the narrower lines near the cut-off.

Although the reionization of helium may not be the only process which can explain the peak in $`T_0`$ at $`z3`$, it appears to be the only process that can simultaneously account for the observed decrease in $`\gamma `$. Galactic winds for example, would be less important in the low density regions, and would therefore result in an increase in the slope of the equation of state. Another possible explanation is a hardening of the ionizing background at $`z3`$, as might be expected because of the increase in the number of quasars. If helium had already been ionized, this would still raise the temperature somewhat, but the effect would be stronger in the high density regions where the gas recombines faster. Hence this would also lead to an increase in $`\gamma `$, contrary to what is observed.

Recently, it was shown \[Theuns et al. 1998, Bryan et al. 1999\] that high resolution simulations of the standard cold dark model, using the ionizing background computed by HM produce a larger fraction of narrow lines than observed at $`z3`$. Different authors have proposed different solutions to this problem. Theuns et al. \[Theuns et al. 1998\] suggested that the gas temperature in the simulations was too low, while Bryan et al. \[Bryan et al. 1999\] argued that the amplitude of primordial fluctuations was too high. For a given reionization history, the temperature in the simulations could for example be increased by increasing the baryon density and the age of the universe \[Theuns et al. 1999\], including Compton heating by the hard X-ray background \[Madau & Efstathiou 1999\] and possibly photo-electric heating by dust grains \[Nath, Sethi & Shchekinov 1999\]. A comparison of our reference model L1, which uses the HM ionizing background (solid lines in Fig. 6), with the data clearly shows that the model underestimates the temperature at $`z3`$. Since the above mentioned mechanisms for increasing the temperature do not change the overall shape of the thermal evolution, a change in the reionization history is required to bring the data and simulations into agreement. Furthermore, the photoheating rates around reionization need to be enhanced to account for the fact that heating by photons with energies significantly above the ionization potential is important when the gas is not optically thin \[Abel & Haehnelt 1999\], as is generally assumed in cosmological simulations.

There are two other lines of evidence for late reionization of He ii. The first are direct measurements of the optical depth from He ii Ly$`\alpha `$ absorption, which have so far been obtained for four quasars \[Jakobsen et al. 1994, Davidsen et al. 1996, Reimers et al. 1997, Anderson et al. 1999, Heap et al. 2000\]. These observations already provide strong evidence for a drop in the mean absorption from $`z3.0`$ to 2.5. The second piece of evidence concerns a change in the spectral shape of the ionizing background. As He ii is ionized, the mean free path of hard UV photons will increase and the spectrum of the UV background will become harder. Songaila & Cowie \[Songaila & Cowie 1996\] and Songaila \[Songaila 1998\] have reported a rapid increase with decreasing redshift of the Si iv/C iv ratio at $`z3`$, which they interpreted as evidence for a sudden reionization of He ii. However, Boksenberg et al. \[Boksenberg, Sargent & Rauch 1998\] found only a gradual change with redshift. The interpretation of this metal line ratio is complicated because local stellar radiation is likely to be important \[Giroux & Shull 1997\].

It should be kept in mind that these three different types of observations probe different physical structures. Our results apply to density fluctuations around the cosmic mean, the effective optical depth depends mostly on the neutral fraction in the voids and the metal line ratios probe the high density peaks. These structures will probably not be ionized simultaneously. After the ionization front breaks through the haloes surrounding the source of He ii ionizing photons, e.g. a quasar, it will propagate quickly into the voids. The filaments, where the recombination rate is much higher, will get ionized more slowly, starting from the outside \[Miralda-Escudé, Haehnelt & Rees 2000, Gnedin 2000\]. The IGM will still be thick to Ly$`\alpha `$ photons when only a small neutral fraction remains in the voids. Hence a drop in the He ii Ly$`\alpha `$ optical depth at $`z3`$ would suggest that the reionization of the voids, which cover most of the volume, is complete \[Miralda-Escudé 1998\]. It is important to note that the peak in $`T_0`$ at $`z3`$ does not imply that the temperature of the general IGM reaches a maximum. In fact, our results imply that the temperature of slightly overdense gas ($`\delta 2`$) is almost constant because the slope of the equation of state, $`\gamma `$, is minimum when $`T_0`$ is maximum.

Detailed modeling, probably in the form of large hydrodynamical simulations, which include radiative transfer, is required to see whether the various observational constraints can be fit into a consistent picture. However, the ingredients necessary to explain our discovery of a peak in the temperature of the low density IGM at $`z3`$ are clear even from our crude model: a softer background at high redshift to delay helium reionization and enhanced heating rates compared to the optically thin limit. Once the evolution of helium heating is understood, the measurements of the temperature at higher redshifts can be used to constrain the epoch of hydrogen reionization.

Finally, we would like to note that because of their hard spectrum, quasars tend to ionize helium shortly after hydrogen, although the delay depends on the clumpiness of the IGM \[Madau, Haardt & Rees 1998\]. It may therefore be difficult to postpone the reionization of helium until $`z3`$, if quasars were responsible for reionizing hydrogen at $`z>5`$. Hence the mounting evidence for helium reionization at $`z3`$ suggests that hydrogen was reionized by stars.

## ACKNOWLEDGMENTS

We would like to thank Bob Carswell for letting us use the spectrum of Q1100$``$264 and for helping us with VPFIT. We are also grateful to Sara Ellison for letting us use the HIRES spectra of APM 08279+5255. We thank Martin Haehnelt, Lam Hui, Jordi Miralda-Escudé and Martin Rees for stimulating discussions. JS thanks the Isaac Newton Trust, St. John’s College and PPARC for support, WLWS acknowledges support from NSF under grant AST-9900733 and GE thanks PPARC for the award of a senior fellowship. Research was conducted in cooperation with Silicon Graphics/Cray Research utilising the Origin 2000 supercomputer COSMOS, which is a UK-CCC facility supported by HEFCE and PPARC. This work has been supported by the TMR network on ‘The Formation and Evolution of Galaxies’, funded by the European Commission.