text:
# Contamination of cluster radio sources in the measurement of the thermal Sunyaev-Zel’dovich angular power spectrum

## 1 Introduction

As a result of gravitationally driven shocks and compression, majority of the baryons in the most massive, collapsed objects in the universe such as groups and clusters of galaxies exist in the form of hot plasma with temperature of $`10^710^8`$ K. These hot baryons manifest themselves by strong diffuse X-ray sources in the sky within the framework of bremsstrahlung emission, which are directly detectable with current X-ray instruments. On the other hand, the hot electrons in the virialized massive dark halos also scatter the passing cosmic microwave background (CMB) photons, giving rise to a subtle change in the CMB spectrum, known as the Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zel’dovich 1972). Unlike the X-ray emission of the hot gas in massive halos, which shows a strong dependence on the detailed structures and enrichness history of the hot gas as well as the distances of the host halos from us, the SZ effect is completely determined by the total thermal energy of the gas intrinsic to the systems. Therefore, it is well suited for studies of distant, massive clusters, which may provide a stringent constraint on cosmological parameters especially for the mean matter density $`\mathrm{\Omega }_\mathrm{M}`$ and the normalization of matter power spectrum $`\sigma _8`$ (e.g. Barbosa et al. 1996; Molnar, Birkinshaw & Mushotzky 2002; Komatsu & Seljak 2002). A more practical and powerful utilization of the SZ effect nowadays is perhaps to statistically measure the angular power spectrum of the SZ signals at various scales. This allows us to gain the weak SZ signals at high statistical significance levels without resolving individual groups and clusters. Indeed, recent detection of the excess power relative to primordial CMB anisotropy at arcminute scales has been attributed to the statistical signals of the thermal SZ effect (Dawson et al. 2001; Bond et al. 2003). A number of ongoing and future experiments will soon be able to provide a more robust measurement of the SZ angular power spectrum at high multipoles or small angular scales.

There are many other effects that should be corrected for before the SZ angular power spectrum is used for cosmological purposes. These include nongravitational effect on the distribution and evolution of the hot gas inside halos (da Silva et al. 2000, 2001; Seljak, Burwell & Pen 2001; Xue & Wu 2001; Holder & Carlstrom 2001; Zhang, Pen & Wang, 2002; Komatsu & Seljak 2002; Zhang & Wu 2003; etc.), possible asymmetric shape of dark halos (Cooray 2000; Lee & Suto 2003; Piffaretti, Jetzer & Schindler 2003), contaminations of radio point sources (Cooray et al. 1998; Lin, Chiueh & Wu 2002; Holder 2002; Rubinõ-Martín & Sunyaev 2002), etc. Despite large uncertainties, many attempts have been made over the past few years towards the evaluation of these effects on the SZ cluster counts and angular power spectrum measurements. In this work, we will explore the contamination of radio sources in clusters in the measurement of the thermal SZ angular power spectrum, namely, cluster radio halos and galaxies.

In theory of hierarchical clustering, massive dark halos like clusters of galaxies form by gravitational aggregation of individual low-mass objects. Sub-halo mergers will generate a considerably large number of energetic electrons that can be accelerated by cluster magnetic fields over cluster regions, giving rise to the radio halos due to synchrotron emission (Buote 2001 and reference therein). Regardless of the debate on the acceleration mechanism for relativistic electrons in clusters, about $`1/3`$ of nearby, rich clusters exhibit radio halos, and this fraction may increase with redshift if mergers were more frequent in the past (Enßlin & Röttgering 2002). Furthermore, there is a strong observational evidence that cluster radio halos show a morphology similar to the diffuse X-ray emission from the hot intracluster gas, and the corresponding radio power is correlated strongly with the X-ray luminosity and temperature of clusters (Liang et al. 2000; Govoni et al. 2001a,b). This enables one to estimate the distribution and evolution of cluster radio halos based on the well-determined X-ray properties of clusters. For example, one can easily derive the radio luminosity function of clusters at different wavebands and redshifts (Enßlin & Röttgering 2002). Together with the radio surface brightness profiles of clusters motivated by empirical formula (e.g. $`\beta `$ model) for the X-ray surface brightness distribution of clusters, one can estimate the radio background from the diffuse radio halos of clusters.

Existence of diffuse radio halos constitutes a major source of contaminations in the measurement of the thermal SZ effect toward clusters. Since radio halos always contribute a positive signal superimposed on the CMB sky, the decrements of the CMB signals due to the thermal SZ effect of the hot intracluster gas at low frequencies would be compensated, while the increments of the SZ signals at high frequencies will be correspondingly enhanced. The crossover point occurs at $`\nu =217`$ GHz. Moreover, the amplitude of such a modification to the thermal SZ signals depends critically on observing frequencies. This arises because the radio flux of cluster halos as a result of synchrotron emission has a steep spectrum of roughly $`S_\nu \nu ^1`$. Consequently, at very high frequencies the cluster radio halos may become hardly visible, and no modification to the SZ signals from the cluster radio halos is needed. However, with the decrease of observing frequencies the SZ decrements and the diffuse radio emission of cluster halos may appear to be comparable, and eventually the latter would become the dominant component in the microwave sky. This will lead to a nonnegligible correction to the SZ measurement (Enßlin 2002) especially in the central regions of clusters, because the radio emission is more centrally concentrated than the thermal SZ effect if the radio emission is assumed to follow the X-ray emission. Therefore, it deserves a close investigation of the influence of the cluster radio halos on the measurement of the SZ angular power spectrum.

Another source of contaminations in the SZ measurements is radio point sources, which include radio galaxies inside clusters and tracing large-scale structures of the universe. In a similar way to cluster radio halos, cluster radio galaxies also compensate the decrements of the thermal SZ intensity in the Rayleigh-Jeans region. On the other hand, the Poisson distribution and clustering of radio galaxies on large scales would enhance the radio fluctuations, resulting in a positive contribution to both primary CMB and SZ power spectra. The role of these field radio galaxies has been well addressed by several authors (e.g. Franceschini et al. 1989; Toffolatti et al. 1998), and the resulting contamination in the measurement of CMB anisotropies is found to be generally small or can be somehow removed. Nevertheless, a significant fraction of the SZ signals could be diluted by the unsubtracted radio sources inside clusters (Cooray et al. 1998; Lin, Chiueh & Wu 2002). Using a toy model of radio point source subtraction, in which a certain central region is removed from all clusters, Holder (2002) showed that the true SZ power at $`l>1000`$ can be underestimated by as much as $`30\%`$ at $`\nu =30`$ GHz. A more realistic model for the radial distribution of radio galaxies in clusters is thus needed to refine this estimation.

In this paper, using a purely analytical approach, we perform a quantitative study of the contamination of cluster radio halos and galaxies in the measurement of the angular power spectrum of the thermal SZ effect. It is hoped that such a study may clarify the issues of to what extent one has to take the correction of cluster radio sources into account in the cosmological applications of the thermal SZ angular power spectrum. Throughout the paper we will work with a flat cosmological model of $`\mathrm{\Omega }_\mathrm{M}=0.35`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$.

## 2 Thermal SZ angular power spectrum

The change in CMB temperature $`T_{\mathrm{CMB}}`$ in direction $`𝜽`$ due to the thermal SZ effect of the hot electrons in massive clusters is simply

$`{\displaystyle \frac{\mathrm{\Delta }T(𝜽)}{T_{\mathrm{CMB}}}}`$ $`=`$ $`g(x)y\left(𝜽\right);`$ (1)
$`y(\theta )`$ $`=`$ $`{\displaystyle n_\mathrm{e}\sigma _\mathrm{T}\left(\frac{kT_\mathrm{e}}{m_\mathrm{e}c^2}\right)𝑑\chi };`$ (2)
$`g(x)`$ $`=`$ $`x{\displaystyle \frac{e^x+1}{e^x1}}4,`$ (3)
where $`x=h_\mathrm{p}\nu /kT_{\mathrm{CMB}}`$, $`n_\mathrm{e}`$ and $`T_\mathrm{e}`$ are the number density and temperature of electrons in clusters, respectively, and the integral is performed along the line of sight $`\chi `$. The angular power spectrum of the temperature fluctuation $`\mathrm{\Delta }T(𝜽)/T_{\mathrm{CMB}}`$ can be separated into the Poisson term $`C_{\mathrm{}}^P`$ and the clustering term $`C_{\mathrm{}}^C`$:

$$C_{\mathrm{}}=C_{\mathrm{}}^P+C_{\mathrm{}}^C.$$
(4)
Since the clustering term $`C_{\mathrm{}}^C`$ is several orders of magnitude smaller than the Poisson term $`C_{\mathrm{}}^P`$ at small angular scales $`\mathrm{}1000`$ where the SZ power becomes dominant in the CMB sky, we will neglect the contribution of source clustering to the SZ angular power spectrum. In the flat sky approximation, the Poisson term can be calculated by (e.g. Cole & Kaiser 1988)

$`C_{\mathrm{}}^P=g^2(x){\displaystyle 𝑑z\frac{d^2V}{dzd\mathrm{\Omega }}𝑑M\frac{d^2N(M,z)}{dMdV}|y_{\mathrm{}}(M,z)|^2},`$ (5)
in which $`d^2V/dzd\mathrm{\Omega }`$ denotes the comoving volume per unit redshift and per steradian, and $`y_{\mathrm{}}(M,z)`$ is the Fourier transform of $`y(\theta )`$:

$$y_{\mathrm{}}=2\pi y(\theta )J_0(\mathrm{}\theta )\theta 𝑑\theta ,$$
(6)
where $`J_0`$ is the Bessel function of zero order. We adopt the standard Press & Schechter (1974) formalism for the mass function of clusters

$$\frac{d^2N}{dMdV}=\sqrt{\frac{2}{\pi }}\frac{\overline{\rho }}{M}\frac{\delta _\mathrm{c}(z)}{\sigma ^2(M)}\frac{d\sigma (M)}{dM}\mathrm{exp}\left(\frac{\delta _c^2(z)}{2\sigma ^2(M)}\right),$$
(7)
where $`\overline{\rho }`$ is the mean mass density of the universe, $`\delta _c(z)`$ is the linear over-density of spherical collapse that virialized at redshift $`z`$, and $`\sigma (M)`$ is the present rms mass fluctuation within the top-hat filter $`M=4\pi \overline{\rho }R^3/3`$:

$$\sigma ^2(M)=\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2P(k)W^2(kR)𝑑k,$$
(8)
in which $`W(x)=3(\mathrm{sin}xx\mathrm{cos}x)/x^3`$ is the Fourier representation of the window function. We parameterize the power spectrum of fluctuation $`P(k)k^{n_\mathrm{s}}T^2(k)`$, and take the fit given by Bardeen et al. (1986) for the transfer function of adiabatic CDM model $`T(k)`$ with the shape parameter $`\mathrm{\Gamma }=0.21`$. The primordial power spectrum is assumed to be the Harrison-Zel’dovich case $`n_\mathrm{s}=1`$. The amplitude in the power spectrum is fixed using the rms fluctuation on an 8 $`h^1`$ Mpc scale, $`\sigma _8`$, which will be taken to be $`\sigma _8=0.9`$ in this study.

## 3 Contribution of cluster radio halos

### 3.1 Properties of cluster radio halos

Over the past few years, diffuse radio halos have been detected in a few tens of nearby, rich clusters. They often extend to a distance of $`1`$ Mpc from cluster centers, and have regular shape, low surface brightness and steep radio spectrum $`S_\nu \nu ^1`$. Moreover, observations suggest that there is a tight correlation between the radio power $`P_{1.4\mathrm{GHz}}`$ and the X-ray luminosity ($`L_\mathrm{X}`$) of host clusters (Liang et al. 2000; Govoni et al. 2001b), which can be well described by the following fitting formula (Enßlin & Röttgering 2002)

$$P_{1.4\mathrm{GHz}}=a_\nu 10^{24}h_{50}^2\mathrm{Watt}\mathrm{Hz}^1\left(\frac{L_\mathrm{X}}{10^{45}h_{50}^2\mathrm{erg}\mathrm{s}^1}\right)^{b_\nu },$$
(9)
in which $`a_\nu =2.78`$, $`b_\nu =1.94`$, and $`h_{50}`$ is the present Hubble constant in units of $`50`$ km s-1 Mpc-1. In order to set up a link between the X-ray luminosity and total mass of the cluster, we first employ the best-fit $`L_\mathrm{X}`$-$`T`$ relation of Xue & Wu (2000), based on the hitherto largest sample of nearby X-ray clusters,

$$\left(\frac{L_\mathrm{X}}{10^{43}h_{50}^2\mathrm{erg}\mathrm{s}^1}\right)=10^{0.032}\left(\frac{kT}{\mathrm{keV}}\right)^{2.79}.$$
(10)
It has been shown that such an $`L_\mathrm{X}`$-$`T`$ relation exhibits no significant evolution at least out to $`z0.8`$ (Rosati, Borgani & Norman 2002 and references therein). We then convert the X-ray temperature into the total cluster mass using virial theorem (Bryan & Norman 1998):

$$\frac{kT}{kT^{}}=\left(\frac{M}{10^{15}M_{}}\right)^{2/3},$$
(11)
in which

$$kT^{}=1.39\mathrm{keV}f_\mathrm{T}[h^2\mathrm{\Delta }_\mathrm{c}(z)E^2(z)]^{1/3},$$
(12)
$`h`$ is the Hubble constant in units of 100 km s-1 Mpc-1, $`f_\mathrm{T}`$ is the normalization factor, which will be taken to be $`f_\mathrm{T}=0.8`$ in this work, and $`\mathrm{\Delta }_\mathrm{c}(z)`$ is the overdensity of the virialized dark halos with respect to the critical density of the universe at redshift $`z`$. Equivalently, one may also employ the empirical $`L_\mathrm{X}`$-$`M`$ relation (Reiprich & Böhringer 1999) incorporating with a simple model of cosmic evolution for dark halos. Consequently, for a given halo of mass $`M`$ one is able to estimate its radio power at 1.4 GHz through

$$P_{1.4\mathrm{GHz}}=3.17\times 10^{20}h_{50}^2\mathrm{Watt}\mathrm{Hz}^1\left(\frac{kT^{}}{\mathrm{keV}}\right)^{5.41}\left(\frac{M}{10^{15}M_{}}\right)^{3.61}.$$
(13)
Radio observations indicate that cluster radio halos usually have a steep power index, which is also well motivated by synchrotron emission process. We adopt a spectral power index of $`\alpha =0.8`$ (e.g. Giovannini et al. 1993) to convert the radio power at 1.4 GHz to the radio emission at other wavelengths.

Furthermore, several observations have also reveled a similarity of the radio and X-ray morphologies in clusters (Deiss et al. 1997; Liang et al. 2000; Govoni et al. 2001a; Feretti et al. 2001). This enables us to model the radio surface brightness of clusters using the empirical fitting formula for X-ray surface brightness profile, i.e., the $`\beta `$ model,

$$\mathrm{\Delta }I_{\mathrm{RH}}(r)=\mathrm{\Delta }I_{\mathrm{RH0}}\left(1+\frac{r^2}{r_\mathrm{c}^2}\right)^{3\beta +0.5},$$
(14)
in which the normalization factor can be determined by the total radio power $`P_\lambda `$ at a given wavelength $`\lambda `$. Meanwhile, the electron number density in clusters follows

$$n_\mathrm{e}(r)=n_{\mathrm{e0}}\left(1+\frac{r^2}{r_\mathrm{c}^2}\right)^{3\beta /2}.$$
(15)
Here we introduce the universal baryon fraction $`f_\mathrm{b}=\mathrm{\Omega }_\mathrm{b}/\mathrm{\Omega }_\mathrm{M}`$ to fix the central electron density $`n_{\mathrm{e0}}`$: $`f_\mathrm{b}=M_{\mathrm{gas}}/M`$, where $`M_{\mathrm{gas}}`$ is the total gas mass within $`r_{\mathrm{vir}}`$. The gaseous halo (also the radio halo) is truncated at the virial radius $`r_{\mathrm{vir}}`$ defined by

$$M=\frac{4\pi }{3}r_{\mathrm{vir}}^3\mathrm{\Delta }_\mathrm{c}(z)\rho _\mathrm{c}(z),$$
(16)
where $`\rho _\mathrm{c}(z)`$ is the critical mass density of the universe at $`z`$. We adopt a constant value of $`2/3`$ for the $`\beta `$ parameter, and assume that the core radius $`r_\mathrm{c}`$ of gas distribution is proportional to the size of the host halo: $`r_\mathrm{c}r_{\mathrm{vir}}`$. We specify the proportionality constant such that $`r_\mathrm{c}=0.15`$ $`h^1`$ Mpc for a massive cluster of $`M=10^{15}h^1M_{}`$ at $`z=0`$ (cf. Komatsu & Kitayama 1999).

Finally, radio halos are preferentially found in clusters that demonstrate violent mergers and substructures, which provide an accelerating mechanism for the relativistic particles responsible for radio emission (e.g. Buote 2001). Yet, the fraction of clusters that possess giant radio halos is still uncertain. Following Enßlin & Röttgering (2002), we adopt a constant radio halo fraction of $`f_{\mathrm{RH}}=1/3`$. Nevertheless, we will also demonstrate the result for a fraction of unity for comparison. Recall that cluster mergers may occur more frequently in the past than at the present.

### 3.2 Modification to the SZ angular power spectrum

In the Rayleigh-Jeans region, the presence of radio halos can dilute the thermal SZ signals towards clusters. The significance of the effect depends critically on working frequencies. Recall that radio emission of clusters has a steep power spectrum of $`\alpha 0.8`$. At very low frequencies, the decrements of the CMB temperature due to the SZ effect could be completely compensated (Enßlin 2002) by cluster radio halos. Here we demonstrate the radio halo contamination using an observing frequency of $`\nu =30`$ GHz. The change in the CMB temperature is related to the intensity fluctuations by

$$\frac{\mathrm{\Delta }T}{T_{\mathrm{CMB}}}=\frac{(e^x1)^2}{x^4e^x}\frac{\mathrm{\Delta }I}{I_0},$$
(17)
where $`I_0=2(kT_{\mathrm{CMB}})^3/(h_\mathrm{p}c)^2`$. $`\mathrm{\Delta }I`$ is composed of two components, the thermal SZ signals $`\mathrm{\Delta }I_{\mathrm{SZ}}`$ and the cluster radio surface brightness $`\mathrm{\Delta }I_{\mathrm{RH}}`$,

$$\mathrm{\Delta }I=\mathrm{\Delta }I_{\mathrm{SZ}}+\mathrm{\Delta }I_{\mathrm{RH}},$$
(18)
or equivalently, we can separated the temperature fluctuation into

$$\frac{\mathrm{\Delta }T}{T_{\mathrm{CMB}}}=\frac{\mathrm{\Delta }T_{\mathrm{SZ}}}{T_{\mathrm{CMB}}}+\frac{\mathrm{\Delta }T_{\mathrm{RH}}}{T_{\mathrm{CMB}}}.$$
(19)
If the fraction of cluster radio halos is $`f_{\mathrm{RH}}`$, then the observed SZ power spectrum at small angular scales is

$$C_{\mathrm{}}=𝑑z\frac{d^2V}{dzd\mathrm{\Omega }}𝑑M\frac{d^2N(M,z)}{dMdV}\left[(1f_{\mathrm{RH}})|g(x)y_{\mathrm{}}(M,z)|^2+f_{\mathrm{RH}}|y_{\mathrm{}}^{}(M,z)|^2\right],$$
(20)
where $`y_{\mathrm{}}^{}(M,z)`$ is the Fourier transform of $`(\mathrm{\Delta }T_{\mathrm{SZ}}+\mathrm{\Delta }T_{\mathrm{RH}})/T_{\mathrm{CMB}}`$.

In order to highlight the effect of cluster radio halos on the SZ angular power spectrum, we demonstrate in Figure 1 the relative correction to the thermal SZ angular power spectrum $`|\mathrm{\Delta }C_{\mathrm{}}|/C_{\mathrm{}}`$, along with the contributions of clusters in different redshift ranges. At $`1000\mathrm{}8000`$, the radio halo contamination in the SZ angular power spectrum arises mainly from clusters at $`0.1<z<0.5`$. High-redshift clusters come into effect only at very small angular scales $`\mathrm{}>10^4`$, while large angular scales are dominated by low-redshift clusters. In Figure 2 we show the contributions of cluster radio halos with different cluster masses to the thermal SZ angular power spectrum. At large angular scales below $`\mathrm{}500`$, the radio contamination is governed by nearby, massive clusters, which can be easily removed from the SZ map because of their rarity. At smaller angular scales beyond $`\mathrm{}=1000`$, the contamination arises mainly from moderately rich clusters at higher redshifts $`z>0.2`$. Contributions from poor clusters and groups of masses less than $`2\times 10^{14}`$ $`M_{}`$ are negligible because of the absence of radio halos in these systems. To provide a quantitative illustration of the modifications due to cluster radio halos, we list in Table 1 the relative corrections $`|\mathrm{\Delta }C_l|/C_l`$ for a set of $`\mathrm{}`$ evaluated by Equation (20), assuming the radio halo fractions of $`f_{\mathrm{RH}}=1/3`$ and $`f_{\mathrm{RH}}=1`$, respectively. We have also shown in Table 1 how our estimates are affected by the uncertainties associated with the measurements of radio power, X-ray luminosities and total masses of clusters. To do this, we simply vary the radio power $`P_\lambda `$, X-ray luminosity $`L_X`$ and normalization factor $`f_T`$ in the virial theorem by $`20\%`$, $`20\%`$ and $`30\%`$, respectively, which are roughly representative of the uncertainties in current determinations of these quantities. It turns out that the modifications to the thermal SZ angular power spectrum from cluster radio halos are only minor, and the maximum correction is no more than $`13\%`$ for $`\mathrm{}>2000`$ even in the extreme case of where $`f_{\mathrm{RH}}=1`$ and $`f_T=1`$. Major uncertainty in our calculation arises from the mass estimate of clusters characterized by the normalization factor $`f_T`$, which could differ by a factor of $`2`$ among different simulations and observations (e.g. Bryan & Norman 1998; Huterer & White 2002).

## 4 Contribution of cluster radio galaxies

### 4.1 Distribution of cluster radio galaxies

A critical issue in the estimate of confusion of cluster radio galaxies in the measurement of the SZ angular power spectrum is to properly handle the distribution of cluster radio galaxies. This includes the luminosity function and radial profile of radio galaxies in clusters. We adopt the univariate radio luminosity function given by Ledlow & Owen (1996), based on a statistical complete sample of 188 radio galaxies in Abell clusters within $`z=0.09`$ and at $`\nu =1.4`$ GHz. The univariate radio luminosity function $`f_{1.4\mathrm{GHz}}`$ can be well fitted by a broken power law

$$\mathrm{log}f_{1.4\mathrm{GHz}}=\mathrm{log}\frac{dN_{\mathrm{radio}}/d\mathrm{log}P_{1.4\mathrm{GHz}}}{N_{\mathrm{gal}}(0.3r_{\mathrm{Abell}})}=a+b\mathrm{log}P_{1.4\mathrm{GHz}},$$
(21)
where $`N_{\mathrm{gal}}(0.3r_{\mathrm{Abell}})`$ is the total number of galaxies brighter than $`R=20.5`$ within 0.3 Abell radius of the cluster center, and $`(a,b)`$ are proportionality constants with ($`0.150,1.766`$) and ($`1.433,33.667`$) for $`P_{1.4\mathrm{GHz}}<10^{24.8}`$ W Hz-1 and $`P_{1.4\mathrm{GHz}}>10^{24.8}`$ W Hz-1, respectively. In order to convert the above radio luminosity function at $`\nu =1.4`$ GHz into the ones ($`f_\nu `$) at other frequencies, we use the spectral index distribution fitted by Lin, Chiueh, & Wu (2002) for a sample of 64 radio sources in 56 clusters observed by Cooray et al. (1998), in which the mean spectral index is $`\alpha =0.71`$. This allows us to take the dispersion of spectral index into account. The luminosity-weighted cluster radio luminosity function takes the form of

$$\mathrm{log}(P_\nu f_\nu )=\{\begin{array}{cc}C_1(\nu )+\gamma _1(\nu )\mathrm{log}P_\nu ;\hfill & \mathrm{log}P<\mathrm{log}P_\mathrm{b}\hfill \\ C_2(\nu )+\gamma _2(\nu )\mathrm{log}P_\nu ;\hfill & \mathrm{log}P>\mathrm{log}P_\mathrm{b}\hfill \end{array},$$
(22)
where

$`C_1(\nu )`$ $`=`$ $`2.680.95e^{0.015\nu }`$
$`\gamma _1(\nu )`$ $`=`$ $`0.85e^{0.0027\nu }+1.49\times 10^3\nu `$
$`C_2(\nu )`$ $`=`$ $`30.42+3.29e^{0.025\nu }`$
$`\gamma _2(\nu )`$ $`=`$ $`0.43e^{0.0022\nu }6.14\times 10^4\nu ,`$ (23)
and the break power $`P_\mathrm{b}`$ has a weak dependence on frequency: log$`P_\mathrm{b}(\nu )=23.860.00396\nu `$.

We assume that the number density profile of galaxies in clusters follows the King model

$$n_{\mathrm{gal}}(r)=n_{\mathrm{gal}}(0)\left(1+\frac{r^2}{r_c^2}\right)^{3/2}.$$
(24)
The surface number density is thus $`\mathrm{\Sigma }_{\mathrm{gal}}(r)=n_{\mathrm{gal}}(r)𝑑\chi `$. A good approximation of the core radius is $`r_c=0.1r_{\mathrm{vir}}`$. In order to fix the central density $`n_{\mathrm{gal}}(0)`$, we introduce the so-called halo occupation distribution, $`N_{\mathrm{gal}}(M)=n_{\mathrm{gal}}(r)4\pi r^2𝑑r`$. We adopt the best fit analytic formulae of Sheth & Diaferio (2001) for the spiral \[$`N_{\mathrm{gal},\mathrm{S}}(M)`$\] and elliptical \[$`N_{\mathrm{gal},\mathrm{E}}(M)`$\] galaxies based on the GIF simulations (Kauffmann et al. 1999):

$`N_{\mathrm{gal},\mathrm{S}}(M)`$ $`=`$ $`(M/M_\mathrm{S})^{\alpha _S}+0.5e^{4[\mathrm{log}(M/10^{11.75}M_{\mathrm{}})]^2}`$
$`N_{\mathrm{gal},\mathrm{E}}(M)`$ $`=`$ $`(M/M_\mathrm{E})^{\alpha _E}e^{(2\times 10^{11}M_{\mathrm{}}/M)^2}`$
$`N_{\mathrm{gal}}(M)`$ $`=`$ $`N_{\mathrm{gal},\mathrm{S}}(M)+N_{\mathrm{gal},\mathrm{E}}(M),`$ (25)
where $`M_\mathrm{S}=7\times 10^{13}h^1M_{\mathrm{}}`$, $`\alpha _\mathrm{S}=0.9`$, $`M_\mathrm{E}=3\times 10^{12}h^1M_{\mathrm{}}`$, and $`\alpha _\mathrm{E}=0.75`$.

Another parameter we need to fix before we proceed to our numerical calculation is $`N_{\mathrm{gal}}(0.3r_{\mathrm{Abell}})`$, the galaxies brighter than $`R=20.5`$ within 0.3 Abell radius $`r_{\mathrm{Abell}}`$. Instead of evaluating $`N_{\mathrm{gal}}(0.3r_{\mathrm{Abell}})`$ directly, we turn to the fraction of galaxies with absolute magnitude brighter than $`R_{\mathrm{limit}}=20.5`$, $`f_{\mathrm{gal}}`$. This can be achieved by employing the luminosity functions of galaxies in clusters, $`\mathrm{\Phi }_S`$ and $`\mathrm{\Phi }_E`$. We adopt the best fit Schechter forms in the $`r^{}`$ band for early and late type galaxies derived from the Sloan Digital Sky Survey (Goto et al. 2002). The best fit parameters in Goto et al. (2002) are based on galaxies within $`0.75`$ $`h_{0.7}^1`$ Mpc, which is only slightly larger than $`0.3r_{\mathrm{Abell}}`$. We convert the $`r^{}`$ magnitude to the standard Johnson-Morgan-Cousins system using $`r^{}=R+0.16(VR)+0.13`$ (Fukugita et al. 1996) and then assign the colors to early and late type galaxies according to the tabulated values given by Fukugita, Shimasaku & Ichikawa (1995). Moreover, we choose a lower magnitude cutoff of $`R_{\mathrm{cut}}=18.5`$ to guarantee the validity of the above adopted halo occupation distribution equation (25) which accounts for galaxies brighter than $`V=17.7+5\mathrm{log}h`$. Taking all these factors into account, we have

$`f_{\mathrm{gal}}={\displaystyle \frac{_{r_{\mathrm{S},\mathrm{limit}}^{}}^{\mathrm{}}\mathrm{\Phi }_S(r^{})𝑑r^{}+_{r_{\mathrm{E},\mathrm{limit}}^{}}^{\mathrm{}}\mathrm{\Phi }_E(r^{})𝑑r^{}}{_{r_{\mathrm{S},\mathrm{cut}}^{}}^{\mathrm{}}\mathrm{\Phi }_S(r^{})𝑑r^{}+_{r_{\mathrm{E},\mathrm{cut}}^{}}^{\mathrm{}}\mathrm{\Phi }_E(r^{})𝑑r^{}}}.`$ (26)
The total radio power of all galaxies in a cluster can be obtained by extrapolating the radio luminosity function of Ledlow & Owen (1996) determined with $`0.3r_{\mathrm{Abell}}`$ to the whole cluster region out to virial radius, which yields

$`P_{\mathrm{tot},\nu }`$ $`=`$ $`f_{\mathrm{gal}}\left[{\displaystyle \frac{N_{\mathrm{gal}}(M)}{N_{\mathrm{gal}}(0.3r_{\mathrm{Abell}})}}\right]{\displaystyle P_\nu (dN_{\mathrm{radio}}/d\mathrm{log}P_\nu )d\mathrm{log}P_\nu }`$ (27)
$`=`$ $`f_{\mathrm{gal}}N_{\mathrm{gal}}(M){\displaystyle P_\nu f_\nu d\mathrm{log}P_\nu }.`$
The total surface brightness of all radio galaxies measured at frequency $`\nu `$ in direction $`𝜽`$ toward a cluster of mass $`M`$ and redshift $`z`$ is thus

$`\mathrm{\Delta }I_{\mathrm{RG}}(𝜽)`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{\Sigma }_{\mathrm{gal}}(𝜽)D_A^2(z)}{N_{\mathrm{gal}}(M)}}\right]\left[{\displaystyle \frac{(1+z)P_{\mathrm{tot},\nu }}{4\pi D_L^2(z)}}\right]`$ (28)
$`=`$ $`{\displaystyle \frac{\mathrm{\Sigma }_{\mathrm{gal}}(𝜽)f_{\mathrm{gal}}P_\nu f_\nu 𝑑\mathrm{log}P_\nu }{4\pi (1+z)^3}},`$
where $`D_A`$ and $`D_L`$ are the angular diameter distance and luminosity distance to the cluster, respectively. This adds to the SZ signals such that the observed radio intensity is $`\mathrm{\Delta }I=\mathrm{\Delta }I_{\mathrm{SZ}}+\mathrm{\Delta }I_{\mathrm{RG}}`$. The corresponding angular power spectrum reads

$$C_{\mathrm{}}=𝑑z\frac{d^2V}{dzd\mathrm{\Omega }}𝑑M\frac{d^2N(M,z)}{dMdV}|y_{\mathrm{}}^{\prime \prime }(M,z)|^2,$$
(29)
where $`y_{\mathrm{}}^{\prime \prime }(M,z)`$ is the Fourier transform of $`(\mathrm{\Delta }T_{\mathrm{SZ}}+\mathrm{\Delta }T_{\mathrm{RG}})/T_{\mathrm{CMB}}`$, and the temperature fluctuation $`\mathrm{\Delta }T_{\mathrm{RG}}`$ is related to $`\mathrm{\Delta }I_{\mathrm{RG}}`$ following equation (17).

### 4.2 Modification to the SZ angular power spectrum

In a similar way to Figure 1, we plot in Figure 3 the relative corrections to the SZ angular power spectrum $`|\mathrm{\Delta }C_{\mathrm{}}|/C_{\mathrm{}}`$ because of the presence of radio galaxies in clusters, together with the contributions of cluster radio galaxies in different redshift ranges. At $`l>1000`$, the distant radio galaxies at $`z>0.2`$, which could be regarded as the unresolved radio sources, dominate the correction term. We list in Table 2 the relative corrections for a set of small scales, together with the uncertainties corresponding to $`20\%`$ variation of the radio luminosity function. It turns out that the confusion of unresolved radio galaxies in clusters in the measurement of the thermal SZ effect is less than $`7.5\%`$ on all scales beyond $`l>2000`$. This is somewhat smaller than the estimate ($`30\%`$) of Holder (2002), who evaluated the effect of radio point sources by removing a certain central region ($`1^{}`$) from all clusters. The latter could simultaneously lead to an underestimate of the thermal SZ effect in the centers of all clusters especially for high-redshift clusters.

## 5 Conclusions

Contamination of radio sources has been a major concern in the measurement of thermal SZ effect, which include radio halos and galaxies of clusters, and clustering of radio galaxies tracing large-scale structures of the universe. The latter comes into effect only on large angular scales and appears to be relatively small as compared with the thermal SZ effect (Toffolatti et al. 1998). We have thus concentrated in this paper on the role of cluster radio halos and galaxies. The effect of cluster radio halos on the measurement of the SZ angular power spectrum has previously remained unknown, while the confusion of the cluster radio galaxies has been recently shown to be very significant (Holder 2002). Using an analytic model for both radio sources and dark matter of clusters, we have performed a quantitative analysis of the extent to which cluster radio halos and galaxies may contaminate the measurement of thermal SZ angular power spectrum. We have found that while the cluster radio sources can indeed cause some confusion in the measurement of SZ angular power spectrum, the relative correction $`|\mathrm{\Delta }C_l|/C_l`$ is well within $`20\%`$ at $`l>2000`$ for an observing frequency of $`\nu =30`$ GHz. The modification becomes even smaller for high-frequency observations. We thus conclude that there is probably no need to take the confusion of cluster radio sources into account in the thermal SZ measurement, if our observing frequencies are chosen to be $`\nu >30`$ GHz, unless the thermal SZ power spectrum is used for the purpose of precision cosmology.

We gratefully acknowledge Torsten Enßlin for useful discussion and help, and an anonymous referee for useful suggestions. This work was supported by the National Science Foundation of China, under Grant No. 10233040, and the Ministry of Science and Technology of China, under Grant No. NKBRSF G19990754.


format: text

text:
# Hot Molecular Cores and High-Mass Star Formation

## 1. Introduction

Named after the Orion Hot Core, hot molecular cores are small ($`d<`$0.1 pc) pockets of warm ($`>100`$ K) dense ($`>10^6`$ cm-3) molecular gas, which show very rich submillimeter spectra. This definition includes both internally and externally heated clumps near stars of all masses, but the focus of this review is on the case of internal heating by young high-mass stars, which form the majority of known hot cores. This review only outlines basic properties of hot cores with an emphasis on recent developments. Recent reviews on the chemistry of star-forming regions, including hot cores, have been given by van Dishoeck & van der Tak (2000) and Langer et al. (2000); high-mass star formation is discussed by Kurtz et al. (2000), Churchwell (2002) and Garay (this volume).

## 2. The embedded phase of high-mass star formation

Stars of masses $`>`$8 M spend a significant fraction of their lifetime, $`>`$10%, embedded in their natal molecular clouds. This embedded phase can be subdivided into four observationally different groups of objects: (1) The initial conditions of high-mass star formation should be large mass reservoirs which are local temperature minima and density maxima. Candidate objects in this phase have been found through infrared (Egan et al. 1998; Carey et al. 1998) and submillimeter (Motte et al. 2003) observations. (2) High-mass protostellar objects, where the central star is surrounded by a massive envelope with a centrally peaked temperature distribution (van der Tak et al. 2000a, Sridharan et al. 2002). Close to the star, ices start evaporating off the grains. (3) Hot cores, with large masses of warm and dense molecular gas, and large abundances of complex organic molecules such as CH3OCH3 and CH3OCHO. These species are thought to be the chemical ‘daughter’ products of reactions of evaporated ice mantle components, their ‘parents’. (4) Small ($`d<`$0.1 pc) pockets of ionized gas, called ultracompact H ii regions. The ordering of these phases from ‘cold’ to ‘hot’ may well be an evolutionary sequence, where the ratio of envelope mass to luminosity decreases, although the mass of the original cloud may also play a role. Phases 2 and 3 are accompanied by massive molecular outflows (Beuther et al. 2002), phases 1–3 by CH3OH masers (Walsh et al. 2001), and H2O masers can be found in all four phases.

The separation between these four phases is not as clean as in the case of low mass star formation. Massive stars form in groups and the four above phases often occur right next to each other, making it hard to assign a single age to the region. An example is G10.47, where VLA images show peaks in both the NH3 (4,4) line, which traces warm molecular gas, and in the 1.3 cm continuum, which traces ionized gas (Cesaroni et al. 1998). In W49, Wilner et al. (2001) imaged CH3CN line and dust continuum emission at 1.4 mm using the BIMA interferometer. The number of hot cores in W49 is about half that of ultracompact H ii regions, indicating that the hot core phase lasts $``$half the lifetimes of ultracompact H ii regions of $``$105 yr. This value agrees well with estimates based on chemical modeling (e.g., Rodgers & Charnley 2001) and with dynamical timescales of molecular outflows. This agreement, together with the lack of spatial ordering of different phases within the region, argues against triggering of star formation by an external event.

## 3. Chemical processes during high-mass star formation

In the earliest, dense, cold phases prior to star formation, atoms and molecules freeze out onto grain surfaces, which act as catalysts for neutral-neutral reactions. The composition of the resulting grain mantles can be probed by mid-infrared spectroscopy. This technique has developed greatly with the ISO mission and the advent of high-quality spectrographs on 8m-class telescopes. The derived abundances (Table 1) exceed by orders of magnitude the values that can be produced by reactions in the gas phase, with CO as the only exception. The data thus make a strong case for an active grain surface chemistry. Additional evidence comes from the observed isotopic composition (D/H ratio) of evaporated H2CO and CH3OH (Charnley et al. 1997).

However, the detailed reaction mechanisms of surface chemistry are not well understood. Laboratory experiments on hydrogenation of CO give conflicting results (Watanabe & Kouchi 2002; Hiraoka et al. 2002), so that the origin of solid CH3OH is not clear. Also, observational limits on solid HDO are inconsistent with the predictions of surface chemistry models (Dartois et al. 2003).

Considering elemental abundances, it is clear that the observational inventory of the solid state is incomplete. The likely reason is that the major carriers of some elements are undetectable through infrared spectroscopy. For example, the only N-bearing species is NH3, which contains only $``$5% of elemental nitrogen. The remaining nitrogen could be in the (undetectable) form of N2. Similarly, OCS does not make up the full sulphur budget. Submillimeter observations suggest trace amounts ($`<`$1%) of solid H2S or SO2 evaporating, but observations with ISO rule out both species, as well as atomic S, as major sulphur carriers (van der Tak et al. 2003).

After evaporation into the gas phase, the molecules of Table 1 start a complex reaction scheme, building up long carbon chains with many different functional groups. Spectroscopy of hot cores at millimeter wavelengths has revealed considerable chemical complexity: the BIMA survey of CH3COOH by Remijan et al. (2003) provides a recent example. Probably even larger molecules are formed, but these are very hard to detect, because of spectroscopic confusion, and because individual spectral lines become weaker as the partition function increases for ever larger species.

The formation of large molecules in hot cores is mostly due to NH3 and CH3OH, which, after acquiring a proton from H$`{}_{}{}^{+}{}_{3}{}^{}`$, H3O+ or HCO+, can form chains. The HCN and C2H5OH molecules could play similar roles, but so far, only upper limits on their solid state abundances have been obtained. An important determinant of the chemistry is therefore the competition for protons between NH3 and CH3OH, which both have a high proton affinity. Rodgers & Charnley (2001) show model results for various values of the NH3/CH3OH ratio, which observations do not yet constrain well.

Chemical diversity can also be the result of high-temperature neutral-neutral reactions. At temperatures $`>`$300 K, the reactions of O and OH with H2 drive all the available oxygen into H2O. Cores at high temperatures are therefore expected to be rich in nitrogen-bearing compounds (Rodgers & Charnley 2001). Although in the models, age differences also cause chemical diversity, observations seem to support a link between temperature and nitrogen-oxygen differentiation. For example, IRAM Plateau de Bure data show that the peak position of nitrogen-bearing species in W3 (OH) has a higher temperature than that of oxygen-bearing species (Wyrowski et al. 1999). Additional support comes from JCMT spectroscopy of the HCN $`J=98`$ line at 797 GHz (Boonman et al. 2001). The data show that the HCN abundance in the envelope of the massive young star AFGL 2591 is not constant with radius, but is enhanced by a factor of 100 at small radii where $`T>300`$ K.

Shocks do not seem as important for the chemistry of high-mass star-forming regions as they are in the low-mass case (e.g., Bachiller et al. 2001). Still, two molecules suggest that the envelopes of young high-mass stars may have been processed by shocks, although neither one conclusively. First, ISO-SWS data show low ratios of gas-phase to solid-state CO2 ($`<0.1`$) and the CO2 gas phase abundance remains low through the 100–300 K temperature regime (Boonman et al. 2003). After evaporating off grains, CO2 must be promptly destroyed, which shocks can do in reactions with H (Charnley & Kaufman 2000) or perhaps H2 (Doty et al. 2002; Talbi & Herbst 2002). Second, Hatchell & Viti (2002) measured NS/CS ratios of 0.02-0.05 in a sample of six hot cores and interpreted these as evidence for shocks. The main reactions to form NS require SH and NH which are produced from OH. Shocks use OH to form H2O and suppress the production of NS. The values of NS/CS=0.001 – 0.01 measured in warm envelopes (van der Tak et al. 2003) may indicate that shocks play a role too. However, for $`t=3\times 10^4`$ yr, the Doty et al. (2002) model of envelope chemistry also predicts NS/CS $`10^3`$, so this ratio cannot be used to demonstrate the influence of shocks.

## 4. Hot cores with cool envelopes

Before the onset of an ultracompact H ii region, young high-mass stars are surrounded by envelopes where ice evaporation occurs in the inner parts. Time-dependent chemical models for this phase were constructed by Doty et al. (2002). These models employ standard reaction networks, and the initial conditions resemble those of dark clouds in the outer envelope ($`T<100`$ K). In the inner envelope ($`T>100`$ K), the initial abundances of the species in Table 1 are enhanced, corresponding to evaporation of ice mantles.

Figure 1 shows results for SO and SO2 as an example. The evaporation of ice mantles leads to pronounced abundance enhancements (‘jumps’) for several molecules. Observational evidence for such jumps is mounting: H2CO, CH3OH, SO, SO2 and OCS show jumps at $`T100`$ K (van der Tak et al. 2000b; Schöier et al. 2002; van der Tak et al. 2003), likely due to ice evaporation. In contrast, HCN is enhanced at $`T300`$ K (Boonman et al. 2001) in the source AFGL 2591. Such an enhancement is a predicted effect of the O$``$H2O reaction and suggests that in this source, high-temperature gas-phase reactions have started, marking the birth of the hot core phase. However, the abundance jumps have so far been inferred from detailed modeling of single-dish spectra. Millimeter-wave interferometers could resolve the jumps spatially, which would be an important test of current models of envelope chemistry.

A subset of these ‘warm envelopes’, where the high mid-infrared brightness allows detection of H$`{}_{}{}^{+}{}_{3}{}^{}`$ rovibrational absorption lines, has been used to derive the cosmic-ray ionization rate (van der Tak & van Dishoeck 2000).111In this paper, the values in the last column of Table 1 were scaled the wrong way around. The conclusions of the paper are not affected though. This process is the main ionization mechanism of dense molecular clouds, regulating for example the influence of magnetic fields on their dynamics. Van der Tak & van Dishoeck compared observations of H$`{}_{}{}^{+}{}_{3}{}^{}`$ absorption and H13CO+ emission to a small chemical network coupled to detailed models of the sources’ structure, to find a mean $`\zeta `$ of $`2.6\times 10^{17}`$ s-1. The source-to-source spread of a factor of 2–3 may be correlated with $`N`$(H2), suggesting that absorption of cosmic rays plays a role. The value of $`\zeta =5.6\times 10^{17}`$ s-1 for one particular source was verified by Doty et al. (2002) using a detailed chemical model.222In this paper, the predicted $`N`$(H$`{}_{}{}^{+}{}_{3}{}^{}`$) should read $`5\times 10^{13}`$ cm-2 rather than $`5\times 10^{14}`$ cm-2.

McCall et al. (2003) infer a much higher value of $`\zeta 1\times 10^{15}`$ s-1 for the diffuse cloud towards $`\zeta `$ Per. More detailed modeling of the chemistry of H$`{}_{}{}^{+}{}_{3}{}^{}`$ and OH and HD, which also depend on $`\zeta `$, give similar results within a factor of 2–3, $``$10$`\times `$ higher than the results for the massive protostars. McCall et al. propose low-energy cosmic rays that penetrate diffuse but not dense clouds. High-resolution observations of HD in a large sample of diffuse clouds with FUSE will help to test this conclusion.

A much lower value of $`\zeta =0.6\times 10^{17}`$ s-1 was obtained by Caselli et al. (2002) for the dense core of L1544. This result needs to be checked in view of the large depletion of CNO-bearing molecules implied by the large abundance of H2D+ towards this source (Caselli et al. 2003). If confirmed, this low value of $`\zeta `$ provides additional evidence for variations in the ionization rates of molecular clouds caused by absorption of cosmic rays. However, more sources need to be studied before local variations in cosmic-ray production can be ruled out.

## 5. Low-mass hot cores

Although the existence of warm ($`>`$100 K) gas in the inner envelopes of low-mass protostars was well-known, it was unclear whether such objects achieved the same chemical complexity as their high-mass counterparts (e.g., Schöier et al. 2002). This situation changed when Cazaux et al. (2003) used the IRAM 30m telescope to detect both ’parent’ and ’daughter’ type species in the envelope of the low-mass protostar IRAS 16293–2422. The abundances are similar to or perhaps even higher than those in the Orion hot core; the assumed $`N`$(H2) of the inner envelope of IRAS 16293 and the size of the Orion hot core are critical to this comparison. The high deuterium fractionation of CH3OH measured by Parise et al. (2002) for IRAS 16293 may also be reproduced by grain surface chemistry models if the atomic D/H ratios from the model of Roberts et al. (2003) are adopted. However, the impact of the rate coefficients of key H$``$D reactions measured by Gerlich et al. (2002) on the Roberts et al. model needs to be studied.

The initial chemical conditions of low- and high-mass hot cores may well be similar. Using ISAAC on the VLT, Pontoppidan et al. (2003) detected CH3OH ice in three low-mass young stellar objects. The solid CH3OH abundances of 15–25% relative to H2O are comparable to those in the most methanol-rich massive sources. However, it is unclear if evaporation of these ice mantles can drive the rich gas-phase chemistry of IRAS 16293. The measured collapse motions in the envelope imply a travel time for gas through the hot core region of only a few hundred years, which is too short to build up high abundances of complex molecules (Schöier et al. 2002; Rodgers & Charnley 2003).

## 6. Future directions

Interferometry at (sub)mm wavelengths allows imaging of dust continuum and molecular lines at 1′′ or slightly higher resolution. The contribution by Beuther in this volume shows a recent example obtained with the Sub Millimeter Array. Such data are useful to constrain the temperature and density structure of high-mass star-forming cores, and to determine the multiplicity of their power sources (Wyrowski et al. 1999). However, current instruments can spatially resolve such structure only in a few nearby objects, and searches for kinematic structure in hot cores have just started (e.g., Cesaroni et al. 1998). Efforts are also hampered by limited sensitivity which at the moment prohibits observation of multiple transitions and of isotopic species, which is necessary to determine the optical depth and the excitation of the lines. In the future, ALMA will provide increases in angular resolution and sensitivity by factors of 10 or more.

Observational progress will also come from the HIFI instrument onboard Herschel, which will perform spectral line surveys at heterodyne resolution. Unhindered by the Earth’s atmosphere, HIFI will not only probe H2O (and maybe O2) directly, but also observe lines that are blocked by atmospheric H2O (especially the ground state transitions of hydride molecules). The calibration will be much better than can be achieved from the ground. However, the success of the HIFI mission depends on the availability of spectroscopic and collisional data for astrophysically interesting molecules at frequencies up to 2 THz. Sideband deconvolution is another critical issue (e.g., Comito & Schilke 2002). These caveats are especially important for hot cores with their rich line spectra.

### Acknowledgments.

The author thanks Fredrik Schöier, Friedrich Wyrowski and Ewine van Dishoeck for useful discussions.

## References

Bachiller, R., Pérez Gutiérrez, M., Kumar, M. S. N., & Tafalla, M. 2001, A&A, 372, 899

Beuther, H., Schilke, P., Sridharan, T. K., Menten, K. M., Walmsley, C. M., & Wyrowski, F. 2002, A&A, 383, 892

Boogert, A. C. A., Helmich, F. P., van Dishoeck, E. F., Schutte, W. A., Tielens, A. G. G. M., & Whittet, D. C. B. 1998, A&A, 336, 352

Boonman, A., Stark, R., van der Tak, F., van Dishoeck, E., van der Wal, P., Schäfer, F., de Lange, G., & Laauwen, W. 2001, ApJ, 553, L63

Boonman, A. M. S., van Dishoeck, E. F., Lahuis, F., & Doty, S. D. 2003, A&A, 399, 1063

Carey, S. J., Clark, F. O., Egan, M. P., Price, S. D., Shipman, R. F., & Kuchar, T. A. 1998, ApJ, 508, 721

Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37

Caselli, P., Walmsley, C. M., Zucconi, A., Tafalla, M., Dore, L., & Myers, P. C. 2002, ApJ, 565, 344

Cazaux, S., Tielens, A.G.G.M., Ceccarelli, C. et al. 2003, ApJ 593, L51

Cesaroni, R., Hofner, P., Walmsley, M., & Churchwell, E. 1998, A&A, 331, 709

Charnley, S. B. & Kaufman, M. J. 2000, ApJ, 529, L111

Charnley, S. B., Tielens, A. G. G. M., & Rodgers, S. D. 1997, ApJ, 482, L203

Churchwell, E. 2002, ARA&A, 40, 27

Comito, C. & Schilke, P. 2002, A&A, 395, 357

Dartois, E., d’Hendecourt, L., Thi, W., Pontoppidan, K. M., & van Dishoeck, E. F. 2002, A&A, 394, 1057

Dartois, E., Schutte, W., Geballe, T. R., Demyk, K., Ehrenfreund, P., & d’Hendecourt, L. 1999, A&A, 342, L32

Dartois, E., Thi, W.-F., Geballe, T. R., Deboffle, D., d’Hendecourt, L., & van Dishoeck, E. 2003, A&A, 399, 1009

Doty, S. D., van Dishoeck, E. F., van der Tak, F. F. S., & Boonman, A. M. S. 2002, A&A, 389, 446

Egan, M. P., Shipman, R. F., Price, S. D., Carey, S. J., Clark, F. O., & Cohen, M. 1998, ApJ, 494, L199

Gerakines, P. A., Whittet, D.C.B., Ehrenfreund, P. et al. 1999, ApJ, 522, 357

Gerlich, D., Herbst, E., & Roueff, E. 2002, PSS, 50, 1275

Gibb, E. L. et al. 2000, ApJ, 536, 347

Hatchell, J. & Viti, S. 2002, A&A, 381, L33

Hiraoka, K., Sato, T., Sato, S., Sogoshi, N., Yokoyama, T., Takashima, H., & Kitagawa, S. 2002, ApJ, 577, 265

Kurtz, S., Cesaroni, R., Churchwell, E., Hofner, P., & Walmsley, C. M. 2000, Protostars and Planets IV, 299

Langer, W., van Dishoeck, E., Bergin, E., Blake, G., Tielens, A., Velusamy, T., & Whittet, D. 2000, Protostars and Planets IV, 29

McCall, B. J., Honeycutt, A.J., Saykally, R.J. et al. 2003, Nature, 422, 500

Motte, F., Schilke, P., & Lis, D. C. 2003, ApJ, 582, 277

Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, ApJ, 479, 839

Parise, B., Ceccarelli, C., Tielens, A.G.G.M. et al. 2002, A&A, 393, L49

Pontoppidan, K. M., Dartois, E., van Dishoeck, E. F., Thi, W.-F., & d’Hendecourt, L. 2003, A&A, 404, L17

Remijan, A., Snyder, L. E., Friedel, D. N., Liu, S.-Y., & Shah, R. Y. 2003, ApJ, 590, 314

Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41

Rodgers, S. D. & Charnley, S. B. 2001, ApJ, 546, 324

Rodgers, S. D. & Charnley, S. B. 2003, ApJ, 585, 355

Schöier, F., Jørgensen, J., van Dishoeck, E., Blake, G. 2002, A&A 390, 1001

Schutte, W., Boogert, A., Tielens, A. et al. 1999, A&A, 343, 966

Sridharan, T. K., Beuther, H., Schilke, P., Menten, K. M., & Wyrowski, F. 2002, ApJ, 566, 931

Talbi, D. & Herbst, E. 2002, A&A, 386, 1139

Tielens, A. G. G. M., Tokunaga, A. T., Geballe, T. R., & Baas, F. 1991, ApJ, 381, 181

Van der Tak, F.F.S., Boonman, A.M.S., Braakman, R., van Dishoeck, E.F. 2003, A&A, in press (astro-ph/0309138)

Van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J., & Blake, G. A. 2000a, ApJ, 537, 283

Van der Tak, F. F. S., van Dishoeck, E. F., & Caselli, P. 2000b, A&A, 361, 327

Van Dishoeck, E. F. 1998, Faraday Discussion, 109, 31

Van Dishoeck, E. F. & van der Tak, F. F. S. 2000, in IAU Symp. 197: Astrochemistry: From Molecular Clouds to Planetary Systems, eds. Y. C. Minh & E. F. van Dishoeck, 97

Van der Tak, F. F. S. & van Dishoeck, E. F. 2000, A&A, 358, L79

Walsh, A. J., Bertoldi, F., Burton, M. G., & Nikola, T. 2001, MNRAS, 326, 36

Watanabe, N. & Kouchi, A. 2002, ApJ, 571, L173

Wilner, D., De Pree, C., Welch, J., & Goss, M. 2001, ApJ, 550, L81

Wyrowski, F., Schilke, P., Walmsley, M., & Menten, K. 1999, ApJ, 514, L43


text:
# A Snapshot Survey for Gravitational Lenses Among 𝑧≥4.0 Quasars: I. The 𝑧>5.7 Sample

## 1 Introduction

Before 1995, only 24 quasars with redshifts larger than four had been published; most had been found in multicolor or grism surveys for high-redshift quasars (e.g., Warren et al. 1987; Schmidt, Schneider, & Gunn 1987). In the past eight years the number of known $`z>4`$ quasars has increased to more than 400, largely due to the discoveries from the Sloan Digital Sky Survey (SDSS; York et al. 2000) and the Digitized Palomar Sky Survey (Kennefick, Djorgovski, & Meylan 1996). The SDSS to date has published the discovery of over 200 quasars with redshifts greater than 3.6 (Fan et al. 2001b; Anderson et al. 2001; Schneider et al. 2001; Schneider et al. 2003). In addition, Fan et al. (2000), Fan et al. (2001a), and Fan et al. (2003) have discovered the seven highest-redshift quasars known, all with $`z>5.7`$ and four of which are the subject of this paper. This large sample provides us with new opportunities to study high-redshift quasars; it has already been used to determine the luminosity function of quasars and its evolution with redshift from $`z=3.6`$ to $`z=6.0`$ (Fan et al. 2001c, 2003). The luminosity function is well-fit with a power-law, while the comoving number density of quasars drops exponentially with redshift over this range. One of the goals of this paper is to determine the extent to which the observed luminosity function might be biased by gravitational lensing.

Understanding the (intrinsic) quasar luminosity function and evolution thereof at high redshifts is important for several reasons. We observe that supermassive, quiescent black holes are ubiquitous at low redshift (Magorrian et al. 1998); comparison of the present-day mass function with the high-redshift luminosity function can constrain models for the duty cycle of black holes and their feeding mechanisms. There is increasing evidence for a population of optically faint AGN at high redshift, which appear in deep X-ray images with the Chandra satellite (Mushotzky et al. 2000; Brandt et al. 2001; Alexander et al. 2001); understanding their nature requires that we understand the nature of the optically luminous objects as well. Finally, the optical (rest-frame UV) luminosity function puts important constraints on the ultraviolet ionizing background (Fan et al. 2001c), which is especially important now that we appear to be observing the first hints that the universe became completely reionized around redshift $`z=6`$ (e.g., Becker et al. 2001; Djorgovski et al. 2001; White et al. 2003).

The continuum luminosities of the $`z>5.7`$ quasars are all very high, with absolute magnitudes at rest frame 1450 Å between $`26.8`$ and $`28`$ for the WMAP cosmology that we use throughout this paper (a flat $`\mathrm{\Lambda }`$ CDM model with $`\mathrm{\Omega }_m=0.29`$ and $`H_0=72\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$; Spergel et al. 2003). This yields a lower limit to the masses of their black holes of order $`4\times 10^9\mathrm{M}_{}`$ by the Eddington argument (Fan et al. 2000). It is a challenge to explain in standard cosmologies how such massive black holes could form at a cosmic epoch less than a billion years after the Big Bang (e.g., Efstathiou & Rees 1988; Turner 1991; Haiman & Loeb 2001).

All these results assume that the apparent magnitudes that we measure for these quasars are intrinsic to the quasar themselves, and represent isotropic fluxes. However, gravitational lensing can systematically brighten the images of quasars by a significant amount. This is a particularly important effect for flux-limited surveys like SDSS, which are sensitive to the most luminous quasars at any given redshift (e.g., Maoz & Rix 1993). The resulting magnification bias (e.g., Pei 1995, and references therein) can systematically change the quasar luminosity function, causing an overestimation of the black hole masses powering these objects.

The expected fraction of multiply imaged quasars at a given redshift depends both on the cosmological model (which sets the distribution and profiles of halos as a function of redshift, as well as the dependence of cosmological pathlength on redshift), and on the luminosity function of quasars, which determines the effect of magnification bias. In the WMAP cosmology (using a halo distribution taken from the large N-body simulations of Jenkins et al. 2001), the fraction of random lines of sight at $`z=4`$ that produce multiple images at all splitting angles (see the discussion of lensing models in § 3 below) is of order 0.2%; this fraction rises to 0.4% at $`z=6`$. However, magnification bias greatly increases the lensing rate in a flux-limited sample. At $`z=4`$, using the canonical cosmology, and a luminosity function $`\mathrm{\Phi }(L)L^{2.5}`$ (Fan et al. 2001b) with a break at $`M_{1450}^{}=24.0`$ (see § 3 below for further discussion), we expect about 1% of quasars in our sample to be lensed (specifically, multiply imaged with magnification of at least two). This fraction increases to 8% if we assume a steeper slope of $`\beta =3.5`$, and to $`2\%`$ if we keep the slope of $`\beta =2.5`$, but assume a break in $`\mathrm{\Phi }`$ at a 10 times lower luminosity ($`M_{1450}^{}=21.5`$). At $`z=6`$, magnification bias can result in lensing probabilities approaching 100% if the luminosity function is particularly steep (Comerford, Haiman, & Schaye 2002; Wyithe & Loeb 2002a, 2002b). However, the detailed shape of $`\mathrm{\Phi }`$ is not well-measured, especially at higher redshifts and luminosities. Such quasars represent the highest peaks in the density field at high redshift, where the mass function of massive dark halos — in which bright quasars likely reside — is expected to be very steep (Fan et al. 2001c).

The median expected splitting of gravitationally lensed quasars is predicted to be somewhat less than $`1\mathrm{}`$ (Turner, Ostriker, & Gott 1984), with flux ratios as large as 20:1; this is in broad agreement with the observed splitting separation found by the Cosmic Lens All Sky Survey (CLASS; Browne et al. 2003), but see the discussion in §3 below on CLASS implications for the details of the lens profiles. The SDSS images themselves have PSF widths of order $`1\mathrm{}1\stackrel{}{\mathrm{.}}5`$ (Abazajian et al. 2003), making them inadequate to look for all but the very rare wide-separation ($`\mathrm{\Delta }\theta >1\stackrel{}{\mathrm{.}}5`$) lenses (Pindor et al. 2003). We are thus carrying out a snapshot imaging survey of 250 SDSS-discovered quasars with $`z>4.0`$ using the Advanced Camera for Surveys (ACS; Ford et al. 2003) on the Hubble Space Telescope (HST), which can easily resolve pairs separated by as little as $`0\stackrel{}{\mathrm{.}}1`$ with flux ratios as large as 10:1. Our observations should therefore be sensitive to essentially all expected lenses (Turner et al. 1984; Hinshaw & Krauss 1987). Fan et al. (2003) presented some preliminary results based on some of these (and other, non-HST) data; here we give a full analysis of the HST observations for all four of the $`z>5.7`$ quasars in our sample.

Throughout this paper, we adopt the WMAP cosmology mentioned above, with an rms mass fluctuation within a sphere of radius $`8h^1`$ Mpc of $`\sigma _8=0.9`$, and power–law index $`n=0.99`$ for the power spectrum of density fluctuations (Spergel et al. 2003). We also adopt the cosmological transfer function from Eisenstein & Hu (1999). Conversions between $`M_B`$ and $`M_{1450}`$ assume $`M_B=M_{1450}0.48`$ (Schmidt, Schneider, & Gunn 1995) with spectral index $`\alpha _\nu =0.5(f_\nu \nu ^{\alpha _\nu })`$.

## 2 The Data

### 2.1 Observations

The full sample of $`z>4.0`$ quasars included in our HST snapshot program will be presented in future papers in this series. In this paper, we present results of imaging the four quasars with $`z>5.7`$ discovered by Fan et al. (2000, 2001a); see Table 1. These four objects were selected from the SDSS imaging database based on their extremely red $`iz`$ color, absence of detection in $`u,g,`$ and $`r`$, and relatively blue $`zJ`$ color (which distinguishes the $`z>5.7`$ quasars from brown dwarfs).

The SDSS images of these quasars have an image scale of $`0\stackrel{}{\mathrm{.}}396\mathrm{pixel}^1`$ and $`z`$-band seeing FWHM between $`1\stackrel{}{\mathrm{.}}4`$ and $`1\stackrel{}{\mathrm{.}}8`$. In order to test the ability of the SDSS photometric pipeline to identify quasar pairs, we created simulated SDSS images of pairs of point sources (see Pindor et al. 2003 for details) in these observing conditions and at the appropriate signal-to-noise ratio. We find that the SDSS star-galaxy separator would have identified pairs of point sources as resolved objects for image separations greater than $`1\mathrm{}`$ and for flux ratios less than $``$ 5:1. Thus the SDSS data are suitable for exploring only a fraction of the parameter space of lensing (in terms of separations and flux ratios) that is of interest; higher resolution HST observations were needed to probe the rest of that parameter space. We did not restrict our sample to point sources since the SDSS star-galaxy separation algorithm (see Scranton et al. 2002) breaks down at these faint magnitudes. The only object of the four that was flagged as an extended source was SDSSp J104433.04-012502.2.

We used the High Resolution Camera (HRC) on ACS to image the four $`z>5.7`$ quasars presented herein. For these four objects, the Ly$`\alpha `$ forest lies in the SDSS $`i`$-band; we therefore observed them in the SDSS $`z`$-band, F850LP, despite its lower sensitivity. The Wide-Field Camera on ACS has higher sensitivity in $`i`$ and $`z`$ than does the HRC, but has has substantially higher overhead and is mildly undersampled. The exposure times were 1200 seconds for each object, 600 seconds in each of two exposures to help in cosmic ray rejection. We did not dither the images, and our cosmic ray rejection is far from perfect (§ 2.2); however, our analysis indicates that imperfect cosmic ray removal is not critical for our search for gravitational lenses.

### 2.2 Data Processing

The raw images were calibrated by the CALACS package in IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. as part of on-the-fly-reprocessing (OTFR) at the time of download. For the most part, we find these initial reductions to be sufficient for our purposes. The images that we present are the “cosmic ray rejected” (CRJ) images that are output by the OTFR algorithms at STScI. These files were re-requested a number of weeks after the observations were completed to ensure that we are using the latest “reference” files and that the OTFR reprocessing is current. The CRJ files have all been reduced in the standard manner, including having been overscan-, bias- and dark-corrected, flat-fielded and photometrically calibrated, in addition to having bad pixels masked and cosmic rays removed (see the ACS manual222http://www.stsci.edu/hst/acs/documents/handbooks/cycle12/cover.html for more details). The plate scale on the HRC of the ACS is not the same in the $`x`$ and $`y`$ directions and the sky axes are not perpendicular on the CCD; these CRJ images have not been corrected for this distortion. A $`5\mathrm{}\times 5\mathrm{}`$ cutout of each of the ACS images is given in Figures 1-4.

With the exception of cosmic ray rejection and hot pixel identification, the images output by the OTFR algorithm are acceptable for our purposes. The CALACS rejects pixels that are significantly different in the two images. This process successfully rejects the majority of the cosmic rays in our images. However, many of the cosmic rays occupy multiple pixels, and we have only two exposures, meaning that some cosmic rays remain. In addition, hot pixels can remain in the corrected images if they have not been properly identified in the masks that are constructed for each period of observations between the annealing procedure that helps to repair hot pixels. Fortunately, the morphologies of the residual cosmic rays and hot pixels are considerably different from that of the ACS PSF (the PSF is spread out over many tens of pixels), so spurious sources due to cosmic rays and most other defects are easy to recognize.

### 2.3 Subtracting the Point Source

The point spread function is narrow enough that any lens with a separation greater than $`0\stackrel{}{\mathrm{.}}2`$ will be obvious by visual inspection (see § 2.5). To search for faint secondary images at separations smaller than this requires that we fit and subtract off a model for the point spread function of each image. We have done this using the v6.0 (pre-launch of ACS) version of the Tiny Tim software (Krist 1995)333http://www.stsci.edu/software/tinytim/, which produces a model PSF for the instruments on HST given the object’s SED and position in the focal plane, the filter curve, and knowledge of the optics of the instrument. Such model PSFs are necessary because it is not always possible to do PSF subtraction using an empirical template (another point source in the field-of-view). That is, given the positional dependence of the PSF across the field-of-view, the Tiny Tim models are often better than using an empirical template. For our ACS/HRC images, the model correctly takes into account the different scales and non-orthogonality of the $`x`$ and $`y`$ axes.

We fit the PSF model to each CRJ image, allowing the location on the CCD and normalization to vary, and using sinc interpolation when shifting the model PSF by fractional pixels. Inputs to Tiny Tim were the camera used (ACS/HRC), the location on the CCD of each quasar to the nearest pixel (as measured by imexam in IRAF) and a Keck spectrum (Becker et al. 2001) of each quasar (for the wavelength dependence of the PSF). We minimized the sum of the square of residuals (since the formal errors per pixel over the image are close to constant). The results are shown in Figures 14. Each of the four images show the familiar first Airy ring; on a very hard stretch, the second Airy ring is faintly visible. Note that HST PSFs are not symmetric as a result of the three primary mirror supports, the signature of which can be seen in the PSF subtracted images. Given inevitably somewhat imperfect PSF subtraction and the expected asymmetry of the residual images, in none of the four cases is there an obvious secondary image; none of them are obviously gravitationally lensed.

### 2.4 Notes on Individual Objects

#### 2.4.1 SDSSp J104433.04-012502.2; $`z=5.74`$

The ACS/HRC image of this quasar is shown in the left panel of Figure 1; the right panel shows the residuals after PSF subtraction. The SDSS $`z`$-band magnitude of this quasar is 19.3; the HST $`z`$-band (ST) magnitude is 19.91, see also Table 1. The offset between the SDSS and ST $`z`$-band magnitude is systematic (see the other objects below) with the ST magnitudes always being fainter. This is probably due to the differences in the filter curves (the SDSS curve having less red response due to atmospheric features) and the differences between a Vega-based (ST) and an AB-based (SDSS) photometric system given that Vega has an AB $`z`$-band magnitude of 0.57 (Fukugita et al. 1996). We find no offset between the SDSS and ST $`i`$-band magnitudes (in the HST images of our lower redshift quasars), where the filter curves appear to be much more similar.

Fan et al. (2000) describe a $`K^{}`$ image of this object taken with the Keck 10m telescope in $`0\stackrel{}{\mathrm{.}}375`$ seeing, which was consistent with a point source. That is consistent with the HST image; there is no apparent splitting in this object.

Shioya et al. (2002) found a faint galaxy ($`m_B25`$) $`1\stackrel{}{\mathrm{.}}9`$ south-west of this quasar. Our 20-minute exposure shows this object faintly after rebinning. Shioya et al. (2002) suggest that this galaxy could amplify the flux of the quasar by up to a factor of two. The lack of a secondary image of the quasar restricts any possible magnification to be no larger than this — in the absence of microlensing and external shear (Wyithe & Loeb 2002a).

#### 2.4.2 SDSSp J083643.85+005453.3; $`z=5.82`$

The ACS/HRC image of this quasar and residuals after PSF subtraction are shown in Figure 2. The SDSS $`z`$-band magnitude of this quasar is 18.79; the HST $`z`$-band (ST) magnitude is 19.49, see also Table 1. One other point source and two other extended sources are in the ACS/HRC field of this quasar, but all are more than $`5\mathrm{}`$ away (out of the field-of-view of Fig. 2), and are therefore unlikely to lens (or significantly magnify) the quasar. The black dots in the field-of-view of Figure 2 are residual cosmic rays that were not properly cleaned by the OTFR reductions.

#### 2.4.3 SDSSp J130608.26+035626.3; $`z=5.99`$

The ACS/HRC image of this quasar is shown in Figure 3. The SDSS $`z`$-band magnitude of this quasar is 19.46; the HST $`z`$-band (ST) magnitude is 20.23, see also Table 1. There is some evidence for a few very faint extended sources in the ACS/HRC field (off the image in Fig. 3), but as with Figure 2, all of the black dots in Figure 3 are residual cosmic rays. Yamada et al. (2003) have taken deep images to $`i25.5`$ in moderate seeing with Subaru in this field, and have found no foreground galaxies close enough to the line of sight to cause gravitational lensing.

#### 2.4.4 SDSSp J103027.10+052455.0; $`z=6.30`$

The ACS/HRC image of this quasar is shown in Figure 4. The SDSS $`z`$-band magnitude of this quasar is 20.02; the HST $`z`$-band (ST) magnitude is 20.54, see also Table 1. There is one other relatively bright point source in the field, nearly due west of the target quasar, but it is at the edge of the field (off the image in Fig. 4) and is not of interest in terms of lensing. As in the previous object, a deep image by Yamada et al. (2003) shows no obvious lensing galaxies in the foreground. However, Petric et al. (2003) have taken deep VLA images at 20 cm in this field, and found a statistically significant overdensity of four radio sources within $`1^{}`$ of the quasar, which may indicate the possibility of a magnifying cluster along the line of sight.

### 2.5 Looking for Multiple Images

Our ACS images have an image scale of $`0\stackrel{}{\mathrm{.}}025\mathrm{pixel}^1`$, and we can visually rule out any multiple imaging at scales larger than about $`0\stackrel{}{\mathrm{.}}2`$. To make this more quantitative, we have carried out simulations of quasar pairs at various separations and flux ratios, by adding a scaled-down and shifted version of our quasar images to itself. This naive procedure produces a slightly incorrect sky level, but the inaccuracy is negligible for our present purposes. Figure 5 shows simulated observations of SDSS J0836+0054 with an additional point source 2.5 and 5.0 mag fainter than the quasar and separated from the quasar by $`0\stackrel{}{\mathrm{.}}3`$. The secondary image is obvious in both cases even without PSF subtraction. Figure 6 shows the result of a $`0\stackrel{}{\mathrm{.}}1`$ offset. With a flux ratio of 10:1, the secondary object is visible as an enhancement in the first Airy ring of the primary object; it becomes clear upon PSF subtraction (Figure 7). With a flux ratio of 100:1, we cannot discern a pair with $`0\stackrel{}{\mathrm{.}}1`$ separation even after subtraction of the PSF; however, even if we were able to do so, it would not significantly improve the constraints that we derive below since most split images will have larger image separations.

## 3 Constraining the Slope of the Quasar Luminosity Function

The shape of the high-redshift quasar luminosity function is not tightly constrained because of the limited size and luminosity range of the available statistical samples. Schmidt et al. (1995), using a set of 90 quasars with $`2.7<z<4.8`$, found a power-law luminosity function slope ($`\beta `$) of $`2`$. Fan et al. (2001c) measured $`\beta 2.5\pm 0.3`$ from 39 quasars in the redshift interval 3.6 to 5.0. These slopes apply to the high-luminosity ($`M26`$) part of the luminosity function. At these high redshifts, direct observations currently yield only poor constraints on the slope of the quasar luminosity function (LF) at $`z5`$, and essentially no constraints at $`z>\mathrm{\hspace{0.33em}6}`$, where the luminosities of the known quasars all lie within $`0.5`$ magnitudes of the detection threshold of SDSS. From the number counts of $`z6`$ quasars, Fan et al. (2001c) place a 3-$`\sigma `$ lower limit on the slope of $`\beta =4.6`$. Gravitational lensing tends to flatten the intrinsic luminosity function, producing an apparent slope of $`\beta =3`$ in the limit of an arbitrarily steep intrinsic LF down to a fixed low–luminosity cutoff (e.g., Schneider, Ehlers, & Falco 1992). As a result, the relatively shallow apparent slope at $`z5`$ does not necessarily reflect the slope of the underlying LF, which could be much steeper, if most quasars were lensed and magnified.

One of the primary goals of this work is to constrain the slope of the quasar LF at high redshift, using the lack of evidence for strong lensing among the four quasars in the HST images (Comerford et al. 2002). In line with previous work, we describe the intrinsic (not necessarily observed) quasar luminosity function as a broken power law (e.g., Boyle, Shanks, & Peterson 1988; Pei 1995):

$$\mathrm{\Phi }_{\mathrm{int}}(L)=\frac{\mathrm{\Phi }_{}/L_{}}{(L/L_{})^{\beta _l}+(L/L_{})^{\beta _h}}.$$
(1)
The LF is described by four parameters: the normalization $`\mathrm{\Phi }_{}`$, the faint-end slope $`\beta _l`$, the bright-end slope $`\beta _h`$, and the characteristic luminosity $`L_{}`$ at which the LF steepens. The lensing probability is very sensitive to the last two parameters, $`\beta _h`$ and $`L_{}`$, while it is only weakly dependent on the faint–end slope below the break, and it is strictly independent of the normalization $`\mathrm{\Phi }_{}`$. Here we apply the lensing model from Comerford et al. (2002), in which lenses are associated with dark matter halos, to compute the total lensing probability, including the effect of magnification bias. In this model, the abundance of halos as a function of potential well depth is adopted from the simulations of Jenkins et al. (2001). While it would be simplest to assume that the halos have a common density profile, such as the singular isothermal sphere (SIS), recent work has suggested that the typical density profile of massive lenses is shallower than that of the SIS (Porciani & Madau 2000; Kochanek & White 2001; Li & Ostriker 2002; Huterer & Ma 2003). The observed dearth of large–separation lenses cannot be explained if all lenses have the same density profile. A prescription that fits the observed steep splitting angle distribution is to assume that all halos below $`M10^{13}\mathrm{M}_{}`$ have SIS profiles (adopting a standard conversion between circular velocity and halo mass), while all halos above this mass follow the dark matter density profile suggested by Navarro, Frenk, & White (1997, hereafter NFW). NFW profiles are much less efficient lenses than are SIS profiles. This prescription is physically motivated in Kochanek & White (2001) based on the efficiency of gas cooling, and we adopt it in this paper. In our calculations, this is essentially equivalent to ignoring all lenses above a halo mass of $`10^{13}\mathrm{M}_{}`$, as the massive halos do not contribute to lensing at small separations.

As an example, we first fix the characteristic luminosity (“knee”) of the quasar LF at a fiducial value of $`M_{1450}^{}=24.0`$, and derive constraints on the bright-end slope, $`\beta _h`$. The choice for this value is somewhat arbitrary: the highest redshift where the quasar LF is measured to faint enough magnitudes to define a reliable break is at $`z2.15`$ (Boyle et al. 2000), where $`M_{1450}^{}25.6(M_B^{}26.1)`$. In pure luminosity evolution models, the break at $`z5`$ would then be at around $`M_{1450}^{}=24.0`$ (e.g. Madau, Haardt, & Rees 1999; Wyithe & Loeb 2002a). Having fixed the value of the break, we next vary the slope $`\beta _h`$ of the LF, and require that the expected probability for no lensing among the four highest redshift quasars is P=0.0026 ($`=3\sigma `$) or greater. Assuming that we can rule out a second image with a flux ratio of 100 down to a splitting angle of $`0\stackrel{}{\mathrm{.}}3`$, we find the constraint on the slope $`\beta _h>4.63`$. We find that resolving high–contrast lensing configurations, with a flux ratio $`>10`$, or resolving angular separations smaller than $`0\stackrel{}{\mathrm{.}}3`$ does not significantly improve this constraint, since in any model, the fraction of lenses with separations less than $`0\stackrel{}{\mathrm{.}}3`$ is small. For example, ruling out flux ratios of 10 down to $`0\stackrel{}{\mathrm{.}}1`$ would only improve the constraint to $`\beta _h>4.58`$.

It is more informative to illustrate joint constraints on the parameters $`(M_{1450}^{},\beta _h)`$. In Figure 8, we show contours for the likelihood at $`P=0.32`$, 0.05, 0.01, and 0.001 that none of the four $`z>5.7`$ quasars are lensed (again, to a separation of 0.3 arcseconds, and to flux ratios of 100). The figure shows, as in Fan et al. (2003), that the derived limit on the slope $`\beta _h`$ depends strongly on the assumed break. Overall, the constraint shown in this figure is similar to that derived from the sample of seven high–redshift quasars in Fan et al. (2003), although the present constraints are somewhat weaker, due to the smaller number of sources. The recent results of Wyithe (2003) are consistent with ours under the assumption of no lensing. As shown in Wyithe (2003), a stronger constraint is available if two of the known $`z>6`$ quasars are indeed magnified, as suggested by Shioya et al. (2002) and White et al. (2003).

## 4 Discussion

The lack of lenses among the four high–redshift quasars allows us to place constraints on the quasar LF at high redshift, shown in Figure 8. It is useful to keep in mind that only models with extreme magnification bias, predicting significant lensing probabilities, can be constrained by the lack of lensing events. For example, for the range of LF slopes considered by Wyithe & Loeb (2002a), $`\beta _h=2.583.43`$, the lensing fraction for the $`z6`$ SDSS sample would be 7-30%. The chance of finding no lenses among the four quasars for their steepest slope of 3.43 is $`(10.3)^4=0.24`$, which translates to a significance of only $`1\sigma `$ for ruling out this slope (consistent with our results shown in Figure 8 at the break of $`M_{1450}=24.0`$).

Accordingly, as Figure 8 shows, the constraints obtained here still allow very steep LF slopes, especially if we assume that the break in the LF lies at relatively bright magnitudes. Nevertheless, these constraints are fairly direct, and are among the most stringent limits we currently have on the shape of the high-redshift quasar LF. We expect that the full $`z>4`$ sample will further tighten the constraints (or yield a measurement should a lens be found). Table 2 shows that the lensing probability (for a slope of $`\beta =3.5`$ and break $`M_{1450}=24.0`$) decreases from 8.2% at $`z=6`$ to 1.2% at $`z=4`$. This implies that seeing no lensing with a single $`z=6`$ quasar provides a constraint that is equivalent to about seven $`z=4`$ quasars, from the equation $`(10.012)^N=(10.082)`$. The lensing probability is a steep function of $`\beta _h`$ and $`M^{}`$, but since our full sample will include $`250/460`$ times more sources than utilized here, we expect that we will be able to significantly tighten the constraints in the $`(\beta _h,M^{})`$ plane with this sample.

Given our assumed shape of the quasar LF, the constraints obtained in the $`(\beta _h,M^{})`$ plane translate directly to the global ionizing emissivity of the $`z6`$ quasar population, $`\nu ϵ_\nu =L\mathrm{\Phi }(L)𝑑L`$. The emissivity is of great interest, since a related quantity — the mean background flux — can be constrained independently from the level of resonant neutral hydrogen absorption detected in high–resolution spectra of the SDSS quasars. Specifically, recent upper limits on this background were obtained from the Gunn–Peterson (1965) trough of the $`z=6.30`$ quasar (e.g., McDonald & Miralda-Escudé 2001; Cen & McDonald 2002; Lidz et al. 2002; Fan et al. 2002), requiring that $`\mathrm{\Gamma }<10^{13}`$ ionizations per hydrogen atom per second at $`z=6`$. The conversion to the ionizing emissivity $`\nu ϵ_\nu `$ is

$$\mathrm{\Gamma }=\frac{\pi \sigma _H}{h_\mathrm{P}}\lambda _{\mathrm{mfp}}\frac{(1+z)^3}{4\pi }\frac{\nu ϵ_\nu }{\nu _\mathrm{H}}\mathrm{s}^1,$$
(2)
where we have assumed $`ϵ_\nu \nu ^1`$, $`\sigma _H=6.3\times 10^{18}\mathrm{cm}^2`$ is the hydrogen ionization cross section at the Lyman limit (and we assumed $`\sigma _H(\nu )\nu ^3`$ above the threshold), $`h_\mathrm{P}`$ is the Planck constant, $`\nu _\mathrm{H}=3.29\times 10^{15}`$Hz is the frequency of a Lyman limit photon, and $`\lambda _{\mathrm{mfp}}`$ is the mean free path of ionizing photons (see the Appendix in Comerford et al. 2002). The last quantity unfortunately has a very large uncertainty. Nevertheless, for a crude estimate, we here assume that the effective mean free path corresponds to a redshift interval of $`\mathrm{\Delta }z=0.17`$ at $`z=3`$ (see Haardt & Madau 1996; Steidel, Pettini, & Adelberger 2001), and scales with redshift as $`\lambda _{\mathrm{mfp}}(1+z)^6`$ towards higher redshift (Cen & McDonald 2002). In this case, the model with $`M_{1450}^{}=24.0`$ and $`\beta _h=4.6`$ corresponds to $`\mathrm{\Gamma }10^{13}\mathrm{sec}^1`$. In other words, the constraint we derive here from the absence of lensing is approximately at the same level as that implied by the upper limit on the ionizing background. Note that the constraints we obtain apply to the total emissivity of quasars alone, but the upper limits on the background from the Gunn-Peterson trough constrains the total emissivity (including galaxies). Future improvements of the lensing constraint can yield tighter limits on the ionizing emissivity of the quasar population than hitherto possible — subject to the uncertainties in (1) the mean evolution of the mean free path with redshift, (2) the faint–end slope $`\beta _l`$, which does affect the total emissivity $`\nu ϵ_\nu `$, and (3) the escape fraction of ionizing photons at $`z6`$.

The fact that none of the four high–redshift quasars appears to be strongly magnified also implies that the masses inferred for the central black holes from the Eddington limit argument are real. If a quasar is strongly beamed towards us, the bolometric luminosity would be overestimated, but Willott, McLure, & Jarvis (2003), following Haiman & Cen (2002) argued that emission line equivalent widths being normal mitigates the beaming hypothesis. The inferred black hole masses in these four quasars are very high, $`24\times 10^9\mathrm{M}_{}`$; how black holes can grow so large so early in the Universe requires explanation (Efstathiou & Rees 1988; Turner 1991). Haiman & Loeb (2001) considered the presence of the $`4\times 10^9\mathrm{M}_{}`$ black hole for the $`z=5.75`$ quasar in simple models for black hole growth. They found that it could be marginally accommodated in simple models of black hole growth by mergers and accretion in an LCDM universe. These considerations depend strongly on the redshift, and on the typical radiative efficiency of black holes. The highest redshift quasars therefore represent a bigger challenge to models, especially if their radiative efficiency during their growth was as high as 20%, the typical value inferred recently for the brightest quasars at lower redshifts (Yu & Tremaine 2002). This is twice the fiducial value of $`10\%`$ adopted in Haiman & Loeb (2001) and most other previous work on black hole growth, doubling the expected growth timescale.

Finally, the lack of strong lensing appears consistent with the large apparent sizes of ionized regions around the $`z6`$ quasars. Haiman & Cen (2002) argue that the $`z=6.30`$ source cannot be magnified by more than a factor of a few, since an intrinsically faint source would not be able to ionize a neutral IGM out to the observed distance away from the source of $`30`$ (comoving) Mpc — the spectrum shows transmission at wavelengths corresponding to this distance. Similar conclusions were drawn for the $`z=6.42`$ quasar (Fan et al. 2003) recently by White et al. (2003).

## 5 Conclusions and Future Work

We have obtained high resolution HST images of the four highest redshift quasars known prior to 2003 to look for the signature of gravitational lensing. We have found no evidence of multiple images, significantly limiting the amount by which these quasars can be magnified by foreground mass concentrations (in the absence of microlensing). Thus the masses of these quasars as derived from the Eddington argument cannot be significantly inflated by magnification.

The lack of any strong lenses puts a $`3\sigma `$ constraint on the bright end slope of the $`z>5`$ luminosity function of $`\beta _h>4.63`$. This limit on the slope of the luminosity function is also consistent with previous limits on the ionizing emissivity of $`z6`$ quasars of $`\mathrm{\Gamma }10^{13}\mathrm{sec}^1`$ ionizations per hydrogen atom per second based on the Gunn-Peterson trough in the highest redshift quasar studied herein.

The next paper in the series will present an analysis of a much larger sample of $`z>4`$ quasar observations with HST and will discuss the constraints that we can place upon the quasar luminosity function from that larger (but lower redshift) sample. Furthermore, there are many more $`z>4`$ (and even a few $`z>5.7`$) quasars known now than were known at the time that these HST observations were approved. Higher resolution images should be taken of as many of these as possible; much of the parameter space can even be covered by ground-based observations since HST resolution splittings are only needed to discern the smallest expected splittings.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. This work was supported in part by HST grant HST-GO-09472.01-A (G. T. R. and M. A. S.), National Science Foundation grants AST-9900703 (D. P. S.), AST-0071091 (M. A. S.), AST-0307384 (X. F.), AST-0307291 (Z. H.) and AST-0307200 (Z. H.). X. F. and D. J. E. acknowledge Alfred P. Sloan Research Fellowships.