text:
# Overconstrained dynamics in galaxy redshift surveys

## 1 Introduction

Galaxy redshift surveys are elements of prime importance in the determination of cosmological parameters. Looking at redshift-space distortions we are able to obtain maximum likelihood estimates of $`\mathrm{\Omega }_\mathrm{m}`$ (Fisher, Scharf & Lahav 1994; Heavens & Taylor 1995; Baker et al. 1998; Hamilton 1998 and references therein; Dekel 1999a,b; Hamilton, Tegmark & Padmanabhan 2000) and the galaxy-galaxy power spectrum (Fisher et al. 1993; Feldman, Kaiser & Peacock 1994; Lin et al. 1996; Sutherland et al. 1999; Hamilton & Tegmark 2000). The dynamics of cosmic velocity flows is also another way to measure $`\mathrm{\Omega }_\mathrm{m}`$ (Willick et al. 1997; Dekel, Burstein & White 1997; da Costa et al. 1998; Kashlinsky 1998; Sigad et al. 1998; Susperregi 2001, hereafter S01) and the power spectrum (Strauss & Willick 1995 and references therein; Kolatt & Dekel 1997; Willick et al. 1997; Freudling et al. 1999; Zehavi & Dekel 1999). In the larger picture, galaxy redshift surveys contribute modestly, alongside a myriad of other datasets (CMB, SN1a data, weak lensing surveys, etc) to the undertaking of estimating the key cosmological parameters ($`\mathrm{\Omega }_\mathrm{m}`$, $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_B`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ or $`\mathrm{\Omega }_Q`$, $`h`$, $`n_s`$, $`n_t`$, etc)(for a review of how this inventory of over a dozen parameters is currently tackled, see e.g. Turner 1999; Primack 2000); all the datasets, including galaxy surveys, can be ideally combined by exploiting cosmic complementarity (Eisenstein, Hu & Tegmark 1998,1999) leading to demonstrably concordant, or at least not inconsistent, predictions.

A number of techniques have been employed to maximize the parameter information extracted from surveys. The Fisher information matrix approach is one example of this (Tegmark 1997; Tegmark, Taylor & Heavens 1997; Goldberg & Strauss 1998; Taylor & Watts 2000). This approach is sound and will lead to very accurate results when applied to forthcoming surveys such as SDSS and 2dF; the computational effort can be successfully minimised in those cases via data-compression techniques, e.g. via the Karnhunen-Loève “signal-to-noise” eigenmodes (Vogeley & Szalay 1996; Tegmark et al. 1998; Matsubara, Szalay & Landy 2000), quadratic compression into high resolution powers (Tegmark 1997; Tegmark, Taylor & Heavens 1998; Tegmark et al. 1998; Padmanabhan, Tegmark & Hamilton 2000), etc. There is, on the other hand, a case to be made for combining different galaxy redshift surveys and extracting information jointly from them, as an ensemble, prior to comparing them with different kinds of datasets. As a step prior to exploiting cosmic complementarity, it can be argued that redshift surveys are best taken into account as an ensemble than individually, however larger and better sampled is any given dataset in comparison to its predecessors. It is to be expected that a comparison of surveys leaves considerable room to manoeuvre to best rid of, or, at least compensate for, the errors inherent in each dataset, stochasticity, unaccounted for contaminants, etc.

The goal of this paper is to set out a procedure to extract an accurate estimate of $`\mathrm{\Omega }_\mathrm{m}`$ from galaxy redshift surveys, when the datasets are studied as an ensemble. The measurement of $`\mathrm{\Omega }_\mathrm{m}`$ will then be the optimal estimate for all the surveys considered. The method will also enable us to estimate the matter and galaxy-galaxy power spectra. The backbone of the procedure is essentially the least-action principle (LAP) reconstruction employed in S01. That paper showed that the LAP method is efficient to measure $`\mathrm{\Omega }_\mathrm{m}`$ from a given galaxy redshift survey and that it conveniently breaks the degeneracy between $`\mathrm{\Omega }_\mathrm{m}`$ and the bias, thus permitting to measure $`\mathrm{\Omega }_\mathrm{m}`$ within fairly arbitrary bias schemes. In this paper we consider a bias relationship that is scale-dependent in order to investigate the dynamics more realistically. Bias is in any case stochastic as well as scale-dependent (Pen 1998; Tegmark & Peebles 1998; Dekel & Lahav 1999; Tegmark & Bromley 1999; Taruya 2000), and the bias relationship is fully determined by the distribution $`P(\delta _\mathrm{g}|\delta _\mathrm{m})`$ (Dekel & Lahav 1999; Sigad, Branchini & Dekel 2000). The stochastic ingredient is largely sample-dependent, and the local scatter that is the measure of stochasticity in $`P(\delta _\mathrm{g}|\delta _\mathrm{m})`$ is a characteristic of each survey, its shot noise and environmental effects. The whole point of studying an ensemble of surveys is to minimise this effect, and therefore for simplicity we neglect the existence of stochasticity in the bias and we focus on its scale dependence.

The galaxy redshift surveys that we examine in this paper are PSC$`z`$, ORS, Mark III and SFI. The LAP method is initially applied to each survey individually in the manner of S01 with a scale-dependent bias model. Hence one obtains the cosmic fields and a measurement of $`\mathrm{\Omega }_\mathrm{m}`$ for each survey. Each survey is regarded as a constraint on the cosmic fields, and as they are all indeed different, the system that results from considering all constraints simultaneously is unavoidably overconstrained. Strictly speaking it is also incompatible. The method to find an optimal solution that is a meaningful representation of the ensemble entails relaxing all the constraints. The errors that this relaxation brings in are dealt with the least-square method, and this yields the best fit. The procedure gives us also an estimate of $`\mathrm{\Omega }_\mathrm{m}`$, which is really a joint measurement for all the surveys in the ensemble. From the cosmic fields and the bias relationship it is then straightforward to compute the power spectra, higher-order correlation functions, etc.

The article is structured as follows. A brief description of the LAP method is given in §2 and the algebra is relegated to the appendices, though the necessary steps to undertake the numerical implementation are set out in §2.1; the characteristics of the surveys that are relevant to this paper are summarized in §3; the procedure for finding LAP solutions for the ensemble of datasets is described in §4 for the most general case; in §5 we discuss the application of the method on the four surveys in question and their results; finally §6 discusses the main results of this paper, in comparison with similar measurements in the literature.

## 2 LAP preliminaries

The three ingredients of the LAP method are: (a) two boundary conditions and (b) the Lagrangian of the self-gravitating matter field. The two boundary conditions largely determine the cosmological model, as they are the endpoints in the evolution of our Universe. Initial homogeneity, as measured by the Sachs-Wolfe constraint on the CMB, suggests the first is

$$\delta (t0,𝐱)0,$$
(1)
where $`\delta `$ is the matter density constrast. The statistics of the primordial $`\delta `$ is thus unconstrained. The second boundary condition is the Universe as is represented in the current epoch by a galaxy redshift survey; this is characterised by the galaxy number-count $`n(𝐬)`$ ($`𝐬`$ denotes redshift coordinates, ($`cz`$,$`\theta `$,$`\phi `$)). This constraint is given by (S01)

$$\rho _s(𝐬)=x^2\frac{N_{\mathrm{gals}}}{V}\left[\frac{1+g_0(𝐱)}{1+\alpha _0(𝐱)^{\prime \prime }}\right],$$
(2)
where $`\rho _s`$ relates to the galaxy number-count via $`\rho _s\mathrm{d}n(𝐬)/\mathrm{d}s\mathrm{d}\Omega `$ (where $`\mathrm{d}\Omega `$ denotes a unit of solid angle), $`V`$ is the volume of the survey, $`x`$ is the radial comoving distance, $`\alpha `$ is the velocity potential (i.e. $`𝐯\alpha `$), the prime denotes the line-of-sight derivative $`\mathrm{d}/\mathrm{d}x`$, $`g\delta _{\mathrm{gals}}`$ is the galaxy number density contrast and the subindex $`0`$ denotes the present time. We use comoving quantities consistently throughout the paper. The bias model that we shall adopt is scale-dependent, such that

$$g_0(𝐤)=b(k)\delta _0(𝐤).$$
(3)
The method is used following the same procedure as that expounded in S01, with the addition of a scale-dependent bias. The revelant LAP equations are summarized in Appendix A for reference.

### 2.1 Numerical resolution

We apply the algorithm as follows.

* The dataset $`𝒟`$ is transformed into a $`z`$-space field $`\rho _s(𝐬)`$, by computing the discrete derivatives of the galaxy number-count $`n(𝐬)`$, binned to the smallest redshift spacing $`\mathrm{\Delta }z`$, and then the result is smoothed choosing an appropriate smoothing length $`r_s`$. We shall only implement Gaussian smoothing, $`W(k)=\mathrm{exp}(k^2r_s^2/2)`$.
* The goal is to find a best fit for the modes $`\delta _y(t)`$, $`\alpha _y(t)`$ (where $`yrlm`$ in the spherical harmonics and Bessel functions decomposition, following the formalism of Appendix A). A consistent method to achieve this is to get progressively closer to the real solution by trying Ansatze of increasing non-linearity in successive iterations. The starting point is the linear solution, which is approximated by inverting the relation $`\delta _s^2\alpha _0`$, where

$$\delta _s(s_{\mathrm{max}}/4\pi N_{\mathrm{gals}})\rho _s1.$$
(4)
For the linear $`\delta `$ field, we estimate

$$\delta _0\delta _s(𝐱+\widehat{𝐱}\alpha _0^{}).$$
(5)
In the surveys where the radial velocities are part of the observational input, in our case Mark III and SFI, the $`v_r`$ data are used to construct the first Ansatz for $`\alpha _0^{}`$, which is then used to compute (5).
* The linear fields $`\delta _{y,0}`$,$`\alpha _{y,0}`$ are the first Ansatz that we input in the LAP system (38). We transform these first into the coefficients $`\delta _y^{(n)}`$, $`\alpha _y^{(n)}`$, using (31), (32). These are the coefficients of the polynomials of order $`N`$ that approximate the fields. This initial Ansatz is equivalent to linear theory and it permits us to solve the homogeneous system. It requires an initial guess of the value of $`\mathrm{\Omega }_\mathrm{m}`$, as the matrix coefficients $`𝒮^\delta `$ depend on it. It is worthwhile to experiment with a ample range of values such as $`0.1\begin{array}{c}<\hfill \\ \hfill \end{array}\mathrm{\Omega }_\mathrm{m}\begin{array}{c}<\hfill \\ \hfill \end{array}1.0`$, to examine the impact of these variations in the convergence of the solutions. The solution obtained from the homogeneous system is least-square fitted to (1),(2). This in turn requires an bias model, and from (3) using Parseval’s theorem (i.e. $`fg=fg`$, where $`f`$ denotes the Fourier transform of $`f`$), we have that $`g(𝐱)=b(x)\delta (𝐱)`$. In the spherical harmonics and Bessel function decomposition this means

$$g(t,𝐱)=\underset{y}{}b_{rl}\delta _y(t)j_l(k_rx)Y_{lm}(\theta ,\phi ).$$
(6)
This assumes that the effect of bias reflects on a superposition of radial modes that are rotationally invariant. In this paper we have fitted $`b_{rl}`$ for $`r5`$ and $`l6`$. Our starting point Ansatz in (2) is $`b_{rl}=1.0`$.
* The resulting $`\delta _y`$,$`\alpha _y`$ are then used to construct the RHS of (38), the inhomogeneous system. These quadratic terms represent the coupling among the normal modes of the perturbations, and their amplitude is a measure of the degree of non-linearity present.
* Successive iterations thus solve (38) and we combine this operation with the least-square fitting the solution to the constraints (1),(2). The procedure eventually yields the correct $`\delta _y`$, $`\alpha _y`$ if the errors in satisfying (1), (2) are monotonously decreasing after each iteration. We use variations around the solution to construct the Fisher information matrix

$$_{ij}\frac{^2(\delta ,\alpha )}{\theta _i\theta _j}=d\rho _s\frac{^2}{\theta _i\theta _j}\mathrm{ln}\rho _s(\delta ,\alpha ),$$
(7)
where $`\theta _0\mathrm{\Omega }_\mathrm{m}`$ and $`\theta _ib_{rl}`$ for $`i>0`$, conveniently ordering the $`b_{rl}`$ in a sequence. The likelihood function $`L`$, is defined as $`L\mathrm{exp}()`$. (7) enables us to compute the relative likelihood of estimates of the parameters, particularly $`\mathrm{\Omega }_\mathrm{m}`$. The working assumption is that $`\delta _s`$ in is a random field, resulting from the galaxy number-count, which is modelled as a 3D Poisson distribution. This enables us to compute $``$. The strength of the LAP method is that it does not require a prior, such as the CDM power spectrum, to estimate the likelihood of the parameters.

## 3 Datasets

### 3.1 PSC$`z`$

The sample contains 14,819 galaxies, within a spherical region of radius $`x_{\mathrm{max}}\begin{array}{c}>\hfill \\ \hfill \end{array}45,000`$ $`\mathrm{km}\mathrm{s}^1`$(its median is at $`8,500`$ $`\mathrm{km}\mathrm{s}^1`$and the galaxy number-count tails off sparsely up to $`\begin{array}{c}>\hfill \\ \hfill \end{array}30,000`$ $`\mathrm{km}\mathrm{s}^1`$), extracted form the IRAS Point Source Catalogue (Saunders et al. 2000), following the selection procedure given by conditions log$`(f60/f25)>0.5`$ and log$`(f100/f60)<0.75`$ (where $`f\alpha `$ denotes the measured flux at $`\alpha `$ microns). The former discards all stars whilst the latter excludes galactic cirruses located away from the Zone of Avoidance (ZoA). The result retains most galaxies plus about 300 mixed contaminants, which are excluded on a one-to-one basis. The source catalogue contains 17,060 sources, of which 1,593 are discarded as they are either within the Milky Way or repeated entries; a further 648 are rejected as very faint galaxies and unidentifiable sources. The resulting 6,500 are a combination of IRAS 1.2 Jy and QDOT galaxies plus 3,000 additional sources. The sky coverage is $`84\%`$ (the ZoA spans a non-axially symmetric region of $`|b|\begin{array}{c}<\hfill \\ \hfill \end{array}5^{}`$). We adopt the unmasked survey, therefore maximising the sky coverage and a selected spherical subvolume of radius $`x_{\mathrm{max}}15,000`$ $`\mathrm{km}\mathrm{s}^1`$. General descriptive information and access to data are available at the PSC$`z`$ webpage http://www-astro.physics.ox.ac.uk/~wjs/pscz.html.

### 3.2 ORS

The sample contains $`8,457`$ optically selected galaxies, of an angular size $`\begin{array}{c}<\hfill \\ \hfill \end{array}1.9`$ arcmin, located within an approximately spherical region $`18,000`$ $`\mathrm{km}\mathrm{s}^1`$(Santiago et al. 1995; Hudson et al. 1995; Baker et al. 1998). The dataset is densely sampled within $`x_{\mathrm{max}}8,000`$ $`\mathrm{km}\mathrm{s}^1`$, which is the spherical subvolume that we shall consider, and it is sparse in the region $`x\begin{array}{c}>\hfill \\ \hfill \end{array}8,000`$ $`\mathrm{km}\mathrm{s}^1`$. The sample is a compilation of optically-selected galaxies extracted from the Upsala General Catalogue (UGC; $`\delta 2.5^{}`$), European Southern Observatory (ESO; $`\delta <17.5^{}`$) catalogue and Extension to the Southern Observatory Catalogue (ESGC; covering the remaining region just south of the Celestial Equator. The sky coverage of these samples is respectively $`34\%`$,$`19\%`$ and $`9\%`$. ORS shares a number of galaxies in common with PSC$`z`$ though it contains a larger fraction of spirals. The predicted velocity fields in ORS and PSC$`z`$ are also in good agreement (Baker et al. 1998)). The sky coverage is 62%, i.e. significantly more limited than in PSC$`z`$ (and even IRAS 1.2 Jy), due to the extinction in the ZoA, and it spans over $`|b|\begin{array}{c}>\hfill \\ \hfill \end{array}20^{}`$. In the sample utilized ORS is filled with 823 IRAS galaxies in order to achieve $`|b|\begin{array}{c}>\hfill \\ \hfill \end{array}5^{}`$ coverage. A drawback of ORS is its poor uniformity, as extinction affects optically-selected galaxies more than it does IRAS galaxies. Further descriptive information and access to data are available at the URL: http://www.astro.princeton.edu/~strauss/ors/.

### 3.3 Mark III

The sample contains radial velocities of $`\begin{array}{c}>\hfill \\ \hfill \end{array}3,400`$ galaxies (Willick et al. 1995,1996,1997), compiled from several sets of S0 spirals and ellipticals. The sky coverage is the entire sky with the exception of the ZoA, at $`|b|\begin{array}{c}>\hfill \\ \hfill \end{array}20^{}30^{}`$. The data are located within a spherical region of radius $`x_{\mathrm{max}}6,000`$ $`\mathrm{km}\mathrm{s}^1`$that is well sampled, though it is fairly anisotropic, as it spans to $`x_{\mathrm{max}}8,000`$ $`\mathrm{km}\mathrm{s}^1`$ in some directions and is consigned to $`x_{\mathrm{min}}4,000`$ $`\mathrm{km}\mathrm{s}^1`$ in others. We consider the spherical subvolume of radius $`x_{\mathrm{max}}8,000`$ $`\mathrm{km}\mathrm{s}^1`$. Distances are inferred via Tully-Fisher and $`D_n\sigma `$ distance indicators and they entail an error $`1721\%`$. The data must be carefully prepared to correct for Malmquist biases (following the procedure set out in Sigad et al. (1998) and used by Zaroubi, Hoffman & Dekel (1999); also in a improved version by Dekel et al. (1999)). This results in the correction of the distances of 1,200 objects. Further descriptive information and access to data are available at the URL: http://redshift.stanford.edu/MarkIII/.

### 3.4 SFI

The sample contains radial velocities of $`1,300`$ Sbc-Sc galaxies, with inclination $`\begin{array}{c}>\hfill \\ \hfill \end{array}45^{}`$ north of $`\delta <45^{}`$ and galactic latitude $`|b|\begin{array}{c}>\hfill \\ \hfill \end{array}10^{}`$ (da Costa et al. 1996,1998; Haynes et al. 1999a,b). The data are distributed within an approximately spherical region of radius $`x_{\mathrm{max}}7,000`$ $`\mathrm{km}\mathrm{s}^1`$. Individual distances are computed in a similar manner as in Mark III, and the distance errors are $`1520\%`$. A comparison with Mark III shows a good agreement (da Costa et al. 1996; ditto for the POTENT reconstruction analysis of Dekel et al. 1999) and it is also in agreement with the predicted velocity field of IRAS 1.2 Jy (da Costa et al. 1998).

## 4 Method

The LAP method applied to one survey is completely determined by the following two constraints:

$$\chi _0(\delta )=0,$$
(8)
$$\chi _{\mathrm{survey}}(\delta ,\alpha )=0,$$
(9)
where (8) is the condition of homogeneity (1), and (9) contains the data given by the survey (2). Any additional dataset adds one constraint of the type (9). Let us suppose we have $`N`$ such datasets. Obviously in practice all samples contain systematic and random errors and not all $`N+1`$ contraints are satisfied exactly by the same solution; therefore, the system $`\chi _i=0`$ with $`i=0,1,\mathrm{},N`$ need not have a consistent solution $`\delta `$,$`\alpha `$. To find a solution, the overconstrained system must be relaxed so that the solution satisfies the constraints (9) within the least margin of error, $`\chi _iϵ_i`$, $`(i>0)`$, and we optimize this by computing the stationary quantity

$$\stackrel{~}{\delta }\left[\underset{i=1}{\overset{N}{}}w_i\chi _i^2(\delta ,\alpha )\right]=0,$$
(10)
where $`\stackrel{~}{\delta }`$ denotes a variation, not a density contrast, and $`w_i`$ is the weight given to each survey in the sum. From the central limit theorem, we assume that $`ϵ_i`$ are normally distributed about zero, i.e. $`ϵ_i=0`$ and $`E(ϵ_i^2)=\sigma _i`$. We will therefore adopt the maximum likelihood values $`w_i=\sigma _i^2`$, where $`\sigma _i^2`$ is the variance of the distribution $`P[\chi _i]`$. The likelihood function is then

$$L=(2\pi )^{N/2}\left(\underset{i=1}{\overset{N}{}}\sigma _i^1\right)\mathrm{exp}\left(\frac{1}{2}\underset{i=1}{\overset{N}{}}\frac{\chi _i^2}{\sigma _i^2}\right).$$
(11)
The solution that satisfies (10) and the LAP equations (38) is the best fit to the $`N`$ surveys that one can find.

### 4.1 Domain of application

The first consideration is to determine the radius of the spherical volume to which the LAP solutions will be consigned. Ideally, one would wish to use the largest selected subvolume possible, in this case that of PSC$`z`$, $`x_{\mathrm{max}}15,000`$ $`\mathrm{km}\mathrm{s}^1`$. This would nonetheless give too little weight to the other datasets, as the constraints $`\chi _i0`$ given by the smaller surveys would have little impact in the outcome of $`\delta ,\alpha `$ (the fractions of the selected volumes of ORS, Mark III and SFI, with respect to PSC$`z`$ are, respectively, 0.15, 0.06 and 0.10). On the other hand, $`x_{\mathrm{max}}`$ should not be too small either, as a large part of data from the larger samples may thus be left unused. An criterion to determine $`x_{\mathrm{max}}^{\mathrm{LAP}}`$ is the following. We assume that the weight each dataset ought to have in determining the volume of the LAP region is given by the density of its sampling, given by the ratio

$$\eta _i\frac{N_{\mathrm{gals}}^i}{V_i}\left(\underset{j=1}{\overset{N}{}}\frac{N_{\mathrm{gals}}^j}{V_j}\right)^1.$$
(12)
This measure of sampling is given in Table 1 for the surveys of our study. It shows that ORS and Mark III are the most densely sampled datasets and SFI is the sparsest. PSC$`z`$ is also on the whole fairly sparsely sampled. The volume of the LAP region is hence the weighed sum

$$V_{\mathrm{LAP}}=\underset{i=1}{\overset{N}{}}\eta _iV_i.$$
(13)
Using the values given in Table 1, for the case at hand this results in $`8,860`$ $`\mathrm{km}\mathrm{s}^1`$. In fact this seems to be a good figure also for PSC$`z`$ as the median redshift of that catalogue is at $`8,500`$ $`\mathrm{km}\mathrm{s}^1`$. We shall adopt the round figure $`x_{\mathrm{max}}^{\mathrm{LAP}}=9,000`$ $`\mathrm{km}\mathrm{s}^1`$.

### 4.2 Solutions

We proceed as follows. First, we apply the LAP method on each sample separately on their own domain of application, following the numerical procedure described in §2.1. Hence, for each survey we obtain $`\delta _i`$,$`\alpha _i`$ within a spherical volume of radius $`x_{\mathrm{max}}^i`$ (using the values given in Table 1). This also yields an estimate of $`\mathrm{\Omega }_\mathrm{m}`$ and an evaluation of the bias coefficients $`b_{rl}`$. The results of the measurement of $`\mathrm{\Omega }_\mathrm{m}`$ in each survey are given in the next section.

Let us consider these solutions, $`\delta _i`$,$`\alpha _i`$, to construct an initial Ansatz of the overall solution that satisfies (10). This Ansatz is given by

$$\delta =\underset{i=1}{\overset{N}{}}\eta _i\delta _i,$$
(14)
$$\alpha =\underset{i=1}{\overset{N}{}}\eta _i\alpha _i.$$
(15)
Thus we adopt the sampling factor as the best weight to average the fields. A few trial and error experiments show this is a good starting point. In order to consign this Ansatz to the domain (13), if $`x_{\mathrm{max}}^{\mathrm{LAP}}<x_{\mathrm{max}}^i`$ then one discards from the survey $`i`$ all the data that lie within the shell $`x_{\mathrm{max}}^{\mathrm{LAP}}<x<x_{\mathrm{max}}^i`$; and if $`x_{\mathrm{max}}^{\mathrm{LAP}}>x_{\mathrm{max}}^i`$ then ones makes the fields identically equal to zero within the range $`x_{\mathrm{max}}^i<x<x_{\mathrm{max}}^{\mathrm{LAP}}`$, effectively “padding” that shell to match the LAP volume. In the example of this paper, PSC$`z`$ data are discarded in the range $`9,000`$ $`\mathrm{km}\mathrm{s}^1`$$`<x<`$15,000 $`\mathrm{km}\mathrm{s}^1`$, and ORS, Mark III and SFI are padded with zeros within the shell $`x_{\mathrm{max}}^i<x<9,000`$ $`\mathrm{km}\mathrm{s}^1`$.

The Ansatz (14),(15) is our starting point to compute the LAP solution over the domain (13). The numerical problem is then to solve (38), introducing variations in the fields $`\delta `$, $`\alpha `$ and in the free parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`b_{rl}`$ to eventually pin down the stationary solution in (10). The numerical resolution is carried out following the procedure set out in §2.1, where all references to the constraints (1), (2) are substituted by (10).

## 5 Results

We apply the method described in §4 to the four datasets in question, PSC$`z`$, ORS, Mark III and SFI. For each survey we construct the galaxy number-count density field $`\rho _s(𝐬)`$ using a high-resolution Gaussian smoothing length of 500 $`\mathrm{km}\mathrm{s}^1`$. First of all, the LAP equations are solved for each survey following §2.1. This results in an estimate of the cosmic fields $`\delta `$, $`\alpha `$. We carry out a likelihood analysis to compute the value of $`\mathrm{\Omega }_\mathrm{m}`$ in each case, and the same analysis also yields an estimate of the parameters $`b_{rl}`$. In order to do so, the underlying assumption is that the errors in the measurements of the fields and the parameters, $`\mathrm{\Delta }\delta _y^{(n)}`$, $`\mathrm{\Delta }\alpha _y^{(n)}`$, $`\mathrm{\Delta }\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Delta }b_{rl}`$, follow a multivariate Gaussian distribution. However, we do not assume that the fields themselves are Gaussian (not at the present time, though we find that they are Gaussian at the initial time, to great accuracy, and therefore it is safe to assume Gaussian initial conditions in the likelihood analysis though our formalism does not require us to do so). The results of $`P(\mathrm{\Omega }_\mathrm{m})`$ for each survey are given in Table 2. The maximum likelihood estimates for $`\mathrm{\Omega }_\mathrm{m}`$ are $`\mathrm{\Omega }_\mathrm{m}^{\mathrm{PSC}z}=0.36`$, $`\mathrm{\Omega }_\mathrm{m}^{\mathrm{ORS}}=0.32`$, $`\mathrm{\Omega }_\mathrm{m}^{\mathrm{MarkIII}}=0.34`$, $`\mathrm{\Omega }_\mathrm{m}^{\mathrm{SFI}}=0.38`$. The 1$`\sigma `$ errors in these measurements are of the order $`_{00}^{1/2}`$, which results in $`2.1\times 10^2`$, $`1.7\times 10^2`$, $`1.3\times 10^2`$ and $`1.9\times 10^2`$, respectively.

### 5.1 Measurement of $`\mathrm{\Omega }_\mathrm{m}`$

The results obtained above are used to construct an initial Ansatz (14), (15) over the domain (13), that we use to solve (38) subject to (10). We follow the method set out in §4 and get the LAP solution that is the best fit to the four surveys. The likelihood analysis of the solution yields an estimate $`\mathrm{\Omega }_\mathrm{m}=0.37\pm 0.01`$ within the 1$`\sigma `$ level. This result is quite surprising in view of Table 2, due to the fact that the two surveys with the largest $`\eta `$, ORS and Mark III, seem to pin down $`\mathrm{\Omega }_\mathrm{m}`$ at a lower value, in the range $`0.320.34`$, and yet PSC$`z`$ and SFI predict a $`\mathrm{\Omega }_\mathrm{m}`$ that is closer to the final value. In this case, the optimization of the morphological features of the overdensity and velocity fields in the LAP solution with those of the four surveys is such that it favours a larger $`\mathrm{\Omega }_\mathrm{m}`$. Indeed the surveys PSC$`z`$ and Mark III are closer to the final solution than ORS and SFI. The latter two yield the highest $`\chi _i^2`$ in (10), though SFI is in better agreement with the LAP solution than ORS and therefore $`\chi _{\mathrm{SFI}}^2`$ is smaller. PSC$`z`$ and Mark III have the smaller variance of the four, and therefore this causes that they have a greater weight in determining the cosmography and the cosmological parameters in the solution.

The likelihood analysis for the bias parameters yields the estimates shown in Table 3. We have computed 42 coefficients, truncating the indices at $`r=5`$ and $`l=6`$. The higher $`l`$ enables us to probe into the nonlinear scales, though a stringent limit is set by the smoothing scale of 500 $`\mathrm{km}\mathrm{s}^1`$. The high $`l`$ coefficients are greater than those in the linear range (low $`l`$), showing greater bias at smaller scales, a trend that is also present in all four surveys when the LAP method is applied to each one individually.

### 5.2 Matter power spectrum

Hereafter we use $`\delta ^{\mathrm{LAP}}`$ and $`\alpha ^{\mathrm{LAP}}`$ evaluated at the maximum likelihood values of the parameters, i.e. $`\mathrm{\Omega }_\mathrm{m}=0.37`$ and $`b_{rl}`$ given in Table 3. We compute the Fourier modes $`\delta (𝐤)`$ of $`\delta ^{\mathrm{LAP}}`$. The matter power spectrum $`P(k)`$ is then given directly by

$$\delta (𝐤)\delta ^{}(𝐤^{})=(2\pi )^3P(k)\delta (𝐤𝐤^{}),$$
(16)
assuming that the modes $`\delta (𝐤)`$ are isotropic and homogeneous, and $`\delta `$ on the RHS is a Dirac delta function. Strictly speaking the fields ought be smooth, though we convolute the solutions with $`W(k)`$, in order to discard residual noise generated in the derivation of the LAP solution. The results for $`P(k)`$ are shown in Table 4. The galaxy-galaxy power spectrum is given by

$$P_{gg}(k)b(k)^2P(k),$$
(17)
and the corresponding results are shown in Table 5. We have consigned the analysis to the range of wavenumbers $`0.02h\mathrm{Mpc}^1\begin{array}{c}<\hfill \\ \hfill \end{array}k\begin{array}{c}<\hfill \\ \hfill \end{array}0.49h\mathrm{Mpc}^1`$ that we have binned in 24 equidistant segments in the log scale. The non-linear part of the power spectrum is located at $`k\begin{array}{c}>\hfill \\ \hfill \end{array}0.3h\mathrm{Mpc}^1`$, and it is in this range that we expect the LAP method to yield significant information. The resulting velocity dispersion is estimated $`\sigma _v286\pm 32`$ $`\mathrm{km}\mathrm{s}^1`$and $`\sigma _8^2=0.64\pm 0.03`$. A COBE-normalized flat $`\mathrm{\Lambda }`$CDM model, untilted and with $`\mathrm{\Omega }_\mathrm{m}=0.37`$ and $`\mathrm{\Lambda }=0.63`$, would fit the measured $`P_{gg}`$ for $`\mathrm{\Gamma }=0.21\pm 0.03`$ within the 1$`\sigma `$ error margin. The errors in Tables 4 and 5 are 1$`\sigma `$. In the case of Table 5, they are larger as they result from the product of the errors in $`P(k)`$ and $`b(k)`$. In addition to the 1$`\sigma `$ errors that are the result of the likelihood analysis, another source of error is due to the limited number of $`\delta _y`$ modes, which is propagated in the coordinate transformations $`\delta _y\delta _𝐤`$ and $`b_{rl}b(k)`$. These amount to $`\begin{array}{c}<\hfill \\ \hfill \end{array}10\%`$ of the total error, and have been accounted for.

Other power spectra of interest are (see e.g. Pen 1998; Tegmark & Peebles 1998; Dekel & Lahav 1999) the galaxy velocity power

$$P_{gv}(k)r(k)b(k)fP(k)$$
(18)
and the velocity-velocity power

$$P_{vv}(k)f^2P(k)$$
(19)
where $`|r(k)|1`$ is a galaxy-velocity correlation coefficient, $`f\mathrm{\Omega }_\mathrm{m}^{0.6}`$ is the linear growth rate. Typically $`r1`$ (Hamilton, Tegmark & Padmanabhan 2000), though it may depart from unity in the nonlinear regime. Using the measurement of $`\mathrm{\Omega }_\mathrm{m}`$ given in §5.1 and Table 4, it is straightforward to compute (19). The resulting errors result from the product of $`\mathrm{\Delta }P(k)`$ and $`\mathrm{\Delta }\mathrm{\Omega }_\mathrm{m}=0.01`$.

## 6 Discussion

The measurement of $`\mathrm{\Omega }_\mathrm{m}`$ that we have obtained by applying the LAP method simultaneously on PSC$`z`$, ORS, Mark III and SFI is $`\mathrm{\Omega }_\mathrm{m}=0.37\pm 0.01`$. This is consistent with the estimates that we have obtained from each survey separately (see Table 2), though ORS and Mark III predict a somewhat lower value, $`\mathrm{\Omega }_\mathrm{m}0.320.34`$. The main factors in the determination of $`\mathrm{\Omega }_\mathrm{m}`$ in our analysis are chiefly two: (1) the variance of the distribution $`P[\chi _i]`$ defined in §4, (2) the errors $`\chi _i^2`$, that reflect the resemblance of cosmographical features in the survey and the LAP soltuion. The ideal combination is a small variance and a minimal value of $`\chi _i^2`$, indicating a close resemblance between the survey and the LAP solution. The surveys that satisfy these requirements best are PSC$`z`$ and Mark III. The variance of ORS and SFI is not much greater than that of PSC$`z`$ and Mark III, so these surveys are consistent with the final LAP solution by an error of amplitude no greater than $`6\%`$.

These results indicate a higher value of $`\mathrm{\Omega }_\mathrm{m}`$ than reported in S01 for IRAS 1.2 Jy ($`\mathrm{\Omega }_\mathrm{m}=0.30`$ with a linear bias of $`b=1.1`$). In order to compare the results of this paper with other measurements reported in the literature, most of which are in the form of the redshift distortions parameter $`\beta \mathrm{\Omega }_\mathrm{m}^{0.6}/b`$, we need to have some information about the optimal value of the linear bias in those measurements. Hamilton (1998) surveys measurements up until mid-1997, most of which fall within the region $`\beta 0.450.75`$. A higher measurement of $`\beta =0.89\pm 0.12`$ is reported by Sigad et al. (1998) from Mark III, and recently a much lower one of $`\beta =0.41_{0.12}^{+0.13}`$ by Hamilton, Tegmark & Padmanabhan (2000) from PSC$`z`$ and a moderately low one of $`\beta =0.5\pm 0.1`$ by Nusser et al. (2000) from a comparison between ENEAR peculiar velocities and the PSC$`z`$ gravity field. Our results are consistent with $`\beta `$ measurements at the lower part of the range, as these allow a density parameter in the region $`0.30.4`$ for reasonable values of the linear bias $`b1.01.3`$. Higher values of $`\beta `$ require a much larger $`\mathrm{\Omega }_\mathrm{m}`$ that is not consistent with our results. Such is the case of Zehavi & Dekel (1999), who predict $`\mathrm{\Omega }_\mathrm{m}\begin{array}{c}>\hfill \\ \hfill \end{array}0.6`$ for both Mark III and SFI, and a value of $`\mathrm{\Omega }_\mathrm{m}0.4`$ that would be consistent with our results has a low likelihood in their analysis.

In Table 4 and 5 we have given the values of the matter and galaxy-galaxy power spectra respectively. The values of $`P_{gg}`$ are consistent, within the 1$`\sigma `$ errors given, with the results of Sutherland et al. (1999) and Hamilton & Tegmark (2000) for PSC$`z`$. The range of wavenumbers considered by Hamilton & Tegmark (2000) is considerable, $`0.01h\mathrm{Mpc}^1\begin{array}{c}<\hfill \\ \hfill \end{array}k\begin{array}{c}<\hfill \\ \hfill \end{array}300h\mathrm{Mpc}^1`$, and they carry out a careful analysis of errors by computing the prewhitened power spectrum. Using similar results, Hamilton, Tegmark & Padmanabhan (2000) fit their measured $`P_{gg}`$ to a COBE-normalized, untilted flat $`\mathrm{\Lambda }`$CDM, $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ (Eisenstein & Hu (1998)), within a range of wavenumbers $`0.01h\mathrm{Mpc}^1\begin{array}{c}<\hfill \\ \hfill \end{array}k\begin{array}{c}<\hfill \\ \hfill \end{array}1.00h\mathrm{Mpc}^1`$. Our results is in agreement with theirs, within the more limited range of wavenumbers that we have estimated, namely $`0.01h\mathrm{Mpc}^1\begin{array}{c}<\hfill \\ \hfill \end{array}k\begin{array}{c}<\hfill \\ \hfill \end{array}0.49h\mathrm{Mpc}^1`$, and our analysis favours a slightly higher value of $`\mathrm{\Omega }_\mathrm{m}`$ that, within the errors, is fairly close to $`\mathrm{\Omega }_\mathrm{m}0.4`$. Furthermore, our analysis predicts a value of the velocity dispersion $`\sigma _v=286\pm 32`$ $`\mathrm{km}\mathrm{s}^1`$, slightly lower than one of the values adopted by Sutherland et al. (1999), $`\sigma _v=300`$ $`\mathrm{km}\mathrm{s}^1`$. The analysis of peculiar velocities in Mark III and SFI, by Zaroubi et al. (1997) and Freudling et al. (1999) respectively (see also Zehavi & Dekel (1999)), yields a high-amplitude estimate of the power spectrum, $`P(k=0.1h\mathrm{Mpc}^1)\mathrm{\Omega }_\mathrm{m}^{1.2}=(4.5\pm 2.0)\times 10^3(h^1\mathrm{Mpc})^3`$. Our estimates are consistent with this (from Table 4, we have $`P(k=0.1h\mathrm{Mpc}^1)(4.63\pm 0.36)\times 10^3(h^1\mathrm{Mpc})^3`$), pinning down the values to greater accuracy. It is also the case that in our estimates of $`P_{gg}(k)`$ the error bars are significantly smaller than in Hamilton & Tegmark (2000) and Sutherland et al. (1999).

## acknowledgement

This research has been funded by the Eusko Jaurlaritza research fellowship BFI99.116 and in part at EHU by research grant UPV172.310-G02/99.

## Appendix A LAP equations

The cosmological perturbations are governed by (S01)

$$=\frac{1}{2}(1+\delta )𝐯^2+\alpha \left\{\dot{\delta }+[(1+\delta )𝐯]\right\}\varphi \delta \frac{1}{3}\mathrm{\Omega }_\mathrm{m}^1|\varphi |^2,$$
(20)
where $`\varphi `$ is the gravitational potential created by the matter density contrast $`\delta `$. The variational equations, namely $`\stackrel{~}{\delta }𝒮/\stackrel{~}{\delta }\delta =0`$ and $`\stackrel{~}{\delta }𝒮/\stackrel{~}{\delta }\alpha =0`$ ($`\stackrel{~}{\delta }`$ denotes a variation), yield:

$$𝐯=\alpha ,$$
(21)
$$^2\varphi =\frac{3}{2}\mathrm{\Omega }_\mathrm{m}\delta ,$$
(22)
and

$$\dot{\delta }+[(1+\delta )𝐯]=0,$$
(23)
$$\dot{\alpha }+\frac{1}{2}|\alpha |^2+\varphi =0.$$
(24)
The radius of the Universe $`a(t)`$ is scaled out alongside $`h`$ with our use of units of km/s for the distance scale. We seek solutions for $`\delta `$,$`\alpha `$ by utilizing the trial fields

$$\delta (t,𝐱)=\underset{y}{}\delta _y(t)j_l(k_rx)Y_{lm}(\theta ,\phi ),$$
(25)
$$\alpha (t,𝐱)=\underset{y}{}\alpha _y(t)j_l(k_rx)Y_{lm}(\theta ,\phi ),$$
(26)
where $`yrlm`$. By substituting (25),(26) into (21),(22) we get

$$𝐯(t,𝐱)=\underset{y}{}\left[\widehat{𝐱}\alpha _y^{}(t)j_l(k_rx)+\frac{1}{x}\widehat{𝐱}𝐉_y(\alpha )\right]Y_{lm}(\theta ,\phi ),$$
(27)
$$\varphi (t,𝐱)=\frac{3}{2}\mathrm{\Omega }_\mathrm{m}\underset{y}{}k_r^2\delta _y(t)j_l(k_rx)Y_{lm}(\theta ,\phi ),$$
(28)
where the coefficients $`\alpha _y^{}`$ and $`𝐉_y(\alpha )`$ are calculated in Appendix B. The time dependence of the coefficients $`\delta _y`$, $`\alpha _y`$ is modelled by a polynomial of degree $`N`$ that is a sum of Chebyshev polynomials, i.e.

$$\delta _y(t)=\underset{n=1}{\overset{N}{}}\delta _y^{(n)}T_n(t),$$
(29)
$$\alpha _y(t)=\underset{n=1}{\overset{N}{}}\alpha _y^{(n)}T_n(t).$$
(30)
The properties of $`T_n`$ are summarised in Appendix C. We do not consider the $`0`$th order term in (29), (30) as it is a constant term. The polynomials (29), (30) that permit us to approximate $`\delta _y`$, $`\alpha _y`$, in such a way that the weighed sum of the errors in the time interval considered is least, are given by the coefficients

$$\delta _y^{(n)}=(\delta _y|T_n)(T_n|T_n)^1,$$
(31)
$$\alpha _y^{(n)}=(\alpha _y|T_n)(T_n|T_n)^1,$$
(32)
where $`(|)`$ is the discrete version of the weighed product (58), i.e.

$$(f|g)\underset{m=1}{\overset{M}{}}w(t_m)f(t_m)g(t_m),$$
(33)
where $`w(t)`$ is (59) and $`M`$ is the number of $`t_m`$ divisions on the interval $`[0,1]`$. The boundary conditions (S01)

$$0=\underset{n=0}{\overset{N}{}}(1)^n\delta _y^{(n)},$$
(34)
$$\rho _s(𝐬)=x^2\left(\frac{N_{\mathrm{gals}}}{V}\right)\left[1+g_0(𝐱)\right]\left[1+\alpha _0^{\prime \prime }(𝐱)\right]^1.$$
(35)
$`g_0`$ is the present galaxy number-count density contrast and $`V`$ is the volume of the survey. The coefficients $`\delta _y`$ enter (35) through the bias relationship (6), that also introduces the coefficients $`b_{rl}`$. The evaluation of the modes at the present time in (35) is simplified by virtue of $`T_n(t_0)=1`$. Substituting the trial fields into equations (23),(24), we get

$$\underset{y}{}(\dot{\delta }_yk_r^2\alpha _y)j_l(k_rx)Y_{lm}(\mathrm{\Omega })=\underset{yy^{}}{}\{\alpha _y^{}\delta _y^{}^{}j_l(k_rx)j_l^{}(k_r^{}x)$$

$$+\frac{1}{x^2}[\widehat{𝐱}𝐉_y(\delta )][\widehat{𝐱}𝐉_y^{}(\alpha )]\}Y_{lm}(\mathrm{\Omega })Y_{l^{}m^{}}(\mathrm{\Omega }),$$
(36)
and

$$\underset{y}{}\left(\frac{3}{2}\mathrm{\Omega }_\mathrm{m}k_r^2\delta _y\dot{\alpha }_y\right)j_l(k_rx)Y_{lm}(\mathrm{\Omega })$$

$$=\frac{1}{2}\underset{y^{}}{}\{\alpha _y^{}\alpha _y^{}^{}j_l(k_rx)j_l^{}(k_r^{}x)$$
(37)
$$+\frac{1}{x^2}[\widehat{𝐱}𝐉_y(\alpha )][\widehat{𝐱}𝐉_y^{}(\alpha )]\}Y_{lm}(\mathrm{\Omega })Y_{l^{}m^{}}(\mathrm{\Omega }).$$

The coefficients $`𝐉_y(\delta )`$ are given in Appendix B, as in $`𝐉_y(\alpha )`$, with the substitution $`\alpha \delta `$. By integrating out all coordinates, it is possible, though arduous, to finally arrange all the coefficients $`\delta _y^{(n)}`$, $`\alpha _y^{(n)}`$ in the following inhomogeneous matrix system for each $`y`$

$$\left[\begin{array}{cc}𝒞_y^\alpha & 𝒞_y^\delta \\ 𝒮_y^\alpha & 𝒮_y^\delta \end{array}\right]\left[\begin{array}{c}\alpha _y^{(n)}\\ \delta _y^{(n)}\end{array}\right]=\left[\begin{array}{c}𝒟_{y^{}y^{\prime \prime }}^y\delta _y^{}^{(m)}\alpha _{y^{\prime \prime }}^{(m^{})}\\ _{y^{}y^{\prime \prime }}^y\alpha _y^{}^{(m)}\alpha _{y^{\prime \prime }}^{(m^{})}\end{array}\right],$$
(38)
where $`𝒞_y^\delta `$, $`𝒞_y^\alpha `$, $`𝒮_y^\delta `$, $`𝒮_y^\alpha `$, $`𝒟_{y^{}y^{\prime \prime }}^y`$ and $`_{y^{}y^{\prime \prime }}^y`$ are computed using the orthogonality relation of the Chebyshev polynomials (58), and the standard orthogonality relations for $`Y_{lm}`$ and $`j_l`$,(65), (71). These terms are estimated numerically. The coefficients $`\alpha _y`$, $`\delta _y`$ in the column matrix on the LHS of (38) are arranged vertically in order of increasing $`n`$, from top to bottom, making a total of $`2N`$ entries. The submatrices $`𝒞`$ result from the integration of the LHS of (36), and the submatrices $`𝒮`$ from integrating the LHS of (37). The RHS of (38) contains cross coefficients coupling different $`y`$ and $`n`$. In conclusion, we have a set of constraints for $`\delta _y^{(n)}`$, $`\alpha _y^{(n)}`$: (34), (35) and (38)(as we have integrated out all coordinates from the dynamical equations, (38) is really an inhomogeneous system of constraints).

## Appendix B Radial derivatives

The radial derivative of the velocity potential can be written as

$$\frac{\mathrm{d}}{\mathrm{d}x}\alpha (t,𝐱)=\underset{y}{}\alpha _y^{}(t)j_l(k_rx)Y_{lm}(\theta ,\phi ),$$
(39)
where, using the equality

$$\frac{\mathrm{d}}{\mathrm{d}u}j_l(u)=(2l+1)^1\left[lj_{l1}(u)(l+1)j_{l+1}(u)\right],$$
(40)
we have

$$\alpha _y^{}=k_r\left[\frac{(l+1)}{(2l+3)}\alpha _{r(l+1)m}\frac{l}{(2l1)}\alpha _{r(l1)m}\right].$$
(41)
Similarly

$$\alpha _y^{\prime \prime }=k_r^2\{\frac{(l+1)}{(2l+3)}\frac{(l+2)}{(2l+5)}\alpha _{r(l+2)m}$$

$$\left[\frac{(l+1)^2}{(2l+3)(2l+1)}+\frac{l^2}{(2l1)(2l+1)}\right]\alpha _y$$
(42)
$$+\frac{l}{(2l1)}\frac{(l1)}{(2l3)}\alpha _{r(l2)m}\}.$$

The coefficients $`𝐉_y(\alpha )`$ in (27) are therefore given by

$$𝐉_y(\alpha )=\frac{\alpha (l,m+1)}{2}\alpha _{rl(m+1)}(i\widehat{𝐱}_1\widehat{𝐱}_2)$$

$$+\frac{\beta (l,m1)}{2}\alpha _{rl(m1)}(i\widehat{𝐱}_1+\widehat{𝐱}_2)+im\alpha _y\widehat{𝐱}_3,$$
(43)
where

$$\alpha (l,m)=\left[l(l+1)m(m1)\right]^{1/2},$$
(44)
$$\beta (l,m)=\left[l(l+1)m(m+1)\right]^{1/2}.$$
(45)
## Appendix C Chebyshev polynomials

The Chebyshev polynomials of the first kind $`T_n`$ satisfy the linear homogeneous second-order differential equation

$$(1t^2)\ddot{y}t\dot{y}+n^2y=0,$$
(46)
where the dot denotes $`\mathrm{d}/\mathrm{d}t`$. Following the convention of Courant & Hilbert (1989), they are explicitly given by

$$T_n(t)=\frac{n}{2}\underset{m=0}{\overset{[n/2]}{}}(1)^m\frac{(nm1)!}{m!(n2m)!}(2t)^{n2m},$$
(47)
that can be rewritten as

$$T_n(\mathrm{cos}\theta )=\mathrm{cos}(n\theta ).$$
(48)
The polynomials satisfy the orthogonality relation

$$T_n|T_m_1^1dtw(t)T_n(t)T_m(t)=\delta _{nm}\frac{\pi }{2}(1+\delta _{n0}),$$
(49)
where $`\delta _{nm}`$ is a Kronecker delta function and

$$w(t)=(1t^2)^{1/2}.$$
(50)
They have the following property

$$T_n(t)=T_n^{}\left(\frac{1+t}{2}\right),$$
(51)
and their recurrence relation is given by

$$2T_n(t)T_m(t)=T_{n+m}(t)+T_{nm}(t),$$
(52)
for $`nm`$. In the asymptotic regime $`t\mathrm{}`$, the zeros of $`T_n`$ are

$$t_m^{(n)}\mathrm{cos}\left(\frac{2m1}{2n}\pi \right),$$
(53)
where $`t_1^{(n)}<t_2^{(n)}<\mathrm{}`$. (49)-(53) are required in the evaluation of the numerical coefficients in (38).

The interval of orthogonality adopted in (49) is $`1t1`$. In our application it is more convenient to scale $`t`$ down to the interval $`[0,1]`$, and we do so by transforming the polynomials via the change of coordinate $`t=2\stackrel{~}{t}1`$. Therefore

$$\stackrel{~}{T}_n(t)=T_n(2t1)=\underset{m=0}{\overset{[n/2]}{}}\stackrel{~}{a}_mt^m.$$
(54)
The coefficients $`\stackrel{~}{a}_m`$ of the new polynomials are found recursively from the old $`a_m`$ through the relations

$$\stackrel{~}{a}_m^{(j)}=2a_m^{(j1)}a_{m+1}^{(j)},$$
(55)
for $`m=n1,n2,\mathrm{},j;`$ $`j=0,1,2,\mathrm{},n;`$

$$\stackrel{~}{a}_m^{(1)}=\frac{a_m}{2},$$
(56)
for $`m=0,1,2,\mathrm{},n;`$

$$a_m^{(j)}=2^ja_n,$$
(57)
for $`j=0,1,2,\mathrm{},n;`$ and $`a_m^{(m)}=\stackrel{~}{a}_m`$ for $`m=0,1,2,\mathrm{},n`$. Thus, the orthogonality relation (49) in the new domain becomes

$$\stackrel{~}{T}_n|\stackrel{~}{T}_m_0^1dt\stackrel{~}{w}(t)\stackrel{~}{T}_n(t)\stackrel{~}{T}_m(t)=\delta _{nm}\frac{\pi }{4}(1+\delta _{n0}),$$
(58)
where

$$\stackrel{~}{w}(t)=\frac{1}{2}t^{1/2}(1t)^{1/2}.$$
(59)
Throughout our calculations we use $`\stackrel{~}{T}_n`$ and the cross-product defined by (58), as opposed to (49), though we have omitted the tildes for simplicity.

## Appendix D Orthogonality relations for $`Y_{lm}`$ and $`j_l`$

### D.1 Spherical harmonics

A spherical harmonic is defined as (following the Condon-Shortley convention)

$$Y_{lm}(\theta ,\phi )=\sqrt{\frac{2l+1}{4\pi }\frac{(l|m|)!}{(l+|m|)!}}P_l^{|m|}(\mathrm{cos}\theta )\mathrm{e}^{\mathrm{i}m\phi }(1)^r,$$
(60)
where $`r=m`$ for $`m0`$ and $`r=0`$ otherwise, and the $`P_l^m(x)`$ are the associated Legendre functions, which are given by

$$P_l^m(x)=\frac{1}{2^ll!}(1x^2)^{m/2}\frac{\mathrm{d}^{l+m}}{\mathrm{d}x^{l+m}}(x^21)^l.$$
(61)
Our convention follows Courant & Hilbert (1969) and differs from Abramowitz and Stegun (1965) by a factor $`(1)^m`$. (61) can be otherwise rewritten as

$$P_l^m(x)=\frac{(1)^m}{\mathrm{\Gamma }(1l)}\left(\frac{x+1}{x1}\right)^{m/2}F(l,l+1;1l;\frac{1x}{2}),$$
(62)
where $`\mathrm{\Gamma }`$ is the gamma function and $`F`$ is the hypergeometric function. The associated Legendre functions satisfy the orthogonality condition

$$_1^1P_l^m(x)P_n^m(x)dx=\frac{2}{2l+1}\frac{(l+m)!}{(lm)!}\delta _{ln}.$$
(63)
The spherical harmonic satisfies the following property

$$Y_l^m(\theta ,\phi )=(1)^mY_l^m(\theta ,\phi ),$$
(64)
and its orthonormal relation is given by

$$_0^{2\pi }d\phi _0^\pi \mathrm{d}(\mathrm{cos}\theta )Y_{lm}(\theta ,\phi )Y_{l^{}m^{}}^{}(\theta ,\phi )=\delta _{ll^{}}\delta _{mm^{}}.$$
(65)
Furthermore, it satisfies the equality

$$P_l(\mathrm{cos}\gamma )=\frac{4\pi }{2l+1}\underset{m=l}{\overset{l}{}}Y_{lm}(\theta ,\phi )Y_{lm}^{}(\theta ^{},\phi ^{}),$$
(66)
where the directions $`(\theta ,\phi )`$ and $`(\theta ^{},\phi ^{})`$ are separated by an angle $`\gamma `$, and $`P_l`$ are the Legendre polynomials, defined as $`P_lP_l^0`$.

### D.2 Spherical Bessel functions

The spherical Bessel function of the first kind, $`j_n`$, is defined as

$$j_l(z)=\sqrt{\frac{\pi }{2z}}J_{l+\frac{1}{2}}(z),$$
(67)
where $`J_\nu (z)`$ is the Bessel function of the first kind

$$J_l(z)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!(l+n)!}\left(\frac{1}{2}z\right)^{l+2n},$$
(68)
that is one of two linearly independent solutions of the second-order differential equation

$$z^2w^{\prime \prime }+zw^{}+(z^2l^2)w=0,$$
(69)
where the prime denotes $`\mathrm{d}/\mathrm{d}z`$. One can rewrite (67) in the so-called Gegenbauer’s generalization

$$j_l(z)=\frac{1}{2}(\mathrm{i})^l_0^\pi \mathrm{e}^{\mathrm{i}z\mathrm{cos}\theta }P_l(\mathrm{cos}\theta )\mathrm{sin}\theta \mathrm{d}\theta .$$
(70)
The spherical Bessel function satisfies the following orthogonality relation

$$_0^1\mathrm{d}xx^2j_l(k_rx)j_l(k_sx)=\frac{\pi }{4k_r}\delta _{rs}\{J_{l+\frac{1}{2}}^{}(k_r)^2$$

$$+[1\left(\frac{l+1/2}{k_r}\right)^2]J_{l+\frac{1}{2}}(k_r)^2\},$$
(71)
where $`\delta _{rs}`$ is a Kronecker delta function and $`k_r0`$, and also a closure relation given by

$$_0^{\kappa _{\mathrm{max}}}d\kappa \kappa ^2j_l(\kappa x)j_l(\kappa x^{})=\frac{\pi }{2x^2}\delta (xx^{}),$$
(72)
where $`\delta `$ is in this case a Dirac delta function.


format: text

text:
# FUSE Limits on FUV Emission from Warm Gas in Clusters of Galaxies1footnote 11footnote 1Based on observations made with the NASA-CNES-CSA Far Ultraviolet Spectroscopic Explorer. FUSE is operated for NASA by the Johns Hopkins University under NASA contract NAS5-32985.

## 1 INTRODUCTION

Recent observations with the Extreme Ultraviolet Explorer (EUVE) provide evidence that a number of clusters of galaxies are strong emitters of EUV radiation (Lieu et al. 1996a, b; Bowyer et al. 1997; Mittaz et al. 1998). This emission is substantially brighter than would be expected from the well-known, X-ray–emitting intracluster medium (ICM; Forman & Jones 1982). Signatures of this “soft excess” are sometimes present in the lowest energy-resolution band of the ROSAT PSPC, but its presence remains controversial (Arabadjis & Bregman 1999; Drake 1999; Krick et al. 2000). At present, two groups agree that an EUV excess exists in both the Virgo and Coma clusters, but they differ on the intensity and morphology of the emitting regions (Lieu et al. 1999; Enßlin et al. 1999; Bowyer & Berghöfer 1998; Bowyer et al. 1999; Berghöfer et al. 2000).

The EUV excess in clusters of galaxies was originally attributed to a diffuse, $`(510)\times 10^5`$ K thermal gas component of the ICM. Gas at these temperatures cools rapidly, however, requiring either a substantial mass of cooling gas or a heretofore unknown energy source to reheat it (Fabian 1996). If a reservoir of warm gas were present in the cores of clusters of galaxies, it would emit strongly in the far-UV (FUV) resonance lines of O VI $`\lambda \lambda 1032,1038`$ and C IV $`\lambda \lambda 1548,1551`$ as it cools through temperatures of a few times $`10^5`$ K (Edgar & Chevalier 1986; Voit et al. 1994), yet Dixon et al. (1996) were unable to detect either doublet in the spectra of five clusters of galaxies obtained with the Hopkins Ultraviolet Telescope (HUT). To probe more deeply for this emission, we have observed the Virgo and Coma clusters with the Far Ultraviolet Spectroscopic Explorer (FUSE). Our limits on the O VI flux of each cluster place tight constraints on the emission integral of any warm component of their ICM.

## 2 OBSERVATIONS AND DATA REDUCTION

FUSE comprises four separate optical systems. Two employ LiF optical coatings and are sensitive to wavelengths from 990 to 1187 Å, while the other two use SiC coatings, which provide reflectivity to wavelengths as short as 905 Å. The four channels overlap between 990 and 1070 Å. For a complete description of FUSE, see Moos et al. (2000) and Sahnow et al. (2000).

The FUSE spectrum of the Coma cluster was obtained in 17 separate exposures on 2000 June 18 and 19. Each was centered on 12h59m49$`\stackrel{\mathrm{s}}{\mathrm{.}}`$0, +27°57′46″ (J2000, or $`l=57.61,b=+87.96`$ in Galactic coordinates), near the cluster center. The total exposure was 28608 s, with 23553 s obtained during orbital night. We use the entire 29-ksec data set in our analysis. Our spectrum of the Virgo cluster is a combination of data from two locations near the cluster center. Two exposures, centered on 12h31m07$`\stackrel{\mathrm{s}}{\mathrm{.}}`$3, +12°23′46″ ($`l=284.03,b=+74.52`$) and totaling 2242 s, were obtained on 2000 June 13. Seven exposures, centered on 12h31m13$`\stackrel{\mathrm{s}}{\mathrm{.}}`$4, +12°22′10″ ($`l=284.2,b=+74.50`$) and totaling 8688 s, were obtained on 2000 June 17. The total integration time is 10,930 s, all of it during orbital night. All observations were made through the 30″ $`\times `$ 30″ (LWRS) aperture.

The standard FUSE data-reduction software are used to confirm that none of our data were obtained during passages through the South Atlantic Anomaly or at low earth-limb angles, to separate day and night observations, and to confirm that the dead-time correction for both observations is negligible. To reduce the detector background, photon events with pulse heights less than 4 or greater than 15 (in standard arbitrary units) are excluded from the data set. Subsequent data reduction is performed by hand. Photon events from regions of the detector with known defects are removed. The standard extraction window defined for the LWRS aperture is centered (in the dimension perpendicular to the dispersion axis) on the diffuse airglow emission lines.

While the resolution of the FUSE spectrograph is approximately 15 km s-1 for a point source, diffuse emission filling the LWRS aperture yields a line profile that is well approximated by a top-hat function with a width of $``$ 106 km s-1. Our spectra are quite faint, so we bin the Coma data by 8 detector pixels, or about 15.5 km s-1. We bin the data from the shorter Virgo observation by 16 pixels. We distinguish between detector pixels and 8- or 16-pixel bins throughout this paper.

We omit the final steps of the standard FUSE data-reduction pipeline for several reasons: first, the pipeline’s corrections for differential Doppler shifts and grating and electronics drifts total less than 10 detector pixels in the dispersion direction. Second, the distortion correction, applied independently for each of the instrument apertures, can cause background photons to “pile up” in the region between two nominal extraction windows, forming a horizontal stripe in the processed detector image. These stripes occasionally overlap our preferred extraction window. For bright point sources, this effect is unimportant, but for faint, diffuse sources, we prefer not to risk including such features in our extracted spectrum. Finally, we wish to include the detector background as a free parameter in our spectral models, so forgo the background subtraction steps included in the standard pipeline.

The FUSE flux calibration, based on theoretical models of white-dwarf stellar atmospheres, is believed accurate to about 10% (Sahnow et al. 2000). Corrections to the nominal FUSE wavelength scale are derived from the measured positions of airglow features in each detector segment and are good to about 0.01 Å. In practice, only the portion of the Coma spectrum obtained during orbital day contains airglow features strong enough to allow an independent wavelength correction. We thus apply the Coma correction to the Virgo data; the resulting wavelength solution is consistent (within the errors) with the positions of measureable airglow lines. Error bars are assigned to the data assuming Gaussian statistics, then smoothed by 9 bins to remove small-scale features in the error spectrum without significantly changing its shape. Segments of the flux- and wavelength-calibrated LiF 1A spectra, showing O VI $`\lambda \lambda 1032,1038`$ at the redshift of each cluster, are presented in 1. We assume $`z=0.0231`$ for Coma (Struble & Rood 1999) and $`z=0.0036`$ for Virgo (Ebeling et al. 1998). Diffuse emission from our own Galaxy is present in both spectra; these results are presented in Dixon et al. (2001).

## 3 SPECTRAL ANALYSIS

33footnotetext: The Image Reduction and Analysis Facility (IRAF) is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc., (AURA) under cooperative agreement with the National Science Foundation.

In our analysis, model spectra are fit to the flux-calibrated data using the nonlinear curve-fitting program SPECFIT (Kriss 1994), which runs in the IRAF22footnotemark: 2 environment, to perform a $`\chi ^2`$ minimization of the model parameters. We find that, because of the wide disparity in the effective area of the various detector segments, combining spectra from different segments does not significantly improve the signal-to-noise ratio. We thus use data only from the segment with the highest effective area at the wavelength of interest.

To set limits on the flux of an emission feature, we first fit a linear continuum to the spectrum. We then add a synthetic emission feature to the data array, raising its flux until $`\chi ^2`$ for the best-fit model with an emission feature differs from $`\chi ^2`$ for the best-fit model without an emission feature by $`\mathrm{\Delta }\chi ^2=9`$ (corresponding to a 3-$`\sigma `$ deviation for one interesting parameter, in this case the flux in the emission line; Avni 1976). We quote as a 3-$`\sigma `$ upper limit the flux of the model emission line that best fits the synthetic emission feature originally added to the data. Upper limits derived in this way are unaffected by small-scale features in the observed continua. The best-fit linear continuum and our 3-$`\sigma `$ upper limits to the flux of the O VI doublet for each cluster are shown in 1.

The observed profile of a diffuse emission feature represents a convolution of its intrinsic profile with the 106 km s-1 top-hat function discussed above. In our models, we assume the intrinsic profile to be a Gaussian with FWHM = 40 km s-1, corresponding to the thermal width of a $`5\times 10^5`$ K gas. Broader lines, perhaps due to bulk motions in the gas, would yield higher upper limits, as the flux is spread over a greater region of the detector. For example, lines with an intrinsic FWHM = 200 km s-1 yield upper limits approximately 70% higher than those of 40 km s-1 lines.

No significant FUV line emission is seen at the cluster redshift in either spectrum. Our limits on the fluxes of various lines predicted to be strong in thermal-gas (Landini & Monsignori Fossi 1990), mixing-layer (Slavin et al. 1993), and shock (Hartigan et al. 1987) models are presented in Table 1. Dereddened intensities are calculated assuming the extinction parameterization of Cardelli et al. (1989, hereafter CCM) with $`E(`$B V) = 0.008 for Coma and 0.030 for Virgo (Schlegel et al. 1998) and $`R_V=3.1`$. Uncertainties in the reddening correction are discussed in Dixon et al. (1996).

At the wavelength of redshifted O VI in the Coma cluster (1055.77, 1061.59 Å), the LiF 1A detector channel has an effective area $`A_{\mathrm{eff}}=24.3`$ cm-2, by far the highest of the four FUSE channels (Sahnow et al. 2000). Unfortunately, a bright band of scattered light contaminates the LiF 1A spectrum from about 1045 to 1058 Å. Present even at night, the band is clearly visible in 1. It shows considerable structure, including an apparent emission feature at 1054.5 Å. Because the intensity of the band varies on small spatial scales, background spectra extracted from detector regions above and below our spectral window cannot be used to predict its shape in our spectrum. The statistical significance of the 1054.5 Å feature is about 3.3 $`\sigma `$; if real, it should be present in the spectrum from the less-sensitive ($`A_{\mathrm{eff}}=15.5`$ cm-2) LiF 2B channel at a significance of $`>`$ 2 $`\sigma `$. No such feature is seen, and we conclude that it is an artifact of the stray-light stripe in the LiF 1A detector segment.

Because the FUSE SiC channels lie on the sun-illuminated side of the spacecraft, their spectra may be contaminated by solar emission (Shelton et al. 2001). We thus use only data obtained during orbital night to set limits on emission features in these spectra. In practice, only the S VI $`\lambda 933.38`$ line in our Coma spectrum is affected, as all other features in the spectrum are redshifted onto the LiF detector channels. All of the Virgo data were obtained during orbital night.

In our Virgo spectrum, the C II\* $`\lambda 1037.02`$ and O VI $`\lambda 1037.62`$ lines are redshifted to 1040.75 and 1041.36 Å, respectively, and, if present, might be contaminated by flux from the O I airglow features at 1040.94 and 1041.688 Å. To set limits on the flux of C II\* and O VI, we include the airglow lines in our model spectrum, but limit their fluxes to that observed in the O I $`\lambda 1027.43`$ line, which is generally much brighter than either of the two longer-wavelength features. The resulting limits are higher than those for the nearby Ly $`\beta `$ and O VI $`\lambda 1031.93`$ lines. Fortunately, all of the Virgo data were obtained during orbital night, when the O I line flux is minimal.

Finally, we set upper limits on the FUV continuum emission from each cluster of galaxies. We consider only the region between 1060 and 1070 Å, which is free of airglow features and for which the LiF 1A detector has $`A_{\mathrm{eff}}>20`$ cm2. For each observation, we determine the mean background-event rate from unilluminated regions of the detector directly above and below the LWRS extraction window. Scaling these rates to the size of the LWRS window, we find that the observed continuum exceeds the detector background by $`110\pm 250`$ and $`200\pm 400`$ photons cm-2 s-1 Å-1 sr-1, respectively, for Coma and Virgo. (We ignore possible structure in the detector background, thought to be present at the 10% level.)

## 4 DISCUSSION

### 4.1 The Coma Cluster

From simultaneous fits to EUVE Deep Survey (DS) and ROSAT PSPC observations of the Coma cluster, Lieu et al. (1996a) derive the temperature and emission integral as a function of radius for a three-component ICM assuming a distance of 139 Mpc and an abundance $`Z=0.21`$ $`Z_{\mathrm{}}`$ (Hughes et al. 1993). They find a temperature of $`8.7\times 10^5`$ K and an emission integral of $`4.8\times 10^{65}`$ cm-3 for the coolest component within $`3\mathrm{}`$ of the cluster center. Combining these parameters with line emissivities estimated with the CHIANTI software package (Dere et al. 1997) using the Feldman et al. (1992) solar abundances scaled to the assumed cluster metallicity, we derive an expected O VI $`\lambda 1032`$ flux of $`3.8\times 10^3`$ photons cm-2 s-1. If the O VI emission were evenly distributed within this region (an assumption consistent with both EUV and soft X-ray maps of the cluster; Bonamente 2000; Briel et al. 1992), then its surface brightness would be 1600 photons cm-2 s-1 sr-1, slightly less than our 3-$`\sigma `$ upper limit to the dereddened O VI surface brightness.

The O VI line emissivity depends critically on the gas temperature. From our limit to the O VI flux of Coma, we derive an upper limit to the emission integral as a function of temperature. In 2, we compare this limit with the best-fit value of the emission integral derived by Lieu et al. (1996a) for the central region ($`r<3\mathrm{}`$) of the cluster, scaled to the area of the LWRS slit. This value of the emission integral, $`EI=4.2\times 10^{63}`$ cm-3, lies just below the FUSE 3-$`\sigma `$ upper limit. If we instead use a 2-$`\sigma `$ limit to the O VI flux, we can nominally exclude the warm-gas model of Lieu et al. (1996a). More recent models of the Coma EUV excess do not require a sub-$`10^6`$ K gas. Bonamente (2000) fits the data with a two-component ICM, the cooler of which has a temperature of $`2.9\times 10^6`$ K. We would not expect significant O VI emission from such a hot gas.

### 4.2 The Virgo Cluster

Two groups (Lieu et al. 1996b; Bonamente et al. 2001) have modeled the EUV excess in the Virgo cluster with a warm component of the ICM. Both find that its temperature rises with radius, while its emission integral falls. The brightest FUV emission would thus be expected from the cluster center. Unfortunately, the center of the Virgo cluster is occupied by the galaxy M87, and its strong stellar continuum complicates the FUV spectrum of this region (cf., Dixon et al. 1996). We thus consider data taken at larger radii, about 4$`\stackrel{}{\mathrm{.}}`$4 (2242 s) and 6$`\stackrel{}{\mathrm{.}}`$0 (8688 s) from the cluster center. Both Lieu et al. (1996b) and Bonamente et al. (2001) model the warm ICM using a series of concentric annuli. Our inner data point samples the 3–5 arcmin annulus, while our outer point falls in the 5–7 arcmin region. Because 80% of our data comes from the outer annulus, we compare our observations with model predictions for the 5–7 arcmin annulus. We estimate that errors in our derived limits to the emission integral resulting from this simplification are on the order of the uncertainties in the FUSE flux calibration.

Our upper limit to the emission integral as a function of temperature is plotted in 3. Also shown are the best-fit values of the temperature and emission integral of the warm ICM, derived from simultaneous fits to the DS and PSPC data from the 5–7 arcmin annulus of the Virgo cluster by Lieu et al. (1996b) and Bonamente et al. (2001), respectively, scaled to the area of the LWRS slit. The Lieu et al. model parameters are $`T=5.8\times 10^5`$ K, $`EI=9.35\times 10^{61}`$ cm-3, and $`Z=0.454`$ $`Z_{\mathrm{}}`$, values firmly excluded by our limit to the O VI flux. The Bonamente et al. model employs a much warmer gas, with $`T=9.2\times 10^5`$ K, $`EI=4.72\times 10^{61}`$ cm-3, and $`Z=0.5`$ $`Z_{\mathrm{}}`$; these parameters fall just below the FUSE 3-$`\sigma `$ upper limit. If we instead use a 2-$`\sigma `$ limit to the O VI flux, we can nominally exclude this model as well.

The redshift of M87 ($`z=0.00436`$; Smith et al. 2000) is slightly greater than that of the Virgo cluster mean. Near the center of the cluster, the X-ray emission is well centered on M87, and certain asymmetrical features seem correlated with radio emission from the nucleus (Belsole et al. 2001). If, in this region, the cooling flow emits at the redshift M87, the O VI $`\lambda 1032`$ line would be strongly attenuated by absorption from Galactic C II $`\lambda 1036.3`$. An upper limit to the emission integral derived from the O VI $`\lambda 1038`$ line would be approximately twice our quoted value.

## 5 CONCLUSIONS

We have obtained FUV spectra of two clusters of galaxies with FUSE, the Far Ultraviolet Spectroscopic Explorer. Neither spectrum shows significant O VI $`\lambda \lambda 1032,1038`$ emission at the cluster redshift. Our 2-$`\sigma `$ upper limits on the O VI $`\lambda 1032`$ flux from Coma and Virgo exclude all published warm-gas models of the EUV excess in these clusters and severely constrain the amount of gas at $`T8\times 10^5`$ K in either cluster.

We thank R. Shelton and E. Murphy for discussions of FUSE data analysis and D. Buote for thoughtful comments on the manuscript. We acknowledge the outstanding efforts of the FUSE P.I. team to make this mission successful. This research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Extragalactic Database (NED), operated by JPL for NASA. This work is supported by NASA grant NAG 5-5197.