text:
# Precision Cosmology from the Lyman-𝛼 Forest: Power Spectrum and Bispectrum

## 1 Introduction

The study of the Ly-$`\alpha `$ forest has been revolutionized in recent years by high resolution measurements using the Keck HIRES spectrograph (Vogt et al., 1994) and by the development of theoretical understanding using hydrodynamical simulations (Cen et al., 1994; Zhang et al., 1995; Hernquist et al., 1996; Theuns et al., 1998) and analytical models (Gnedin & Hui, 1998). The picture that has emerged from these studies is one in which the neutral gas responsible for the absorption is in a relatively low density, smooth environment, which implies a simple connection between the gas and the underlying dark matter. The neutral fraction of the gas is determined by the interplay between the recombination rate (which depends on the temperature of the gas) and ionization caused by ultraviolet photons. Photoionization heating and expansion cooling cause the gas density and temperature to be tightly related, except where mild shocks heat up the gas. This leads to a tight relation between the absorption and the gas density. Finally, the gas density is closely related to the dark matter density on large scales, while on small scales the effects of thermal broadening and Jeans smoothing have to be included. In the simplest picture described here all of the physics ingredients are known and can be modeled.

Within the model described above, the relation between the observed flux and the underlying dark matter is completely defined, enabling one to study the dark matter correlations with the help of the Lyman-$`\alpha `$ forest. This was proposed first by Croft et al. (1998) and has been subsequently investigated and applied to the real data by several groups (Croft et al., 1999, 2002b; McDonald et al., 2000; Hui et al., 2001). This work established the simple flux power spectrum as the statistic of choice, although other statistics such as the flux probability distribution have also been investigated (McDonald et al., 2000).

This paper addresses in more detail the question of how much information about cosmological parameters can be extracted from the analysis of the Lyman-$`\alpha `$ forest and what are the best statistics to use. We focus on N-point statistics (where N=2, 3, 4), since these provide the simplest parametrization of the long-range correlations. We use a Fisher matrix analysis to assess the uncertainties in cosmological parameters for a given uncertainty in Lya forest measurements. To calculate the Fisher matrix, we must calculate the derivatives of the observables (flux power spectrum and bispectrum) with respect to cosmological parameters, centered on a reasonable fiducial model. Using this approach, we show that higher order statistics add to and independently confirm the information gained from the flux power spectrum, and we address practical issues of such a study. For the data sample, we assume a sample with noise and resolution characteristics of the Sloan Digital Sky Survey, which currently contains a few thousand QSO spectra with measured Lyman-$`\alpha `$ forest ($`z>2.2`$).

Statistical correlations in the Lyman-$`\alpha `$ forest are sensitive to several parameters of cosmological interest. Broadly, they are sensitive to the linear power spectrum of matter fluctuations on scales around 1 comoving $`h^1`$Mpc. In this paper, this power spectrum $`P_L(k,a)`$ is expressed as

$$P_L(k,a)=AD^2(a)k^{n+(\alpha /2)\mathrm{ln}(k/k_p)}T^2(k).$$
(1)
The parameters that are being studied are the power spectrum amplitude ($`\kappa \mathrm{ln}\sqrt{A}`$), the primordial slope $`n`$, and primordial curvature $`\alpha `$. We also investigate the sensitivity to the linear growth factor $`D(a)`$, where $`a=(1+z)^1`$ is the expansion factor.

Several additional parameters affecting absorption in the Lyman-$`\alpha `$ forest must be included, such as those of the gas temperature-density relation, which for our redshift range is assumed to have the form

$$T(\rho )=T_{1.4}\left(\frac{\rho }{1.4\overline{\rho }}\right)^{\gamma 1}.$$
(2)
$`T_{1.4}`$ is used because McDonald et al. (2001) determined the temperature at this density most precisely. $`T_{1.4}`$ and $`\gamma `$ are two of the parameters studied in this paper. In addition, the mean transmitted flux $`\overline{F}(z)`$ is also assumed to be a free parameter (which can be related to the UV background and baryon density though the equation of ionizing equilibrium). Thus six parameters were varied at each redshift. Our approach is to use hydrodynamical and N-body simulations, and vary all the parameters at each redshift. We are interested in the sensitivity to cosmological parameters, so we marginalize over all the other parameters when presenting the results.

The outline of the paper is the following. In section 2 we describe our HPM simulations and the hydrodynamic simulations that are used for comparison, introduce the statistics used in the course of the analysis, and discuss our choice of simulations, filters and instrumental effects for this study. Section 3 describes the procedure of calculating the Fisher matrix, and the determination of errorbars for the parameters of interest. In section 4, we address the possibility of using the Lyman-$`\alpha `$ forest to determine quintessence density $`\mathrm{\Omega }_q`$ or equation of state $`w_q`$, and more general deviations from the Einstein-de Sitter growth factor $`D(a)=a`$. After discussing the implications of the results in section 5, we conclude with a discussion of future work in this area.

## 2 Simulations

To assess the sensitivity of the Ly-$`\alpha `$ forest to cosmological parameters we performed a Fisher matrix analysis (e.g. Tegmark et al. (1997)). We conducted a redshift-dependent analysis using 9 redshift bins, centered at $`z=2.2,2.4,\mathrm{},3.8`$ and using a realistic redshift distribution of quasars as found in SDSS data. Since we only had simulations with values of $`a`$ spaced by 0.04, we chose the closest $`a`$ available.

We use Hydro-PM (HPM) simulations (Gnedin & Hui, 1998; McDonald, 2001; McDonald et al., 2002; Cen & McDonald, 2002) rather than more accurate hydrodynamic simulations due to the very large number of simulations needed. The boxes used for the main analysis had length 80 $`h^1`$Mpc, with $`N=512^3`$ particles, and spectra were averaged over all three axes. Pixelization, smoothing, and noise of the data were imitated in the simulations, with $`S/N=5`$ and Gaussian smoothing $`\sigma `$ of length 0.625 comoving $`h^1`$Mpc (same size as the pixels), as will be described in more detail in section 2.7. $`64^2`$ lines of sight into each face of the box were used to reduce computation time after finding that using $`128^2`$ did not yield more information.

The fiducial model used for the full Fisher analysis is $`\gamma 1=0.3`$, $`T_{1.4}=17000`$ K, $`n=0.75`$, $`\alpha =0.0`$, $`\mathrm{\Omega }_{m,0}h=0.26`$, $`\sigma _8=0.75`$ and $`\overline{F}(z)`$ from

$$\overline{F}(z)=e^{e^{0.95+3.77\mathrm{ln}[(1+z)/4.0]}},$$
(3)
a fit to the data in McDonald et al. (2000). Note that $`n`$ and $`\alpha `$ here are the primordial slope and curvature. The values measured are $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$, which include a contribution from the transfer function. The value chosen for $`\gamma 1`$ is the center of the range of theoretical values, and two-sided derivatives were done in steps of $`\pm 0.03`$ (in other words, using simulations with $`\gamma 1=0.0`$ and 0.6). Likewise, the value chosen for $`T_{1.4}`$ is observationally justifiable, and derivatives were taken in steps of $`\pm 5000`$ K. Derivatives in $`n`$ about 0.75 and $`\alpha `$ about 0.0 (a value consistent with observation as well) were both done in steps of $`\pm 0.05`$. Derivatives in $`\overline{F}(z)`$ were in steps of $`\pm 0.02`$. Instead of varying the actual power spectrum amplitude, which requires a large number of different simulations, we varied $`a`$. In relating derivatives such as $`\mathrm{d}P/\mathrm{d}a`$ to the desired derivatives $`\mathrm{d}P/\mathrm{d}\kappa `$, we used the approximation (shown in McDonald (2001) to be quite accurate) that $`z`$-evolution of the power spectrum can be treated as a rescaling of $`\kappa `$ and $`T_{1.4}`$. The change in growth factor can be treated primarily as a change in $`\kappa `$, along with a small change in $`T_{1.4}`$ (fixed in velocity coordinates, it consequently differs in simulation coordinates). Two-sided derivatives with respect to $`a`$ were calculated using steps of $`\pm 0.04`$. Then, by subtracting off this temperature term from the $`a`$-derivatives, we obtain the portion of the $`a`$-derivative that mimics a change in power spectrum amplitude evolution. Then we use $`\mathrm{d}\kappa =\mathrm{d}a/a`$ to form the amplitude derivative (this form assumes an Einstein-de Sitter universe with $`D(a)=a`$).

While the bulk of our simulations are HPM, we use the hydrodynamic simulations for comparison whenever possible. These are described in Cen et al. (2001); McDonald et al. (2002); Cen & McDonald (2002); and for more detail, Cen et al. (2002). They have box size 25 $`h^1`$Mpc, with $`768^3`$ Eulerian cells for baryons and $`384^3`$ dark matter particles, and the cosmology is $`\mathrm{\Lambda }`$CDM, with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_b=0.035`$, $`\sigma _8=0.9`$, $`h=0.67`$, and $`n=1`$. Median temperature at the mean density is around 15,000 K, slightly decreasing with redshift for $`2<z<4`$, consistent with observational constraints in McDonald et al. (2000). The three outputs used have redshifts 1.9, 2.45, and 3.0. The HPM simulations used for the comparison with hydrodynamic simulations had a smaller box size than those used for the full analysis. The parameters were box size 25 $`h^1`$Mpc, $`384^3`$ particles, the same cosmology as the hydrodynamic simulations, and $`\gamma 1`$=0.6 (the temperature-density relation in the hydrodynamic simulations is not precisely a power law, but the slope is near 0.6). The two had nearly identical initial conditions, minimizing the sampling variance. The missing mass that was concentrated into stars in the hydrodynamic simulations was neglected in the HPM simulations. In all cases, $`\overline{F}=0.7`$ was used for the comparison, with 2-sided derivatives using increments of $`\pm 0.02`$. The simulations were binned so that the hydrodynamic and HPM simulations had the same size pixels.

We now introduce the statistics used for analysis, and then use them to address several questions about these simulations: accuracy in comparison with hydrodynamic simulations, sufficiency of the box size and resolution, the magnitude of random fluctuations due to finite box size, and the effects of chunking on analysis. We discuss all but the last of these issues only to determine the sufficiency of our simulations for a Fisher matrix study, rather than requiring the higher precision necessary for their use in a data analysis. The final issue, chunking, is discussed since it may be done in a data analysis, so it is necessary to show that it does not affect the results.

### 2.1 Statistics

In general, the statistics were computed using the flux residual $`\delta _F=F/\overline{F}1`$ rather than $`F`$ itself. Besides the flux power spectrum $`P`$ with normalization

$$_{\mathrm{}}^+\mathrm{}\frac{dk}{2\pi }P(k)=\delta _F^2,$$
(4)
several higher order statistics were studied. Rather than investigate the full bispectrum and trispectrum information, we decided to use a more restricted form of these general statistics. In this sense our results will be conservative. The higher order statistics are created using the following procedure: $`\delta _F(\lambda )`$ is Fourier transformed to find $`\stackrel{~}{\delta }_F(k)`$. Next, it is band-pass filtered to create the field

$$\stackrel{~}{\delta }_H(k)=\stackrel{~}{\delta }_F(k)W(k)$$
(5)
where $`W(k)`$ is some filtering function. Two types of filters used in this analysis are a square window,

$$W_s(k)=\{\begin{array}{cc}0,\hfill & \text{if }k<k_1\text{ or }k>k_2\hfill \\ 1,\hfill & \text{if }k_1kk_2\hfill \end{array}$$
(6)
and a Gaussian window characterized by some $`\overline{k}`$ and $`\sigma _k`$.

Next, the filtered field $`\stackrel{~}{\delta }_H(k)`$ is inverse transformed to create the real-space counterpart, $`\delta _H(\lambda )`$. The field $`H(\lambda )`$ is created by squaring:

$$H(\lambda )=(\delta _H(\lambda ))^2.$$
(7)
We define $`T`$ as the power spectrum of $`H`$. Another statistic, $`B`$, is the cross-spectrum between $`\delta _F`$ and $`H`$. $`B`$ and $`T`$ are related to the three- and four-point functions, respectively. However, because $`T`$ is not the reduced four-point function, there is a significant covariance between it and $`P`$.

In addition, the cross-correlation coefficient

$$R=\frac{B}{\sqrt{PT}}$$
(8)
was studied by Zaldarriaga et al. (2001b) due to its simple physical interpretation. As shown there, it is expected that $`R`$ computed from the flux residual should be significantly negative for low $`k`$, since perturbations in high density (low flux) regions tend to grow faster. As shown in this paper, one utility of these higher order statistics is that they break degeneracies between the power spectrum amplitude and $`\overline{F}`$ as measured by $`P`$ alone. Furthermore, as shown in Zaldarriaga et al. (2001b), they may also be useful in discriminating between gravitational processes (which give negative $`B`$ and $`R`$) and continuum fluctuations (which give positive $`B`$ and $`R`$) or other extraneous effects such as metal lines, star-formation induced outflows of gas, inhomogeneity in the UV background, and the limitations of simulations.

Fig. 1 shows $`P`$, $`B`$, $`T`$, and $`R`$ for the central point in parameter space around which variation occurs for $`z=3.0`$ ($`\gamma 1=0.3`$, $`T_{1.4}=17000`$ K, $`a=0.24`$, $`n=0.75`$, $`\alpha =0.0`$, $`\overline{F}=0.68`$). Note that we vary $`a`$ rather than the power spectrum amplitude itself (which involves generating many more simulations). By Eq. 1, we have $`\delta \kappa =\delta a/a`$ for an Einstein-de Sitter cosmology ($`D(a)=a`$). Errorbars shown are computed using the variances of the statistics for a single simulation and the length of spectrum for this redshift in the SDSS data available to us.

Our study showed that $`T`$ and $`R`$ do not add much information, since $`T`$ is so closely related to $`P`$, and $`R`$ is just a combination of the other statistics. Consequently, in the rest of this paper, only $`P`$ and $`B`$ will be used. Figs. 2 and 3 show the relative variation of these statistics with all 6 parameters. Note that in this section, all plots of the values of statistics or their derivatives were created using the same filter, a square filter from 0.2-1.1 $`h`$Mpc-1 (unless noted otherwise), and derivative plots have the same scale for easy comparison.

We briefly explain the effects of changing the various parameters on $`P`$ and $`B`$ that are illustrated in Figs. 2 and 3. Raising the overall power spectrum amplitude causes an increase in $`P`$ and $`B`$ on all scales shown here, because the increase in power in the matter power spectrum carries directly over to the flux power spectrum. Increasing $`n`$ lowers $`P`$ and $`B`$ on large scales and raises it on small scales, directly reflecting its effects on the matter power spectrum (for even smaller scales than those shown here, increasing the power spectrum amplitude and $`n`$ actually causes a decrease in $`P`$ because non-linear peculiar velocities suppress power, as in McDonald (2001)). The statistics are nearly insensitive to $`\alpha `$ because it does not have much effect on these scales. $`P`$ and $`B`$ show only a very slight increase with $`T_{1.4}`$. The decrease of these statistics with an increase in $`\gamma 1`$, consistent with McDonald (2001), can be explained by considering the optical depth $`\tau \rho _b^{20.7(\gamma 1)}`$, so that the exponent on $`\rho _b`$ is decreasing, and hence the fluctuations in optical depth for a given density fluctuation also decrease. The decrease in power with an increase in the mean flux is consistent with the results in McDonald (2001), and can be understood by the proportionality between the optical depth and $`\rho _b^{20.7(\gamma 1)}`$. Because the increased mean flux is equivalent to a lower optical depth, we can see that this should have the same effect as lowering the density fluctuations, and consequently the flux power spectrum (hence the degeneracy between power spectrum amplitude and mean flux).

### 2.2 Comparison of HPM and Hydro simulations

Fig. 4 shows the relative difference between HPM and hydrodynamic simulations at three redshifts, Fig. 5 shows their flux derivatives at $`z=2.45`$, and Fig. 6 shows their amplitude derivatives at $`z=2.45`$.

For the derivatives, the relative change (e.g., $`\delta P/P`$) was used for comparison, because differences in the magnitude of the derivatives could be simply artifacts of the differences in the magnitudes of the statistics themselves. Since this kind of change will cancel out of a Fisher matrix calculation (see equation 11), the relative change is plotted to emphasize changes in the $`k`$-dependent shape of the derivatives rather than their amplitude.

As illustrated by Figs. 4-6, the HPM simulations yield statistics and derivatives fairly similar to those of the hydrodynamic simulations. The worst discrepancies in Fig. 4 occur for $`z=1.9`$, already below the range of this study. A Fisher matrix analysis using just amplitude and $`\overline{F}`$ was completed for these simulations, and found that the error bounds for both parameters were 6% lower when calculated with the HPM than with the hydro simulations, a tolerable difference for our purposes.

### 2.3 Resolution and box size

Another question about the HPM simulations is the box size and resolution necessary for a realistic study of these statistics. To check whether the 80 $`h^1`$Mpc box with $`512^3`$ particles is sufficient, we computed the statistics and their derivatives for a 40 $`h^1`$Mpc box with $`512^3`$ particles (higher resolution), and a 40 $`h^1`$Mpc box with $`256^3`$ particles (same resolution, different box size). The values of the statistics and their derivatives with respect to $`\overline{F}`$ and $`\kappa `$ with $`\gamma 1=0.3`$, $`T_{1.4}=17000`$ K, $`a=0.24`$, $`n=0.95`$, $`\alpha =0.0`$, and $`\overline{F}=0.7`$, with typical instrumental effects, are in Figs. 7, 8, and 9, respectively. The filters are noted on the plot.

As shown in Fig. 7, box size and resolution differences are generally within tolerance for $`P`$ and $`B`$ with the higher filter shown. Any small box size effects may be attributed to the difference in random seed in the 80 $`h^1`$Mpc boxes versus the 40 $`h^1`$Mpc boxes. There is a noticeable difference for $`B`$, the lower filter, between the higher and lower resolution 40 $`h^1`$Mpc boxes. This difference cannot be attributed to a seed difference since the seed is same for these boxes. The difference is not evident between the 40 and 80 $`h^1`$Mpc 5123 boxes, because of either a box-size effect or random fluctuations that oppose this resolution effect. For the most part, the $`k`$-dependent shapes of the derivatives agree, with the difference being one of magnitude. A smaller version of the Fisher analysis done in section 3 was completed for these boxes of different resolutions to check the effects of the derivative differences. The errors on amplitude are about 13% higher and on flux are about 15% lower with the lower resolution box than with the higher resolution box. However, since we are only looking at error bounds, which should hopefully be within about 25% accuracy, this error is tolerable, and the difference in Figure 7 is (by this criterion) not significant.

### 2.4 Random fluctuations due to finite box size

Another question is the magnitude of the random fluctuations due to finite box size, which can be determined by using different random seeds to generate the initial random Gaussian field. A study was done using three seeds, with $`\gamma 1=0.3`$, $`T_{1.4}=17000`$ K, $`a=0.24`$, $`n=0.95`$, $`\alpha =0.0`$, $`\overline{F}=0.68`$, and standard resolution and noise effects, on the values of $`P`$, $`B`$, and their $`\overline{F}`$\- and $`\kappa `$-derivatives. Plots of the relative difference between the statistics and their $`\overline{F}`$-derivatives computed with individual seeds and averaged over three seeds are in Figs. 10 and 11 respectively. A plot of the relative difference between the statistics computed with individual seeds and averaged over four seeds for a 40 $`h^1`$Mpc box is in Fig. 12.

As shown in these figures, the statistics and their derivatives for an 80 $`h^1`$Mpc box do not have very large fluctuations due to finite box size, so only one random seed (the same each time, for consistency) was used for this analysis rather than averaging over the statistics obtained from several seeds. To check concretely that this is not a problem, a Fisher analysis for just $`\kappa `$ (derivatives not shown, but they have a similar pattern as the flux derivatives) and $`\overline{F}`$ at one $`z`$ was completed using three different seeds, and using the spectra and covariance matrices computed by averaging over lines of sight from all three seeds. The $`\overline{F}`$ and $`\kappa `$ error bounds for individual seeds were within 0.015 (relative difference) of the bounds for the averaged spectra. Note that the 40 $`h^1`$Mpc box has significantly higher fluctuations, as expected, so it is possible that apparent box size effects are actually due to these fluctuations.

### 2.5 Chunking

The simulations were also used to study the effects of chunking, a procedure often used for data analysis. In the data, the spectra can be much longer than 80 $`h^1`$Mpc, and are cut into smaller chunks for analysis. The fields $`\delta _F(\lambda )`$ and $`H(\lambda )`$ are both created from the full spectrum in real space before chunking.

To identify the effects of chunking, we compared the values of the statistics for an 80 $`h^1`$Mpc box with $`512^3`$ particles; and the same 80 $`h^1`$Mpc box, computed with half-lines of sight and averaged (see results in Fig. 13). All other simulation parameters were identical: $`\gamma 1=0.3`$, $`T_{1.4}=17000`$ K, $`a=0.24`$, $`n=0.95`$, $`\alpha =0.0`$, $`\overline{F}=0.7`$, and typical noise and resolution effects.

As shown in Fig. 13, chunking has little effect on $`P`$ and $`B`$, with the values of the statistics computed using chunking almost perfectly matching those for the 80 $`h^1`$Mpc box without chunking. Thus, the chunks of spectrum as small as 40 $`h^1`$Mpc used in a data analysis can accurately represent the results in an infinite-spectrum limit.

### 2.6 Filters

As stated in section 2.1, we tried two filter shapes and several filter sizes and numbers of filters. The square window, in Eq. 6, was used in Zaldarriaga et al. (2001b). The Gaussian window was also tried since the Fourier transform of the square filter is not well-localized. The sizes of the filters were chosen based on the following consideration: the $`k`$ values in these 80 $`h^1`$Mpc boxes range from 0.078-5.0 $`h`$Mpc-1. Assuming $`\mathrm{\Omega }_m=0.3`$, at $`z=3`$, this range corresponds to 0.0007-0.05 s km-1. Below about 0.001 s km-1, continuum fluctuations become important, and above about 0.02 s km-1, resolution and noise effects become very important. Consequently, for square filters, the largest range of filter used is 0.002-0.02 s km-1 (0.2-2.0 $`h`$Mpc-1), and several smaller filters were chosen. For Gaussian filters, the question of filter width is a bit more difficult because there is some small contribution from beyond even $`2\sigma `$. Thus, we hoped to find the balance between the competing considerations of filter smoothness (compactness of the Fourier transform) versus its ability to eliminate $`k`$-values outside the range we want to consider. For the data, other possibly useful filters that can balance the smoothness and noise considerations are a triangular filter, or a Gaussian filter with a sharper low-$`k`$ cutoff. We show results for several sizes of square and Gaussian filters.

First we investigated how does an increase in number of filters with different sizes improve the information from $`B`$. A Fisher analysis using varying numbers of square filters was done to determine the number of filters to use in the range 0.2-2.0 $`h`$Mpc-1. We tried using one or two filters by subdividing that interval accordingly. No filters overlapped, as analysis showed that overlap did not increase the amount of information at all. The error bars computed from one filter versus from two filters in the same range were within 5% of each other. Thus one filter is sufficient for analysis of a given range of $`k`$. However, for this analysis we used two non-overlapping filters for the total range, allowing us to analyze separately the information contribution of lower and higher filter ranges.

We also investigated filter sizes and shapes as shown in Table 1. Filter B1 is the lower filter of the pair referred to as filter B, and B2 is the higher filter. The numbers shown are actually for $`z=3`$; the filter parameters are constant in s km-1, and so their value in $`h`$Mpc-1 is $`z`$-dependent. While the magnitude of the statistics are rather filter-dependent, the shape of the statistics is essentially filter-independent and so the results were almost independent of the filter size and shape used, at least within the range of filters explored here.

### 2.7 Instrumental dependencies

Noise and resolution play an important role in determining the sensitivity of a given data set. We parameterize the noise with two numbers, $`f_c`$ (the fraction of noise whose amplitude is correlated with the signal) and $`S/N`$ (the overall signal to noise ratio relative to the continuum, including all noise). Noise whose amplitude is correlated with the signal is characterized by its variance $`\sigma _c^2=c\overline{F}`$ (where $`\overline{F}`$ is the mean flux from the quasar). A typical example is photon shot noise, where the variance is proportional to the flux. Noise whose amplitude is uncorrelated with the signal is characterized by its variance $`\sigma _u^2`$ (for example, Gaussian readout noise and photon shot noise from the sky) These variables are related as follows, defining a unique $`c`$ and $`\sigma _u^2`$ for each $`f_c`$ and $`S/N`$:

$`f_c`$ $`={\displaystyle \frac{c\overline{F}}{c\overline{F}+\sigma _u^2}}`$ (9)
$`S/N`$ $`={\displaystyle \frac{1}{\sigma _f}}={\displaystyle \frac{1}{\sqrt{\sigma _u^2+c\overline{F}}}}`$ (10)
The distinction between correlated and uncorrelated noise amplitude is made because they have different effects on $`B`$ and $`T`$. Uncorrelated noise contributes to $`P`$ and $`T`$, but not to $`B`$, whereas correlated noise contributes to all three. The noise’s contribution to $`P`$ can be calculated analytically whether its amplitude is correlated or uncorrelated with the signal. The correlated noise contribution to $`B`$ depends not only on the noise amplitude, but also on the signal power spectrum $`P`$.

In this study, we tried to match the pixelization, resolution, and noise to that of the SDSS data. The pixels in SDSS data are roughly 1 Å wide (so the smoothing $`\sigma `$ is approximately the size of one pixel). Assuming $`z=3`$ and $`\mathrm{\Omega }_m=0.3`$ this corresponds to 0.6 $`h^1`$Mpc. To match the 0.6 $`h^1`$Mpc pixels, we used pixels of size 0.625 $`h^1`$Mpc in the simulations. This small difference from the real value makes little effect in the value of the statistics, since even reducing the pixel size by a factor of four changes $`B`$ and $`P`$ by at most 5%. We used Gaussian smoothing $`\sigma `$ the same size as our pixels. SDSS data have $`S/N`$ (averaged over redshift bins in our range) ranging from 1 to 30, with median values in each redshift bin ranging from 4 to 5. We used the value $`S/N=5.0`$ for all redshift bins. The noise amplitude was entirely uncorrelated with the signal ($`f_c=0`$), since we could not evaluate the amount of correlated noise in the data.

## 3 Fisher Matrix Calculation and Results

The construction of the Fisher matrices for each individual redshift bin is straightforward. A vector containing all values $`x_i`$ that were to be used (for example, $`P`$ and $`B`$ for all $`k`$) was created, as was its full covariance matrix (computed using all 12288 lines of sight) and its derivatives with respect to the 6 parameters $`\gamma `$, $`T_{1.4}`$, $`\kappa `$, $`n`$, $`\alpha `$, and $`\overline{F}`$. The derivatives with respect to $`P`$ were created from the values with the power due to noise subtracted off. These covariance matrices and derivatives were then used to create the $`6\times 6`$ Fisher matrices for each redshift bin by computing

$$F_{kl}=\left(\frac{\mathrm{d}𝐱}{\mathrm{d}y_k}\right)^T𝖢^1\left(\frac{\mathrm{d}𝐱}{\mathrm{d}y_l}\right)$$
(11)
so that $`𝖥`$ is the inverse covariance matrix for $`𝐲`$, the vector of all the parameters. In general, there is another term related to $`\mathrm{d}𝖢/\mathrm{d}y_i`$ that we neglect in the expectation that it is small.

Next, a block-diagonal $`54\times 54`$ Fisher matrix was created, using the $`6\times 6`$ Fisher matrices for the nine redshift bins. Each block was multiplied by the length of spectrum in that redshift bin in the SDSS data; the total spectrum length is $`7.7\times 10^7`$ km/s. Fisher matrices in this form were created for all seven sets of filters listed in Table 1, for various combinations of statistics. Further information about error-bounds is extracted from these large Fisher matrices by projecting from them down to a smaller number of parameters.

### 3.1 Projection to fewer parameters

To reduce to the smallest number of parameters, we used the following procedure. The $`54\times 54`$ Fisher matrices are the inverse covariance matrices $`𝖢^1(𝐲)`$ for the values of the parameters in each redshift bin. We projected down to a smaller set of parameters that includes:

* $`\gamma _0`$ and $`\gamma _1`$, the constant and linear terms of a linear expansion of $`\gamma (z)`$ about $`z=3.0`$,
* $`T_{1.4}`$, an overall redshift-independent value,
* $`\kappa _{\mathrm{eff}}`$, an overall effective amplitude at some $`z`$ and $`k`$,
* $`n_{\mathrm{eff}}`$, an overall effective slope, where the evolution of the pivot point (fixed in velocity coordinates) is expressed via $`n(z)=n(k_p)+\alpha \mathrm{ln}(k_p(z)/k_p)`$,
* $`\alpha _{\mathrm{eff}}`$, an overall effective curvature, and
* The nine values of $`\overline{F}`$, one per redshift bin.

This projection includes 15 parameters, and can be computed using Eq. 11.

Here we present results for the analysis done with just filter B, using $`P`$ and $`B`$. All references to “error bounds” implies that we have marginalized over other parameters (by using $`\sqrt{(F^1)_{ii}}`$ rather than $`1/\sqrt{F_{ii}}`$).

First, the error bounds for $`\gamma _0`$ and $`\gamma _1`$, defined by \[$`\gamma 1](z)=\gamma _0+\gamma _1(z3)`$, are 0.054 and 0.10, respectively. We are not very sensitive to $`T_{1.4}`$, with a relative error of 0.08.

Table 2 shows the relative errors for the nine values of $`\overline{F}`$. The errors increase with redshift mainly because the amount of data decreases, and at low redshifts they are extremely small.

The error on $`\kappa =(\mathrm{\Delta }\sqrt{A})/\sqrt{A}`$ is .017. The errors on $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ (defined at pivot point $`k=0.009`$ s km-1 at $`z=2.6`$) are $`\mathrm{\Delta }n_{\mathrm{eff}}=0.013`$, and $`\mathrm{\Delta }\alpha _{\mathrm{eff}}=0.029`$. The error ellipse in terms of $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ computed from their joint covariance matrix is in Fig. 14. To answer the question of whether $`B`$ alone provides independent information competitive wih $`P`$, or whether $`B`$ and $`P`$ together is more sensitive than either of the two alone, we carried out the analysis using $`P`$ alone, $`B`$ alone, and $`P`$ and $`B`$ together. The error bounds for all parameters for various combinations of statistics can be found in Table 3.

Fig. 14 shows the 68% confidence level error ellipses for $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ for $`P`$, $`B`$, and both together. It is clear that for this combination of parameters $`P`$ provides most of the information. While $`B`$ can act as a useful independent confirmation, adding its information by itself does not decrease the statistical errors significantly.

Fig. 15 shows the 68% CL error ellipses for $`\kappa _{\mathrm{eff}}`$ and $`\overline{F}(z=3.0)`$ for $`P`$, $`B`$, and both together. In contrast with Fig. 14, combining the two significantly improves upon $`P`$ or $`B`$ alone, reducing the errors on $`\kappa `$ and $`\overline{F}`$ by almost a factor of 3. The primary reason for this is the partial degeneracy between the two parameters in $`P`$ alone, as shown in Fig. 15. This degeneracy is easy to understand from the derivatives in Figs. 2 and 3. Changing $`\kappa `$ or $`\overline{F}`$ has a relatively flat effect on the power spectrum $`P`$, so the two parameters are significantly degenerate from $`P`$ alone, especially when all the other parameters are included. In contrast, the effect on $`B`$ is much larger for $`\kappa `$ than for $`\overline{F}`$. Combining $`P`$ and $`B`$ thus allows one to break the degeneracy between the two parameters and determine both with a much higher precision.

We also investigated how the errors vary with filter, finding that they did not vary significantly for filters given in Table 1. In particular, the error bounds on $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ were virtually independent of filter, and those on $`\kappa _{\mathrm{eff}}`$ were nearly so. Among the Gaussian filters, there was a slight trend towards lower error bounds on $`\kappa _{\mathrm{eff}}`$ for the larger filters (0.015 for the largest versus 0.018 for the smallest) but even here the difference is not very large.

### 3.2 Expansions of $`\kappa _{\mathrm{eff}}`$, $`n_{\mathrm{eff}}`$, and $`\alpha _{\mathrm{eff}}`$

In order to determine how much variation in the values of $`n_{\mathrm{eff}}`$, $`\alpha _{\mathrm{eff}}`$, and $`\kappa _{\mathrm{eff}}`$ (at a particular $`z`$ and $`k`$) is induced by deviations in $`n`$, $`\alpha `$, $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_b`$, and $`h`$, we used CMBFAST (Seljak & Zaldarriaga, 1996) outputs to study the derivatives with respect to these parameters, using spectra that were normalized to the same amplitude at $`k=0.05`$ $`h`$Mpc-1 at $`z=20`$. The derivatives were computed at $`z=2.6`$, $`k=0.009`$ s km-1, the pivot point for SDSS data sample. We assume the fiducial model $`\mathrm{\Omega }_m=0.3`$, $`h=0.7`$, $`\mathrm{\Omega }_b=0.04`$ (or $`x\mathrm{\Omega }_b/\mathrm{\Omega }_m=0.13`$). The pivot point for $`\alpha `$ is $`k=0.05`$ $`h`$Mpc-1. The expansion for $`\kappa _{\mathrm{eff}}`$ is denoted $`\mathrm{\Delta }\kappa _{\mathrm{eff}}`$ to emphasize that this expansion is to be used to calculate what relative change in power spectrum amplitude would result from changing the cosmological parameters from our chosen values (hence it gives $`\mathrm{\Delta }\kappa _{\mathrm{eff}}=0`$ for our cosmological model). The resulting expansions are

$`\mathrm{\Delta }\kappa _{\mathrm{eff}}`$ $`=\mathrm{\Delta }\kappa +1.5(n0.75)0.032\alpha `$ (12)
$`+0.28\mathrm{ln}\left({\displaystyle \frac{\mathrm{\Omega }_m}{0.3}}\right)+0.65\mathrm{ln}\left({\displaystyle \frac{h}{0.7}}\right)+0.14\mathrm{ln}\left({\displaystyle \frac{x}{0.13}}\right)`$
$`n_{\mathrm{eff}}`$ $`=3.25+n+3.0\alpha +`$ (13)
$`0.16\mathrm{ln}\left({\displaystyle \frac{\mathrm{\Omega }_m}{0.3}}\right)+0.23\mathrm{ln}\left({\displaystyle \frac{h}{0.7}}\right)+0.034\mathrm{ln}\left({\displaystyle \frac{x}{0.13}}\right)`$
$`\alpha _{\mathrm{eff}}`$ $`=0.22+\alpha 0.11\mathrm{ln}\left({\displaystyle \frac{\mathrm{\Omega }_m}{0.3}}\right)`$ (14)
$`0.17\mathrm{ln}\left({\displaystyle \frac{h}{0.7}}\right)+0.034\mathrm{ln}\left({\displaystyle \frac{x}{0.13}}\right)`$
The effect on the amplitude is primarily determined by the primordial amplitude and slope, but not the derivative of the slope (given that the pivot point $`k=0.05`$ $`h`$Mpc-1 is at much lower $`k`$ than the wavevectors the forest is sensitive to). Today’s matter density $`\mathrm{\Omega }_m`$ and Hubble constant $`h=H_0/100`$ $`\mathrm{km}/\mathrm{s}/\mathrm{Mpc}`$ change the transfer function and so affect the amplitude, slope and derivative of the slope. However, $`\mathrm{\Omega }_m`$ also changes the relation between velocities and comoving coordinates as given by the Hubble parameter $`H(z)`$. This partially cancels the transfer function effect of $`\mathrm{\Omega }_m`$, so that the coefficients in front of $`\mathrm{\Omega }_m`$ are reduced relative to $`h`$. Baryon density has a smaller effect than other parameters, at least around the small value of $`\mathrm{\Omega }_b/\mathrm{\Omega }_m`$ we expanded here.

It is clear that since we can only determine precisely two numbers, the amplitude and the slope, and somewhat less precisely also the derivative of the slope, we cannot determine all of the parameters with high precision. Hence it is best if the Ly-$`\alpha `$ forest is combined with other cosmological probes, such as cosmic microwave background (CMB) anisotropies. Since one expects CMB to provide very strong constraints on $`h`$ and $`\mathrm{\Omega }_m`$ (at least in the context of spatially flat models), this then allows the Ly-$`\alpha `$ forest to determine the primordial slope and its derivative with high accuracy.

## 4 Power Spectrum Evolution and Quintessence

In this section we try to determine the sensitivity of the Lyman-$`\alpha `$ forest to deviations from the Einstein-de Sitter (EdS) growth factor $`D(a)=a`$. We analyzed both specific quintessence models and also more general deviations of the growth factor from EdS. Sensitivity to quintessence model parameters can arise from the change in the growth factor $`D(z)`$ and from the Hubble parameter $`H(z)`$. In our tests the latter was subdominant and most of the sensitivity was from the growth factor.

First, we studied the sensitivity of the Lyman-$`\alpha `$ forest to quintessence parameters, $`\mathrm{\Omega }_q`$ and the equation of state $`w_q`$, for nine different quintessence models. Several models considered had static equation of state $`w_q=w_0`$, while the other models were dynamic, $`w_q=w_0+w_1(a1)`$. The models studied are shown in Table 4.

All models had $`\mathrm{\Omega }_b=0.04`$ and $`h=0.7`$, with $`\mathrm{\Omega }_{\mathrm{CDM}}`$ determined by requiring $`\mathrm{\Omega }_b+\mathrm{\Omega }_{\mathrm{CDM}}+\mathrm{\Omega }_q=1`$.

### 4.1 Projection with quintessence parameters

To find constraints on all parameters once quintessence parameters are included, we started from the $`54\times 54`$ Fisher matrices created in section 3, and projected down to the parameters in section 3.1 plus $`\mathrm{\Omega }_q`$ and $`w_q`$ (or $`\mathrm{\Omega }_q`$, $`w_0`$, and $`w_1`$ for dynamical models). All error bounds are computed from the Fisher matrix for filter B using $`P`$ and $`B`$. To do so, we computed derivatives of $`\kappa _{\mathrm{eff}}`$, $`n_{\mathrm{eff}}`$, and $`\alpha _{\mathrm{eff}}`$ with respect to quintessence parameters at some pivot point fixed in s km-1. These derivatives were computed from CMBFAST outputs for the 9 $`z`$-bins for all 9 quintessence models studied, normalized to have the same amplitude for $`z=20`$ at the same pivot point at which derivatives were taken.

Results for all models indicated that $`w_q`$ and $`\mathrm{\Omega }_q`$ are highly degenerate, with a correlation coefficient of -0.9. The well-determined direction for model 1, which is fairly representative, in $`(w_q,\mathrm{\Omega }_q)`$ space is $`(0.60,0.80)`$. The strong correlation is not surprising, since the redshift range probed is too small to determine two parameters from the growth factor or $`H(z)`$ evolution. Some information on $`\mathrm{\Omega }_q`$ also comes from the power spectrum shape, which is why the correlation is not perfect. In the following we thus fix one parameter, for example $`\mathrm{\Omega }_q`$, and present errors on equation of state only. The motivation for this is that other tests, most notably CMB combined with large scale structure tests (e.g. cluster counts, galaxy clustering, weak lensing), will be able to determine $`\mathrm{\Omega }_q`$ at $`z=0`$ very accurately. If these tests also provide independent constraints on $`w_0`$ then for the time dependent $`w`$ one can use our results to constrain $`w_1`$.

The errors for fixed $`w_q`$ (marginalized over all other variables) can be found in Table 5. The 68% CL error ellipses show how $`\mathrm{\Omega }_q`$ and amplitude are degenerate (Fig. 16). Analysis of the covariance matrices for these model parameters (or the fact that the error bounds for $`\kappa _{\mathrm{eff}}`$, $`n_{\mathrm{eff}}`$, and $`\alpha _{\mathrm{eff}}`$ have risen from their value computed without including quintessence) indicated a very high degree of degeneracy between $`\mathrm{\Omega }_q`$ and $`\kappa _{\mathrm{eff}}`$, $`n_{\mathrm{eff}}`$, and $`\alpha _{\mathrm{eff}}`$. If $`\alpha `$ is fixed then the errors on $`\mathrm{\Omega }_q`$ consequently decrease (though this effect is much greater for the models to which we are least sensitive). This indicates that at least some of our information about quintessence parameters comes from shape information. The degeneracy with $`\kappa _{\mathrm{eff}}`$ indicates that some information is from the growth factor information as well, since the magnitude of the growth factor (as opposed to its shape) is degenerate with $`\kappa _{\mathrm{eff}}`$. The errors with $`w`$ and $`\alpha `$ fixed are in Table 6.

Error bounds on $`w_q`$ for static models, with $`\mathrm{\Omega }_q`$ fixed (as may be done by measurements other than the Lyman-$`\alpha `$ forest), are shown in Table 7. As shown there, now that $`\mathrm{\Omega }_q`$ is fixed, $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ have error bounds roughly equal to those computed in the previous section with no quintessence. This fact shows that these parameters contain information about $`\mathrm{\Omega }_q`$ but not $`w_q`$. $`\kappa _{\mathrm{eff}}`$, on the other hand, is degenerate with $`w_q`$, indicating that $`w_q`$ is detected via its influence on the growth factor. Consequently, the model with least negative $`w_q`$ has the best error bounds. Table 8 shows error bounds for $`w_1`$ with both $`\mathrm{\Omega }_q`$ and $`w_0`$ fixed. Errors for $`n_{\mathrm{eff}}`$ and $`\alpha _{\mathrm{eff}}`$ are not shown since they are the same as in the previous section, without quintessence. Not surprisingly, the error bounds are best for the model with the most rapidly-varying equation of state.

### 4.2 More general growth

Besides considering specific quintessence models, we also considered more general power spectrum evolution, modeled as

$$D(a)=D_0(a/a_0)^s,$$
(15)
where $`s=1`$ corresponds to an Einstein-de Sitter universe. Thus, we can determine if our sensitivity to quintessence is due to growth factor deviations. This analysis was completed by starting with the $`54\times 54`$ Fisher matrix for $`P`$ and $`B`$ (filter B), and projecting down to the same set of parameters as in Section 3.1 along with $`s`$. The projection for $`\kappa `$ and $`s`$ used $`Pe^{2\kappa }D^2(a)`$ to relate the necessary derivatives $`\mathrm{d}P/\mathrm{d}\kappa `$ and $`\mathrm{d}P/\mathrm{d}s`$ to the computed values $`\mathrm{d}P/\mathrm{d}a`$. Using the value of $`a_0`$ determined by requiring that amplitude and $`s`$ have the minimum degeneracy, and inputting a value of $`s=1`$, yields an error $`\mathrm{\Delta }s=0.11`$. Note that if only information from the power spectrum $`P`$ is used (rather than including $`B`$ as done here), the error rises to $`\mathrm{\Delta }s=0.47`$.

Values of $`\mathrm{\Delta }s`$ for arbitrary static quintessence models were computed using an integral expression for $`D(a)`$ obtained from Hamilton (2001), Taylor expanding to several orders, which gave a form $`D(a)=a\left(1+_{n=1}c_n(\mathrm{\Omega }_q,w_q)a^{3nw_q}\right)`$. Expanding this expression about $`a_0`$ gave an expression for $`s`$, and hence for $`\mathrm{\Delta }s=s1`$, in terms of the model parameters. Requiring $`|\mathrm{\Delta }s|<0.11`$ defines the region in $`w_q`$-$`\mathrm{\Omega }_q`$ parameter space that is detectable using the Lyman-$`\alpha `$ forest (growth factor information only, not curvature). This curve is shown in Fig. 17.

To check the validity of this expansion, the growth factor was computed directly from CMBFAST outputs for the 9 quintessence models studied above. This also allowed the computation of $`\mathrm{\Delta }s`$ for dynamic models, for which an analytic expression was not obtained. Values of $`\mathrm{\Delta }s`$ for the 9 models studied earlier, and values of $`\mathrm{\Omega }_q(z=2.6)`$, are shown in Table 9. As shown, $`\mathrm{\Delta }s`$ and $`\mathrm{\Omega }_q(z=2.6)`$ are highly correlated.

### 4.3 Change in $`H(z)`$

So far the only tests of dark energy proposed have used either the growth factor or the redshift luminosity distance. If one wishes to extract $`w(z)`$, then both of these involve a double integral over this quantity and degeneracies arise. A more direct and still in principle observable way is to measure the Hubble parameter $`H(z)`$, which is related to $`w(z)`$ through a single integral. One way to measure it is to have a characteristic feature fixed in comoving coordinates which is observed in redshift space. The relation between redshift space and comoving space is determined by $`H(z)`$ and so observing the feature as a function of redshift determines its evolution. The problem of course is that there are no characteristic features imprinted, since the distribution of structure in the universe is stochastic in nature.

One must therefore look for a characteristic scale in correlations between structures. In principle such features could be provided by baryonic oscillations imprinted in the matter power spectrum, but in practice this is a weak effect limited to very large scales and so cannot be made very precise. One is thus left with the variations in the correlation function slope as a function of scale. The slope varies from $`n1`$ on large scales to $`n3`$; on Ly-$`\alpha `$ forest scales we find $`\alpha _{\mathrm{eff}}=0.22`$ in CDM models. Hence the scale at which the slope takes a specific value can be viewed as a standard ruler and can be traced with redshift. If this slope is measured in redshift space, then one is measuring directly $`H(z)`$. To be able to do this one has to detect the curvature in the slope over the dynamic range of observations. This is challenging, since the dynamic range is narrow and the error on the slope will be large. We find that such a detection should be possible with the current sample, but is not expected to improve the constraints on $`w(z)`$ significantly.

## 5 Discussion and Conclusions

Our results show that the statistical precision one can achieve from Ly-$`\alpha `$ forest using SDSS data is truly impressive. The overall amplitude of the power spectrum, its slope and its curvature can all be determined with 1-3% precision. Our results are significantly more optimistic than those in Zaldarriaga et al. (2001a), which focused on a higher resolution but smaller data sample available from Keck. This is primarily a consequence of the large sample size of SDSS data, which increases the statistical precision on larger scales and moves the pivot point from $`0.03`$ to $`0.01`$ s km-1. On larger scales the sensitivity of the power spectrum to the gas temperature-density relation is reduced so one is measuring more directly the underlying matter power spectrum. In addition, we have shown here that higher order correlations break the degeneracy between the mean flux and the amplitude (this degeneracy is broken to some extent already by the power spectrum alone). Using the power spectrum alone and the bispectrum alone are two independent ways to confirm this determination. Using them together gives even lower error bounds due to the breaking of degeneracies. The statistical precision of such a data set is competitive with the one from the CMB using the MAP and Planck satellites. More importantly, combining the two probes will allow one to determine the primordial power spectrum to a high precision over 3 decades in scale ($`0.001<k<2`$ $`h`$Mpc-1).

Furthermore, since one is measuring the power spectrum over a range of redshifts ($`2<z<4`$), one can also study the growth factor evolution. While for standard cosmological constant models the universe is effectively Einstein-de Sitter for $`z>2`$, dynamical models with rolling scalar fields often produce an equation of state increasing with redshift (this is especially true for many of the so-called tracker models). In this case the dark energy or quintessence is still dynamically important for $`z>2`$ and can affect the growth rate of perturbations and the Hubble parameter $`H(z)`$. While there is no simple single parameter combination that describes the sensitivity to dark energy, it is clear that the precision is correlated with $`\mathrm{\Omega }_q`$ at $`z>2`$. Our results show that if $`\mathrm{\Omega }_q(z=2.6)>0.2`$ then the deviations in the growth factor are sufficiently large to be detected in Ly-$`\alpha `$ forest spectra using the current SDSS sample. With the full SDSS sample this limit can be improved further and models with $`\mathrm{\Omega }_q(z=2.6)>0.1`$ should be detectable. This sensitivity can increase further if one can extend the available data set to even lower redshifts by using space based observations.

While it is clear that the statistical power of upcoming data sets is impressive, the main remaining issue is whether the simplest model of the QSO absorption adopted here is valid, or whether there are other more complicated astrophysical processes that can spoil the picture. Since the expected statistical precision is so high one must investigate processes that affect the statistics even at a very small (1%) level. There are several possible complications that may have an effect and we mention some here: metal line absorption along the line of sight (e.g. C-IV, Si-III…), an inhomogeneous UV background, kinematic gas outflows generated by supernovae from the galaxies, temperature fluctuations, etc.

Fortunately, while there are several possible contaminations, there are also several possible tests one can apply to identify and remove the contamination. Some of these processes can be investigated with the statistics used in this paper, such as the power spectrum and bispectrum. In the simplest model the main driver of correlations is gravity, which induces a very specific relation between low and higher order statistics. Any astrophysical contamination is likely to destroy these relations. Other tests such as the probability distribution of the flux (as a function of scale) will also yield useful information (McDonald et al., 2000). Another way to investigate these effects is through the cross-correlation between the forest and galaxies at the same redshift (Adelberger et al., 2003); for theoretical attempts to interpret this result, see McDonald et al. (2002), Croft et al. (2002a), Kollmeier et al. (2002) and Bruscoli et al. (2002).

If the simplest model passes this and other high precision tests it will receive an important confirmation that will strengthen the credibility of the results. In the opposite case one can still identify the contaminant and apply the corrections assuming the contamination is sufficiently well understood. This parallels the situation in the CMB, where secondary processes and foregrounds are searched for and, if identified, subtracted from the primary CMB. It is clear however that a validation of the interpretation set forth by the simplest model will require a lot of coordinated effort from several groups and more effort should be put into this. Our results suggest that this is well worth the effort and that the Ly-$`\alpha `$ forest could be our next cosmological gold-mine.

R.M. is supported by an NSF Graduate Research Fellowship. This work was supported by NASA, NSF CAREER, David and Lucille Packard Foundation and Alfred P. Sloan Foundation. We thank Nick Gnedin for the HPM code. The simulations were performed at the National Center for Supercomputing Applications.


format: text

text:
# Roche tomography of cataclysmic variables – II. Images of the secondary stars in AM Her, QQ Vul, IP Peg and HU Aqr

## 1 Introduction

The secondary, Roche-lobe-filling stars in CVs are key to our understanding of the origin, evolution and behaviour of this class of interacting binary. To best study the secondary stars in CVs, we would ideally like direct images of the stellar surface. This is currently impossible, however, as typical CV secondary stars have radii of 400 000 km and distances of 200 pc, which means that to detect a feature covering 20 per cent of the star’s surface requires a resolution of approximately 1 microarcsecond, 10 000 times greater than the diffraction-limited resolution of the world’s largest telescopes. Rutten & Dhillon (1994) and Watson & Dhillon (2001, hereafter referred to as Paper I) described a way around this problem using an indirect imaging technique called Roche tomography which uses phase-resolved spectra to reconstruct the line intensity distribution on the surface of the secondary star.

Obtaining surface images of the secondary star in CVs has far-reaching implications. For example, a knowledge of the irradiation pattern on the inner hemisphere of the secondary star in CVs is essential if one is to calculate stellar masses accurately enough to test binary star evolution models (see Smith & Dhillon 1998). Furthermore, the irradiation pattern provides information on the geometry of the accreting structures around the white dwarf (see Smith 1995).

In this paper we present new Roche tomograms of the secondary star in the dwarf nova IP Peg, and the magnetic CVs AM Her and QQ Vul in the light of the Na I $`\lambda \lambda `$8183,8195 Å absorption doublet. In addition, we also present a Roche tomogram of the magnetic CV HU Aqr in the light of the He II $`\lambda `$4686 Å emission line component known to originate from the secondary star. These tomograms allow a study of the irradiation pattern on the secondary stars as well as providing measurements of the binary parameters for all four CVs.

## 2 Observations and reduction

The spectra of QQ Vul and HU Aqr were taken on the 3.5-m telescope at Calar Alto and full details of the observations and subsequent data reduction can be found in Catalán, Schwope & Smith (1999) and Schwope, Mantel & Horne (1997), respectively.

The observations of IP Peg and AM Her were carried out over two nights on 1994 August 11–12 using the ISIS dual-beam spectrograph and 1200 lines mm-1 gratings on the 4.2m William Herschel Telescope (WHT). The EEV3 CCD chip in the red arm covered the wavelength range 7985–8437Å at a resolution of 0.37Å ($``$27 km s-1). The TEK2 CCD chip on the blue arm covered the wavelength range 4560–4968Å at a resolution of 0.4Å ($``$50 km s-1). To reduce readout noise, the pixels were binned by a factor of 2 in the spatial direction. Each object was observed using exposure times of 300 seconds and a nearby field star was located on the slit to correct for slit losses. Wide slit spectra of the companion star were also taken to permit absolute spectrophotometry. Sky-flats and tungsten flats were taken on both nights and a CuAr+CuNe arc spectrum was taken approximately every 50 minutes in order to wavelength calibrate the spectra. Narrow and wide slit spectra at a range of different air-masses were taken of flux standards for the purposes of flux calibration and water removal. We also observed a range of M-dwarf spectral-type templates in order to obtain a Na I intrinsic line profile for use in Roche tomography.

During the reduction process, bias subtraction, division by a master flat-field frame and sky subtraction were carried out and then the spectra were optimally extracted (Horne, 1986). The arc spectra were used to wavelength calibrate the data, with an rms error of 0.01Å. Slit loss correction was then applied by dividing the object spectra by the average value per pixel (after masking out strong lines, cosmic rays and telluric features) in the corresponding comparison star spectra, and later multiplying by the average value per pixel in the wide-slit comparison star spectrum. Telluric absorption and absolute flux calibration were then performed using the standard star spectra.

In total we obtained 79 spectra of IP Peg and 51 spectra of AM Her. We present only the red spectra in this paper, as we are only concerned with the Na I doublet at $`\lambda \lambda `$8183,8195Å for the purposes of this work. A full journal of observations is presented in Table LABEL:tab:journal.

## 3 Roche tomography

### 3.1 Techniques

Roche tomography (Rutten & Dhillon 1994; Dhillon & Watson 2001; Paper I) is analogous to the Doppler imaging technique used to map rapidly rotating single stars (e.g. Vogt & Penrod 1983). In Roche tomography the secondary star is assumed to be Roche-lobe-filling, locked in synchronous rotation and to have a circularized orbit. The secondary star is then modelled as a series of quadrilateral tiles or surface elements, each of which is assigned a copy of the local (intrinsic) specific intensity profile. These profiles are then scaled to take into account the projected area of the surface element, limb darkening and obscuration, and then Doppler shifted according to the radial velocity of the surface element at that particular phase. Summing up the contributions from each element gives the rotationally broadened profile at that particular orbital phase.

By iteratively varying the strengths of the profile contributed from each element, the ‘inverse’ of the above procedure can be performed. Due to the variable and unknown contribution to the spectrum of the accretion regions in CVs, the data must be slit-loss corrected and continuum subtracted prior to mapping with Roche tomography. Thus, Roche tomograms present images of the distribution of line flux on the secondary star. The data is fit to an aim $`\chi ^2`$, and a unique map is selected by employing the maximum entropy (MEMSYS) algorithm developed by Skilling & Bryan (1984). The definition of image entropy, $`S`$, that we use is

$$S=\underset{j=1}{\overset{k}{}}m_jd_jm_j\mathrm{ln}\left(m_j/d_j\right)$$

where $`k`$ is the number of tiles in the map, and $`m_j`$ and $`d_j`$ are the map of line flux and the default map, respectively. We employ a moving uniform default map, where each element in the default is set to the average value in the map. For further details of the Roche tomography process see the review by Dhillon & Watson (2001).

For those CVs where we wish to map the Na I doublet, Roche tomography requires the intrinsic line profile to be input as a template. Ideally, this should be the same line taken from the spectrum of a slowly rotating star of the same spectral class as the star being mapped, preferably obtained with the same instrumental set-up. Unfortunately, only one of the spectral-type template stars that were observed is useful (GL806 – see Table LABEL:tab:templates); the others show signs of broadening of the Na I doublet and have been identified as either flare stars, and are therefore probably rapidly rotating, or binary stars, where $`V_{rot}\mathrm{sin}i`$ measurements using the Na I absorption doublet have been found to be systematically greater than those given by the narrow metal lines (Bleach et al. 2000). We used the GL806 template in all of our reconstructions and, even though the spectral class of the template differs from those of the stars being mapped, the exact shape of the template profile used has little effect on the reconstructions due to the large degree of rotational broadening present in CV secondary stars.

For the CVs mapped using the Na I doublet, it is necessary to correct the derived systemic velocities for the systemic velocity of the spectral-type template used in the reconstructions. The systemic velocity of GL806 was determined by fitting gaussians to the Na I doublet, which yielded a value of $`\gamma `$ = –24.3 $`\pm `$ 1.1 km s-1. For comparison, Wilson (1953) found the systemic velocity of GL806 to be –15 $`\pm `$ 5 km s-1, whereas Gliese (1969) quote a value of –24.4 km s-1.

It is important to note that the maximum-entropy algorithm we use requires that the input data are positive. As a result, it is necessary to invert absorption-line profiles and, in the case of low signal-to-noise data, add a positive constant in order to ensure that there are no negative values in the continuum-subtracted data. If this positive constant is not added, and any negative points are either set to zero or ignored during the reconstruction, then the fit will be positively biased in the wings of the line profile wings. The exact value of this constant is unimportant, so long as it prevents the occurrence of negative values in the data and is effectively removed. We account for this positive constant during the reconstruction by employing a virtual image element which contributes a single value to all data points, effectively cancelling out the constant. For the systems where the Na I doublet has been used in the reconstructions, we have employed this technique of adding a positive constant to each data value. We find that this significantly improves the reliability of the technique to determine the component masses when noisy data are being used. As the data for HU Aqr are from Gaussian fits to the He II $`\lambda `$4686Å emission line, there are no errors in the continuum and hence we do not apply this technique.

In addition, we have accounted for the effects of velocity smearing due to the motion of the secondary star over the duration of an exposure for all four systems. This is done by calculating the profile at evenly spaced intervals over one exposure and then averaging them. Although computationally more time-consuming, this reduces the effects of smearing, and can also take into account features that may cross the limb of the secondary star during an exposure.

### 3.2 Determining the system parameters

Roche tomograms are greatly influenced by parameters such as the component masses, inclination and systemic velocity of the binary. Errors in the assumed values of these parameters degrades the quality of the surface map and, in general, results in additional structure in the Roche tomograms. The correct parameters are, therefore, those that produce the map containing the least artefacts, corresponding to the map of highest entropy (Paper I). By carrying out reconstructions for many pairs of component masses (iterating to the same $`\chi ^2`$ on each occasion) we can construct an ‘entropy landscape’ (e.g., Fig. 3), where each point corresponds to the entropy value obtained in a reconstruction for a particular pair of component masses.

In order to determine the binary parameters of the CVs studied in this paper, we constructed a series of entropy landscapes assuming a wide range of different values for the orbital inclination and systemic velocity. For each combination of $`i`$ and $`\gamma `$, we picked the maximum entropy value in the corresponding entropy landscape, and plotted this value in the $`i\gamma `$ plane, which is shown for each object on the left-hand side of Fig. 1. In all cases we found that a unique or optimal systemic velocity could be selected that consistently gave the map of maximum entropy. In addition, the value of the optimal systemic velocity was found to be largely independent of the orbital inclination assumed during the reconstructions, as demonstrated in the left-hand plots of Fig. 1. For example, in reconstructions carried out over a range of inclinations spanning 30, the systemic velocity we obtain for both IP Peg and AM Her only varies by 1 km s-1.

The inclination, however, is not as well constrained for any of the CVs. The right-hand plots on Fig. 1 show cuts through the plots in the $`i\gamma `$ plane at the optimal systemic velocity, and the scatter clearly demonstrates the difficulty in selecting an optimal inclination, although the deterioration in the quality of the reconstruction for some inclinations allow limits to be set.

The difficulty in constraining the inclination can be explained by considering that there are two ways of determining the inclination of non-eclipsing CVs from the spectra of the secondary star alone. First, from the variation in the projected radius of the Roche-lobe shaped secondary star as it is viewed at different aspects, corresponding to a variation in the shape (or measured $`v\mathrm{sin}i`$) of the line profiles (see Shahbaz 1998). The second is from the variation in the strength of the profiles, since a high inclination system will exhibit a much greater variation in the line profile strength than a low inclination system with the same surface intensity distribution.

Unfortunately, the variation in $`v\mathrm{sin}i`$ is expected to be, at most, around 20–25 km s-1 for high inclination systems. Although this is comparable to the velocity resolution of our data, it is buried in noise and hence the mapping technique is largely blind to any such variations. Our only constraint on the inclination, therefore, comes from the variation of the line strength, which also has its limitations. Although, rapidly varying features can only be mapped onto high inclination systems, slowly varying features produced by a low inclination system can also be mapped onto the polar regions of a high inclination system with little loss in the quality of the reconstruction. Thus we only expect to constrain the inclination in CVs which have high inclinations and exhibit strong asymmetries in their surface intensity distributions. In addition, for the same reason outlined above, we also expect our determinations to be biased towards higher inclinations.

As we cannot reliably constrain the inclination of the non-eclipsing CVs in this work, and hence the component masses, Fig. 2 shows the derived component masses as a function of inclination. The top section of each panel in Fig. 2 also shows the mass ratio derived for each inclination. Since $`v\mathrm{sin}i`$ is purely a function of the mass ratio $`q`$ (=$`M_1/M_2`$) and the radial velocity semi-amplitude of the secondary star $`K_r`$, reconstructions over a range of inclinations should yield the same value of the mass ratio on each occasion. This is confirmed in Fig. 2, with the scatter in the values serving as an indication of the noise in the technique. Indeed, the scatter can be accounted for by the fact that the component masses are only varied in 0.02M increments in the entropy landscapes.

Entropy landscapes representative of the mass estimates for each system are shown in Fig. 3, and a summary of the adopted system parameters can be found in Table LABEL:tab:estimates. The results for each CV are discussed in more detail in the relevant sections below.

## 4 Results for AM Her

### 4.1 The binary geometry of AM Her

From Fig. 1 the maximum entropy value yields a value of –17 km s-1 for the systemic velocity of AM Her, taking into account the systemic velocity of the template star. This is within 2$`\sigma `$ of the range of values found by Young, Schneider & Shectman (1981) of –14 $`\pm `$ 4 km s-1 and –40 $`\pm `$ 15 km s-1 from circular fits to the H$`\gamma `$ and He I $`\lambda `$4471Å emission lines, respectively.

Fig. 1 also shows the maximum entropy value as a function of inclination obtained in entropy landscapes constructed assuming an optimum systemic velocity of –17 km s-1. It appears as though we cannot determine a unique inclination for AM Her, though our work suggests that the inclination may lie around 60–80. Given that we expect our results to be biased towards higher inclinations, the inclination is most likely to be around 60. Although this estimate is far greater than the $``$35 determined by Brainerd & Lamb (1985), it is consistent with the findings of Davey & Smith (1996) who used radial-velocity curve fitting to constrain the inclination to be between 45 and 60. Greeley et al. (1999) observed a partial eclipse of the He II $`\lambda `$1640 Å emission line region by the secondary star and hence inferred an inclination $``$45. In addition to this, Wickramasinghe et al. (1991) obtained $`i=52^{}`$ from polarimetric observations.

Fig. 3 shows the entropy landscape for AM Her constructed using $`i`$ = 60, $`\gamma `$ = –17 km s-1 and taking into account velocity smearing introduced by the 300 second exposure times. The brightest point on this plot indicates the map of highest entropy, and corresponds to $`M_1`$ = 0.36 M and $`M_2`$ = 0.18 M, which we adopt as our best estimate of the component masses. The optimum component masses decrease to $`M_1`$ = 0.24 M and $`M_2`$ = 0.12 M for $`i`$ = 80. These masses appear suspiciously low since AM Her lies above the period gap and we would not, therefore, expect the secondary mass to be significantly less than $`M_2`$ = $``$0.25 M. This partly confirms our speculation that our inclination determinations are biased towards high inclinations (Section 3.2) and that the inclination is not likely to be greater than 60.

The mass of the secondary star derived from the entropy landscape agrees well with the value of $`M_2`$ = 0.20 – 0.26 M estimated by Southwell et al. (1995). This also agrees with the main-sequence mass-period relation derived by Smith & Dhillon (1998), which gives a secondary star mass of 0.23$`\pm `$0.02 M.

Meaningful comparison of the primary mass derived in our work with other published white dwarf masses is, unfortunately, more difficult as they cover a wide range of values, including: 0.39 M (Young et al., 1981), 0.69 M (Wu, Chanmugam & Shaviv, 1995), 0.75 M (Mukai & Charles, 1987), 0.91 M (Mouchet, 1993) and 1.22 M (Cropper, Ramsay & Wu, 1998). Gänsicke et al. (1998), however, concluded that the white dwarf mass in AM Her was between 0.35 – 0.53 M and ruled out a white dwarf mass as high as 1.22 M. In addition to their secondary star mass estimate, Southwell et al. (1995) derived $`q`$ = 0.47$`\pm `$0.08, in agreement with the value of $`q`$ = 0.53 obtained in this work (Fig. 2).

### 4.2 The surface map of AM Her

The Roche tomogram for AM Her is presented in Fig. 4 and shows a decrease in the amount of Na I absorption around the inner Lagrangian ($`L_1`$) point on the trailing hemisphere. We tested the significance of this feature using the following procedure:

1. Fit the ‘true’ map using the observed data (Fig. 4).
2. Compute 200 ‘trial’ maps from simulated data sets constructed using bootstrap re-sampling. From our observed trailed spectrum containing $`n`$ data points we formed a simulated trailed spectrum by selecting, at random and with replacement, $`n`$ data values and placing these at their original positions in the new simulated trailed spectrum. For points that are not selected, the associated error bar was set to infinity and hence they were effectively omitted from the fit. For points selected more than once, the error bars were divided by the square root of the number of times they were picked (see Paper I for more details).
3. Calculate 200 ‘difference’ maps by subtracting each ‘trial’ map from the ‘true’ map. These ‘difference’ maps give the scatter in each pixel. (It would be incorrect to now calculate a summary error statistic like the standard deviation, as the distribution of pixel values is often found to be non-normal; see Paper I).
4. Subtract some mean or comparison level from the ‘true’ map. For instance, if one wished to determine the significance of a starspot, an appropriate comparison level to subtract from the ‘true’ map would be a region of immaculate photosphere. Subtract the corresponding level from each of the ‘difference’ or ‘scatter’ maps.
5. Rank each pixel in the ‘true’ map computed in step (iv) against the corresponding pixel in the ‘scatter’ maps. This will show the fraction of fluctuation values that the observed value in the ‘true’ map exceeds.

For our first test we took the average intensity of the outer (and presumably non-irradiated) hemisphere of the star as our comparison level, thereby testing the reality of the asymmetry between the outer and inner hemispheres. The left-hand panel in Fig. 5 shows the resultant significance map, calculated using the procedure listed above. The bright greyscale depicts regions where the ‘true’ map ranks above 100 per cent of the ‘scatter’ maps, and therefore represents regions where the intensity is significantly different from the outer hemisphere. Dark greyscales represent regions of the map that do not rank highly against the scatter map. Features in these regions are not, therefore, significantly different from the average intensity of the outer hemisphere.

The left-hand panel in Fig. 5 shows that the irradiation seen on the inner hemisphere in the Roche tomogram is real, and appears more significant on the trailing hemisphere (i.e., $`\varphi `$ = 0.25). Note that the significance of the regions near the south pole (SP) is an artefact. Due to the inclination of AM Her, latitudes $`>i`$ in the southern hemisphere are never seen and hence the reconstructions always assume the (same) default level, with very little scatter in the intensity values. As a result, regions near the south pole will appear significantly different from the comparison level unless, of course, the comparison level is very similar to the default level.

The right-hand panel in Fig. 5 shows another significance map, but taking the comparison level as a region around the $`L_1`$ point on the leading hemisphere ($`\varphi `$ = 0.25). This shows that the asymmetry between the leading and trailing hemispheres is real, and that the effects of irradiation are primarily found on the trailing hemisphere. The asymmetry in the irradiation pattern between the northern and southern hemisphere is again an artefact, since regions on the southern hemisphere are less visible and the data does not constrain the reconstructions in these regions. Note that in this case, the comparison level selected was similar to the default level of the map, and therefore regions near the south pole (SP) no longer appear to be significant.

Our results are in agreement with the surface maps of AM Her obtained using the method of radial-velocity curve fitting (Davey & Smith 1992, 1996). They proposed that the lack of irradiation on the leading hemisphere could be explained if the incoming accretion column was able to block the radiation produced close to the white dwarf and found that the accretion geometry given by Cropper (1988) supported this argument.

## 5 Results for QQ Vul

### 5.1 The binary geometry of QQ Vul

The systemic velocity and inclination of QQ Vul were determined in the manner outlined in Section 3.2, and the results are shown in Fig. 1. We find that the systemic velocity is well constrained with a value of –14 km s-1, comparable with the value of 17 $`\pm `$ 23 km s-1 obtained by Mukai & Charles (1987). Catalán et al. (1999) found $`\gamma `$ = –0.2 $`\pm `$ 2.2 km s-1 from the results of elliptical fits to the Na I absorption lines used in this work. In addition, using the same method they also obtained $`\gamma `$ = –8.3 $`\pm `$ 1.6 km s-1 using data obtained on the WHT in 1993. One explanation for this difference in the systemic velocity is if a change in the intensity distribution between the two epochs can alter the measured value of $`\gamma `$ (we explore this possibility in Section 7).

It is clear from Fig. 1 that it is difficult to select a single value for the inclination given the scatter in the distribution of the points. The quality of the reconstructions does, however, begin to deteriorate towards lower inclinations allowing us to set a lower of limit of $``$72 for the inclination of QQ Vul.

Previous estimates of the orbital inclination include $`i`$ = 50–70 (Schwope et al., 2000), $`i`$ = 65$`{}_{}{}^{}\pm `$7 (Catalán et al., 1999), $`i`$ = 60$`{}_{}{}^{}\pm `$ 14 (Mukai & Charles, 1987) and 46$`{}_{}{}^{}<i<`$ 74 (Nousek et al., 1984). In addition, Nousek et al. (1984) also found evidence of a brief eclipse of the cyclotron region combined with the simultaneous disappearance of emission lines. They attributed this to an eclipse by the secondary star, forcing the orbital inclination to the upper limit of 46$`{}_{}{}^{}<i<`$ 74. From EUVE observations Belle, Howell & Mills (2000) also suggested that the X-ray emitting region is eclipsed by the secondary star, constraining the inclination to $`>`$ 60. Certainly, in this work we can rule out inclinations $`<`$ 60 due to the low entropy of the reconstructions at these inclinations. These claims of eclipse features, however, are outside phase 0 and thus cannot be attributed to the secondary star. The lack of an X-ray eclipse at this phase means that an inclination $`>`$ 72 can be excluded.

Fig. 3 shows the entropy landscape for QQ Vul, constructed assuming an inclination of 72, from which we derive optimum masses of $`M_1`$ = 0.66 M and $`M_2`$ = 0.42 M. We adopt these masses and inclination for the reconstruction in Section 5.2. The derived masses, however, decrease smoothly as the inclination is increased to 90 and cover the range $`M_1`$ = 0.58–0.66 M, $`M_2`$ = 0.34–0.44 M at a constant $`q`$ = 0.63. This agrees well with, and is better constrained than, the masses of $`M_1`$ = 0.54$`{}_{0.16}{}^{}{}_{}{}^{+0.21}`$ M and $`M_2`$ = 0.30$`\pm `$0.10 M estimated by Catalán et al. (1999).

The white-dwarf mass we derive agrees with the values of $`M_1`$ = 0.58$`{}_{0.09}{}^{}{}_{}{}^{+0.44}`$ M (Mukai & Charles, 1987) and $`M_1`$ = 0.59 M (Mouchet, 1993), but is almost half that of Cropper et al. (1998) and Wu et al. (1995) who obtained $`M_1`$ = 1.22 M and $`M_1`$ = 1.1–1.3 M, respectively. The two latter mass determinations make use of the X-rays emitted from the accreting white dwarf and involve fitting model spectra to the observed X-ray spectra, and are therefore subject to underlying assumptions in the model.

The mass ratio of $`q`$ = 0.63 obtained in this work is in excellent agreement with previous determinations, including $`q`$ = 0.54$`\pm `$0.14 (Catalán et al., 1999) and $`q`$ = 0.45–0.67 determined by Schwope et al. (2000) from Doppler tomography and polarimetry.

### 5.2 The surface map of QQ Vul

The Roche tomogram of QQ Vul is presented in Fig. 6 and shows a distinct reduction in the Na I absorption around the $`L_1`$ point which can be explained by irradiation from the white dwarf and accretion regions ionising the Na I on the inner face of the secondary, as in AM Her. The effects of irradiation also appear to be slightly stronger on the trailing hemisphere, suggesting that the leading hemisphere may be partially shielded by the accretion stream or curtain.

Another notable feature is a bright patch on the leading hemisphere, corresponding to an increased Na I flux deficit around phase 0.75. A similar feature can also be seen in the tomogram of IP Peg (Fig. 8) and larger flux deficits at $`\varphi `$ = 0.75 than at $`\varphi `$ = 0.25 have also been reported for the Na I doublet in HT Cas (Catalán, 1999) and in the TiO light curves of Z Cha (Wade & Horne, 1988).

In their study of QQ Vul using the same data used in this work, Catalán et al. (1999) found that their surface maps derived from radial-velocity curve fitting did not give satisfactory fits to the velocities and line fluxes between orbital phases $`\varphi `$ = 0.6 and $`\varphi `$ = 0.8. The bright patch seen around phase 0.75 in Fig. 6 is almost certainly the cause of their unsatisfactory fits. Catalán et al. (1999) attributed this discrepancy to the presence of starspots. This interpretation is almost certainly wrong, however, as an increase in the Na I flux deficit is inconsistent with a star-spot, as the Na I flux deficit is known to decrease with later spectral type (Brett & Smith, 1993) and star-spots are generally cooler than the surrounding photosphere. Therefore, a starspot should appear as a dark feature in the tomograms. (Although spots hotter than the photosphere have also been imaged, e.g. Donati et al. 1992).

The reality of the spot feature can be assessed in Fig. 7, which shows the significance map when we take the comparison level as the average intensity on the outer hemisphere, avoiding both the region that appears to be irradiated and the bright patch. This shows that the bright patch is not significant. The significance map does, however, show that the irradiated inner hemisphere is real and there is evidence for stronger irradiation of the trailing hemisphere.

As with AM Her, the reduction in absorption around the $`L_1`$ point can be explained by irradiation from the white dwarf and accreting regions ionising the Na I on the inner face of the secondary, with possible shielding by the accretion stream/curtain reducing the effect on the leading hemisphere.

## 6 Results for IP Peg

### 6.1 The binary geometry of IP Peg

The optimum systemic velocity, as shown in Fig. 1, is –9 km s-1. This is significantly different from the values of 56 $`\pm `$ 4 km s-1 obtained by Martin, Jones & Smith (1987) and 31 $`\pm `$ 7 km s-1 obtained by Martin et al. (1989). Fig. 1 also shows that we have been unable to determine the inclination of IP Peg. This is due to the small variations in the line strength over the orbital cycle, which provides our only constraint on the inclination before taking into consideration the eclipse width. The reconstructions do, however, yield a consistent result of 0.43 for the mass ratio of IP Peg (Fig. 2). This is in disagreement with $`q`$ = 0.55–0.63 and $`q`$ = 0.55–0.62 found by Martin et al. (1989) and Marsh (1988), but agrees with 0.35$`<`$q$`<`$0.49 (Wood & Crawford, 1986), $`q`$ = 0.39$`\pm `$0.04 (Catalán, 1999) and $`q`$ = 0.32$`\pm `$0.08 (Beekman et al., 2000). All the authors, except Beekman et al. (2000), made corrections for the effects of irradiation where appropriate.

Although we cannot determine the inclination using the surface maps alone, IP Peg is an eclipsing system, and its geometry is further constrained by the observed eclipse width, which is dependent upon the inclination and size of the secondary star’s Roche lobe. Fig. 3 shows the entropy landscape constructed by varying the inclination at each different mass pairing in order to match the eclipse width of $`\mathrm{\Delta }\varphi `$ = 0.0863 (Wood & Crawford, 1986). This results in optimum masses of $`M_1`$ = 1.18 M and $`M_2`$ = 0.50 M at an inclination of 82.3. Altering the eclipse width to the other extreme of $`\mathrm{\Delta }\varphi `$ = 0.0918 given by Wood & Crawford (1986) changes this determination to $`M_1`$ = 1.16 M and $`M_2`$ = 0.50 M and the optimum inclination to 84.4. Our best estimates are therefore $`M_1`$ = 1.16 – 1.18 M, $`M_2`$ = 0.50 M and $`i`$ = 82–85.

Our secondary star mass determination appears to lie between the values of $`M_2`$ = 0.59–0.75 M (Martin et al., 1989), $`M_2`$ = 0.55–0.73 M (Marsh, 1988) and $`M_2`$ = 0.33$`{}_{0.05}{}^{}{}_{}{}^{+0.14}`$ M (Beekman et al., 2000), $`M_2`$ = 0.33$`\pm `$0.07 M (Catalán, 1999).

The primary mass is consistent with the values of $`M_1`$ = 1.05–1.25 M determined by Martin et al. (1989), $`M_1`$ = 0.99–1.19 M (Marsh, 1988) and $`M_1`$ = 1.05$`{}_{0.07}{}^{}{}_{}{}^{+0.14}`$ M (Beekman et al., 2000), but slightly larger than $`M_1`$ = 0.85$`\pm `$0.09 M determined by Catalán (1999).

### 6.2 The surface map of IP Peg

The Roche tomogram of IP Peg in the light of the Na I absorption doublet is shown in Fig. 8. The most noticeable feature is the weakened absorption flux deficit on the inner hemisphere of IP Peg, which can be attributed to irradiation. The effects of irradiation appear to be strongest, however, towards the leading hemisphere of the secondary star. We also see a bright patch of enhanced Na I absorption on the leading hemisphere similar to that seen in QQ Vul (Section 5.2).

The significance map constructed taking the comparison level to be the average intensity of the outer hemisphere, avoiding the bright patch, is shown in Fig. 9. From this we can see that the bright patch is not significant, whereas the irradiated region on the leading hemisphere is.

In their radial velocity curve-fitting studies, Davey & Smith (1992) also found that the effects of irradiation were strongest on the leading hemisphere, and that regions of IP Peg appeared to be irradiated some distance onto the back half of the star. They found that this feature could not be explained by irradiation from the bright spot, as the illuminating source would have to originate from a direction about 45 away from the line of centres and instead attributed this asymmetry to circulation currents (e.g. Martin & Davey 1995).

It is clear in Fig. 8 that we cannot identify any significant irradiation effects on the back half of the star and, as such, cannot provide any support for the evidence of circulation currents on IP Peg. It seems likely, instead, that the irradiation of the leading hemisphere can be explained by the bright spot, which is located on the correct side of the secondary star to cause this effect, and is a strong illuminating source (as seen in, for example, Szkody 1987).

Our Roche tomogram, however, bears little resemblance to the surface map that Davey & Smith (1992) computed from a model for a secondary star heated equally by the bright spot and the white dwarf (their Fig. 15), in that there does not appear to be any significant signs of irradiation on the trailing hemisphere. This suggests that irradiation from both the white dwarf and the bright spot is shielded from this region. An obvious way to block irradiation from the surface of the white dwarf is through shielding by the accretion disc. Irradiation from the bright spot may be blocked by the accretion stream, especially if much of the emission from the bright spot is located downstream from the impact point, as suggested in the eclipse maps of IP Peg by Bobinger et al. (1999). This would result in shielding of the trailing hemisphere, but still allow the bright spot to irradiate the leading hemisphere.

One further explanation for the features seen around the $`L_1`$ point is that they could be caused by the effects of an eclipse of the secondary star by the accretion disc. This would cause darker regions to appear in the Roche tomogram around the $`L_1`$ point and a reduction in the quality of the fit. Reconstructions carried out neglecting phases between 0.45–0.55 show no improvement in the quality of the fit, and there is negligible difference in the Roche tomograms. We are therefore confident that these features are not the result of an eclipse.

## 7 Results for HU Aqr

### 7.1 The binary geometry of HU Aqr

Fig. 1 shows the systemic velocity and inclination of HU Aqr as a function of entropy. We obtain a value of $`\gamma `$ = 7 km s-1 and, although it is difficult to assign a value to the inclination, the reconstruction quality begins to decrease rapidly for inclinations less than 80, and we therefore adopt a best estimate of $`i`$ 80.

Our estimate of the systemic velocity is compatible at the 2$`\sigma `$ level with $`\gamma `$ = 83 $`\pm `$ 44 km s-1 derived by Glenn et al. (1994) who fitted a sinusoid to the narrow H$`\alpha `$ component thought to originate from the secondary star. Schwope et al. (1997) derived a systemic velocity of –7.8 $`\pm `$ 1.4 km s-1 by fitting radial velocity curves to the same data used in this work. Although this result disagrees with the value obtained using Roche tomography, the sine-fitting method takes no account of any non-uniform intensity distribution present on the secondary star, whereas Roche tomography does. Since there is clearly an asymmetrical intensity distribution on the secondary star in HU Aqr (see Section 7.2) we believe that this is the cause of the difference between the two results.

In order to test the reliability of the sine-fitting method in determining the systemic velocity, we generated a model secondary star with the same binary geometry as we have found for HU Aqr in this work ($`M_1`$ = 0.61 M, $`M_2`$ = 0.15 M and $`i`$ = 84.4, see later) but with $`\gamma `$ = 0 km s-1. The inner hemisphere of the model was then uniformly irradiated except for a totally dark shadow covering 40 per cent of the leading hemisphere. The outer hemisphere of the secondary star was also set to be totally dark. From this model, trailed spectra of infinite signal-to-noise were generated, assuming 100 equally spaced phase bins and 30 second exposure lengths.

By fitting radial-velocity curves to this fake data, we obtained a systemic velocity of –14 km s-1. Fig. 10 shows the result of the radial-velocity curve fit (solid line) to the synthetic trailed spectrum, and the true motion of the centre-of-mass of the secondary star (dashed line) for comparison. The value of $`\gamma `$ = –14 km s-1 is significantly different from the true value of $`\gamma `$ = 0 km s-1 used in the model, and brings our systemic velocity of 7 km s-1 into agreement with the value of –7.8 $`\pm `$ 1.4 km s-1 obtained by Schwope et al. (1997). This shows that a knowledge of the intensity distribution across the secondary is required to accurately determine the systemic velocity (as also found by Schwope et al. 1997). This may have implications for the accurate determinations of $`\gamma `$’s used in CV age tests (e.g., van Paradijs, Augusteijn & Stehle 1996 and North et al. 2002).

Since HU Aqr is an eclipsing system, its inclination has been reasonably constrained by previous studies, including estimates of $`i`$ = 80$`{}_{}{}^{}\pm `$5 (Glenn et al., 1994) obtained from fits to polarisation curves, and $`i`$ = 85.6 (Schwope et al., 2001) from modelling the soft X-ray light-curves. Our determination of $`i`$ 80 does not provide any further constraint to those values already published.

Previous estimates of the mass ratio, however, are somewhat more uncertain and range from roughly 0.18 – 0.4 (e.g. Schwope et al. 1997, Hakala et al. 1993 and Heerlein, Horne & Schwope 1999). Fig. 2 shows the component masses and mass ratio obtained when reconstructions are carried out over different inclinations. This yields a mass ratio for HU Aqr of $`q`$ = 0.25, which is in excellent agreement with $`q`$ = 0.256 derived by Heerlein et al. (1999) assuming $`i`$ = 84, and the value of $`q`$ = 0.25 that Schwope et al. (2001) found.

Fig. 3 shows the entropy landscape for HU Aqr which has been constructed by changing the inclination at each different mass pairing in order to match the eclipse width of $`\mathrm{\Delta }\varphi `$ = 0.0763 (Heerlein et al., 1999), as described for IP Peg in Section 6.1. This results in a best estimate of the component masses of $`M_1`$ = 0.61 M and $`M_2`$ = 0.15 M. This agrees with the result of Schwope et al. (2001), who determined $`M_2`$ = 0.15 M using the mass-radius relation of Caillault & Patterson (1990), and $`M_2`$ = 0.17 M according to the relation by Neece (1984). Our derived masses, combined with the eclipse width, give our best estimate of the inclination as 84.4.

### 7.2 The surface map of HU Aqr

The Roche tomogram of HU Aqr in the light of the He II 4686Å emission line is presented in Fig. 11. The two most prominent features in the Roche tomogram are the strong asymmetry between the inner and outer hemispheres of the secondary star and the weaker asymmetry between the leading and trailing hemispheres. The emission from the inner hemisphere is due to irradiation from the primary and the accretion regions, and the tomogram is compatible with 90 per cent of the emission arising from the inner hemisphere alone. (Note that the remainder of the emission that arises on the outer hemisphere is most likely due to the smearing effects of the default map).

The reality of the asymmetry between the leading and trailing hemispheres can be assessed from the left-most panel in Fig. 12, which shows a slice passing through the leading hemisphere (LH), north pole (NP), trailing hemisphere (TH) and south pole (SP) near the $`L_1`$ point of the secondary star. The triangular points represent the ‘true’ map values; 67 per cent of the bootstrapped ‘trial’ maps (measured relative to the mode of the distribution) lie within the solid line, and 100 per cent within the dotted line111We have not generated a significance map for HU Aqr since the low scatter in the pixel values derived from the bootstrap reconstructions mean that each region of the map would be significant when compared with any other region.. It is clear from this that the asymmetry between the leading and trailing hemispheres is not due to statistical errors. This asymmetry was also found by Schwope et al. (1997) who attributed it to shielding by the accretion curtain.

Closer examination of the tomogram reveals more subtle intensity variations across the secondary star. In particular there appears to be a region of more intense emission located towards the north pole but shifted slightly towards the leading hemisphere. The plot on the right hand-side of Fig. 12 shows a vertical slice cutting through the secondary star close to the polar regions but slightly towards the $`L_1`$ point. This again shows the shadow on the leading hemisphere (LH) as discussed earlier, but another feature of note is the drop in emission from the trailing hemisphere (TH) compared to the polar regions. This feature is not apparent near the $`L_1`$ point and only becomes significant for regions halfway between the centre-of-mass of the secondary and the $`L_1`$ point, disappearing again when the effects of irradiation are diminished near the terminator. This cannot be explained by shielding of the secondary star by the accretion curtain, since this would only block irradiation on the leading hemisphere.

HU Aqr does, however, have a strong stream component with as much as 40 per cent of the flux in the U band due to accretion stream emission (Hakala et al., 1993). Doppler maps by Schwope et al. (1997) also show that the accretion stream is a dominant emission-line feature. Presumably the accretion stream can also irradiate the secondary star and, since the stream is deflected in the direction of the orbital motion, the stream will be invisible to regions of the secondary star on the trailing hemisphere away from the $`L_1`$ point. This obscuration of the accretion stream by the secondary star itself may explain the observed shadow on the trailing hemisphere.

This feature could also be explained if systematic errors were able to produce a greater line flux in the centres of the lines, corresponding to artificially brighter polar regions in the Roche tomograms and giving the illusion of a shadow at lower latitudes. It was first thought that the assumption of a Gaussian profile for the emission from the secondary star could introduce such an artefact, since the actual profile from the secondary star is expected to be flatter, especially around quadrature. We found, however, that at low spectral resolutions the emission profile from the secondary is well approximated by a Gaussian. Artefacts could be generated if the emission from the various components of HU Aqr (e.g. from the secondary star and the accretion stream) are incorrectly disentangled, which is especially difficult at phases where the individual components are crossing each other. This situation is more difficult to assess and emphasises the importance of mapping lines that are known to originate solely from the secondary star.

## 8 Conclusions

The Roche tomograms show an asymmetric irradiation pattern on all 4 CVs. This can be explained if the accretion column/stream shields the leading hemisphere in AM Her, QQ Vul and HU Aqr. There is also evidence for a fainter shadow present on the trailing hemisphere of HU Aqr. We propose that the stream itself contributes significantly to the irradiation, and that the shadow on the trailing hemisphere occurs as these regions on the secondary star have no direct view of the majority of the stream. Further observations using absorption lines originating from the secondary star are needed to confirm this.

A different scenario occurs in IP Peg, where irradiation affects the leading hemisphere only. We suggest that the bright spot, which is on the correct side of the secondary star to cause such an irradiation pattern, may be the principal illuminating source in this system. Although we still expect irradiation of the trailing hemisphere by the white dwarf (as seen in the magnetic CVs), there is no evidence for this in the tomogram. It may be the case that the accretion disc in IP Peg is shielding the secondary star from irradiation by the white dwarf. Unlike Davey & Smith (1992), we find no evidence for irradiation extending onto the back half of the star, and hence cannot support their idea that circulation currents operate in IP Peg.

We have used the entropy landscape technique to derive the system parameters. Although we cannot reliably constrain the inclination, and therefore the masses, in the non-eclipsing CVs, the technique does provide a consistent result for the mass ratio that is independent of the inclination. For the two eclipsing CVs (IP Peg and HU Aqr), we believe that the masses derived in this work are the most accurate to date. In addition, it is found that measurements of the systemic velocity using circular orbit fits to the radial-velocity variations are unreliable, and are dependent upon the intensity distribution across the secondary star. We might therefore expect such measurements to be both line and time-dependent. Since the intensity distribution across the secondary star is known in Roche tomography, systemic velocity measurements using this technique should prove far more accurate.

## Acknowledgements

We thank Stuart Littlefair for many useful conversations. CAW is supported by a PPARC studentship. This project was supported in part by the Bundesministerium für Bildung und Forschung through the Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR) under grant number 50 OR 9706 8. We would also like to thank the referee, Andrew Cameron, for his comments, which substantially improved the clarity of the paper.