text:
# Reionization of Hydrogen and Helium by Early Stars and Quasars

## 1 Introduction

Following cosmological recombination at a redshift $`z10^3`$, the baryonic gas filling up the universe became predominantly neutral. Given that this gas is known to be mostly ionized today, one arrives at two of the major questions in current extragalactic astronomy: (i) when were the cosmic hydrogen and helium re-ionized? and (ii) which sources dominated this reionization process? The answers to both questions are likely to be different for hydrogen and helium (see review by Barkana & Loeb 2001). Recent observations provide preliminary answers to the first question. The absorption spectra of SDSS quasars at $`z6`$ indicate that the neutral fraction of hydrogen increases significantly at $`z6`$ (Becker et al. 2001; Fan et al. 2002), and the UV spectrum of quasars implies that helium is fully ionized only at $`z3`$ (Jacobsen et al 1994; Tytler 1995; Davidsen et al. 1996; Hogan et al. 1997; Reimers et al. 1997; Heap et al. 2000; Kriss et al. 2001; Smette et al. 2002). The latter observations also indicate, through a cross-correlation between the hydrogen and helium forest of absorption lines (Kriss et al. 2001; Smette et al. 2002), that the ionization of helium at $`z3`$ had significant contributions from both quasars and stars. No analogous evidence exists for hydrogen at $`z6`$.

In this paper we make theoretical predictions for the reionization histories of hydrogen and helium. We calculate the contribution from quasars using a model that matches all the existing observational data on the quasar luminosity function at high redshifts (see Wyithe & Loeb 2002 for details). We model the less constrained stellar contribution using two free parameters, which are in turn constrained by existing observational data. Since the recombination times of both hydrogen and helium are shorter than the age of the universe in overdense regions of the intergalactic medium (IGM), it is possible that these species experienced more than one epoch of reionization. One of the goals of our detailed study is to identify the parameter values for which multiple reionization epochs are possible.

The reionization history has important implications for the temperature and polarization anisotropies of the cosmic microwave background (CMB). Anisotropies on scales smaller than the angular size of the horizon at hydrogen reionization are suppressed by a factor $`e^{\tau _{\mathrm{es}}}`$ where $`\tau _{\mathrm{es}}`$ is the line-of-sight optical depth for electron scattering (Hu & White 1997; Haiman & Loeb 1997), and secondary anisotropies are added (Hu 2000 and references therein). Reionization is the primary source of polarization anisotropies on large angular scales (Zaldarriaga & Seljak 1997; Kamionkowski et al. 1997; Hu 2000). Polarization anisotropies might be detected in the near future by MAP on large angular scales (Kaplinghat 2002) or by ground-based experiments on small angular scales. So far, the available data on the temperature anisotropies of the CMB provides an upper limit of $`\tau _{\mathrm{es}}0.18`$ (Wang, Tegmark & Zaldarriaga 2002; Bond et al. 2002). We will use this upper limit to constrain the free parameters of our model.

Throughout the paper we assume density parameters values of $`\mathrm{\Omega }_m=0.35`$ in matter, $`\mathrm{\Omega }_b=0.052`$ in baryons, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$ in a cosmological constant, and a Hubble constant of $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (or equivalently $`h=0.65`$). For calculations of the Press-Schechter (1974) mass function (with the modification of Sheth & Tormen 1999) we assume a primordial power-spectrum with a power-law index $`n=1`$ and the fitting formula to the exact transfer function of Cold Dark Matter, given by Bardeen et al. (1986). Unless otherwise noted we adopt an rms amplitude of $`\sigma _8=0.87`$ for mass density fluctuations in a sphere of radius $`8h^1`$Mpc.

## 2 Reionization in a Clumpy Universe

The simplest estimate of the epoch of reionization is based on the following simple considerations. Given a co-moving density of ionizing photons $`n_\gamma `$ in a homogeneous but clumpy medium of comoving density $`n_0`$ (where the size of the HII region is much larger than the scale length of clumpiness), the evolution of the volume filling factor $`Q_\mathrm{i}`$ of ionized regions is (Haiman & Loeb 1997; Madau et al. 1999; Barkana & Loeb 2001)

$$\frac{dQ_\mathrm{i}}{dz}=\frac{1}{n_0}\frac{dn_\gamma }{dz}\alpha _\mathrm{B}\frac{C}{a^3}Q_\mathrm{i}n_\mathrm{e}\frac{dt}{dz},$$
(1)
where $`\alpha _\mathrm{B}`$ is the case B recombination coefficient, $`a=1/(1+z)`$ is the scale factor, $`n_\mathrm{e}`$ is the comoving electron density, and $`Cn_0^2/n_0^2`$ is the clumping factor. This equation describes statistically the transition from a fully neutral universe to a fully ionized one, and yields reionization redshifts for hydrogen of around 7-12. Large uncertainties arise in both the source term and in the value of the clumping factor (because more rapid recombinations lead to a slower evolution of $`Q_\mathrm{i}`$).

A more realistic description of reionization in a clumpy medium is provided by the model of Miralda-Escude et al. (2000). In what follows, we draw primarily from their prescription and refer the reader to the original paper for a detailed discussion of its motivations and assumptions. The model assumes that reionization progresses rapidly through islands of lower density prior to the overlap of individual cosmological ionized regions. Following overlap, the remaining regions of high density are gradually ionized. It is therefore hypothesized that at any time, regions with gas below some critical overdensity $`\mathrm{\Delta }_\mathrm{i}\rho _i/\rho `$ are ionized while regions of higher density are not. The assumption of homogeneity in equation (1) implies that the volume filling factor equals the mass filling factor. Therefore, within the model of Miralda-Escude et al. (2000) we replace $`dQ_\mathrm{i}/dz`$ by $`dF_\mathrm{M}(\mathrm{\Delta }_\mathrm{i})/dz`$, where

$$F_\mathrm{M}(\mathrm{\Delta }_\mathrm{i})=_0^{\mathrm{\Delta }_\mathrm{i}}𝑑\mathrm{\Delta }P_\mathrm{V}(\mathrm{\Delta })\mathrm{\Delta }$$
(2)
is the fraction of mass in regions with overdensity below $`\mathrm{\Delta }_\mathrm{i}`$, and $`P_\mathrm{V}(\mathrm{\Delta })`$ is the volume weighted probability distribution for $`\mathrm{\Delta }`$. The product $`Q_\mathrm{i}C`$ can be rewritten

$$R(\mathrm{\Delta }_\mathrm{i})Q_\mathrm{i}CQ_\mathrm{i}\frac{\rho ^2}{\rho ^2}=F(\mathrm{\Delta }_\mathrm{i})\frac{_0^{\mathrm{\Delta }_\mathrm{i}}𝑑\mathrm{\Delta }P_\mathrm{V}(\mathrm{\Delta })\mathrm{\Delta }^2}{_0^{\mathrm{\Delta }_\mathrm{i}}𝑑\mathrm{\Delta }P_\mathrm{V}(\mathrm{\Delta })}=_0^{\mathrm{\Delta }_\mathrm{i}}𝑑\mathrm{\Delta }P_\mathrm{V}(\mathrm{\Delta })\mathrm{\Delta }^2,$$
(3)
where $`F(\mathrm{\Delta }_\mathrm{i})`$ is the fraction of the volume with $`\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{i}`$. The quantity $`\alpha _\mathrm{B}(1+z)^3n_\mathrm{e}_0^{\mathrm{\Delta }_\mathrm{i}}𝑑\mathrm{\Delta }P_\mathrm{V}(\mathrm{\Delta })\mathrm{\Delta }^2`$ is therefore the number of recombinations per atom in the IGM per second. Note that the term analogous to the clumping factor is calculated from the volume weighted distribution. The mass fraction $`F_\mathrm{M}(\mathrm{\Delta }_\mathrm{i})`$ (or equivalently $`\mathrm{\Delta }_\mathrm{i}`$) therefore evolves according to the equation

$$\frac{dF_\mathrm{M}(\mathrm{\Delta }_\mathrm{i})}{dz}=\frac{1}{n_0}\frac{dn_\gamma }{dz}\alpha _\mathrm{B}\frac{R(\mathrm{\Delta }_\mathrm{i})}{a^3}n_\mathrm{e}\frac{dt}{dz}.$$
(4)
This equation assumes that all ionizing photons are absorbed shortly after being emitted, so that there is no background ionizing field, and no loss of ionizing photons due to redshift. We therefore implicitly assume that the mean free path of ionizing photons is much smaller than the Hubble length. This should be valid at redshifts not too much smaller than the overlap redshift.

The integration of equation (4) requires knowledge of $`P_\mathrm{V}(\mathrm{\Delta })`$. Miralda-Escude et al. (2000) found that a good fit to the volume weighted probability distribution for the density as seen in N-body simulations has the functional form

$$P_\mathrm{V}(\mathrm{\Delta })d\mathrm{\Delta }=A\mathrm{exp}\left[\frac{(\mathrm{\Delta }^{2/3}C_0)^2}{2(2\delta _0/3)^2}\right]\mathrm{\Delta }^\beta d\mathrm{\Delta },$$
(5)
with $`\delta _0=7.61/(1+z)`$ and $`\beta =2.23`$, 2.35 and 2.48, and $`C_0=0.558`$, 0.599 and 0.611 at $`z=2`$, 3 and 4. At $`z=6`$ they assume $`\beta =2.5`$, which corresponds to the distribution of densities of an isothermal sphere, and solve for $`A`$ and $`C_o`$ by requiring the mass and volume to be normalized to unity. We repeat this procedure to find $`P_\mathrm{V}(\mathrm{\Delta })`$ at higher redshifts. The proportionality of $`\delta _0`$ to the scale factor is expected for the growth of structure in an $`\mathrm{\Omega }_\mathrm{m}=1`$ universe or at high redshift otherwise, and its amplitude should depend on the amplitude of the power-spectrum. The simulations on which the distribution in Miralda-Escude et al. (2000) was based assumed $`\mathrm{\Omega }_m=0.4`$ in matter, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$ in a cosmological constant and $`\sigma _8=0.79`$, close to the values used in this paper.

Equation (4) provides a good description of the evolution of the ionization fraction following the overlap of individual ionized bubbles, because the ionization fronts are exposed to the mean ionizing radiation field. However prior to overlap, the prescription is inadequate, due to the large fluctuations in the intensity of the ionizing radiation. A more accurate model to describe the evolution of the ionized volume prior to overlap was suggested by Miralda-Escude et al. (2000). In our notation the appropriate equation is

$$\frac{d[Q_\mathrm{i}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})]}{dz}=\frac{1}{n^0}\frac{dn_\gamma }{dz}\alpha _\mathrm{B}(1+z)^3R(\mathrm{\Delta }_{\mathrm{crit}})n_\mathrm{e}Q_\mathrm{i}\frac{dt}{dz}.$$
(6)
or

$$\frac{dQ_\mathrm{i}}{dz}=\frac{1}{n^0F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})}\frac{dn_\gamma }{dz}\left[\alpha _\mathrm{B}(1+z)^3R(\mathrm{\Delta }_{\mathrm{crit}})n_\mathrm{e}\frac{dt}{dz}+\frac{dF_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})}{dz}\right]\frac{Q_\mathrm{i}}{F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})}.$$
(7)
In this expression, $`Q_\mathrm{i}`$ is redefined to be the volume filling factor within which all matter at densities below $`\mathrm{\Delta }_{\mathrm{crit}}`$ has been ionized. Within this formalism, the epoch of overlap is precisely defined as the time when $`Q_\mathrm{i}`$ reaches unity. However, we have only a single equation to describe the evolution of two independent quantities $`Q_\mathrm{i}`$ and $`F_\mathrm{M}`$. The relative growth of these depends on the luminosity function and spatial distribution of the sources. The appropriate value of $`\mathrm{\Delta }_{\mathrm{crit}}`$ is set by the mean free path of the ionizing photons. More numerous sources can attain overlap for smaller values of $`\mathrm{\Delta }_{\mathrm{crit}}`$. Assuming $`\mathrm{\Delta }_{\mathrm{crit}}`$ to be constant with redshift, we find that results do not vary much (less than 10% in the optical depth to electron scattering) for values of $`\mathrm{\Delta }_{\mathrm{crit}}`$ ranging from a few to a few tens. At high redshift, these $`\mathrm{\Delta }_{\mathrm{crit}}`$ correspond to mean free paths comparable to the typical separations between galaxies or quasars. We assume $`\mathrm{\Delta }_{\mathrm{crit}}=20`$ (which lies between the values for galaxies and quasars) throughout the remainder of this paper.

## 3 The Reionization of Hydrogen and Helium

Next we describe the joint evolution of the filling factors and ionized mass fractions for hydrogen and helium using generalizations of equations (4) and (7). Let $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ be the filling factors of ionized hydrogen, singly ionized helium and doubly-ionized helium, and $`F_{\mathrm{H}^+}`$, $`F_{\mathrm{He}^+}`$ and $`F_{\mathrm{He}^{++}}`$ be the mass-fractions of ionized hydrogen, singly ionized helium and doubly ionized helium, within the volumes $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ respectively. For later use we also define111Note that $`Q^\mathrm{m}`$ is the mass-filling factor. The mass filling factors for H+ and He++ equal $`Q_{\mathrm{H}^+}F_{\mathrm{H}^+}`$ and $`Q_{\mathrm{He}^{++}}F_{\mathrm{He}^{++}}`$ respectively. Prior to the overlap of HeIII regions, $`Q_{\mathrm{He}^+}^\mathrm{m}=Q_{\mathrm{He}^+}F_{\mathrm{He}^+}`$. However following the overlap of HeIII regions, $`Q_{\mathrm{He}^+}^\mathrm{m}=F_{\mathrm{He}^+}`$. the fractions of ionized mass in the universe $`Q_{\mathrm{H}^+}^\mathrm{m}`$, $`Q_{\mathrm{He}^+}^\mathrm{m}`$ and $`Q_{\mathrm{He}^{++}}^\mathrm{m}`$.

We assign the comoving densities of hydrogen and helium; $`n_\mathrm{H}^0=1.88\times 10^7(\mathrm{\Omega }_\mathrm{b}h^2/0.022)\mathrm{cm}^3`$ and $`n_{\mathrm{He}}^0=0.19\times 10^7(\mathrm{\Omega }_\mathrm{b}h^2/0.022)\mathrm{cm}^3`$, respectively, and adopt the case B recombination coefficients222At high redshift the absorbtion mean-free-path divided by the speed of light is much smaller than the Hubble time. The appropriate recombination coefficient is therefore case B. for $`\mathrm{H}^+`$, $`\mathrm{He}^+`$ and $`\mathrm{He}^{++}`$ at a temperature of $`10^4`$K (Osterbrock 1974; Barkana & Loeb 2001); $`\alpha _\mathrm{B}^{\mathrm{H}^+}=2.6\times 10^{13}\mathrm{cm}^3\mathrm{s}^1`$, $`\alpha _\mathrm{B}^{\mathrm{He}^+}=2.73\times 10^{13}\mathrm{cm}^3\mathrm{s}^1`$ and $`\alpha _\mathrm{B}^{\mathrm{He}^{++}}=13\times 10^{13}\mathrm{cm}^3\mathrm{s}^1`$. The ionizing radiation field is described by the comoving densities of ionizing photons for ionized hydrogen, singly ionized helium and doubly ionized helium, namely $`n_\gamma ^{\mathrm{H}^+}`$, $`n_\gamma ^{\mathrm{He}^+}`$ and $`n_\gamma ^{\mathrm{He}^{++}}`$, respectively. These quantities are computed through integration of the corresponding spectra in the frequency intervals $`3.29\times 10^{15}<\nu <5.94\times 10^{15}`$Hz for HI, $`5.94\times 10^{15}<\nu <1.31\times 10^{16}`$Hz for HeI (He+), and $`\nu >1.31\times 10^{16}`$Hz for HeII (He++). Two classes of sources, quasars and stars, contribute to the ionizing radiation field. We discuss these in turn in the following subsections.

### 3.1 Ionizing Photons from Quasars

To calculate the comoving density of ionizing photons emitted by quasars, $`n_\gamma `$, we integrate over both the quasar luminosity function and over the mean quasar spectrum (specified in Schirber & Bullock 2002). Given the B-band quasar luminosity function $`\mathrm{\Phi }(L_\mathrm{B},z)`$ at redshift $`z`$ in units of $`\mathrm{Mpc}^3L_\mathrm{B}^1`$, we find for example

$$\frac{dn_\gamma ^{\mathrm{H}^+}}{dz}=\frac{dt}{dz}_0^{\mathrm{}}𝑑L_\mathrm{B}\mathrm{\Phi }(L_\mathrm{B},z)_{3.29\times 10^{15}\mathrm{Hz}}^{5.94\times 10^{15}\mathrm{Hz}}𝑑\nu \frac{L_\nu (L_\mathrm{B})}{h_\mathrm{p}\nu },$$
(8)
where $`h_\mathrm{p}`$ is the Planck constant, and $`L_\nu (L_\mathrm{B})`$ is the luminosity in units of $`\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$ of a quasar with a B-band luminosity $`L_\mathrm{B}`$. At redshifts below $`z=2.5`$ we use the luminosity function of Boyle et al. (2000). At $`z>2.5`$ we use the theoretical luminosity function of Wyithe & Loeb (2002); the reader is referred to the original paper for more details. This quasar evolution model successfully describes all known properties of the high redshift quasar luminosity function (Fan et al. 2001a,b), and reproduces measurements of the black-hole/dark matter halo mass relation (Ferrarese 2002) as well as estimates of the quasar duty cycle (Steidel et al. 2002). Since quasars at high redshift appear similar to typical examples at low redshift (Fan et al. 2001b), the luminosity function fully determines the ionizing radiation field. We therefore have almost no freedom in our calculation of the ionizing radiation field from quasars at high redshifts. This allows us to explore the effect of different stellar radiation fields on the reionization history of the universe.

### 3.2 Ionizing Photons from Stars

To compute the comoving densities of stellar ionizing photons, $`n_\gamma `$, we integrate over model spectra for stellar populations in the frequency ranges quoted above to find (for each species) the total number of ionizing photons, $`N_\gamma `$, emitted per baryon incorporated into stars. We find

$$\frac{dn_\gamma ^{\mathrm{H}^+}}{dz}=N_\gamma f_{\mathrm{esc}}f_{}\frac{dF_\mathrm{b}}{dz}n_\mathrm{b},$$
(9)
where $`n_\mathrm{b}`$ is the local number density of baryons, $`f_{\mathrm{esc}}`$ is the escape fraction for the ionizing radiation, $`f_{}`$ is the fraction of baryons within a galaxy that are incorporated in stars, and $`F_\mathrm{b}(z)`$ is the fraction of baryons in the universe at redshift $`z`$ that collapsed and cooled inside galaxies. To compute $`F_\mathrm{b}(z)`$ we require the dark matter collapse fraction in halos $`F_{\mathrm{col}}(z,T_{\mathrm{min}})`$ with virial temperatures larger than $`T_{\mathrm{min}}`$, given by

$$F_{\mathrm{col}}(z,T_{\mathrm{min}})=\frac{1}{\rho _\mathrm{m}}_{M_{\mathrm{min}}(T_{\mathrm{min}})}^{\mathrm{}}𝑑MM\frac{dn_{\mathrm{ps}}}{dM},$$
(10)
where $`\rho _m`$ is the co-moving mass density in the universe, and $`dn_{\mathrm{ps}}/dM`$ is the Press-Schechter (1974) mass function for dark matter halos. In regions of the universe that are neutral, the critical virial temperature (the temperature at which atomic cooling becomes efficient) is $`T_{\mathrm{min}}=10^4`$K (a circular velocity of $`v_{\mathrm{cir}}10`$km s-1). In reionized regions, infall of gas is suppressed below a virial temperature of $`T_{\mathrm{min}}2.5\times 10^5`$K ($`v_{\mathrm{cir}}50`$km s-1, Thoul & Weinberg 1996). Given $`F_{\mathrm{col}}(z,10^4K)`$ and $`F_{\mathrm{col}}(z,2.5\times 10^5K)`$ we can take the baryonic collapse fraction

$$F_\mathrm{b}(z)=Q_{\mathrm{H}^+}^\mathrm{m}F_{\mathrm{col}}(z,2.5\times 10^5\text{K})+(1Q_{\mathrm{H}^+}^\mathrm{m})F_{\mathrm{col}}(z,10^4\text{K}).$$
(11)
This prescription for calculating the stellar ionizing field will be referred to hereafter as case A. Note that $`f_{\mathrm{esc}}`$ could be different for hydrogen and helium; however, we adopt a single value for it assuming that this value is dictated by the geometry of the optically-thick gas that allows limited escape routes which are common for all ionizing photons (Wood & Loeb 2000 and references therein).

We also consider a second evolution for the stellar ionizing field, denoted hereafter as case B. A recent study (Kauffmann et al. 2002) has shown that in a large sample of local galaxies, the ratio $`ϵ=M_{}/M_{\mathrm{halo}}`$ (where $`M_{}`$ and $`M_{\mathrm{halo}}`$ are the total stellar and dark matter halo masses respectively) scales as $`ϵM_{\mathrm{halo}}^{2/3}`$ for $`M_{}<3\times 10^{10}M_{}`$, but is constant for larger stellar masses. Note that the star-formation efficiency is proportional to $`ϵ`$. Since star-formation is thought to be regulated by supernova feedback (Dekel & Silk 1986), the important quantity is the depth of the galactic potential well, or equivalently the halo circular velocity. Using the stellar mass Tully-Fisher relation of Bell & De Jong (2001), we find the threshold circular velocity $`v_{}=176\mathrm{km}\mathrm{s}^1`$ that at $`z=0`$ corresponds to a stellar mass of $`3\times 10^{10}M_{}`$. In this case we define $`f_{}`$ as the star-formation efficiency in galaxies with circular velocities larger than $`v_{}`$, and calculate the effective product of the dark matter collapse fraction and star-formation efficiency

$$f_{}F_{\mathrm{col}}(z)=\frac{1}{\rho _\mathrm{m}}_{M_{\mathrm{min}}(T_{\mathrm{min}})}^{\mathrm{}}𝑑Mϵf_{}M\frac{dn_{\mathrm{ps}}}{dM},$$
(12)
where $`ϵ=1`$ for $`M>M_{\mathrm{halo}}^{}`$ and $`ϵ=\left(M/M_{\mathrm{halo}}^{}\right)^{2/3}`$ for $`M<M_{\mathrm{halo}}^{}`$. The value of $`M_{\mathrm{halo}}^{}(z)`$ is calculated from

$$M_{\mathrm{halo}}^{}(z)=4.3\times 10^{10}h^1\left(\frac{v_{}}{176\mathrm{k}\mathrm{m}\mathrm{s}^1}\right)^3\left[\frac{\mathrm{\Omega }_\mathrm{m}}{\mathrm{\Omega }_\mathrm{m}^\mathrm{z}}\frac{\mathrm{\Delta }_{\mathrm{crit}}}{18\pi ^2}\right]^{1/2}\left(\frac{1+z}{10}\right)^{3/2}M_{},$$
(13)
where $`\mathrm{\Delta }_{\mathrm{crit}}=18\pi ^2+82d39d^2`$, $`d=\mathrm{\Omega }_\mathrm{m}^\mathrm{z}1`$ and $`\mathrm{\Omega }_\mathrm{m}^\mathrm{z}=\mathrm{\Omega }_\mathrm{m}(1+z)^3/[\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]`$. Finally, the value of $`f_{}F_\mathrm{b}`$ required for calculation of $`dn_\gamma /dz`$ may then be computed in analogy with case A, and substituted into equation (9).

In both cases A and B we make the distinction between the ionizing radiation field due to a possible early population of zero-metallicity stars, and the metal enriched stars observed at lower redshifts. It is thought that the primordial initial mass function favored massive stars (Bromm, Copi & Larson 1999, 2001; Abel, Bryan & Norman 2000; Mackey, Bromm, Hernquist 2002). The possible existence of this population is very important for reionization because the spectrum of these stars would result in an order of magnitude more ionizing photons per baryon incorporated into stars (Bromm, Kudritzki, & Loeb 2001). The formation of the very massive stars was suppressed as the material out of which stars form was enriched with metals. The fraction of the ionizing photons produced by metal-free stars depends on several unknown parameters, including the mixing efficiency of metals, the environments in which new stars form, and most importantly, the threshold metallicity above which star formation is transformed from being dominated by massive stars to a Scalo (1998) initial mass function (IMF). The threshold metallicity is believed to be small; Bromm et al. (2001) argue for a threshold $`\frac{Z_{\mathrm{thresh}}}{Z_{}}10^3`$ of the solar metallicity value. The efficiency of mixing of metal enriched outflows from star forming galaxies to the surrounding IGM is even more uncertain; Scannapieco, Ferrara & Madau (2002) find that the mass weighted mean metallicity can reach values greater than $`10^3`$ of the solar value at redshifts as high as 20, and note that the average metallicity scales with star formation efficiency, supernovae rate, and the fraction of supernovae energy that is channeled into outflows. The average metallicity increases roughly exponentially with redshift as it is modulated by the exponential growth in the collapse fraction of baryons at high redshifts.

Since the first metals were produced by supernovae, it is reasonable to suppose that the enrichment of the IGM with metals at redshift $`z`$ is proportional to the mass of baryons that has formed stars by that redshift. The average metallicity at a redshift $`z`$ can then be written as

$$\frac{Z(z)}{Z_{}}=C_{\mathrm{metal}}\frac{Z_{\mathrm{thresh}}}{Z_{}}_{\mathrm{}}^z𝑑z\frac{dF_\mathrm{b}}{dz},$$
(14)
where $`C_{\mathrm{metal}}`$ is a constant. Note that the baryonic collapse fraction $`F_\mathrm{b}`$, which is computed from equation (11) is affected by the reionization history. The formation of metal-rich stars is assumed to be proportional to the average metalicity, with the fraction reaching unity when the average metalicity of the IGM reaches $`\frac{Z_{\mathrm{thresh}}}{Z_{}}`$. Before $`\frac{Z(z)}{Z_{}}`$ reaches $`\frac{Z_{\mathrm{thresh}}}{Z_{}}`$ we may approximate the fraction of star formation at redshift $`z`$ that is in zero-metallicity stars to be $`(1\frac{Z(z)}{Z_{}}/\frac{Z_{\mathrm{thresh}}}{Z_{}})`$, while after $`\frac{Z(z)}{Z_{}}`$ reaches $`\frac{Z_{\mathrm{thresh}}}{Z_{}}`$ we may assume all star formation to have a Scalo (1998) IMF333The situation in reality is more complicated since mixing of metals is incomplete and the formation sites of new stars are correlated with the enriched regions. Our discussion uses the minimum number of free parameters to describe this complicated process..

For a given star formation efficiency and escape fraction, the above prescription yields a characteristic redshift for the build-up of the threshold metal enrichment. The exponential growth in the collapse fraction with redshift implies that the evolution in the mode of star formation resembles a step-function. We therefore define a transition redshift $`z_{\mathrm{tran}}`$ below which metal-rich stars with a Scalo (1998) IMF dominate the production rate of ionizing photons.

Through most of the paper we present results in terms of $`z_{\mathrm{tran}}`$, but return to a justification of this choice in § 7. At redshifts above $`z_{\mathrm{tran}}`$, we use values of $`N_\gamma `$ calculated for each of the three species $`\mathrm{H}^+`$, $`\mathrm{He}^+`$ and $`\mathrm{He}^{++}`$ assuming massive ($`100M_{}`$) zero-metallicity stars and the generic spectrum calculated by Bromm, Kudritzki & Loeb (2001). The resulting values of $`N_\gamma `$ are $`14020`$ for H+, $`25200`$ for $`\mathrm{He}^+`$ and $`4480`$ for $`\mathrm{He}^{++}`$. At lower redshifts we assume metal enriched stars (1/20th solar metallicity) with a Scalo (1998) mass-function, and use spectral information from the stellar population model of Leitherer et al. (1999)444Model spectra of star-forming galaxies were obtained from http://www.stsci.edu/science/starburst99/.. This results in values of $`N_\gamma 3250`$ for H+, $`N_\gamma 1020`$ for $`\mathrm{He}^+`$ and $`N_\gamma 0.2`$ for $`\mathrm{He}^{++}`$.

### 3.3 Constraints on the Evolution of $`\mathrm{H}^+`$, $`\mathrm{He}^+`$ and $`\mathrm{He}^{++}`$ Ionization Fronts

The typical distance that an ionizing photon of neutral helium can travel through neutral hydrogen having a cross-section $`\sigma _\mathrm{H}(\nu )`$ before the optical depth reaches unity is the absorption mean-free-path $`\lambda `$, where

$$\lambda ^1=\frac{_{5.94\times 10^{15}\mathrm{Hz}}^{1.31\times 10^{16}\mathrm{Hz}}𝑑\nu (L_\nu /\nu )\sigma _\mathrm{H}(\nu )n_H^0(1+z)^3}{_{5.94\times 10^{15}\mathrm{Hz}}^{1.31\times 10^{16}\mathrm{Hz}}𝑑\nu L_\nu /\nu }.$$
(15)
We find $`\lambda 4\left[\left(1+z\right)/10\right]^3`$kpc. This distance is much shorter than $`1`$Mpc which is the typical size of an HII region around a $`10^8M_{}`$ halo prior to reionization (Barkana & Loeb 2000). Therefore, at $`z10`$ the $`\mathrm{He}^+`$ propagation front can only lead the $`\mathrm{H}^+`$ propagation front by a small fraction of the size of the HII region, and the excess $`\mathrm{He}^+`$ ionizing photons photo-ionize additional hydrogen. On the other hand, most ionizing photons for hydrogen ($`\nu >5.94\times 10^{15}`$Hz) are not absorbed by helium, and so the expansion of the HII region is not inhibited by the presence of helium. Obviously, due to the absence of singly ionized helium the $`\mathrm{He}^{++}`$ front cannot propagate beyond the $`\mathrm{H}^+`$ front.

Following these considerations, we impose two restrictions on the evolution of $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$. First, $`Q_{\mathrm{He}^{++}}Q_{\mathrm{H}^+}`$. Second, $`Q_{\mathrm{He}^+}Q_{\mathrm{H}^+}Q_{\mathrm{He}^{++}}`$. Following helium overlap, we impose the same restrictions on the evolution of the mass fraction, i.e. $`F_{\mathrm{He}^{++}}F_{\mathrm{H}^+}`$ and $`F_{\mathrm{He}^+}F_{\mathrm{H}^+}F_{\mathrm{He}^{++}}`$. The full equations governing the evolution of the filling factors and mass fractions are presented in the Appendix. The equations are split into three cases describing three different epochs. The evolution of the filling factors prior to hydrogen overlap are described in case I. Case II describes the evolution of the filling factors and of the ionized hydrogen mass fraction following the overlap of HII regions, but prior to the overlap of HeII regions. Finally, case III describes the evolution of the ionized mass-fractions following the overlap of cosmological HeII regions. Figure 1 provides a schematic representation of the three epochs in reionization history.

## 4 Reionization Histories

In this section we derive reionization histories for different values of $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$. These histories are shown in figures 2 (for case A star-formation) and 3 (for case B star-formation).

For case A we show results for combinations of $`f_{\mathrm{esc}}f_{}=0.0005`$ and 0.005, and $`z_{\mathrm{tran}}=10.0`$ and 20.0. We note the following features. Large values of $`f_{\mathrm{esc}}f_{}`$ allow the massive zero metallicity stars to reionize both hydrogen and helium in the universe at very early times. If the zero metallicity stars do not reionize helium, quasars result in the overlap of HeII regions by $`z4`$, consistent with observations showing transmission just blueward of the helium Ly$`\alpha `$ line at $`z3`$ (Jacobsen et al 1994; Tytler 1995; Davidsen et al. 1996; Hogan et al. 1997; Reimers et al. 1997; Heap et al. 2000; Kriss et al. 2001; Smette et al. 2002) as well as evidence for a temperature rise in the IGM at $`z3.4`$ indicating the reionization of helium at that time (Schaye et al. 2000; Theuns et al. 2002). Interestingly, we also find that it is possible for hydrogen and/or helium to have been reionized twice. This situation arises if the zero-metallicity stars reionize the universe prior to $`z_{\mathrm{tran}}`$, but below $`z_{\mathrm{tran}}`$ the recombination rate is sufficiently high that the softer spectra of normal stars cannot maintain overlap. Overlap of HII regions is then re-achieved at a later time as the collapse fraction grows and the density of the IGM is lowered. A second overlap of HeII regions results from the rise in quasar activity at $`z10`$. Similar results are found in case B, for values of $`f_{\mathrm{esc}}f_{}`$ that are a factor of $`10`$ larger than in case A.

We would like to reiterate the point that multiple reionization epochs occur for some sets of parameters as a result of multiple peaks in the emissivity of ionizing photons as a function of cosmic time. In such cases the first emissivity peak occurs at $`z_{\mathrm{tran}}`$, at which point it drops due to the softer spectra of metal enriched stars. Later, when the collapsed fraction of baryons has reached a sufficiently high value, the resulting increase in emissivity further reionizes the IGM. As we will show in § 7, this double peaked behavior is still seen if a smoothly varying fraction of collapsed objects at metallicities below the threshold produce the hard spectra rather than sources emitting prior to a fixed $`z_{\mathrm{tran}}`$. A double peaked emissivity will always arise so long as the relevant metal enrichment time is less than the time required for the collapse fraction to change by a factor equal to the ratio between the emissivities of the first generation of stars and subsequent metal enriched stars.555Note that even if the mixing of metals is partial and the intergalactic filling factor of metals is smaller than unity (e.g., Scannapieco et al. 2002), new stars are likely to form in regions that were already enriched because galaxies tend to be clustered on large scale sheets and filaments. This effect needs to be included when evaluating the effective $`z_{\mathrm{tran}}`$ from numerical simulations.

There are two small, sharp jumps in the reionization histories that are artifacts, but warrant explanation. The first is in the evolution of $`Q_{\mathrm{He}^+}^\mathrm{m}`$ immediately following the overlap of HeIII regions. This small jump (note that the curve is smooth, not discontinuous) arises because our formalism demands $`Q_{\mathrm{He}^+}^\mathrm{m}=0`$ at this epoch. However following the overlap of He++, the He+ front is free to rapidly expand, powered by the strong ionizing radiation field, and the mass fraction $`F_{\mathrm{He}^+}`$ therefore grows until the He+ front reaches the H+ front. The second jump occurs at $`z=2.5`$ where we switch from computation of the ionizing flux due to quasars using the model quasar luminosity function, to calculation using the empirical luminosity function. Neither is a perfect description around this redshift, and there is a discontinuous jump in the luminosity density, and hence in the ionizing radiation density, resulting in a jump in the reionization histories.

At early times the rise in temperature associated with reionization curtails star-formation. The effect of this is clearly seen in panel (a) of figure 2. Here the zero metallicity stars rapidly increase the ionized fraction, but as $`Q_{\mathrm{H}^+}`$ approaches unity, the ionizing flux is reduced, and the rate of reionization slows. This feature is not seen in cases where overlap occurs at lower redshift (e.g. panel (d) in figure 2). Many combinations of the parameters $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$ predict significant reionization without overlap due to the early presence of metal-free stars, resulting in a peak in the ionized fraction centered on the transition redshift, followed by recombinations before further reionization due to normal stars and quasars. We term these peaks failed overlaps. In section § 6 we show that these failed overlaps leave a significant imprint on the visibility function for electron scattering of CMB photons.

### 4.1 The Optical Depth to Electron Scattering

Different reionization histories result in different densities of electrons as a function of redshift. One simple probe of the reionization history is the optical depth to electron scattering $`\tau _{\mathrm{es}}`$, which depend on the mass filling factors $`Q^\mathrm{m}`$

$$\tau _{\mathrm{es}}=_0^{1000}𝑑z\frac{cdt}{dz}\sigma _T\left[Q_{\mathrm{H}^+}^\mathrm{m}n_\mathrm{H}^0+Q_{\mathrm{He}^+}^\mathrm{m}n_{\mathrm{He}}^0+2Q_{\mathrm{H}^{++}}^\mathrm{m}n_{\mathrm{He}}^0\right](1+z)^3,$$
(16)
where $`\sigma _T=6.652\times 10^{25}\mathrm{cm}^2`$ is the Thomson cross-section. In the following section we plot contours of $`\tau _{\mathrm{es}}`$ as a function of $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$.

## 5 Joint Constraints on $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$

In this section we discuss the constraints that current observations place on $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$. First we discuss the variation of $`\tau _{\mathrm{es}}`$ with $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$. The optical depth for electron scattering is constrained from recent CMB anisotropy experiments to a value $`\tau _{\mathrm{es}}<0.18`$ (Wang & Tegmark 2002). An upper limit on $`\tau _{\mathrm{es}}`$ implies that there is maximum redshift beyond which the universe was neutral. As a second constraint, we note that the spectra of quasars at $`z6`$ do not show a Gunn-Peterson (1965) trough, indicating that HII regions had achieved overlap earlier than that time.

As we have seen in the previous section, the universe is reionized earlier for higher values of $`f_{\mathrm{esc}}f_{}`$, which result in higher values of $`\tau _{\mathrm{es}}`$. Hence $`f_{\mathrm{esc}}f_{}`$ cannot be too large so as not to violate the limit $`\tau _{\mathrm{es}}<0.18`$. On the other hand, if the product $`f_{\mathrm{esc}}f_{}`$ is too small, then the ionizing radiation field will be too weak for the HII regions to overlap by $`z6`$. If the transition redshift $`z_{\mathrm{tran}}`$ is large, then $`f_{\mathrm{esc}}f_{}`$ may obtain higher values since the harder spectra of metal free stars are only available to the ionizing radiation field at high redshifts when the IGM is dense and the collapse fraction is small. However if $`z_{\mathrm{tran}}`$ is low, then $`f_{\mathrm{esc}}f_{}`$ may take smaller values and still result in overlap of HII regions by $`z6`$. In addition to the above constraints, observations suggest that the IGM was metal enriched by $`z6`$ (Songaila 2001) to a level that would lead to a Scalo stellar mass function (Bromm et al. 2001) if new galaxies formed out of it. As a result, $`z_{\mathrm{tran}}`$ is limited to be larger than six.

Figures 4 and 5 show contours of $`\tau _{\mathrm{es}}`$ (solid lines) as a function of $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$ for case A and B star-formation respectively. Also shown are contours of $`Q_{\mathrm{H}^+}`$ at $`z=6`$ (dashed lines), and a vertical line at $`z_{\mathrm{tran}}=6`$ (dot-dashed line). The regions excluded by the aforementioned constraints are shaded grey. For our case A star-formation model we find that $`f_{\mathrm{esc}}f_{}0.004`$ if $`z_{\mathrm{tran}}10`$. For a higher $`z_{\mathrm{tran}}`$, we find that larger values of $`f_{\mathrm{esc}}f_{}`$ are possible. If $`z_{\mathrm{tran}}6`$, $`f_{\mathrm{esc}}f_{}`$ may be as low as $`10^4`$, though for $`z_{\mathrm{tran}}15`$ we find $`f_{\mathrm{esc}}f_{}0.001`$. We find that the lower limits on $`f_{\mathrm{esc}}f_{}`$ are larger by a factor of $`3`$ for our case B star-formation model. However since the intensity of the stellar radiation field is curtailed by feedback at high redshift in this model, we find that $`f_{\mathrm{esc}}f_{}`$ can be as large as $`0.1`$ for $`z_{\mathrm{tran}}15`$, or even larger for $`z_{\mathrm{tran}}15`$. Interestingly, the contours for $`\tau _{\mathrm{es}}`$ and $`z_{\mathrm{reion}}`$ are not parallel. A significant range of $`\tau _{\mathrm{es}}`$ exist for different $`z_{\mathrm{tran}}`$ (particularly where case A star-formation is assumed). This indicates that measurement of optical depth does not determine the overlap redshift. We discuss the reason for this puzzling result in § 6.

In § 4 we mentioned that hydrogen and/or helium could be reionized twice, and showed examples of reionization histories that exhibit this behavior. In figures 6 and 7 we show contours of the overlap redshift for hydrogen (left panels), and helium (right panels), for case A and case B star-formation models. The redshift of the most recent overlap is described by the dark contours. At high values of $`z_{\mathrm{tran}}`$ and low values of $`f_{\mathrm{esc}}f_{}`$, the overlap redshift is almost independent of star-formation since then quasars dominate the ionizing radiation field. This is particularly true for helium, which is ionized by quasars at $`z4`$ in the absence of any stellar ionizing radiation. Cases where there was an additional earlier overlap are denoted by the grey contours. The region of the $`[z_{\mathrm{tran}},f_{\mathrm{esc}}f_{}]`$ plane that results in multiple overlap epochs is shaded grey. We see that this area covers a significant portion of parameter space, and therefore that it is quite possible the universe was reionized twice. The contours representing the upper limit of $`\tau _{\mathrm{es}}<0.18`$, and the lower limit for the overlap redshift $`z>6`$ are also shown (thick dashed lines). As mentioned earlier, parameters in the region between these limits are not excluded by current observations.

## 6 The Visibility Function

The visibility function $`g(z)d\eta `$ describes the probability that an observed photon was scattered between conformal times $`\eta (z)`$ and $`\eta +\mathrm{\Delta }\eta `$

$$g(z)=\frac{d\tau _{\mathrm{es}}}{d\eta }e^{\tau _{\mathrm{es}}}=\sigma _\mathrm{T}\left[Q_{\mathrm{H}^+}^\mathrm{m}n_\mathrm{H}^0+Q_{\mathrm{He}^+}^\mathrm{m}n_{\mathrm{He}}^0+2Q_{\mathrm{He}^{++}}^\mathrm{m}n_{\mathrm{He}}^0\right](1+z)^2e^{\tau _{\mathrm{es}}},$$
(17)
where $`d\eta =(1+z)cdt`$, $`\tau _{\mathrm{es}}`$ is the optical depth from redshift 0 to $`z`$ and the $`Q^\mathrm{m}`$ are mass filling factors. It is commonly assumed that the visibility function is peaked during the recombination era, with a second, lower amplitude, but broader peak coinciding with the epoch of overlap of HII regions. However we have found that overlap may have occurred on multiple occasions, and that many combinations of the parameters $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$ predict significant reionization without overlap (failed overlap) at high redshift due to the early presence of metal-free stars. During these failed overlaps, examples of which are shown in figures 2 and 3, the ionization fraction can be significant. Following a failed overlap at $`z_{\mathrm{tran}}`$, recombinations dominate the evolution and the ionization fraction drops before further reionization ensues due to normal stars and quasars. The higher cosmic density at early times might more than compensate for the smaller ionization fraction, so that the peak of the visibility function does not coincide with the epoch of reionization.

This is indeed what we find. Figures 7 and 8 show visibility functions corresponding to the sample reionization histories shown in figures 2 and 3. Generally, if the early metal-free stars succeed in achieving reionization, then the peak in the visibility function coincides with this overlap epoch (e.g. case A in figures 2 and 3). However the optical depth $`\tau _{\mathrm{es}}`$ is large and many of these cases are already excluded by CMB measurements. Of the 8 cases shown in figures 2 and 3, those that exhibit failed overlap have visibility peaks that coincide with this failed overlap, rather than with the actual reionization epoch. Thus, the visibility function may probe the nature of the early generation of stars rather than the reionization epoch itself.

The reionized gas during a failed overlap will be in bubbles. Gruzinov & Hu (1998) have pointed out that a universe which is ionized in patches induces additional CMB anisotropies on arc-minute scales. Thus, patchy reionization due to a failed overlap at high redshifts will modify the shape of the small-scale power spectrum of CMB anisotropies. However, the amplitude of the added secondary anisotropies was predicted to be at a level well below the feature in the spectrum detected more recently by the Cosmic Background Imager (Padin et al 2001). This feature was explained by Bond et al. (2002) as being due to the foreground Sunyaev-Zeldovich effect from X-ray clusters.

## 7 Justification for Step-Function Approach to the Evolution in the Mode of Star Formation

As mentioned in § 3.2 we have approximated the evolution in the mode of star formation as a step function. That is, we have assumed that the production rate of ionizing photons in the universe is transformed sharply at $`z_{\mathrm{tran}}`$ from being dominated by massive, metal-free stars to being dominated by metal-rich stars with a Scalo IMF. In this section we justify this approach by comparing its results to those obtained with a scheme involving gradual metal enrichment of the IGM as described in § 3.2 by equation (14).

Figure 10 shows two sample reionization histories computed using the gradual enrichment scheme, together with the resulting visibility function. The examples correspond to panels a) and c) for Case B star formation in Figure 3. The thick grey line in these panels shows the fraction of star formation resulting in a Scalo (1998) mass function as a function of redshift. In each of these cases the value of $`C_{\mathrm{metal}}`$ (listed in the figure) produces around 50% Scalo (1998) IMF star formation at the previously used redshift $`z_{\mathrm{tran}}`$. Note the different values of $`C_{\mathrm{metal}}`$, and the different histories of the collapsed mass-fraction. The features discussed for the corresponding examples in Figures 3 and 9 are still present in Figure 10. These include the double overlap of He III in a), and the non-coincidence of the overlap redshift of HII and the peak in the visibility function in b). Furthermore, we find unchanged values of $`\tau _{\mathrm{es}}`$. These calculations appear to justify our use of a single parameter $`z_{\mathrm{tran}}`$ to characterize the reionization histories and visibility functions.

## 8 Conclusion

We have explored the reionization histories of hydrogen and helium due to stars and quasars. The results were analyzed as a function of the two free parameters in our model, namely: (i) the transition redshift, $`z_{\mathrm{tran}}`$, above which the stellar population is dominated by massive, zero metallicity stars; and (ii) the product of the escape fraction of ionizing photons and the star-formation efficiency, $`f_{\mathrm{esc}}f_{}`$. The quasar model was not varied since it provides an excellent fit to all existing data on the luminosity function of quasars up to redshift $`z6`$ (Wyithe & Loeb 2002).

Figures 2 and 3 show sample reionization histories for different choices of $`z_{\mathrm{tran}}`$ and $`f_{\mathrm{esc}}f_{}`$. We find that a wide range of $`z_{\mathrm{tran}}6`$ is allowed for $`f_{\mathrm{esc}}f_{}2\times 10^3`$ if the star-formation efficiency does not depend on galaxy mass (figure 4) or $`f_{\mathrm{esc}}f_{}10^2`$ if the efficiency is suppressed in low mass galaxies (figure 5). This wide range satisfies the constraints that overlap of the HII regions must be achieved by $`z6`$ (Fan et al 2002) and that the optical depth for electron scattering must be limited to $`\tau _{\mathrm{es}}<0.18`$. The allowed range leads generically to $`\tau _{\mathrm{es}}7\%`$ (see figures 4 and 5). The MAP satellite is expected to have sufficient sensitivity to detect values of $`\tau _{\mathrm{es}}`$ as small as $`5\%`$ (Kaplinghat et al. 2002), well below this range of expected values.

A major fraction of the allowed range of $`f_{\mathrm{esc}}f_{}`$ and $`z_{\mathrm{tran}}`$ leads to an early peak in the ionized fraction due to metal-free stars at high redshifts. Often this peak results in a small but significant filling factor that is subsequently reduced (temporarily) due to recombination. In a restricted range of the allowed parameter values, we find that either hydrogen or helium experience two overlap epochs, separated by recombination (see shaded regions in figures 6 and 8). That helium might have been reionized twice due to the presence of early very massive stars was previously pointed out by Oh et al. (2001). The first overlap phase is caused by the population of zero-metallicity, massive stars and the second is dominated by the quasars for helium or by stars and quasars for hydrogen. Even if early overlap is not achieved, the peak in the visibility function for scattering of the CMB often coincides with the early partial ionization peak rather than with the actual reionization epoch (see figures 7 and 9). The resulting value of $`\tau _{\mathrm{es}}`$ is therefore larger than expected based the reionization redshift alone. Thus, the CMB visibility function may be probing the nature of the early generation of stars rather than the reionization epoch itself.

Future CMB experiments, such as MAP and Planck, will provide tighter constraints on $`\tau _{\mathrm{es}}`$ and will reduce the range of allowed $`f_{\mathrm{esc}}f_{}`$ in our model (see figures 4 and 5). Any additional observational data on the composition of the stellar population or the abundance of HeI or HeII at $`z4`$ will provide stronger lower limits on $`z_{\mathrm{tran}}`$.

We thank Rennan Barkana and an anonymous referee for useful detailed comments on the manuscript. JSBW is supported by a Hubble Fellowship grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. AL acknowledges generous support from the Institute for Advanced Study at Princeton and the John Simon Guggenheim Memorial Fellowship. This work was also supported in part by NSF grant AST-0071019.

## Appendix A Equations Describing Reionization of Hydrogen and Helium in a Clumpy Universe

The filling factors $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$, and ionized mass fractions $`F_{\mathrm{H}^+}`$, $`F_{\mathrm{He}^+}`$ and $`F_{\mathrm{He}^{++}}`$ evolve according to six coupled ODEs. It is convenient to write the ODEs and constraints in terms of the following quantities:

$`S_{\mathrm{CF}}^{\mathrm{H}^+}`$ $`=`$ $`{\displaystyle \frac{1}{n_\mathrm{H}^0}}{\displaystyle \frac{dn_\gamma ^{\mathrm{H}^+}}{dz}}`$
$`S_{\mathrm{CF}}^{\mathrm{He}^+}`$ $`=`$ $`{\displaystyle \frac{1}{n_{\mathrm{He}}^0}}{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}`$
$`S_{\mathrm{CF}}^{\mathrm{He}^{++}}`$ $`=`$ $`{\displaystyle \frac{1}{n_{\mathrm{He}}^0}}{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}`$
$`S_\mathrm{C}^{\mathrm{H}^+}`$ $`=`$ $`{\displaystyle \frac{1}{n_\mathrm{H}^0F_{\mathrm{H}^+}}}{\displaystyle \frac{dn_\gamma ^{\mathrm{H}^+}}{dz}}`$
$`S_\mathrm{C}^{\mathrm{He}^+}`$ $`=`$ $`{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}`$
$`S_\mathrm{C}^{\mathrm{He}^{++}}`$ $`=`$ $`{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}}}{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}`$
$`S_{\mathrm{KF}}^{\mathrm{H}^+}`$ $`=`$ $`\alpha _\mathrm{B}^{\mathrm{H}^+}{\displaystyle \frac{R^{\mathrm{H}^+}}{a^3}}n_\mathrm{e}^{\mathrm{H}^+}{\displaystyle \frac{dt}{dz}}`$
$`S_{\mathrm{KF}}^{\mathrm{He}^+}`$ $`=`$ $`\alpha _\mathrm{B}^{\mathrm{He}^+}{\displaystyle \frac{R^{\mathrm{He}^+}}{a^3}}n_\mathrm{e}^{\mathrm{He}^+}{\displaystyle \frac{dt}{dz}}`$
$`S_{\mathrm{KF}}^{\mathrm{He}^{++}}`$ $`=`$ $`\alpha _\mathrm{B}^{\mathrm{He}^{++}}{\displaystyle \frac{R^{\mathrm{He}^{++}}}{a^3}}n_\mathrm{e}^{\mathrm{He}^{++}}{\displaystyle \frac{dt}{dz}}`$
$`S_\mathrm{K}^{\mathrm{H}^+}`$ $`=`$ $`\left(\alpha _\mathrm{B}^{\mathrm{H}^+}{\displaystyle \frac{R^{\mathrm{H}^+}}{a^3}}n_\mathrm{e}^{\mathrm{H}^+}{\displaystyle \frac{dt}{dz}}+{\displaystyle \frac{dF_{\mathrm{H}^+}}{dz}}\right){\displaystyle \frac{Q_{\mathrm{H}^+}}{F_{\mathrm{H}^+}}}`$
$`S_\mathrm{K}^{\mathrm{He}^+}`$ $`=`$ $`\left(\alpha _\mathrm{B}^{\mathrm{He}^+}{\displaystyle \frac{R^{\mathrm{He}^+}}{a^3}}n_\mathrm{e}^{\mathrm{He}^+}{\displaystyle \frac{dt}{dz}}+{\displaystyle \frac{dF_{\mathrm{He}^+}}{dz}}\right){\displaystyle \frac{Q_{\mathrm{He}^+}}{F_{\mathrm{He}^+}}}`$
$`S_\mathrm{K}^{\mathrm{He}^{++}}`$ $`=`$ $`\left(\alpha _\mathrm{B}^{\mathrm{He}^{++}}{\displaystyle \frac{R^{\mathrm{He}^{++}}}{a^3}}n_\mathrm{e}^{\mathrm{He}^{++}}{\displaystyle \frac{dt}{dz}}+{\displaystyle \frac{dF_{\mathrm{He}^{++}}}{dz}}\right){\displaystyle \frac{Q_{\mathrm{He}^{++}}}{F_{\mathrm{He}^{++}}}}.`$ (A1)
In the above definitions, the $`S_\mathrm{C}`$ and $`S_{\mathrm{CF}}`$ are source terms describing the ionizing radiation field, while the $`S_\mathrm{K}`$ and $`S_{\mathrm{KF}}`$ are sink terms describing the recombination rate. In the recombination terms, the $`n_\mathrm{e}`$ are the electron densities available for recombination (equal to the mean cosmological density multiplied by the volume filling factor and the mass fraction). These are therefore

$`n_\mathrm{e}^{\mathrm{H}^+}`$ $`=`$ $`\left[n_\mathrm{H}^0+n_{\mathrm{He}}^0{\displaystyle \frac{Q_{\mathrm{He}^+}F_{\mathrm{He}^+}+2Q_{\mathrm{He}^{++}}F_{\mathrm{He}^{++}}}{Q_{\mathrm{H}^+}F_{\mathrm{H}^+}}}\right]`$
$`n_\mathrm{e}^{\mathrm{He}^+}`$ $`=`$ $`\left[n_\mathrm{H}^0+n_{\mathrm{He}}^0\right]`$
and
$`n_\mathrm{e}^{\mathrm{He}^{++}}`$ $`=`$ $`\left[n_\mathrm{H}^0+2n_{\mathrm{He}}^0\right].`$ (A2)
We write the generalization of equations (4) and (7) for the three species in terms of the above quantities. There are 3 cases corresponding to the period before HII regions overlap, the period following the overlap of HII regions but preceding the overlap of HeII regions, and the period following the overlap of HeII regions:

$``$CASE I) Pre Hydrogen Overlap: If $`Q_{\mathrm{H}^+}<1`$ or $`F_{\mathrm{H}^+}<F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$

Prior to the overlap of HII regions, the mass fractions $`F`$ and their derivatives are set by the critical value of $`\mathrm{\Delta }_{\mathrm{crit}}`$ chosen as the ionization density threshold. The filling factors $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ evolve according to a generalization of equation (4). Hence

$`{\displaystyle \frac{dF_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`{\displaystyle \frac{d}{dz}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})F_{\mathrm{H}^+}=F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$
$`{\displaystyle \frac{dF_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`{\displaystyle \frac{d}{dz}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})F_{\mathrm{He}^+}=F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$
$`{\displaystyle \frac{dF_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`{\displaystyle \frac{d}{dz}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})F_{\mathrm{He}^{++}}=F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$
$`{\displaystyle \frac{dQ_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{H}^+}+{\displaystyle \frac{1}{n_\mathrm{H}^0F_{\mathrm{H}^+}}}\left[\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}+\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}\right]S_\mathrm{K}^{\mathrm{H}^+}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^+}{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_\mathrm{K}^{\mathrm{He}^+}{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^{++}}{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}S_\mathrm{K}^{\mathrm{He}^{++}}.`$ (A3)
In the above, $`\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^+}}{dz}`$ and $`\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{H}^{++}}}{dz}`$ are the excess ionizing photon rates having frequencies $`5.94\times 10^{15}<\nu <1.31\times 10^{16}`$Hz for HeI ionization, and $`\nu >1.31\times 10^{16}`$Hz for HeII ionization. The values of these excess source terms are determined from the constraints on the relative evolution of the different ionization fronts. As these are excess rates, their values must be positive.

If $`Q_{\mathrm{He}^{++}}<Q_{\mathrm{H}^+}`$ then $`Q_{\mathrm{He}^+}`$ is limited to be smaller than $`Q_{\mathrm{H}^+}Q_{\mathrm{He}^{++}}`$ since the $`\mathrm{He}^+`$ front cannot propagate beyond the H+ front. Therefore, if $`Q_{\mathrm{He}^+}=Q_{\mathrm{H}^+}Q_{\mathrm{He}^{++}}`$ we require $`\frac{dQ_{\mathrm{He}^+}}{dz}=\frac{dQ_{\mathrm{H}^+}}{dz}\frac{dQ_{\mathrm{He}^{++}}}{dz}`$, and since $`\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}=0`$ in this case we find

$`\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}`$ $`=`$ $`\left({\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}+{\displaystyle \frac{1}{n_\mathrm{H}^0F_{\mathrm{H}^+}}}\right)^1\left[(S_\mathrm{C}^{\mathrm{He}^+}S_\mathrm{C}^{\mathrm{H}^+})(S_\mathrm{K}^{\mathrm{He}^+}S_\mathrm{K}^{\mathrm{H}^+})\right].`$ (A4)
In addition to the above constraint, the $`\mathrm{He}^{++}`$ filling factor $`Q_{\mathrm{He}^{++}}`$ is limited to be less than $`Q_{\mathrm{H}^+}`$. Therefore, if $`Q_{\mathrm{He}^{++}}=Q_{\mathrm{H}^+}`$ we require $`\frac{dQ_{\mathrm{He}^{++}}}{dz}=\frac{dQ_{\mathrm{H}^+}}{dz}`$ so that the fronts propagate at the same rate, and $`\frac{dQ_{\mathrm{He}^+}}{dz}=0`$ since the He+ front cannot propagate if the HeII region is filling all of the available volume. These constraints yield

$`\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}`$ $`=`$ $`{\displaystyle \frac{n_{\mathrm{He}}^0F_{\mathrm{He}^+}n_\mathrm{H}^0F_{\mathrm{H}^+}}{2n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}+n_\mathrm{H}^0F_{\mathrm{H}^+}}}\left[(S_\mathrm{C}^{\mathrm{He}^{++}}S_\mathrm{C}^{\mathrm{H}^+})(S_\mathrm{K}^{\mathrm{He}^{++}}S_\mathrm{K}^{\mathrm{H}^+})\right]`$
$``$ $`n_{\mathrm{He}}^0F_{\mathrm{He}^+}{\displaystyle \frac{n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}+n_\mathrm{H}^0F_{\mathrm{H}^+}}{2n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}+n_\mathrm{H}^0F_{\mathrm{H}^+}}}\left[(S_\mathrm{C}^{\mathrm{He}^{++}}S_\mathrm{C}^{\mathrm{He}^+})(S_\mathrm{K}^{\mathrm{He}^{++}}S_\mathrm{K}^{\mathrm{He}^+})\right]`$
and
$`\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}\left[(S_\mathrm{C}^{\mathrm{He}^{++}}S_\mathrm{C}^{\mathrm{He}^+})(S_\mathrm{K}^{\mathrm{He}^{++}}S_\mathrm{K}^{\mathrm{He}^+})+{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}\right].`$ (A5)
$``$CASE II) Post Hydrogen Overlap: If $`Q_{\mathrm{H}^+}=1`$ and $`F_{\mathrm{H}^+}>F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$ and ($`Q_{\mathrm{He}^{++}}<1`$ or $`F_{\mathrm{He}^{++}}<F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$).

Following the overlap of HII regions, the mass fraction of hydrogen $`F_{\mathrm{H}^+}`$ is free to evolve according to equation (7) (with modified source term to account for the excess ionizing photons). The filling factors $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ again evolve according to a generalization of equation (4). Note that if $`Q_\mathrm{i}=1`$, substitution of equation (4) into equation (7) yields $`\frac{dQ_\mathrm{i}}{dz}=0`$. Thus, post HII overlap the value of $`Q_{\mathrm{H}^+}`$ remains at unity. However by continuing to follow all the equations we allow for the possibility of a recombination epoch following an early reionization. We have

$`{\displaystyle \frac{dF_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`S_{\mathrm{CF}}^{\mathrm{H}^+}+{\displaystyle \frac{1}{n_\mathrm{H}^0}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_{\mathrm{KF}}^{\mathrm{H}^+}`$
$`{\displaystyle \frac{dF_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`{\displaystyle \frac{d}{dz}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})F_{\mathrm{He}^+}=F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$
$`{\displaystyle \frac{dF_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`{\displaystyle \frac{d}{dz}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})F_{\mathrm{He}^{++}}=F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$
$`{\displaystyle \frac{dQ_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{H}^+}+{\displaystyle \frac{1}{n_\mathrm{H}^0F_{\mathrm{H}^+}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_\mathrm{K}^{\mathrm{H}^+}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^+}{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_\mathrm{K}^{\mathrm{He}^+}{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^{++}}S_\mathrm{K}^{\mathrm{He}^{++}}.`$ (A6)
As before, the values of the excess ionizing photon rates are determined from constraints on the co-evolution of the ionization fronts. In the case where $`Q_{\mathrm{He}^{++}}<Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ is limited to be smaller than $`Q_{\mathrm{He}^+}=Q_{\mathrm{H}^+}Q_{\mathrm{He}^{++}}`$. However as mentioned above, if $`Q_{\mathrm{H}^+}=1`$ then the above equations imply $`\frac{dQ_{\mathrm{H}^+}}{dz}=0`$. We therefore have $`\frac{dQ_{\mathrm{He}^+}}{dz}=\frac{dQ_{\mathrm{He}^{++}}}{dz}`$, and since $`\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}=0`$ in this case we find

$$\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^+}}{dz}=n_{\mathrm{He}}^0F_{\mathrm{He}^+}\left[S_\mathrm{C}^{\mathrm{He}^+}S_\mathrm{K}^{\mathrm{He}^+}\right].$$
(A7)
$``$CASE III) Post Helium Overlap: If $`Q_{\mathrm{He}^{++}}=1`$ and $`F_{\mathrm{He}^{++}}F_\mathrm{M}(\mathrm{\Delta }_{\mathrm{crit}})`$.

Finally, following the overlap of HeII regions, the mass fractions $`F_{\mathrm{H}^+}`$, $`F_{\mathrm{He}^+}`$ and $`F_{\mathrm{He}^{++}}`$ are free to evolve according to equation (7) (with modified source terms to account for the excess ionizing photons). The filling factors $`Q_{\mathrm{H}^+}`$, $`Q_{\mathrm{He}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ still evolve according to generalizations of equation (4). At the point of He++ overlap our formalism demands that there is no He+ in the universe since $`Q_{\mathrm{He}^+}=0`$. However, following He++ overlap the He+ front (which must lead the He++ front) is free to propagate into the denser IGM. Therefore, following the overlap of He++, $`F_{\mathrm{He}^+}`$ becomes the fraction of helium in the universe (equal to zero initially) that is singly ionized. The recombination rate $`R^{\mathrm{He}^+}`$ is calculated

$$R^{\mathrm{He}^+}=_{\mathrm{\Delta }_\mathrm{i}^{\mathrm{He}^{++}}}^{\mathrm{\Delta }_\mathrm{i}^{\mathrm{He}^+}}𝑑\mathrm{\Delta }\frac{dP}{d\mathrm{\Delta }}\mathrm{\Delta }^2,$$
(A8)
where $`F_{\mathrm{He}^{++}}=_0^{\mathrm{\Delta }_\mathrm{i}^{\mathrm{He}^{++}}}𝑑\mathrm{\Delta }\frac{dP}{d\mathrm{\Delta }}\mathrm{\Delta }`$ and $`F_{\mathrm{He}^{++}}+F_{\mathrm{He}^+}=_0^{\mathrm{\Delta }_\mathrm{i}^{\mathrm{He}^+}}𝑑\mathrm{\Delta }\frac{dP}{d\mathrm{\Delta }}\mathrm{\Delta }`$. As before, if $`Q_\mathrm{i}=1`$ substitution of equation (4) into equation (7) yields $`\frac{dQ_\mathrm{i}}{dz}=0`$. The filling factors $`Q_{\mathrm{H}^+}`$ and $`Q_{\mathrm{He}^{++}}`$ therefore remain at unity. We have

$`{\displaystyle \frac{dF_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`S_{\mathrm{CF}}^{\mathrm{H}^+}+{\displaystyle \frac{1}{n_\mathrm{H}^0}}\left[\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}+\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}\right]S_{\mathrm{KF}}^{\mathrm{H}^+}`$
$`{\displaystyle \frac{dF_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`S_{\mathrm{CF}}^{\mathrm{He}^+}{\displaystyle \frac{1}{n_{\mathrm{He}}^0}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_{\mathrm{KF}}^{\mathrm{He}^+}{\displaystyle \frac{dF_{\mathrm{He}^{++}}}{dz}}`$
$`{\displaystyle \frac{dF_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`S_{\mathrm{CF}}^{\mathrm{He}^{++}}{\displaystyle \frac{1}{n_{\mathrm{He}}^0}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}S_{\mathrm{KF}}^{\mathrm{He}^{++}}`$
$`{\displaystyle \frac{dQ_{\mathrm{H}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{H}^+}+{\displaystyle \frac{1}{n_\mathrm{H}^0F_{\mathrm{H}^+}}}\left[\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}+\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}\right]S_\mathrm{K}^{\mathrm{H}^+}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^+}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^+}{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^+}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}S_\mathrm{K}^{\mathrm{He}^+}{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$
$`{\displaystyle \frac{dQ_{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`S_\mathrm{C}^{\mathrm{He}^{++}}{\displaystyle \frac{1}{n_{\mathrm{He}}^0F_{\mathrm{He}^{++}}}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}S_\mathrm{K}^{\mathrm{He}^{++}}.`$ (A9)
The excess ionizing photon rates are determined as follows. If $`F_{\mathrm{He}^{++}}<F_{\mathrm{H}^+}`$ then $`F_{\mathrm{He}^+}`$ is limited to be smaller than $`F_{\mathrm{He}^+}=F_{\mathrm{H}^+}F_{\mathrm{He}^{++}}`$ since the He+ front cannot propagate beyond the H+ front. Therefore, if $`F_{\mathrm{He}^+}F_{\mathrm{H}^+}F_{\mathrm{He}^{++}}`$ we require $`\frac{dF_{\mathrm{He}^+}}{dz}=\frac{dF_{\mathrm{H}^+}}{dz}\frac{dF_{\mathrm{He}^{++}}}{dz}`$, and since $`\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}=0`$ in this case, we find

$$\mathrm{\Delta }\frac{dn_\gamma ^{\mathrm{He}^+}}{dz}=\left(\frac{1}{n_{\mathrm{He}}^0}+\frac{1}{n_\mathrm{H}^0}\right)^1\left[(S_{\mathrm{CF}}^{\mathrm{He}^+}S_{\mathrm{CF}}^{\mathrm{H}^+})(S_{\mathrm{KF}}^{\mathrm{He}^+}S_{\mathrm{KF}}^{\mathrm{H}^+})\right].$$
(A10)
The $`\mathrm{He}^{++}`$ mass fraction $`F_{\mathrm{He}^{++}}`$ is limited to be less than $`F_{\mathrm{H}^+}`$. Therefore, if $`F_{\mathrm{He}^{++}}=F_{\mathrm{H}^+}`$ we require $`\frac{dF_{\mathrm{He}^{++}}}{dz}=\frac{dF_{\mathrm{H}^+}}{dz}`$ so that the fronts propagate at the same rate, and $`\frac{dF_{\mathrm{He}^+}}{dz}=0`$ since the He+ front cannot grow if the $`\mathrm{He}^{++}`$ is filling all of the available volume. These constraints yield

$`\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}`$ $`=`$ $`{\displaystyle \frac{n_{\mathrm{He}}^0n_\mathrm{H}^0}{2n_{\mathrm{He}}^0+n_\mathrm{H}^0}}\left[(S_{\mathrm{CF}}^{\mathrm{He}^{++}}S_{\mathrm{CF}}^{\mathrm{H}^+})(S_{\mathrm{KF}}^{\mathrm{He}^{++}}S_{\mathrm{KF}}^{\mathrm{H}^+})\right]`$
$``$ $`n_{\mathrm{He}}^0{\displaystyle \frac{n_{\mathrm{He}}^0+n_\mathrm{H}^0}{2n_{\mathrm{He}}^0+n_\mathrm{H}^0}}\left[(S_{\mathrm{CF}}^{\mathrm{He}^{++}}S_{\mathrm{CF}}^{\mathrm{He}^+})(S_{\mathrm{KF}}^{\mathrm{He}^{++}}S_{\mathrm{KF}}^{\mathrm{He}^+})\right]`$
and
$`\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^{++}}}{dz}}`$ $`=`$ $`n_{\mathrm{He}}^0\left[(S_{\mathrm{CF}}^{\mathrm{He}^{++}}S_{\mathrm{CF}}^{\mathrm{He}^+})(S_{\mathrm{KF}}^{\mathrm{He}^{++}}S_{\mathrm{KF}}^{\mathrm{He}^+})+{\displaystyle \frac{1}{n_{\mathrm{He}}^0}}\mathrm{\Delta }{\displaystyle \frac{dn_\gamma ^{\mathrm{He}^+}}{dz}}\right].`$ (A11)


format: text

text:
# Deconfining transition in Full QCD

## 1 Introduction

A schematic phase diagram of full QCD with 2 dynamical flavours ($`m_u=m_d=m_q`$) is shown in Fig. 1. The upper part of the diagram is the deconfined phase, the lower part is confined. The line is determined by the maxima of a number of susceptibilities , including the susceptibility $`\chi _L`$ of the Polyakov line

$$\chi _L=d^3xL(\stackrel{}{x},0)L^{}(\stackrel{}{0},0)$$
(1)
and the susceptibility $`\chi _{\mathrm{ch}}`$ of the chiral order parameter

$$\chi _{\mathrm{ch}}=d^3x\overline{\psi }\psi (\stackrel{}{x},0)\overline{\psi }\psi (\stackrel{}{0},0).$$
(2)
All of them have a maximum at the same value of $`T`$, for a given $`m_q`$, which defines the line in Fig. 1. For $`m_q>`$ 3 GeV, the maxima of $`\chi _L`$ diverge proportionally to the volume $`V`$, indicating a first order transition. At $`m_q=0`$ there are theoretical reasons and numerical indications that the transition is second order . At intermediate values of $`m_q`$ (tiny part of the line in Fig. 1) the susceptibilities do not diverge with $`V`$, and this is interpreted as absence of a phase transition: the line would correspond to a crossover.

Across the transition the density of the free energy $`F`$ is a function of the order parameters. The singularities of derivatives of $`F`$ are related to susceptibilities of the order parameters. In our case $`L`$ is the order parameter only at $`m_q=\mathrm{}`$, $`\overline{\psi }\psi `$ only at $`m_q=0`$. A good order parameter in the whole range of $`m_q`$’s could lead to a different assignment for the order of the transition.

Such a parameter could be the disorder parameter $`\mu `$ which describes condensation of magnetic charges . That parameter has been constructed and tested in quenched theory  and its definition can be extended to full QCD . In the spirit of $`N_c\mathrm{}`$ arguments one expects that the mechanism of confinement is the same in quenched and full QCD, quark loops being non leading in $`1/N_c`$ expansion. The operator $`\mu `$ creates a magnetic charge, as defined by some given abelian projection. $`\mu 0`$ signals dual superconductivity. In the quenched theory the specific choice of the abelian projection proves to be immaterial . For the details about the definition of $`\mu `$ we refer to . For the quenched theory it is found that $`\mu 0`$ for $`T<T_c`$ and that $`\mu =0`$ for $`T>T_c`$, $`\mu \tau ^\delta `$ as $`TT_c^{}`$, with $`\tau =1T/T_c`$. A finite size scaling analysis of the infinite volume limit yields $`\delta =0.20(3)`$ for $`SU(2)`$; $`\delta =0.50(2)`$ for $`SU(3)`$. The critical index $`\nu `$ has the values $`\nu =0.62(1)`$ for $`SU(2)`$, in agreement with Ref. , $`\nu =0.33(1)`$ for $`SU(3)`$, corresponding to a first order transition . We repeat the same analysis for 2 staggered flavours, using the Wilson action for the pure gauge sector, and $`12^3\times 4`$, $`16^3\times 4`$,$`32^3\times 4`$ lattices. The machines used are APEmille crates. Part of the results are already published in Ref. .

## 2 Numerical Results

As usual instead of $`\mu `$ itself we determine the quantity $`\rho =\frac{}{\beta }\mathrm{log}\mu `$ in terms of which $`\mu =\mathrm{exp}\left(_0^\beta \rho (\beta ^{})𝑑\beta ^{}\right)`$. For $`T<T_c`$ we find that $`\rho `$ is practically size independent, i.e. that $`\mu 0`$ as $`N_s\mathrm{}`$ (see Fig. 2). For $`T>T_c`$ we find (Fig. 3)

$$\rho =kN_s+\mathrm{const}.$$
(3)
i.e. that $`\mu `$ is strictly zero as $`N_s\mathrm{}`$. More extended checks of superselection of magnetic charge in the deconfined phase have been done, showing that for different magnetic charges created by $`\mu `$, $`\rho \mathrm{}`$ if the net magnetic charge is non-zero, $`\rho \mathrm{const}.`$ if it is zero. Some examples are given in Fig. 3.

Around $`T_c`$ a finite size scaling analysis goes as follows. By dimensional arguments one can parametrize $`\mu `$ as

$$\mu =\tau ^\delta \varphi (\frac{a}{\xi },\frac{N_s}{\xi },m_qN_s^\gamma )$$
(4)
where $`a`$ is the lattice spacing and $`\tau =1T/T_c`$. If the correlation length goes large, $`\xi \tau ^\nu `$, $`a/\xi 1`$ and the dependence on $`a`$ can be neglected. The variable $`N_s/\xi `$ can be traded for $`\tau N_s^{1/\nu }`$. In the quenched theory $`m_q`$ is absent and $`\mu =\tau ^\delta \overline{\varphi }(0,\tau N_s^{1/\nu })`$, which gives for $`\rho `$ the scaling behaviour $`\rho /N_s^{1/\nu }=f(\tau N_s^{1/\nu })`$, that can be tested and allows the determination of $`T_c`$, $`\nu `$ and $`\delta `$. In particular it implies that the height of the peak scales as $`N_s^{1/\nu }`$.

In presence of dynamical quarks, $`m_q`$ sets another scale, and we face a two scale problem. Simulations done at fixed ($`m_\pi /m_\rho `$ show that the height of the peak $`\rho `$ roughly scales as $`N_s^3=V`$ (see Fig. 4). A more refined analysis can be done by choosing $`m_q`$ and $`N_s`$ such that $`m_qN_s^\gamma `$ in Eq. (4) is kept constant, and the problem is again reduced to a one scale problem, allowing to determine $`\nu `$, which gives information on the order of the transition. The index $`\gamma `$ is known to be $`\gamma 2.49`$ , so that the three sets ($`m_q=0.075,N_s=16`$), ($`m_q=0.043,N_s=20`$) and ($`m_q=0.01335,N_s=32`$) keep $`m_qN_s^\gamma `$ constant. Preliminary results indicate that $`\nu =1/3`$, or that the transition is compatible with first order. This is clearly visible from Fig. 5, where the height of the peak is plotted versus $`N_s^3`$. More detailed numerical analyses involving more values of $`N_s`$, $`m_q`$ and the use of an improved action, to improve the quality of scaling, are needed to draw a definite conclusion.

In any case it is demonstrated that also in full QCD confinement is produced by monopole condensation.


text:
# BL Lac Identification for the Ultraluminous X-ray Source Observed in the Direction of NGC 4698

## 1 Introduction

In recent years, the improved imaging capabilities and increased sensitivity of *ROSAT*, *Chandra* and *XMM-Newton* have allowed us to effectively study discrete sources in nearby galaxies beyond the Local Group. Particularly intriguing is the discovery of off-nuclear X-ray sources with luminosities well above the Eddington limit for a typical neutron star, $`10^{38}`$ erg s-1, and up to $`2\times 10^{40}`$ erg s-1 (e.g., Read et al. 1997; Colbert & Mushotzky 1999; Roberts & Warwick 2000; Makishima et al. 2000; Fabbiano et al. 2001). These sources are typically called ultraluminous X-ray sources (ULXs). Despite much effort, little is presently known about these sources. The identification of their optical counterparts is often problematical when the sources are superposed against regions of high surface brightness in the host galaxy. To date, only two ULXs appear to have a clear optical identification (Roberts et al. 2001; Wu et al. 2002), while for others it has been possible to study only the nearby environment (Pakull & Mirioni 2002; Wang 2002).

We have started a search for ULXs in a sample of nearby galaxies (Foschini et al. 2002). This paper concerns follow-up observations of the ULX in NGC 4698. We discuss the radio and optical counterparts and give a redshift estimate. We identify the source with a background source, most likely a BL Lac object at $`z=0.43`$.

## 2 X-ray data (*XMM-Newton*)

NGC 4698 is an Sab spiral galaxy located in the Virgo cluster ($`d=16.8`$ Mpc). It hosts an active nucleus, classified as a Seyfert 1.9 by Ho et al. (1997). The galaxy was observed on 16 December 2001 using the European Photon Imaging Camera (EPIC) on board the *XMM-Newton* satellite. EPIC is composed of two instruments: the PN-CCD camera (Strüder et al. 2001) and two MOS-CCD detectors (Turner et al. 2001). The effective exposure time was $`9.2`$ ks.

We find only one ULX apparently associated with NGC 4698 (Foschini et al. 2002), located $`73\mathrm{}`$ from the optical nucleus of the galaxy, at $`\alpha =12^h48^m25.9^s`$ and $`\delta =+08^{}30\mathrm{}20\mathrm{}`$ (J2000). The point source centroid is measured to precision of better than $`0\stackrel{}{.}1`$, and the absolute pointing uncertainty is $`<4\mathrm{}`$ (Jansen et al. 2001). We detected $`198`$ photons with MOS1, in a circle of radius $`30\mathrm{}`$, $`249`$ with MOS2, and $`414`$ with PN. These statistics are sufficient to perform spectral fitting with simple models. The ones we use are the power law (PL), black body (BB), bremsstrahlung (BR), and the multicolor black body disk (MCD, *diskbb* in Xspec) by Mitsuda et al. (1984). Fluxes were corrected according to the energy encircled fraction (Ghizzardi 2001).

For the processing, screening, and analysis of the data we used the standard tools in the XMM-SAS software (v. 5.2) and HEAsoft Xspec (v. 11.0.1). Correction for vignetting has not been applied because the source is close to the center of the field of view ($`<2\mathrm{}`$) and most of the detected photons have energies less than $`5`$ keV (see Lumb 2002).

The best fit (Fig. 3) is found using the power-law model with $`\mathrm{\Gamma }=2.0\pm 0.2`$ ($`\chi ^2=18.3`$, $`\nu =24`$), giving a flux of $`8.6\times 10^{14}`$ erg cm-2 s-1. If the source is located in NGC 4698, for which we assume a distance of 16.8 Mpc, the corresponding luminosity, assuming isotropic emission, is $`3\times 10^{39}`$ erg s-1. For $`z=0.43`$, as discussed below, the luminosity becomes $`3.9\times 10^{43}`$ erg s-1. Additional information on the fits using the other models is given in Table 1.

## 3 Radio data (VLA)

The VLA observations of NGC 4698 were performed by Ho & Ulvestad (2001) on 29 August 1999 (20 cm) and 31 October 1999 (6 cm). The source detection at 6 cm was reported by Ho & Ulvestad (2001); the flux density is $`1.13`$ mJy and the background noise is $`0.072`$ mJy beam-1. The source, which appears largely unresolved at a resolution of $`1\mathrm{}`$, is located at $`\alpha =12^h48^m25.9^s`$ and $`\delta =+08^{}30\mathrm{}21\mathrm{}`$ (J2000); the radio position is accurate to $`0\stackrel{}{.}1`$.

A reanalysis of the VLA data detected the source also at 20 cm. It is unresolved with the 1$`\mathrm{}`$ beam, and its 20 cm position is consistent with that measured at 6 cm. We find a flux density of $`0.34`$ mJy, with a background noise of $`0.032`$ mJy beam-1. The spectral index, defined as $`S_\nu \nu ^\alpha `$, is $`\alpha _r1.0`$.

## 4 Optical data

### 4.1 Digitized Sky Survey (DSS)

The optical counterpart of the source is visible on the DSS and the HST observations of NGC 4698. We are certain of the identification from the relative offset between the nucleus and the optical counterpart. From the Automatic Plate Measurement (APM111http://www.ast.cam.ac.uk/$``$apmcat/.) catalog, the point source has $`R=18.9`$ mag and $`B=21.9`$ mag, indicating a very red color ($`BR=3.0`$ mag). The original observation was made with the 48 inch Schmidt telescope at Palomar Observatory on 17 February 1950.

### 4.2 Hubble Space Telescope (HST)

We retrieved the calibrated F450W, F606W, and F814W WFPC2 images of NGC 4698 from the HST Archive222Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Institute. STScI is operated by the association of Universities for Research in Astronomy, Inc. under the NASA contract NAS 5-26555. to perform photometry of the optical counterpart of the source. The observations with the F450W and F814W filters were performed on 10 August 2001, while the observation with F606W was taken on 31 March 2000. In all the filters, the source is located on the WF detector, which has a scale of 0$`\stackrel{}{.}`$1 per pixel. The total exposure time is 460 s ($`2\times 230`$ s) in the F450W and F814W filters and 600 s (400 s and 200 s) in the F606W filter. We combined the separate images taken in each filter using the crrej task in IRAF to reject cosmic rays. Aperture photometry was done using the tasks available in the apphot package. The integrated magnitudes, determined using the aperture growth curve method, were transformed to the Vega magnitude system using the zeropoints given in the HST Data Handbook. We find $`m_{\mathrm{F450W}}=21.1\pm 0.1`$ mag, $`m_{\mathrm{F606W}}=19.9\pm 0.1`$ mag, and $`m_{\mathrm{F814W}}=18.6\pm 0.1`$ mag. In the Johnson-Kron-Cousins system these values translate into $`B=21.3\pm 0.1`$ mag, $`R=19.6\pm 0.1`$ mag, and $`I=18.5\pm 0.1`$ mag.

### 4.3 Las Campanas (Magellan)

XMMU J$`124825.9+083020`$ was observed for 900 s ($`3\times 300`$ s) on 6 April 2002 with the LDSS-2 spectrograph (Allington–Smith et al. 1994) on the 6.5 m Baade telescope at Las Campanas Observatory. We observed through a $`1\stackrel{}{.}25`$ longslit in $`0\stackrel{}{.}9`$ seeing, with the slit oriented at the parallactic angle. We used the medium red grism, which gives a dispersion of $`5.3`$ Å pixel-1. The FWHM spectral resolution was $`16`$ Å, and the spatial resolution was $`0\stackrel{}{.}38`$ pixel-1.

Basic data reduction was performed using the IRAF package. The individual spectroscopic frames were corrected for overscan, flat-fielded using domeflats and summed to obtain the final object frame. The spectral extraction was done by summing the counts within an aperture of $`6`$ pixels ($`2\stackrel{}{.}3`$). Wavelength calibration was achieved using a polynomial fit to the lines in the sky spectrum.

The spectrum in the range $`50007000`$ Å is approximately described by the form $`f_\lambda \lambda ^4`$.

### 4.4 Very Large Telescope (VLT)

VLT observations333Based on observations made with ESO Telescopes at the Paranal Observatory under programme ID $`269.D5014`$. were performed with VLT-UT3 (Melipal) plus the FORS1 spectrograph. FORS1 is equipped with a $`2048\times 2048`$ pixel Tektronix CCD, which covers a $`6\stackrel{}{.}8\times 6\stackrel{}{.}8`$ field in the standard resolution imaging mode with a scale of $`0\stackrel{}{.}2`$ pixel-1.

Two spectra, both with an exposure time of 570 s, were acquired on 21 April 2002, starting at 4:52:50 UT. The spectra were acquired using Grism #150I plus order separator GG435, which avoids overlapping of spectral orders over a given wavelength; this limited the spectral range to $`45009000`$ Å. The slit width was $`1\mathrm{}`$ for both spectra, and this setup secured a final dispersion of $`5.5`$ Å pixel-1, corresponding to a FWHM resolution of 13 Å. Before the spectroscopic observation, a one minute $`R`$-band acquisition image was obtained on the same night starting at 04:42:52 UT under very good seeing conditions ($``$0$`\stackrel{}{.}`$6). The object appeared clearly elliptical and extended, with a bright core and a fuzzy halo. Moreover, an irregular spot located around 1$`\mathrm{}`$ north of the object core is apparent in the VLT image.

After correction for flat-field and bias, the spectra were background subtracted and optimally extracted (Horne 1986) using IRAF. He-Ne-Ar and Hg-Cd lamps were used for wavelength calibration. The wavelength calibration was checked against the position of night sky lines; the typical error was 0.5 Å. Finally, the two spectra were stacked together in order to increase the signal-to-noise ratio. We encountered problems during observation of the spectrophotometric standard, and hence our spectra are not flux calibrated.

Fig. 4 shows the normalized VLT spectrum, with several absorption lines identified \[see Laurent-Muehleisen et al. (1998) for a discussion on optical identification of BL Lacs\]. The detected lines include Ca II H&K $`\lambda \lambda 3933,3968`$, the G band $`\lambda 4304`$, H$`\beta `$ $`\lambda 4861`$, and Fe I $`\lambda 5270`$: they all indicate $`z0.43`$. The line strengths are heavily diluted by the featureless continuum, but the spectrum shares close similarity to that of the BL Lac object $`0548322`$ (Barth et al. 2002). We shifted the spectrum of $`0548322`$, which has a redshift of 0.069, to $`z=0.43`$.

## 5 Discussion

### 5.1 Classification

The absorption features detected in the VLT spectrum clearly places XMMU J124825.9+083020 at $`z=0.43`$. We identify the source as a BL Lac object, for the following reasons. (1) We detect no emission lines; the upper limit to any emission feature is $``$5 Å. (2) The stellar features superposed on the featureless continuum have strengths generally consistent with those of other BL Lac objects. We demonstrate this concretely by comparing a redshifted spectrum of the BL Lac object $`0548322`$ (Fig. 4). The break contrast at 4000 Å is $``$29%, again similar to other BL Lac objects (Laurent-Muehleisen et al. 1998). (3) The optical continuum slope measured in the Magellan spectrum is consistent with those of high-frequency peaked BL Lac objects (e.g., Stickel et al. 1993). And (4) the multiwavelength spectral energy distribution (Fig. 5), although not assembled from simultaneous observations, is highly reminiscent of those of BL Lac objects.

By using two-point spectral indices, namely the radio-to-optical $`\alpha _{\mathrm{ro}}`$ and optical-to-X-ray $`\alpha _{\mathrm{ox}}`$, it is possible to show that different objects populate different regions of the $`\alpha _{\mathrm{ro}}\alpha _{\mathrm{ox}}`$ plane (e.g., Brinkmann et al. 1997; Laurent-Muehleisen et al. 1999). Another method has been suggested by Maccacaro et al. (1988), who proposed a nomograph to link the X-ray flux in the energy band $`0.33.5`$ keV and the visual magnitude. The values of $`\alpha _{\mathrm{ro}}`$ and $`\alpha _{\mathrm{ox}}`$ for the present source are $`0.42`$ and $`0.95`$, respectively, thus placing it in the region of X-ray selected BL Lacs (Brinkmann et al. 1997), or high-energy peaked BL Lacs in the $`\alpha _{\mathrm{ro}}\alpha _{\mathrm{ox}}`$ diagram of Laurent-Muehleisen et al. (1999). The nomograph of Maccacaro et al. (1988) gives a ratio $`f_\mathrm{x}/f_\mathrm{v}`$ between $`1.3`$ and $`3.8`$ (depending on whether we use the HST optical magnitudes through the F606W or F450W filter, respectively), in the regime of AGNs and BL Lacs.

The spectral indices can be used to deduce some general properties of the dominant radiation mechanism. If $`\alpha _\mathrm{x}\alpha _\mathrm{r}`$, the source may exhibit relativistic beaming, while if $`\alpha _\mathrm{x}>\alpha _\mathrm{r}`$, it may fulfill the conditions of the homogeneous synchrotron model (Harris & Krawczynski 2002). In our case, we have $`\alpha _\mathrm{x}\alpha _\mathrm{r}1`$, so that we cannot clearly discriminate between these two cases. If the source is a high-frequency peaked BL Lac, however, it is likely that the homogeneous synchrotron model is more applicable. This is confirmed by the $`\alpha _{\mathrm{xox}}`$ test of Sambruna et al. (1996): in this case, the difference by $`\alpha _{\mathrm{ox}}\alpha _\mathrm{x}=\alpha _{\mathrm{xox}}`$ is approximately equal to zero, so avoiding a clear discrimination between the physical mechanism of the source.

The spectral energy distribution (Fig. 5) is comparable to those of BL Lac objects peaked in the $`\gamma `$-ray domain (e.g., Fossati et al. 1998), but the third EGRET catalog (Hartman et al. 1999) does not have any source within several degrees of XMMU J$`124825.9+083020`$.

### 5.2 Search for ULXs and Contamination with Background Objects

This research is a useful demonstration of just how difficult it is to identify the physical nature of ULXs. It shows the vital importance of redshift determinations. Although we had a lot of photometric data on XMMU J$`124825.9+083020`$, by themselves they were insufficient to clearly establish whether the source belongs to NGC 4698 or is a background object. In the first case, the source would have been something similar to a microquasar (e.g., Mirabel & Rodríguez 1999), probably located in a globular cluster. Indeed, the optical image of the source appeared to be slightly extended (angular extent $`3\mathrm{}`$), and its position with respect to NGC 4698 suggested that it could be a globular cluster, albeit an unusually large one. The possibility of an accreting black hole in a globular cluster is not so remote: *Chandra* observations of the globular cluster system of NGC 4472 show that about $`40`$% of the bright low mass X-ray binaries are associated with optically identified globular clusters (Kundu et al. 2002). In addition, the X-ray luminosity function shows a break near $`3\times 10^{38}`$ erg s-1, suggesting that the brightest X-ray binaries are accreting black holes (Kundu et al. 2002), perhaps microquasars.

*XMM-Newton* observations of the Lockman Hole (Hasinger et al. 2001) show that the number of background sources in the energy band $`0.52`$ keV with flux greater than $`4.0\times 10^{14}`$ erg cm-2 s-1 (the best fit value from Tab. 1) is $`15`$ deg-2. In the energy band $`210`$ keV, there are 40 sources deg-2 with flux higher than $`4.6\times 10^{14}`$ erg cm-2 s-1. Assuming the same $`\mathrm{log}N\mathrm{log}S`$ relation, and considering that the $`D_{25}`$ area of NGC 4698 is about $`7.9`$ arcmin2, we expect $`0.08`$ background objects in the $`210`$ keV energy band and $`0.03`$ in the energy band $`0.52`$ keV. However, despite these low values, we have found, in the present case, that the only ULX is a background AGN.

The above calculations could be underestimated in the present case because the VLT images show an unknown concentration of galaxies north-east of NGC 4698 (see Fig. 1), thus suggesting the possibility of a statistically meaningful excess of background sources. However, no additional X-ray sources is seen in the present *XMM-Newton* observation. Perhaps, a longer exposure may reveal soft X-ray emission or additional sources, if NGC 4698 lies along the line of sight to a galaxy cluster. The three X-ray sources identified to date, however, have three redshifts: XMMU J$`124825.9+083020`$ has $`z=0.43`$, NGC 4698 has $`z=0.0033`$, and the ROSAT source 1RXS J$`124828.1+083103`$ has been recently identified with a Seyfert nucleus at $`z=0.12`$ (Xu et al. 2001). Therefore, these three sources are not members of a single cluster.

## 6 Final Remarks

We have presented a multiwavelength analysis of the ultraluminous X-ray source XMMU J124825.9+083020 observed in the direction of NGC 4698. The availability of X-ray, optical, and radio data makes it possible to construct the spectral energy distribution, but the crucial information to decipher the true nature of the source was provided by the optical spectrum. We show that the absorption features are unambiguously associated with a background source at $`z=0.43`$. We argue that the source is most likely a BL Lac object, most likely of the high-frequency peaked category, based on the absence of emission lines, the dilution of the absorption features, the slope of the optical continuum, and the overall shape of the spectral energy distribution.

###### Acknowledgements.

We would like to thank A. Pizzella for useful discussion about NGC 4698, the Service Mode Observations personnel at ESO, and in particular M. Romaniello, for the help in the preparation of VLT observations. This work is based on observations obtained with *XMM-Newton*, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). This research has made use of the NASA Astrophysics Data System Abstract Service and of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We acknowledge the partial support of the Italian Space Agency (ASI) to this research. L. C. H. is grateful for financial support from the Carnegie Institution of Washington and from NASA grants.