text:
# Galaxy spectral parameterisation in the 2dF Galaxy Redshift Survey

## 1 Introduction

Galaxy redshift surveys are now probing the galaxy distribution of the local Universe more accurately than ever before, and in so doing they are establishing many fundamental properties of the galaxy population and its large-scale structure. The 2dF Galaxy Redshift Survey (2dFGRS) is one such ambitious project, conceived with the aim of mapping the galaxy distribution to an extinction corrected magnitude limit of $`b_\mathrm{J}19.5`$. This particular survey is now essentially complete with approximately 230,000 galaxy spectra having already been acquired, and has already started to yield significant results (e.g. Percival et al. 2001; Madgwick et al. 2002 and others). In addition to this survey, the Sloan Digital Sky Survey (Strauss et al., 2002) is also underway, and once complete will obtain up to 1,000,000 individual galaxy spectra.

Apart from the main scientific goals of quantifying the large-scale structure of the Universe, one of the most significant contributions of galaxy redshift surveys is to our understanding of the galaxy population itself, through the information about galaxy properties contained in the observed spectra. Having a data set of 230,000 galaxy spectra — as in the case of the 2dFGRS — allows us to test the validity of galaxy formation and evolution scenarios with unprecedented accuracy. However, the sheer size of the data set presents its own unique problems. Clearly, in order to make the spectral data set more ‘digestible’, some form of data compression is necessary. Familiar statistics such as equivalent width measurements, morphological types, and broad-band colours are really just compression techniques in some sense. These quantities can be compared with theoretical predictions and simulations, and hence can set constraints on scenarios for galaxy formation and biasing.

The approach that has been adopted in Madgwick et al. (2002) is to define a spectral indicator, $`\eta `$, based on Principal Component Analysis (PCA). Because of instrumental limitations, the flux calibration of the 2dF spectra is unreliable, and hence robust measurements of the shape of the spectral continua are not possible. For this reason, the $`\eta `$ statistic was designed to (in a loose sense) preserve the most discriminatory information in the spectra, while being robust to the instrumental uncertainties. Effectively, the parameter $`\eta `$ measures the strength of absorption and emission lines, while remaining insensitive to the slope of the continuum (or broad band colour).

The motivation of the definition of the $`\eta `$ statistic was in part pragmatic. The goal of this paper is to address the fundamental question of the *physical significance* of this parameter, in terms of how it relates to the physical properties of galaxies. A substantial amount of work has been carried out in the past on relating the observed spectra of galaxies to the physical processes occurring in them (see e.g Kennicutt, Tamblyn & Congdon 1994; Ronen, Arágon-Salamanca & Lahav 1999; Carter, Fabricant, Geller & Kurtz 2001 and references therein), and the advent of a large and uniform survey such as the 2dFGRS will result in many more advances in this field. However, for the present analysis we restrict ourselves to relating only the 2dFGRS spectral classification parameter, $`\eta `$, to these processes and leave a fuller development of this subject for future works. The investigation presented here is particularly timely as many new results from the 2dFGRS are being presented which make exclusive use of this parameter to characterise and partition the galaxy population (e.g. Madgwick et al. 2002; Norberg et al. 2002; Martinez et al. 2002).

We make use of two kinds of models to pursue this question. In the first approach, we adopt a simple parameterisation to describe the star formation history of each galaxy in terms of an exponentially decreasing function of time. For a given redshift of observation, the star formation history is then characterised by a single parameter: the timescale of the exponential decline. These star formation histories are then convolved with spectral synthesis models, using the PEGASE (Projet d’Etude des GAlaxies par Synthese Evolutive) code, developed by Fioc & Rocca-Volmerange (1997). A particular advantage of this package compared to others available in the literature (e.g. Bruzual & Charlot 1993) is that it is possible to include modelling of nebular emission from star forming regions in each galaxy — an important ingredient in the calculation of the $`\eta `$ parameter. We refer to these models as the “simple” models, because of the rather simplified nature of the parametrised star formation histories.

The limitation of this simple method is that we must explicitly assume some family of star-formation histories for our galaxy population. We therefore extend this approach to incorporate the cosmological framework of the Cold Dark Matter (CDM) scenario, by using semi-analytic techniques (see e.g. Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; Somerville & Primack 1999). We make use of a mock catalogue of synthetic spectra which was designed to match the selection criteria of the 2dFGRS, analysed previously in Slonim et al. (2001), and created using an updated version of the code described in Somerville & Primack (1999). There are several advantages of using a semi-analytic model to create spectra for this exercise. The star formation history of each galaxy will reflect at least some of the complexity of the interconnected processes of dark matter clustering, gas cooling, galaxy merging, supernova feedback, etc., modelled self-consistently within a specific cosmological framework. In addition, as we can select galaxies so as to reproduce the same redshift and luminosity distributions displayed by the observed 2dF galaxies, the resulting ensemble should contain a mix of galaxies with different sorts of star formation histories that is similar to the actual observed sample. Finally, the semi-analytic model also yields many other physical parameters for each galaxy, giving us the potential to further probe the physical processes that may be related to $`\eta `$.

We find that in the simple models, there is essentially a one-to-one relationship between $`\eta `$ and the ratio of the present rate of star formation to its past time-averaged value. This quantity is sometimes known as the “birthrate” parameter, $`b`$ (Scalo 1986). Even more surprisingly, we find that there is a strong correlation between $`\eta `$ and $`b`$ even in the spectra produced by the semi-analytic models, despite the complex and diverse nature of the underlying star formation histories. Empirically, the $`b`$ parameter is known to be tightly correlated with morphological type and colour. We therefore argue that the $`\eta `$ parameter presents a practical and robust means of extracting a fundamental measure of galaxy type from the relatively low-quality spectra typical of present-day large redshift surveys.

Section 2 of this paper briefly reviews the operational definition of and motivation behind the spectral classification parameter $`\eta `$, used in the 2dFGRS. In Section 3 we make use of the “simple” models based on exponentially declining star formation histories to show how $`\eta `$ is related to the star-formation history of a galaxy. Section 4 then generalises these results using formation histories generated with a semi-analytic model. We conclude with a discussion of our results and ways to build upon the results presented here.

## 2 Spectral Classification in the 2dFGRS

### 2.1 Principal Component Analysis

Principal component analysis (PCA) is a well established statistical technique which has proved very useful in dealing with high dimensional data sets (see e.g. Murtagh & Heck 1987; Connolly et al. 1995). In the particular case of galaxy spectra we are typically presented with approximately 1000 spectral channels per galaxy, however when used in applications this is usually compressed down to just a few numbers, either by integrating over small line features — yielding equivalent widths — or over broad band filters. The key advantage of using PCA in our data compression is that it allows us to make use of all the information contained in the spectrum in a statistically unbiased way, i.e. without the use of such ad hoc filters.

In order to perform the PCA on our galaxy spectra we first construct a representative volume limited sample of the galaxies. We compute the “mean spectrum” from this ensemble, and subtract this mean from all of the galaxy spectra. When we apply the PCA to this sample an orthogonal set of components (eigenspectra) is constructed, which span the wavelength space occupied by the galaxy spectra. These components have been specifically chosen by the PCA in such a way that as much information (variance) is contained in the first eigenspectrum as possible, and that the amount of the remaining information in all the subsequent eigenspectra is likewise maximised. Therefore, if the information contained in the first $`n`$ eigenspectra is found to be significantly greater than that in the remaining eigenspectra we can significantly compress the data set by swapping each galaxy spectrum (described by 1000 channels) with its projections onto just those $`n`$ eigenspectra.

In the case of the 2dF galaxy spectra we find that approximately two thirds of the total variance (including the noise) in the spectra can be represented in terms of only the first two projections ($`pc_1`$, $`pc_2`$). So, at least to a first approximation, galaxy spectra can be thought of as a two dimensional sequence in terms of these two projections.

In Fig. 1 we show these first two eigenspectra. It can be seen from this figure that whilst the first eigenspectrum contains information from both lines and the continuum, the second is dominated by absorption and emission lines. Because of this it is possible to take two simple linear combinations which isolate either the continuum or the emission/absorption line features. In effect what we are doing when we utilise these linear combinations is rotating the axes defined by the PCA to make the interpretation of the components more straightforward. In so doing we can see that a parameterisation in terms of $`pc_1`$ and $`pc_2`$ is essentially equivalent to a two dimensional sequence in colour (continuum slope) and the average emission/absorption line strength. This ability to isolate the continuum and line components of the galaxy spectra turns out to be very useful, as will be discussed in Section 2.2.

### 2.2 The $`\eta `$ Parameter

The 2dF instrument makes use of up to 400 optical fibres with a diameter of 140$`\mu `$m (corresponding to $`2.0^{\prime \prime }`$ on the sky, depending on plate position. See Lewis et al. 2002). The instrument itself was designed to measure large numbers of redshifts in as short an observing time as possible. However, in order to optimise the number of redshifts that can be measured in a given period of time, compromises had to be made with respect to the spectral quality of the observations. Therefore if one wishes to characterise the observed galaxy population in terms of their spectral properties, care must be taken in order to ensure that these properties are robust to the instrumental uncertainties.

The quality and representativeness of the observed spectra can be compromised in several ways and a more detailed discussion of these issues is presented in previous work (see e.g. Madgwick et al. 2002). The net effect is that the uncertainties introduced into the fibre-spectra predominantly affect the calibration of the continuum slope and have relatively little impact on the emission/absorption line strengths. For this reason any given galaxy spectrum which is projected onto the plane defined by ($`pc_1`$,$`pc_2`$) will not be uniquely defined in the direction of varying continuum but will be robust in the orthogonal direction (which measures the average line strength).

The linear combination of the first two eigenspectra which is robust to these uncertainties is shown in Fig. 1 (b) and denoted by $`\eta `$ (ETA\_TYPE in the 2dFGRS catalogue111http://www.mso.anu.edu/2dFGRS/). It is simply

$$\eta =apc_1pc_2,$$
(1)
where $`a`$ is a constant which we find empirically to be $`a=0.5\pm 0.1`$.

We have now identified the single statistically dominant component of the galaxy spectra which is robust to the known instrumental uncertainties. We have therefore chosen to adopt this (continuous) variable as our measure of spectral type. Having defined $`\eta `$ in this formal and pragmatic manner, however, we are left with the question of whether it is physically meaningful, and how it is to be interpreted. This is the issue that we now address.

### 2.3 Correlation of $`\eta `$ with traditional galaxy classifiers

We show the observed distribution of $`\eta `$ for the 2dF galaxies in Fig. 2. It is rather intriguing that this distribution is strongly bimodal. Fig. 2 also shows the $`\eta `$ projections of spectra from galaxies in the Kennicutt Atlas (Kennicutt 1992), which have known morphologies. There is a correspondence between morphology and $`\eta `$ in the expected sense: the value of $`\eta `$ increases as one moves towards later type objects in the Hubble sequence. Larger values of $`\eta `$ indicate more dominant emission lines, and it is well-known that later type galaxies have stronger emission lines. From this diagram we see that $`\eta `$ seems to be correlated with morphological type, at least based on this small sample of galaxies. Norberg et al. (2002) showed the distribution of $`\eta `$ for a larger sample of galaxies that had been classified into four morphological types (E, S0, Sp and Irr), and showed that while there is a correlation between morphological type and $`\eta `$, there is a large scatter. The correspondence between $`\eta `$ and morphological type is investigated in more detail in Madgwick (2002).

We can also compare the value of $`\eta `$ derived for each galaxy to the equivalent width of the H$`\alpha `$ emission line, for a sample of high signal-to-noise ratio emission line spectra. This is shown in Fig. 3, from which it can be seen that there is a strong correlation between these quantities. It is well-known that the equivalent width of H$`\alpha `$ is a measure of star formation; indeed it has been used to estimate the birthrate parameter $`b`$ (the ratio of the present rate of star formation to the past averaged value — see e.g. Figure 3 of Kennicutt, Tamblyn & Congdon 1994). We return to this point in Section 3.2, where we use galaxy formation models to investigate how $`\eta `$ relates to the star-formation history.

## 3 Comparison with Simple Models

The results shown above suggest that $`\eta `$ is related both to galaxy morphology and to star formation history. In this section we further explore and test this hypothesis with models produced by combining simple parameterised star formation histories with spectral synthesis models. In this way, we can make a direct connection between a galaxy’s star formation history and the value of our spectral classification statistic, $`\eta `$.

### 3.1 Model Ingredients

Stellar population models combine theoretical temperature-luminosity tracks for stars of various masses with model stellar atmospheres and an assumed Initial Mass Function (IMF) in order to produce synthetic spectra for a “single burst” population of a uniform age and metalicity. In order to use these models to obtain predictions for a composite stellar population (i.e. a galaxy), the “single burst” models must be combined with an assumed star formation history, which provides the distribution of stellar ages and metalicities within the galaxy. Since the classical work on spectral synthesis by Tinsley (1980), many authors have adopted a simple parametrisation of the star formation history in terms of an exponential function of time (e.g. Bruzual 1983):

$$SFR(t)=\frac{1}{\tau }\mathrm{e}^{(tt_f)/\tau }.$$
(2)
Here, $`t_\mathrm{f}`$ is the time at which star formation first commences, and $`\tau `$ is the characteristic timescale of star formation. Most observed spectra can be modelled rather well using this simple function with the appropriate choice of $`\tau `$. A typical spiral galaxy spectrum might be well-fit with $`\tau 1`$ Gyr, while later type galaxies require larger values and early type galaxies require smaller values of $`\tau `$.

Nebular emission lines from ionised HII regions and extinction due to dust also contribute to the appearance of a galaxy spectrum. Modelling these processes is complicated, requiring assumptions about the ionization state, the geometry, the metalicity, the dust mass and the composition of HII regions.

The models presented in this Section are based on this simple but well-defined picture, and are realized using version 2.0 of the PEGASE software package (Fioc & Rocca-Volmerange) In all of the models shown we have adopted a Kennicutt (1983) initial-mass function with stellar masses in the range $`0.1120M_{\mathrm{}}`$. In addition, the metalicity of each galaxy is traced in a self-consistent manner using the prescription of Woosley & Weaver (1995), to model the enrichment of the inter-stellar medium.

Nebular emission lines result from light re-emitted by the ionised gas in star-forming regions. A prescription for determining the strengths of these lines is implemented in the PEGASE code. This process involves the absorption of Lyman continuum photons (below 912Å) by the nebular gas, which gets ionised, and reaches recombination equilibrium (Osterbrock 1989). It is assumed that 70% of these photons are absorbed by the gas at solar metalicity. In this approach, the strength of the nebular emission lines is a function of the age of the stellar population only, and metalicity and geometric effects are neglected.

The PEGASE code provides a simple way to model the effects of dust extinction on the synthetic spectra, in which the optical depth is estimated from the mass of gas and the metalicity. The absorption is then estimated using observational data for a mixture of graphites and silicates as in the Milky Way and the Magellanic Clouds (Draine & Lee 1984 and Pei 1992). In making this calculation it is also necessary to make an assumption about the geometry of the galaxy. We investigated models in which we assumed a spheroidal geometry and also an inclination-averaged disk geometry.

The spectra generated by the PEGASE code over the optical wavelength range are generally given in 10Å bins, and the emission line strengths are specified separately in terms of their peak fluxes over the continuum. In order to make comparisons between these spectra and the 2dFGRS we therefore interpolate the synthetic spectra onto the 4Å binned wavelength range of the 2dFGRS. Because the given emission lines are not resolved in this binning we create lines by superposing Gaussian profiles with the specified peak flux and FWHM corresponding to that calculated for the 2dF spectrograph using arc line measurements ($`2.3`$ pixels). We convert the given synthetic spectra to units of counts/bin by multiplying the flux (erg/s/Å) by the wavelength. We then normalise each spectrum to have mean counts of $`1`$ over our entire wavelength range. The spectra processed in this way are directly comparable to those of the 2dFGRS.

Spectral histories of a given galaxy are compiled for a range of formation times ($`t_f=020`$ Gyr), and a grid of values of $`\tau `$: 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 1.0, and 2.0 Gyr. We also calculate spectral histories for an instantaneous burst of star-formation at $`t=t_f`$. Because this grid does not comprise a well-defined sample of galaxies, to compare with the results from the 2dFGRS, we first subtract the mean spectrum from the 2dFGRS volume limited sample (described in the previous section) from each spectrum. We then compute the projection of each spectrum on the eigenspectra derived from the 2dFGRS volume limited sample. The first two projections ($`pc_1`$, $`pc_2`$) derived in this manner are shown for our grid of models in Fig. 4, along with the corresponding projections for a random subsample of galaxies in the $`M_{b_\mathrm{J}}5\mathrm{log}_{10}(h)<18`$ volume-limited 2dFGRS sample. Model grids are shown both with and without including the effects of dust extinction. It can be seen that dust does not have a dramatic impact on the results, but does seem to improve the agreement between the model tracks and the locus of the observed galaxies.

Already it is interesting that these simple models seem to cover the same range of the ($`pc1`$, $`pc2`$) parameter space occupied by the observed galaxy population. Evidently, at least in this simple scheme, galaxies are “born” in the lower right-hand corner of the plot, and progress towards the upper left-hand corner as they age. As well, galaxies in the upper right part of the diagram tend to have smaller values of the characteristic timescale for star formation, $`\tau `$, associated with early type galaxies, while as one moves diagonally towards the lower left corner, galaxies have the more extended timescales for star formation associated with late type galaxies. Perhaps unsurprisingly, the vast majority of the 2dF galaxies (which are fairly luminous) are consistent with times since formation $`\mathrm{\Delta }t`$ between 2–10 Gyr. It should be kept in mind, however, that the instrumental effects discussed in Section 2.2 can introduce random scatter in this diagram. We therefore turn now to the more robust $`\eta `$ parameter.

### 3.2 The Physical Significance of $`\eta `$

We argued based on a visual inspection of the eigenspectrum used to define the $`\eta `$ projection that $`\eta `$ essentially measures the strength of stellar and interstellar absorption line features and the strength of nebular emission line features. The strength of the absorption features mainly indicates the age and metalicity of the stellar population that dominates the optical luminosity of the galaxy, and the emission lines are strongly correlated with star formation activity. The comparison with classical galaxy classification methods such as visual morphology or equivalent width, shown in Section 2.3, suggested that $`\eta `$ may be connected with a measure of the star formation *relative* to the existing older stellar population, such as the birthrate parameter, $`b`$ (Scalo 1986). The birthrate parameter is defined as the ratio between the current and past-averaged star-formation rate, which in the case of the simple exponential star formation law is given analytically by

$$\frac{SFR(t)}{SFR(t)_{\mathrm{past}}}=\frac{(tt_f)}{\tau }\frac{\mathrm{e}^{(tt_f)/\tau }}{(1\mathrm{e}^{(tt_f)/\tau })}.$$
(3)
In this section we investigate our hypothesis that $`\eta `$ might be correlated with the birthrate parameter, $`b`$.

The relationship between $`\eta `$ and $`b`$ for the grid of models with simple star formation histories (as before) is shown in Fig. 5. For a fixed formation time, $`t_f`$, there is a well-defined relationship between the star formation history, as characterised by the birthrate parameter $`b`$, and the $`\eta `$ spectral parameter. This helps in the interpretation of $`\eta `$ for the observed population of galaxies. At least in this simple scheme, galaxies with small values of $`\eta `$ formed most of their stars in the past, and are currently evolving rather passively. This fits in well with the classical picture of early type galaxy formation. Galaxies with larger values of $`\eta `$ have had a significant amount of recent star formation, characteristic of spiral and irregular galaxies. Recalling from Fig. 1 that the 2dFGRS sample spans a range of $`\eta `$ values of about -5 to 10, this suggests again that most of the observed galaxies in the 2dF require times since formation greater than about 1 Gyr. From this diagram, we arrive at a prediction that the distribution of $`\eta `$ for galaxies observed at larger look-back times should exhibit a shift towards higher values. We find that the effect of extinction corrections is to slightly lower the calculated $`\eta `$ for a given star formation scenario, however, this is only a small effect and does not alter the qualitative interpretation that we have suggested.

A limitation of the simple characterisation of the star formation histories that we have adopted here is manifested by the asymptote of the $`b`$ values towards unity: because the star formation rate was assumed to be monotonically decreasing, $`b`$ is always less than one. It is clear that we must be cautious in over-interpreting these results, as real galaxies presumably may have much more complex star formation histories than we have assumed here. We now pursue a similar investigation for models with much more complex star formation histories, generated using a cosmological semi-analytic model of galaxy formation.

## 4 Comparison with Semi-Analytic Models

In this section, we again extract the $`\eta `$ parameter from model galaxy spectra created using stellar population synthesis models. The difference is that instead of using a simple parameterisation for the star formation history, the formation history of each galaxy is now modelled based on physical recipes set within the framework of hierarchical structure formation. The ensemble of model spectra that we use here was created using the Somerville et al. semi-analytic models (e.g. Somerville & Primack 1999; Somerville, Primack & Faber 2001) and has been described previously in Slonim et al. (2001). We refer the reader to those works for details, and here sketch the ingredients of the models very briefly.

In the hierarchical picture, present-day galaxies such as the ones observed by the 2dFGRS formed by the merging and accretion of smaller objects over time. In this picture, the star formation history of a galaxy is determined by the mass accretion and merging histories of its host dark matter halo, and the efficiency of gas cooling within those halos. As the hot gas is enriched with heavy elements by metal-rich winds from massive stars and supernovae, the cooling efficiency is increased. Galaxy mergers may trigger powerful bursts of star formation. In turn, violent star formation events may inhibit future star formation by heating the interstellar medium or driving winds that blow it out of the galaxy.

These star formation histories are convolved with stellar population models to produce model spectra in much the same way as in the “simple” models discussed in the previous section. In the models used here, we have used the multi-metalicity GISSEL models (Bruzual & Charlot, 1983) with a Salpeter IMF to calculate the stellar part of the spectra. These models are very similar to the PEGASE stellar population models used in the previous section, and for the purposes of this investigation, this difference should not significantly affect our results. Nebular emission lines are added to the spectra in the same way as before, using the empirical library from PEGASE.

Dust extinction is included using an approach similar to that of Guiderdoni & Rocca-Volmerange (1987), which is also very similar to the approach implemented in the PEGASE package. Here, the mass of dust is assumed to be proportional to the gas fraction times the metalicity of the cold gas. We then use a standard Galactic extinction curve and a “slab” model to compute the extinction as a function of wavelength and inclination.

The recipes for star formation, supernova feedback, and chemical evolution involve free parameters, which we set by requiring a halo with a virial velocity of 220 km/s to host a galaxy with an average I-band luminosity of about $`21.5+5\mathrm{log}h_{100}`$, and with an average gas fraction of $`0.1`$ to $`0.2`$, consistent with observations of local spiral galaxies. If we assume that mergers with mass ratios greater than $``$ 1:3 form spheroids, we find that the models produce the correct morphological mix of spirals, S0s and ellipticals at the present day. It has previously been shown in numerous works (e.g. Somerville & Primack 1999) that this approach leads to fairly good agreement with numerous key galaxy observations, such as the local luminosity function, colours, and clustering properties.

We construct a “mock 2dF catalogue” of $`2611`$ model galaxies with the same magnitude limit, wavelength coverage and spectral resolution, and redshift range as the 2dF survey. The synthetic spectra are expressed in terms of photon counts and the total number of counts in each spectrum is normalised to unity, as in the prepared observed spectra.

The calculation of the strength of the nebular emission features is one of the less precise aspects of the spectral synthesis packages. For this reason in Fig. 6 we show the first two principal components ($`pc_1`$,$`pc_2`$) for the semi-analytic model (SAM) ensemble of synthetic spectra. It can be seen that the distribution of the mock catalogue principal components is quite similar to that of the observed 2dFGRS galaxies. To guide the eye, the same set of evolutionary tracks shown in Fig. 4, derived from the “simple” models, is also shown. The most noticeable difference between the observed distribution and that of the mock, is that the former is substantially broader (most likely due to observational effects such as noise, reddening and evolution combined with the fact that the observed galaxies form a much larger sample). It is also interesting that the simple model tracks span the same locus of the ($`pc1`$, $`pc2`$) parameter space as the SAM galaxies, despite the fact that as we have emphasised the star formation histories are very different in the two kinds of models.

We can now discover whether the strong connection between $`\eta `$ and the birthrate parameter $`b`$ that we demonstrated for the simple models hold true for galaxies with more arbitrary star formation histories. In Fig. 7, we show $`bSFR/SFR_{\mathrm{past}}`$ vs. $`\eta `$ for the mock-2dF ensemble of SAM galaxy spectra. Again we see a strong relationship between $`\eta `$ and the $`b`$ parameter. Also, it is encouraging to note that the range of $`\eta `$ spanned by the SAM galaxies (-5 to 15) is similar to that spanned by the 2dF galaxies. The fact that the spectral properties of the galaxies in the SAM ensemble appear similar to those of the real 2dF galaxies enables us to use the results of the SAMs to draw a correspondence between the actual numerical value of $`b`$ and $`\eta `$. This leads to the interesting conclusion that the dividing line between the two “bumps” in the bimodal distribution of $`\eta `$ seen in Fig. 1, $`\eta 1.4`$, corresponds to galaxies that are forming stars at about 1/10th of their past-averaged rate. Similarly, the other two $`\eta `$ thresholds adopted by Madgwick et al. (2002) to calculate luminosity functions per spectral type may be matched up with $`b`$ values using Fig. 7.

We have also investigated correlations between $`\eta `$ and other physical properties of galaxies in the SAMs, such as bulge-to-total ratio or mean stellar age. While other correlations exist, none are as tight as the correlation shown in Fig. 7 between $`\eta `$ and the birthrate parameter $`b`$. Of all the correlations we investigated, the connection between $`\eta `$ and star formation history seems to present the most straightforward interpretation of the $`\eta `$ parameter. We defer a more detailed investigation of the other physical correlations in the SAMs and comparison with the 2dF and other data sets to a future work.

## 5 Discussion

In this paper we have studied the physical interpretation of a PCA-based spectral parameter, $`\eta `$, defined for the analysis of the 2dFGRS spectra (Madgwick et al. 2002). Although the definition of $`\eta `$ was motivated by the need to compromise between the desire to extract the maximum amount of statistical information from the spectra (in the sense of PCA’s maximum variance) and the limitations imposed by the instrumental uncertainties of the 2dF, we have argued here that this parameter has a straightforward and physically meaningful interpretation. We find that $`\eta `$ is an indicator of the star formation history of the galaxy, and is tightly correlated with the birthrate parameter $`b=SFR/SFR_{\mathrm{past}}`$, which characterises the ratio of present to past-averaged star formation.

A first indication of this correspondence is the correlation we find between $`\eta `$ and and the equivalent width of the H$`\alpha `$ emission line, EW(H$`\alpha `$), which has been used as a direct measure of the birthrate parameter in previous studies (e.g. Kennicutt, Tamblyn & Congdon 1994). We show that in models with simple star formation histories, parameterised in terms of an exponential function of time, there is a one-to-one correspondence between $`\eta `$ and the birthrate parameter $`b`$ (for a fixed formation time). Perhaps this is not too surprising in the context of classical spectral synthesis work — it has long been known (e.g. Tinsley 1980) that a simple SF history can be used to interpret spectral appearance, i.e. $`\eta `$ in this context. More surprisingly, we show also that a strong relationship between $`\eta `$ and $`b`$ is also exhibited by model galaxies with much more complex star formation histories, created within a hierarchical cosmological framework using a semi-analytic model.

The explanation of this surprisingly simple result must lie in the physical processes that determine the appearance of a galaxy spectrum. First, it is important to remember that at optical wavelengths, where all of this analysis has been carried out, the spectrum is dominated by the most luminous stars, and is therefore biased towards the most recent significant star formation activity. The slope of the optical continuum and the strength of stellar absorption lines evolve on the timescale of the lifetimes of intermediate type main sequence stars, roughly several Gigayears. The strengths of nebular emission features, however, evolve on much shorter timescales, as they require the presence of very hot, short lived O and B type stars. These features therefore depend on the star formation history on the timescale of 5–10 Megayears. The $`\eta `$ parameter has been defined so as to be insensitive to the continuum slope (for practical reasons), and hence represents a sequence from strong absorption lines to strong emission lines. It is well-known that the strengths of stellar absorption lines (especially Hydrogen recombination lines such as the Balmer series) are good indicators of stellar age, while the strength of the nebular emission lines are indicators of present star formation rate. With all this in mind, we can deduce that the rapid evolution in the value of $`\eta `$ over the first several hundred Myr of a galaxy’s life (see Fig. 5) reflects the fading of the emission lines as the star formation rate declines and the young stars burn out. Galaxies then tend to “pile up” at low values of $`\eta `$ as the age of their dominant stellar population exceeds a few Gyr.

In this sense it is straightforward to understand the relationship between $`\eta `$ and $`b`$ for our “simple” models, of monotonically decreasing star formation activity. However, the star formation histories of galaxies created by the semi-analytic models can be quite complex. The star formation rate fluctuates dramatically and non-monotonically over time as the galaxy exhausts its gas and then accretes a new gas supply, or experiences bursts of star formation triggered by mergers. The issue we must now address is how an optical diagnostic such as $`\eta `$, which is dominated by only the most recent events of star-formation, can be so closely related to the birth-rate parameter, $`b`$, which incorporates the entire formation history of a galaxy.

The basis of the hierarchical picture of galaxy formation — from which our mock catalogue of SAM galaxies has been derived — is that galaxy properties are determined by the merger histories of their host dark matter halos. In essence, the history of a specific halo (or galaxy) can be understood in terms of how the density on the scale of that structure compares to the background density. A region with density much above the average (a “many sigma” peak, in the language of Gaussian random fields) collapses early, while lower density peaks collapse later. It can therefore be conjectured that the early collapse of a dense dark matter halo can be associated with early star formation, early consumption of all available gas, and low present-day star formation rates. Conversely, late-forming objects will have a sustained gas supply and ongoing star formation. This suggests that the birthrate parameter $`b`$ reflects perhaps the most theoretically fundamental property of a galaxy (or of its host halo): small-$`b`$ galaxies represent rare, many-sigma peaks in the primordial density field, while large-$`b`$ galaxies are formed in more common, lower density peaks. This trend is manifested in the numerous empirical correlations between galaxy “type” (as characterised by morphology, colour, or spectral type), luminosity, and environment. We expect this simple picture to be complicated by the details of gas cooling, star formation, feedback, etc., but the results from the semi-analytic models (which include all of these effects at some level, although of course the real Universe is likely to be even more complicated) suggest that this only introduces a moderate amount of scatter on top of the general trend.

We therefore conclude that the $`\eta `$ parameter represents a promising candidate for galaxy classification in large modern redshift surveys. It was adopted in previous analyses of the 2dFGRS because of its practical advantages: it is straightforward and efficient to compute in an automated fashion, and it is robust to the instrumental uncertainties commonly associated with fibre-based multi-object spectrographs. Here we have shown that $`\eta `$ also has a straightforward physical interpretation which can be intuitively connected both with traditional classifiers such as morphology, and that it seems to be connected with theoretically fundamental properties of galaxies within the modern hierarchical structure formation paradigm.

## Acknowledgments

The efforts of the 2dFGRS collaboration, in preparing and compiling the data used in this analysis, are greatly appreciated. In particular we thank Karl Glazebrook for providing many useful comments and suggestions. DSM was supported by an Isaac Newton studentship from the Institute of Astronomy and Trinity College, Cambridge. RSS thanks the IoA for hospitality during completion of this work, and acknowledges support from the IoA visitor program.


format: text

text:
# Infrared Observations During the Secondary Eclipse of HD 209458b I. 3.6-Micron Occultation Spectroscopy Using the VLT 1footnote 11footnote 1Based on observations collected at the European Southern Observatory, Chile (proposal no. 67.C-0196)

## 1 INTRODUCTION

The discovery of the transiting extrasolar planet HD 209458b (Charbonneau et al. 2000; Henry et al. 2000) has provided a unique opportunity to deduce the physical characteristics of one example of the so-called ‘hot Jupiter’ class of exoplanets. The transit geometry has allowed an accurate derivation of the stellar and planetary mass and radius (Brown et al. 2001b) and also provides an opportunity to measure the composition of the exoplanet atmosphere. The scale height of an atmosphere in hydrostatic equilibrium is proportional to $`T/g`$, where $`T`$ is the atmospheric temperature and $`g`$ is the surface gravity. Because stellar irradiation of the planet is intense at $`0.046`$ AU (Seager & Sasselov 1998), the atmosphere could be as hot as $`T1400`$ K or more. Coupled with the relatively low surface gravity ($`g=848`$ cm s-2), the atmospheric scale height is potentially as large as $`H750`$ km. Moreover, when viewed tangent to the limb, the atmosphere can be opaque over several scale heights at the wavelength of a strong absorption line. The effective absorbing area of the planet can thus be greater in a strong line by a detectable amount, producing a slightly deeper eclipse at wavelengths close to the line core (Seager & Sasselov 2000; Brown 2001; Hubbard et al. 2001). That is, the opacity of the exoplanet atmosphere will impose weak, but potentially detectable, absorption lines on the stellar spectrum as the planet passes in front of the star. Several investigators have attempted to detect this effect (Bundy & Marcy 2000; Moutou et al. 2001). Using the Hubble Space Telescope to observe the primary eclipse, Charbonneau et al. (2002) succeeded in detecting the sodium doublet in the atmosphere of HD 209458b.

In addition to sodium, absorption in molecular features of water, methane, and carbon monoxide may also be observable (Hubbard et al. 2001; Brown 2001). Brown et al. (2002) attempted to detect absorption by carbon monoxide using ground-based observations. From their attempt and our own experience, we recognize that the terrestrial atmosphere can potentially frustrate such efforts, but we believe that good observing conditions and appropriate observational and analytical techniques will enable successful measurements from the ground. The known timing of the transit (Brown et al. 2001b) is a boon to this technique, because it discriminates against false signals by providing a temporal modulation on the signal.

In this paper we introduce a new observational approach involving the secondary eclipse, when the planet passes behind the star. The importance of the secondary eclipse was first pointed out by Charbonneau et al. (2000), but no measurement attempt has been reported in the literature. We have performed moderate resolution ‘occultation spectroscopy’ to measure the modulation (i.e., disappearance and reappearance) of the potential exoplanet line spectrum as the planet is hidden and revealed by the star. We have also performed lower resolution spectroscopy during and near secondary eclipse using the NASA Infrared Telescope Facility (IRTF), attempting to isolate the broad infrared (IR) flux peaks from the planet’s presumably hot atmosphere at wavelengths where the opacity is minimal. This paper introduces the concept of occultation spectroscopy for extrasolar planets and reports the Very Large Telescope (VLT) spectroscopic results; lower resolution occultation spectroscopy of the IR flux peaks using the IRTF observations will be the subject of a second paper.

## 2 OCCULTATION SPECTROSCOPY

In this section, we briefly explain our rationale for attempting occultation spectroscopy at secondary eclipse, as opposed to more conventional transmission spectroscopy during primary eclipse. At visible and near-IR wavelengths, the self-emitted flux of the planet is negligible, making occultation spectroscopy of thermal emission impossible. However, many of the strongest molecular features of interest in exoplanets (e.g., the very strong fundamental $`\nu _3`$ band of methane) occur at longer IR wavelengths. With increasing wavelength, the total signal detectable using occultation spectroscopy becomes as large, or larger, than that from transit spectroscopy. To see this, consider a spectral line having sufficient strength to be opaque $`N`$ scale heights higher in the atmosphere compared to the nearby continuum when viewed tangent to the limb. Cases of interest are for large $`N`$, so we can safely assume that the line is also optically thick when viewed at normal incidence to the exoplanet atmosphere. Transit spectroscopy looks for this absorption against a stellar continuum flux proportional $`2\pi RNHB_\lambda ^{}`$, where $`H`$ is the atmospheric scale height of the exoplanet, $`2\pi RNH`$ is the area of the annulus of the planetary atmosphere, and $`B_\lambda ^{}`$ is the Planck function for the stellar continuum temperature at wavelength $`\lambda `$. Occultation spectroscopy, on the other hand, looks for the modulation of the exoplanet spectrum itself, with continuum flux proportional to $`\pi R^2B_\lambda ^p`$, where $`B_\lambda ^p`$ is the Planck function for the exoplanet continuum temperature. With increasing wavelength, the increase in $`B_\lambda ^p`$ and the large area of the exoplanet ($`\pi R^2`$) dominates the area of the atmospheric annulus ($`2\pi RNH`$) in spite of the greater stellar intensity ($`B_\lambda ^{}`$). Evaluating the ratio of these fluxes for $`N=3`$, $`H=750`$ km, $`R=1.0\times 10^5`$ km, and for stellar and exoplanet temperatures of $`6000`$ K and $`1200`$ K, respectively, we find that the ratio equals unity for $`\lambda =2.6\mathrm{\mu m}`$. A rigorous treatment would of course include the line profiles and many other factors, but the inevitable tendency is that, for wavelengths longward of the K photometric band, occultation spectroscopy becomes increasingly favorable over transmission spectroscopy in terms of the amplitude of the potential signal.

Occultation spectroscopy is essentially direct IR spectroscopy of the exoplanet, and in principle this could be done at most orbital phases, and also done for planets which do not transit. However, it requires detection of the small-amplitude exoplanet spectrum (the precise nature of which is unknown), superposed on a strong stellar background. The measurement becomes much easier if the signal is modulated predictably. The disappearance of the exoplanet spectrum during secondary eclipse is therefore an essential part of our technique. Note, however, that modulation by variations in Doppler velocity could also be used for all hot Jupiter systems, even ones without transits (Wiedemann et al. 2001). Note also that the concept of occultation spectroscopy is not new; it was successfully used to obtain the IR spectra of Pluto and its moon Charon (Sawyer et al. 1987).

## 3 OBSERVATIONS

We observed two secondary eclipses, on UT 2001 July 8 and UT 2001 July 15, with the ISAAC spectrometer on the VLT (Antu) at Cerro Paranal in Chile. The observations were performed in visitor mode, with real-time decisions regarding nodding frequency, integration time, etc. Both nights were clear, with $`1020\%`$ relative humidity, and seeing in the range from $`0.50.7`$ arcsec. The zenith column of precipitable water vapor, as measured from telluric water lines in our spectra, was $`6`$ mm on July 8 and $`1.6`$ mm on July 15.

ISAAC is capable of imaging and spectroscopy in the wavelength region between 1 and 5 µm (Moorwood 1997). We used spectroscopic mode with a resolving power of 3300, covering the 3.5 to 3.75 µm range ($`26672852`$ cm-1). Both eclipses occurred within an hour of transit across the local meridian, where the star reached a minimum airmass of $`1.38`$. The spectra maintained good quality to the largest airmass we observed ($`2.51`$).

At approximately half-hour intervals, we observed a comparison star, HD 210483. The comparison star is within a degree of HD 209458, with nearly the same B and V magnitudes and photometric colors (Table 1). We observed both stars using a conventional ‘nod’ technique, placing them alternately at ‘a’ and ‘b’ positions on the slit, in the order ‘abba’, with $`60`$-second integrations at each slit position. On July 8 we obtained $`36`$ ‘abba’ sets of HD 209458, and $`20`$ ‘abba’ sets of HD 210483, while on July 15 we recorded $`41`$ and $`24`$ sets, respectively. About once per hour, we obtained spectra of a continuum lamp for flat-fielding. We did not rely on the standard flats; instead, we recorded flat field exposures immediately after each block of HD 209458 observations without moving the telescope or instrument, in order to avoid small but non-negligible flexure effects.

## 4 MODELED SPECTRUM

We do not expect to achieve sufficient precision to detect the planetary spectrum in each individual spectral resolution element; instead, we look for the candidate signal by performing a least-squares fit of the observations to a model ‘template’. We computed the template as the emergent flux from a model atmosphere for the exoplanet, using a series of specific intensity spectra at different values of $`\mu `$ (the cosine of the zenith angle). These were computed from the formal integral of the radiative transfer equation, using an LTE source function; we neglected all scattering terms, since the long-wavelength region is of interest and we are assuming a cloud-free model. Flux was computed by a quadrature integration of the intensity spectra. The spectrum synthesis was done at very high spectral resolution, and convolved to the $`0.8`$ cm-1 resolution of our ISAAC data using a Gaussian instrument profile. We adopted the parameters for the star and planet given by Cody & Sasselov (2002).

Model atmospheres for the planet can vary widely, with one major difference being the treatment of clouds and aerosols. The maximum thermal contrast in IR spectra is generally obtained from clear atmospheres; one reason for this is the fact that scattering processing in clouds tend to wash out strong absorption features. Moreover, the presence of clouds high in the atmosphere effectively limits the detection of methane by allowing only the small fraction of the column above the cloud to be probed by the observer. We have therefore used the most optimistic case of the clear atmosphere as a reference point. Our fiducial model follows the calculation of Goukenleuque et al. (2000), with the surface gravity appropriate for HD 209458b. This model includes irradiation by the star, but uses only gaseous opacities; condensates and particulates are assumed to have settled out of the atmosphere forming a deep cloud layer. Figure 1 compares the pressure/temperature structure from our fiducial model with a recent model from Sara Seager (2002, private communication). The temperature vs. optical depth structures of these models are similar, but the fiducial model requires a much higher pressure (due to lower opacity) to attain the same optical depth as the Seager model. (We also inspected a model for HD 209458b kindly sent to us by David Sudarsky (2002, private communication); it is slightly hotter than the Seager model at all depths.)

We calculated the number densities of methane in the fiducial model, using the simple analytic formulae given by Burrows & Sharp (1999). The methane mixing ratio peaks at $`7\times 10^6`$ at $`25`$ mbar pressure, similar to the values for the 51 Peg model (Goukenleuque et al. 2000). In computing the template, we used the current best-available experimental methane line parameters from Brown, Dulick, & Devi (2001a), now incorporated into HITRAN (Rothman et al. 1998), and the calculation followed Wiedemann et al. (2001) in other respects (partition functions, etc.). Since line parameters are critical to calculation of the template, we calculated spectra of hot methane ($`773`$ K) at ISAAC resolution for laboratory conditions, using both our adopted line data and also theoretical line parameters from the spherical top data system theory (Wenger & Champion 1998). We compared these calculated spectra to laboratory measurements of hot methane at $`773`$ K, measured by Steyert & Reuter (2003, private communication) at the same resolution. Based on this comparison, we have confidence in the HITRAN methane line data used to compute our model spectrum. Since our spectral region also exhibits absorption due to water, we experimented using water lines in our template, also based on HITRAN. Because the HITRAN database is incomplete at high temperatures, we also experimented with using HITEMP, a HITRAN line database extension for higher temperatures, to obtain water line data. This results in abundant weak absorption lines that do not significantly change the template. We have therefore neglected water absorption in the analysis until more complete information is available for individual hot water lines in this spectral region.

We expressed the exoplanet flux spectrum as a ratio to the stellar flux. We computed the stellar flux from a gray atmosphere at $`T=5800`$ K, reducing the adopted $`6000`$ K stellar temperature by $`200`$ K to allow for the decrease in brightness temperature from the visible to $`3.6\mathrm{\mu m}`$, as documented for the solar spectrum by Vernazza et al. (1976).

Figure 2 shows the modeled template spectrum, with and without water lines from HITRAN included. In the results quoted and discussed below, we have used the ‘methane only’ version of the model template, for reasons given above. However, our results and associated errors do not change significantly if ‘HITRAN water’ or ‘HITEMP water’ is included in the model template.

We computed methane spectra from models intermediate between our fiducial model and the Seager model. If we increase the continuous opacity in the fiducial model, while preserving the same temperature-optical depth relation, the methane bands weaken rapidly (again, see Figure 2). The greater continuous opacity pushes all line formation to lower pressure layers, giving lower column densities ‘above the continuum.’ The mixing ratio of methane decreases at low pressures (Burrows & Sharp 1999), further lowering the methane column density. Our fiducial model is a limiting case of an exceptionally clear atmosphere, and many other models for the exoplanet will not exhibit significant methane features, thereby escaping the test of our present analysis.

## 5 DATA ANALYSIS

Analysis of the HD 209458 observations took place in four stages: 1) extraction of spectra from the 2-D data frames, 2) removal of the telluric (and the few intrinsically stellar) absorption features from the spectra, yielding residuals which potentially contain the exoplanet spectrum, 3) estimation of the ‘amplitude’ of the exoplanet’s spectrum in each residual spectrum, and 4) fitting of the aggregate results to a secondary eclipse curve. Since our present (and future) results depend critically on the fidelity of our data analysis, we describe this process in some detail, and we also describe checks made to insure the integrity of the analysis.

### 5.1 Spectral Extraction

The spectral frames were first cleaned of ‘hot pixels’ and temporary anomalies such as energetic particle events. Temporal sequences of frames at a single nod position were examined to locate hot pixels based on their variability, and these pixels were given zero weight in subsequent analysis. Energetic particle events were identified and repaired separately using a median filter applied to the temporal sequence. The cleaned ‘abba’ set was combined by adding the ‘a’ frames and subtracting the ‘b’ frames to produce a difference image, which was then divided by a flat-field frame. The flat-field frame was computed as an average of the individual flats during the night. We verified that temporal variations in the flats were negligible for purposes of our analysis.

At an excellent IR site such as Paranal, changes in the thermal background are approximately linear with time during the several minutes required to record an ‘abba’ set. Consequently, most of the background emission is removed by construction of the difference image a-b-b+a. However, we found that a second order correction was necessary to remove the background completely. We compute the residual background spectrum by summing the region between the a and b slit positions in the difference image. The relatively long slit used in ISAAC facilitates the precision of this correction by providing a $`17`$-arcsec length of background. This residual background spectrum is subtracted from each row of the two-dimensional difference image prior to the extraction of the spectra. Note that our technique differs from usual IR practice, which would be to subtract the (negative) ‘b’ spectrum from the ‘a’ spectrum after their extraction from the difference image, which would automatically subtract the background to second order. However, this would not maintain separate ‘a’ and ‘b’ spectra, which we desired in order to verify that our results are independent of the position of the spectra on the detector array.

We extracted the ‘a’ and ‘b’ spectra from the difference image using the optimal extraction formalism of Horne (1986). The spatial profile needed in the optimal extraction was constructed by fitting polynomials in the intensity vs. wavelength distribution at each slit position, and then sampling these curves in the spatial direction at each wavelength. In the process of extracting the spectra, we noted that the approximate signal-to-noise ratio of each spectrum ($`200`$) was close to the value expected based on the noise from Poisson statistics of photoelectrons from the background.

At this stage of the analysis, the spectra were separated into sets corresponding to July 8 or July 15, HD 209458 or HD 210483, and ‘a’ or ‘b’ slit positions. Each of the 8 sets was analyzed independently of the other sets. The individual spectra in each set were shifted slightly using spline interpolation, so that the telluric absorption features were all coincident. Frequencies in wavenumber for the telluric lines were derived from high-resolution solar Fourier transform spectra (Livingston & Wallace 1991), convolved down to the resolution of the ISAAC data. (This accounted for any blending of telluric lines, which can change their effective wavenumbers.) Wavenumbers for each point in our ISAAC spectra were derived by spline interpolation using the telluric lines as standards; the error of this procedure was $`0.1`$ cm-1, much less than the $`0.8`$ cm-1 spectral resolution. A single accurate wavenumber scale was thereby inferred for all the spectra in each set. Sample spectra from both nights are shown in Figure 3.

### 5.2 Telluric Correction

The intensities in each spectrum fluctuate due to variable losses at the ISAAC slit, as well as changes in the telluric line and continuum absorption. For each set of spectra, we made an intensity normalization, followed by a correction to remove the telluric absorption. A pseudo-continuum was used as an aid to normalization; it was determined as a spline fitted to a set of points in each spectrum having the greatest intensity in pre-determined wavenumber intervals. We then used a wavelength-independent scale factor to normalize each spectrum so that the integral intensity under the pseudo-continuum was constant for all spectra. That is, we enforce the condition that the total intensity under each pseudo-continuum is equal to the total intensity under the average pseudo-continuum; the normalization factor $`f_i`$ is calculated from

$$f_i=\frac{_\lambda \overline{c}}{_\lambda c_i}$$
(1)
where $`\overline{c}`$ represents the average pseudo-continuum, $`c_i`$ represents the fitted pseudo-continuum to spectrum $`i`$, and $`f_i`$ is the factor by which each spectrum $`i`$ is multiplied. Thus, the normalized spectra are calculated by

$$n_i=f_is_i$$
(2)
where $`n_i`$ is the normalized spectrum and $`s_i`$ is the original spectrum.

In the normalized spectra, at each wavelength element, we fit a least-squares line to the natural log of intensity with airmass, $`\mathrm{ln}n_i(x,\lambda )`$. The fitted slope $`b_\lambda `$ was used to correct each intensity value to the minimum airmass in that set of spectra (this corrects for airmass-dependent telluric absorption). That is, the airmass-corrected spectrum $`n_i^{}`$ is obtained from

$$n_i^{}=\mathrm{exp}\left(\mathrm{ln}n_i(x,\lambda )b_\lambda (xx_{min})\right)$$
(3)
Each set of corrected spectra was then averaged. Residual spectra were calculated by subtracting the average spectrum and dividing the difference by a continuum fitted to the average spectrum. This produced residual spectra given by

$$r_i=\frac{n_i^{}\overline{n^{}}}{c_{\overline{n^{}}}}$$
(4)
which presumably contain the candidate exoplanet signal.

It is important to note several properties of the residual spectra. Division by the pseudo-continuum of the average spectrum leaves the residuals expressed in units of the stellar continuum. However, since we divided by only one continuum curve per data set, temporal changes in intensity were not masked by fluctuating errors in continuum fits. Also, the subtraction of the average spectrum means that the sum of the residuals at each wavelength is identically zero. Thus, only variations in the spectra survive this process. In particular, the exoplanet spectrum will survive (since it varies due to the eclipse), but the intrinsically stellar lines (Farmer 1989) are subtracted.

### 5.3 Higher-order Corrections

Several additional corrections were needed before the residuals could be compared with the model spectrum. First, we noticed that a broad absorption feature was evident in some of the residual spectra, and it corresponded to a blend of water lines near $`2720\text{cm}^1`$. The depth of this feature relative to the stellar continuum was $`30\%`$ greater on the more humid night (July 8), and variations of the depth of the feature on either night were about $`\pm 2\%`$. To correct these variations, we obtained a low-noise residual water spectrum by differencing the natural logarithms of the average spectrum for July 8 and July 15 (for both stars, and for both the ‘a’ and ‘b’ spectra, separately). This exploits the (fortuitous) fact that the precipitable water vapor column was appreciably different on the two nights, but telluric methane absorption was the same. (Unlike the concentration of water in the troposphere, which is highly variable, methane is effectively constant because it is long-lived and therefore well-mixed.) We tested each residual spectrum to see whether a significant correlation existed with the water vapor spectrum. If it did, then we removed it from the residuals by least-squares fitting.

We also noticed that the residuals sometimes exhibit variations at low spatial frequencies (see upper panel of Figure 4). These baseline effects (typically $`0.20.3\%`$) are not surprising, because our analysis treats each wavelength point independently. (We do not introduce artificial baselines - e.g., by continuum fitting - but neither do we suppress real baseline variations.) A Fourier analysis also revealed a spike of excess power at the Nyquist frequency ($`0.5`$ cycles per pixel), which is commonly seen in IR data, usually indicating a difference in the two output amplifiers reading alternate columns of the detector array. We removed these effects by computing the Fourier transform for each set of residuals; then we zeroed-out an appropriate number of low frequency elements, as well as the DC offset (zero frequency) and the Nyquist frequency components. The number of low frequencies removed (up to $`0.0079`$ cycles per pixel) was taken to be the maximum possible without impinging on frequencies where the exoplanet model spectrum exhibited significant power. The same filter operations were performed on the model template spectrum as on the observed residuals. Then we inverse-transformed the filtered residuals and template, which were then suitable for comparison. The effect of this Fourier filtering on the observed residuals is seen the lower panel of Figure 4.

Even after the removal of variable water vapor absorption, and Fourier filtering, the residuals needed an additional correction. We discovered that they contained an excess of ‘outliers’, i.e., points more than several standard deviations from zero. We exploited our comparison star observations to find those wavelengths which tended to become outliers and reject them from the analysis. The procedure was applied to the ‘a’ and ‘b’ spectra separately, since we suspected low-amplitude ‘hot pixels’ as a likely cause of outliers. We applied a wavelength-dependent ‘weighting mask’ to the residuals, where the weight of the mask (normally unity) was set to zero when a particular wavelength was found to deviate from zero by more than $`2.5\sigma `$ in two or more spectra for both HD 209458 and HD 210483. The requirement that points are zero-weighted only if that wavelength tends to be discrepant in spectra from both stars assures that this procedure does not bias the results. About 7% of the points were zero-weighted using this method; we found that our final results were remarkably insensitive to the limits for the rejection. The reason for this is that the total number of outliers was much less than $`7\%`$, because a given point was always zero-weighted when it met the above criterion, even if it was ‘well-behaved’ in the majority of spectra. The standard deviation of the Fourier-filtered and masked residual spectra was in the range $`0.003`$ to $`0.006`$.

### 5.4 Fit to the Exoplanet Spectrum

Prior to fitting the model template and Fourier filtering, the residuals are corrected for the heliocentric Doppler shift due to the relative motion of HD 209458 with respect to the Sun (Nidever et al. 2002), as well as the geocentric Doppler shift due to the relative motion of the Earth around the Sun. We have also corrected for the non-negligible Doppler shift due to the orbital motion of the exoplanet with respect to HD 209458. The residuals are shifted in wavenumber to the rest frame of the exoplanet, and the model spectrum (already in the exoplanet rest frame) is interpolated onto the same wavenumber grid. Note that although the correction is small, we have corrected each residual (each frame) for the Doppler effect, based on its time of observation. The residuals calculated for HD 210483 are analyzed in an identical manner. Although HD 210483 is not known to host a planet, we apply the same Doppler correction as we did to HD 209458, but we calculate the correction based on the time of observation of each HD 210483 residual.

The comparison between each residual spectrum and the model exoplanet spectrum is made using linear least-squares. For this purpose we use the ‘methane only’ version of the model (see Section 4), since the state of the spectroscopic line data gives us the most confidence in this version. (Our overall results do not change with the ‘methane plus water’ version.) The result of a least-squares regression of the residuals versus the model spectrum gives a quantitative estimate of the degree to which each residual spectrum ‘contains’ the exoplanet model spectrum. We refer to these estimates as ‘model amplitudes’, since they represent the degree to which the model spectrum appears in each residual spectrum. The least-squares fit also estimates the random error in the model amplitudes, typically $`1.52.5`$. In other words, the signal-to-noise ratio for detecting the modeled exoplanet spectrum in a single set of residuals is about $`0.5`$. This is sufficient precision to expect detection of the exoplanet spectrum in the average of $`50`$ residual spectra, provided that the signal-to-noise ratio increases as the square root of the number of spectra averaged.

### 5.5 Checks on the Analysis

We have checked our analysis procedures in several ways. Since we use linear least-squares, the analysis should have the property that averaging the best-fit model amplitudes for a set of residual spectra should give the same result as fitting to the average of those residual spectra. We verified that our results satisfy this identity exactly, if the time-variable Doppler shift applied to the residuals is neglected (with Doppler shifts included, the equality is approximate.) We have also confirmed that our numerical procedures do not attenuate potential signals. We added a synthetic signal, identical to the modeled spectrum, at the earliest practical stage of the analysis (immediately after the extraction of the spectra from the 2-D data frames). The least-squares solutions recover this signal (as a difference, out of eclipse minus in eclipse) with a best-fit model amplitude of near unity.

Another important check on our analysis is to examine the nature of the noise, specifically the distribution of noise as compared to a normal error distribution. This is important because our error estimates implicitly assume that the data reflect a normal error distribution. Figure 5 shows the error distribution for all of the points (after masking) in the ‘a’ residuals of HD 209458 from July 8. Apart from an excess of points within $`1\sigma `$ of zero, the distribution closely approximates a normal error curve. Similar distributions are found for the other sets of data.

Figure 6 shows the distribution of ‘model amplitudes.’ In this case we construct the distribution from the least-squares fits to all residual spectra for both HD 209458 and HD 210483 from both nights; this comprehensive inclusion is needed to provide sufficient points to expect a reasonable approximation to the normal error curve. Considering the number of residual spectra represented (242), the data are in good agreement with the normal error curve. As a further check on the errors, we computed the variation in fitted model amplitudes which result from fitting Gaussian random noise whose standard deviation varied from $`0.003`$ to $`0.006`$, the same as the range of our filtered and masked residuals. The results of this simulation were in close agreement with the error estimates derived from the least-squares fits to the actual residual data.

### 5.6 Secondary Eclipse Timing

The timing of secondary eclipse depends on the orbital eccentricity. If the eccentricity is identically zero (as might be expected based on tidal circularization arguments) the secondary eclipse will occur exactly mid-way between primary eclipses. If the eccentricity is non-zero, it can occur up to $`90`$ minutes earlier or later (depending also on $`\omega `$). Doppler observations to date (G. W. Marcy 2002, private communication) give $`e=0.011\pm 0.015`$ with $`\omega =156^{}`$, which would imply $`\delta t=31`$ minutes, i.e., a slightly earlier secondary eclipse. We computed the time of secondary eclipse as a function of the eccentricity, using $`P=3.52474`$ days (G. W. Marcy 2002, private communication), and the zero-point from primary eclipse at $`T_0=2451659.93675`$ HJD (Brown et al. 2001b), and we corrected for light travel time. We checked our calculations against a primary eclipse ephemeris by David Charbonneau (2002, http://www.astro.caltech.edu/$$dc/frames.html), obtaining essentially identical results, and we did a similar check for the time difference introduced by using a non-zero eccentricity (D. Charbonneau 2002, private communication). Given a calculated time for secondary eclipse, we fit a simple eclipse curve to our model amplitude data by linear least-squares. Our eclipse curve uses two levels (in- and out-of-eclipse) connected by straight lines. The time from first to fourth contact was taken to be $`184.3`$ minutes, and from second to third contact $`132.2`$ minutes (Charbonneau et al. 2002). Note that, due to the subtraction of an average spectrum in our data analysis, only the amplitude of the eclipse curve is significant; the ‘zero-point’ in the fitted eclipse curve is not meaningful.

## 6 RESULTS AND DISCUSSION

Given the precision achieved by our analysis, we would expect to detect the exoplanet spectrum if it is represented by the irradiated, cloudless, low-albedo atmosphere, with thermochemical equilibrium abundances of methane (and water), and also provided that the timing of the secondary eclipse is known and the stellar and exoplanet parameters (radii, etc.) are exactly as adopted. We first conclude that this precision can indeed be achieved using ground-based observations. The recent results and discussion of Brown et al. (2002) imply that ground-based detection is feasible, since these authors came within a factor of $`3`$ of the required precision using data from only one night in poor weather. Our experience, using a somewhat different technique, confirms that the requisite precision can be obtained if good observing conditions prevail.

Our results for the secondary eclipse are illustrated in Figures 7 and 8, showing the derived ‘model amplitudes’ for each HD 209458 spectrum (‘a’ and ‘b’ for both nights) versus time from the assumed mid-point of secondary eclipse. Recall that the error bars for individual spectra were derived from the least-squares fits to the residuals after telluric corrections, Fourier filtering, and masking. We fit an eclipse curve to the aggregate results at each assumed eclipse time. (We fold the curve about the mid-point of secondary eclipse for the figures.) Also, there is no evidence that the results from the ‘a’ and ‘b’ spectra analyzed separately are significantly different. Figure 7 shows the result for $`e=0.011`$, i.e., taking the slightly non-zero eccentricity from the Doppler data at face value, again for all spectra. The solid line is the least-squares estimate of the eclipse amplitude, $`0.1\pm 0.3`$, and the dashed line shows the levels corresponding to an eclipse amplitude of unity (i.e., if the model template spectrum correctly represents the exoplanet, and the timing is correct.) Figure 8 shows the corresponding result ($`0.5\pm 0.4`$) under the assumption that the eccentricity is identically zero. In each case the reduced chi-squared of the fit is $`1.5`$, indicating that the scatter of the data is only slightly larger than the independently-assigned error bars. Note that each plot for HD 209458 is accompanied by a similar fit to the comparison star, shown in the lower panels. In both cases, the comparison star shows no significant change, as expected.

Adopting $`e=0.011`$ from the Doppler data, our results for non-detection of the secondary eclipse ($`0.1\pm 0.3`$) exclude the model template spectrum at the $`>3\sigma `$ level. Given our error distributions, this case is firmly rejected. Assuming an eccentricity of zero, non-detection of the secondary eclipse is ambiguous ($`0.5\pm 0.4`$). However, even if we explore eccentricities over the entire plausible range from $`0.0`$ to $`0.03`$, and assume $`\omega `$ values in the first and second quadrants, the fitted secondary eclipse amplitudes never reach unity. Instead, the amplitudes vary from $`0.2`$ to $`+0.9`$, and the amplitudes for the comparison star similarly range from $`0.5`$ to $`+0.9`$. At a given eclipse time, the amplitudes for the two stars are uncorrelated, as expected when sampling the noise envelopes resulting from fitting independent random data to eclipse curves of variable timings. We therefore conclude that the secondary eclipse is not detected using our ‘methane template.’ Given the similarity between HD 209458 and the comparison star in terms of the noise envelope, and noting the nearly identical reduced chi-squared values, we conclude that HD 209458b probably does not exhibit significant methane absorption features, in agreement with the models of Sudarsky et al. (2002). Methane absorption in the combined light spectrum (star plus planet) has a total equivalent width of no more than $`0.025`$ cm-1 within our absorption bandpass ($`26672852`$ cm-1).

Deviation of the exoplanet atmosphere from our fiducial model is not difficult to explain. As discussed in Section 4, models of this exoplanet must be exceptionally clear if significant methane absorption is to appear in their spectrum. In this respect our fiducial model represents a limiting case. Our results are certainly consistent with the suggestion of a cloudy atmosphere based on the observed low sodium abundance (Charbonneau et al. 2002). However, even in a clear atmosphere there may be other viable explanations for both the sodium result (Barman et al. 2002), and our present result (e.g., photochemical depletion of methane). Nevertheless, it is significant that we can now begin to limit the range of parameters of viable exoplanet models using ground-based observations.

## 7 Acknowledgements

Jeremy Richardson is supported in part by a NASA Graduate Student Researchers Program Fellowship, funded by the Office of Space Science at NASA Headquarters (grant number NGT5-50273). Other aspects of this research were supported by the NASA Origins of Solar Systems program. We extend special thanks to Gianni Marconi, the VLT night astronomer, who contributed substantially to the success of these technically challenging observations. We thank Dave Charbonneau for several discussions concerning the secondary eclipse, and for sharing his insights generally. We benefitted from discussions with Sara Seager on the characteristics of exoplanet model atmospheres, and we thank her for providing her recent temperature/pressure relation. We are grateful to Geoff Marcy for communicating the results of the latest solution for orbital eccentricity, and to Dave Sudarsky for sending us his temperature/pressure profile for HD 209458b. We also thank Linda Brown for discussions on the status of methane line data. Finally, we thank the referee for helpful comments and suggestions.