text:
# Infrared Observations During the Secondary Eclipse of HD 209458 b II. Strong Limits on the Infrared Spectrum Near 2.2 µm

## 1 INTRODUCTION

The discovery of the first transiting extrasolar planet, HD 209458 b, (Charbonneau et al., 2000; Henry et al., 2000) has led to several new observations designed to characterize the physical properties of the planet. These observations have provided a determination of the planetary and stellar radii and the true planetary mass (Brown et al., 2001; Cody & Sasselov, 2002). Furthermore, given the geometry of the orbit, we are now beginning to learn about the structure of the planet’s atmosphere. The atmosphere was first probed conclusively by Charbonneau et al. (2002), who reported a detection of the sodium doublet in transmission as the planet crossed in front of the star. Although typical models account for the strong stellar irradiation and estimate that the effective temperature of the planet is $`T`$ 1100–1800 K (e.g., Seager & Sasselov, 1998; Charbonneau et al., 2000), no measurements of the temperature are available from actual observations of the planet. An attempt to detect the reflected starlight from the planet during *secondary eclipse* (i.e., the time when the planet disappears behind the star) has recently been reported (Kenworthy & Hinz, 2003), but these observations were performed in the visible region and were limited by instrumental and atmospheric effects. Pioneering infrared observations (Lucas & Roche, 2002) did not achieve sufficient sensitivity to detect realistic planetary models of HD 209458 b and similar extrasolar planet systems. The infrared observations reported here, however, have sufficient sensitivity to detect the thermal emission spectrum from HD 209458 b at the level predicted by several models for the planetary atmosphere.

In paper I (Richardson et al., 2003), we applied the method of occultation spectroscopy observations of the secondary eclipse from the Very Large Telescope (VLT), and using this technique over a narrow bandpass near 3.6 µm, we were able to place limits on the abundance of methane in the planetary atmosphere. The most stringent limit applied only to an exceptionally clear model atmosphere, and only for certain values of the eclipse timing, which is uncertain by as much as $`30`$ minutes, given the 1$`\sigma `$ error in the eccentricity of the orbit. In this paper, we report a further attempt to detect the secondary eclipse of HD 209458 b using the SpeX instrument at the NASA Infrared Telescope Facility (IRTF). We took a different approach for these observations; we used a broader wavelength range, in order to look for the infrared continuum of the planet near 2.2 µm. By observing the combined light from the star and planet, we search for a change in the shape of the spectrum as the planet disappears behind the star and later reemerges. The nature of this technique makes our observations quite sensitive to the temperature gradient in the planetary atmosphere.

## 2 OBSERVATIONS

We obtained a total of four nights of data from the SpeX instrument (Rayner et al., 2003) at the NASA IRTF located on Mauna Kea in Hawaii. These are UT 19, 20, 26, and 27 September 2001, and secondary eclipse events were predicted for UT 20 and 27 September.

The SpeX instrument is a cross-dispersed echelle spectrometer capable of imaging and spectroscopy in the wavelength region between 0.8 and 5.5 µm with low to moderate spectral resolution ($`\lambda /\mathrm{\Delta }\lambda 2500`$). We operated the instrument in 1.9 to 4.2 µm mode using the 0.5 arcsec slit, giving a spectral resolution of 1500; this mode gives nearly-continuous wavelength coverage in this region, over six orders of the echelle. The weather was excellent (low water vapor), and the seeing conditions were good—less than 1.0 arcsec on all four nights. Both eclipses occurred within an hour of the transit of the star across the local meridian, meaning that the eclipse was observed with the object directly overhead.

The use of a comparison star is an important aspect of our observational technique, and we selected HD 210483, a close spectral match to HD 209458 and nearby in the sky (see Table 1 in Paper I). The comparison star was observed immediately following each observation of the HD 209458 system. We nodded the telescope between the ‘a’ and ‘b’ positions on the slit to remove the terrestrial atmospheric background, and we recorded spectra in the sequence ‘abba.’ A given spectrum at either slit position was recorded with an integration time of 6 seconds per coadd and a total of 3 coadds; since the two objects were nearly the same visual magnitude, we used the same integration time for observing both objects. Typically we would record two successive ‘abba’ sets of HD 209458 and then switch to the comparison star for two more ‘abba’ sets, giving approximately a 50% duty cycle. It required only about 15 seconds to switch between the two stars. Over the four nights of observations, we obtained 1036 individual spectra of HD 209458 and 868 spectra of HD 210483, with a typical signal-to-noise ratio of $`100`$ near 2.2 µm. Finally, we take calibration spectra approximately once per hour. To account for flexure and other changes over the night, important for the large SpeX instrument, we record flats using an IR continuum lamp, and we record arc lamp spectra for wavelength calibration.

## 3 ANALYSIS

In this section we discuss the details of the analysis process used to interpret the observations. This process consists of rejecting outlying points (or ‘hot pixels’), extracting the spectra from the raw data, and subtracting a suitable comparison spectrum to remove the telluric features. We search for changes in the resulting ‘difference spectrum’ that are synchronous with the secondary eclipse and are therefore due to the spectrum of the planet.

### 3.1 Spectral Extraction

The first step in the analysis was the rejection of energetic particle events and other intermittently bad pixels from the raw data. This was accomplished using a median filter over a set of eight raw data frames, corresponding to two ‘abba’ sets.

For extracting the spectra from the raw data frames, we used the Spextool software package (version 3.0), written by Mike Cushing and Bill Vacca for analysis of SpeX data (Vacca et al., 2003). The IDL widgets allow the user to enter an ‘ab’ pair of raw images, with the corresponding calibration files (flats and arcs). Given an individual ‘ab’ pair, the program extracts the sky background using the region of the slit between the ‘a’ and ‘b’ object positions, and removes this from the extracted ‘a’ and ‘b’ spectra. Further, it checks for the residual sky background in the pair-subtracted image, and removes it if necessary. Currently, the program outputs a standard extraction of the ‘a’ and ‘b’ spectra; that is, it performs a sum over the spatial dimension to calculate the flux at each wavelength point. Spextool also performs a wavelength calibration of the extracted spectra using the arc frames (recorded by taking an image of an Argon lamp) as well as tabulated information on telluric lines and comparing with the extracted spectra. Finally, Spextool corrects for slight detector non-linearities in the SpeX array.

### 3.2 Corrections to the Extracted Spectra

The spectra of each object for a given night are first interpolated onto a constant and uniform wavelength scale. At this point we perform a second quality control check to find and correct any remaining outlying points due to uncorrected ‘hot pixels.’ For each wavelength point we perform a median filter over the time series of all values of that point for a given object during the night. A point is considered to be an outlier if the difference between the value of the point and the median value is greater than $`4\sigma `$, and it is then replaced by the median value. This process effectively removes outlying points from the spectra, but it is only necessary to correct about 0.3–0.4% of the total number of points for each object.

Next we remove from the stack any spectra that are clearly discrepant. The rejected spectra correspond in most cases to observations that were listed in our notes as questionable, either because of short-term changes in terrestrial atmospheric conditions or poor focus of the telescope. The stacks of spectra for both objects now contain only those spectra that will later be used in the calculation of the difference spectra.

With the individual spectra (again, of a given object) now interpolated onto the same wavelength scale and cleaned of outlying points, we then correct the data at each wavelength point for air mass. We fit a line to the log of the intensity as a function of air mass, as described in Richardson et al. (2003, Equation 3), although in this case, we correct to air mass of unity. In order to ensure that we do not remove the effect of the planet from the HD 209458 spectra, we experimented with using the air mass correction calculated for the comparison star to correct the HD 209458 spectra. We found that this did not have a large effect on the resulting difference spectrum, showing that correcting the HD 209458 and HD 210483 spectra separately did not introduce a bias in the resulting residual spectrum. An example of extracted spectra (only order 4, $``$2.0–2.4 µm) of HD 209458 and HD 210483 after the air mass correction is shown in Figure 1.

We then calculate an average air mass-corrected spectrum of each object for that particular night. By comparing the average spectra of HD 209458 and HD 210483 we can identify stellar lines; at this resolution ($`\lambda /\mathrm{\Delta }\lambda =1500`$), only strong stellar lines are important. Because of the large difference in the radial velocities of the two stars, $`|\mathrm{\Delta }v|55.2`$ km s<sup>-1</sup> (Nidever et al., 2002; Wilson, 1953), a given stellar line will appear at a slightly different wavelength in one star relative to the other. These show up clearly by subtracting the two average spectra. We remove them by linearly interpolating between the points on either side of the line. In comparing HD 209458 and HD 210483, we found a total of 16 stellar lines, identified as Mg and Si, as well as a single H line and a single Na line.

### 3.3 Difference Spectra

With the groups of spectra cleaned and adjusted for air mass, we are now ready to consider individual spectra. First, we remove the known stellar lines from the individual spectra in the same way that they were removed from the average spectra. Note that at this point in the analysis, the two sets of spectra (corresponding to the two stars) have been extracted and interpolated onto separate wavelength scales. Next, an optimum shift value in wavelength is calculated (for both the average as well as the individual spectra separately) by minimizing the standard deviation of the difference in intensity between a target and a normalized comparison spectrum for each order separately. Using the calculated wavelength shift, we can then interpolate each comparison spectrum onto the wavelength scale of each HD 209458 spectrum. This ensures that the telluric absorption features line up in the average and individual spectra of both objects.

Next, we calculate difference spectra by comparing individual spectra of HD 209458 and HD 210483. As described, we record spectra of HD 209458, typically in groups of 8, or two ‘abba’ sets, and then switch to the comparison star. The subsequent observations of the comparison star are recorded at nearly the same air mass as those of HD 209458. We then calculate the normalized ‘difference spectrum’ $`d_i`$ from

$$d_i=\frac{(t_if_ic_i)}{g_i\overline{c}}$$
(1)
where $`t_i`$ represents a single spectrum of HD 209458, $`c_i`$ is the corresponding comparison star spectrum, $`f_i`$ and $`g_i`$ are normalization factors, and $`\overline{c}`$ is the average comparison star spectrum for a given night of observations. We calculate a difference spectrum for each individual spectrum of HD 209458 by subtracting the corresponding subsequent comparison spectrum; for example, the first ‘a’ spectrum of HD 209458 in a given ‘abba’ set would be compared to the first ‘a’ spectrum in the subsequent comparison set. The normalization factors are used to place the comparison spectra (average and individual) on the same scale as the spectra of HD 209458. Although the two stars are nearly the same brightness, we nevertheless expect that they will not have exactly the same intensity levels, as seen in the example in Figure 1, due to slit losses, variable atmospheric absorption, and changes throughout the night. Both normalization factors are calculated for the entire spectrum, rather than for each order independently, to ensure that any overall slope in the planetary spectrum is not removed. The factors are computed by enforcing the condition that the total of all orders in the comparison spectrum is equal to the total of all orders in the spectrum of HD 209458, as in

$$f_i=\frac{_\lambda t_i}{_\lambda c_i}$$
(2)
and similarly for $`g_i`$.

This process effectively and completely removes the terrestrial absorption lines that dominate the object spectra. Furthermore, the process also has two desirable side-effects. First, it removes any time-varying changes in the detector response that may not have been corrected by the flat-fielding. Second, it removes many of the variations in the continuum and line absorption that were not corrected by the fit to air mass, because these variations may also be time-variable. Note also that the resulting normalized difference spectrum contains the candidate planetary spectrum in flux units relative to the stellar spectrum, or the ‘contrast’ (Sudarsky et al., 2003).

### 3.4 Averaging over Wavelength

As described above, the individual difference spectra $`d_i`$ have been calculated by subtracting an individual (i.e., at a given slit position ‘a’ or ‘b’) HD 210483 spectrum from the corresponding individual HD 209458 spectrum, as indicated by Equation 1. We then proceed to average the individual difference spectra $`d_i`$ over wavelength in order to improve the signal-to-noise ratio. We separate the spectra into sections, or ‘bins’, and average the points within each bin. The bins are defined specifically for each region, such that a bin boundary does not fall on a telluric or stellar feature. This ensures that we do not introduce a bias in the bin average. The width of each bin is therefore variable but is typically $`0.025`$ µm.

In looking at the series of binned difference spectra, we see that most of them exhibit an overall slope. The variation in intensity is a slowly-varying function of wavelength, amounting to typically 1–5% from 2–4 µm; a slight slope is evident in the example difference spectrum from UT 27 September shown in Figure 2. We attribute this effect to image motion on the slit due to guiding errors, coupled with the $`\lambda ^{0.2}`$ wavelength dependence of seeing (Linfield et al., 2001). Our interpretation is based on three considerations: 1) the extreme slopes are clearly related to image quality, as judged by the slit losses through total intensity; 2) we constructed a simple model of the process and found that we expect to see such slopes, based on the wavelength dependence of seeing combined with guiding errors of plausible magnitude; and 3) the IRTF staff also see the effect (M. Cushing, private communication, 2002). We have investigated several methods for removing this baseline effect, but we found that our results for the short wavelength region (2–2.5 µm) are robust and do not change significantly based on whether a correction is made or on the details of the correction. This makes sense, given that the flux peak in the planetary spectrum predicted by most models (e.g., Sudarsky et al., 2003) is a relatively sharp peak and therefore insensitive to the removal of a low-order polynomial over the entire 2–4 µm region. We therefore present our results with no attempt to remove these gradual slopes.

Finally, the uncertainties have been carefully propagated through the entire analysis. At the beginning of the analysis, after the spectra have been extracted from the raw frames using Spextool, we calculate the noise level in each data set. (Each data set consists of a time-series of spectra of a given object for a given night of observations.) We take the error in each point to be the standard deviation in the time series at each wavelength point independently. This error value is propagated through the calculation of the difference spectra and the wavelength binning using the standard error propagation formulae. Furthermore, we tested our analysis technique by adding a synthetic planetary signal to the data and verifying that our algorithms could extract it with the correct amplitude.

### 3.5 Fit to Eclipse Curve

The final step in the analysis is the comparison of the processed stack of residual difference spectra to the eclipse timing curve. The eclipse curve is calculated based on the predicted time of center of secondary eclipse (corrected for the Earth-Sun light travel time), using known values of the period and time of center eclipse from Schultz et al. (2003) ($`P=3.5247542\pm 0.0000044`$ days and $`T_c=2452223.896173\pm 0.000086`$ HJD) and assuming a circular orbit (eccentricity $`e=0`$). We also note that updated ephemeris data ($`P=3.5247501\pm 0.0000004`$ days and $`T_c=2452618.66888\pm 0.00010`$ HJD) is available from Wittenmyer et al. (2002). The eclipse curve is constructed in a simple way; zero represents during eclipse, unity represents out of eclipse, and we perform a linear interpolation between the ingress and egress points. Then we perform a linear least squares fit of the difference spectra to the eclipse curve, at each wavelength bin. The amplitude from the least squares fit effectively produces a final spectrum that represents the out-of-eclipse spectra minus the in-eclipse spectra, and therefore represents the candidate planetary spectrum.

Finally, we consider the possibility of a shift in the predicted time of center of secondary eclipse, caused by a non-zero value of the orbital eccentricity (Charbonneau, 2003). The predicted times are based on the assumption of zero eccentricity, which implies that the orbit is circular and that the secondary eclipses occur exactly halfway between primary eclipses. However, the current value of the eccentricity, based on a single-planet fit to the Doppler velocity measurements, is $`e=0.0281\pm 0.0120`$ with $`\omega =69.78\pm 1.2^{}`$ (G. Laughlin, private communication, 2003). If this eccentricity is taken at face value, it implies a shift in the timing of the secondary eclipse of $`\delta t=31.4`$ minutes, a slightly later eclipse. As shown below in Section 5, we find certain models for the planetary atmosphere to be inconsistent with the data at $`e=0`$. We have investigated the effect of changing the time of center eclipse by as much as 120 minutes in either direction. At no eclipse times do we find evidence for the rejected model spectra at full strength in the data.

The final residual spectrum shown in Figure 3 represents the resulting fit to the eclipse curve, averaged for both of the in-eclipse nights, assuming zero orbital eccentricity. We show only the result for the short wavelength region (1.9–2.5 µm). The larger error bars in the 3.0–4.0 µm region (suggested by Figure 2) are caused by the terrestrial background and prevent a conclusive interpretation of the planetary spectrum in this region. Analysis of these long-wavelength data is continuing but is beyond the scope of this paper. Also plotted in Figure 3 is a model spectrum by Sudarsky et al. (2003)—their baseline model for HD 209458 b; this model is excluded by our data, as described in detail in Section 5.

We note that a flat line would be consistent with our data shown in Figure 3, implying that any model that does not exhibit the peak shown in the Sudarsky et al. (2003) baseline model would escape detection in the present analysis. A flat line drawn through the data overlaps 16 of the 27 points within the individual 1$`\sigma `$ error bars, with a reduced chi-squared value of 2.35 (suggesting that the error bars are slightly underestimated and smaller than the scatter in the data). Note that for a normal distribution of errors, 19 of 27 points would be expected to fall within $`\pm 1\sigma `$; therefore, our data are roughly consistent with random deviations from a flat line. We verified that the errors in each wavelength bin are normally distributed. In order to explore the range of models that are rejected by our data, we do some diagnostic model calculations.

## 4 MODEL CALCULATIONS

Our model calculations are intended to help interpret the null result for the planetary spectrum in the region between 2.0–2.5 µm. We have implemented a simple spectral synthesis code to calculate model spectra for HD 209458 b. The algorithm consists of adopting a temperature-pressure profile and computing the emergent thermal emission spectrum. We then scale the temperature-pressure profile in an ad hoc fashion and re-compute the spectrum. Although we do not enforce radiative equilibrium when scaling the profile, the method is nonetheless consistent with our intent, which is to provide a diagnostic of the planetary atmosphere.

The three fiducial temperature-pressure profiles, described below, were used as input to the spectral synthesis code and were calculated in radiative equilibrium using improvements on the method given by Seager & Sasselov (1998) and Seager et al. (2000). We devised a simple way of scaling the fiducial temperature-pressure profiles to see how these scalings changed the resulting spectrum. The prescription we used was

$$T_i^{}=T_i+\gamma (T_iT_0)+T_{\text{offset}},$$
(3)
where $`T_i`$ is the temperature of layer $`i`$, and $`T_0`$ is the boundary temperature, i.e., the temperature at small optical depth. Two parameters control the shape of the synthetic profile: $`\gamma `$ determines the temperature gradient, or atmospheric heating, and $`T_{\text{offset}}`$ is used to adjust the overall temperature of the profile. We characterize the shape of the resulting profile using the temperature at optical depth unity for $`\lambda =2.15`$ µm (the center of the bandpass) and subtracting the boundary temperature:

$$\mathrm{\Delta }T=T(\tau _{\nu _c}=1)T_0.$$
(4)
Armed with a temperature-pressure profile, we continue by calculating the continuous opacity due to collision-induced absorption (CIA) for H<sub>2</sub>-H<sub>2</sub> (Borysow, 2002) and H<sub>2</sub>-He (Jørgensen et al., 2000). We include line opacities for H<sub>2</sub>O, CO, and CH<sub>4</sub>. The water lines are taken from the extensive line database calculated by Partridge & Schwenke (1997) and compiled by Kurucz,<sup>1</sup><sup>1</sup>1http://kurucz.harvard.edu/molecules/h2o/ the CO lines are from Goorvitch (1994), and the CH<sub>4</sub> lines were obtained from HITRAN (Rothman et al., 1998). The relative mixing ratios of the three species were determined using the simple thermochemical equilibrium formulae provided by Burrows & Sharp (1999). We then compute emergent intensity as the formal solution to the radiative transfer equation, with the source function equal to the Planck function, from

$$\mu \frac{dI_\nu }{d\tau }=I_\nu B_\nu (T).$$
(5)
We then solve for $`I_\nu `$ and calculate the flux density by integrating over $`\mu `$, as in

$$F_\nu =_\mu _\tau B_\nu (T(\tau ))e^{\tau /\mu }𝑑\tau 𝑑\mu .$$
(6)
Finally, we calculate the contrast ratio by dividing by a Kurucz model (Kurucz, 1992) for HD 209458. The parameters used in the Kurucz model were $`T=6000`$ K, $`\mathrm{log}g=4.25`$ (cgs), and \[Fe/H\]$`=0.0`$, as calculated by Cody & Sasselov (2002) using the parameters in Mazeh et al. (2000).

The assumption that the source function is given by the Planck function amounts to ignoring scattering processes in the model atmosphere. Only clouds and aerosols will produce significant scattering at this wavelength. Any attempt to specify the distribution of clouds in our suite of diagnostic models would be problematic, since the structure of these models has been perturbed in an ad hoc fashion. Fortunately, the amplitude of the flux peak at this wavelength is determined primarily by the wavelength distribution of CO and H<sub>2</sub>O opacity and the temperature profile of the atmosphere.

It is useful at this point to clarify our usage of the Sudarsky et al. (2003) baseline model. We use their tabulated flux values<sup>2</sup><sup>2</sup>2http://zenith.as.arizona.edu/$``$burrows/ for the planet, divided by the Kurucz model described above. We also multiply by the area ratio $`0.016`$, obtained by taking $`R_p=1.42R_J`$ (Cody & Sasselov, 2002) and $`R_{}=1.146R_{{\scriptscriptstyle }}`$. This yields approximately twice their tabulated contrast values, because their contrast is phase-averaged (Sudarsky, private communication, 2002). Using their temperature-pressure profile as input, we calculate a flux peak of similar but slightly larger magnitude with our simplified spectral synthesis code, and the differences could be explained by the fact that we ignore their cloud opacities in our calculation. Note that our result is based on the Sudarsky et al. (2003) baseline flux divided by the Kurucz model, *not* on our re-calculation of the spectrum.

The assumed redistribution of incident stellar radiation is also different among the models, and we consider this factor in comparing with our data. Recalling the equation for the effective temperature of the planet (Guillot et al., 1996), we have

$$T_{\text{eff}}=T_{}\sqrt{\frac{R_{}}{2D}}\left[f(1A)\right]^{1/4},$$
(7)
where $`T_{}`$ is the effective temperature of the star, $`R_{}`$ is the stellar radius, $`D`$ is the orbital radius of the planet, and $`A`$ is the Bond albedo. The parameter $`f`$ quantifies the redistribution of incident stellar radiation; the $`f=1`$ case represents the situation in which the incident radiation is evenly redistributed over the entire planet, while the $`f=2`$ case represents the one in which the incident radiation is completely absorbed and re-emitted on the day side of the planet. The Sudarsky et al. (2003) baseline model assumes $`f=2`$ (or equivalently, $`f=0.5`$ using their definition of $`f`$ (Saumon et al., 2002)). As noted, our first case for comparing with the observational results is the baseline Sudarsky et al. (2003) flux divided by the Kurucz model. We choose three other cases as fiducial models for comparison with the data. These have been selected in order to test a wide range of the possible physical conditions in the atmosphere. We choose two $`f=1`$ models based on the formulation by Seager & Sasselov (1998); Seager et al. (2000), one with no clouds and one with cloud opacity due to MgSiO<sub>3</sub>, Fe, and Al<sub>2</sub>O<sub>3</sub>. Since the Sudarsky et al. (2003) baseline model considers cloud opacities and assumes $`f=2`$, we further select a cloudless model with $`f=2`$ (based on the Seager & Sasselov (1998); Seager et al. (2000) formulation). The temperature profiles of these four cases are shown in Figure 4, and the results for the four flux peaks (the Sudarsky et al. (2003) baseline result as well as the cloudless and cloudy fiducial cases) are shown in Figure 5.

## 5 RESULTS AND INTERPRETATION

We now fit these four modeled flux peaks to the data using linear least squares, as in

$$r_i=a+bm_i,$$
(8)
where $`r_i`$ represents the resulting difference spectrum (e.g., Figure 3), $`m_i`$ represents the contrast ratio from a given model, $`i`$ represents the wavelength bin, and $`a`$ and $`b`$ are the coefficients of the fit. The slope of this fit, $`b`$, which we call the ‘model amplitude,’ therefore represents the ‘amount’ of the model signal that appears in the data. A model amplitude of unity would mean that the model clearly fits the data (depending on the errors), while a value of zero represents no correlation between the model and the data. Note that the $`y`$-offset between the model and the data is irrelevant, since we are searching only for the evidence of the shape of a given model within the data. The results for the four fiducial models are shown in Table 1. For each model, we indicate the model amplitude $`b`$, the uncertainty in $`b`$, and the reduced chi-squared of the fit. We further list the reduced chi-squared for the case in which we *force* the model amplitude to be unity; that is, we assume that the model appears in the data at full strength. For all four fiducial models, the slope of the fit is roughly zero, suggesting that none of them fit the data. However, the rejection is much stronger for the both $`f=2`$ models, as seen by the reduced chi-squared values for forcing the model amplitude to be unity. These high values indicate that it is extremely unlikely that the $`f=2`$ models fit the data. Given that the observations and analysis are highly differential, we are furthermore confident that a peak in the planetary spectrum is not being fortuitously canceled by some systematic error in the data.

The least squares fit described above is difficult to interpret for models with small amplitude flux peaks, as in some of our scaled models (see Equations 3 and 4). It is also desirable to explore another way to quantify the degree to which the four specific models are rejected by the data. We use a Monte Carlo technique to accomplish this. To quantify the statistical significance of the rejection of any particular model, we constructed ‘fake data’ by adding normally-distributed noise to the modeled contrast values, with $`\sigma `$ equal to the size of the error bar in each bin of Figure 3. We created 10<sup>5</sup> fake data sets for each model and then fit each one to the original model using linear least squares (in the same way in which the real data were fit to the model to calculate the model amplitude). Next, we determine the number of resulting model amplitudes which are greater than the value obtained by fitting the original model to the actual data; this is a liberal criterion for concluding that the model is consistent with a given fake data set. The percentage of fake data sets consistent with a given model, using this criterion, is the ‘detection efficiency’, $`Q`$. A high value of $`Q`$, say $`>99\%`$, implies that a given model consistently produced at least a minimal signature in the fake data. The results for the four specific models discussed thus far are shown in Table 2, and are also plotted in Figure 6 (filled symbols). A model with a high detection efficiency using the Monte Carlo technique implies that it is detectable. Thus, the two $`f=2`$ models—having model amplitudes near zero and high detection efficiencies—are again concluded to be rejected. These results are consistent with the chi-squared analysis described above.

We tested our suite of scaled models, with the prescription given by Equation 3, using this technique. The parameters and results of these models, as well as the model amplitudes, are shown in Table 3 for the cloudless $`f=1`$ cases and in Table 4 for the cloudy $`f=1`$ cases. The results of $`Q`$ vs. $`\mathrm{\Delta }T`$ are plotted in Figure 6 for both sets of models. We reject a model if it has a value of $`Q>99\%`$ and also has a model amplitude near zero. Looking at Tables 3 and 4, the rejected models correspond to the steeper T/P profiles, with $`\gamma 1.25`$. Thus, we can see that the temperature profiles with a higher temperature gradient, causing a larger peak near 2.2 µm in the emergent spectrum, can be eliminated as inconsistent with our observations. The more isothermal profiles, on the other hand, are consistent with the data. Figure 6 suggests that only models with $`\mathrm{\Delta }T1100`$ K are consistent with our data. Although the rejection limit also depends on $`T_{\text{offset}}`$ from Equation 3, it is strongly dependent on $`\mathrm{\Delta }T`$.

## 6 CONCLUSION

Our results are extremely sensitive to the temperature-pressure profile, which is as expected, given that the Planck function varies strongly with temperature in this wavelength region. Besides the shape of the profile, there are a variety of other possible processes in the planetary atmosphere that could render the planet undetectable with these observations. For example, the broad peak in the spectrum near 2.2 µm is caused by the large number of weak water lines on either side of this feature. Absorption by these lines lowers the continuum level on either side, thus creating the peak. One obvious process by which the peak would not appear is a lowering of the water abundance in the planetary atmosphere. We tested this idea and found that only when the water abundance was lowered by a factor of 10 did the peak become small enough to escape detection using our method. This value for the water abundance ($`9.8\times 10^5`$) seems unrealistically small.

We suggest three ways in which the 2.2 µm flux peak would be lowered sufficiently to be consistent with our observations. The first is the efficient redistribution of the incident radiation over the entire planet, which can be achieved through atmospheric dynamics, i.e., winds. Depending on how efficiently the incident radiation is redistributed, the effect could be a temperature asymmetry between the ‘day’ and ‘night’ sides of the planet. Showman & Guillot (2002) suggest that the day-night temperature difference may be as much as 300 K; even a difference of this magnitude would lead to strong winds of 1 km s<sup>-1</sup> or more. The $`f=2`$ situation would tend to steepen the T/P profile in the atmosphere, since the strong stellar irradiation would have to be absorbed relatively deep in the atmosphere and re-radiated. The two specific $`f=2`$ models considered above, which assume that the heat is absorbed and re-radiated only on the day side, are inconsistent with our data (as shown in Table 1). This suggests that the large temperature asymmetry scenario may not be viable and that some sort of transport mechanism may be at work to move the incident energy from the hot to the cold side. We caution that we have observed only two secondary eclipses of HD 209458 b; recent simulations by Cho et al. (2003) indicate that the temperature asymmetry between the two hemispheres may be time-variable. This effect can be so strong as to render the day side to be colder than the night side at times.

A second mechanism by which the 2.2 µm flux peak could be lowered is additional opacity in the upper atmosphere, which would flatten the T/P profile and thus reduce the relative strength of the peak. For example, the T/P profiles shown by Barman et al. (2001) are nearly isothermal, and will certainly produce 2 µm spectra consistent with our observed limit. Barman (2003, private communication) attributes this behavior to their treatment of the opacity, specifically the inclusion of strong gaseous opacity due to metals, which differs from the opacity treatment in other models. However, consistency with our data does not necessarily indicate that a given model is correct, and we leave it to the theorists to sort out the differences among competing models. We suggest that a high-altitude cloud (first suggested by Seager & Sasselov (2000) and favored by Charbonneau et al. (2002) to explain their low sodium result) would also make the atmosphere more isothermal if the cloud is sufficiently absorbing and vertically extended. Note, however, that the Charbonneau et al. (2002) sodium result requires placing such a cloud at the limb of the planet, whereas our results require the clouds to be broadly distributed on the day side.

A high cloud might also affect the emergent spectrum through a third mechanism. For a sufficiently high cloud, the atmosphere above the cloud deck would be thin, giving a small column density in which a spectral signature could be formed. Very little incident radiation would penetrate an optically thick cloud, meaning that the temperature structure would be dominated by the blackbody reflection or emission of the cloud itself. The emergent spectrum would be determined by the small part of the atmosphere above the cloud deck, rather than the details of the T/P profile.

Finally, we note that our results have been obtained using the 3-meter IRTF, indicating that valuable observations of extrasolar planets, leading to new information on the atmospheres of these objects, can be conducted using modest-aperture telescopes from the ground, given favorable observing conditions and careful analysis.

This paper is based upon observations at the IRTF, supported by the National Aeronautics and Space Administration under Cooperative Agreement no. NCC 5-538 issued through the Office of Space Science, Planetary Astronomy Program. We thank Dave Sudarsky for providing us with his temperature-pressure profile and for clarifying discussions of these results. We thank the IRTF staff, in particular Mike Cushing for his lengthy discussions of the Spextool extraction software and John Rayner for his insights into the SpeX instrument. Finally, we thank the referee for helpful comments and suggestions. Jeremy Richardson is supported in part by a NASA Graduate Student Researchers Program Fellowship, funded by the Office of Space Science at NASA Headquarters (grant number NGT5-50273). Other aspects of this research were supported by the NASA Origins of Solar Systems Program. Sara Seager is supported by the Carnegie Institution of Washington.


format: text

text:
# Effects of Electroweak Instantons In High-Energy Neutrino Telescopes

## I Introduction

Instantons are classical solutions of non-Abelian gauge theories in Euclidean space-time instant , that represent tunneling transitions between topologically inequivalent vacua. The processes induced by electroweak instantons could be extremely interesting since they violate baryon+lepton number $`(B+L)`$ conservation thooft . The transition rate for the tunneling processes are exponentially suppressed at low energies compared to the the energy barrier separating the vacua, referred as the “sphaleron” energy ($`E_{sph}\pi M_W/\alpha _\mathrm{w}8`$ TeV) sphaleron . However, at high temperatures highT , or high energies highE , the transition rate may be unsuppressed. It has thus been of great interest to explore the impact of the sphalerons on generating the baryon number asymmetry of the Universe baryon , and to ask if these effects can be observed in high-energy reactions such as at future colliders sphaleron ; pheno .

It is a challenge from a theoretical point of view to reliably estimate the two-body scattering rate for the instanton-induced processes near and above $`E_{sph}`$ Mattis:1991bj . In view of the close analogy of QCD, parton scattering amplitudes induced by electroweak instantons have been calculated using a perturbation method in the instanton background, with the saddle-point approximation state . It was argued that, due to the very large pre-exponential factor in the cross section formula prefactor

$$\sigma _0\frac{(2\pi /\alpha _\mathrm{w})^{7/2}}{M_W^2}5.3\times 10^3\mathrm{mb},$$
(1)
the scattering cross section may be sufficiently enhanced at parton center-of-mass energies of about $`E4\pi M_W/\alpha _\mathrm{w}30`$ TeV, leading to possibly observable effects at a Very Large Hadron Collider (VLHC) and at future ultra-high energy cosmic ray experiments Fodor:2003bn . Recently, a new calculation based on a generalized semi-classical approach has been performed to determine the exponent function rubakov . Their numerical results demonstrated a severe exponential suppression, extended to energies as high as $`30\times E_{sph}250`$ TeV, under the simplification of the $`S`$-wave dominance. Beside the phenomenological and observational excitement motivating searches for electroweak instanton effects, it would be of significant importance to obtain a quantitative description of such interactions and to test the different theoretical approaches, such as the two mentioned above. An initial comparison has just been made prefactor  as we were completing this work.

We are approaching the era of neutrino astrophysics in which it will become possible to observe neutrinos of extra-galactic origin, at energies of PeV to EeV scales. Present experiment, such as AMANDA amanda , Baikal baikal  and RICE rice  have successfully measured the atmospheric neutrino spectrum to energies up to $`10\mathrm{TeV}`$ and have placed limits on any diffuse high-energy neutrino flux beyond this energy. IceCube icecube , with a cubic kilometer effective volume, will be sensitive to charged and neutral current neutrino interactions from $``$TeV to the highest observed energies. Furthermore, planned experiments such as ANITA anita  and EUSO/OWL euso  will have even greater sensitivities to ultra-high energy neutrinos. For a review of high-energy neutrino astronomy, see Ref. review .

In this Letter, we study in detail the instanton or sphaleron-induced processes in high-energy neutrino telescopes, following some early proposals duncan . Neutrino telescopes can measure neutrino scattering cross sections at energies far greater than those accessible to colliders. At energies around $`E_\nu `$100 TeV, the Earth becomes opaque to neutrinos with the standard model (SM) cross sections. Taking advantage of this feature, comparing the flux of neutrinos observed traveling through the Earth (up-going events) to the flux from above the detector (down-going), the neutrino-nucleon cross section can be inferred measure . More generally, the zenith angle distribution of events in an underground detector can, given a sufficient flux, reveal this cross section. Satellite-based experiments and air shower arrays cannot perform this task as effectively, as they are not designed to observe neutrino events over as great an angular distribution, or at as great a range of energies. A large volume, underground detector is ideal for this task. For this reason, we focus on the IceCube detector, which is presently under construction and will be completed well before other post-LHC collider projects can be realized.

In our analysis, we consider both of the calculations for the cross section of instanton-induced neutrino-nucleon interactions described earlier. For the calculation by Ringwald, we perform the numerical evaluation by closely following the formalism in Ref. state . As for the calculation by Bezrukov et al. (BLRRT), only the exponential function is given rubakov . We have thus assumed a constant pre-factor for the cross section given by Eq. (1), multiplied by the exponential function parameterized from rubakov . We note that this may have been a very crude estimate for the cross section since the full pre-factor as in state  is energy-dependent and it may decrease at higher energies.

The cross sections predicted by these calculations are shown in Fig. 1 by the dotted and dashed curves for the Ringwald’s and BLRRT calculations, respectively. As expected, the cross section grows sharply above the threshold, with a $`\nu N`$ center-of-mass energy about 30 TeV for the Ringwald’s calculation, and as high as about 250 TeV for the BLRRT calculation. For comparison, the standard model neutrino-nucleon total cross section of neutral and charged currents uhecr2  is also shown (solid line), where the three active flavors of neutrinos have been included and the CTEQ-5 parton distribution functions cteq  are used. Also illustrated as horizontal dashed lines are the cross sections which correspond to interaction lengths equal to the vertical down-going and horizontal depths of IceCube (D Depth and H Depth). It is in the range between these lines that distinctive features in the angular distribution of down-going events appear. For a neutrino flux $`\mathrm{\Phi }`$, the number of events $`N_\nu `$ observed as hadronic or electromagnetic showers in a neutrino detector of effective volume, $`V`$, is given by the convolution over energy of the quantity $`V\times \mathrm{\Phi }\times n\times \sigma _\nu `$. Here $`n`$ is the density of the target that interacts with a neutrino with cross section $`\sigma _\nu `$ to produce a shower.

In Fig. 2, the energy spectrum of down-going shower events predicted in the IceCube experiment is shown. In the Ringwald’s calculation (dotted line), as the center-of-mass neutrino-nucleon energy exceeds the sphaleron energy threshold near 30 TeV, the number of events increases dramatically above the standard model prediction. Even farther above this energy, however, more of the neutrinos are absorbed in the ice before reaching the detector and the event rate is suppressed. This drastic “bump” structure in the spectrum indicates the sharply enhanced cross section at the sphaleron energy threshold position. The peak of this bump occurs at the associated neutrino energy and is mainly generated by charged current electron neutrino interactions. The “shoulder” slightly to the left of the bump is from neutral and charged current interactions which generate showers less energetic than the incident neutrino. These features occur at considerably higher energies for cross sections found using the calculation by Bezrukov et.al. rubakov . In this approach, due to the exponential suppression of instanton-induced interactions well above the sphaleron energy, observations of these interactions will be considerably more difficult to study in neutrino telescopes, or colliders.

The left frame of Fig. 2 considers a flux of neutrinos equal to the upper bound found by Waxman and Bahcall wbbound . This choice of flux represents neutrinos from compact engines, such as gamma ray bursts or hadronic blazars. The limit of Waxman and Bahcall is somewhat controversial and it has been argued that larger high-energy neutrino fluxes may be possible from such sources mpr . The right frame considers the neutrino flux from the interactions of ultra-high energy protons with the cosmic microwave background, called the cosmogenic neutrino flux cosmogenic . We have used the cosmogenic flux as calculated in Ref. cosmogenic2 . The cosmogenic neutrino flux can be reliably calculated from the observed flux of ultra-high energy protons. Therefore, the choice of this flux is quite conservative. It is interesting to note that at EeV energies, the Waxman-Bahcall bound is not far above the conservative cosmogenic prediction. For this reason, in our discussion the detected neutrino flux plays a secondary role.

Figure 3 shows the zenith angle distribution of showers above 1 EeV in a kilometer-scale detector. Up-going events correspond to $`1<\mathrm{cos}\theta _{\mathrm{zenith}}<0`$, whereas down-going events correspond to $`0<\mathrm{cos}\theta _{\mathrm{zenith}}<1`$. For standard model interactions, the distribution (solid curve) is nearly flat for down-going events, and essentially no up-going events occur due to very efficient neutrino absorption by the Earth at these energies. For models with larger cross sections from instanton-induced interactions, vertical down-going events become more frequent, producing more events near $`\mathrm{cos}\theta _{\mathrm{zenith}}1`$. At zenith angles near the horizon, $`\mathrm{cos}\theta _{\mathrm{zenith}}0`$, more of the neutrinos are absorbed and the rate can be suppressed. We have chosen the shower threshold of 1 EeV to optimize the angular effects for the case of Ringwald’s calculation.

Another interesting characteristic feature of instanton-induced processes is the large multiplicity of final state particles and the violation of $`B+L`$. The basic operators involving quark and lepton fields are of the form $`(qqq\mathrm{})^{n_g}`$ thooft , where $`n_g=3`$ is the number of fermion generations. It has been argued that the processes involving multiple gauge bosons and Higgs bosons, like $`(qqq\mathrm{})^{n_g}W^nH^m`$, can be significantly enhanced highE . A typical neutrino-induced event could thus be

$$\nu _eu\overline{d}\overline{d}+\overline{c}\overline{c}\overline{s}\mu ^++\overline{t}\overline{t}\overline{b}\tau ^++nW+mH.$$
(2)
With both quarks and leptons of all three generations involved simultaneously at the primary production, this type of events should look quite unique. It is, however, difficult to predict how the events would exactly look like in the IceCube detector given the fact that the particles are highly collimated and will be challenging to separate.

Shown in Table 1 are the event rates predicted in a kilometer-scale high-energy neutrino telescope, such as IceCube. Rates are shown for the standard model charged plus neutral current (CC+NC), as well as for the calculations by Ringwald and by BLRRT (including CC+NC). Rates are shown for three (shower) energy ranges, chosen to illustrate the features in the energy spectrum associated with each approach (see Fig. 2). For the Ringwald’s calculation, the event rate predicted for the cosmogenic neutrino flux is about $`0.7`$ per year in the narrow range of 0.4 to 2 EeV, a factor of ten above the SM (CC+NC) prediction. This rate is somewhat higher (1.2) for the flux given by the Waxman-Bahcall bound. Either of these cases provide a possible signature for observations over several years. The predictions of BLRRT are considerably more difficult to test, however. In the peaked region of the energy spectrum, at 10 to 20 EeV, only on the order of 0.03 events per year are expected. Although this is much larger than the SM (CC+NC) prediction, it will be very challenging for kilometer-scale instruments to observe this signature. An order of magnitude larger detector in effective volume would be needed to probe such a scenario.

It has been proposed that ultra-high energy neutrinos undergoing instanton-induced interactions may have generated many of the observed cosmic ray events above the GZK cutoff Fodor:2003bn . Although it would be extremely interesting to establish this interpretation, it will be difficult to determine with confidence that this is the case with air shower experiments. Neutrinos with mb scale cross sections have similar experimental signatures to protons and, therefore, will be difficult to distinguish.

Future satellite-based cosmic ray experiments, such as EUSO/OWL euso , may be able to observe similar numbers of ultra-high energy neutrino events compared to IceCube halzen20 . Much like ground based cosmic ray experiments, EUSO/OWL does not have the phenomenological advantages of a deeply buried neutrino telescope described in this paper. Such experiments will have a more difficult time identifying neutrino-nucleon cross section enhancements.

In summary, the SM electroweak instantons may provide observable signatures in kilometer-scale high-energy neutrino telescopes. The large deviations of the neutrino-nucleon cross section from the charged and neutral current prediction found in these calculations provide distinctive features in both the spectrum and angular distribution of events in the detector. These variations reflect both the increased probability of a neutrino interacting in the detector and the increased probability of being absorbed in the Earth. Another qualitative feature would be presented by the particle content of the events (all three generations of quarks and leptons with $`B+L`$ violation).

The calculation of Ringwald’s should be unambiguously tested by IceCube before post-LHC colliders can be realized. The prediction of BLRRT will be considerably more difficult to probe due to the suppression of the rate in the relevant energy region.

Acknowledgments: This work was supported in part by a DOE grant No. DE-FG02-95ER40896 and in part by the Wisconsin Alumni Research Foundation.


text:
# History of Hydrogen Reionization in the Cold Dark Matter Model

## 1 INTRODUCTION

The epoch of reionization of the universe started with the emission of the first ionizing photons into the intergalactic medium (hereafter, IGM), and ended when all the low-density regions of the universe were ionized. During this period of time, the ionizing photons reaching the IGM had to be sufficient to ionize every atom in the universe and to balance any recombinations in the IGM (in addition, many more ionizing photons may have been emitted in regions of dense, self-shielded gas, which were locally absorbed and did not contribute to ionizing the IGM). At the end of reionization, the mean free path of ionizing photons increased to scales much larger than the typical separation of the collapsed halos hosting the sources. Regions of dense, self-shielded gas (observed as Lyman limit systems) shrank in size as the ionizing background increased in intensity, due to the growing mean free path (Miralda-Escudé, Haehnelt, & Rees 2000). If the fraction of baryons in regions of intermediate density which were ionized as the mean free path increased was small, the rise in the background intensity could be relatively fast, as seen in the numerical simulations of Gnedin (2000).

Precisely this rapid change of the ionizing background intensity has been inferred by Fan et al. (2002) to occur at $`z6`$, by analyzing the change of the mean transmitted flux of the Ly$`\alpha `$ forest in the highest redshift quasars known. The transmitted flux apparently drops abruptly to very low levels near $`z=6`$. This observation most likely implies that the end of reionization occurred, in fact, at $`z6`$, although as a larger number of quasars at $`z>6`$ are discovered, more detailed investigations of this issue will need to be done to understand precisely how the mean free path grew. In particular, the number of small gaps of transmitted flux, such as the one seen in a $`z=6.37`$ quasar by White et al. (2003), provides information on the size of the H II regions before they overlapped (Miralda-Escudé 1998; Barkana 2002).

The results of the Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2003 and references therein) have strengthened the observational support of the Cold Dark Matter (CDM) scenario in a flat model that contains a vacuum energy component (in addition to the ordinary baryonic matter and dark matter). The model parameters are being measured with increasing accuracy via combined observations of the cosmic microwave background (CMB), galaxy clustering, weak lensing, and the Ly$`\alpha `$ forest. The theory then makes an increasingly reliable prediction for the abundance of collapsed dark matter halos as a function of mass and redshift. It is in these halos that the process of radiative gas dissipation and star formation can take place, leading to the emission of ionizing radiation. However, the uncertainties in the efficiency of star formation, the stellar initial mass function (IMF), the emission due to gas accretion into black holes (be they in X-ray binaries or active galactic nuclei), and the fraction of the ionizing photons that escape into the IGM from the dense regions in which they are emitted, prevent a clear prediction for the history of reionization.

Even if the time at which reionization ends is determined to be at $`z6`$ by the Gunn-Peterson trough (Gunn & Peterson 1965) observations, there can be a long time interval in which the universe was partially ionized, with cosmological H II regions around the sources filling only some fraction of the volume. A very good observational probe for the characteristic time at which the bulk of the baryons in the universe were ionized is the electron scattering optical depth to the CMB, $`\tau _e`$. This has been measured for the first time by WMAP (Kogut et al. 2003; Spergel et al. 2003) using the large-scale polarization-temperature anisotropy correlation. At present, the result still has a large error, both observational and model-dependent: Kogut et al. find a model-independent value of $`\tau _e=0.16\pm 0.04`$ from the polarization-temperature correlation of the WMAP data, but Spergel et al. find $`\tau _e=0.17\pm 0.06`$ from a fit of a CDM model with a running spectral index to a combined data set of CMB, galaxy clustering and Ly$`\alpha `$ forest observations, and $`\tau _e=0.12\pm 0.06`$ when the fit is forced to the usual CDM model with a primordial power-law spectrum. A value of $`\tau _e=0.17`$ would imply that most baryons were reionized as early as $`z20`$.

The announcement of the measurement of $`\tau _e`$ by WMAP has led to a large number of papers discussing the possibility of this early reionization in the CDM model (Whyithe & Loeb 2003; Haiman & Holder 2003; Holder et al. 2003; Ciardi, Ferrara, & White 2003; Somerville & Livio 2003; Sokasian et al. 2003; Cen 2003; Chiu, Fan, & Ostriker 2003). These authors have generally found that, if one is willing to assume a high efficiency in the production of ionizing photons per unit mass in the low-mass halos collapsing at early times, then early reionization and a high optical depth is possible in CDM models. Naturally, the lower the power spectrum amplitude on small scales, the higher the efficiency that is required, and some models can be ruled out based on their small power on small scales, such as the warm dark matter model (Spergel et al. 2003). But to a large extent, the ability of the CDM model to accommodate a wide range for the epoch at which the bulk of the IGM was reionized simply reflects our ignorance of the efficiency in producing ionizing photons that escape.

However, a reasonable question we can ask is how the required efficiency compares to any observationally determined values of the rate at which ionizing photons have been emitted over the history of the universe, taking into account the predicted abundances of halos of different masses that form at each redshift. In this paper, we adopt a very simple model of the ionizing emissivity from CDM halos, but which still captures the essential physical ingredient for the process of reionization: as the medium is ionized, the gas is heated and the Jeans mass is raised, implying that some low-mass halos where gas can collapse when the medium is neutral are unable to continue to accrete gas after the medium is ionized. We divide halos into two populations, the “high-mass” ones (which will be called Population A) that are able to accrete gas and form new stars irrespective of the state of the intergalactic medium, and the “low-mass” ones (Population B) that can form stars and emit ionizing radiation only when the medium is neutral. Low-mass halos collapse first in CDM, but as they ionize the medium around them the rate of star formation can be suppressed, until the high-mass halos start to collapse (Couchman & Rees 1986; Shapiro, Giroux, & Babul 1994). A similar separation of halo classes was made by Haiman & Holder (2003); our Population A corresponds to their Type II and Type Ia halos (Haiman & Holder separated these halos into the ones that can or cannot cool by atomic cooling after they have accreted the gas; we will not distinguish between these two classes in our more simple model), and our Population B corresponds to their Type Ib population.

Our calibration of the efficiency to produce ionizing photons is based on observations of the Ly$`\alpha `$ forest during the post-reionization era. The mean transmitted flux of the Ly$`\alpha `$ forest, combined with CDM simulations of Ly$`\alpha `$ forest spectra and measurements of the baryon density of the universe, allows us to infer the intensity of the ionizing background; and the abundance of Lyman limit systems tells us the mean free path of ionizing photons. The ratio of the background intensity to the mean free path then yields the volume-averaged emissivity. This emissivity was inferred to be, at most, seven ionizing photons per baryon and per Hubble time at $`z=4`$ by Miralda-Escudé (2003). Remaining uncertainties on the value of this emissivity are related mainly to the temperature of the IGM and the CDM power spectrum amplitude. At $`z>4`$, the abundance of Lyman limit systems becomes highly uncertain, so the emissivity cannot be inferred. The key question for reionization is, of course, how this emissivity changes with increasing redshift. If the comoving emissivity stays constant at $`z>4`$, the optical depth $`\tau _e`$ is not larger than $`0.09`$ (Miralda-Escudé 2003). In CDM models, the fraction of mass in halos that can form stars decreases with redshift at $`z>4`$, so to obtain $`\tau _e>0.09`$ (which requires an increase of the emissivity with redshift), the low-mass halos (or Population B) formed at $`z4`$ must have a much higher efficiency to emit ionizing radiation than the halos of higher mass at $`z4`$.

The question we are addressing in this paper is: how much higher does this efficiency have to be in order to produce a certain value of the optical depth $`\tau _e`$, and what constraints are imposed by requiring reionization to end at $`z=6`$? We describe our model for the efficiency to emit ionizing radiation in §2. Results are presented in §3, and a general discussion is given in §4.

## 2 MODELS

### 2.1 CDM Power Spectrum

We use three CDM models to predict the halo abundances. All the models assume flat space with a cosmological constant, with $`\mathrm{\Omega }_{m0}=0.27`$, $`\mathrm{\Omega }_{b0}=0.044`$, $`h=0.71`$. The three models vary in the normalization and slope of the power spectrum: Model 1 has $`\sigma _8=0.84`$, $`n_s=0.93`$; Model 2 has $`\sigma _8=0.9`$, $`n_s=0.96`$; and Model 3 has $`\sigma _8=0.9`$, $`n_s=1`$ (here, $`\sigma _8`$ is the linear rms mass fluctuation at present on spheres of radius $`8h^1\mathrm{Mpc}`$, and $`n_s`$ is the primordial scalar spectral index). The power spectrum is computed with the fitting formula of Eisenstein & Hu (1999). Model 1 has the smallest amplitude of fluctuations on the small scales on which the first halos that emit ionizing radiation collapse ($`0.03`$ to 1 Mpc), and Model 3 has the largest amplitude. Model 3, with the largest power on small scales, has halos forming at the highest redshift, and therefore predicts an earlier start of reionization and higher optical depth $`\tau _e`$ than the other models for a fixed ionizing emissivity per unit mass.

Model 1 has the same parameters as the best fit model found by WMAP (Bennett et al. 2003; Spergel et al. 2003) to the combined data set of CMB observations, the 2dF Galaxy Redshift Survey (Colless et al. 2001), and Ly$`\alpha `$ forest power spectrum data (Croft et al. 2002). However, we do not include the running spectral index introduced in Spergel et al. (inclusion of the running spectral index would reduce the power on small scales). Model 2 has higher values of $`\sigma _8`$ and $`n_s`$ by $`1.5`$ and $`1`$ standard deviations, respectively, from the best fit values of Spergel et al. Model 2 also has a slightly larger power than the best fit CDM model without a running spectral index found by Spergel et al. (which has parameters $`\sigma _8=0.8`$ and $`n_s=0.96`$). Finally, Model 3 is very close to the best fit of Spergel et al. to CMB data alone, and is the same model that was adopted in various recent studies on the reionization history (Ciardi et al. 2003; Somerville & Livio 2003; Sokasian et al. 2003; Cen 2003).

Measurement of the Ly$`\alpha `$ forest has allowed a determination of the power spectrum amplitude on the smallest scales. We therefore compare the amplitude of our models with the determination of the Ly$`\alpha `$ forest power spectrum by Croft et al. (2002), who find a normalization at redshift $`z=2.72`$ of

$$\mathrm{\Delta }^2(k_p)\frac{k_p^3}{2\pi ^2}P(k_p)=0.74_{0.16}^{+0.20},$$
(1)
where $`H(z)^1k_p=0.03(\mathrm{km}/\mathrm{s})^1`$. The normalization of our Models (1, 2, 3) are $`\mathrm{\Delta }^2(k_p)=(0.85,1.06,1.27)`$, respectively. Hence, Models 2 and 3 have a power spectrum amplitude higher than the Croft et al. (2002) measurement by $`1.6\sigma `$ and $`2.6\sigma `$, respectively. The high amplitude of Model 3 would seem to be ruled out by the result of Croft et al.; note, however, that the quoted error on this measurement does not include some potentially important systematic effects, as discussed in Croft et al., and the amplitude can be higher if the mean Ly$`\alpha `$ transmitted flux is higher than the value used by Croft et al. (see Seljak, McDonald, & Makarov 2003).

### 2.2 Halo Abundances

We use the Sheth & Tormen (1999) prescription for computing the abundances of halos as a function of mass and redshift. Halos are divided into two populations: those which do not emit any ionizing photons after the medium around them has been reionized (Population B), and those which are massive enough to continue forming stars and producing ionizing photons after reionization (Population A). In this paper we place the division between the two populations at a velocity dispersion $`\sigma _0=35\mathrm{km}\mathrm{s}^1`$, based on the suppression of gas infall and cooling in smaller halos caused by photoionization (Efstathiou 1992; Thoul & Weinberg 1996; Navarro & Steinmetz 1997). The lower limit to the velocity dispersion of the Population B of low-mass halos is placed at $`\sigma _{min}=3.68\mathrm{km}\mathrm{s}^1`$, because molecular hydrogen cooling below a temperature of $`2000`$ K is ineffective (Yoshida et al. 2003).

The fraction of baryons, $`F_{A,B}`$, which occupy each of the halo populations is shown in Figure 1 for Models 1 to 3. We also show in Figure 2 the difference between using the Sheth-Tormen and the Press-Schechter (Press & Schechter 1974) mass functions, for Model 1. The Sheth-Tormen mass function predicts a slightly larger number of objects at the highest redshifts.

### 2.3 Emissivity

For Population A, i.e., halos with $`\sigma >\sigma _0=35\mathrm{km}\mathrm{s}^1`$, the ionizing photon emissivity per baryon is assumed to be of the form

$$ϵ_A(z)=F_A(z)\frac{ϵ_4}{F_A(4)}\left(\frac{1+z}{5}\right)^\alpha ,$$
(2)
where $`F_A(z)`$ is the fraction of baryons in Population A halos at redshift $`z`$, and $`ϵ_4`$ is the total emissivity at $`z=4`$, which is assumed to be equal to 7 ionizing photons per baryon and per Hubble time at $`z=4`$, the upper limit obtained in Miralda-Escudé (2003). We will use two values of $`\alpha `$ in this paper: for $`\alpha =0`$ the energy emitted per unit physical time and per unit mass remains constant, and for $`\alpha =1.5`$ the total energy emitted per Hubble time at redshift $`z`$ remains constant. Note that gas cooling is inefficient in very massive halos when the virialized gas temperature is $`>10^7K`$, but at $`z=4`$ the abundance of these halos, which arise from $`>4\sigma `$ fluctuations, are negligibly small, so our assumption of a constant emissivity per unit mass for Population A is reasonable.

For Population B, we allow an emissivity per unit mass different from the one in Population A halos, parameterizing it as

$$ϵ(\sigma )=\left[1+\gamma \mathrm{ln}\frac{\sigma _0}{\sigma }\right]\frac{ϵ_4}{F_A(4)}\left(\frac{1+z}{5}\right)^\alpha ,$$
(3)
where $`\gamma `$ is a scaling factor used to adjust the variation in efficiency between the low-mass and high-mass halos. This results in a smooth variation of the emissivity with $`\sigma `$. The Population B emissivity is then

$$ϵ_B(z)=_{\sigma _{min}}^{\sigma _0}f_B(\sigma ,z)ϵ(\sigma )𝑑\sigma ,$$
(4)
where $`f_B`$ is the fraction of mass in halos in the interval $`d\sigma `$ around $`\sigma `$, and we set $`\sigma _{min}=3.68\mathrm{km}\mathrm{s}^1`$, corresponding to $`2000`$ K for the mean particle mass of the primordial atomic gas.

### 2.4 Analytical Framework

The equation for the ionization of the intergalactic medium is (Madau, Haardt, & Rees 1999; Miralda-Escudé 2003)

$`{\displaystyle \frac{dy}{dt}}=`$ $`ϵRy`$ (5)
$`=`$ $`ϵ_A+ϵ_B(1y)Ry,`$ (6)
where $`y`$ is the volume fraction of the IGM that is ionized at any given time $`t`$, and $`R`$ is the mean number of recombinations per baryon in the ionized regions per unit time. This equation makes the approximation that the photon mean free path is very short, so emitted photons are instantaneously used to ionize atoms. It also hides all the complicating physics of reionization in the averaged recombination rate $`R`$, which depends on the clumping factor of the photoionized gas. Note that Population B halos emit ionizing radiation only when they are located in the atomic medium, so the emissivity $`ϵ_B`$ is multiplied by $`1y`$. We assume a clumping factor of unity throughout this paper, computing the recombination coefficient at $`T=10^4`$ K, which gives $`R=1.16\times 10^{17}\mathrm{s}^1`$ at $`z=4`$. The solution of equation (6) is

$`y(t)=`$ $`{\displaystyle _{t_i}^t}\mathrm{exp}\left\{{\displaystyle _t^{}^t}\left[R(t^{\prime \prime })+ϵ_B(t^{\prime \prime })\right]𝑑t^{\prime \prime }\right\}`$ (7)
$`\times \left[ϵ_B(t^{})+ϵ_A(t^{})\right]dt^{},`$
where $`t_i`$ is an arbitrarily chosen initial time (before any appreciable halo formation).

### 2.5 Optical Depth

The Thomson scattering optical depth, $`\tau _e`$, is

$$\tau _e(z)=\frac{n_e\sigma _ec}{H_0}_0^z\frac{F_i(1+z)^2y(z)dz}{[\mathrm{\Omega }_{\mathrm{\Lambda }0}+\mathrm{\Omega }_{m0}(1+z)^3]^{1/2}},$$
(8)
where $`\mathrm{\Omega }_{\mathrm{\Lambda }0}=1\mathrm{\Omega }_{m0}`$, $`n_e`$ is the comoving electron number density, $`\sigma _e`$ is the electron cross section, and $`F_i`$ is the fraction of the baryons in the ionized regions that are actually in the ionized IGM (with the rest being in self-shielded atomic or molecular clouds, or in stars). We assume $`F_i=0.9`$ (as in Miralda-Escudé 2003) for $`z<6`$, but $`F_i=1`$ at $`z>6`$; lowering $`F_i`$ at $`z>6`$ would demand starting reionization at higher redshift to obtain a fixed optical depth, increasing the required emissivity from the first star-forming halos (note that equation (5) needs to be modified if $`F_i<1`$ during reionization). For our value of $`\mathrm{\Omega }_bh^2=0.022`$, the electron number density is $`n_e=2.7\times 10^7(1f_YY)\mathrm{cm}^3`$, where $`Y=0.24`$ is the helium abundance by mass, and $`f_Y=(0.5,0.75)`$ when the helium is (doubly, singly) ionized. We assume that helium becomes doubly ionized at $`z=3.5`$ (e.g., Heap et al. 2000), and that it is singly ionized at higher redshift wherever hydrogen is ionized. This results in an optical depth up to $`z=6`$ of $`\tau _e(z=6)=0.032`$.

## 3 RESULTS

We now present the results on the fraction of the volume that is ionized as a function of redshift, $`y(z)`$, with the emissivity that arises from our two populations of halos given by equations (2) and (3). To summarize the emissivity model, the Population A (high-mass halos) emissivity per unit mass is fixed at $`z=4`$ and is allowed to vary as $`(1+z)^\alpha `$, while the Population B (low-mass halo) emissivity per unit mass can also be increased with the parameter $`\gamma `$.

Figure 3 shows $`y(z)`$ for each of the three CDM models and for both $`\alpha =0`$ and $`\alpha =1.5`$. The case $`\gamma =0`$, where the emissivity per unit mass from Population A and B is the same, is shown as solid lines. The cases $`\gamma =1`$, $`10`$ and $`100`$ are shown as dotted, dashed, and dash-dot lines, respectively. The models with high $`\gamma `$ values have an early reionization as the first Population B halos collapse; then, $`y(z)`$ grows slowly when it approaches unity because of the decrease in the emission from Population B as most of the IGM is ionized, and finally $`y`$ reaches unity when the Population A halos start forming and emitting radiation in the ionized regions of the universe.

The optical depth to electron scattering for all the models shown is listed in Table 1.

The emissivity as a function of redshift is shown in Figure 4 for the same models. At redshift just above $`4`$, the emissivity declines with redshift because the abundance of Population A halos (which are the only ones contributing when the medium is entirely ionized) drops rapidly. The emissivity increases again at redshifts higher than the end of reionization, when Population B halos contribute. Models with higher $`\gamma `$ can obviously have a higher emissivity at very high redshift, but their emissivity declines at intermediate redshifts when $`y`$ is close to unity and the emission from Population B halos is suppressed, as assumed in equation (3). The models with higher $`\gamma `$ result in higher optical depths because the IGM is reionized earlier.

## 4 DISCUSSION

For a fixed CDM halo population, the most simple model for the emissivity is to say that the rate of photon emission per unit mass is constant for all halos in which cooling can take place, i.e., with $`\sigma >\sigma _{min}`$. This model, using the emissivity per unit mass calibrated with the observational determination at $`z=4`$ and a clumping factor of one, predicts a redshift at which reionization ends of $`5`$ to $`7`$, in very good agreement with the observation of the appearance of the Gunn-Peterson trough at $`z=6`$. The emissivity history in this model shows a decline with redshift from $`z=4`$ to $`6`$, then a slight increase up to $`z10`$ as more Population B halos contribute, and then a rapid decline beyond $`z`$ 10 to 15 as the number of all halos decreases.

Obviously, the reionization history should also depend on the clumping factor affecting the recombination rate. We have assumed a clumping factor of unity, and a higher clumping factor would delay the end of reionization. A moderate clumping factor of 2 or 3 would not greatly affect the result, because the recombination time at the mean baryonic density is equal to the Hubble time at $`z=6`$ and the increase in $`y(z)`$ is relatively fast for the $`\gamma =0`$ models (solid lines in Fig. 3; this is simply due to the rapid increase in the mass fraction in collapsed halos shown in Fig. 1), so the number of recombinations is not very large. For example, a clumping factor of 3 for Model 1 with $`\alpha =1.5`$ (panel b) delays the end of reionization from $`z=6`$ to $`z=5`$. As discussed by Miralda-Escudé et al. (2000), and Miralda-Escudé (2003), the clumping factor is unlikely to be large during reionization. Therefore, we have the remarkable result that the simplest model of a constant emissivity per unit mass in CDM halos, calibrated with the observed emissivity at $`z=4`$, roughly predicts the right epoch for the end of reionization.

This same model predicts an optical depth to Thomson scattering, $`\tau _e`$, of $`0.05`$ to $`0.09`$, depending on $`\alpha `$ and on the allowed variations on the power spectrum parameters determining the amplitude of fluctuations on the small scales of the first halos. Although this is smaller than the value found by the WMAP mission in a CDM model fit according to Spergel et al. (2003), $`\tau _e=0.17\pm 0.06`$, it is not ruled out with the present errorbar: the upper value of $`\tau _e=0.09`$ of our $`\gamma =0`$ models is consistent with this measurement at the $`1.3\sigma `$ level. The likelihood function of this range of $`\tau _e`$ shown in Figure 8 of Spergel et al. (2003) is also seen to be not very small. Moreover, if the running spectral index in the CDM model is not allowed, the best fit then yields $`\tau _e=0.12\pm 0.06`$ (see Table 7 of Spergel et al.). The smaller errorbar obtained by the analysis in Kogut et al. (2003) is puzzling, because the addition of observational data and the use of a parameterized model instead of a model-independent analysis should not result in an increased error.

If a high value of $`\tau _e`$ is confirmed, then the ionizing emission efficiency in the earliest low-mass halos to form in CDM is required to be much higher than in the sources present in the more massive halos at $`z=4`$. The required factor by which the efficiency must increase can be inferred from Table 1. The increase in emissivity per unit mass is $`(1+2\gamma )(1+z)^\alpha `$ (see eq. 3), for a typical halo contributing to the emissivity at high redshift with $`\sigma 5\mathrm{km}\mathrm{s}^1`$. For the case $`\alpha =0`$, the factor $`1+2\gamma `$ gives the increase in the rate of emission per unit mass, while for the case $`\alpha =1.5`$, the factor $`1+2\gamma `$ gives the increase in the total energy per unit mass emitted in ionizing photons, taking into account the shorter Hubble time at higher redshift. Interpolating from Table 1, we find that to reach $`\tau _e=0.17`$, the total ionizing energy per unit mass emitted by Population B halos at $`z15`$ needs to be larger by a factor $`1+2\gamma 10`$ compared to the Population A halos at $`z=4`$, and the rate of emission needs to be higher by a factor $`100`$. For Model 1 (more consistent with the Croft et al. 2002 measurement of the Ly$`\alpha `$ forest power spectrum), the energy emitted per unit mass needs to increase by a factor $`25`$, and the emission rate per unit mass needs to increase by $`300`$.

Are these very high efficiencies for producing ionizing photons in the earliest low-mass halos physically realistic? The emissivity of Population A halos at $`z=4`$ is about 60 photons per baryon per Hubble time, as obtained by dividing the mean emissivity $`ϵ_4`$ of 7 photons per baryon per Hubble time (see §2.3) by the fraction of baryons in Population A halos at $`z=4`$ (see Fig. 1). From the emission ratios mentioned above required to produce $`\tau _e=0.17`$, Population B halos at $`z20`$ would need to emit $`1500`$ and $`600`$ photons per Hubble time for Models 1, and 3, respectively. Since in these models with high optical depth the universe is reionized at $`z20`$ over a time interval corresponding to $`\mathrm{\Delta }z5`$ (see Fig. 3), or one quarter of a Hubble time, the total emission during this time interval coming from the first halos needs to be about $`400`$ and $`150`$ photons per baryon in Population B halos, for Models 1 and 3.

To compare this to the number of photons that can be produced by star formation, we use the models of Leitherer et al. (1999; see their Fig. 78) to find that for a Salpeter IMF with solar metallicity in the stellar mass range $`1M_{}<M<100M_{}`$, 6000 ionizing photons per baryon are emitted. Most of the ionizing photons are produced by stars with $`M>20M_{}`$, and in the case of a Salpeter IMF with stars forming only in the range $`20M_{}<M<100M_{}`$, 25000 photons per baryon can be emitted. Note that this number is equal to the ratio of the energy obtained by hydrogen fusion (7 MeV) to the average energy of an ionizing photon (20 eV) times the fraction of stellar mass that is fused to helium ($`30`$ %) times the fraction of energy emitted as ionizing photons ($`25`$%). Very massive stars with zero metallicity can further increase this emission to $`10^5`$ photons per baryon (Tumlinson & Shull 2000; Bromm, Kudritzki, & Loeb 2001; Schaerer 2002), because their convective cores during the main-sequence include almost all the mass of the star (zero metallicity stars also have higher temperature than metal-rich ones during the main-sequence, which increases the fraction of energy they emit in ionizing radiation, but also increases the average energy of their ionizing photons). Thus, if $`f_{}`$ is the fraction of the gas in a halo that can be converted to stars with $`M>20M_{}`$, and $`f_{esc}`$ is the fraction of ionizing photons that escape to the IGM, then an optical depth $`\tau _e=0.17`$ requires $`f_{}f_{esc}>(10^{2.4},10^{2.8})`$ for Models (1,3) for zero-metallicity stars, and $`f_{}f_{esc}>(10^{1.8},10^{2.2})`$ for Models (1,3) for metal-enriched stars.

It is clear from these values that an optical depth $`\tau _e=0.17`$ can be produced with physically realistic models. Even $`\tau _e=0.3`$ is possible in Model 3 with $`\gamma =100`$ and $`\alpha =1.5`$, which implies an emission of 3000 ionizing photons per baryon (or $`f_{}f_{esc}=10^{1.5}`$ for zero-metallicity stars) from the first star-forming halos over the redshift range $`25z32`$. However, if a high optical depth is confirmed by future measurements, an explanation for the dramatic change in ionizing emission from the low-mass halos that formed the first stars to the more massive halos at lower redshift will need to be found.

A difficulty in making star formation very efficient in low-mass halos is that their shallow potential wells makes them more vulnerable to lose their gas as a result of stellar winds, photoionization and supernovae. The escape fraction might be higher than in more massive halos at lower redshift, but only after gas has been ionized and expelled from the vicinity of stars, in which case the gas can easily escape from the halo. The typical halo where stars form at $`z20`$ with velocity dispersion $`5\mathrm{km}\mathrm{s}^1`$ and total mass of $`10^6M_{}`$ (e.g., Yoshida et al. 2003) has a gas mass of $`2\times 10^5M_{}`$ within a virialized radius $`200`$ parsecs, and cooling by molecular hydrogen leads to the formation of a dense central gas clump containing only $``$ 1% of the gas mass. The rest of the gas must take more than a halo dynamical time ($`10^{7.5}`$ years) to reach the center and cool before it has any possibility to form stars, by which time supernovae would have taken place. The hydrodynamic simulations of these first star-forming halos (Abel, Bryan, & Norman 2002; Bromm, Coppi, & Larson 2002) suggest that the first stars may all be very massive. If the $`2000M_{}`$ of the core of cooled gas can efficiently form one or a few massive stars, and most of the emitted ionizing photons escape, one may barely be able to achieve the fraction $`f_{}f_{esc}10^{2.5}`$ required for $`\tau _e=0.17`$. After the formation of these stars, the supernovae explosions would expel the remaining gas from the halo (Bromm, Yoshida, & Hernquist 2003).

In halos of higher mass (with $`\sigma `$ of 10 to 30 $`\mathrm{km}\mathrm{s}^1`$), a larger number of supernovae would be required to expel the gas. Nevertheless, the effects of photoionization and supernova explosions should still self-regulate the rate of massive star formation, and it is not clear why the earliest galaxies at high redshift could have higher values of $`f_{}f_{esc}`$ than the more massive galaxies formed at lower redshift. While the mean gas density in galaxies should increase with redshift, increasing gas cooling rates, the shallower potential wells of low-mass halos allow for gas dispersal after a smaller fraction $`f_{}`$ of gas has turned to massive stars.

Another possibility is that black holes formed by the first massive stars are able to accrete gas. Gravitational accretion can produce radiation with an efficiency $`20`$ times greater than stellar fusion (and stars generally fuse only a fraction of their mass). In this case, requiring a fraction as low as $`10^{3.5}`$ of the halo gas to accrete on a central black hole could produce the $`300`$ photons per baryon in halos required for $`\tau _e=0.17`$, depending on the fraction of the accretion energy emitted as ionizing photons (we have assumed 20% here). However, whether one speculates about a high efficiency of massive star formation or of gas accretion into black holes, the question that remains unanswered is why such high efficiency to produce ionizing radiation would occur only in the first halos at high redshift, and not in similar or more massive halos at later times.

If it were nevertheless possible to produce these high emissivities when the first stars formed, a second possible problem of an early reionization might be that the high emissivity leads to a rapid end of reionization, in contradiction with the presence of the Gunn-Peterson trough at $`z=6`$. Our results suggest a way to solve this particular problem: as seen in Figure 3, once the fraction of the IGM that is ionized becomes close to unity, the emissivity can be decreased because star formation in Population B halos is suppressed. This can regulate the rate of emission to a level that just balances recombinations, and it is only when the more massive Population A halos reach a critical abundance that the emissivity can rise again, completing reionization and causing the fast increase in the mean free path of ionizing photons that clears the Gunn-Peterson trough. Thus, in all our models the end of reionization occurs between redshift 5 and 8, even though the models with high $`\gamma `$ can ionize half of all the volume as early as $`z20`$. In the models with the earliest reionization in Figure 3, the redshift for the end of reionization could easily be lowered to fit the observed value of 6 by moderately increasing the clumping factor after $`y`$ is close to unity or increasing the value of $`\sigma _0`$ that separates our Populations A and B.

Note that models with high $`\gamma `$ have an emissivity history that is double peaked: there is a first maximum due to Population B halos, then the emissivity declines as most of the IGM is ionized, and a second maximum is then caused by Population A halos. The possibility of a double-peaked history of star formation due to the effect of reionization raising the Jeans mass was proposed by Couchman & Rees (1986) and Shapiro et al. (1994). However, this does not necessarily imply a “double reionization”, as other authors have argued (Cen 2003; Whyithe & Loeb 2003). The fraction of ionized volume $`y(z)`$ rises monotonically with time, and it is only “stalled” by the suppression of the star formation due to ionization, in a self-regulated process. A double reionization in which the universe recombines again after having been ionized would require a long delay of the feedback ionization effects on the star formation rate. If this delay were present, the suppression of star formation might start to operate only after stars have already emitted more than enough photons to ionize all the universe; later, when the IGM recombines again, star formation might not resume at the same level for some time, allowing for partial recombination to occur. This feedback delay might be caused by a reservoir of gas in Population B halos that can keep forming stars after the surrounding low-density medium is ionized, and by the finite time required for new recombined and cooled gas to accrete in the Population B halos after the IGM recombines again; but it is not clear if this would be sufficient to cause the global instability in the star formation rate that would be required for a double reionization.

Finally, comparing our results with other recent work, we note that other authors (e.g., Whyithe & Loeb 2003; Haiman & Holder 2003; Ciardi et al. 2003; Somerville & Livio 2003; Sokasian et al. 2003; Cen 2003) have generally emphasized that the central value of the WMAP result $`\tau _e0.17`$ can be obtained in CDM models with a physically allowed emissivity of ionizing photons. Although we agree with this conclusion, we note that a drastic change in the emissivity per unit mass of star-forming halos from $`z15`$ to $`z=4`$ is required. For example, Ciardi et al. (2003) find, using our Model 3, that $`\tau _e=0.16`$ can be reached in models where the dominant emission comes from halos corresponding roughly to our Population A, but the emissivity they assume would vastly overpredict the intensity of the ionizing background at $`z=4`$. At the same time, models with no such change in the emissivity yield a lower optical depth that is still within the plausible error range of the WMAP result. Our approach is closest to that of Chiu et al. (2003), who use an ionizing emission efficiency chosen to yield reasonable models for the ionizing background, and our conclusions are generally in agreement with theirs.

## 5 CONCLUSION

The first detection of the Thompson optical depth to the CMB by the WMAP mission introduces a new era in the study of reionization. The observations of the Gunn-Peterson absorption in the highest redshift quasars, which are sensitive to the time when reionization ended and the mean free path of ionizing photons increased, can be confronted with the total ionized column to the surface of last scattering, which is sensitive to the epoch when about half of the baryons in the universe were ionized. The simplest assumption we can make about the emissivity is that it is roughly constant in all halos in which gas can cool. When we calibrate the emissivity to the upper limit allowed by observations at $`z=4`$, we find that the end of reionization is close to $`z=6`$, as suggested by observations, but that the optical depth is smaller than $`0.09`$. To obtain an optical depth as high as $`\tau _e=0.17`$ (the central value of the WMAP result), an increase of the rate of ionizing emission per unit mass in collapsed halos with redshift, of a factor $`100`$ to $`300`$ (depending on the CDM power spectrum model) from $`z=4`$ to $`z20`$, is required.

Further progress on our understanding of the reionization epoch will come, among other things, from an increased accuracy in the measurement of $`\tau _e`$, as well as more robust measurements of the power spectrum on small scales from Ly$`\alpha `$ forest observations. The discovery of quasars at $`z>6`$ and the search for any small fraction of Ly$`\alpha `$ transmitted flux in gaps along their Gunn-Peterson trough (White et al. 2003) may also allow us to test if the fraction of the ionized medium was close to unity for a wide redshift range above $`z=6`$ (as in our models with high $`\gamma `$ in Fig. 3), or if it increased rapidly over a short redshift range. If an optical depth as high as the central value of the WMAP measurement is confirmed, then either a physical explanation for the high efficiency of the ionizing emissivity of the first objects will need to be found, or some more fundamental modification of the CDM model, such as non-Gaussian fluctuations on small scales (e.g., Avelino & Liddle 2003), may be implied.

This work was supported in part by NSF grant NSF-0098515.