text:
# Topology of the Las Campanas Redshift Survey

## 1 Introduction

At the present time the topology of the large-scale structure (LSS) of the Universe is not well known. Some surveys show pronounced coherent planar structures such as the “Great Wall” (see e.g. de Lapparent et al. 1986), others show evidence of filaments (see e.g. Haynes and Giovanelli 1986), and yet others appear dominated by superclusters and voids (see e.g. Kirshner et al. 1981). In many cases, however, these conclusions are reached by visual inspection rather than by a rigorous statistical analysis.

The topology of the LSS is very important for cosmology. Simple inflation-based models predict a Gaussian distribution of the density perturbations (see e.g. Mukhanov et al. 1992, for a comprehensive review), whereas topological defect-based models predicted a non-Gaussian density field (see e.g. Brandenberger 1994, for a review). These models to a first approximation all predict a scale-invariant spectrum of primordial perturbations; and therefore the predictions of the models cannot be differentiated by means of power spectra and two-point functions. However, defect models can be differentiated amongst one-another as well as from inflation-based models in terms of the topology of LSS which they yield.

The currently favored simple inflation-based theories predict a Gaussian distribution of the density fluctuations. Nonlinear evolution and biasing are expected to produce secondary non-Gaussianities, but of a type very different from the primary non-Gaussian signals predicted in defect models. For example, the genus curve (one of the Minkowski functionals used in our work) of a simulated SLOAN redshift survey is indistinguishable from that of a Gaussian distribution (Colley et al. 1999).

The nontrivial topological signatures of most defect models are impossible to detect in the current CMB experiments (and hard to detect even in future experiments such as MAP) since due to the central limit theorem the non-Gaussian effects disappear at the large angular scales accessible (see e.g. Moessner et al. 1994, for a specific study). Conversely, the topology of the LSS may be less sensitive to complicated biasing effects than the topology on smaller scales, thus rendering the study of the topology of the LSS an interesting endeavor.

Recently, the use of the Minkowski functional statistics has been suggested as a useful way to study the topology of structure in the Universe (Mecke et al. 1994, Schmalzing and Buchert 1997). It has been shown (Mitsouras 1998, Mitsouras et al. 1999) that Minkowski functionals can differentiate between models with identical power spectra such as defect toy models and inflation-based models. Amongst the four Minkowski functionals $`M0M3`$, the $`M2`$ statistic (integrated mean curvature) was shown to be a better discriminator than the frequently used genus statistic $`M3`$. Minkowski functionals have been applied to Abell galaxy clusters surveys (Kerscher et al. 1997) and to the IRAS catalog (Kerscher et al. 1998), among others, and to numerical simulations (Schmalzing et al. 1999).

Another set of useful statistics to explore the morphology of the large-scale structure are the structure functions (Babul and Starkman 1992, Sathyaprakash et al. 1996) which have e.g. been applied to the 1.2 Jy IRAS catalog (Sathyaprakash et al. 1998). These statistics are based on ratios of eigenvalues of the inertia tensor of the mass distribution and carry information different from that of the Minkowski functionals. On the other hand, the Shapefinder statistics (Sahni et al. 1998) are ratios of Minkowski functionals. They have very recently (Bharadwaj et al. 2000) been used to analyze the LCRS data set (without taking into account the important selection effects), finding evidence for filamentarity.

The deepest wide field redshift survey to date is the Las Campanas Redshift Survey (LCRS) (Shectman et al 1996). An initial analysis (Doroshkevich et al 1996) indicated that the LCRS field was dominated by planar structures, with higher density filaments at their intersections. In contrast, the analysis of Mitsouras (Mitsouras 1998) using the Minkowski functional statistics detected no deviations from a Gaussian distribution. Neither analysis, however, took into account the complicated selection effects intrinsic to the LCRS survey.

In this paper, we calculate the Minkowski functionals for the LCRS, taking careful account of the nontrivial selection effects of that survey (except for the limited surface brightness selection effect). We compare the resulting statistics with the predictions of some toy models of structure formation: the Standard Cold Dark Matter (SCDM) model as a representative of the class of inflation-based models, and three toy defect models (based on string wakes, string filaments, and textures, respectively). Our aim in this paper is to focus on the potential of the LCRS survey to differentiate between different models of structure formation with identical power spectra rather than on studying specific COBE normalized models. We find that no model comes close to the data. In work in progress we are studying a class of COBE-normalized inflation-based models with the goal of discovering whether some members of this class come closer to explaining the observed topology than the SCDM model.

## 2 Data Analysis

The Las Campanas Redshift Survey (LCRS) is a fiber-optic survey consisting of 26418 redshifts of galaxies. $`1.5^{}\times 1.5^{}`$ spectroscopic fields were observed to pave 700 deg2 of sky in six strips, each roughly $`1.5^{}\times 80^{}`$. Three strips centered on $`\delta =3^{}`$, $`6^{}`$, and $`12^{}`$ are in the North galactic cap (NGC) and three centered on $`\delta =39^{}`$, $`42^{}`$, and $`45^{}`$ are in the South galactic cap (SGC). The LCRS probed to redshifts $`z0.2`$ with mean redshift $`\overline{z}0.1`$ or equivalently to distances of 600 $`h^1`$Mpc. Approximately 70% of galaxies having magnitudes in the range $`15.0m_r<17.7`$ were measured spectroscopically. For our analysis, we only consider galaxies out to a redshift $`z<0.15`$ or distances of 450 $`h^1`$Mpc or less.

The selection criteria for the LCRS are more complex than a simple magnitude limit and have their origin in the technique of the fiber-optic measurements. Galaxies were selected from a CCD-based catalog based on four selection criteria. For the first 20% of survey galaxies, redshifts were measured using a 50-fiber spectrograph system and the final 80% using an improved 112-fiber system. Galaxies were first selected by applying both faint and bright isophotal magnitude limits, with nominal values

$$\begin{array}{c}m_1m<m_2:\hfill \\ \{\begin{array}{c}m_1=16.0,m_2=17.3,50\mathrm{fiber}\mathrm{data}\hfill \\ m_1=15.0,m_2=17.7,112\mathrm{fiber}\mathrm{data}\hfill \end{array}.\hfill \end{array}$$
(1)
The second criterion selects against low surface brightness (LSB) galaxies by imposing the nominal central surface brightness cutoff

$$\begin{array}{c}m_c<m_{cen}0.5(m_2m):\hfill \\ \{\begin{array}{c}m_{cen}=18.15,50\mathrm{fiber}\mathrm{data}\hfill \\ m_{cen}=18.85,112\mathrm{fiber}\mathrm{data}\hfill \end{array},\hfill \end{array}$$
(2)
where the central magnitude $`m_c`$ corresponds to the flux within a 2 pixel radius of the centroid of the galaxy. There are small field-to-field variations from the nominal photometric limits which were made in response to variations in sky brightness and in seeing. The third criterion restricts the maximum number of galaxies in a given field to the number of fibers available during observation of that particular field. Lastly, the fourth criterion imposes a minimum separation of 55′′ between galaxies because mechanical constraints prevent the fibers from being placed any closer together.

## 3 Models

There are two main classes of theories for structure formation which may possibly be able to account for the origin of primordial perturbations and observed large-scale structure in the universe. The first are inflation-based models in which quantum fluctuations are generated at an early stage of inflation on microphysical scales, and the wavelength of these perturbations subsequently grows exponentially and thus gives rise to density perturbations on cosmological scales (Mukhanov and Chibisov 1981, see also Press 1980, Lukash 1980, Chibisov & Mukhanov 1980, Sato 1981, Chibisov & Mukhanov 1982 for related work). Simple inflation-based models predict a scale-invariant spectrum of adiabatic perturbations with random phases, each having a Gaussian distribution. The second class is based on topological defects, formed during a symmetry breaking phase transition in the early universe, giving rise to seeds for structure formation (Kibble 1976, Zel’dovich 1980, Vilenkin 1981). Topological defects represent regions in space with trapped energy density and they result in a roughly scale-invariant spectrum of non-adiabatic perturbations with non-Gaussian phases. In this paper, we compare the inflation-based “standard” cold dark matter (SCDM) model with topological defect models like cosmic strings and textures. 111After this work was completed, the results of the BOOMERANG (de Bernardis et al 2000) and MAXIMA-I (Hanany et al. 2000) CMB experiments were announced, which with high statistical significance demonstrate the existence of a well-defined and narrow first acoustic peak in the angular power spectrum, thus ruling out topological defects as the main source of density fluctuations (see e.g. Albrecht et al 1999 and references therein). It is, however, important to keep in mind that defects may still contribute a significant amount to the primordial power spectrum (see e.g. Bouchet et al. 2000), and that the resulting non-Gaussian signals in the large-scale structure of the Universe may well be important.

### 3.1 Standard Cold Dark Matter Model (SCDM)

In the SCDM model, the dominant constituent of the universe is cold dark matter and the standard model is described by a density parameter $`\mathrm{\Omega }=1`$, Hubble constant $`H_{}=100h`$ km/s/Mpc with $`h=0.5`$, and cosmological constant $`\mathrm{\Lambda }=0`$. This specific model conflicts with observations if we compare the normalization of the power spectrum of density perturbations with that of the CMB fluctuations. Other inflationary models like tilted CDM, $`\mathrm{\Lambda }`$CDM, and mixed dark matter models have the potential to account for observations. Since we are primarily concerned with the power of the Minkowski functionals to differentiate between the signatures of primordial Gaussian and non-Gaussian phases, we have chosen to work with the standard model. We have also modified the standard value of the Hubble constant from $`h=0.5`$ to $`h=0.7`$ to reflect some of the more recent measurements.

In the SCDM simulation, 1923 particles are evolved in a cube of comoving length of 630 $`h^1`$Mpc using the Bertschinger & Gelb (1991) $`P^3M`$ N-body code. The normalization of the primordial power spectrum was chosen such that $`\sigma _8=0.95`$ to agree with the observed clustering (Pen 1998). A COBE normalized model would require $`\sigma _8=1.22`$.

### 3.2 Cosmic Strings

Cosmic strings are linear topological defects which arise if the vacuum manifold $``$ is not simply connected. At the time of the phase transition, a network of strings is formed. The phase transition takes place on a short time scale $`\mathrm{\Delta }t<H^1`$, and will lead to correlation regions of radius $`\xi <t`$. The strings are either infinite in length or closed loops. Infinite strings give rise to filaments or wakes while loops form seeds for spherical accretion of matter and radiation. In our study, we are concerned only with the infinite strings, but realistic cosmic string models should contain a mixture of both.

In the string wake model, strings are straight on the scale of the Hubble radius and the string tension $`|p|`$ equals the mass per unit length $`\mu `$. As a result, the strings will move relativistically and exert no local gravitational force. The strings will have a large transverse velocity $`v`$ and give rise to a planar velocity perturbation. A string at time $`t_i`$ produces a wake of comoving length $`v\xi (t_i)z(t_i)`$ and width $`\xi (t_i)z(t_i)`$. The wake thickness can be calculated using the Zel’dovich approximation (Zel’dovich 1970) and will depend on the nature of the dark matter since free streaming has an effect. For hot dark matter, the thickness $`h(t_i)`$ at the present time $`t_{}`$ of a wake created at time $`t_i>t_{eq}`$ is (Perivolaropoulos et al. 1990, Brandenberger 1991)

$$h(t_i)=\frac{24\pi }{5}G\mu v\gamma (v)z(t_i)^{1/2}t_{},$$
(3)
where $`\gamma (v)`$ is the relativistic factor associated with $`v`$. The cosmic string model of structure formation is viable even if dark matter is hot since free streaming will not erase the non-adiabatic primordial perturbations on a small scale. The wakes will have dimensions

$$v\xi (t_i)z(t_i)\times \xi (t_i)z(t_i)\times h(t_i),$$
(4)
with mass

$$m(t_i)=\frac{24\pi }{5}G\mu v^2\gamma (v)\nu ^2\rho _{}t_i^{1/3}t_{}^{8/3}.$$
(5)
The string correlation length is taken to be $`\xi (t)=t/\nu `$, where $`\nu `$ could in principle be determined from numerical simulations. Hot dark matter wakes created before $`t_{eq}`$ are diluted by the radiation pressure and, therefore, do not contribute any significant amount to structure formation. The string wake model predicts that planar wakes are the dominant structures in the large-scale galaxy distribution with the most numerous and thickest wakes created at $`t_{eq}`$.

In the string filament model, strings have a lot of small-scale structure and the string tension $`|p|`$ will be smaller than $`\mu `$. The strings will have a small transverse velocity and will exert local gravitational attraction. This induces a filamentary pattern of accretion. Zanchin et al. (1996) have analyzed the string filament model for moderate transverse velocities and we will here present an extreme case in which the transverse velocity is negligible (Mitsouras et al. 1999). Using the Zel’dovich approximation, the comoving radius $`q(t_i)`$ of a filament created at time $`t_i`$ is (Aguirre & Brandenberger 1995)

$$q(t_i)=\left[\frac{12}{5}G\lambda \mathrm{ln}\left(\frac{t_{}}{t_i}\right)\right]^{1/2}t_{},$$
(6)
where $`\lambda =\mu |p|`$. The mass in an individual filament is then

$$m(t_i)=\frac{12\pi }{5}G\lambda \nu ^1\rho _{}t_i^{1/3}t_0^{8/3}\mathrm{ln}\left(\frac{t_{}}{t_i}\right),$$
(7)
where $`\nu `$ is defined as previously.

In the toy topological defect simulations, galaxies are just laid down and no actual evolution takes place. The procedure for setting up a string filament or string wake model first required dividing the time interval between $`t_{eq}`$ and $`t_{}`$ into Hubble times. At each time step $`t_i`$, a network of strings with correlation length $`\xi (t_i)`$ is laid down to simulate the evolution of the strings. For each time step $`t_i`$, beginning with $`t_{eq}`$, we divided the simulation cube into $`N(t_i)`$ Hubble volumes $`t_i^3`$. A fixed number $`c_s^3`$ of strings per Hubble volume is laid down, choosing string centers and directions at random in the cube. Thus at each time step, $`c_s^3\times N(t_i)`$ is the number of strings being laid down in the cube. For given values of $`G\lambda `$ and $`\nu `$ in the case of the string filament simulation and $`G\mu `$, $`\nu `$ and $`v`$ in the case of the string wake simulation, the total number $`c_s`$ of strings is determined by demanding that the total mass satisfies critical density. Galaxies are then laid down at random in the volume of the strings, with a uniform density amongst the strings. For our string filament simulation, we have taken $`G\lambda =10^6`$ to agree with the COBE normalization and $`\nu =1`$, such that the length of the string is equal to the Hubble radius at its time of creation. The string wake simulation has values $`G\mu =10^6`$ and $`v\gamma =1`$.

### 3.3 Textures

Textures are topological defects whose vacuum manifold $``$ is a sphere $`S^3`$. In 4-dimensional space-time, textures are unstable (Turok 1989) and can collapse to form seeds for structure formation. While the energy in cosmic strings is trapped in potential energy, in textures, the energy is in spatial gradients with possibly some in kinetic energy.

The simplest texture configurations are spherically symmetric (Turok 1989). At each time $`t_i`$, there is a finite probability $`p`$ that in a volume $`t_i^3`$, there will be a texture configuration entering the horizon (Prokopec 1991, Leese & Prokopec 1991). The initial texture configuration can be realized by a spherically symmetric velocity perturbation which results in a spherically symmetric density perturbation (Gooding et al. 1991). The comoving radius $`q(t_i)`$ of this perturbation, calculated (Aguirre 1995) using the Zel’dovich approximation, is

$$q(t_i)=\frac{6}{5}ϵt_i^{1/3}t_{}^{4/3},$$
(8)
where $`ϵ=G\eta ^2`$ and $`\eta `$ is the scale of symmetry breaking.

In the toy texture model, $`p`$ spherical texture configurations are laid down per Hubble volume $`t_i^3`$ with centers chosen at random in the box. Galaxies are laid down within each texture with a Gaussian radial density profile and random angular distribution. The number of galaxies per texture was chosen to satisfy critical density for the entire volume. A value of $`ϵ=3\times 10^5`$ is an appropriate value for a COBE-normalized texture model (Pen et al. 1997).

### 3.4 Construction of mock LCRS catalogs

We construct mock catalogs by applying the LCRS selection criteria to the simulations. The first step involved randomly assigning luminosities and corresponding absolute magnitudes to all the galaxies in each cube, according to the Schechter luminosity function

$$\varphi (L)=\varphi ^{}\left(\frac{L}{L^{}}\right)^\alpha e^{L/L^{}}.$$
(9)
where $`\varphi ^{}=0.019h^3\mathrm{Mpc}^3`$, $`\alpha =0.70`$, and $`L^{}`$ corresponds to the absolute magnitude $`M^{}=20.29+5\mathrm{log}h`$ for the LCRS (Lin et al. 1996). The faint end slope $`\alpha =0.70`$ results from fitting the LCRS luminosity function and not from the absence of intrinsically faint galaxies in the survey. The next step was to assign a surface brightness to the galaxies, but this presented a significant problem since the distribution of surface brightness is nontrivial. While LSB galaxies are found in the same large scale structure as other galaxies, they tend to be less strongly clustered (Mo et al. 1994). Since biasing may occur in one or more of the models, the LSB selection criterion has been omitted.

To apply the selection criteria, slices of appropriate geometry are first cut from the simulation cubes. We have chosen to consider the three strips in the North galactic cap as a composite slice and the three strips in the South galactic cap as another composite slice. The 3-dimensional geometry of the NGC and SGC slices should reveal more information about the topology of the galaxy distribution. To avoid boundary problems in the 630 $`h^1`$Mpc simulation cube, slices are cut from an interior cube with side-lengths of 490 $`h^1`$Mpc. After converting the absolute magnitudes into apparent magnitudes, the exact faint and bright apparent magnitude limits in each spectroscopic field (Shectman et al. 1996) are applied. Next, the restriction on the maximum number of galaxies in each spectroscopic field is imposed, depending on the fiber system used for that particular field (random elimination of galaxies if the number of galaxies in the simulation field exceeds the maximal number of fibers used for that field). Finally, to apply the minimum separation criterion without introducing bias, the galaxies are ordered in RA to form the set $`\{g_i|g_1,\mathrm{},g_N\}`$. Starting with galaxy $`g_1`$ and choosing to keep it, any galaxy within $`55^{\prime \prime }`$ is immediately removed from the set. This step is iteratively carried out until the end of the set is reached.

## 4 Minkowski Functionals

Hadwiger’s theorem states that in d spatial dimensions, the global morphological description of any spatial distribution can be completely characterized by d+1 descriptors comprising the set of Minkowski functionals (Hadwiger 1957). Given an isodensity surface $`K`$ in 3d, the morphology is completely determined by the following four statistics. The first Minkowski functional

$$M_{}=V,$$
(10)
is simply the volume enclosed by the surface $`K`$ and the second functional

$$M_1=\frac{A}{8},$$
(11)
is related to the surface area $`A`$. The third functional

$$M_2=\frac{H}{2\pi _2},$$
(12)
provides information about shape and is proportional to the integrated mean curvature

$$H=\frac{1}{2}_K\left(\frac{1}{R_1}+\frac{1}{R_2}\right)𝑑A,$$
(13)
where $`R_1`$ and $`R_2`$ are the principal radii of the surface $`K`$. The fourth Minkowski functional

$$M_3=\frac{3\chi }{4\pi },$$
(14)
which is related to the Euler characteristic

$$\chi =_K\frac{1}{R_1R_2}𝑑A.$$
(15)
is purely a topological quantity and can be geometrically defined as

$$\chi (\mathrm{\#}\mathrm{of}\mathrm{clusters})(\mathrm{\#}\mathrm{of}\mathrm{tunnels})+(\mathrm{\#}\mathrm{of}\mathrm{cavities}).$$
(16)
The Minkowski functionals are realized for a galaxy distribution by placing spheres of radius $`r`$ centered on the galaxies. The union of these spheres specifies an ensemble of clusters, where two spheres are said to belong to the same cluster if they are connected by a chain of intersecting spheres. The Minkowski functionals will be used to measure the topology of the coverage. We use the method of Mecke et al. (1994) and Kerscher et al. (1997) and display the functionals as a function of the length scale r, which acts as a smoothing length of the density field.

## 5 Results

For each model studied, the simulations were performed in the following way. First, a cubical box of comoving length 630 $`h^1`$Mpc was populated with galaxies according to the prescriptions discussed in Section 3. Next, slices with the appropriate geometry of the LCRS survey were cut out from the simulation cube. To avoid boundary effects, the slices were extracted from an interior cube with side length of 490 $`h^1`$Mpc. For comparing with the individual LCRS slices, four galaxies in the interior box were chosen as centers of slightly oversized slices of $`90^{}`$ in length and $`1.6^{}`$ in width. For comparing with the composite LCRS slices, four larger slices of $`90^{}`$ in length and $`10.6^{}`$ in width were cut out.

Luminosities were assigned randomly in the initial simulation box according to the Schechter luminosity function. For each slice, this then allows us to determine the apparent magnitude of all galaxies contained in the slice. Next, the LCRS selection criteria were applied as detailed in Section 3.4.

For each slice, the Minkowski functionals were calculated using a slightly modified version (Mitsouras 1998) of the program package publically available (Kerscher et al. 1998). Averages over the four centers were taken. Under the ergodic hypothesis, the volume averages (averages over the slices) will be equal to the ensemble averages over different realizations of the theories. The error bars indicated in our plots correspond to these statistical errors.

Figures 1 - 4 show the resulting volume-averaged Minkowski functionals for the NGC slice. The results of our simulations are compared with the LCRS data. As mentioned in the Introduction, none of the models considered in this paper is able to explain the data. The SCDM model gives the least bad fit, but even it is ruled out with high statistical significance. 222The oscillations seen at large values of $`r`$ are artifacts of the boundary conditions.

As is apparent from Figures 1 - 4, the Minkowski functionals provide a very powerful tool to distinguish between models with identical power spectra but different topology. Amongst the Minkowski functionals, the integrated mean curvature statistic M2 appears to be the most powerful discriminant, followed by the genus statistic M3. This confirms the results seen in previous studies (Mitsouras 1998, Mitsouras et al. 1999, Kerscher et al. 1997).

The predictions of the texture toy model differ dramatically from those of the other models studied, and also from the data. The reason is that in the texture toy model, structure is dominated by well separated spherical clumps, whereas in the other models there is a higher degree of connectedness of the individual structures. In addition, the central density in the texture toy model is larger than in the other theories studied. Therefore, M0, the volume enclosed by the “isodensity” surface, is significantly smaller for any value of $`r`$ than in the other models. The surface area (M1) itself continues to increase with $`r`$ since to a first approximation the number of disconnected components remain the same while the size of each object increases. In the other models studied, the value of M1 for small $`r`$ is higher since the shape of the objects does not minimize the surface area for fixed enclosed volume. However, as $`r`$ increases, the growing interconnectedness of the components leads, for large $`r`$, to a gradual decrease in the value of M1. The radius corresponding to the maximum of M1 corresponds to the mean separation between clusters.

For a single component surface without holes and with fixed volume, the sphere minimizes the integrated mean curvature. Hence, for small values of $`r`$, the value of M2 is lower in the texture model than in the other models. In particular, in the string filament and wake toy models the individual filaments and wakes have a small curvature radius in one direction, and hence a large mean curvature. Once the individual filaments and wakes merge by overlapping, the overall shape of the surface becomes concave and hence M2 turns negative. This feature is shared by the SCDM model, and by the data. Similarly, the overlapping of filaments and wakes leads to a structure with many holes for larger values of $`r`$. This leads to negative genus for these scales, a feature which is again shared by the SCDM model and by the data, but is absent in the texture toy model. The radius corresponding to the minimum of M3 gives another good indicator of the mean separation of clusters.

By comparing the LCRS data with the simulations, it appears that the data has less connectivity than what would be seen in the SCDM model, and significantly less than what would be predicted by the toy string theories. These results are confirmed by the analysis of the SGC slice (Figures 5 - 8). Note that the differences between the NGC and SGC slices are due to the different selection effects.

The number of galaxies in the model slices is not the same as in the LCRS slices. In order to check that the disagreement between the theories and the data is not due to this difference, we have repeated the analysis with “randomly sampled” model slices, model slices from which in a random way galaxies were removed in order to obtain the same number of galaxies as in the actual LCRS slices. Figures 9 - 12 give the results for the “randomly sampled” NGC slice, Figures 13 - 16 the corresponding results for the SGC slice. As is apparent, the conclusions are unchanged.

As mentioned in Section 3.4, we did not apply the LSB selection criterion of the observational survey. To investigate whether this omission could account for the disagreement between the data and the SCDM model, we repeated the calculation of the Minkowski functionals of the observational data with the lowest $`20\%`$ of the galaxies (in terms of surface brightness) removed. There were no appreciable changes in the results.

## 6 Discussion

The main conclusion of our study is that Minkowski functionals provide a very powerful statistical tool to differentiate between theories with identical power spectra but different topologies. We have compared the predictions of the SCDM model and of three topological defect toy models with the LCRS data, taking careful account of most of the selection criteria of the LCRS survey.

None of the models studied comes close to giving a satisfactory match to the data. This is not surprising since the SCDM model is known not to correctly describe the observed clustering of matter in the Universe, and since the defect models considered here are all unrealistic toy models. The texture model produces isolated and well-separated dense clusters and gives a picture of the distribution of galaxies which already at first glance is inconsistent with any observational data set. The string models are based on ideal straight filaments or wakes, and the filaments and walls seen in observations are not straight.

It would be of great interest to explore whether more realistic models can provide a better fit to the data. In the class of inflation-based models yielding a spectrum of roughly scale-invariant primordial perturbations it would be of interest to explore the effects of reducing $`\mathrm{\Omega }_m`$ and introducing a cosmological constant $`\mathrm{\Lambda }`$. In particular, it would be of great interest to see if the Minkowski functionals as applied to the LCRS geometry can differentiate between the set of inflation-based models which are compatible with the present CMB data and with the amplitude of the matter power spectrum on large scales. Work on this issue is in progress (Merali et al. 2000).

From the point of view of particle cosmology, the defect models studied in this paper are too idealized to be realistic. In texture models the resulting structures are not all spherically symmetric, and there are important non-topological fluctuations which also lead to density inhomogeneities (Borrill 1994). Similarly, the structure seeded in a cosmic string model will also include effects of string loops (see e.g. Brandenberger 1994 for a review), and the structures seeded by the long string network will be given by a combination of the filament and wake effects. In addition, the structures will not be straight. It would also be of interest to use the Minkowski functionals to study the predictions of numerically simulated more realistic defect models, and to see if there are parameters for which good agreement with the data is obtained.

It is also possible that the differences between the theoretical predictions and the LCRS data are due to biasing effects. This issue could be investigated by using hydrodynamical simulations of structure formation. Redshift distortions also complicate the interpretation of clustering in redshift surveys. The peculiar velocities of galaxies cause them to appear displaced along the line of sight in redshift space, thus causing the observed galaxy distribution to be anisotropic. It would be interesting to compare the real and redshift space results from N-body simulations using Minkowski functionals.

In conclusion, we hope to have convinced the reader that Minkowski functionals provide a very powerful tool for quantitative cosmology.


format: text

text:
# The cluster environments of the 𝐳∼𝟏 3CR radio galaxies

## 1 Introduction

Clusters of galaxies are the largest, most massive, collapsed structures in the Universe, and as such are of fundamental importance for many cosmological studies. They provide a unique probe of large–scale structure in the early Universe and, as the systems which separated from the Hubble flow at the earliest epochs, they contain the oldest galaxies known. These can strongly constrain the first epoch of the formation of ellipticals, and hence set a lower limit to the age of the Universe; the effectiveness of using such old galaxies at high redshifts towards this goal has been well demonstrated by Dunlop et al. \[Dunlop et al. 1996\]. Further, because clusters contain large numbers of galaxies at the same distance, they are important testbeds for models of galaxy evolution.

The cores of optically–selected clusters are dominated by a population of luminous early–type red galaxies which occupy a narrow locus in colour–magnitude relations. Stanford, Eisenhardt and Dickinson \[Stanford et al. 1998\] showed that, out to redshifts $`z1`$, the evolution of the colours of early–type cluster galaxies on these relations is completely consistent with passive evolution of an old stellar population formed at high redshift. The small intrinsic scatter of the galaxy colours, and particularly the fact that it remains small at redshifts $`z\mathrm{}>0.5`$ \[Stanford et al. 1998\], implies that the star formation of the ellipticals comprising a cluster must have been well synchronized, and sets tight limits on the amount of recent star formation that might have occurred (e.g. Bower, Kodama and Terlevich 1998).

Cluster ellipticals also show a tight relationship between their effective radius, effective surface brightness, and central velocity dispersion (the ‘fundamental plane’; c.f Dressler et al. 1987, Djorgovski and Davies 1987). The location of these ellipticals within the fundamental plane has been shown to evolve with redshift out to $`z=0.83`$, in a manner which implies that the mass–to–light ratio of the galaxies evolves as $`\mathrm{\Delta }\mathrm{log}(M/L_\mathrm{r})0.4\mathrm{\Delta }z`$, roughly in accordance with passive evolution predictions (e.g. van Dokkum et al 1998 and references therein).

These results are in qualitative agreement with ‘monolithic collapse’ models of galaxy formation, in which an elliptical galaxy forms the majority of its stars in a single short burst of star formation at an early cosmic epoch. On the other hand, the hierarchical galaxy formation models favoured by cold dark matter cosmologies (e.g. White and Frenk 1991) predict a later formation of galaxies with star formation on–going to lower redshifts. These are supported by the appearance of a population of bluer galaxies in many clusters at redshifts $`z\mathrm{}>0.3`$ (the Butcher–Oemler effect; Butcher and Oemler 1978) and by the high fraction of merging galaxies seen in the $`z=0.83`$ cluster MS1054$``$03 \[van Dokkum et al. 1999\]. Since finding clusters at high redshift is strongly biased towards the very richest environments, in which early–type galaxies will have formed at the very highest redshifts, hierarchical models can still explain the apparent passive evolution and small scatter of the colour–magnitude relation out to redshifts $`z1`$.

It is clearly important to extend cluster studies out to still higher redshifts, but the difficulty lies in the detection of clusters. At optical wavelengths, the contrast of a cluster above the background counts is minimal at these redshifts: the deep wide–area ESO Imaging Survey (EIS) has found 12 ‘good’ cluster candidates above redshift 0.8, but none above $`z=1.2`$ \[Scodeggio et al. 1999\]. Addition of near–infrared wavebands helps (e.g. Stanford et al. 1997), but is still a relatively inefficient method. Selection using X–ray techniques is more reliable, but X–ray surveys are currently sensitivity limited: the ROSAT Deep Cluster Survey found about 30 clusters with $`z\mathrm{}>0.5`$, but none above redshift 0.9 \[Rosati et al. 1998\]. Chandra and XMM will make a big improvement here. An alternative approach is to use powerful radio galaxies as probes of distant clusters: these can be easily observed out to the highest redshifts, and there is growing evidence that at high redshifts they lie in rich environments.

At low redshifts powerful double radio sources (FR II’s; Fanaroff & Riley 1974) are associated with giant elliptical galaxies that are typically the dominant members of galaxy groups; the only nearby radio source of comparable radio luminosity to the powerful high redshift radio galaxies is Cygnus A, and this source lies in a rich cluster \[Owen et al. 1997\]. At a redshift $`z0.5`$, analysis of the galaxy cross–correlation function around FR II radio galaxies (Yates, Miller and Peacock 1989) and an Abell clustering classification \[Hill & Lilly 1991\] have shown that about 40% of radio sources are located in clusters of Abell richness class 0 or greater. At $`z1`$, the circumstantial evidence that at least some powerful radio sources are located at the centres of clusters is overwhelming and includes the following:

* Detections of luminous X–ray emission from the fields of the radio galaxies (e.g. Crawford and Fabian 1996), sometimes observed to be extended, indicating the presence of a relatively dense intracluster medium.
* Large over-densities of galaxies in the fields of some distant radio sources, selected by infrared colour \[Dickinson 1997\].
* Direct detections of companion galaxies with narrow–band imaging (McCarthy, Spinrad and van Breugel 1995) and spectroscopic studies \[Dickinson 1997\].
* The observation that powerful radio galaxies are as luminous as brightest cluster galaxies at $`z0.8`$, and have radial light profiles which are well–matched by de Vaucouleurs law with large (10 to 15 kpc) characteristic radii (Best, Longair and Röttgering 1998b).
* The radio sources display large Faraday depolarisation and rotation measures (e.g. Best et al. \[Best et al. 1998a\]; see also Carilli et al. \[Carilli et al. 1997\] for redshift $`z>2`$ radio sources), requiring a dense, ionised surrounding medium.
* Theoretical arguments that to produce such luminous radio sources requires not only a high AGN power, due to a very massive central black hole being fueled at close to the Eddington limit \[Rawlings & Saunders 1991\], but also a dense environment to confine the radio lobes and convert the jet kinetic energy efficiently into radiation (e.g. see Barthel and Arnoud 1996).

At still higher redshifts, radio sources have been detected out to redshift $`z=5.2`$ \[van Breugel et al. 1999\], and some well–studied sources are known to lie in cluster environments. For example, towards the radio galaxy 1138$``$215 ($`z=2.2`$) extended X–ray emission has been detected (Carilli et al. 1998) and narrow–band imaging reveals over 30 nearby Ly-$`\alpha `$ emitters \[Kurk et al. 2000\]. Radio sources may therefore offer a unique opportunity to study dense environments back to the earliest cosmic epochs. As yet, however, there have been no systematic studies of radio galaxy environments much beyond $`z0.5`$, and so it is important to investigate in detail the nature of the environments of the general population: do all powerful distant radio galaxies lie in cluster environments, or do only a minority but these have ‘grabbed the headlines’?; what are the properties of any clustering environments (richness, radius, shape, etc.) surrounding these objects?; what is the nature of the constituent galaxies (morphological composition, segregation, etc.) of any detected cluster?

In this paper, deep near–infrared observations of the fields of a sample of 28 powerful radio galaxies with $`0.6<z<1.8`$ are analysed, in conjunction with Hubble Space Telescope (HST) images of the same fields, to investigate the ubiquity and richness of the environments of the radio galaxies. The observations, data reduction, and source extraction are described in Section 2. In Section 3 the integrated galaxy counts are considered, in Section 4 the angular and spatial cross–correlation amplitudes are derived, and in Section 5 an investigation of the colour–magnitude and colour–colour relations is carried out. The implications of the results are discussed in Section 6. Throughout the paper values for the cosmological parameters of $`\mathrm{\Omega }=1`$ and $`H_0=50`$ km s-1Mpc-1 are assumed.

## 2 Observations and Data Reduction

### 2.1 The dataset

The data used for this research was presented and described by Best, Longair and Röttgering \[Best et al. 1997\]. In short, the sample consists of 28 radio galaxies with redshifts $`0.6<z<1.8`$ drawn from the revised 3CR sample of Laing, Riley and Longair \[Laing et al. 1983\]. The fields of these radio galaxies were observed at optical wavelengths using the Wide-Field Planetary Camera II (WFPC2) on the HST generally for one orbit in each of two different wavebands. They were also observed at near–infrared wavelengths using IRCAM3 on the UK Infrared Telescope (UKIRT) in the K–band for approximately 54 minutes, and in 20 of the 28 cases also in the J–band.

For 5 sources (3C13, 3C41, 3C49, 3C65, 3C340) a further 3 to 4 hours of K–band observations have subsequently been taken using IRCAM3 in September 1998 (Best et al., in preparation), and these were combined with the original data to provide much deeper images. These further data were taken in the same manner as the original IRCAM3 data except for the use of the tip-tilt system that was available on UKIRT then, but wasn’t available for the original runs; use of this system provided a significant reduction in the effective seeing.

The HST data were reduced according to the standard Space Telescope Science Institute (STScI) calibration pipeline \[Lauer 1989\], following the description of Best et al. \[Best et al. 1997\]. The UKIRT data were reduced using IRAF and mosaicked using the dimsum111dimsum is the ‘Deep Infrared Mosaicking Software’ package developed by Eisenhardt, Dickinson, Stanford and Ward. package, again following the general procedure outlined by Best et al. \[Best et al. 1997\]. During the mosaicking process, the individual galaxy frames were block replicated by a factor of 4 in each dimension to allow more accurate alignment. The final images were then block–averaged in 2 by 2 pixels, resulting in frames with a pixel size of 0.142 arcseconds. The seeing varied from image to image, from about 0.7 to 1.2 arcsec. dimsum also produces an output exposure map indicating the exposure for each pixel on the combined image; this is useful for weighting the different regions of the images, which reach different limiting sensitivities owing to the dithering process used in taking the data (see Best et al. 1997).

This dataset has a number of advantages and disadvantages for use in studies of clustering around distant radio galaxies. On the negative side, since these data were originally taken with the goal of studying the host galaxies of the radio sources, no comparison fields were taken with which the data can be compared to statistically remove background counts. Further, IRCAM3 is only a 256 by 256 array with a field of view per frame of just over 70 by 70 arcseconds. After mosaicking this provided a final image of nearly 100 by 100 arcsec, which corresponds to a little over 800 kpc at a redshift of one for the adopted cosmology, but the highest sensitivity is only obtained in the central 50 by 50 arcsec region. This is a significantly smaller size than a cluster would be expected to have at this redshift, and so only those galaxies in the central regions of any prospective clusters will be investigated (cf. the much larger extent of the structure around the radio galaxy 3C324 at redshift 1.206, as demonstrated in Figure 2 of Dickinson \[Dickinson 1997\], although just over half of the associated cluster members found there would lie within the IRCAM3 field of view).

Although these two factors limit the usefulness of the dataset, they are far outweighed by the benefits. First, the large size and essentially complete nature of the sample make it ideal for surveying the average environments of these sources. Second, the multi–colour nature of the dataset provides important information on any cluster membership through the creation of colour–magnitude diagrams and the determination of photometric redshifts. Third, the availability of near–infrared data is essential for such studies at these high redshifts since it continues to sample the old stellar populations longwards of the 4000Å break; in the K–waveband the K and evolutionary corrections are small and relatively independent of morphological type, even at redshifts $`z\mathrm{}>1`$, in contrast to the situation at optical wavelengths. Fourth, the availability of the high resolution HST data allows accurate star–galaxy separation down to the faintest magnitudes studied. For these reasons, this dataset is well–suited to investigating the environments of $`z1`$ radio galaxies.

### 2.2 K–band image detection and photometry

Throughout these analyses, the K–band data were used as the primary dataset. Image detection and photometry on the K–band frames was carried out using sextractor version 2.1 \[Bertin & Arnouts 1996\]. As a result of the dithering technique used to obtain near–infrared data, the exposure time varies with position across the image meaning that the rms background noise level also varies. sextractor allows the supply of an input weight map by which the local detection threshold is adjusted as a function of position across the image to compensate the varying noise levels, thus avoiding missing objects in the most sensitive regions of the image or detecting large numbers of spurious features in the noisiest regions. The output exposure map of each field produced by dimsum was used as the weight map for sextractor. The source extraction parameters were set such that, to be detected, an object must have a flux in excess of 1.5 times the local background noise level over at least $`N`$ connected pixels, where $`N`$ was varied a little according to the seeing conditions, but in general was about 20 (equivalent to 5 pixels prior to the block replication and averaging during the mosaicking procedure).

To test the validity of this extraction method, a search for negative holes was carried out using the same extraction parameters, and resulted in a total of only 43 negative detections throughout the 28 images, that is, about 1.5 negative detections per frame. The fluxes associated with these negative detections all correspond to positive features below the 50% completeness limit and so it is expected that there are essentially no false positive detections above that limit.

The output from sextractor was examined carefully, with minor (typically 1–2 per field) modifications being made to the catalogue if necessary: objects coinciding with the spikes of bright stars were removed, the data for occasional entries that had been separated into two objects by sextractor but were clearly (on comparison with data at other wavebands) a single object were combined, and very occasionally it was necessary to add to the catalogue an obvious object which had been missed due to its proximity to a brighter object. Further, all objects that lay within 21 pixels (3 arcsec) of the edge of an image were removed from the catalogue in case their magnitudes were corrupted by aperture truncation. sextractor provides a flag parameter to indicate the reliability of its measured magnitudes. After this trimming of objects close to the edges of the fields, less than 2% of all objects had high flag values indicative of truncated apertures, saturated pixels, or corrupted data. These objects were also removed from the catalogue.

sextractor’s mag\_best estimator was used to determine the magnitudes of the sources; this yields an estimate for the ‘total’ magnitude using Kron’s \[Kron 1980\] first–moment algorithm, except if there is a nearby companion which may bias the total magnitude estimate by more than 10% in which case a corrected isophotal magnitude is used instead. The determined magnitudes were also corrected for galactic extinction, using the extinction maps of Burstein and Heiles \[Burstein & Heiles 1982\].

To investigate the accuracy of these total magnitudes, and the completeness level of the source extraction as a function of position on the image, Monte–Carlo simulations were carried out using the following 4–step process.

Step 1: The point–spread function (PSF) of each K–band image was determined using objects which were unsaturated and which were unresolved on the HST images.

Step 2: A series of model galaxies was made. 30% of the galaxies were assumed to be ellipticals, with radial light profiles governed by de Vaucouleurs’ law, $`I(r)\mathrm{exp}\left[7.67\left((a/a_\mathrm{e})^2+(b/b_\mathrm{e})^2\right)^{1/8}\right]`$, where $`a`$ and $`b`$ are the distances along the projected major and minor axes, and $`a_\mathrm{e}`$ and $`b_\mathrm{e}`$ are the characteristic scale lengths in those directions. Lambas, Maddox and Loveday \[Lambas et al. 1992\] provide a parameterization for the ellipticities of elliptical galaxies in the APM Bright Galaxy Survey in terms of the parameter $`p=b_\mathrm{e}/a_\mathrm{e}`$. An ellipticity was drawn from each galaxy at random from this distribution, and the position angle of the major axis was also chosen randomly. Songaila et al. \[Songaila et al. 1994\] found that for galaxies with $`18<K<19`$ the median redshift is about 0.6, which means that for $`K\mathrm{}>18`$ (where detection and completeness are being tested) the majority of the galaxies will be distant enough that their angular sizes are relatively insensitive to redshift. High redshift ellipticals have typical characteristic sizes ranging from 2 to 10 kpc (e.g. Dickinson 1997), and so the apparent characteristic size, $`r_\mathrm{e}`$ ($`=(a_\mathrm{e}b_\mathrm{e})^{1/2}`$), of each elliptical was chosen randomly from the range 0.2 to 1.2 arcsec.

The remaining 70% of the galaxies were built using a single exponential profile appropriate for galaxy disks, $`I(r)\mathrm{exp}[r/r_\mathrm{d}]`$, where $`r_\mathrm{d}`$ is the characteristic scale length of the disk. Mao, Mo and White \[Mao et al. 1998\] show that the characteristic scale lengths of disks decreases with redshift approximately as $`(1+z)^1`$, but Simard et al. \[Simard et al. 1999\] shows that a magnitude–size correlation somewhat counter–balances this: at a given apparent magnitude, higher redshift objects must be intrinsically brighter and so are larger. From their data, the typical scale length of galaxy disks at around our completeness limit will be 1–4 kpc. The disk scale length of each galaxy was therefore chosen at random from the range 0.1 to 0.5 arcsec. The disk inclination was chosen at random from 0 to 90 degrees, with inclinations greater than 75 degrees being replaced by 75 degrees to account for the non–zero thickness of the disk. The projected orientation of the galaxy was also chosen randomly.

The model galaxies were then convolved with the stellar PSF derived in step 1.

Step 3: Five stellar or galactic objects were added to each frame with a random position and scaled to a random magnitude in the range $`15<K<22`$. sextractor was then run on the new image using the same input parameters as for the original source extraction, to see if the added objects were detected, and if so to determine the difference between the total magnitude measured by sextractor and the true total magnitude.

Step 4: Step 3 was repeated until 25000 stellar objects and 25000 model galaxies had been added to each image.

From these results, the mean completeness fraction was determined as a function of magnitude for both stars and galaxies over the entire set of images. The results are shown in Figure 1. The 50% completeness limits are $`K=19.75`$ for galaxies and $`K=20.4`$ for stars. As discussed in Section 2.1, new deeper observations are available for five fields, and considering these deeper fields alone, the 50% completeness limits are 20.2 and 20.9 for galaxies and stars respectively.

The input and measured magnitudes of each detected object can also be compared. Typical results for an individual field are shown in Figure 2. Note the generally low scatter, particularly at bright magnitudes, but the occasional source with a large deviation. These are caused by the proximity of the model object to a bright source in the field; large errors such as this are avoided in the observed data by careful examination of the output sextractor catalogues, as discussed above. The mean difference between the input and measured magnitudes, and the scatter of the measured magnitudes around this mean, are shown as functions of input magnitude for both model stars and galaxies in Figure 3.

### 2.3 Multi–colour source extraction

The J–band frames were aligned with the K–band frames by using a number of objects that appeared unresolved on the HST images. This alignment often involved a small shift in position of the frames, but generally much less than a degree of frame rotation. sextractor was then run in its double image mode using the J and K–band images. In double image mode sextractor uses one image to detect objects and define the apertures to be used for flux determination, and then the fluxes and magnitudes are measured from the second image. In this way J–band magnitudes or upper limits were determined for the sources detected in the K–band, through exactly the same apertures.

For the HST data, sextractor was first run on the HST frames obtained at the end of the calibration procedure. Then, the four separate WFPC2 frames were overlaid individually with the K–band data, using between 3 and 15 unresolved objects visible at both wavelengths. In this process the HST data was re-pixelated to match the UKIRT data, essential for a combined running of sextractor on the K–band and HST images; the re-pixelated and aligned HST images were then convolved to the angular resolution of the K–band data using a Gaussian convolving function. sextractor was run in double image mode on the convolved HST data, thus measuring accurate fluxes and magnitudes for the objects through the same apertures as the K–band data. The HST fluxes were corrected for the small differences in gain ratios between the different WFPC2 chips and using a 4% linear correction ramp for the charge transfer efficiency effects \[Holtzman et al. 1995\]. They were also corrected for galactic extinction using the extinction maps of Burstein and Heiles \[Burstein & Heiles 1982\]. Finally, for each object the ‘stellaricity index’ (see Section 2.4), used for the separation of stars and galaxies, was replaced by that calculated in the run of sextractor on the original HST frame: using the highest angular resolution data provides a much more accurate determination of this parameter.

For both the J–band and HST frames, only objects which were detected with fluxes of at least 3$`\sigma `$ were considered, where $`\sigma `$ is the uncertainty on the flux measurement provided by sextractor; this flux error estimate includes the uncertainty due to the Poisson nature of the detected counts and that from the standard deviation of the background counts. An additional source of flux error arises from the uncertainty in the subtraction of the background count level as a function of position across the image. This value was estimated as the product of the area of the extraction aperture and the rms variation of the subtracted background flux across the image. This background subtraction error estimate was combined in quadrature with the flux error given by sextractor to determine the uncertainties on the magnitudes of the extracted objects.

### 2.4 Star–galaxy separation

sextractor provides a ‘stellaricity index’ for each object, which is an indication of the likelihood of an object being a galaxy or a star, based on a neural network technique \[Bertin & Arnouts 1996\]. In the ideal cases a galaxy has a stellaricity index of 0.0 and a star has 1.0. In practice, at faint magnitudes, low signal–to–noise and galaxy sizes smaller than the seeing lead to an overlap in the calculated stellaricity indices for the two types of object.

For sources which were present on the HST frames (typically only about 80% of objects, the precise percentage depending upon the sky rotation angle at which the HST frame was taken), the stellaricity index for each object determined from the un-convolved HST data was adopted: the high resolution of these images meant that, apart from the very reddest objects, star–galaxy separation222Of course, star–galaxy separation using the resolved or unresolved nature of sources means that any quasars will be placed into the star category, but their contribution is negligible. was relatively unambiguous right down to the completeness limit of the K–band frames. For the objects for which no HST data was available, stellaricity indices were taken from the K–band data. At magnitudes $`K\mathrm{}<17.5`$, star–galaxy separation could be carried out directly from these stellaricity values. At fainter magnitudes it was still possible to provide a fairly accurate segregation of galaxies from stars by combining the stellaricity index with the $`JK`$ (or F814W$`K`$) colour of the object, as shown in Figure 4: the bluer average colour of stars can be used to separate out the stars and galaxies in the ambiguous range of stellaricity indices from about 0.5 to 0.9.

Figure 5 shows the total star counts as a function of magnitude derived from all of the frames, fit with the function $`\mathrm{log}N(K)=0.247K0.830`$. This good fit to the differential star counts using a single power–law distribution clearly demonstrates that the star–galaxy separation is working well.

## 3 Galaxy Counts

Using the differences between true and measured magnitudes determined for the simulated galaxies in Section 2.2 (and Figure 3), the measured ‘total’ magnitudes of the galaxies were converted to true ‘total’ magnitudes, and these were binned in 0.5 magnitude bins over the magnitude range $`14<K<20`$ to determine the raw galaxy counts, $`n_{\mathrm{raw}}`$. These were corrected for completeness using the observed mean completeness fraction as a function of magnitude from the simulations (Figure 1). A possible source of systematic error in the number counts at faint magnitudes must also be taken into account: the increased scatter in the photometric measurements at faint magnitudes, and the fact that there are more faint galaxies than bright galaxies, mean that it is more likely for faint galaxies to be brightened above the $`K=20`$ limit (or be brightened and move up one magnitude bin) than for brighter galaxies to appear erroneously in a fainter bin or fall out of the catalogue completely. This leads to an apparent increase in the number counts at the faintest magnitudes. The results from the model galaxy simulations were used to correct for this effect: the required corrections were smaller than 10% in all cases. Combining these effects produced corrected galaxy counts ($`n_\mathrm{c}`$), and scaling by the observed sky area, produced final counts per magnitude per square degree, $`N_\mathrm{c}`$.

Number counts were also derived using the same technique considering just the five fields with deeper images (see Section 2.1). This provided determinations of galaxy counts for $`17.5<K<20.5`$ which may be more reliable at the fainter magnitudes due to smaller completeness corrections.

These galaxy counts are tabulated in Table 1 and plotted in Figure 6. The error on the number counts in each bin was calculated from two factors: the Poissonian error on the raw galaxy counts, and an error in the completeness correction which is generously assumed to be 30% of the number of counts added in the correction. The number counts are compared on the plot with counts from various K–band field surveys in the literature. The final column in Table 1 shows the mean galaxy counts per magnitude bin determined from these literature field surveys.

The galaxy counts derived in this paper show an excess of counts relative to the literature counts, which for magnitudes $`K>15.5`$ is at greater than the $`1\sigma `$ significance level. For comparison, the K–magnitudes of the radio galaxies themselves span the range $`15.4<K<18.0`$. The integrated excess counts over the magnitude range $`15.5<K<20`$ corresponds to, on average, 11 galaxies per field. Unfortunately, as described in Section 2.1, the original goals of these observations was to carry out detailed studies of the radio source host galaxies, and so no blank sky frames were taken. Therefore no comparison can be made between radio galaxy and blank sky frames to determine whether the excess counts found here are associated with structures surrounding the radio galaxies or whether there is some systematic offset between the counts calculated here and the literature counts.

We believe that the majority of the excess counts are associated with the presence of the radio galaxies, for a number of reasons. First, Roche, Eales and Hippelein \[Roche et al. 1998\] found a similar excess in the K–band galaxy counts in a study of the fields surrounding 6C radio galaxies at similar redshifts. Second, as reviewed in the introduction, there are numerous lines of evidence suggesting that there are clusters around at least some of these objects. Third, later in this paper it will be shown that cross–correlation analyses and colour–magnitude relations indicate an overdensity of galaxies comparable to that observed.

## 4 Cross–correlation analyses

### 4.1 Angular cross-correlation estimators

The clustering of galaxies around the radio galaxies can be investigated using the angular cross–correlation function $`w(\theta )`$, which is defined from the probability of finding two sources in areas $`\delta \mathrm{\Omega }_1`$ and $`\delta \mathrm{\Omega }_2`$ separated by a distance $`\theta `$:

$$\delta P=N^2[1+w(\theta )]\delta \mathrm{\Omega }_1\delta \mathrm{\Omega }_2$$

where $`N`$ is the mean surface density of sources on the sky. From this, it follows that for a given survey:

$$DD(\theta )=\frac{1}{2}n_\mathrm{D}(n_\mathrm{D}1)[1+w(\theta )]\frac{\delta \mathrm{\Omega }_\mathrm{D}(\theta )}{\mathrm{\Omega }},$$
(1)
where $`DD(\theta )`$ is the number of data–data pairs with angular separation of between $`\theta `$ and $`\theta +\delta \theta `$, $`n_\mathrm{D}`$ is the total number of sources in data catalogue, $`\mathrm{\Omega }`$ is the total angular area of sky sampled and $`\delta \mathrm{\Omega }_\mathrm{D}(\theta )`$ is the mean angular area of sky accessible at a distance $`\theta `$ to $`\theta +\delta \theta `$ around the data points. Clearly $`\delta \mathrm{\Omega }_\mathrm{D}(\theta )`$ is extremely difficult to calculate due to boundary effects; various estimators for $`w(\theta )`$ have therefore been derived using comparisons between the data points and catalogues of randomly distributed points (e.g. see Cress et al. 1996 for a discussion). For the estimator adopted in this paper, it is considered that if $`n_\mathrm{R}`$ random points are added to the image (where $`n_\mathrm{R}n_\mathrm{D}`$ to minimise errors introduced by the random catalogue) then the number $`DR(\theta )`$ of data–random pairs between $`\theta `$ and $`\theta +\delta \theta `$ will be given by:

$$DR(\theta )=n_\mathrm{D}n_\mathrm{R}\frac{\delta \mathrm{\Omega }_\mathrm{D}(\theta )}{\mathrm{\Omega }}$$
(2)
Combining equations 1 and 2 gives the following estimator for $`w(\theta )`$:

$$w(\theta )=\frac{2n_\mathrm{R}}{(n_\mathrm{D}1)}\frac{DD(\theta )}{DR(\theta )}1$$
(3)
In addition, it is possible to consider the angular cross–correlation function just around an individual source, in this case the radio galaxy. For this, the parameters become $`RGD(\theta )=(n_\mathrm{D}1)[1+w_{\mathrm{rg}}(\theta )]\delta \mathrm{\Omega }_{\mathrm{RG}}(\theta )/\mathrm{\Omega }`$ and $`RGR(\theta )=n_\mathrm{R}\delta \mathrm{\Omega }_{\mathrm{RG}}(\theta )/\mathrm{\Omega }`$, where $`RGD(\theta )`$ and $`RGR(\theta )`$ are respectively the number of radio galaxy–data and radio galaxy–random pairs between $`\theta `$ and $`\theta +\delta \theta `$, $`w_{\mathrm{rg}}(\theta )`$ is the angular cross–correlation function for galaxies around the radio galaxy, and $`\delta \mathrm{\Omega }_{\mathrm{RG}}(\theta )`$ is the angular area of sky accessible at a distance $`\theta `$ to $`\theta +\delta \theta `$ from the radio galaxy. These equations combine to give:

$$w_{\mathrm{rg}}(\theta )=\frac{n_\mathrm{R}}{n_\mathrm{D}}\frac{RGD(\theta )}{RGR(\theta )}1$$
(4)
### 4.2 Calculating $`w(\theta )`$ and $`w_{\mathrm{rg}}(\theta )`$

One problem with estimating $`w(\theta )`$ and $`w_{\mathrm{rg}}(\theta )`$ from the current data is that, due to the K–band dithering technique, the detection limit is a function of position across each image. More objects are detected in the central regions of the image, which would produce a spurious peak in the angular cross–correlation statistics. This non-uniform noise level can be accounted for using the simulations which were discussed in Section 2.2. From these simulations, for each field the completeness fraction can be calculated as a function of both magnitude and position, and then applied to the random catalogue in the following way.

Step 1: Excluding those objects determined to be stars, $`DD(\theta )`$, $`RGD(\theta )`$, and $`n_\mathrm{D}`$ were determined for each frame down to a chosen limiting magnitude $`K_{\mathrm{lim}}`$, for a set of bins in $`\theta `$.

Step 2: 25000 random objects were each assigned a random position on the image, and a magnitude $`KK_{\mathrm{lim}}`$ drawn at random from a distribution matching the observed number counts (Figure 6).

Step 3: For each random object, the completeness fraction of objects of that magnitude and position (averaged in 50 by 50 pixel bins) was determined from the simulations. The object was accepted or rejected at random with a likelihood of acceptance based upon the derived completeness fraction.

Step 4: For those objects which remained in the catalogue, $`DR(\theta )`$, $`RGR(\theta )`$, and $`n_\mathrm{R}`$ were calculated. Hence $`w(\theta )`$ and $`w_{\mathrm{rg}}(\theta )`$ were derived.

This process was repeated for all images and for two different limiting magnitudes, $`K_{\mathrm{lim}}=19`$ and $`K_{\mathrm{lim}}=20`$. $`K_{\mathrm{lim}}=19`$ corresponds to the limit at which the data are still almost 100% complete, whilst $`K_{\mathrm{lim}}=20`$ requires a significant completeness correction but provides improved statistics for the galaxy counts. The results were averaged over all 28 fields; investigating the values for individual fields was not attempted since for these small fields of view the measured amplitudes of single images are too strongly affected by variations in the background (and foreground) counts to provide results with a high statistical significance; only by combining the data of the 28 fields is a sufficiently robust measurement obtained. An estimation for the uncertainty in the value of the combined $`w(\theta )`$ in each bin was provided by the scatter in the values between the different fields. The results are shown in Figure 7.

$`w(\theta )`$ is usually assumed to have a power–law form; $`w(\theta )=A(\theta /\mathrm{deg})^\delta `$. If so, the observed $`w(\theta )`$ will follow a form $`w(\theta )=A(\theta /\mathrm{deg})^\delta C`$, where $`C`$ is known as the integral constraint and arises from the finite size of the field of view. Its value can be estimated by integrating $`w(\theta )`$ over the area of each field, and for our data corresponds to $`C_{\mathrm{gg}}=41.1A_{\mathrm{gg}}`$ and $`C_{\mathrm{rg}}=48.0A_{\mathrm{rg}}`$ for, respectively, all galaxy-galaxy pairs and just radio galaxy-galaxy pairs. It is not possible to determine both the amplitude and the slope of the fit from the current data, and so the canonical value of $`\delta =0.8`$ (which also seems to be appropriate at high redshifts; Giavalisco et al. 1998) has been adopted. The observed data in Figure 7 has been fit with functions of this form, and the resulting angular cross–correlation amplitudes are provided in Table 2.

The results at $`K<19`$ and $`K<20`$ are in approximate agreement. For the galaxy–galaxy pairs there may be a decrease of $`A_{\mathrm{gg}}`$ as fainter magnitude limits are used (as expected, e.g. see Roche et al. 1998), but the errors are too large to determine this with any degree of confidence. The galaxy–galaxy amplitudes are similar to those derived by Roche et al. to the same magnitude limits in the fields of radio galaxies at $`z0.75`$. A more obvious feature is that at both magnitude limits the cross–correlation amplitude around the radio galaxy is significantly larger ($`A_{\mathrm{rg}}A_{\mathrm{gg}}`$). Figure 7 shows that much signal originates at small angular separations ($`\theta \mathrm{}<10^{\prime \prime }`$) and so this cannot be related to any problems with completeness in the outer regions of the frames; the similarity of the results for $`K<19`$ and $`K<20`$ also demonstrates this.

### 4.3 The spatial cross–correlation amplitude

As has been described by many authors (e.g. Longair & Seldner 1979, Prestage and Peacock 1988) it is possible to convert from an angular cross–correlation amplitude to a spatial cross–correlation amplitude if the galaxy luminosity function is known. The spatial cross–correlation function is usually assumed to have a power–law form:

$$\xi (r)=B_{\mathrm{rg}}\left(\frac{r}{\mathrm{Mpc}}\right)^\gamma $$

where the power–law slope is related to the slope of the angular cross–correlation function by $`\gamma =\delta +1`$ ($`=`$1.8 for the $`\delta =0.8`$ adopted here). The spatial cross–correlation amplitude, $`B_{\mathrm{rg}}`$, is then related to that of the angular cross-correlation by $`A_{\mathrm{rg}}=H(z)B_{\mathrm{rg}}`$ (see Longair & Seldner \[Longair & Seldner 1979\] for a full derivation) where:

$$H(z)=\frac{I_\gamma }{N_\mathrm{g}}\left(\frac{D(z)}{1+z}\right)^{3\gamma }\varphi (m_0,z)\left(\pi /180\right)^{(\gamma 1)}$$
(5)
Here, $`I_\gamma `$ is a definite integral which for $`\gamma =1.8`$ has a value of 3.8, $`N_\mathrm{g}`$ is the measured sky density of objects above the magnitude limit (per steradian), $`D(z)`$ is the proper distance to a source at redshift $`z`$, $`\varphi (m_0,z)`$ is the number of galaxies per unit comoving volume which at redshift $`z`$ are more luminous than apparent magnitude $`m_0`$, and the factor $`\left(\pi /180\right)^{(\gamma 1)}`$ is required to convert $`A_{\mathrm{rg}}`$ from degrees to radians. The value of $`H(z)`$ is only a very weak function of $`\gamma `$ provided $`\gamma `$ is of order 2 \[Prestage & Peacock 1988\], and so the fixing of $`\gamma `$ at 1.8 will introduce no significant errors.

$`\varphi (m_0,z)`$ requires assumptions for both the K–correction of galaxies and the local luminosity function, as well as being cosmology dependent, and so this conversion will always be somewhat uncertain. In this paper a Schechter \[Schechter 1976\] form is adopted for the local luminosity function, that is:

$$\varphi (L)\mathrm{d}L=\varphi ^{}\left(\frac{L}{L^{}}\right)^\alpha \mathrm{exp}\left(\frac{L}{L^{}}\right)\mathrm{d}\left(\frac{L}{L^{}}\right)$$

where $`\varphi (L)\mathrm{d}L`$ is the number of galaxies per comoving cubic Mpc with luminosity $`L`$ to $`L+\mathrm{d}L`$, $`\varphi ^{}`$ is the density normalization factor, $`L^{}`$ is the characteristic luminosity and $`\alpha `$ is the faint–end slope. A pure luminosity evolution model is also adopted; that is, values for $`\varphi ^{}`$, $`\alpha `$ and $`L^{}(z=0)`$ are taken, and the only evolution of this function is then evolution of the value of $`L^{}`$ with redshift in accordance with passive evolution predictions. As discussed in the introduction, passive evolution models provide a good description of cluster galaxy properties back to $`z1`$, and they can also provide a reasonable fit to the observed K–band field number counts (see below). These models are, however, undoubtedly a simplification in view of hierarchical galaxy formation theories. In hierarchical models, fewer luminous galaxies are predicted to exist at high redshifts (e.g. see Figure 4 of Kauffmann and Charlot 1998), resulting in the values of $`\varphi (m_0,z)`$ and consequently $`H(z)`$ being lower, and hence $`B_{\mathrm{rg}}`$ will be increased. The pure passive luminosity approach therefore provides a conservative lower estimate for the spatial cross–correlation amplitude around high redshift radio galaxies.

For the passive evolutionary models, the galaxy population was split into four different galaxy types: ellipticals and S0’s (E’s), Sa and Sb types (Sab’s), Sc types (Sc’s) and Sd types and irregulars (Sdm’s). Galaxies of these types were built up using the Bruzual and Charlot \[Bruzual & Charlot 1993\] stellar synthesis codes (1996 version), assuming a Scalo initial mass function, solar metallicity, a formation redshift $`z_\mathrm{f}=10`$, and four different star formation histories (cf. Guiderdoni and Rocca–Volmerange 1990). The E’s were assumed to form their stars in a rapid early burst with the star formation decreasing exponentially on a $`\tau _{\mathrm{sfr}}=0.5`$ Gyr timescale. The Sab’s had 50% of their stars in a bulge component formed in the same manner as the E’s, and the remaining 50% in a disk–like component with a much longer star formation timescale ($`\tau _{\mathrm{sfr}}=6`$ Gyr). The Sc’s were modelled using just a single long star formation timescale ($`\tau _{\mathrm{sfr}}=8`$ Gyr), and the Sdm’s using a still longer timescale ($`\tau _{\mathrm{sfr}}=50`$ Gyr) to produce significant star formation at the current epoch. These star formation histories approximately reproduce the colours of the various morphological types at the current epoch. It was assumed that all of these four morphological types displayed the same luminosity function (that is, $`M^{}`$ and $`\alpha `$), with the weightings of the different types (ie. the relative contributions to $`\varphi ^{}`$) being 30%, 30%, 30% and 10% respectively. In fact, at wavelengths as long as the K–band, the K and evolutionary corrections are small and relatively independent of morphological type, and so the results are not strongly dependent upon any of these assumptions. Changing the values adopted above by a sufficiently large amount that they produce unacceptable distributions of galaxy colours at the current epoch produces only 10-20% changes in the final K–band luminosity function.

Local K–band luminosity functions have been derived by Gardner et al. \[Gardner et al. 1993\], by Mobasher, Sharples and Ellis \[Mobasher et al. 1993\] and by Loveday \[Loveday 2000\]. All three of these datasets are approximately consistent with a characteristic absolute magnitude $`M_\mathrm{K}^{}=25.1`$ and a slope $`\alpha =1.0`$. To investigate whether the combination of this luminosity function and the assumption of passive evolution models produces acceptable results, and to calculate the appropriate value of $`\varphi ^{}`$, the K–band number counts as a function of apparent magnitude that would be expected in this model were derived. The value of $`\varphi ^{}`$ was adjusted to provide a good fit to the observed data points. The resulting fit, for a value $`\varphi ^{}0.004`$, is shown as the solid line in Figure 6, and demonstrates that this simple model can provide a reasonable fit to the observed number counts.

With this model luminosity function, it is possible to calculate $`\varphi (m_0,z)`$ and, using the values of $`N_\mathrm{g}`$ from Table 1 ($`N_\mathrm{g}=9.2\times 10^7\mathrm{sr}^1`$ for $`K<19`$ and $`2.0\times 10^8\mathrm{sr}^1`$ for $`K<20`$), hence to determine $`H(z)`$. For $`K<19`$, $`H(z)`$ ranges from $`15.4\times 10^7`$ to $`1.8\times 10^7`$ over the redshift range 0.6 to 1.8. The mean value of $`1/H(z)`$ averaged over all 28 radio source redshifts is $`\overline{1/H(z)}1.7\times 10^6`$, and using this value to convert from $`A_{\mathrm{rg}}`$ to $`B_{\mathrm{rg}}`$ gives $`B_{\mathrm{rg}}=600\pm 240`$. For $`K<20`$ the corresponding range of $`H(z)`$ is from $`12.2\times 10^7`$ to $`3.1\times 10^7`$, with a mean value of $`\overline{1/H(z)}1.4\times 10^6`$, corresponding to $`B_{\mathrm{rg}}=510\pm 120`$.

These values can be interpreted physically by comparing with the equivalent values for Abell clusters calculated between the central galaxy and the surrounding galaxies ($`B_{\mathrm{cg}}`$). This has been calculated independently as a function of Abell cluster richness by a number of authors \[Prestage & Peacock 1988, Hill & Lilly 1991, Andersen & Owen 1994, Yee & López–Cruz 1999\]. Converting their values to $`H_0=50`$ km s-1 Mpc-1 and $`\gamma =1.8`$ they are all approximately consistent with each other, and average to $`B_{\mathrm{cg}}350`$ for Abell class 0 and $`B_{\mathrm{cg}}710`$ for Abell class 1. The environments surrounding the redshift one radio galaxies are therefore comparable, on average, to those of clusters of between Abell Classes 0 and 1 richness.

Hill and Lilly \[Hill & Lilly 1991\] further showed that there is a correlation between the value of $`B_{\mathrm{cg}}`$ and the parameter $`N_{0.5}`$, where $`N_{0.5}`$ is an Abell–type measurement defined as the net excess number of galaxies within a radius of 0.5 Mpc of the central galaxy with magnitudes between $`m_1`$ and $`m_1+3`$, $`m_1`$ being the magnitude of the central galaxy. Using a larger dataset, Wold et al. \[Wold et al. 2000\] calibrated this relation as $`B_{\mathrm{cg}}=(37.8\pm 10.9)N_{0.5}`$. The average $`B_{\mathrm{rg}}`$ value for $`K<20`$ then implies that an average net excess of 13.5 galaxies should be found around each radio galaxy within a radius of 0.5 Mpc and with a magnitude down to three magnitudes fainter than the radio galaxy. The data presented in this paper cover about 80% of this sky area and, since the radio galaxies have typical magnitudes of $`K17`$, the galaxy counts to $`K=20`$ do sample approximately 3 magnitudes below the radio galaxy. Therefore, the net excess counts of 11 galaxies per field down to $`K=20`$ (see Section 3) is fully consistent with the determined value of $`B_{\mathrm{rg}}`$.

## 5 Colour–magnitude and colour–colour relations

### 5.1 Near–infrared colour–magnitude relations

For each field the magnitudes of all galaxies were determined through all of the different filters, as described in Section 2.3. After correction for galactic reddening, these were used to construct colour–magnitude relations for each field. The near–infrared $`JK`$ vs $`K`$ relation for each field (if $`J`$–band data were not taken, the longest wavelength HST filter available was used instead) are shown in Figure 9 in order of ascending redshift. The stars (open diamonds) and galaxies (asterisks) have been separated on these diagrams using the technique described in Section 2.4, and the radio galaxy is plotted using a large symbol. The uncertainties on the measured colours are not shown to avoid cluttering up the diagrams; at bright magnitudes ($`K\mathrm{}<17`$) they are dominated by calibration uncertainties ($`\mathrm{}<0.1`$ mags), but increase to about 0.3 magnitudes by $`K20`$.

The colour–magnitude relation for the Coma cluster has been redshifted and evolved according to the passive evolutionary models for elliptical galaxies described in the previous section; this redshifted relation is over-plotted on the figures. Note that the radio galaxies generally have near–infrared colours similar to this theoretical line, indicating that these are old elliptical galaxies (e.g. Best et al. 1998b), but not in all cases: some of the radio galaxies appear redder because a heavily reddened nuclear component contributes to the K–band emission (e.g. 3C22, 3C41; Economou et al. 1995, Rawlings et al. 1995, Best et al. 1998), whilst for others the bluer filter is significantly contaminated by excess optical–UV emission induced by the radio source (the ‘alignment effect’, e.g. McCarthy et al. 1987).

A number of features are apparent from these near–infrared colour–magnitude diagrams. Considering initially those radio galaxies with redshifts $`z\mathrm{}<0.9`$, many of the sample show reasonably convincing evidence for associated clusters; here it should be born in mind that for an Abell Class 1 cluster only 10 to 15 associated cluster galaxies are expected in the observed region of sky down to three magnitudes below the magnitude of the radio galaxy (see Section 4.3). Some fields (e.g. 3C34, 3C337) clearly show at least this number (the ‘background’ counts at these K–magnitudes and colours are small, as can be seen by a comparison with the colour–magnitude relations of the higher redshift sources), whilst other fields (e.g. 3C217) show few if any associated galaxies. There are clearly large source–to–source variations in environmental richness.

Where a clear colour–magnitude relation is observed, the mean colour of this relation lies close to that calculated theoretically by just passively evolving the Coma colour–magnitude relations back in redshift, and the scatter around the colour–magnitude sequence is small. These results have been shown before for optical and X–ray selected clusters at redshifts out to $`z0.8`$ (e.g. Stanford et al. 1998). There are some small deviations in the colour of the observed sequence from the passive evolution relation (e.g. see 3C265), but none much larger than a couple of tenths of a magnitude. The near–infrared colours of the galaxies in the radio galaxy fields are therefore consistent with those observed in other clusters at the same redshift, implying that the excess galaxy counts are associated with a structure at the radio galaxy redshift. This is important because of previous suggestions that powerful distant radio galaxies may be systematically amplified by foreground lensing structures (Le Fèvre et al. 1987; Benítez, Martínez–Gonzáles and Martin–Mirones 1997); were such structures to be present, they could account for both the excess K–band counts and the peak in the cross–correlation statistics, but the colour–magnitude relations argue against this.

At higher redshifts the evidence for clear colour–magnitude relations is much poorer. This is mainly because the combination of the greatly increased contribution from field galaxies at these fainter magnitudes and the increased scatter in the colours due to photometric uncertainties, results in any colour–magnitude sequence appearing much less prominent. The difficulty of selecting cluster candidates at these redshifts on the basis of a single colour can be gauged by examining 3C324; for this field, Dickinson \[Dickinson 1997\] has confirmed the presence of a poor cluster of galaxies, but this is barely apparent from its colour–magnitude relation.

Although no prominent colour–magnitude sequences are seen in the higher redshift fields, there remains a net excess of K–band counts: if the fields of the radio galaxies are divided into two redshift bins, no significant differences in the faint galaxy counts are seen between the high and low redshift fields. A simple analysis also shows that the excess in the fields of the higher redshift radio galaxies is comprised of red ($`JK1.75`$) galaxies. In Table 3 are given the mean number of galaxies per field with magnitudes $`17<K<20`$ and colours $`JK1.75`$ or $`1.25<JK<1.75`$, in the low and high redshift bins. There are more galaxies with bluer $`JK`$ colours in the fields of the lower redshift radio galaxies than in those at higher redshifts, due to the associated cluster galaxies in the lower redshift fields which have colours $`1.25\mathrm{}<JK\mathrm{}<1.75`$ (cf. Figure 8). On the other hand, there are more galaxies with $`JK1.75`$ colours in the high redshift than low redshift fields. The excess K–band galaxy counts in the high redshift field appear to be predominantly associated with red galaxies, with colours similar to those expected for old cluster galaxies at these redshifts. This again indicates that the excess number counts are associated with a structure at the redshift of the radio galaxy rather than a foreground structure.

### 5.2 Multi–colour relations

In Figure 9 are plotted a complete set of colour–magnitude relations for six galaxies in the sample, chosen to be those which show amongst the best examples of near–infrared colour–magnitude relations for their redshift333Equivalent plots for the other fields are available from the author on request, but are not reproduced here to save space.. From these it is apparent that, although the scatter around the near–infrared colour–magnitude sequence remains small even at these high redshifts, colours that reach shortward of the rest–frame 4000Å break show a dramatic increase in the scatter of the relation. (e.g. compare the various relations for 3C34); these colours can be strongly influenced by small amounts of recent or on–going star formation, indicating that this may be common in these high redshift clusters. Whether such star formation is in some way connected with the presence of a powerful radio source in these clusters cannot be distinguished from these data.

To properly investigate the nature of the galaxies in these fields, all of the colour information must be used simultaneously to derive photometric redshifts and investigate star formation activity. This is beyond the scope of the current paper but will be addressed later (Kodama & Best, in preparation). Here, the use of the all multi–colour information simulataneously is merely demonstrated in Figure 10 through colour–colour plots for these six fields. For each field, the near–infrared $`JK`$ colour of each galaxy is plotted against its optical-infrared colour. These data are then compared against theoretical evolutionary tracks for the four different passively evolving galaxy models considered in Section 4.3 (E’s, Sab’s, Sc’s, Sdm’s).

For the lower redshift radio galaxies there is clearly a large concentration of galaxies with colours very close to those of the model elliptical galaxy at the redshift of the radio source. Further, there is a distribution of galaxies with colours between this and the colours of the model spiral galaxies at that redshift. In contrast, for the higher redshift radio galaxies, no strong concentration of galaxies is seen close to the elliptical galaxy prediction, and the number of cluster candidates lying between the locations of the elliptical and spiral model galaxies is smaller than that found in the lower redshift cases. Clearly, despite the excess $`K`$–counts and red $`JK`$ galaxies in the fields of these high redshift radio sources, to accurately investigate cluster membership requires the construction of photometric redshifts using several colours measured with high photometric accuracy.

## 6 Discussion

In this paper a number of pieces of evidence for clustering around distant 3CR radio galaxies have been presented. These can be summarised as follows:

* The K–band number counts show an overdensity of faint galaxies in the fields of the radio galaxies, with a mean value of 11 excess galaxies per field.
* This excess is comparable to the galaxy overdensity expected for a field of view of this size centred on a cluster of approximately Abell Class 0 richness.
* Cross–correlation analyses show a pronounced peak in the angular cross–correlation function around the radio galaxies.
* Assuming that the galaxy luminosity function undergoes pure passive luminosity evolution with redshift, the corresponding spatial cross–correlation amplitude lies between those determined for Abell Class 0 and Abell Class 1 clusters.
* The galaxies in the fields of most of the lower redshift radio galaxies in the sample show clear near–infrared colour–magnitude relations with only small scatter. The colours of these sequences are in agreement with those of other clusters at these redshifts, indicating that the excess number counts and cross–correlation peak are both associated with a structure at the redshift of the radio galaxy.
* There is considerably more scatter in the relations involving shorter wavelength colours, suggesting low levels of recent or on–going star formation in many of the galaxies.
* At higher redshifts the colour–magnitude relations are less prominent due to increased background contributions, but there is a clear excess of galaxies with very red infrared colours.

These features all provide strong evidence that distant radio galaxies tend to reside in rich environments. The number counts, the cross–correlation statistics, and the colour–magnitude relations all complement the previous results from X–ray imaging, narrow–band imaging, spectroscopic studies and radio polarisation studies discussed in Section 1. A coherent picture that most, but not all, high redshift radio sources live in poor to medium richness clusters has now been built.

Taking the results at face value, the environmental richness around these $`z1`$ radio galaxies is higher than that around powerful radio galaxies at $`z0.5`$ calculated by Yates et al. \[Yates et al. 1989\] and by Hill and Lilly \[Hill & Lilly 1991\]. This suggests that the increase between $`z=0`$ and $`z=0.5`$ in the mean richness of the environments surrounding FR II radio galaxies, found by those authors, continues to higher redshifts.

Some notes of caution must be added to this conclusion. First, from the variations in richness of the colour–magnitude relation at any given redshift, (e.g. compare 3C217 and 3C226), it is apparent that there is a wide spread in the density of the environments in these fields. Although most show some evidence of living in at least group environments, not all powerful distant radio galaxies lie in clusters. Second, a simple visual comparison with the extremely rich high redshift cluster MS1054$``$03 ($`z=0.83`$, van Dokkum 1999) is sufficient to demonstrate that even at high redshifts, powerful radio galaxies still avoid the most extreme richness clusters.

Further, although the galaxy count excesses and the cross–correlation amplitudes have been compared with those of Abell clusters at low redshifts, in hierarchical galaxy formation models such comparisons will always be somewhat ambiguous. On–going mergers between galaxies mean that more sub–clumps are seen at higher redshifts, whilst the general galaxy cross–correlation length also evolves with redshift in a manner dependent upon both cosmological parameters and the method of selecting the galaxy populations \[Kauffmann et al. 1999\]. Therefore, quantitative interpretations of either of these parameters at high redshift must be considered with some care. On the other hand, in hierarchical growth models, the structures in which the radio galaxies reside will also continue to grow and evolve into much richer structures by a redshift of zero, meaning that the qualitative result that the high redshift radio galaxies lie in rich environments for their redshift is secure.

Finally, it is important to consider the consequences for our understanding of the onset and nature of powerful radio sources of a change in the environments of FR II radio galaxies from groups at low redshifts to clusters at high redshift. As discussed by Best et al. \[Best et al. 1998b\], if this result holds then the standard interpretation of the tightness and slope of the Hubble $`Kz`$ relation, that of ‘closed–box passive evolution’ of radio galaxies at $`z\mathrm{}>1`$ into radio galaxies at $`z0`$, is no longer valid. It is not possible that the environments can become less rich with progressing cosmic time. Instead, Best et al. propose that powerful radio galaxies selected at high and low redshift have different evolutionary histories but must contain a similar mass of stars, a few times $`10^{11}M_{}`$, so conspiring to produce the observed ‘passively evolving’ K$`z`$ relation. In their model, powerful FR II radio sources are seldom formed in more massive galaxies (that is, in central cluster galaxies at low redshifts) because of the difficulties in supplying sufficient fueling gas to the black hole: in rich low redshift clusters the galaxies and gas have been virialised and take up equilibrium distributions within the cluster gravitational potential, the galaxies have high velocity dispersions greatly reducing the merger efficiency, and there is a dearth of gas–rich galaxies close to the centre of the clusters which might merge with, and fuel, the central galaxy. Thus, the formation of a powerful radio source in these environments is a rare event (but can still happen, e.g. Cygnus A).

At high redshifts, radio galaxies can be found in (proto) cluster environments because these are not yet virialised, have frequent galaxy mergers, and have a plentiful supply of disturbed intracluster gas to fuel the central engine and confine the radio lobes. The central cluster galaxies will be amongst the most massive galaxies at these redshifts and so, from the correlation between black hole mass and bulge mass (e.g. Kormendy and Richstone 1995), will have the most massive black holes. The kinetic energy of the relativistic radio jets of distant 3CR radio galaxies corresponds to the Eddington luminosity of a black hole with $`M10^810^9M_{}`$ \[Rawlings & Saunders 1991\], implying that these sources are fueled close to the Eddington limit. Therefore the most powerful radio sources will tend to be powered by the most massive central engines and hence be hosted by the most massive galaxies, which tend to be found at the centres of forming rich clusters. The significant scatter in the black hole versus bulge mass correlation \[Kormendy & Richstone 1995\] would, however, result in some scatter in the richness of the environments of the radio galaxies. The data presented in this paper are in full agreement with this model.

## Acknowledgements

This work was supported in part by the Formation and Evolution of Galaxies network set up by the European Commission under contract ERB FMRX– CT96–086 of its TMR programme. The United Kingdom Infrared Telescope is operated by the Joint Astronomy Centre on behalf of the U.K. Particle Physics and Astronomy Research Council. This work is, in part, based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA Inc., under contract from NASA. I thank Huub Röttgering for useful discussions, and the referee for helpful comments.