text:
# The evolution of faint AGN between z ≃ 1 and z ≃ 5 from the COMBO-17 survey

## 1 Introduction

The luminosity function of quasi-stellar objects (QSOs) and its evolution with redshift provides one of the most important tools for the cosmic demography of active galactic nuclei (AGN). It constrains physical models for QSOs, particularly those for the growth of supermassive black holes in galaxies within the context of hierachical collapse of structure in the universe (Haehnelt & Rees, 1993; Haiman & Loeb, 1998; Kauffmann & Haehnelt, 2000). It is also relevant for understanding the extragalactic UV background (Meiksin & Madau, 1993; Boyle & Terlevich, 1998).

Previous studies of the QSO luminosity function (QLF) have established a strong rise in the activity with look back time from the local universe to redshifts $`z2`$ (Boyle et al., 1988). The main debate at these low to intermediate redshifts is now about the question whether the shape of the QLF changes with redshift (Boyle et al., 2000; Wisotzki, 2000; Miyaji et al., 2000; Cowie et al., 2003). If so, this would mean that low-luminosity AGN evolve differently from high-luminosity QSOs, possibly calling for a substantial revision of our understanding of the cosmic duty cycle of nuclear activity in galaxies.

At higher redshift $`z>2`$, the obvious expectation that the rise of the QSO activity with redshift must turn over at some point was satisfied by observations at $`z>3`$ (Warren et al., 1994; Schmidt et al., 1995, hereafter WHO and SSG). However, optical surveys in this redshift regime were often plagued by selection effects and small number statistics. Moreover, studies of X-ray-selected (Miyaji et al., 2000) and radio-selected QSOs (Jarvis & Rawlings, 2000) did not strongly support claims of a very steep drop in nuclear activity for $`z>3`$, keeping the issue unresolved whether the detected turnover was physical reality or not. A rather robust observation of the space density decline beyond $`z>3.6`$ has been established recently, albeit only for very luminous QSOs, by the Sloan Digital Sky Survey (SDSS, Fan et al., 2001).

Most optically selected QSO samples are expected to be more or less complete at either $`z2.2`$ or $`z3.6`$, where QSOs show conspicuous colours in broad-band searches, while within this redshift band, confusion arises with stars and compact low-redshift galaxies. Any practical approach of following up QSO samples in this intermediate redshift range is forced to avoid strong contamination and observed mostly extreme objects, resulting in a high degree of incompleteness that is difficult to account for. Thus, while the fact that cosmic QSO activity shows a maximum around $`z2`$–3 is not doubted as such, this maximum has very rarely been detected in a single survey.

This selection issue was one of the original motivations for the COMBO-17 survey. This survey is based on photometric classification and redshift estimation using a set of 17 filters, of which 12 are medium-band filters with $`10,000`$ km/sec FWHM, well matched to optimum detection of QSO emission lines and spread across the range of wavelengths observable by modern CCD detectors.

Another major driver for this project was to probe fainter regimes of the AGN luminosity function, a task which has always beed limited by the overwhelming need for spectroscopic follow-up telescope time. Using medium-band spectrophotometry, a reasonably complete and clean AGN sample could be obtained across a wide range of redshifts down to $`R24`$. Based on such a deep sample, we were able to derive luminosity functions down to $`M_B23`$ at redshifts above 3, entering the domain of Seyfert galaxies even at high redshift.

In this paper, we present a new determination of the luminosity function of faint optically selected QSOs, aimed at a broad redshift range around the elusive turnover. The analysis is based on a sample of 192 objects between $`z=1.2`$ and $`z=4.8`$, all selected from the COMBO-17 survey. The survey is briefly described in Sect. 2, followed by a discussion of technical aspects such as completeness and redshift quality in Sect. 3. In Sections 4 through 7 we present and discuss the results.

## 2 COMBO-17 observations

The COMBO-17 project (“Classifying Objects by Medium-Band Observations in 17 Filters”) was designed to provide a sample of $``$50,000 galaxies and several hundred AGN with precise photometric redshifts and spectral energy distributions (SEDs). As shown below, the filter set provides a redshift accuracy of $`\sigma _z0.03`$ for quasars (and similarly for galaxies; cf. Wolf et al., 2003), smoothing the unknown true redshift distribution of the sample only very mildly and certainly allowing the derivation of luminosity functions.

The first step of the COMBO-17 data analysis was to convert the photometric observations into a very low resolution ‘fuzzy spectrum’, allowing for simultaneous spectral classification into stars, galaxies and QSOs, as well as for accurate redshift and SED estimation for the latter two categories. The full survey catalogue will contain about 50,000 objects with classifications and redshifts covering a solid angle of 1.0 $`\text{deg}^2`$. This fuzzy spectroscopy consciously compromises on redshift accuracy in order to obtain large samples of quasars with a reasonable observational effort.

While the photometric redshift technique has already been applied to galaxy samples about 40 years ago (Baum, 1962; Butchins, 1983), we have modified and improved the approach by increasing the number of filters and narrowing their bandwidth to obtain better spectral resolution and more spectral bins. This way, COMBO-17 provides identifications and reasonably accurate redshifts not only for galaxies but also for quasars, a novelty pioneered in CADIS (Wolf et al., 1999) and more recently applied in the SDSS (Richards et al., 2001) and to observations with superconducting tunnel junctions (deBruijne et al., 2002).

All observations presented here were obtained with the Wide Field Imager (WFI, Baade et al., 1998) at the MPG/ESO 2.2-m telescope on La Silla, Chile. The WFI provides a field of view of $`34\mathrm{}\times 33\mathrm{}`$ on a CCD mosaic consisting of eight 2k $`\times `$ 4k CCDs with a scale of $`0\stackrel{}{.}238`$/pixel. The total observing time for COMBO-17 was about $``$175 ksec per field, including a $``$20 ksec exposure in the R-band with seeing below 0$`\stackrel{}{.}`$8.

Observations and data analysis have been completed on three fields including the Chandra Deep Field South (CDFS), covering an area of 0.78 $`\text{deg}^2`$ (Tab. 2). The deep $`R`$-band images have 5-$`\sigma `$ point source limits of $`R26`$ and provide the highest signal-to-noise ratio for object detection and position measurement among all data in the survey, except for L-stars and quasars at $`z>5`$, which can be extremely faint in R.

Using SExtractor (Bertin & Arnouts, 1996), we obtained a catalogue of $``$ 200,000 objects with positions and the SExtractor photometry MAG–AUTO (for details on the construction of the catalogue see Wolf et al. 2001b, ). This $`R`$-band selected source catalogue was then used to obtain SEDs across 17 passbands with a photometric technique specifically tailored for measuring colour indices with high signal-to-noise and as little as possible interference by changing observing conditions. To this end, we projected the object coordinates into the frames of reference of each single exposure and measured the object fluxes at the given locations using with the well-established MPIAPHOT approach of a seeing-adaptive weighted aperture (Röser & Meisenheimer, 1991).

The COMBO-17 photometric calibration is based on a system of faint standard stars in the COMBO-17 fields, which we tied to spectrophotometric standard stars during photometric nights. Our standards were selected from the Hamburg/ESO survey database (Wisotzki et al., 2000) of digital objective prism spectra. By having standard stars within each survey exposure, we were independent from photometric conditions for imaging. More details of the data reduction will be provided in a forthcoming technical survey paper (Wolf et al., in prep). In the present paper, all magnitudes are quoted with reference to Vega as a zero point.

## 3 The quasar catalogue

### 3.1 Sample definition

The quasar sample is extracted from the full survey catalogue purely on the basis of spectral information, thus deliberately ignoring morphological evidence. We believe that the classification of the ‘fuzzy spectra’ based on 17 filters, outlined below, allows to clearly differentiate between stars, galaxies and quasars, providing a safer separation between the object classes than morphological criteria. This is particularly important in the context of our interest in low-luminosity AGN which may well appear extended on a deep $`R`$-band image. Our purely spectroscopic classification approach ensures that the sample will not be heavily biased against such objects.

The sample finally used for all analyses in this paper is defined by limits in magnitude and redshift. Objects are selected to have a magnitude of $`R>17`$ to avoid the saturation regime in the individual frames, and $`R<24`$ which is where the completeness of the quasar identification has dropped below 30% (see Sect. 3.3). Since quasars show brightness variations, the sample selection will depend on the epoch of observation at the faint end. Here, we have used our R-band photometry from January 2000.

The sample is further limited to redshifts of $`z>1.2`$, since host galaxies may contribute significantly to the spectra at lower redshifts, where the 4000 $`\mathrm{\AA }`$-break is still contained within the filterset. Our templates currently contain only pure quasar and pure galaxy spectra, but no mixed templates with contributions from both. Therefore, identification and redshift estimation of low-luminosity quasars at $`z<1.2`$ is not straightforward at this point, and complicated completeness issues arise in the low-$`z`$ domain.

In order to ensure minimum contamination from non-active galaxies, we have set our probability threshold for an object to be classified as quasars quite high. As a result, we have eliminated many trustworthy low-luminosity AGN from the sample, but keep a well-controlled selection of higher-luminosity QSOs. We have chosen this conservative approach because we have not yet obtained spectroscopic confirmation redshifts for a sufficiently large sample to understand the selection function of low-luminosity Seyfert galaxies well enough.

The catalogue resulting from this selection contains 192 quasars between $`z=1.2`$ and $`z=4.8`$, with a median redshift of $`z=2.23`$. Examples of quasar filter spectra are shown in Fig. 2, and the Hubble diagram for the full sample is presented in Fig. 3.

### 3.2 Classification and redshift estimation

The photometric measurements from 17 filters provide low-resolution spectra for each object which are analysed by a statistical technique for classification and redshift estimation based on spectral template matching (for details see Wolf et al. 2001a, ). Meanwhile, we have improved the template library for quasars by deriving it from the recent SDSS QSO template spectrum (vanden Berk et al., 2001), rather than from the earlier used pure emission line contour by Francis et al. (1991). This leads to a better detection and redshift estimation of quasars at $`z<2.5`$ where the ‘little blue bump’ may render the spectral shape between the prominent emission lines as quite different from the power-law which we assumed for the quasar templates previously. Since the CADIS work (Wolf et al., 1999), we have also learned that photometric redshifts for quasars are strongly influenced by flux variability.

Since the multi-colour observations were collected over a period of two years, and quasars as well as some stars show variability, it was necessary to correct for variability when constructing the 17-filter spectra for classification. For this purpose we used the $`R`$-band observations which are available for each observing epoch, allowing us to completely map the variability at least in this band. When constructing the SEDs of variable objects, we related the measurements of observed photometric band to the $`R`$-band magnitude obtained in the same observing run. As a result, the SEDs are not distorted by long-term magnitude changes. This variability correction works well only as long as flux variations are not accompanied by major changes in spectral shape. Furthermore, we also cannot account for short-term variability on a time scale of days, as within each observing run we have obtained only one $`R`$-band image.

### 3.3 Completeness correction

The subsequent analysis includes only objects with successful $`z`$ estimates. It is therefore critical to understand for which quasars the data permit such a classification. Given the photometric properties of the survey, we can produce mock catalogues of stars, galaxies and quasars with realistic photometric errors, and investigate the classification performance as a function of object type, redshift and magnitude. The 17-filter-spectra in the mock catalogues are constructed from the library templates, using the empirically derived observational error distributions. With this approach, the completeness of the classification and redshift estimation can be derived from Monte-Carlo simulations (see Wolf et al. 2001a, ). We implicitly assume that the spectral templates truly resemble observed objects, and if that is not the case, our simulation will be too optimistic.

In fact, we can test whether our completeness maps are realistic, given that within our fields and selection limits 12 broad-line AGN are known from spectroscopic observations in the CDFS (Hasinger 2002, priv. comm.; Szokoly et al, in prep.). The map predicts that 10 out of 12 objects should be identified, while in fact we recover 8. Two objects are missing from our sample although they lie in regions of high expected completeness: A QSO with broad absorption lines (BALs) at $`z=3.6`$ and a Sy-1 galaxy at $`z=1.2`$. The BAL QSO was misclassified as a star because of its unusually star-like colours. The Sy-1 galaxy resides just at our low redshift limit which we adopted to avoid incompleteness arising from the pure AGN spectra getting contaminated by host galaxy contributions. We conclude that on the whole our map is realistic, but that (a) rare QSOs with unusual colours could still escape our attention, and (b) right at the low-redshift limit the level of completeness might be reduced compared to our maps.

The product of the simulations for quasars is a completeness map (selection function), providing a formal probability $`C(R,z,SED)`$ that a quasar of given intrinsic properties is recovered, in bins of observed $`R`$-band magnitude and redshift. Our quasar template library contains a range of emission line strengths and a range in spectral indices or continuum colours. Effective spectral indices depend on redshift as the quasar continuum is not a pure power law. However, for a quasar at $`z=2.0`$ the range of template $`BI`$ colours runs from about $`+0.35`$ to $`+1.75`$, corresponding roughly to power law indices of $`\alpha =1.66\mathrm{}+0.4`$. After we recognized that the simulated completeness depends very little on the spectral index, we collapsed the completeness function into $`C(R,z)`$ removing the explicit SED dependence (see Fig. 4).

At most redshifts, our classification algorithm is more than 90% complete for quasars with $`R23`$. The contour line for 50% completeness ranges around $`R23.5`$. The redshift dependence of the completeness shows conspicuous oscillations above redshift 2, which are caused by the strong signature of the Lyman-$`\alpha `$ emission line migrating through the filter set and alternating between visibility in a medium-band filter and invisibility when it falls between two neighboring medium-band filters. Whenever the line is visible quasars can be identified to fainter levels due to its contrast.

At certain redshifts, the selection function reaches values above 100 %. This occurs when redshift aliasing creates local overdensities in the estimated $`z`$-distribution based on a flat underlying simulated distribution. Such ‘overcomplete’ regions are always accompanied by adjacent ‘undercomplete’ zones, so that the total number of objects is conserved.

### 3.4 Influence of limited redshift accuracy

For the analysis it is crucial to check to which extent the limited redshift accuracy of $`\sigma _z0.03`$ (as determined from the simulations) could affect inferences about luminosity functions. Two major aspects need to be explored, redshift aliasing and catastrophic mistakes, which are both irrelevant for interpreting well-exposed data from an aperture spectrograph, but could both play a role in our case:

1. Aliasing results from structures which are finer than the redshift resolution, violating the sampling theorem; it appears as fake structure on a scale that is typically slightly larger than the resolution. This effect is doubtlessly present in our data, as is also visible from the $`z`$-dependent oscillations in Fig. 4. In this paper, we avoid dealing with the problem altogether by only considering redshift bins of typically $`\mathrm{\Delta }z=0.6`$, an order of magnitude larger than the resolution limit.
2. Catastrophic mistakes occur in certain regions of colour space where different interpretations can be assigned to the same colour vector and probabilistic assumptions are used to make a final redshift assignment. Obviously, a number of cases could lead to the assignment of a wrong redshift, but their effect on luminosity functions should be relatively small unless a major fraction of objects were affected. The only real impact would be in a situation in which many low-redshift objects of medium luminosity were wrongly assumed to reside at high redshift, boosting the abundance of the rare luminous quasars. We see no evidence for any such effect to be important.

A large-scale performance check of our spectrophotometric redshifts based on a substantial sample of spectroscopic redshifts is still pending. However, we have been able to test our classifications and redshifts against two samples of albeit limited sizes that have spectroscopic data. In the S11 field, the COMBO-17 sample overlaps at its bright end with the ‘2QZ’ quasar survey (Croom et al., 2001), and in the CDFS we found that a number of spectroscopically identified X-ray sources coincide with faint COMBO-17 AGN (Hasinger 2002, priv.comm.; Szokoly et al, in prep.). Fig. 5 shows a comparison of spectroscopic and multi-colour redshifts for all 22 objects from the COMBO-17 QSO sample with $`R<24`$ for which spectra exist. Only one of these 22 objects, classified as a QSO at $`z0.6`$ by COMBO-17 (i.e. outside the redshift range considered in this paper), was actually a misclassification; a 2QZ spectrum identified this object correctly as a White Dwarf/M dwarf-pair. The other 21 objects were found to be QSOs, and their redshifts agree all within $`\pm 0.05`$. This comparison suggests that the accuracy of COMBO-17 redshifts is certainly no worse than the assumed $`\sigma _z=0.03`$, and the rate of catastrophic mistakes is probably well below 10 %.

The CDFS contains several further faint Seyfert galaxies, particularly at $`z<1`$, most of which were recently identified as optical counterparts to faint X-ray sources. These are very challenging objects for the multi-colour redshift estimation, as their host galaxies will contribute significantly to their overall SED. Clearly, more work is needed to tackle these objects, which are therefore explicitly excluded from the scope of this paper.

## 4 Number counts

From the AGN sample and the completeness map described above it is straightforward to compute surface densities as a function of apparent magnitude, using the relation

$$𝒩(<R)=\underset{i}{}\frac{1}{A_i}=\frac{1}{A_0}\underset{i}{}\frac{1}{C(R_i,z_i)}$$
(1)
where the summation runs over all AGN brighter than $`R`$. $`A_0`$ is the formal total survey area of 0.78 deg<sup>2</sup>; multiplying this with the value of the completeness map $`C(R_i,z_i)`$ for object $`i`$ yields the estimated ‘effective survey area’ $`A_i`$ for that particular source. In other words, the decreased probability of detecting very faint AGN, particularly relevant for certain redshift ranges, is simply interpreted as being equivalent to a reduced survey area for such objects.

Results are shown in Fig. 6 and listed in Table 3. Besides the full sample of all $`z>1.2`$ AGN, we also show the number counts for subsamples split at $`z<2.2`$ and $`z>2.2`$, the typical high-redshift limit of UV excess surveys. The two subsamples show very different trends towards the faint limit: While the low-$`z`$ objects dominate at brighter magnitudes, their surface density increases only slowly with decreasing flux level. On the other hand, the cumulative number of AGN with $`z>2.2`$ is a strong function of magnitude, and quite well described by a power law $`𝒩(<R)R^{0.31}`$, corresponding to a power law index of $`1.75`$ for the differential number-flux relation. Towards the faint end it even seems as if the slope might be increasing, but this is not formally significant. Nevertheless, for $`R>23`$ it is clear that high-redshift AGN with $`z>2.2`$ become the dominant population.

## 5 Estimation of rest frame luminosities

For each object in the sample we have individual SEDs from the 17-filter spectrophotometry at our disposal. In order to derive absolute magnitudes for the subsequent analysis, we decided to tie all luminosities to the UV continuum level at $`\lambda _{\mathrm{rest}}=145`$ nm. This was realised by integrating the best-fitting redshifted template over a synthetic narrow rectangular passband at 143 nm–147 nm, thus avoiding the Si iv/O iv emission line at 140 nm. This procedure enabled us to directly measure rest-frame luminosities over a redshift range $`1.4z5`$ without any need for extrapolation.

While such individually determined luminosities may suffer less from biases arising in uncertain assumptions on AGN spectral energy distributions than some previous samples, they considerably complicate the statistical exploitation. Although most of our objects are detected in all or nearly all of the 17 filter bands, the sample is presently *defined* only in the $`R`$-band: (1) The flux limit is homogeneously truncated at $`R<24`$. (2) The completeness map has the $`R`$ band magnitude as one of its two independent parameters. It is also assumed in the algorithms used for luminosity function estimation that the survey is defined by just one flux limit. We have therefore explored whether it is possible to estimate absolute UV magnitudes $`M_{145}`$ from the measured $`R`$ band flux, using a modified but basically traditional $`K`$ correction approach.

For this purpose we computed for each quasar its absolute magnitude directly from the distance modulus, $`M_0=R(mM)`$, equivalent to assuming a $`K`$ correction for a power-law spectrum with slope $`\alpha =1`$. We then plotted the differences $`\mathrm{\Delta }M`$ between $`M_0`$ and $`M_{145}`$ (the spectrophotometric estimates) against source redshift; this is shown as the distribution of points in Fig. 7. This plot shows that at least at redshifts $`z2`$ there is a well-defined relation between $`z`$ and $`\mathrm{\Delta }M`$. When fitting a constant line through these points, as shown in Fig. 7, we get an overall rms scatter of 0.24 mag (0.13 mag for $`z>2.2`$ and 0.26 mag for $`z<2.2`$), which can only insignificantly be improved by fitting e.g. a low-order polynomial. This line defines then an empirical correction relation $`K(z)`$, so that for the redshift range of interest in this paper we can use $`R`$ band magnitudes to estimate $`M_{145}`$ with an accuracy of 0.24 mag, even without using any SED information.

Most earlier studies of the optical AGN luminosity function, especially those focussing on lower redshifts, have expressed their results in terms of blue magnitudes $`M_B`$. In order to facilitate a statistical comparison, we obtained a crude estimate for the offset $`M_BM_{145}`$ by repeating the above procedure for the rest-frame $`B`$ band. Note that for all $`z>1.2`$ objects this involves a good deal of extrapolation, and the diagram corresponding to Fig. 7 shows more than 0.6 mag of scatter around the mean trend. Nevertheless, the offset of 1.75 mag thus determined is actually identical (to 0.01 mag) to the result when the same quantity is measured in the mean quasar energy distribution of Elvis et al. (1994). In conclusion this means that all quoted absolute magnitudes $`M_{145}`$ can be converted into $`M_B`$ using the relation $`M_B=M_{145}+1.75`$.

In the subsequent analysis we use two sets of cosmological parameters, an Einstein-de Sitter universe with $`H_0=50`$ km s<sup>-1</sup> Mpc<sup>-1</sup>, $`\mathrm{\Omega }_m=1.0`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, and a flat ‘concordance model’ with $`H_0=65`$ km s<sup>-1</sup> Mpc<sup>-1</sup>, $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. While the former is now physically almost obsolete, it still has some relevance as a ‘yardstick’ model, as most earlier studies of AGN evolution were expressed preferentially in these terms. Actually we found that the results change very little when switching between the two models, except for small shifts in both axes, and our displayed results always refer to the ‘concordance universe’ unless explicitly stated otherwise. Figure 8 shows the distribution of absolute magnitudes vs. redshift for the entire AGN sample. In terms of the conventional (if arbitrary) distinction between high-luminosity quasars and low-luminosity Seyferts around $`M_B23`$, our sample covers just the region close to this dividing line. It is therefore the first optically selected AGN sample of substantial size to probe this important regime of nuclear activity in galaxies. It also matches very well the luminosity range covered by low-$`z`$ AGN surveys that provide a $`z0`$ reference (e.g., Köhler et al., 1997).

## 6 Evolution of the AGN luminosity function

### 6.1 Non-parametric estimates

#### 6.1.1 Method of computation

We employ the usual $`1/V_{\mathrm{max}}`$ estimator Schmidt (1968) to give the space density contributions of individual objects. Luminosity functions are then readily obtained by forming the appropriate sums: Denoting the luminosity-binned differential LF as $`\varphi (M)`$ and the cumulative LF as $`\psi (M)`$ we have

$`\varphi (M,z)`$ $`=`$ $`{\displaystyle \underset{M\mathrm{\Delta }M}{\overset{M+\mathrm{\Delta }M}{}}}{\displaystyle \underset{z\mathrm{\Delta }z}{\overset{z+\mathrm{\Delta }z}{}}}{\displaystyle \frac{1}{V_i}},`$ (2)
$`\psi (M,z)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{M}{}}}{\displaystyle \underset{z\mathrm{\Delta }z}{\overset{z+\mathrm{\Delta }z}{}}}{\displaystyle \frac{1}{V_i}}`$
where each value of $`V_i`$ is the total comoving volume within which object $`i`$ would still be included in the sample. The value of $`V_i`$ follows from integrating the survey selection function $`C(R,z)`$, introduced in Sect. 3.3, over all relevant redshifts:

$$V_i(R,z)=A_0_z^{z_{\mathrm{max}}}C(R,z)\frac{dV}{dz}𝑑z$$
(3)
(cf. Wisotzki, 1998) where $`A_0`$ is again the formal total survey area, $`A_0=0.78`$ deg<sup>2</sup> in the present case. The assigned error bars are based purely on Poissonian shot noise due to the limited AGN counts within each bin, i.e. we ignored errors in magnitude or redshift estimation. The error bars for a given sum are then given by

$$\sigma _\varphi =\sqrt{\underset{i}{}1/V_i^2(M,z)}.$$
(4)
It is well known that binned estimates of luminosity functions suffer from a number of biases. In particular, two effects are worth being recalled, both of which can be seen as two different variants of Malmquist bias and lead to the observed LF appearing flatter than the intrinsic one:

* Binning in magnitude and a steep LF create a positive bias in the derived space densities, due to a shift in effective bin centres towards lower luminosities.
* Differential evolution within redshift shells leads to an overestimation of space densities where the LF is steep.

While the first effect is absent in the cumulative LF representation, the latter is unavoidable as long as no model assumptions for the evolution of space densities are made. It should thus be kept in mind that the nonparametric LFs discussed in the following section do not necessarily show the data ‘as they are’, but just form the most conventional way of presenting luminosity function data.

#### 6.1.2 Results

Inside the redshift interval $`1.2<z<4.8`$ covered by the COMBO-17 sample we defined six redshift shells of equal sizes to trace the evolution of the QLF. Figure 9 gives a synopsis of the results; for the benefit of interested readers who wish to use our data for their own computations we present the tabulated binned LFs in Table 4. This table shows also how many objects contributed to each bin. The numbers given in this table show that the differences between adjacent redshift shells are not huge, i.e. evolution is relatively moderate over this redshift range. We therefore have deviated from the usual practice to plot all LFs into one frame, as any real trend would be very hard to discern in such a diagram. Instead, we created individual subpanels for each redshift shell and plot each luminosity function separately, but together with the $`1.8<z<2.4`$ LF taken as a reference. This visualisation shows that the measured values of the LF are below the corresponding reference values in nearly all data points. In other words, we detect an unambiguous maximum in the comoving AGN space density near these redshifts, and a significant drop towards both lower and higher $`z`$. The existence of such a maximum becomes even clearer in Fig. 11 where integrated space density at given lower bound in luminosity (taken directly from the cumulative LFs in each redshift shell) is plotted directly against $`z`$. We return to a more detailed discussion of the evolution of space densities and LF shape properties below.

### 6.2 Parametric analysis of the luminosity function

#### 6.2.1 Ansatz

Within certain bounds and accuracy limits, an observed luminosity function can usually be described quite well by some simple analytic expression, allowing one to compress the results into a few well-determined numbers, with the positive side effect that the above mentioned Malmquist-type biases due to data binning can be avoided by obtaining luminosity function and evolution parameters simultaneously from an observed sample. This was first demonstrated by Marshall et al. (1983) for a very simple power law model of the evolving QLF. Since the parametric forms employed for the analysis in this paper are somewhat non-standard, we use the following paragraphs to spell out the adopted ansatz.

##### Shape of the luminosity function:

The most common analytic description of the QLF is a double power law, involving a bright-end slope $`\gamma _1`$, a faint-end slope $`\gamma _2`$, and a smooth turnover at a characteristic ‘break’ luminosity $`M^{}`$. We have explored this ansatz and decided not to use it for the present sample, mainly because we fail to identify a well-defined break luminosity in the data. Direct fitting of a double power law LF to any of the redshift shells subsamples shows that $`M^{}`$ is a very ill-constrained quantity that often takes a value outside the luminosity range covered by the sample; furthermore, even minor modifications in the subsample definition can cause major changes in $`M^{}`$. Since the luminosity functions in Fig. 9 undoubtedly show *some* indication of curvature, we decided to parametrise the LF, at given $`z`$, as a polynomial in $`M_{145}`$:

$$\mathrm{log}\varphi =\underset{i=0}{\overset{n}{}}A_i\mu ^i.$$
(5)
where $`\mu =M_{145}M^{}`$ and $`M^{}`$ is a constant which can be chosen freely but may also contain the effects of luminosity evolution (see below). In this representation, $`A_0=\mathrm{log}\varphi _0`$ is the space density normalisation at $`M_{145}=M^{}`$, which in our case has mostly been set to $`M^{}=25`$ (roughly corresponding to $`M_B=23.25`$). For $`n=1`$, the LF is a single power law with slope $`\gamma =2.5A_11`$.

##### Evolution parameters:

The last decade has seen a sometimes hot debate over the question whether ‘pure luminosity’ or ‘pure density’ evolution were more appropriate. We have implemented both parametrisations as well as mixed forms. As independent variable we use

$$\zeta =\mathrm{log}\left(\frac{1+z}{1+z_{\mathrm{ref}}}\right)$$
(6)
which vanishes for $`z=z_{\mathrm{ref}}`$. For our analysis we avoid extrapolating to redshifts outside the sample definition interval and adopt $`z_{\mathrm{ref}}=2`$. Pure density evolution is then expressed as polynomial expansion of the space density normalisation:

$$\mathrm{log}\varphi _0(z)=\underset{j=0}{\overset{m_D}{}}C_j\zeta ^j.$$
(7)
For $`m_D=1`$, the parameter $`C_1`$ is equal to the density evolution index $`k_D`$ used in several earlier analyses.

On the other hand, pure luminosity evolution involves a redefinition of the LF parameter $`M^{}`$:

$$M^{}(z)=M^{}(\zeta =0)+\underset{j=1}{\overset{m_L}{}}B_j\zeta ^j.$$
(8)
(Note that unlike for a ‘broken power law’, $`M^{}(\zeta =0)`$ is here an arbitrary constant and not a free fit parameter.) Again, the case of $`m_L=1`$ corresponds to a standard form already used by several previous authors, $`M^{}(z)(1+z)^{k_L}`$. In this case the luminosity evolution index is $`k_L=0.4B_1`$.

##### Fitting method:

We followed a maximum likelihood approach (Marshall et al., 1983) to search for optimal parameter combinations, using a modified downhill-simplex minimisation scheme (Press et al., 1992). Several goodness-of-fit tests were applied to check whether the best-fit model was adequate. These included one- and two-dimensional Kolmogorov-Smirnov tests for the distributions over the $`(M,z)`$ plane, as well as a $`\chi ^2`$ test that directly compares observed and predicted numbers of AGN in redshift/luminosity bins. A fit was considered satisfactory only when it passed all of the statistical tests with a probability for the validity of the null hypothesis of $`>10`$ %. This highly conservative approach ensures that the missing distribution-free properties of the KS test (see Wisotzki 1998 for discussion on Lilliefors et al., 1968) can be neglected.

#### 6.2.2 Results

Starting with the simplest possible luminosity function, we found that a single power law LF as an overall shape model had to be rejected as there is weak but significant curvature in the LF, especially for the lower redshift ranges, $`z2.4`$. But already a 2nd order poynomial gives a statistically adequate description of the LF at all redshifts. Note that this form involves one free parameter less than the case of a double power law LF.

As a first evolution mode we explored pure density evolution (PDE). A global PDE fit over the entire redshift range $`1.2<z<4.8`$ was achieved with again just a 2nd-order polynomial to reproduce the maximum of comoving space density around $`z2`$. Alternatively, the sample can be described equally well by pure luminosity evolution (PLE), albeit requiring one more parameter than PDE. Our best-fit PLE model features a 3rd order polynomial for the function $`M^{}(z)`$. The predicted parametric luminosity functions are compared with the binned nonparametric estimates in Fig. 10. Both fits provide excellent descriptions of the data, with goodness-of-fit probabilities of $`p>50`$ % for all statistical tests (see Tab. 5 for best-fit models). The differences between PDE and PLE luminosity functions become significant only outside the luminosity range sampled by the COMBO-17 data.

In Fig. 11 we show the cumulative luminosity function as a function of redshift, for different limiting luminosities, again together with the corresponding nonparametric estimates. For simplicity, only the PDE result is plotted. Since the fit is presently unconstrained below $`z=1.2`$, we have marked the extrapolation into this region by a dotted line. The fit is completely consistent with the individual binned datapoints.

## 7 Discussion

### 7.1 The maximum of AGN activity

The dominant feature in Fig. 11 is the peak of comoving AGN space densities around $`z2`$. Although its existence was beyond doubt for a long time, rarely has it been possible to locate the maximum *within a single survey*. UV excess-selected samples typically reach up to $`z=2`$–2.5 and mainly trace the rise at low $`z`$; conversely, most dedicated high-redshift surveys start being effective only beyond $`z3`$ and just show the decreasing branch. The COMBO-17 sample combines properties of both search techniques and covers enough redshift range that the peak of AGN activity is clearly bracketed. The PDE and PLE fits both place the maximum at a redshift of $`z_{\mathrm{max}}=2.1`$.

An important issue for the astrophysical interpretation of AGN evolution is a possible luminosity dependence of the peak location. For example, in hierarchical structure formation one expects objects of higher mass objects to form later, which in a simple scenario of correlated masses and luminosities should manifest in high-$`L`$ AGN to show a peak at lower redshifts. This is certainly not the case in our data; on the contrary, one might be tempted to speculate from Fig. 11 that there could be a gradual shift of the maximum towards higher redshifts when the luminosity limit is increased.

We have tried to model such a trend by including explicit luminosity-dependent density evolution parameters, e.g. in the form of a linear dependency between absolute magnitude and $`z_{\mathrm{max}}`$. While the data are statistically *consistent* with a moderate shift, the best-fit model was always very close to simple PDE, despite the additional degrees of freedom. We conclude that there is no indication for a dramatic difference between the space density peaks of AGN of different luminosities, within the range covered by COMBO-17.

A similar but much stronger effect, a shift of maximum space density towards *lower* redshift for very low luminosity AGN, has recently been reported for deep X-ray selected samples (Cowie et al., 2003; Hasinger et al., 2003). While a detailed comparison between the properties of X-ray and optically selected AGN samples is beyond the scope of this paper, we just note two aspects which need to be taken into account: Firstly, the new X-ray selected samples probe even much deeper than COMBO-17 into the population of AGN with very low luminosities. Assuming an Elvis et al. (1994) AGN spectrum, our low-luminosity limit of around $`M_{145}=23`$ corresponds to a 2–8 keV luminosity of roughly $`L_x43.6`$; this is just the point where the difference becomes visible, according to Cowie et al. and Hasinger et al.

But secondly, the X-ray and optical luminosity functions of high-redshift AGN also show significant dissimilarities. In particular, the X-ray LF of type 1 (broad-line) AGN near $`\mathrm{log}L_x44`$ is very nearly constant (cf. Cowie et al.; Hasinger et al.), while the optical LF around the corresponding $`M_{145}24.6`$ is still significantly rising (cf. Fig. 9). This may be indicative of different spectral shapes between low- and high-luminosity, or between low- and high-redshift AGN, all leading to nontrivial $`L_x`$ to $`M_{145}`$ conversions. Thus, while on one hand our $`R<24`$ sample may still be not quite deep enough to see similar effects as found in the X-ray domain, it is not even obvious that exactly the same effects should be expected for a yet deeper sample.

### 7.2 Evolution of the AGN luminosity density

The contribution of AGN to the metagalactic radiation field is probably relevant for several wavebands. Recent satellite missions were very successful in resolving the extragalactic X-ray background as mainly due to distant low-luminosity AGN (e.g. Miyaji et al., 2000). Likewise, the diffuse ionising UV background is believed to be strongly influenced by AGN, but a quantitative synthesis still relies on several assumptions and extrapolations. One of the principal uncertainties has always been the shape of the low-luminosity end of the AGN luminosity function. In order to estimate the AGN luminosity density $`\varrho _L(z)`$, one has to evaluate the integral

$$\varrho _L(z)=_LL\varphi (L,z)𝑑L.$$
(9)
The integral diverges only for a luminosity function slope $`\gamma <2`$, but the relative contributions of AGN at different luminosity levels depend critically on the overall shape of the LF.

The COMBO-17 AGN sample goes deep enough that, for the first time, the quantity $`L\varphi `$ can be safely integrated without depending on heavy extrapolation into the unobserved range. We have computed both binned nonparametric as well as parametric estimates of $`\varrho _L(z)`$, which we show in Fig. 12. We present the monochromatic luminosity density $`\varrho _{L(145)}`$, based on the value of $`\lambda L_\lambda `$ evaluated at $`\lambda =145`$ nm (which in turn is directly derived from $`M_{145}`$). Assuming a typical QSO spectrum such as the one given by Elvis et al. (1994), the quantity $`\varrho _{L(145)}`$ can also be extrapolated into a frequency-integrated luminosity density. Since we are particularly interested in the AGN contribution to the cosmic production of hydrogen-ionising photons, we provide the necessary conversion:

$$L_{\lambda <91.2\mathrm{nm}}=\mathrm{\hspace{0.25em}1.45}\times \lambda L_{\lambda =145\mathrm{nm}}.$$
(10)
Notice that we have individual SED information on each of the COMBO-17 objects, which in principle could be used to derive a more accurate estimate of $`L_{\mathrm{ion}}`$. This will be done in a later paper dedicated to the synthesis of the UV background from the COMBO-17 AGN sample.

For the nonparametric estimates we computed $`\varrho _L`$ by summing over the luminosity-weighted inverse volumes $`L/V_i`$ up to the survey limit; the resulting numbers should set lower limits to the full integral over all luminosities. The two adopted parametric forms of PDE and PLE can be directly integrated to provide smooth functions $`\varrho _L(z)`$. It is worth noting that the integral effectively converges within the COMBO-17 survey limits; the difference between setting the low-luminosity integration boundary at $`M_{145}=23`$ (approximately the COMBO-17 limit) and at $`M_{145}=10`$ (or $`L=0`$) is only 2 %. In other words, the constraints on the faint end of the AGN luminosity function from COMBO-17 are already sufficient to estimate the total UV radiative output from AGN. This is also illustrated by the relatively good agreement in Fig. 12 between the binned estimates (which, as said above, are formally just lower limits to $`\varrho _L`$) and the PDE model.

On the other hand, the flat slope of the LF makes the contribution of brighter AGN to $`\varrho _{L(145)}`$ anything but negligible. In fact, the inclusion or exclusion of individual high-luminosity objects makes a noticeable difference for the binned estimate. (This is much less so for the parametric estimates, because of the uniform weighting inherent in the maximum likelihood procedure.) The importance of higher luminosity objects becomes particularly apparent when comparing the PDE- and PLE-derived relations. The substantial and significant differences seen in Fig. 12 are *not* due to faint-end extrapolation effects, but entirely originate in the much flatter LF slope of the PLE model at *brighter* magnitudes, already documented in Fig. 10 above. We reiterate that these discrepancies are due to the limited power of the COMBO-17 dataset to discriminate between evolution models. They will be resolved as soon as constraints from other, brighter AGN surveys are combined with the COMBO-17 sample.

### 7.3 Comparison with other surveys

#### 7.3.1 Redshift $`1.5<z<4.5`$

In terms of the luminosity range, our quasar sample pushes to fainter limits than any previous survey. Our targets mainly have observed magnitudes of $`R=20\mathrm{}24`$, while previous surveys either observed down to $`R20`$ if they derived luminosity functions, or observed to $`R22.5`$ but produced object lists which did not constrain the luminosity functions very much further. Therefore, this work enters a new regime of studying low-luminosity AGN at intermediate to high redshifts.

A first attempt to conduct an AGN survey with a strategy similar to COMBO-17 was performed, albeit on a much smaller scale, in the course of the CADIS survey (Wolf et al. 1999), where a sample of 12 QSOs at $`z>2`$ and $`R<22`$ was identified within an area of 250 arcmin<sup>2</sup>. Their resulting surface density is twice as high as for COMBO-17, but because of the small sample size, the results of CADIS and COMBO-17 are still formally compatible. The discrepancy underlines, however, the importance of cosmic variance and hence the need to obtain large samples.

Major previous work producing luminosity functions include the 2QZ at $`z2.3`$ (Boyle et al., 2000), the work by WHO at $`2.0<z<4.5`$, and by SSG at $`z>2.7`$ as well as the results from the high-redshift SDSS sample at $`z>3.6`$ by Fan et al. (2001). As we have little overlap with all these previous surveys in luminosity, we cannot test directly to what extent our and their luminosity functions coincide. However, at the bright end our LF is completely consistent with an extrapolation of the SDSS-based LF given by Fan et al., with similar slope and similar normalisation. SSG cover a slightly wider range of redshifts down to $`z=2.7`$, but there still is not much overlap, so as before we need to extrapolate. The SSG slope is steeper than that of Fan et al., so that the SSG prediction for our regime is even further above the COMBO-17 results.

On the lower-redshift side, our results smoothly connect to the 2QZ results. At intermediate redshift again, WHO used a broken power-law to characterize the QLF, where they found a steep end slope of $`\alpha =1.67`$ and a faint end slope of $`\beta =0.45`$, which is mostly constrained by the probably non-optimal assumption of a broken power-law. The bright end of our LF varies around 0.4 to 0.7 in their units and can be considered consistent with the WHO faint-end.

#### 7.3.2 Quasars at $`z>4.5`$

The redshift range above 4.5 is subject to several studies from dedicated surveys to find faint high-redshift quasars. At highest luminosities the SDSS has delivered samples of six objects in 180 $`\text{deg}^2`$ at $`i^{}<20`$ (Fan et al., 2001) and 29 objects in $`700\text{deg}^2`$ at $`i^{}<20.5`$ (Andersen et al., 2001). Medium-deep surveys across cover smaller areas, among which the BTC40 has so far reported two $`z>4.5`$ QSOs at $`I<21.5`$ across 36 $`\text{deg}^2`$ and still have a list of fainter candidates at $`I22`$ to follow up (Monier et al., 2002). The Oxford-Dartmouth Thirty Degree Survey (ODT) also aims at finding $`z>4.5`$ quasars down to $`I<22`$ on 30 $`\text{deg}^2`$ (Dalton, priv. comm.). Monier et al. (2002) demonstrate the consistency of their currently available small numbers with the SDSS result and predict on the basis of the SDSS-LF a surface density of 0.026 $`/\text{deg}^2`$ at $`z>4.5\mathrm{}5.0`$ and $`I<19.9`$, and much less at $`z>5`$.

COMBO-17 covers only 0.78 $`\text{deg}^2`$ currently, but reaches a bit deeper. An extrapolation of the SDSS luminosity function suggests we should find 2.0 quasars per $`\text{deg}^2`$ at $`I<23`$ and $`z>4.5`$. Our observation of some curvature in the fainter domain of the luminosity function reduces that prediction to $`1.25/\text{deg}^2`$, so one or two objects is the total number to be expected in the COMBO-17 dataset.

The present sample contains three objects at $`z>4.5`$, all in the S11 field (see Fig. 14 for filter spectra). They all have $`I<22`$, although our selection should be complete to $`I<23`$. Basically, this observation is statistically consistent with an extrapolation of both the SDSS LF and our own. As we can not draw strong conclusions about the cosmic abundance of these objects from COMBO-17, they have been excluded from the LF results presented above.

### 7.4 Predictions for future deep surveys

Given a parametrised model for the evolving AGN luminosity function, it is straightforward to predict surface densities (numbers of AGN per unit solid angle) for other survey specifications. In Table 6 we provide a set of such numbers, computed for grid of $`R`$ band magnitudes and redshifts. We deliberately stretch these predictions to the limit of credibility, e.g. in the lowest ($`R<25`$) row or in the rightmost ($`z>5`$) column, in order to enable a direct comparison with possible future dedicated ultra-deep or very high redshift surveys. It should be understood that some of these numbers involve a considerable amount of extrapolation outside the areas covered by our data. For the same reason we list at each grid point the predictions of both PDE and PLE models, hoping that the difference between these approximately bracket the level of uncertainty, especially in the extrapolation regions.

At intermediate redshifts and flux levels, well sampled by COMBO-17 sources, the agreement of PDE and PLE is excellent, simply reflecting the fact that both are valid descriptions of the observed data. On the other hand, the two evolution modes predict substantially different numbers of faint high-redshift AGN, the PLE prediction exceeding the PDE prediction by a factor of several. This is only partly due to the uncertainties imposed by the limited dataset. Figure 10 shows why this effect is actually expected from the properties of the two chosen evolution modes: A luminosity function which fits the data well at $`z2`$, which then is displaced either vertically or horizontally to account for the substantial negative evolution beyond $`z>3`$, will result in dramatically different space densities of low-luminosity AGN.

## 8 Conclusions

We have presented work on the evolution of the quasar luminosity function which has novel aspects by addressing two of the main unresolved problems in quasar research: The location of the peak of quasar activity, and the evolution of low-luminosity objects, the bulk of the quasar population:

* Using the medium-band approach of COMBO-17 it was possible to break the colour degeneracies between stars and quasars around $`z3`$ which had prevented earlier optical surveys from mapping out the shape of the turnover in quasar activity from a single, homogeneous survey. As a result, we have measured a broad maximum around $`z2`$.
* By reaching faint magnitudes of $`R<24`$, we were probing the quasar LF and measured its slope down to rather low luminosity. The constraints on the faint end of the LF are fully sufficient to estimate the total UV radiative output. Although the two investigated evolution modes predict drastically different numbers of very faint ($`R25`$), very high redshift ($`z5`$) AGN, these have not much effect on the luminosity density integral. The main contribution to the UV background comes indeed from intermediate luminosity objects within $`M_{145}[25,28]`$.

However, a number of questions and details remain still open. Among these are:

1. The analysis of COMBO-17 catalogues themselves has so far omitted redshifts $`z<1.2`$, because low-luminosity Seyfert galaxies, also those of type 1, may remain undetected due to their complex SEDs consisting of a mix of contributions from the host galaxy and the active nucleus. It will be a technical issue for the COMBO-17 classification to improve on that in the near future.
2. At $`z>4.5`$ COMBO-17 reaches to sufficient depth, but does not cover sufficient area to make a significant contribution. Wider area surveys such as the BTC40, the ODT or a deep sub-area of SDSS are needed to identify a larger sample of faint high-redshift quasars.
3. Our sample alone fits well to both of the simplest evolutionary models for the quasar LF, the PDE and PLE models. Thus, further improvements in determining the overall shape of the LF and understanding its evolution on a broader basis will require the combination of the COMBO-17 sample with other samples at higher luminosity and lower redshift (to be addressed by a forthcoming paper).
4. While our survey is conducted in the rest-frame UV domain, one could wish to estimate the *bolometric* luminosity density by applying a bolometric correction factor. Using the Elvis et al. (1994) template SED, the relation would be $`L_{\mathrm{bol}}=4.5\times L_{145}`$. However, we have little reason to assume *a priori* that the quoted bolometric correction is universally applicable to all AGN including our faint sample, so at this stage we refrained from following this up further. For similar reasons we do not want to derive accretion rates here.
5. Probably, our sample breaks new ground in terms of depth and redshift coverage, however, it is clear that the results apply only to the evolution of broad-line type-1 AGN.
6. Deep X-ray surveys with the Chandra and XMM/Newton observatories have started to charter into new territory of AGN research. Comparing the first results of X-ray and deep optical samples such as COMBO-17, we see similarities but also a number of differences. Exploring the interplay between X-ray and optical domain for AGN will be a highly rewarding exercise, which will lead to a more complete understanding of the cosmic evolution of accretion-powered galactic nuclei.

In summary, understanding the accretion history of black holes in the nuclei of galaxies will certainly require some more quite fundamental work. However, establishing the evolution of faint quasars in the turnover epoch is one step forward.

###### Acknowledgements.

CW was supported by the PPARC rolling grant in Observational Cosmology at University of Oxford and by the DFG–SFB 439. We thank Prof. Hasinger for providing spectroscopic redshifts of Chandra X-ray sources to facilitate the cross check with multi-colour redshifts and foster our trust in their quality. We thank Lance Miller, Steve Warren and an anonymous referee for helpful comments improving the manuscript.


format: text

text:
# To appear in The Astrophysical Journal Scattering by Interstellar Dust Grains. I. Optical and Ultraviolet

## 1 Introduction

Interstellar grains scatter electromagnetic radiation. Reflection nebulosities are conspicuous at optical and UV wavelengths when dust is brightly illuminated by a nearby star. Dust clouds which are not unusually close to a star are illuminated by the general interstellar radiation field. Finally, the starlight scattered by dust in the diffuse interstellar medium constitutes the so-called “diffuse galactic light”. Observations of reflection nebulae, dust clouds, and the diffuse galactic light, provide a means of determining the scattering properties of interstellar grains, thereby testing models for interstellar dust.

Photoionization and photodissociation of molecules play a major role in interstellar chemistry, and the chemical structure of molecular clouds is therefore directly linked to the ability of ultraviolet starlight to penetrate into these dusty regions. Knowledge of dust scattering properties in the ultraviolet is therefore required for realistic modeling of interstellar clouds.

The nature of interstellar grains remains uncertain (see Draine 2003a, and references therein). This paper will examine the scattering properties for a grain model consisting of two separate grain populations – carbonaceous grains and silicate grains. With the grains approximated by homogeneous spheres with the size distributions found by Weingartner & Draine (2001; hereafter WD01), this grain model is consistent with the observed interstellar extinction, the observed infrared emission from interstellar dust (Li & Draine 2001, 2002), and the X-ray scattering halo observed around Nova Cygni 1992 (Draine & Tan 2003). The carbonaceous grains are assumed to be primarily carbon when the grains are large, but to extend down to very small sizes with the smallest grains being individual polycyclic aromatic hydrocarbon molecules. The scattering is dominated by grains with radii $`a100\mathrm{\AA }`$, containing $`10^6`$ atoms; carbonaceous grains in this size range are modeled using the optical properties of graphite.

The primary objective of this paper is to calculate the scattering and extinction properties of this dust model at infrared, optical, and ultraviolet wavelengths, and to make these results available for use in radiative transfer calculations and for comparison with observations. The X-ray scattering and absorption properties of this grain model are the subject of Paper II (Draine 2003b).

The adopted dielectric functions are presented in §2. The scattering properties of interstellar dust at optical and ultraviolet energies, as calculated for the carbonaceous-silicate grain model, are presented in §3. We show the scattering phase function at selected wavelengths from the SDSS z band ($`\lambda =8930\mathrm{\AA }`$) to the vacuum ultraviolet ($`\lambda =1820\mathrm{\AA }`$). Scattering properties are calculated for Milky Way dust with $`R_V=3.1`$, and also for models for dust in the LMC and SMC. In §3.2 we show the degree of polarization as a function of scattering angle for selected wavelengths.

The Henyey-Greenstein phase function has often been used to approximate the anisotropic scattering properties of interstellar dust. In §4.1 we show that the Henyey-Greenstein function has r.m.s. error $`<10\%`$ for $`0.47\mathrm{\mu m}<\lambda <0.94\mathrm{\mu m}`$, but has larger errors outside this range. We present a new analytic phase function (equation 5) with a wider range of applicability, with r.m.s. error $`<10\%`$ for $`\lambda >0.27\mathrm{\mu m}`$.

In §5 we collect observational determinations of the scattering albedo and $`\mathrm{cos}\theta `$ for dust in reflection nebulae, in dense clouds, and in the diffuse interstellar medium. Discrepancies among these determinations are noted, and possible reasons for this are discussed.

The principal results are summarized in §6

## 2 Dielectric Function

As discussed by WD01, the grain population is assumed to include very small grains with the optical properties of polycyclic aromatic hydrocarbon molecules (PAHs), plus larger grains which are approximated as carbonaceous or silicate spheres. The PAHs produce negligible scattering.

From the observed 3.4$`\mathrm{\mu m}`$ C-H stretch feature, Pendleton & Allamandola (2002) estimate that $``$85% of the C is aromatic, and $``$15% is aliphatic (chainlike). The graphite dielectric function will be used to approximate the optical and ultraviolet response of interstellar carbonaceous grain material. Scattering and absorption by the carbonaceous spheres is calculated using the dielectric tensor of graphite, using the usual “1/3-2/3 approximation” (Draine & Malhotra 1993).

The dielectric functions used here are taken from Paper II, which constructs self-consistent dielectric functions extending from microwave to X-ray energies, including realistic structure near X-ray absorption edges. The adopted dielectric functions for graphite and “astronomical silicate” are shown in Figures 1 and 2. These dielectric functions are close to the dielectric functions obtained previously by Draine & Lee (1984), although differing in detail.

## 3 IR-Optical-UV Scattering by Interstellar Dust

### 3.1 Angular Distribution of Scattered Light

Weingartner & Draine (2001; hereafter WD01) obtained size distributions of spherical carbonaceous and silicate grains which reproduce the observed extinction curve both in the local Milky Way and in the Large and Small Magellanic Clouds. Here we calculate the scattering properties of these dust mixtures.

The scattering properties of a particular dust mixture $`X`$ are characterized by the differential scattering cross section per H nucleon

$$\left(\frac{d\sigma _{\mathrm{sca}}(\lambda ,\theta )}{d\mathrm{\Omega }}\right)_X\underset{j}{}𝑑a\left(\frac{1}{n_\mathrm{H}}\frac{dn_j}{da}\right)_X\left(\frac{dC_{\mathrm{sca}}}{d\mathrm{\Omega }}\right)_{j,a,\lambda },$$
(1)
where $`n_\mathrm{H}^1(dn_j/da)da`$ is the number of grains of type $`j`$ per H nucleon with radii in $`(a,a+da)`$, and $`(dC_{\mathrm{sca}}/d\mathrm{\Omega })_{j,a,\lambda }`$ is the differential scattering cross section for grain type $`j`$, radius $`a`$, at wavelength $`\lambda `$, for scattering angle $`\theta `$, for a grain illuminated by unpolarized light. The grains are assumed to be spherical, and the differential scattering cross sections are calculated using Mie theory (Bohren & Huffman 1984), using the code developed by Wiscombe (1980, 1996).

In Figure 3 we show the differential scattering cross section per H nucleon for $`R_V=3.1`$ Milky Way dust, at the central wavelengths of SDSS z (8930Å), i (7480Å), r (6165Å), g (4685Å), and u (3550Å), Cousins I (8020Å) and R (6492Å), V (5470Å), and the F250W (2696Å), F220W (2220Å), and F25CN182 (1820Å) filters for the Space Telescope Imaging Spectrograph (STIS). The scattering becomes stronger and more forward-throwing at shorter wavelengths.

To see the sensitivity to variations in the dust mixture, Fig. 3 also shows $`d\sigma _{\mathrm{sca}}/d\mathrm{\Omega }`$ calculated for the WD01 dust mixtures for the “average LMC” and the SMC bar. The calculated scattering closely resembles the Milky Way scattering, but with an overall reduction of about a factor of 4 for the LMC, and a factor of 8 for the SMC, in line with the overall dust and metal abundance relative to the Milky Way.111 Ne/H is $`30\%`$ of solar in the LMC and $`14\%`$ of solar in the SMC (Dufour 1984; Kurt & Dufour 1998). In the LMC, $`E(BV)/N_\mathrm{H}`$ is $`24\%`$ of the local Milky Way value (Koorneef 1982; Fitzpatrick 1986) while in the SMC it is only $`13\%`$ of the local value (Martin et al. 1989).

Figure 4 shows the wavelength dependence of $`\sigma _{\mathrm{ext}}`$ (the total extinction cross section per H), the scattering albedo $`\sigma _{\mathrm{sca}}/\sigma _{\mathrm{ext}}`$, and the first and second moments222The second moment $`\mathrm{cos}^2`$ will be used in the new analytic phase function proposed in §4.2.  $`\mathrm{cos}\theta `$ and $`\mathrm{cos}^2\theta `$ for the scattered radiation. Six different grain models are shown, fitted to different observed extinction curves (see WD01 for details). There are considerable differences in albedo, $`\mathrm{cos}`$, and $`\mathrm{cos}^2`$ among the models. In particular, the SMC bar model shows a high albedo near $`4.6\mathrm{\mu m}^1`$, in contrast to the other five models which have a local minimum in the albedo at this frequency. This is because the SMC bar model differs from the other 5 models in lacking PAHs and $`a0.02\mathrm{\mu m}`$ graphite grains, as these are not allowed by the observed absence of a $`4.6\mathrm{\mu m}^1`$ extinction “bump”.

### 3.2 Polarization of Scattered Light

Even when a grain is illuminated by unpolarized light, the scattered radiation is generally polarized. The degree of polarization depends upon both the scattering angle and the wavelength of the radiation. The fractional polarization $`P(I_{}I_{})/(I_{}+I_{})`$, where $`I_{},I_{}`$ are the intensities of scattered light in polarization modes perpendicular or parallel to the scattering plane.

The polarization $`P`$ of the scattered light is shown in Fig. 5 as a function of scattering angle $`\theta `$ for Milky Way dust with $`R_V=3.1`$, at 11 different wavelengths. Rayleigh scattering would have $`P=(1\mathrm{cos}^2\theta )/(1+\mathrm{cos}^2\theta )`$, with $`P=1`$ for $`\theta =90^{}`$. At long wavelengths, the polarization has a distinct peak near $`90^{}`$, but even for $`\lambda 0.9\mathrm{\mu m}`$ the peak polarization is only $`0.45`$. As the wavelength is reduced, the peak polarization declines. For $`\lambda 0.6\mathrm{\mu m}`$, the polarization becomes negative at large scattering angles (see Figure 5), with large negative polarizations in the 120–150 region for $`0.2\mathrm{\mu m}\lambda 0.4\mathrm{\mu m}`$. Note, however, that this negative polarization occurs for scattering angles where the scattering is very weak (see Figure 3) and therefore could be masked by scattering at other points on the sightline where the scattering contributes a positive polarization. Observations of the predicted negative polarization will probably require simple scattering geometries, such as dust in a thin disk, illuminated by a single source.

Similar results are found for the LMC and SMC dust mixtures – see Fig. 5. Note the very large negative polarizations found for the LMC mixture for $`0.2\mathrm{\mu m}\lambda 0.55\mathrm{\mu m}`$. From the variation of the ultraviolet polarization signature between the different grain size distributions in Figure 5 it is apparent that the ultraviolet polarization is sensitive to the details of the grain size distribution.

## 4 Analytic Approximations for the Phase Function

### 4.1 Henyey-Greenstein Phase Function

Scattering of unpolarized incident light by a dust mixture at wavelength $`\lambda `$ is characterized by the total scattering cross section per H nucleon, $`\sigma _{\mathrm{sca}}(\lambda )`$, and a “phase function”

$$\mathrm{\Phi }(\theta ,\lambda )\frac{1}{\sigma _{\mathrm{sca}}(\lambda )}\frac{d\sigma _{\mathrm{sca}}(\theta ,\lambda )}{d\mathrm{\Omega }}$$
(2)
characterizing the angular distribution of the scattered light, where $`d\sigma _{\mathrm{sca}}/d\mathrm{\Omega }`$ is the differential scattering cross section, at scattering angle $`\theta `$, for unpolarized incident light. The definition (2) implies the normalization $`𝑑\mathrm{\Omega }\mathrm{\Phi }=1`$. Isotropic scattering would have $`\mathrm{\Phi }=1/4\pi `$. The first moment

$$\mathrm{cos}\theta =𝑑\mathrm{\Omega }\mathrm{cos}\theta \mathrm{\Phi }(\theta ,\lambda )$$
(3)
of the phase function is a measure of the asymmetry between forward and backward scattering.

Henyey & Greenstein (1941) proposed an analytic function to model anisotropic scattering for dust grain mixtures:

$$\varphi _0(\theta )=\frac{1}{4\pi }\frac{1g^2}{(1+g^22g\mathrm{cos}\theta )^{3/2}},g=\mathrm{cos}\theta .$$
(4)
With a single parameter $`g`$, this is a convenient analytic form that has been widely used to represent dust scattering properties in radiative models of dusty regions.

By construction, $`\varphi _0(\theta )`$ has the correct first moment $`\mathrm{cos}\theta =𝑑\mathrm{\Omega }\mathrm{cos}\theta \varphi _0(\theta )=g`$, but of course $`\varphi _0`$ does not perfectly reproduce the angular dependence of the real phase function $`\mathrm{\Phi }(\theta )`$. As will be seen below, $`\varphi _0`$ is a poor approximation at both $`\lambda 0.4\mathrm{\mu m}`$ and $`\lambda >1\mathrm{\mu m}`$.

### 4.2 A New Phase Function

At long wavelengths, Rayleigh scattering prevails, with $`\mathrm{cos}\theta 0`$ and $`\mathrm{\Phi }(3/16\pi )(1+\mathrm{cos}^2\theta )`$. When $`\mathrm{cos}\theta 0`$, the Henyey-Greenstein phase function $`\varphi _01/4\pi `$; $`\varphi _0`$ is 33% low at $`\theta =0,\pi `$ and 33% high at $`\theta =\pi /2`$.

Consider the phase function

$$\varphi _\alpha (\theta )=\frac{1}{4\pi }\left[\frac{1g^2}{1+\alpha (1+2g^2)/3}\right]\frac{1+\alpha \mathrm{cos}^2\theta }{(1+g^22g\mathrm{cos}\theta )^{3/2}},$$
(5)
with two adjustable parameters, $`\alpha `$ and $`g`$. For $`\alpha =0`$, equation (5) reduces to the Henyey-Greenstein phase function $`\varphi _0`$ (equation 4). For $`g=0`$ and $`\alpha =1`$ this reduces to the phase function for Rayleigh scattering. For $`\alpha =1`$ this corresponds to the phase function proposed by Cornette & Shanks (1992). Analytic results for this phase function are given in Appendices A \- C.

To determine the parameters $`\alpha `$ and $`g`$ we can require that $`\varphi _\alpha (\theta )`$ have the same first and second moments $`\mathrm{cos}\theta `$ and $`\mathrm{cos}^2\theta `$ as $`\mathrm{\Phi }`$. The parameters $`g`$ and $`\alpha `$ are then given by equations (B4-B6). However, although the resulting phase function has correct first and second moments, the fit is poor when the dust is strongly forward-scattering ($`\mathrm{cos}\theta 0.6`$).

We find that an improved fit is obtained if we obtain $`g`$ and $`\alpha `$ from equations (B4-B6) only when the resulting $`\alpha 1`$; for values of $`\mathrm{cos}\theta `$ and $`\mathrm{cos}^2\theta `$ for which equations (B4-B6) lead to $`\alpha >1`$, we instead set $`\alpha =1`$ and obtain $`g`$ from equation (C2). We will refer to this new analytic phase function as $`\varphi _{\alpha 1}`$.

Figure 6 shows the scattering phase function $`\mathrm{\Phi }(\theta )`$ for the WD01 Milky Way dust model, together with the Henyey-Greenstein phase function $`\varphi _0`$ and our new phase function $`\varphi _{\alpha 1}`$, at wavelengths $`\lambda =1.2\mathrm{\mu m}`$, $`0.80\mathrm{\mu m}`$, and $`0.40\mathrm{\mu m}`$. For these three wavelengths, the new phase function provides an improved fit to the actual phase function $`\mathrm{\Phi }`$. For $`\lambda =0.40\mathrm{\mu m}`$, we also show $`\mathrm{\Phi }`$ approximated by a sum over Legendre polynomials $`P_l`$ up to $`l=8`$. This 9 term expansion does not reproduce $`\mathrm{\Phi }`$ as well as the new phase function $`\varphi _{\alpha 1}`$.

At ultraviolet wavelengths the grains become more forward-throwing, and the analytic phase functions $`\varphi _0`$ and $`\varphi _{\alpha 1}`$ no longer provide a good fit, as seen in Figure 7 for $`\lambda =0.20\mathrm{\mu m}`$ and $`0.10\mathrm{\mu m}`$. At $`\lambda =0.20\mathrm{\mu m}`$, for example, both the Henyey-Greenstein phase function $`\varphi _0`$ and the new phase function $`\varphi _{\alpha 1}`$ underestimate the forward scattering intensity by a factor of $``$2 for $`\mathrm{cos}\theta 0.98`$ ($`\theta 10^{}`$).

To quantify the error associated with using an analytic phase function $`\varphi (\theta )`$ to approximate an actual phase function $`\mathrm{\Phi }(\theta )`$, we define the r.m.s. relative error

$$h_{\mathrm{rel}}\left[\frac{d\mathrm{\Omega }}{4\pi }\left[\frac{\varphi (\theta )\mathrm{\Phi }(\theta )}{\mathrm{\Phi }(\theta )}\right]^2\right]^{1/2},$$
(6)
and r.m.s. absolute error

$$h_{\mathrm{abs}}\left[\frac{d\mathrm{\Omega }}{4\pi }\left[\frac{\varphi (\theta )\mathrm{\Phi }(\theta )}{\mathrm{\Phi }}\right]^2\right]^{1/2}.$$
(7)
where, of course, $`\mathrm{\Phi }=1/4\pi `$. The r.m.s. relative error $`h_{\mathrm{rel}}`$ would appear to be the best measure of the overall quality of fit. $`h_{\mathrm{rel}}`$ and $`h_{\mathrm{abs}}`$ will differ substantially only when the phase function is very asymmetric; large fractional errors in directions where the scattering is weak then make only a small contribution to $`h_{\mathrm{abs}}`$, while modest fractional errors in directions where the scattering is very strong will make large contributions to $`h_{\mathrm{abs}}`$.

In the limit of Rayleigh scattering ($`\lambda a`$), $`\mathrm{\Phi }(\theta )\frac{3}{4}\left(1+\mathrm{cos}^2\theta \right)`$, and $`\varphi _0`$ has errors

$`h_{\mathrm{abs}}`$ $``$ $`1/\sqrt{20}0.22,`$ (8)
$`h_{\mathrm{rel}}`$ $``$ $`(134\pi )^{1/2}/30.22.`$ (9)
Figure 8 compares the error of different analytic approximations to the phase function, using as an example $`\mathrm{\Phi }(\theta ,\lambda )`$ for the WD01 model for $`R_V=3.1`$ Milky Way dust. The errors $`h_{\mathrm{rel}}`$ and $`h_{\mathrm{abs}}`$ are shown in Figure 8 for the Henyey-Greenstein function $`\varphi _0`$, the Cornette-Shanks phase function $`\varphi _1`$, and our new phase function $`\varphi _{\alpha 1}`$. At long wavelengths $`\lambda 3\mathrm{\mu m}`$, the Henyey-Greenstein approximation $`\varphi _0`$ has $`h_{\mathrm{rel}}h_{\mathrm{abs}}0.22`$, as expected from eq. (8,9). Between $`0.5`$ and $`1\mathrm{\mu m}`$, the relative errors for $`\varphi _0`$ are modest, $`h_{\mathrm{rel}}<10\%`$, so $`\varphi _0`$ provides a good approximation to the actual scattering properties. However, $`\varphi _0`$ has $`h_{\mathrm{rel}}>10\%`$ for $`\lambda <0.5\mathrm{\mu m}`$, rising to shorter wavelengths.

Figure 8 shows the fractional errors $`h_{\mathrm{rel}}`$ and $`h_{\mathrm{abs}}`$ for the new phase function $`\varphi _{\alpha 1}`$. As expected, the new phase function provides an excellent fit at long wavelengths, with fractional errors $`<3\%`$ for $`\lambda >1\mathrm{\mu m}`$. As the wavelength is decreased, the new phase function $`\varphi _{\alpha 1}`$ (as measured by $`h_{\mathrm{rel}}`$) remains preferable to the Henyey-Greenstein phase function for $`\lambda 0.19\mathrm{\mu m}`$. $`\varphi _{\alpha 1}`$ has $`h_{\mathrm{rel}}<10\%`$ for $`\lambda >0.27\mathrm{\mu m}`$.

In the ultraviolet the phase function becomes very strongly forward-scattering, and all of the the analytic approximations have significant errors. More complicated parameterizations – such as the use of multiple Henyey-Greenstein phase functions, as by Witt (1977) and Hong (1985) – may be considered. Alternatively, a given phase function $`\mathrm{\Phi }(\theta )`$ can be approximated by summing over Legendre polynomials, although the sum needs to include a large number of terms to provide a good approximation – 16 terms are required for $`\lambda =0.20\mathrm{\mu m}`$, and even 51 terms are insufficient for $`\lambda =0.10\mathrm{\mu m}`$. At short wavelengths it may be best to simply tabulate $`\mathrm{\Phi }(\theta )`$ and interpolate as required.

## 5 Discussion

The scattering properties of interstellar dust have been determined observationally by comparing the observed surface brightness of reflection nebulae with model nebulae computed with different dust scattering properties, selecting the model which provides the best match to the observations. The usual approach has been to try to determine only two grain properties – the albedo and $`\mathrm{cos}\theta `$ – by finding radiative transfer models which appear to be consistent with the observed surface brightness of scattered light.

It has been customary to assume that the scattering phase function $`\mathrm{\Phi }(\theta )`$ can be approximated by the Henyey-Greenstein function $`\varphi _0(\theta )`$, with $`g=\mathrm{cos}\theta `$

This approach has been used with (1) observations of high surface brightness reflection nebulae (e.g., NGC 2023, NGC 7023, IC 435) illuminated primarily by a single star, (2) observations of individual dust clouds illuminated by ambient starlight, and (3) observations of the much fainter diffuse galactic light (the entire Galaxy as a reflection nebulosity). Figures 9 and 10 show the results so obtained from a number of independent studies, at wavelengths from the optical to vacuum ultraviolet. Also shown are the albedo and asymmetry factor calculated for WD01 grain models representing average Milky Way dust with $`R_VA_V/E(BV)=3.1`$, and dust from denser regions with $`R_V=4.0`$ and 5.5. While there are differences among the different observational determinations, particularly in the ultraviolet, many of the observational results are in general agreement with the albedo and asymmetry factor calculated for the WD01 grain model.

For $`\lambda ^1<4\mathrm{\mu m}^1`$, the observational studies appear to be in general agreement with one another, but at shorter wavelengths the observational results are sometimes in conflict. For example, from observations of the diffuse galactic light at $`6.25\mathrm{\mu m}^1`$ ($`\lambda =0.16\mathrm{\mu m}`$) Hurwitz et al. (1991) found $`(a,g)=(0.185\pm .055,0.2\pm 0.2)`$, while Witt et al. (1997) found $`(a,g)=(0.45\pm 0.05,0.68\pm 0.10)`$.

Aside from genuine regional variations in dust properties, there are many factors which could contribute to such discrepancies. Foremost may simply be the difficulty of observations in the vacuum ultraviolet, followed by uncertainties concerning the nebular geometry and the illuminating radiation. However, given that use of the Henyey-Greenstein phase function $`\varphi _0`$ introduces an absolute error $`h_{\mathrm{abs}}70\%`$ in the scattering phase function at $`\lambda <0.16\mathrm{\mu m}`$ (see Fig. 8), it seems possible that reliance on the Henyey-Greenstein phase function may contribute to the disagreements among different observational determinations of the albedo and scattering asymmetry in the ultraviolet.

Unfortunately, the new phase function $`\varphi _{\alpha 1}`$ is inaccurate for $`\lambda <0.27\mathrm{\mu m}`$. For $`\lambda ^14\mathrm{\mu m}^1`$, radiative transfer models should use a more accurate representation of the phase function $`\mathrm{\Phi }`$. As discussed above, one possibility is expansion in Legendre polynomials, although many terms are required. Alternatively, it may be simplest to use tabulated values of $`\mathrm{\Phi }`$ with interpolation. “Inversion” of the observational data to infer the grain properties may not be feasible; the best approach may be to assume the dust properties \[i.e., $`\sigma _{\mathrm{ext}}(\lambda )`$, albedo($`\lambda `$), and $`\mathrm{\Phi }(\theta ,\lambda )`$\] based on a grain model, and then seek a dust spatial distribution which maximizes consistency with the observed surface brightness and whatever other constraints are available. Failure to find agreement would suggest that the grain model may be incorrect.

Applying this approach to reflection nebulae is complicated by the possibility of clumpy dust distributions, which Mathis et al. (2002) show can produce a range of ratios of nebular surface brightness to unscattered stellar flux, depending on the detailed dust clump distribution and the viewing direction. Diffuse clouds (i.e., the diffuse galactic light) and externally-illuminated clouds may prove easier to interpret than bright reflection nebulae, which generally have a bright star embedded in the dust distribution.

## 6 Summary

The following are the principal results of this work:

1. Differential scattering cross sections have been calculated at selected wavelengths from infrared to ultraviolet for dust mixtures appropriate to the diffuse interstellar medium in the Milky Way, the LMC, and the SMC. These scattering functions, which can be used for modeling reflection nebulae, are available at http://www.astro.princeton.edu/$``$draine .
2. The polarization of scattered light, as a function of scattering angle, is calculated at selected wavelengths.
3. The calculated phase function $`\mathrm{\Phi }(\theta )`$ for the WD01 model for $`R_V=3.1`$ MW dust is compared to the widely-used Henyey-Greenstein phase function $`\varphi _0(\theta )`$. At wavelengths $`0.48\mathrm{\mu m}<\lambda <0.96\mathrm{\mu m}`$, $`\varphi _0`$ provides a good fit to $`\mathrm{\Phi }`$, with r.m.s. relative error $`h_{\mathrm{rel}}<10\%`$. In the ultraviolet, however, $`\varphi _0`$ provides a poor fit, with r.m.s. absolute error $`h_{\mathrm{abs}}>50\%`$ at $`\lambda <2400\mathrm{\AA }`$. The Henyey-Greenstein phase function is not suitable for accurate modeling of reflection nebulae in the ultraviolet.
4. A new phase function, $`\varphi _{\alpha 1}(\theta )`$ is proposed \[see equation (5)\]. This phase function provides a good approximation to our calculated phase functions for $`\lambda >0.27\mathrm{\mu m}`$, with r.m.s. relative error $`h_{\mathrm{rel}}<10\%`$. At shorter wavelengths, however, $`\mathrm{\Phi }`$ becomes strongly forward-throwing, and neither the Henyey-Greeenstein phase function $`\varphi _0`$ nor the new phase function $`\varphi _{\alpha 1}`$ provide good approximations.
5. There are significant discrepancies among observational determinations of the albedo and $`\mathrm{cos}\theta `$ in the ultraviolet for the diffuse galactic light and dust in discrete clouds and reflection nebulae. One possible cause for these discrepancies may be reliance on the Henyey-Greenstein phase function in the radiative transfer models. It is recommended that the new phase function $`\varphi _{\alpha 1}`$ be used at wavelengths $`\lambda 0.27\mathrm{\mu m}`$; at shorter wavelengths accurate radiative transfer models should use some other procedure to obtain $`\mathrm{\Phi }`$.

I thank Karl Gordon, Aigen Li, Kalevi Mattila, Jonathan Tan, Adolf Witt, and Michael Woolf for valuable comments, the anonymous referee for helpful suggestions, and Robert Lupton for making available the SM software package. This work was supported in part by NSF grant AST-9988126.

## Appendix A Phase Function Moments

A phase function of the form

$$\varphi _\alpha (\theta )=\frac{1}{4\pi }\left[\frac{1g^2}{1+\alpha (1+2g^2)/3}\right]\frac{1+\alpha \mathrm{cos}^2\theta }{(1+g^22g\mathrm{cos}\theta )^{3/2}}$$
(A1)
has first and second moments given by

$`\mathrm{cos}\theta `$ $`=`$ $`g{\displaystyle \frac{1+\alpha (3+2g^2)/5}{1+\alpha (1+2g^2)/3}},`$ (A2)
$`\mathrm{cos}^2\theta `$ $`=`$ $`{\displaystyle \frac{1+2g^2+(3\alpha /35)\left(7+20g^2+8g^4\right)}{3+\alpha (1+2g^2)}}.`$ (A3)
## Appendix B Reproduce $`\mathrm{cos}\theta `$ and $`\mathrm{cos}^2\theta `$

If we specify $`\mathrm{cos}\theta `$ and $`\mathrm{cos}^2\theta `$, the adjustable parameters $`g`$ and $`\alpha `$ are determined from equations (A2,A3). It can be shown that $`g`$ satisfies the cubic equation

$$g^3\frac{17}{3}\mathrm{cos}\theta g^2+7\mathrm{cos}^2\theta g\frac{7}{3}\mathrm{cos}\theta =0.$$
(B1)
Let

$`a`$ $``$ $`{\displaystyle \frac{7}{3}}\mathrm{cos}^2\theta {\displaystyle \frac{289}{81}}\mathrm{cos}\theta ^2,`$ (B2)
$`b`$ $``$ $`{\displaystyle \frac{119}{18}}\mathrm{cos}\theta \mathrm{cos}^2\theta {\displaystyle \frac{4913}{729}}\mathrm{cos}\theta ^3{\displaystyle \frac{7}{6}}\mathrm{cos}\theta .`$ (B3)
If $`a^3+b^2>0`$, the solution is

$$g=\left[\left(a^3+b^2\right)^{1/2}b\right]^{1/3}\left[\left(a^3+b^2\right)^{1/2}+b\right]^{1/3}+\frac{17}{9}\mathrm{cos}\theta ,$$
(B4)
whereas if $`a^3+b^2<0`$, the appropriate root is

$$g=2|a|^{1/2}\mathrm{cos}(\psi /3)+\frac{17}{9}\mathrm{cos}\theta ,\psi \mathrm{arccos}(b/|a|^{3/2}).$$
(B5)
The parameter $`\alpha `$ is then obtained from equation (A2):

$$\alpha =\frac{25(\mathrm{cos}\theta g)}{3(3+2g^2)g5(1+2g^2)\mathrm{cos}\theta }.$$
(B6)
## Appendix C Fix $`\alpha `$, Reproduce $`\mathrm{cos}\theta `$

If we choose to fix the value of $`\alpha `$ and $`\mathrm{cos}\theta `$, then $`g`$ must satisfy

$$g^3\frac{5}{3}\mathrm{cos}\theta g^2+\left(\frac{3\alpha +5}{2\alpha }\right)g\frac{5}{6\alpha }(3+\alpha )\mathrm{cos}\theta =0.$$
(C1)
For $`\alpha 0`$, equation (C1) has the real root

$`g`$ $`=`$ $`\left[(a_\alpha ^3+b_\alpha ^2)^{1/2}+b_\alpha \right]^{1/3}\left[(a_\alpha ^3+b_\alpha ^2)^{1/2}b_\alpha \right]^{1/3}+{\displaystyle \frac{5}{9}}\mathrm{cos}\theta ,`$ (C2)
$`a_\alpha `$ $``$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{5}{6\alpha }}{\displaystyle \frac{25}{81}}\mathrm{cos}\theta ^2,`$ (C3)
$`b_\alpha `$ $``$ $`{\displaystyle \frac{125}{729}}\mathrm{cos}\theta ^3+{\displaystyle \frac{5}{9\alpha }}\mathrm{cos}\theta .`$ (C4)
For $`\alpha 0`$ (the Henyey-Greenstein phase function) we see immediately that equation (C1) requires $`g=\mathrm{cos}\theta `$.


text:
# 1. Introduction

## 1. Introduction

The last few years have brought rapid advances in the understanding and application of quantum entangled photons. Exploitation of quantum correlations is important in various applications, including teleportation , quantum lithography , clock synchronization , and optical communications . To date, most of this work has focused on characterizing quantum correlation effects at low intensities due to the lack of bright sources of correlated photon pairs. However, Lamas-Linares et al.’s  recent report of the laser-like production of polarization entangled photons now opens up the possibility that one can take advantage of correlated effects for communications. In addition, the demonstration of efficient parametric down conversion using a diode laser and a pair of solid state parametric conversion chips , demonstrates that parametric generation technology is rapidly approaching the fabrication thresholds where economies of scale can be achieved. Earlier attempts, by Mandel , to apply correlated photons to improving the signal to noise ratio in a communications link depends on coincidence counters, which, because of the detector dead time, are limited to low rates. Our approach depends on the development of a special detector that is intrinsically sensitive to two-photons . This single detector replaces the coincidence counters, gating electronics, amplifiers, and the computer interface employed in that earlier effort. Consequently, it is important to do a link analysis to determine whether there is an advantage to using correlated photons and two photon detectors for conveying information in an optical communications channel.

The answer is yes, but only for the special situation where the signal strength is swamped by in-band background noise. We refer specifically to situations where free space high bandwidth of optical communications are desired, but in-band solar background radiation hampers the signal to noise. The other situation of interest would be in fiber networks operating under conditions of a large background of incoherent in-band scattering, which adversely impacts recovery of the signal at the receiver. For shot noise limited telecom links there is generally no advantage to using correlated photons as opposed to classical photons for the carrier. Our approach is predicated on the development of a two photon sensitive photodetector which eliminates the need for separate coincidence electronics. Consequently, one can define the required figure of merit that must be met for a performance advantage to be realzied of the corrrelated 2-photon link over a classical uncorrelated photon link in terms of the two photon and single photon quantum efficiencies of the detectors. Besides introducing a new method for performing the coincidence measurement, this paper defines the conditions under which one can expect to the two-photon correlated detection to realize and advantage when compared to classical uncorrelated signal recovery. Finally, it briefly examines the effect of diffraction on the 2-photon collection efficiency in the free space application.

## 2. Near-Field Optical Communications Links

The near field communication link is defined as an architecture for which the full output power of the transmitter is subtended and collected by the detector. The schematics in Figures 1(a) and 1(b) depict the key differences between the classical and the correlated-photon near field telecom links. In this analysis, a classical link, using photons of frequency, $`\omega _1`$, is compared with a correlated-photon link transmitting photon pairs of approximate frequency, $`\omega _1`$, but actual frequency $`\omega _a+\omega _b`$ = $`2\omega _1`$ = $`\omega _2`$. The power received at the detector is given by the following expression

$$P_r(\omega _i,t)=\mu P_t(\omega _i,t)L(\omega _i).$$
(1)
where $`P_t(\omega _i,t)`$ is the modulated output power of the communications laser at carrier frequency $`\omega _i`$, $`\mu `$ represents conversion losses suffered after emission from the transmitter, and L represents any losses that occur during transit through the communications channel. This expression assumes that the collimating optics at the transmitter and the collecting optics at the receiver are sufficiently large that the footprint of the transmitted beam is completely subtended by the receiver aperture. This would be a reasonable assumption at optical frequencies over free space transmission ranges up to 1000 km. To facilitate comparison of the performance of both links, we assume that the initial laser output power of the classical link is identical to the laser output power of the correlated photon link,

$$P_t(\omega _1,t)=P_t(\omega _2,t).$$
(2)
The link designs differ by the conversion loss factor, $`\mu `$, and the carrier frequency used to transmit the signal information through the channel. In the classical telecom link, $`\mu =1`$ and the carrier frequency is $`\omega _1`$. Thus the signal current generated at the receiver for the classical link is given by:

$$I_0=\eta _{det}P_r(\omega _1,t)=\eta _{det}P_t(\omega _1,t)L(\omega _1),$$
(3)
where $`\eta _{det}`$ is the receiver efficiency for converting input photons to carriers in the signal current.

In contrast, the value of $`\mu `$ in the correlated photon link is determined by the method used to produce the correlation. Furthermore, we assume that the modulated output of the laser transmitter is passed through a nonlinear parametric down conversion crystal and correlated photons are generated with efficiency, $`\mu =\eta _{PDC}`$. The following additional properties apply.

* Both photons were created at the same point in time.
* Both photons were created at the same spatial point.
* Energy is conserved $`\omega _2=\omega _a+\omega _b`$.
* Momentum is conserved $`\stackrel{}{k_2}=\stackrel{}{k_a}+\stackrel{}{k_b}`$.

where one assumes that the daughter photons from the down conversion process are degenerate with $`\omega _a+\omega _b=2\omega _1=\omega _2`$. After down conversion, the output is modulated with the signal waveform before being sent through the transmission medium for collection at the receiver. Therefore in the correlated link, the signal current generated at the receiver is proportional to the signal power and is given by

$$I_0^{}=\eta _{2ph}(2\omega _1)[P_t^{}(\omega _a,t)L(\omega _a,t)+P_t^{}(\omega _b,t)L(\omega _b)]=\eta _{2ph}(2\omega _1)\eta _{PDC}P_t(\omega _2,t)L(\omega _1).$$
(4)
Even though we do not know exactly when or where any pair of twin photons are born within the down conversion crystal, the fact that they are simultaneously created means that standard geometrical imaging optics can be exploited in a straightforward matter to reunite these photon pairs in coincidence at the receiver. The figure of merit for comparing the ultimate link performance is given by the signal to noise ratio of each approach;

$$SNR_{classical}=\frac{I_0}{\sqrt{_i\sigma _i^2}},$$
(5)
$$SNR_{correlated}=\frac{I_0^{}}{\sqrt{_i\sigma _i^2}},$$
(6)
where

$$\underset{i}{}\sigma _i^2=\sigma _{thermal}^2+\sigma _{shot}^2+\sigma _{laserRIN}^2+\sigma _{background}^2$$
(7)
is the sum of the variances of all noise contributions, thermal noise, shot noise, relative intensity noise due to the laser, and in-band background noise arriving with the signal at the detector. These factors are defined for each link in Table 1 below. The expression for the correlated SNR is optimized when the source has been imaged onto the detector so that the correlated pairs, that are collected by the imaging lens, arrive at the same time and overlap at the same point of the detector. This means that correlated telecom links are practical for short range links where the foot print of the transmitter s output beam are fully subtended by the receiver collection aperture.

Here $`k`$ is Boltzman s constant, $`B`$ is the circuit bandwidth, $`T`$ is the ambient absolute temperature, $`R_i`$ is the detector input impedance, $`e`$ is the electric charge, $`F_{RIN}`$ is the relative intensity noise factor, and $`P_{rB}`$ is the power level of the background noise at the receiver. To compare the relative efficiencies of classical signal recovery with the signal recovery process of correlated photons, we assume similar links are established for the classical photon signal recovery and correlated photon signal recovery. $`P_t`$, $`A_t`$, $`R_i`$, $`L`$, $`A_r`$, and $`\omega _1`$ are chosen to be identical for both links, where $`A_t`$ is the area of he transmiter aperture and $`A_r`$ is the area of the receiver collection aperture. They differ in the choice of down conversion efficiencies, detector efficiencies, and the background noise contribution to the signal to noise ratio (SNR).

The differences between these two types of links are summarized a follows. The correlated photon transmitter design can never be as efficient as the classical transmitter due to the parametric down conversion factor, $`\eta _{PDC}`$ $`<1`$. But an advantage can be obtained if one employs a detector that is specifically designed to have a poor single photon detection efficiency, and a two-photon detection efficiency that is significantly larger such that $`\eta _{2ph}(2\omega _1)`$ $`>>`$ $`\eta _{1ph}(\omega _1)`$. To illustrate this point, we will compare the correlated photon signal to noise to the classical signal to noise. The ideal communications link is often designed so that the shot noise is dominant. Which yields,

$$\frac{SNR_{correlated}}{SNR_{classical}}=\sqrt{\frac{\eta _{2ph}(2\omega _1)\eta _{PDC}L(\omega _1)}{\eta _{det}}}<1.$$
(8)
Clearly the ratio in the expression above will always be less than 1, largely because the down conversion efficiency will never be 100$`\%`$. Alternatively, if one examines the case where the background noise exceeds the signal current, $`I_B`$ $`>`$ $`I_0`$, the ratios change to

$$\frac{SNR_{correlated}}{SNR_{classical}}=\frac{\eta _{2ph}(2\omega _1)\eta _{PDC}L(\omega _1)}{\sqrt{\eta _{det}\eta _{1ph}(\omega _1)}}>1.$$
(9)
Here we assume that $`\eta _{accident}(2\omega _1)P_{rB}`$ shown in Table 1 is negligible compared to $`\eta _{1ph}(\omega _1)`$. Equation 9 says that the correlated photon link has a performance advantage as long as $`\eta _{2ph}(2\omega _1)`$ $`>>`$ $`\eta _{1ph}(\omega _1)`$, in accordance with this relation.

## 3. Conclusions

Certain communications channels need to operate in the presence of severe background noise, such as when the Sun is positioned near the line of sight (LOS) between the transmitter and the receiver. In such situations noise-immune coding techniques are of limited help due to saturation of the detector. Narrow-band filters can limit the background but with a sufficiently intense source the detectors can still be saturated by noise. In this situation two-photon correlated detection can avoid the noise in an entirely different way, where a two-photon detector does not see it, or only sees “accidental” coincidences which are small for incoherent sources like the Sun. If this is the case a quantum-correlated communication channel can out-perform other techniques because it eliminates the background before a detector signal is generated.

In telecom links for which the footprint of the of the beam is larger than the collection aperture, quantum-correlated telecom links are essentially non competitive against classical communications links because the correlation cross section of the product, $`P_t(\omega _a,t)P_t(\omega _b,t)`$, falls off at a rate proportional to $`1/R^4`$ while the classical link has a $`1/R^2`$ dependence. However, over ranges of about 1000 km realizable communications links using 1 to 2-meter mirror optics permit the design of links that collect most of the emitted photons, thereby making free space links under this range superior if the conditions of Equation 9 are met. In addition, fiber optic links, which intrinsically permit one to image all the incident photons onto the receiver, may find an advantage in using correlated photons when the in-band background noise exceeds the signal level.

In summary, one will observe a performance advantage for correlated photon links when the classical link is limited by an intense source of uncorrelated background noise. The requirements for using this technique are (1) that an intrinsically two-photon detector be developed (this rejects noise without being saturated by it) and (2) that most of the transmitted photons be collected (which can be done in free space over distances up to about 1000km).

## Acknowledgements

The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, and was supported by a contract with the Office of Naval Research.

Figure 1, D.J. Jackson et. al.

$``$ $`\eta _{1ph}(2\omega _1)`$ is the single photon detection efficiency.

$``$ $`\eta _{accidental}(2\omega _1)`$ is the two photon absorption efficiency for statistically random detection events.