text:
# Using Perturbative Least Action to Reconstruct Redshift Space Distortions

## 1 Introduction

The current generation of galaxy redshift surveys is producing an avalanche of data about the structure of the universe. The Sloan Digital Survey (SDSS; York et al. 2000) has already measured the redshifts of $`30,000`$ galaxies (X. Fan, private communication) of the $`10^6`$ galaxy redshifts over 10,000 $`\mathrm{deg}^2`$ which the survey will ultimately cover. The Two Degree Field redshift survey (2dF; Colless 1999) has measured redshifts for $`10^5`$ galaxies, and will eventually measure a quarter of a million galaxies over 2000 $`\mathrm{deg}^2`$ in the southern hemisphere. Meanwhile, the IRAS 0.6 Jy Point Source Catalog redshift Survey (PSCz; Saunders et al. 2000), with about $`15,000`$ galaxies over virtually all of the sky, provides a fertile testbed for cosmological models and methods.

As impressive as these surveys are, they are limited to providing a somewhat distorted snapshot of the universe. For example, it is both suggested observationally (Hubble 1936; Oemler 1974; Davis & Geller 1976; Kaiser 1984; Santiago & Strauss 1992; Blanton 2000 and references therein) and predicted (Davis et al. 1985; Bardeen et al. 1986; Blanton 1999 and references therein; Dekel & Lahav 1999) that the luminous structure, that which the surveys record, is biased with respect to the underlying matter density of the universe. Additionally, though the Hubble relation can be used to give an approximate 3-dimensional picture of structure, peculiar velocities (e.g. Strauss & Willick 1995 and references therein for a review) cause a distortion of the structure along the line of sight. Since the peculiar velocity field is a function of the underlying mass field and cosmology, velocity distortions and bias are intimately related.

Perturbative Least Action (PLA; Goldberg & Spergel 2000, hereafter GS) was shown to be an excellent technique for the reconstruction of nonlinear structure in real space. In this paper, we extend PLA into redshift space, and show that it is an excellent tool for extracting information from redshift surveys even into the nonlinear regime.

However, before going too far afield, it will be useful to review some of the basic issues involved in redshift space distortions of the density field, and define some of the symbols which will be used throughout this paper. This discussion is not meant to be comprehensive, however, and the interested reader will certainly benefit from some of the excellent reviews on the subject (Hamilton 1998; Hatton & Cole, 1998; Zaroubi & Hoffman 1996; Strauss & Willick 1995; Sahni & Coles 1995; Kaiser 1987).

### 1.1 Definitions and Conventions

Let us consider an observer in an expanding universe. Hubble’s law states that in a Friedman-Robertson-Walker universe, the distance, $`d`$ to a test particle with redshift, $`z`$, will be, to first order in $`z`$:

$$czH_0d^{FRW},$$
(1)
where $`H_0`$ is the Hubble constant at the present epoch, $`d`$ is the distance to the particle, and the approximation comes from the fact that relativistic effects become important at high redshifts. However, we will confine our discussion to the non-relativistic case and this approximation throughout this discussion.

Hubble’s Law assumes that a particle (galaxy) is at rest in comoving coordinates. A particle with a local, peculiar velocity, $`𝐯`$, will have a redshift of:

$$cz=H_0d+\widehat{𝐱}𝐯,$$
(2)
where $`\widehat{𝐱}`$ points along the line of sight of the test particle. Since it is actually this redshift that we observe, and not the position of the particle, it is worthwhile to construct a comoving coordinate which reflects the observations of the observer at the origin. We define the comoving redshift space coordinate:

$$𝐬=𝐱+\widehat{𝐱}\left(\frac{a\widehat{𝐱}\dot{𝐱}}{\dot{a}}\right),$$
(3)
where we have substituted $`a\dot{𝐱}`$ for $`𝐯`$ and $`a𝐱`$ for $`d\widehat{𝐱}`$ in our implicit definition of the peculiar velocity. Comparison with equation (2) yields the relation, $`s(t)=cz/\dot{a}(t)`$.

Equation (3) shows that the mapping of $`𝐱`$ to $`𝐬`$ is inherently non-invertible. Any trajectory of $`𝐱`$ will yield a single trajectory in redshift space, but the converse does not necessarily hold.

Indeed, even in the Zel’dovich regime (defined below), the infall of matter from both sides (front and back) of a structure can give rise to a “triple-valued zone” (see Strauss & Willick 1995, §5 for a discussion), as matter at different physical distances from the observer appear to have identical redshifts due to conspiracy between the Hubble flow and the peculiar velocity. We will attempt to disentangle this degeneracy in § 2.

### 1.2 The Linear Regime

For now, let us consider the linear regime. By convention, particle $`i`$ sits at position, $`𝐪_i`$, at $`a=0`$, and that the ensemble of $`𝐪_i`$ forms a uniform grid. If a field is linear, that is, if

$$\delta (𝐱,t)\frac{\rho (𝐱,t)}{\overline{\rho }(t)}1$$
(4)
remains small at all times, and if the velocity field initially has no curl, then the Zel’dovich approximation (Zel’dovich 1970) can be used to give the trajectory of a particle as:

$$𝐱_i(t)=𝐪_i+D(t)𝐩_i,$$
(5)
where $`D(t)`$ is a cosmology dependent, monotonically increasing growth factor, normalized to unity at the present, and $`𝐩_i`$ is the final displacement of particle, $`i`$, from its initial position.

Fields for which equation (5) well approximates the trajectories of all particles at all times will be referred to as Zel’dovich fields. This is to be distinguished from fields in the “linear regime” for which:

$$\delta (𝐪_i,t)=D(t)_q𝐩(𝐪_i)$$
(6)
holds at all times. When perturbations are very small, both equalities hold. However, as perturbations become larger, equation (6) breaks down first, and densities evolve according to a more complex function of time.

In the Zel’dovich regime, the redshift space coordinate evolves as:

$$𝐬(t)^{(zel)}=𝐪+D(t)𝐩+\widehat{𝐱}\left[\widehat{𝐱}\left(\frac{a(t)\dot{D}(t)}{\dot{a}(t)}𝐩\right)\right]=𝐪+D(t)𝐩+\widehat{𝐱}\left[D(t)f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t)\widehat{𝐱}𝐩\right]$$
(7)
where we define $`f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t)`$ such that:

$$f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t)\frac{a(t)\dot{D}(t)}{\dot{a}(t)D(t)}.$$
(8)
For $`\mathrm{\Omega }_M=1`$, this function is a constant in time. At $`t=t_0`$, a good analytic fit can be given by $`f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)\mathrm{\Omega }_M^{0.6}`$ (Peebles 1980). This function is normally used in discussions of bias, and is generally combined with the linear bias parameter, $`b`$, to relate the divergence of the velocity field to the overdensity of the galaxy field at the present day, via the parameter, $`\beta f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)/b`$, where $`\delta _g=b\delta `$ for a linear biasing model.

In this section, we will be treating only unbiased fields, and have introduced $`\beta `$ as a means of comparing this discussion of redshift space distortions to the standard approach (e.g. Strauss & Willick 1995). In §3, we’ll return to the degeneracy in $`\beta `$ and show how PLA might be used to break it.

### 1.3 The Distant Observer Approximation

Up to this point, we’ve treated linear redshift space distortions with more or less full generality. However, since our ultimate goal is to apply these distortions in the context of PLA, we will want to make some simplifying assumptions. For example, PLA (GS) uses a Particle Mesh (PM) Poisson solver (Hockney & Eastwood 1981). This method takes advantage of Fast Fourier Transforms (FFTs), which assume Cartesian coordinates.

The general form of the comoving redshift coordinate, $`𝐬`$, above (equation 3), is not separable in Cartesian coordinates. If this form were to be applied in general, one would wish to describe coordinates with spherical harmonics (e.g. Susperregi 2000). The primary purpose of the current discussion, however, is to examine the underlying dynamics, and while treating redshift space distortions of nearby systems is undoubtedly of cosmological interest, the matter at hand is greatly simplified by assuming the distant observer approximation (d.o.a.).

In the d.o.a. we essentially assume that the system of interest is sufficiently far away that the $`\widehat{𝐱}_i`$ are parallel for all particles. For convenience, we will label the third orthonormal coordinate (the z-axis), as the line of sight. Using this definition, we redefine the comoving redshift coordinate:

$$s^\alpha =x^\alpha +\delta _{\alpha 3}^K\frac{a\dot{x}^\alpha }{\dot{a}},$$
(9)
where $`\alpha `$ denotes the index of the direction vector, and $`\delta _{\alpha \beta }^K`$ is the Kronecker-$`\delta `$ function.

In the linear regime, this becomes

$$s^\alpha (t)=q^\alpha +D(t)p^\alpha [1+\delta _{\alpha 3}^Kf(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t)].$$
(10)
Thus, if a particle is observed at $`𝐬_i^F`$ and its initial position, $`𝐪_i`$, and the cosmology are known, this expression may be inverted and combined with equation (5) to give:

$$x^\alpha (t)=q^\alpha +D(t)p^\alpha =q^\alpha +\frac{D(t)(s^{\alpha F}q^\alpha )}{[1+\delta _{\alpha 3}^Kf(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)]}.$$
(11)
Equations (10,11) may be combined to extend linear theory into redshift space. One must keep in mind that we have assumed that we know both the initial and final position of a particle in the comoving redshift space coordinate. As pointed out above, without both constraints, inverting the redshift coordinate becomes ill-posed.

## 2 Method: Least Action in Redshift Space

Even if we have the idealized set of observations discussed above, and have a complete, unbiased mapping of the density field in redshift space, as perturbations become large, complications will arise both in mapping the redshift space field to a real space one, and in reconstructing an initial density field. A number of researchers have attempted worked on the dual problems of reconstructing an underlying real space CDM density field from a velocity field, and the calculation of a velocity field field from a redshift survey.

For the former problem, the POTENT algorithm (Bertschinger & Dekel 1989; Dekel, Bertschinger & Faber 1990; Dekel et al. 1999) uses the Zel’dovich approximation to relate a redshift/distance survey, in which one component of the peculiar velocity can be directly computed, to an underlying mass density field. Nusser et al. (1991) uses a nonlinear generalization to extend this into the nonlinear regime. Others (Kudlicki et al. 1999; Chodorowski et al. 1998; Chodorowski & Łokas 1997; Bernardeau 1992) use higher order perturbation theory to compute the relationship between the velocity divergence and real space density field. In particular, Chodorowski & Łokas (1997) point out that application of these methods may be used to break the degeneracy between the linear bias constant, $`b`$, and cosmology.

A related problem concerns the calculation of the peculiar velocity field from a galaxy redshift survey. Nusser & Davis (1994) use a quasi-linear correction to the Zel’dovich approximation to relate the redshift space density field to the peculiar velocity field. Like the methods listed above which use perturbation theory to related the velocity and density field, Chodorowski (2000) uses 2nd- and 3rd- order perturbation theory to expand and compare the real space and redshift space density fields.

Others have taken a slightly different approach, which attempts to essentially solve these two problems simultaneously. Giavalisco et al. (1993) suggested that the Least Action approach described by Peebles (1989) could be used in redshift space with only a canonical transform of the coordinates. Schmoldt and Saha (1998) consider the difficulties of running least action reconstruction in redshift space, and test this by reconstructing the velocity field of the Local Group.

Susperregi & Binney (1994) modified this approach and described a technique whereby one could expand the density and velocity fields in Fourier space, and write down the Least Action equations in Eulerian form. Susperregi (2000) took this a step further, and applied a similar code (using Spherical Harmonic transforms) to the IRAS 1.2 Jy redshift survey (Fisher et al. 1995). Each of these techniques use a smoothing filter on the density field such that they are in the mildly nonlinear regime at the present.

PLA has several distinct advantages over these approaches. First, since they are inherently Eulerian, as perturbations becomes large, they no longer fairly sample the matter field. PLA, on the other hand, is Lagrangian in the sense that it performs the time integral over the particle field, rather than the density field. Moreover, the Eulerian PLA approaches assume a locally curl-free velocity field at all times, by the velocity-density relationship. While PLA generally assumes curl-free initial conditions, vorticity is permitted to develop.

Ultimately, we want to reconstruct the underlying real space CDM density field and evolution from a set of observed galaxy redshifts. The approach taken in this paper differs from those discussed in GS in that we now deal with quasi-linear structure, rather than the highly nonlinear constraints. While previously, we were content to get a “realistic” set of initial conditions, here, our aim is reconstruct the details of the observations exactly. By doing this, we hope to disentangle degeneracies in bias, get a handle on the true cosmological power spectrum of perturbations, and examine the growth of large scale structure.

In this section, we describe applying the PLA approach to redshift space constraints. In particular, we will deal with two main issues:

1. Given some observed redshift space density field, $`\delta ^F(𝐬)`$, find final particle constraints, $`\{𝐬_i^F\}`$, which satisfy the density field, which maps to a corresponding initial uniform field, $`\{𝐪_i\}`$, such that the constraints can be most easily satisfied.
2. Given an initial particle constraints, $`\{𝐪_i\}`$, and final particle constraints in redshift space, $`\{𝐬_i^F\}`$, find the trajectories of particles which self-consistently satisfy the boundary conditions.

### 2.1 Computing Particle Constraints

We must first determine the boundary conditions for the particles’ positions at $`a=0`$ and $`a=1`$. Let’s consider a set of idealized observations, in which a smoothed, complete, and unbiased density field in redshift space, $`\delta ^F(𝐬)`$, has been observed.

In order to determine the final boundary condition, we generate a set of particle redshift positions which yield this observed density field. We begin by assuming by distributing particle positions uniformly distributed on a grid, $`\stackrel{~}{𝐬}_i=𝐪_i`$, where the tilde over the redshift coordinate will be explained shortly.

From here, we iterate in the following way. In each iteration, we take the density field of the current value of the particle positions yielding $`\delta (𝐬)`$. We then transform this density field into the target field by using the Jacobian determinant to create a laminar flow. That is, one can adjust the particle positions in the former grid by using a coordinate transformation:

$$\stackrel{~}{𝐬}_i^{}=\stackrel{~}{𝐬}_i+\psi (\stackrel{~}{𝐬}).$$
(12)
In that case, the density field as measured in the primed frame compared to the unprimed frame will be:

$$\left[1+\delta (𝐬^{})\right]d^3𝐬^{}=\left[1+\delta (𝐬)\right]d^3𝐬.$$
(13)
From the form of the transform, the determinant of the Jacobian is easily computed. For small perturbations:

$$\left[1+\delta (𝐬^{})\right]\left[1+\delta (𝐬)\right]\left(1+^2\psi \right).$$
(14)
Assigning $`\delta (𝐬^{})=\delta ^F(𝐬)`$, and using standard Fourier techniques, one can compute the scalar field, $`\psi (𝐬)`$, and taking the gradient, one can compute the coordinate transform, and therefore, can iteratively produce a particle map which satisfies the target density field.

Since we assume throughout that structure evolves out of an initially uniform density field, the initial constraints on these particles must be the uniform grid, $`𝐪_i`$. This is exactly analogous to mapping from an unperturbed simulation since a uniform particle field run through an N-body will remain uniform, and hence serve as a perfectly legitimate unperturbed simulation.

But which final positions correspond to the values of $`𝐪_i`$? In other words, given any set of particle redshift positions, $`\stackrel{~}{𝐬}_i^F`$ which satisfy the density field, $`\delta ^F(𝐬)`$, what is the “best” permutation matrix, $`𝐌_{ij}`$, such that

$$𝐬_i^F=\underset{j}{}M_{ij}\stackrel{~}{𝐬}_j^F\mathrm{?}$$
(15)
We define $`𝐌_{ij}`$ as an $`N_p\times N_p`$ matrix which has exactly one “1” in each row and column, and zeros elsewhere.

We need to define what we mean by “best”. In general, we mean that we wish to compute boundary constraints which most naturally provide us with physically well-motivated orbits. Since the laminar flow method of generating final constraints necessarily assumes no shell crossings, for small perturbations, the permutation matrix will simply be the identity matrix. In this ultra-linear case, in which perturbations are assumed to be so small that there are no orbit crossings even in redshift space, we estimate the physical particle displacement from the uniform field as:

$$p_i^\alpha =\frac{\stackrel{~}{s}_i^{\alpha F}q_i^\alpha }{1+f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t)\delta _{\alpha 3}^K},$$
(16)
The final physical position of particle, $`i`$, is thus related by equation (5). Taking the smoothed density field of $`𝐱_i`$ gives an estimate of the density field in real space.

However, in many cases of interest, even where there are no shell crossings in real space, there are orbit crossings in redshift space. These are the famous “triple value zones,” (Strauss & Willick 1995 §5.9) so named because a particle at a given redshift is a triply degenerate function of distance. Even if a system is dynamically in the linear regime and the real space density field is known with great accuracy, this degeneracy can arise.

In order to illustrate this, we have run a simulation in which a plane wave density field is laid down along the line of sight. Its amplitude is such that there are triple-valued zones at $`a=1`$. The solid lines in Figure 1 show the initial and final density field of this distribution in both real and redshift space. The solid line in Figure 2 demonstrates the existence of triple value zones. There, we plot the relationship between real and redshift particle coordinates for our simulated Zel’dovich pancake.

Given the “observed” redshift space density field, $`\delta ^F(𝐬)`$ in panel d) of Fig. 1, we have already described how to generate a set of redshift space coordinates $`𝐬_i^F`$. Since the simulation is approximately in the Zel’dovich regime, it is assumed that if we can estimate the real space coordinates, we can use that information to estimate the position of a particle on the uniform grid,$`𝐪_i`$. However, given the degeneracy of the mapping from redshift to real space, this is no simple task.

Willick et al. (1997; also Sigad et al. 1998) suggest a likelihood approach to breaking this degeneracy, called VELMOD. Part of VELMOD relates the observed redshift space density field to a test value of the real space density field. We use a similar approach here. We first assume that the real space particle field is the one generated from the assumption of no orbit crossings. That is, equations (16) is used to approximate the positions of the particles. One may then compute a real space density field from the particle field approximations, and from that the potential field, $`\varphi (𝐱)`$ can then be computed. The assumption of linearity gives the following relation:

$$𝐩(𝐱)=\frac{\varphi (𝐱)}{(a_0\ddot{D}_0+\dot{a}_0\dot{D}_0)}.$$
(17)
This makes $`\stackrel{~}{𝐬}^F`$ a straightforward function of $`𝐱^F`$ via equation (5). For triply valued $`\stackrel{~}{𝐬}_i^F`$, we can then determine a posterior probability that the particle is at real space position, $`𝐱^F`$:

$$P(𝐱_i^F|\stackrel{~}{𝐬}_i^F)[1+\delta ^{OLD}(𝐱_i)]\delta ^D\left[𝐱_i𝐱(\stackrel{~}{𝐬}_i^F)\right]\frac{x(𝐬_i)}{\stackrel{~}{s}_i}$$
(18)
where $`\delta ^{OLD}(𝐱)`$ is the previous iteration of the estimated real space density field. This is simply the discrete form of the continuous distribution function used in the VELMOD approach. Using this distribution function, a real space coordinate is randomly assigned to each triple-valued particle. The real space density field and potential are then recomputed and the process is repeated until satisfactory convergence is reached.

We apply this VELMOD-like approach to our Zel’dovich pancake. The dotted and short dashed lines in the lower-left panel of Figure 1 show the initial guess and final estimate of the real space density field. While the fit between the true and estimated real space density fields are good, they are not perfect. One way of thinking about this is that equation (17) assumes that the final velocity of a particle is linearly proportional to the force on that particle. In the Zel’dovich approximation, this is a good assumption. However, as the limits of that approximation are approached, the relation between acceleration and velocity may become more complex. On scales on which shell crossings occur, the two may even be of opposite signs. One approach that people have historically used is to simply smooth density fields until all structures are linear. It should be noted this density field is not our final estimate of the real-space field. Rather, we are using it as a working model to set up redshift space constraints for PLA, which makes no assumptions about the linearity of orbits.

This artificial steepening is also apparent in the estimated relationship between real and redshift space coordinates, as is illustrated by the open square points in Figure 2. Notice, however, that we have qualitatively reproduced the features of the tripled value zones.

At the final iteration of VELMOD, not only are assumed final positions of the particles computed, but so are their velocity/displacement vectors, $`𝐩_i`$. If linear theory approximately holds, then each of those particles ought to have originated at:

$$\stackrel{~}{q}_i^\alpha \stackrel{~}{s}_i^{F\alpha }[1+\delta _{\alpha 3}^Kf(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)]p_i^\alpha .$$
(19)
We may now return to the problem posed at the beginning of this section: How do we compute the “best” permutation matrix, $`𝐌_{ij}`$? We find the matrix which minimizes:

$$\chi ^2[𝐌_{ij}\stackrel{~}{𝐪}_j,𝐪_i]\underset{i}{}\left[𝐪_i𝐌_{ij}\stackrel{~}{𝐪}_j\right]^2.$$
(20)
In order to actually perform this minimization, we apply a simulated annealing method (Press et al. 1992) Even for systems with a number of triple-value zones, the first guess of final and initial particle position pairings produces rapid convergence. One can then permute $`\stackrel{~}{𝐬}_i^F`$ into $`𝐬_i^F`$, resulting in a well-motivated set of boundary constraints.

### 2.2 Computing Trajectories in Redshift Space

We will now consider the simultaneous determination of the orbits of interacting particles when the boundary constraints are given in redshift space. Let us say that we have an isolated, uniform grid of particles at $`a=0`$, with positions given by $`\{𝐪_i\}`$, and at $`t=t_0`$, those particles are “observed” at redshift coordinates, $`\{𝐬_i^F\}`$. We write down the trajectories of the particles as the sum of a part given by linear perturbation theory, and coefficients times basis functions. However, unlike the discussion of PLA in real space (GS), redshift space gives us heterogeneous constraints on our basis functions. In real space, the basis functions were constrained such that:

$$f_n(t_0)=0;\underset{a0}{lim}a^2\dot{f}_n(t)=0$$
(21)
In GS, we showed that these constraints can be satisfied by using

$$f_n=\underset{m=n}{\overset{m_{max}}{}}b_{nm}D(t)^m\left[D_0D(t)\right],$$
(22)
where $`b_{1m}=\delta _{1m}^k`$, and the higher order coefficients are based on fitting to an N-body simulation.

In redshift space, however, we need to define a slightly different set of basis functions, $`\stackrel{~}{f}_n(t)`$. In this case the basis functions along the line of sight must satisfy

$$\stackrel{~}{f}_n(t_0)+\frac{a_0\dot{\stackrel{~}{f}}_n(t_0)}{\dot{a}_0}=0;\underset{a0}{lim}a^2\dot{\stackrel{~}{f}}_n(t)=0,$$
(23)
in order that the varying the coefficients of $`\stackrel{~}{f}_n`$ do not change the corresponding radial redshift space coordinate.

To satisfy these constraints, we introduce a complementary set of basis functions to those introduced in real space:

$$\stackrel{~}{f}_n(t)=f_n(t)\frac{\dot{f}_n(t_0)}{\dot{a_0}}\frac{a^n}{1+n}.$$
(24)
It can be shown that these functions satisfy equation (23) if the real space basis functions described above are used. As with the unaccented basis functions, only the first function goes linearly or slower at early times. Using these basis functions, we are able to describe the trajectory of any particle as:

$`x_i^\alpha (t)=x_i^{(0)\alpha }+D(t)(s_i^{\alpha F}x_i^{(0)\alpha })+{\displaystyle \underset{n}{}}C_{in}^\alpha f_n(t)`$ ; $`\alpha =1,2`$ (25)
$`x_i^\alpha (t)=x_i^{(0)\alpha }+{\displaystyle \frac{D(t)(s_i^{\alpha F}x_i^{(0)\alpha })}{[1+f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)]}}+{\displaystyle \underset{n}{}}C_{in}^\alpha \stackrel{~}{f}_n(t)`$ ; $`\alpha =3,`$ (26)
where $`x_i^{(0)\alpha }`$ are the set of some physically self-consistent orbits, as output from an N-body code. Of course, we may also set $`x_i^{(0)\alpha }=q_i`$.

In GS, we showed that by simultaneously minimizing the action, $`S`$, for each coefficient, $`C_{i,n}^\alpha `$, then the equations of motion of the particles are necessarily satisfied, and hence, we may find these orbits. That is,

$$\frac{S}{C_{i,n}^\alpha }=_0^{t_0}𝑑t\left[f_n(1\delta _{\alpha 3}^K)+\stackrel{~}{f}_n\delta _{\alpha 3}^K\right]\left[\frac{}{t}(a^2\dot{x}_i^{(1)\alpha })+\frac{\varphi _i^{(0)}}{x_i^\alpha }\frac{\varphi _i}{x_i^\alpha }\right],$$
(27)
We will solve equation (27) for all possible basis functions simultaneously, by determining the coefficients such that the kernel

$$𝐠_i(t)\frac{(a^2\dot{𝐱}_i^{(1)})}{t}+\left(\varphi _i^{(0)}\varphi _i\right)$$
(28)
vanishes at all times, and for all particles.

In order to do this, we can use the metric:

$$X^2\underset{i}{}𝑑t|W(t)𝐠_i(t)|^2,$$
(29)
where $`W(t)`$ is an arbitrary weighting function. By minimizing $`X^2`$, we find the set of trajectories which come closest to satisfying the equations of motion.

Perturbations of the basis functions may be approximated by:

$$\delta C_{i,n}^\gamma \left[𝑑t\mathrm{\Delta }_{i,m,\beta }^\alpha (t)\mathrm{\Delta }_{i,n,\gamma }^\alpha (t)\right]_{m\gamma ,n\beta }^1\left[𝑑tW(t)^2g_i^{\alpha (old)}(t)\mathrm{\Delta }_{i,m,\gamma }^\alpha (t)\right]_{n\beta },$$
(30)
where

$$\mathrm{\Delta }_{i,n,\beta }^\alpha (t)\frac{g_i^\alpha (t)}{C_{i,n}^\beta }=\delta _{\alpha \beta }\frac{[a^2x_i^{\alpha (1)}f_n(t)]}{t}\frac{^2\varphi _i}{x_i^\alpha x_i^\beta }f_n(t).$$
(31)
Where $`\beta =3`$, the basis functions, $`f_n(t)`$ should be replaced by $`\stackrel{~}{f}_n(t)`$.

We illustrate the results of this method on the Zel’dovich pancake in Figure 1 and Figure 2. For this, we have used 4 basis functions and two iterations. For each iteration, we took the constraint pairs, $`\{𝐪_i,𝐬_i^F\}`$, and used PLA to compute the best fit full trajectory. We then evaluated the positions and velocities of the particles at $`a0.01`$, and ran the particles through the PM code again.

The long dashed lines in panels c) and d) of Figure 1 represent the density field of the second iteration in redshift and real space, respectively. While by the nature of the constraints we would necessarily expect the redshift space density field to converge to the “true” field, we have no such guarantee in real space. Nevertheless, the real space density field does seem to give a somewhat better fit than the initial estimate given by the VELMOD-like approach, especially around the edge of the peak, PLA gives a smoother edge.

Additionally, even though the peak is nonlinear in the sense that $`\delta >1`$ in real space, PLA is able to very successfully generate an initial density field. Figure 1a) and b) show the redshift and real space density field determined by PLA, as well as the true initial density field. PLA is able to determine the amplitude of the initial peak to about $`10\%`$. However, it should be noted that PLA may erroneously generate too much small scale power. Given some a priori knowledge of the power spectrum, however, one may use a power-preserving filter like the one described in GS in order to appropriately smooth the field. In this 1-d case, in which we only expect a single mode, using such a filter would be gratuitous.

Finally, in Figure 2, we show the relationship between real and redshift space coordinates as determined from the output of the PLA code. Note that the PLA points turn over more smoothly than our initial guess points. This is due to the fact that PLA assumes the field to be evolving from an initially uniform particle field, while the VELMOD-like approach makes no such assumption.

## 3 A High Resolution Test of the Code

To illustrate PLA’s success as a reconstruction scheme somewhat into the nonlinear regime, we have run a high resolution simulation, with $`N_p=128^2`$, $`N_g=256^3`$, and with $`L=800h^1`$ Mpc. The cosmology used is $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`\sigma _8=0.7`$. We then apply PLA under the assumption that the cosmology was known to reconstruct the field.

In Figure 3 we show a density contour of $`\delta =0.7`$ for the observed redshift space density field, smoothed with a Gaussian filter with a radius of $`8h^1`$ Mpc, and a similar plot, but for the field resulting from our PLA reconstruction. Recall that the underlying particle field for the latter is based on running the reconstructed initial conditions through a PM code and taking the smoothed density field. A visual inspection demonstrates that the two fields are virtually identical.

We plot a similar comparison in Figure 4, except with a smoothing radius of only $`4h^1`$ Mpc, and a density contour of $`\delta =2.1`$. On these scales, too, the fields are very similar. The reconstructed field, however, seems to differ somewhat on small scales. If we had an a priori model of the power spectrum, this small scale power could be suppressed numerically.

More quantitatively, Figure 5 shows the fit between the reconstructed and true field as a function of scale. Narayanan & Croft (1999) provide the goodness of fit metric:

$$\mathrm{\Delta }^2(k,t)=\frac{[\delta _1(𝐤)\delta _2(𝐤)]^2}{[\delta _1(𝐤)^2+\delta _2(𝐤)]^2}.$$
(32)
In this case, $`\delta _1(𝐤)`$ and $`\delta _2(𝐤)`$ are the Fourier transforms of the true and reconstructed redshift density fields. For perfect matching on a particular scale, this metric goes to zero. For uncorrelated fields, it goes to one.

As Figure 5 illustrates, the fit for all four comparisons is very good even into the nonlinear regime ($`k0.4h\mathrm{Mpc}^1`$). The best fit was found for the final conditions in redshift space, since this was the actual set of observations to be matched. However, it is shown that the real space field at z=0 is also reconstructed quite well, as are the corresponding initial conditions. In particular, the relevant scale is that for which $`\mathrm{\Delta }^2=0.5`$ for each comparison. In this case, the real and redshift initial conditions are well matched down to a scale of $`18.6`$ and $`19.5h^1`$ Mpc, respectively, and the real and redshift space final conditions are matched down to scales of $`15.1`$ and $`12.7h^1`$ Mpc, respectively.

Another measure of the quality of the reconstruction is the isotropy of the reconstructed real space density field. We may do this by decomposing the power spectrum of the field into Legendre Polynomials (e.g. Hamilton, 1992):

$$P(𝐤)=\underset{l\mathrm{even}}{}𝒫_l(\mu )P_l(k),$$
(33)
where $`\mu `$ is the cosine of the angle between the direction vector, $`𝐤`$, and the line of sight, $`𝒫_l(\mu )`$ are the Legendre polynomials, $`P_l(k)`$ are the azimuthally-averaged multipole expansion of the power spectrum. Taking the inverse Legendre transform, we find (Cole, Fisher and Weinberg 1995):

$$P_l(k)=\frac{2l+1}{4\pi }_1^1𝑑\mu _0^{2\pi }𝑑\varphi P(𝐤)𝒫_l(\mu ).$$
(34)
Here, we will only be using the monopole and quadrupole moments, which, as a reminder are $`𝒫_0(\mu )=1`$ and $`𝒫_2(\mu )=(3\mu ^21)/2`$. We may thus measure the isotropy of the distribution by computing the ratio:

$$Q_x(k)\frac{P_2(k)}{P_0(k)}.$$
(35)
Since the real universe is isotropic, an accurately computed reconstruction should have a quadrupole equal to zero on all scales.

Figure 6 shows the quadrupole ratio as a function of scale for both 1 and 2 iterations. Two things are clear from this plot, however. First, the overall isotropy does improve with additional iterations. Secondly, there is a systematic effect in generating these anisotropies which is almost certainly caused by an anisotropic noise term. In the appendix, we use the quadrupole and hexadecipole moments to illustrate that this effect is dominated by noise, rather than by a systematic underestimate of the radial velocity term, for example.

This noise term comes out of the reconstruction scheme, itself. In the previous section, we discussed how one goes about approximately rewinding the trajectory of a particle in order to determine its initial constraint. However, there was an assumption that linearity approximately held. As structure gets more and more nonlinear, this assumption will fail to hold, and the particle matching technique will break down. Future advances in reconstruction methods will have to take this into account. It may be possible to apply PLA in an iterative and statistical way in order better do this matching.

## 4 Applications: Breaking the Bias Degeneracy

### 4.1 Motivation

Though the reconstruction of nonlinear fields is interesting in its own right, we have begin this investigation into redshift space for the purpose of finding out something about cosmology. As a motivation for this sort of reconstruction analysis of redshift surveys, we address the $`\mathrm{\Omega }_M`$-$`b`$ degeneracy in redshift surveys, and discuss how the degeneracy might be broken without recourse to outside dynamical estimates of $`\mathrm{\Omega }_M`$. We will test the effectiveness of using PLA to break the degeneracy.

We begin by introducing the problem as it appears in the linear regime. Excellent recent reviews of this topic is given by Hamilton (1998) and by Strauss & Willick (1995), and we will therefore present only an overview of the linear biasing problem.

Let us first suppose that we have the full velocity field information about a group of almost uniformly distributed particles. Let us further suppose that these particles are biased with respect to some true underlying CDM field, such that:

$$\delta _g(𝐱,t_0)=b\delta (𝐱,t_0),$$
(36)
where the unsubscripted $`\delta `$ is the CDM density field, and $`\delta _g`$ (for galaxies) represents the density field of some biased tracer of the mass. Note that a straight linear biasing model may be replaced with any deterministic, local, and monotonically increasing function of $`\delta (𝐱,t_0)`$ without changing the essence of the discussion or the biasing problem in general. Blanton (1999) provides an excellent review of various types of biasing models.

With linear bias greater than unity the velocity field is of lower amplitude than one would directly infer from a measurement of $`\delta _g`$. From (36), equation (6) may be replaced by:

$$𝐯(𝐱,t_0)=\beta \delta _g(𝐱,t_0),$$
(37)
where

$$\beta \frac{f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)}{b}\frac{\mathrm{\Omega }_M^{0.6}}{b}$$
(38)
and

$$𝐯=f(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda },t_0)𝐩.$$
(39)
Thus, the divergence of the velocity field and density field are related by the same proportionality constants for all combinations of $`\mathrm{\Omega }_M`$ and $`b`$ which yield the same $`\beta `$. This is the crux of the degeneracy problem.

In a redshift survey, we do not actually know the divergence of the velocity field, but rather must infer it through anisotropies in the redshift density field. Kaiser (1987) shows that in the linear regime, there is a straightforward relationship between the real and redshift space density fields in k-space in the d.o.a.:

$$\delta _g^s(𝐤)=\widehat{𝐒}\delta ^x=\delta ^x(𝐤)[1+\beta \mu ^2],$$
(40)
where $`\widehat{𝐒}`$ is the linear redshift distortion operator, $`\delta _g^s(𝐤)`$ is the Fourier transform of the galaxy density field in redshift space, and $`\delta ^x(𝐤)`$ is the Fourier transform of the underlying CDM density field in real space.

From this relationship, a redshift space power spectrum field may be compared to the real space power spectrum:

$$P_g^s(𝐤)|\delta _g^s(𝐤)|^2=P(𝐤)[1+\beta \mu ^2]^2.$$
(41)
We can further decompose the power spectrum field into Legendre polynomials (equation 33; Hamilton 1992). We then take the inverse Legendre transform (equation 34; Cole, Fisher and Weinberg 1995) and compute the quadrupole moment.

Taking the Legendre expansion of the form of the linear density field (equation 33) above, we find that the moments of the biased, redshift space distribution may be related to the underlying real space distribution as follows (Hamilton 1998):

$`P_{g,0}^s(k)`$ $`=`$ $`\left(1+{\displaystyle \frac{2}{3}}\beta +{\displaystyle \frac{1}{5}}\beta ^2\right)P(k)`$ (42)
$`P_{g,2}^s(k)`$ $`=`$ $`\left({\displaystyle \frac{4}{3}}\beta +{\displaystyle \frac{4}{7}}\beta ^2\right)P(k).`$ (43)
Though the underlying power spectrum is not known, this quadrupole ratio may be computed as:

$$Q_s(k)\frac{P_{g,2}^s(k)}{P_{g,0}^s(k)}=\frac{\frac{4}{3}\beta +\frac{4}{7}\beta ^2}{1+\frac{2}{3}\beta +\frac{1}{5}\beta ^2}.$$
(44)
Since the multipole expansion may be computed directly from the observed redshift space density field, it is clear that in the linear regime, the quadrupole ratio is degenerate for a particular value of $`\beta `$.

In addition, $`\beta `$ may be estimated through the use of distance-velocity comparisons. In practice, measurements of $`\beta `$ are still quite difficult due to the noisy data involved. Willick (2000) gives a review of current estimates, and finds a value of $`\beta =0.5\pm 0.04`$ from the IRAS velocity field and the Mark III Catalog used to estimate distances. This is relatively unchanged from the earlier estimate based on measured redshift-space anisotropy by Cole, Fisher, & Weinberg (1995) of $`\beta =0.54\pm 0.3`$. Ballinger et al. (2000) estimate $`\beta =0.4\pm 0.1`$ for the IRAS 0.6 Jy PSCz survey from analysis of anisotropies. Similar values are found from the Optical Redshift Survey catalog (Baker et al. 1998).

After estimating $`\beta `$, how does one estimate $`\mathrm{\Omega }_M`$ without recourse to direct mass estimates? We answer this by first pointing out that a correctly reconstructed real space density field will be perfectly isotropic. That is, $`Q_x(k)=0`$ on all scales. Nonlinearities is the evolution of the density field will mean that the linear redshift distortion operator will no longer be valid on all scales. By reconstructing fields using PLA for different combinations of $`\mathrm{\Omega }_M`$ and $`b`$, and measuring the quadrupole moment for the reconstructed field, we can determine the “true” cosmology as that which minimizes the anisotropy.

### 4.2 Simulations

#### 4.2.1 Constraining Bias

In order to test this approach, we have run two simulations, each extending only into the mildly nonlinear regime. We have found that simulations which contain highly nonlinear structure do not effectively differentiate between different cosmologies due to the excessive noise and difficulty of doing the particle matching on small scales.

The two simulations were each run with $`\beta =0.5`$ and a boxsize of $`L=1000h^1`$ Mpc, the first with $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`b=1`$, and the second with $`\mathrm{\Omega }_M=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, and $`b=2`$. Each simulation was run with $`64^3`$ particles, and $`128^3`$ gridcells, and the PLA reconstruction used 4 basis functions. The “observations” of these simulations were the redshift-space density fields in the d.o.a. After reconstructing the initial density field for a particular assumed bias, we ran the initial conditions through a PM code, and measured the quadrupole moment in real-space, where it should vanish. The “mean” quadrupole moment is estimated as:

$$Q_x^2=\frac{𝑑kQ_x^2(k)k^2e^{k/k_l}}{𝑑kk^2e^{k/k_l}},$$
(45)
where $`k_l`$ is a limiting scale at which point grid and/or nonlinear effects will become important. We have assumed $`k_l=20\times 2\pi /L_{box}`$, but found similar results for $`k_l=\mathrm{}`$. The peak of this contribution occurs around $`3k_l`$, or on a physical scale of $`16h^1`$ Mpc, the nonlinear scale. The results of each model tested are shown in Figure 7 for the $`\mathrm{\Omega }_M=0.3`$ simulation, and Figure 8 for the $`\mathrm{\Omega }_M=1`$ simulation.

In each case, we find that the best fit value of $`\mathrm{\Omega }_M`$ corresponds to the actual value of $`\mathrm{\Omega }_M`$ used in the simulation. That is, by tracing the detailed evolution of a mildly nonlinear field, PLA can effectively break the $`\beta `$ degeneracy.

In estimating this effect from a real survey, we would have to Monte Carlo observations based on the survey geometry, selection function, and the like, in order to estimate $`\mathrm{\Omega }_M`$ and its errors. While direct error estimation is difficult with only two realizations, the results are quite suggestive that this will be an effective way to constrain $`\mathrm{\Omega }_M`$ directly. Susperregi (2000) uses an Eulerian least action code to similarly show that the bias degeneracy may be broken through accurate reconstruction.

As a final test of PLA as a redshift space reconstruction scheme, we show that one may obtain somewhat better estimates of $`\beta `$, itself, from the assumed isotropy of the reconstructed real space density field than from the redshift space anisotropy. This is illustrated in Figure 9, in which we show a comparison between $`Q_x(k)`$ of the reconstructed field for simulation 1, assuming the correct value of $`\mathrm{\Omega }_M=0.3`$, and the corresponding quadrupole moment residuals, $`Q_s(k)Q_s^L`$ of the observed redshift space field. Note that the term $`Q_s^L`$ is simply that which one would estimate from an assumption of the correct value of $`\beta `$.

On large scales, these two statistics are almost identical. However, on smaller scales, when nonlinearities begin to become important, the two statistics both exhibit anisotropies. The redshift space field, however, becomes anisotropic on larger scales. Since the corresponding uncertainties in $`\beta `$ are approximately proportional to one over the square root of the number of modes probed, we may relate the expected uncertainties from the reconstructed field to that from the redshift field as:

$$\frac{\sigma _\beta ^{PLA}}{\sigma _\beta ^Z}\frac{\sigma _Q^{PLA}\frac{\beta }{Q^{PLA}}(k_{\mathrm{max}}^{PLA})^{3/2}}{\sigma _Q^Z\frac{\beta }{Q^Z}(k_{\mathrm{max}}^Z)^{3/2}},$$
(46)
where the superscript “Z” refers to the estimate from the redshift space field and the superscript “PLA” refers to the estimate from the reconstructed real space field. For $`\mathrm{\Omega }_M=0.3`$, all of the partial derivatives are almost exactly one. Moreover, estimates of the scatter in the quadrupole estimates show that $`\sigma _\beta ^{PLA}0.9\sigma _\beta ^Z`$. Finally, since Figure 9 shows $`k_{\mathrm{max}}^{PLA}1.25k_{\mathrm{max}}^Z`$, the approximate relation between the uncertainty in the bias between the two errors is $`\sigma _\beta ^{PLA}0.65\sigma _\beta ^Z`$. Thus, we are able to generate a somewhat better constrained estimate of $`\beta `$ from the reconstructed density field, as compared to the observed field.

## 5 Future Prospects

This paper has discussed the problem of reconstructing the underlying real space CDM density field and its evolution from galaxy redshift surveys under rather idealized conditions. We have assumed full sampling, a constant linear bias relation with no morphological segregation, the distant observer approximation, no errors in measurement, and a very regular geometry. We have shown that theoretically even a mildly nonlinear field can be reconstructed using PLA to break the bias degeneracy. In actual observations almost none of these assumptions will hold. We would like to end this work with a brief discussion of how these effects might be appropriately modeled, and mention a few possible candidates of real surveys to which the PLA method might be applied.

Throughout, we have assumed a cubic geometry. In simulating a realistic survey, a mask must be applied such that statistics may be correctly computed for the true survey volume. Additionally, for many observational samples of interest, the distant observer approximation no longer describes the system adequately. Over the course of its lifetime, a galaxy may have traversed a significant angle in the sky, making the use of Cartesian coordinates difficult. An obvious solution to this problem is to write the PLA equations in spherical coordinates.

An additional concern in the application of PLA to realistic observations is that the biasing model that we have used here is a strictly linear one, and observational evidence suggest that biasing may be much more complex (e.g. Blanton 1999 and references therein). In fact, we have only used linear bias in this work for its simplicity. PLA would be applicable to any model of deterministic bias, even one which had an explicit model of morphological segregation. It is not clear how one might effectively reconstruct a field under the assumption of a significantly stochastic bias model.

Even under the simplest of approximations, realistic surveys are not cubic, not necessarily contiguous, and generally have nonuniform selection functions. The issues of contiguity and geometry are related, in that in both we need to approximate the density field outside the survey volume in order to correctly estimate the potential field within.

Lahav et al. (1994) apply one such technique to the IRAS 1.2 Jy survey (Strauss et al. 1992; Fisher et al. 1995) in which expansion of the density field in spherical harmonics is used to reconstruct the field outside the survey volume. In essence, this is very similar to assuming a particular autocorrelation function, and generating an outside field based on a truncated form of that function and the observed field near the edges. In many respects, this is quite similar to the sort of reconstruction done using the “Constrained Initial Conditions” technique by Hoffman & Ribak (1991,1992), since both use observed the observed autocorrelation function to build realistic external fields around observed structure.

The final issue in real surveys concerns non-uniformity and noise within the survey volume. However, for flux limited surveys such as the IRAS 1.2 Jy and PSCz (Saunder et al. 2000) surveys, at large distances, shot noise begins to dominate calculations of the density field. Since the observed density field is derived from an incomplete sampling of a finite number of discrete points, the uncertainties in the observed density field behaves like a Poisson statistic, with $`\sigma 1/\sqrt{N_{gal}}`$. In order to correctly anticipate the effects of shot noise when running simulations of a particular survey, a random component needs to be added to the observed density field.

Moreover, a non-uniform selection function, coupled with observations in redshift space results in Malmquist bias. That is, the “true” selection function is based on observed fluxes, and hence is a function in real space. The observations, however, are of densities in redshift space. Strauss & Willick (1995) give an excellent review of how these issues are dealt with in reconstructing density fields.

These considerations are all with an eye toward learning about cosmology from redshift surveys. For example, the IRAS 0.6 Jy PSCz Survey (Saunders et al. 2000) is an especially promising recent candidate for analysis, as it is publicly available and has nearly full sky coverage. Hamilton, Tegmark, & Padmanabhan (2000) have already estimated $`\beta =0.41_{.12}^{+.13}`$ from this survey, and Ballinger et al. (2000) find a similar result of $`\beta =0.4\pm 0.1`$. Both methods used only linear theory, however. Valentine, Saunders, & Taylor (2000) do a somewhat higher order reconstruction, by using the PIZA method, and find a best fit to the survey with a slightly higher value of $`\beta 0.5`$. In testing our code, we have found that one can get a better measure of $`\beta `$ by using PLA to reconstruct a density field under some fiducial cosmology, and comparing cosmologies to see which produce the minimum anisotropy in the real space field. We have showed that ideally, PLA can be used to discriminate between “degenerate” pairs of bias and $`\mathrm{\Omega }_M`$. Finally, we showed that PLA produces a means by which uncertainties in the measurement of $`\beta `$ itself can be reduced.

Another interesting prospect is the application to PLA to large distance redshift surveys, since one of the byproducts of PLA is the real space density field. Nusser et al. (2000) use the $`D_n\sigma `$ relation of the ENEAR redshift-distance survey (da Costa et al. 2000) as test particles within the PSCz survey, much as we would wish to do using PLA. A reconstruction was done using the method described by Nusser & Davis (1995). Predicted distances from the reconstructed field can then be compared with the estimated distances from the ENEAR survey. They also estimated $`\beta 0.5`$.

While the PSCz survey contains $`15,000`$ redshifts, the current generation of redshift surveys is producing an even greater opportunity to measure statistical and global properties of the universe. When complete, the SDSS redshift survey (York et al. 2000), will produce $`10^6`$ galaxy redshifts, and will cover a quarter of the sky out to Petrosian magnitude, $`r_p^{}=17.7`$ to $`z0.15`$. The 2dF survey (Colless 1999) will ultimately measure redshifts over a quarter of a million galaxies, out to $`b_J=19.5`$.

Reconstruction of fields from these enormous datasets will prove a significant computational challenge. However, it is well worth it, as PLA can yield insight into the underlying power spectrum, bias, and cosmological parameters.

I would like to gratefully acknowledge many helpful comments by David Spergel, Michael Strauss, Vijay Narayanan, and Jim Peebles, as well as a superior visualization tool by Michael Blanton. This work was supported by an NSF Graduate Research Fellowship and NASA ATP grant NAG5-7154.

## Appendix A The Effects of Noise on Multipole Moments

In §3, we found that the reconstructed real space density field of the high-resolution simulation had a significant and seemingly systematic quadrupole moment at small scales. The question which arises out of this is, does this quadrupole represent a systematic error in the reconstruction along the line of sight (such as would occur, for example, if the assumed value of $`\beta `$ were incorrect), or does it represent a stochastic term?

To examine this question, let us consider the following form of a reconstructed field with a very simple error term:

$$\stackrel{~}{\delta }(𝐤)=\delta (𝐤)\left[1+\mu ^2B(𝐤)+\mu ^2ϵ(𝐤)\right],$$
(A1)
where $`𝐤`$ is a Fourier space component, $`\delta (𝐤)`$ is the true real space field we are attempting to reconstruct, $`\stackrel{~}{\delta }(𝐤)`$ is the reconstructed form of the field, $`\mu ^2`$ is the cosine of the angle between the $`𝐤`$ and the line of sight, $`B`$ is a systematic, “bias” term, and $`ϵ(𝐤^\mathrm{𝟐})`$ is an anisotropic random error drawn from a $`N(0,\sigma ^2(𝐤))`$ distribution.

If the corresponding real space density field of redshift space observations have been perfectly reconstructed and contain no noise, the reconstructed field is simply equal to the true underlying field, and thus will be perfectly isotropic in $`k`$-space. However, let us imagine that a field contains no noise, but the reconstruction is such that we have assumed that $`\stackrel{~}{\delta }(𝐤)=\delta _s(𝐤)`$, or the real space density field is the same as the redshift space field. Under those circumstances, $`B(𝐤)=\beta `$, hence our terminology.

Finally, let us consider a more general case, one in which we wish to test for a systematic form of $`B(𝐤)`$ and for the existence of a random noise component. In that case, the reconstructed three-dimensional power spectrum may be written as:

$$\frac{\stackrel{~}{P}(𝐤)}{P(𝐤)}=1+2\mu ^2B(𝐤)^2+\mu ^4B(𝐤)^2+\mu ^4\sigma ^2(𝐤)$$
(A2)
If we then decompose these terms into multipole moments (see §5.4.1) and assume that the noise and systematic terms are simply scale dependent, we find:

$`\stackrel{~}{P}_0(k)`$ $`=`$ $`P(k)\left[1+{\displaystyle \frac{2}{3}}B(k)+{\displaystyle \frac{1}{5}}B(k)^2+{\displaystyle \frac{1}{5}}\sigma ^2(k)\right]`$ (A3)
$`\stackrel{~}{P}_2(k)`$ $`=`$ $`P(k)\left[{\displaystyle \frac{4}{3}}B(k)+{\displaystyle \frac{4}{7}}B(k)^2+{\displaystyle \frac{4}{7}}\sigma ^2(k)\right]`$ (A4)
$`\stackrel{~}{P}_4(k)`$ $`=`$ $`P(k)\left[{\displaystyle \frac{8}{35}}B(k)^2+{\displaystyle \frac{8}{35}}\sigma ^2(k)\right]`$ (A5)
We may then look at the relationship between the quadrupole ratio, $`\stackrel{~}{Q}(k)\stackrel{~}{P}_2(k)/\stackrel{~}{P}_0(k)`$ and the hexadecipole ratio, $`\stackrel{~}{H}(k)\stackrel{~}{P}_4(k)/\stackrel{~}{P}_0(k)`$. We have two extremes. In the case of no anistropic noise component, each ratio is simply a parametric function of $`B(k)`$, and thus, a straightforward relation between the two may be plotted.

If, on the other hand, the anisotropic noise term dominates, we find the relation:

$$\stackrel{~}{H}^{NOISE}(k)=\frac{2}{5}\stackrel{~}{Q}^{NOISE}(k).$$
(A6)
We illustrate this in Figure 10. We show that for an observed redshift density field, the deterministic bias term dominates. Though we do not necessarily have the correct form of the anisotropic error term, this error which would seem to be at the root of the systematic small scale quadrupole in the high-resolution simulation.

Since we may take as a prior that the universe is inherently isotropic, regularization or iterative techniques might be employed in future reconstruction schemes which explicitly find a fully isotropic real space solution.


format: text

text:
# Testing the Hypothesis of the Morphological Transformation from Field Spiral to Cluster S0

## 1 Introduction

The morphologies of galaxies remain one of the key observational characteristics used to classify and differentiate classes of galaxies. It is known from observations of the local Universe that galaxy morphology depends on environment, with the galaxy population within clusters dominated by early-types and spirals being found predominately in lower-density environments. Even within the clusters, the morphological mix changes with local galaxy density (and hence with cluster-centric distance) – producing a morphology-density relation – as shown by Dressler (1980, D80). Moreover, similar segregation behaviour has been seen in Hubble Space Telescope (HST) imaging of ten clusters at higher redshifts, $`z=0.37`$–0.56 by Dressler et al. (1997, D97). These results confirm galaxy morphologies as one of the most sensitive tracers of environmental influences currently known.

One of the most interesting results to come out of the morphological studies of distant clusters has been the claim by Dressler et al. (1997) that the S0-to-elliptical (S0/E) ratio is a strongly decreasing function of redshift: by a factor of $`5`$ between $`z0`$ and 0.5 (although see Andreon 1998; Fabricant, Franx & van Dokkum 2000). This contrasts with the increasing number of blue, spiral galaxies seen in these regions at higher redshifts; the origin of the Butcher-Oemler effect (Butcher & Oemler 1984; Dressler et al. 1994; Couch et al. 1994, 1998; D97). D97 therefore proposed that most of the cluster S0’s are formed relatively recently by the morphological transformation of spirals (see also Poggianti et al. 1999; Fabricant et al. 2000). These spiral galaxies are continuously supplied by accretion from the surrounding field in the course of the assembly of the cluster. Couch et al. (1998, C98) and Fasano et al. (2000) have added estimates of the S0/E ratios for three clusters at $`z=0.31`$ and nine clusters at $`0.09<z<0.26`$, respectively, to bridge the gap in the D97 analysis between $`z=0.37`$ and the present-day. These studies confirm the trend of a decline in the cluster S0 population from $`z=0`$ out to $`z>0.5`$. Proposed mechanisms for transforming the morphologies of bulge-strong spirals to S0 galaxies include dynamical interactions, such as harassment (Moore et al. 1996, 1999a), while gas dynamical processes, e.g. ram-pressure stripping (Abadi, Moore & Bower 1999), will also aid in truncating star formation and hence fading the disk component of any spiral galaxy within the cluster environment.

The luminous galaxy populations of rich clusters at $`z0`$ are dominated by passive galaxies with early-type morphologies and apparently old stellar populations, as shown by their tight colour-magnitude relation (e.g. Bower, Lucey & Ellis 1992; Terlevich et al. 1999). The transformation of star-forming spiral galaxy to passive cluster member must occur, as the blue, spiral galaxies we see in clusters at $`z0.5`$ cannot escape and hence their descendents have to reside in these regions at the present-day (Bower, Kodama & Terlevich 1998; Smail et al. 1998; van Dokkum et al. 1998). Following the work by Smail et al. (1998) and Poggianti et al. (1999) on the relationship between the blue galaxies in distant clusters and nearby cluster S0’s, Kodama & Bower (2001) quantified the analysis by sketching the flow of galaxies across the colour-magnitude diagram down to the present-day, taking into account the galaxy accretion from the surrounding fields.

Motivated by this success in connecting the photometric evolution of the galaxy populations in distant clusters to the present-day counterparts within the framework of the assembly of rich clusters, we now try to seek to explain the evolution of the morphological mix in clusters, and especially the decline in the S0 population, using a simple model for the growth of clusters and the evolution of their galaxy populations.

Since the proportion of spiral to elliptical galaxies in the field is much higher than that in clusters at any time from $`z0.5`$ down to the present-day (e.g. Driver et al. 1998; van den Bergh et al. 2000), continuous accretion of galaxies from the field should be an effective process in changing the morphological mix in clusters. If field spiral galaxies are turned into S0’s upon accretion on to clusters, then the accretion process effectively increases the fraction of S0 galaxies with time, and might explain the rapid increase of S0/E ratio presented in D97. The aim of this Paper is to test this hypothesis quantitatively and thereby constrain both the galaxy accretion rate and the timescale of the morphological transformation. Unless otherwise stated, we use H$`{}_{0}{}^{}=50`$ km s-1 Mpc-1 and q$`{}_{0}{}^{}=0.1`$.

## 2 The Morphological Mix in Distant Clusters and the Field

### 2.1 The Morphological Mix in Clusters

To investigate the morphological evolution of the galaxy population in clusters we use observations of galaxy morphologies in three samples of rich clusters spanning $`z=0`$–0.5. These come from the survey of local clusters by Dressler (1980), where we use the same ten high concentration clusters as in D97, and two different HST–based studies: Couch et al. (1998)’s work on three clusters at $`z=0.31`$ and the morphological catalogue of ten rich clusters at $`z=0.37`$–0.56 from Smail et al. (1997b), which is the basis of the D97 analysis.

To reliably compare the galaxy samples at different redshifts we must apply the same selection criteria to each survey. Specifically, when attempting to trace the evolution of galaxies between high redshifts and the present day, we have to take into account the various degrees of fading they will suffer after their star formation ceases, which depends upon their star formation histories. For example, the most vigorously star forming late-type spirals in the distant cluster samples may fade substantially ($`<1`$ mag) after they cease star formation (see Kodama & Bower 2001) and therefore only the brighter examples ($`M_V<21`$) will appear in a local cluster sample limited at $`M_V20`$. For the simple model discussed here we have chosen to ignore the added complication of a star-formation dependent magnitude cut and simply adopt a fixed magnitude limit for all morphological types. The D80 sample is limited at an absolute magnitude equivalent to $`M_V20.4`$ (D97). We have therefore also chosen to limit the more distant samples from D97 and C98 at an absolute magnitude of $`M_V=20`$ (we note that the magnitude limit quoted in D97 for their analysis is incorrect and actually corresponds to $`M_V=19`$). By adopting a fixed absolute magnitude cut we will tend to over-estimate the proportion of descendents of late-type spirals predicted in the local clusters. However, we note that the fraction of late spirals (Scdm) is small (7 per cent) and therefore our results do not significantly change if we use, for example, a $`M_V=21`$ cut for this population.

In addition, because of the morphology–density relation (D97), it is important to use equivalent areas in the nearby and distant clusters when comparing the galaxy populations. The morphological mix in the high redshift clusters in D97 and C98 is calculated within the HST fields, which roughly correspond to the inner 130′′ or 1.0 Mpc ($`z0.5`$) diameter regions. For the nearby sample from D80 therefore we use the same physical radius to define the morphological sample.

We show in Figs. 1 and 2 the observed S0/E and Sab/S0 ratios, respectively, from D80, C98 and D97. For the local sample of D80 we plot both the individual clusters as well as the average and the scatter for the whole sample. We note that observations of all of the ten clusters are consistent with the mean value given the large uncertainties of the individual cluster measurement arising from sampling errors.

### 2.2 The Morphological Mix in the Field

To successfully model the evolution of the morphological mix in clusters we need to know the distribution of morphologies for the galaxy population they accrete from the field. Ideally the morphological mix for this field sample would be based on a large, volume-limited sample culled from a morphologically-classified faint field galaxy redshift survey. Unfortunately, there are no large redshift surveys with adequate morphological information and so we have instead chosen to use a magnitude-limited morphological survey at a depth which should place the majority of the galaxies within the redshift range we require. We therefore use the morphological classifications from the HST–based Medium Deep Survey (MDS, Griffiths et al. 1994) which reaches $`I=22`$ and will have a median redshift of $`z0.5`$ with 75 per cent of the galaxies at $`z0.65`$, based upon the redshift distribution at this depth in Lilly et al. (1995). This means that the majority of the sample represent the population which will be accreted onto the clusters in the redshift range covered by our analysis. The morphological mix of field galaxies for this sample is then: E, $`0.10\pm 0.01`$; S0, $`0.10\pm 0.01`$; Sab, $`0.25\pm 0.01`$; Scdm, $`0.30\pm 0.01`$; and Irr, $`0.25\pm 0.01`$ (as used by Smail et al. 1997b). The errors are based on Poisson statistics. This distribution is consistent with the smaller sample of van den Bergh et al. (2000) in the Hubble Deep Field with spectroscopic redshifts which place them at $`z=0.25`$–0.6, as well as the small number of field galaxies with morphological and spectroscopic information from Dressler et al. (1999). This then provides us with a morphological distribution for the field galaxies which are accreted onto clusters at $`z<0.5`$.

In a similar manner to that discussed above, when comparing the cluster populations at different epochs, we also need to consider a variable magnitude limit when estimating the properties of the field galaxies which are accreted onto the clusters, depending upon their previous star formation histories (or crudely the morphological type). However, the lack of any large morphologically-classified faint galaxy redshift surveys also means that we do not have the information available to reliably apply a variable apparent magnitude cut to the different morphological subsamples in the field population. To quantify the likely error which will arise from adopting a single magnitude cut across all morphological types we turn to the type-dependent luminosity functions available for local galaxy samples111While there is no suitable high redshift sample to use, we note that there is no evidence for differences in the evolution of the bright end of the luminosity functions for colour-selected (red and blue) galaxy populations out to $`z0.5`$, at least down to our limit of $`M_V=20`$ (Lilly et al. 1995). from the 2dF Galaxy Redshift Survey (Folkes et al. 1999). Correcting the $`B`$-band luminosity functions to $`V`$-band using for type-dependent $`(BV)`$ colours from RC3 (Buta et al. 1994), and adopting limits of $`M_V=20`$ for field E/S0s and $`M_V=20`$ or $`M_V=21`$ for spirals (Sab+Scdm) we can estimate the likely variation in the fractions of E/S0, Sab and Scdm. Using a limit of $`M_V=21`$ for the spiral galaxies, we estimate relative proportions of 1:0.90:0.95 for E:Sab:Scdm, compared to 1:1.25:1.5 for the $`M_V=20`$ cut (we note that this latter ratio is very similar to the value derived above for the distant field, using $`M_V=20`$). We conclude therefore that by adopting a fixed magnitude limit across all morphological types in the field we increase the proportion of spiral galaxies in the field by factor of $`1.4`$–1.5 for both early- and late-type systems. We will discuss the consequence of this effect in §4.

## 3 Modelling the Evolution of Galaxy Morphology in Clusters

Using the morphological distributions for the cluster and field populations derived in the previous section we devise a series of models for the evolution of the morphological mix in clusters from $`z0.5`$ to $`z=0`$. These models are by necessity simple and ignore many of the details of the dynamical processes which may be involved in the morphological transformation, such as galaxy harassment (Moore et al. 1996; 1999a) and/or ram-pressure stripping (Abadi et al. 1999). However, the relative importance of these processes is still a matter of debate (e.g. Poggianti et al. 1999) and the details of the processes are highly uncertain. We therefore adopt a phenomenological approach to test the basic foundations of the morphological transformation hypothesis for the evolution of the cluster S0 population.

We compare the predictions of the various models for the evolution of the S0/E and Sab/S0 ratios within the clusters with the observations from D80, C98 and D97. As we show below, these comparisons constrain both how large a fraction of the cluster population must be accreted since $`z0.5`$ and how long the morphological transformation takes.

There are three key parameters in the models: (1) the total galaxy accretion rate onto the cluster; (2) the range of spiral types (e.g. Sab, Scdm) which transform into cluster S0’s, and (3) the timescale of this morphological transformation. We note that the accretion rate and morphological range have almost the same effect on the morphological evolution in the clusters, since both determine how many of the accreted field galaxies are turned into S0’s. In contrast the transformation timescale controls how long the accreted spirals are classified as having ‘spiral’ morphology after entering the cluster. The parameters for the models we consider are summarised in Table 1.

The cumulative galaxy accretion rate onto a cluster, $`A_{\mathrm{cum}}(t)`$, is normalised to unity at $`z=0`$ (equivalent to $`t=t_{\mathrm{univ}}`$, the age of the present-day Universe — 16.6 Gyr in our cosmology) as:

$$A_{\mathrm{cum}}(t_{\mathrm{univ}})=_0^{t_{\mathrm{univ}}}A(t)𝑑t=1,$$
(1)
where $`A(t)`$ is the normalised accretion rate at a cosmic time $`t`$. To quantify the amount of recent accretion we define $`A_{\mathrm{recent}}(t)`$, the integrated accretion rate from time $`t`$ (or redshift $`z`$) to the present-day ($`t=t_{\mathrm{univ}}`$, or $`z=0`$):

$$A_{\mathrm{recent}}(t)=A_{\mathrm{cum}}(t_{\mathrm{univ}})A_{\mathrm{cum}}(t).$$
(2)
The actual form for the accretion rate is taken from the Extended Press-Schechter theory (Bower 1991), assuming a power spectrum with $`n=1.5`$. Based upon the weak lensing estimates of the masses of the clusters used in D97’s analysis (Smail et al. 1997a) we adopt a mass for the present-day cluster of $`10^{15}`$ M, or 100 M. Since the galaxy accretion history is strongly dependent on the threshold mass for the substructure within the clusters (Bower 1991), we take the fiducial threshold mass of 10, 30 and 50 M for the low, high and super-high accretion rate models, respectively.

For a particular accretion history, the evolution of the morphological mix (fraction of E, S0 and Sab galaxies) in clusters is calculated as follows:

$`f_\mathrm{E}(t)`$ $`=`$ $`{\displaystyle \frac{1}{A_{\mathrm{cum}}(t)}}\{A_{\mathrm{cum}}(t_1)\times f_\mathrm{E}(t_1)+`$ (3)
$`(A_{\mathrm{cum}}(t)A_{\mathrm{cum}}(t_1))\times f_\mathrm{E}^{\mathrm{field}}\},`$
$`f_{\mathrm{S0}}(t)`$ $`=`$ $`{\displaystyle \frac{1}{A_{\mathrm{cum}}(t)}}\{A_{\mathrm{cum}}(t_1)\times f_{\mathrm{S0}}(t_1)+`$ (4)
$`(A_{\mathrm{cum}}(t)A_{\mathrm{cum}}(t_1))\times f_{\mathrm{S0}}^{\mathrm{field}}+`$
$`(A_{\mathrm{cum}}(t\tau _{\mathrm{trans}})A_{\mathrm{cum}}(t_1\tau _{\mathrm{trans}}))`$
$`\times f_{\mathrm{Sp}\mathrm{S0}}^{\mathrm{field}}\},`$
$$f_{\mathrm{Sab}}(t)=\frac{1}{A_{\mathrm{cum}}(t)}\{(A_{\mathrm{cum}}(t)A_{\mathrm{cum}}(t\tau _{\mathrm{trans}}))\times f_{\mathrm{Sp}\mathrm{S0}}^{\mathrm{field}}\},$$
(5)
where $`t_1`$ is the cosmic time at $`z=0.54`$ (9.8 Gyr in our cosmology) the mean redshift of the four most distant clusters in the D97 sample, the epoch at which we normalise the fraction of ellipticals and S0’s ($`f_\mathrm{E}(t_1)=0.40\pm 0.08`$ and $`f_{\mathrm{S0}}(t_1)=0.16\pm 0.04`$, see below). $`\tau _{\mathrm{trans}}`$ is the timescale for the morphological transformation from a spiral to S0 after accretion, and $`f_{\mathrm{Sp}\mathrm{S0}}^{\mathrm{field}}`$ is the fraction of field galaxies that are eventually turned into cluster S0’s. Based on §2.2 this is taken as:

$$f_{\mathrm{Sp}\mathrm{S0}}^{\mathrm{field}}=\{\begin{array}{ccc}f_{\mathrm{Sab}}^{\mathrm{field}}\hfill & =0.25\pm 0.01\hfill & \mathrm{Model}(\mathrm{c}),\hfill \\ f_{\mathrm{Sab}}^{\mathrm{field}}+f_{\mathrm{Scdm}}^{\mathrm{field}}\hfill & =0.55\pm 0.02\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(6)
In the following we also assume that any bulge-strong galaxies observed in the clusters are in the process of morphological transformation into S0’s and we ignore merging of galaxies, since it is unlikely to be important for typical galaxies in these rich cluster environments at $`z<0.5`$, where the relative velocities are high, $`>1000`$ km s-1.

To normalise the model for the evolution of the morphological mix within distant clusters we start from the morphologically-classified galaxy population seen within clusters at the highest redshifts covered by our analysis, $`z0.5`$. For this we use the average of fractions for the four most distant clusters in the Smail et al. (1997b) sample used by D97: Cl 0412$``$65 ($`z=0.51`$), Cl 1601+42 ($`z=0.54`$) Cl 0016+16 ($`z=0.55`$) and Cl 0054$``$27 ($`z=0.56`$). These are: E, $`0.40\pm 0.08`$; S0, $`0.16\pm 0.04`$; Sab, $`0.32\pm 0.10`$; Scdm, $`0.07\pm 0.06`$; and $`0.05\pm 0.03`$ Irr/compact/unclassified. The errors show the scatter between the individual clusters. The Poisson errors are always small; about half of the errors shown. These numbers should represent the morphological mix for a ‘typical’ $`z0.5`$ cluster, since the fractions are simply averaged across the four clusters rather than weighted by richness.

## 4 Results

We plot in Figs. 1 and 2 the model S0/E and Sab/S0 ratios as a function of redshift calculated from our model, normalised to the average cluster population at $`z=0.54`$. Fig. 1 shows that the evolution of the S0/E ratio depends largely on the recent accretion rate, $`A_{\mathrm{recent}}`$, (compare Models (a),(b) and (f)) and/or what proportion of the accreted spiral population are turned into S0’s (Models (a) and (c)), and that it is relatively insensitive to the transformation timescale ,$`\tau _{\mathrm{trans}}`$ (Models (a),(d) and (e)). To match the strong evolution seen in the observations we must have both a high accretion rate, typically $`>50`$ per cent of the cluster population added since $`z0.5`$, and also transform the majority of the accreted spirals into S0s (not only bulge-strong Sab galaxies, but also bulge-weak Scdm’s). We discuss the possibilities for transforming bulge-weak spirals into S0’s in §5. Based upon Fig. 1 we conclude that Models (a), (d), (e) and (f) are preferred. We note that if we applied the type-dependent magnitude cut for the accreted field galaxies to take into account the fading effect discussed above, the necessary accretion rate would need to be even higher since there would be fewer bright spirals available to be added from the field. This would strengthen the conclusions presented here. This result confirms quantitatively the earlier discussion by Fabricant et al. (2000), who suggested that a high accretion rate of field galaxies was needed to explain the strong evolution in the morphological mix in clusters between $`z=0.33`$ and the present-day.

A further constraint on the models comes from the spiral fraction (Fig. 2). As our model assumes that most spiral galaxies observed in the cluster cores are undergoing a morphological transformation into S0’s, the fraction of spiral galaxies in clusters can be used to estimate the timescale of this process (Poggianti et al. 1999). However, the spiral fraction also depends on the accretion rate and/or the proportion of the total spiral population which can transform into S0’s. As the S0 fraction depends upon the same parameters, we can reduce this dependence by normalising the Sab fraction using the S0 fraction. The Sab/S0 ratio thus depends primarily on the morphological transformation timescale (Fig. 2) and can be used to constrain this timescale.

We show in Fig. 2 the predicted Sab/S0 ratios for our models and compare these to the observations. To illustrate the dependence on the transformation timescale we show three models: (d), (a) and (e) for which $`\tau _{\mathrm{trans}}=3`$, 2, 1 Gyrs respectively. As Fig. 2 shows, if $`\tau _{\mathrm{trans}}`$ is too short (1 Gyr in Model (e)), the model struggles to retain enough accreted galaxies with spiral morphologies in the clusters before they are transformed into S0’s. On the other hand, Model (d) ($`\tau _{\mathrm{trans}}=3`$ Gyr) tends to over-produce the proportion of cluster spirals compared to the observations. Equally, adopting a very high accretion rate model, Model (f), which is consistent with the S0/E ratio evolution (Fig. 1), and $`\tau _{\mathrm{trans}}`$ of 1 Gyr tends to underestimate the Sab/S0 ratios (Fig. 2). Taking into account the systematic uncertainty in the normalisation of the models at $`z0.5`$ using the fractions of E and S0 galaxies in the four distant clusters (Fig. 1 and §3), our model suggests that the timescale for the morphological transformation should be in the range $`\tau _{\mathrm{trans}}1`$–3 Gyr.

Note that we have included field Scdm’s into the cluster Sab category in the models, except for Model (c) (Eqs. (5) and (6)), since it is possible that the accreting Scdm’s are transformed into Sab’s either prior to, or immediately after, their entry into the cluster cores (§5). If we did not take into account this ‘pre-transformation’ from bulge-weak spirals to bulge-strong ones, the Sab/S0 ratios in the models would go down by $`50`$ per cent. In this case, only Model (d) would reproduce the observed trend, suggesting $`\tau _{\mathrm{trans}}3`$ Gyr.

The relatively long timescale required for the morphological transformation in the model is consistent with a number of recent results on the properties of galaxies in distant clusters. Based upon the paucity of post-starburst signatures in the spectra of the luminous S0 population in distant clusters and the prevalence of red, passive galaxies with late-type morphologies, Poggianti et al. (1999) concluded that the morphological transformation of S0’s had to occur on a longer timescale than the stellar evolution lifetimes ($`>1`$ Gyrs) of any A-stars formed in the most recent star formation event in these galaxies (either while they were in the field, or any activity which was induced on entering the cluster). Comparisons between the properties of the luminous S0 and elliptical galaxies in distant clusters have been used to place limits on the relative evolution of the two populations: work on three clusters at $`z=0.31`$ – using spectral (Jones, Smail & Couch 2000) and photometric (Couch et al. in preparation) analysis concluded that the S0 and elliptical populations are remarkably homogeneous, and similar conclusions were reached using precise photometric analysis of the colours of S0/ellipticals in the cluster cores at $`z=0.18`$ (Smail et al. 2001), $`z=0.33`$ (van Dokkum et al. 1998) and $`z0.5`$ (Ellis et al. 1997). All of these studies suggest that the luminous S0 galaxies in the cluster cores ceased their star formation at least 2 Gyr prior to the observed epoch (but see the discussion in Smail et al. 2001 for less luminous S0’s). Observations of the evolution of the colour distribution of cluster galaxies (Kodama & Bower 2001) and the radial gradients in star formation, measured using \[Oii\], (Balogh, Navarro & Morris 2000) both suggest that the star formation in a galaxy declines relatively slowly after it has been accreted by a cluster, $`1`$ Gyr.

Taken together these results would indicate that an accreted star-forming spiral galaxy transforms into a passive S0 within $`>2`$ Gyrs of entering the cluster environment. Our analysis confirms this suggestions. This is encouraging as our analysis relies purely on modelling the evolution of the morphological mix within the clusters, whereas the previous results are based on the stellar populations within individual cluster galaxies and so these two approaches are completely independent.

## 5 Discussion and Conclusions

We illustrate a simple phenomenological model to attempt to reproduce the broad properties of the morphological evolution of cluster populations across $`7`$ Gyrs from $`z0.5`$ to the present-day.

In the framework of this model, we find that the observed evolution of the S0/E ratio requires high rates of accretion onto the clusters and also that a significant fraction of all accreted spiral galaxies, both bulge-strong Sab and bulge-weak Scdm, are transformed into S0’s.

The need to transform not only Sab, but also weaker-bulged galaxies into S0’s is a concern as there is no obvious and effective mechanism for transforming an Scdm into an S0. The bulges of late-type spirals are just not luminous enough to evolve into the bright bulges of S0 galaxies simply as a result of stopping the star formation in the disk or even by ‘harassment’ in the cluster environment (Moore et al. 1996). One possibility would be that the accreted field populations are not the same as ‘pure’ field population we assumed from the MDS. They may be ‘pre-processed’ in group or sub-cluster units. or intrinsically biased towards earlier types. Most of the bulge-weak Scdm’s may have been already turned to bulge-strong galaxies by merging or through bar-formation due to strong galaxy harassment/interaction, before entering the clusters. This is likely, given that the outskirts of a rich cluster, where accreting spirals are located, also originate from high density regions in the early Universe (e.g. Moore et al. 1999b). Therefore, the morphological mix in the field local to the cluster may be skewed towards earlier Hubble types than that in the ‘pure’ field through galaxy–galaxy interactions in these denser regions (Moore et al. 1999a, 1999b; Springel 1999; Diaferio et al. 2001).

An alternative solution is to adopt a high density Universe. With $`\mathrm{\Omega }_o=1`$ (q$`{}_{0}{}^{}=0.5`$) for example the accretion rate is higher than in the current models $`A_{\mathrm{recent}}(z=0.5)0.7`$ (Bower 1991). And if this is really the case, the large S0/E ratio at $`z0`$ can be reached without the need to transform Scdm galaxies into S0’s. Furthermore, the S0/Sab ratio would also be consistent even with $`\tau _{\mathrm{trans}}=1`$ Gyr. However, recent observations suggest slower cluster evolution than expected in an $`\mathrm{\Omega }_o=1`$ cosmology, with little evolution seen in X-ray luminosity function of clusters out to $`z0.3`$ (Ebeling et al. 1997) and beyond (Fairley et al. 2000), indicating the recent accretion rate onto rich clusters cannot be that high.

Finally, we note that there is the possibility of significant scatter in the morphological mix in clusters even at the same epoch (Figs. 1 and 2). However, with the present observations this scatter is consistent with the observational errors and so we have chosen to model the evolution of the morphological mix in an ‘average’ cluster. If real differences do exist between the morphological distributions in similar clusters at a fixed redshift this would provide a powerful tool for investigating the nature of galaxy morphology and the mechanisms responsible for its transformation.

In summary, we have simulated the evolution of the morphological mix in clusters with the aim of applying a phenomenological test for the hypothesis of the morphological transformation. We therefore use a simple model in which the accreted field spirals are transformed into S0’s in the cluster environment to attempt to match the observed evolution of the distribution of morphologies for cluster galaxies out to $`z0.5`$. To reproduce the rapid increase in the S0/E ratio between distant clusters and their present-day counterparts, it is found that at least nearly half of the present-day cluster population must have accreted since $`z0.5`$, consistent with the earlier estimate by Fabricant et al. (2000). Moreover, we found that for this hypothesis to be correct, not only the bulge-strong field spirals, but also relatively bulge-weak field spirals (Scdm) must have eventually transformed into cluster S0’s. This requirement may be difficult to achieve unless the field population is pre-processed to produce stronger-bulged galaxies prior to their arrival in the cluster core. Finally, in order to explain the relatively high fraction of cluster spirals we find that the process responsible for the morphological transformation from spiral to S0s upon entry to clusters appears to be relatively slow ($``$1–3 Gyr).

## Acknowledgements

We thank Richard Bower, Michael Balogh, Warrick Couch, Alan Dressler and Bianca Poggianti for useful discussions. We also thank an anonymous referee for their comments which significantly improved the content and presentation of this paper. TK thanks the Japan Society for the Promotion of Science for support through its Research Fellowships for Young Scientists. IRS acknowledges support from the Royal Society.