text:
# The spatial clustering of radio sources in NVSS and FIRST; implications for galaxy clustering evolution

## 1 Introduction

In striking contrast with the extremely high level of isotropy observed in the temperature of the cosmic microwave background (see e.g. de Bernardis et al. 2000), galaxies are not distributed throughout the Universe in a random manner. According to the gravitational theory of instability the present structures originated from tiny fluctuations in the initial mass density field. This has shaped the large-scale structure of the Universe, which consists of vast empty regions (voids), and strings of dark and luminous matter (walls) where billions of galaxies are found.

The clustering properties of galaxies can be quantified using statistical techniques, such as methods of nearest neighbours, counts in cells, power spectra, and correlation functions (see Peebles 1980, for an in-depth mathematical review). In particular the two-point correlation function is a simple, but powerful tool that has become a standard for studying large-scale structure. The clustering of cosmological objects can be characterized by their spatial correlation function, which has the form $`\xi (r)=(r/r_0)^\gamma `$ where $`r_0`$ is the present-day correlation length and $`\gamma 1.8`$ for objects ranging from clusters to normal galaxies (see Bahcall & Soneira 1983, for a review). The local population of galaxies is a relatively unbiased tracer of the underlying matter distribution, with $`r_0=5.4`$ $`h^1`$ Mpc derived from galaxies in the early CfA redshift survey by Davis & Peebles (1983), however more recent low-redshift surveys show that the clustering of galaxies depends strongly on luminosity and/or morphological type. For example, local $`LL_{}`$ ellipticals represent spatial structures that are much more strongly clustered with $`r_0712`$ $`h^1`$ Mpc (e.g. Guzzo et al. 1997; Willmer et al. 1998; Norberg et al. 2002). From deep, magnitude-limited redshift samples it has been found that the comoving correlation length of galaxies declines with redshift, roughly as expected from simple gravitational theory (e.g. CFRS, Le Fèvre et al. 1996; Hawaii K, Carlberg et al. 1997; CNOC2, Carlberg et al. 2000; CFDF, McCracken et al. 2001). In contrast to this, the clustering strength of quasars appears to vary little over $`0z2.5`$. Croom et al. (2001) found an approximately constant amplitude of $`5`$ $`h^1`$ Mpc from $`10,000`$ quasars in the 2dF QSO Redshift Survey. Likewise, Daddi et al. (2001, 2002) found that the (comoving) correlation length of massive elliptical galaxies also shows little evolution with redshift. They find $`r_0=12\pm 3`$ $`h^1`$ Mpc for a population of extremely red objects (EROs) at $`z1`$ (see also McCarthy et al. 2001; Roche et al. 2002; Firth et al. 2002), which are consistent with being the passively evolving progenitors of local massive ellipticals (e.g. Dunlop et al. 1996; Cimatti et al. 1998, 2002; Dey et al. 1999; Liu et al. 2000). Color selection methods such as Lyman-break (Steidel et al. 1995) and narrow-band imaging techniques are providing statistical samples of very high redshift galaxies, allowing us to study large-scale structure at even earlier epochs. Lyman-break galaxies have correlation lengths as high as $`r_03`$ $`h^1`$ Mpc even at $`z34`$, and are thought to be associated with (mildly) biased star-forming galaxies (e.g. Adelberger 2000; Ouchi et al. 2001; Porciani & Giavalisco 2002).

Studying clustering as a function of redshift and galaxy type may provide important constraints on some long-standing problems in cosmology concerning galaxy formation and evolution. For example, which of the galaxies observed at high redshift are the progenitors of local galaxy populations, and which of the local galaxies host the remnant black holes that once powered high redshift active galactic nuclei (AGN)? Two common views on how structures observed at high redshifts may be related to structures observed today are represented by (i) the galaxy conservation model (e.g. Fry 1996; Tegmark & Peebles 1998) in which it is assumed that galaxies formed very early in a monolithic collapse (e.g. Eggen et al. 1962) and have evolved passively with a decreasing star formation rate since $`z2`$, and (ii) the hierarchical merging model (e.g. Mo & White 1996) in which it is assumed that the most luminous galaxies formed more recently in massive dark matter haloes that have grown hierarchically by the merging of less massive galaxies and their haloes. Kauffmann & Charlot (1998) computed the evolution of the observed K-band luminosity function for both the monolithic case and the hierarchical case, and found that by a redshift of $`1`$ these models differ greatly in the abundance of bright galaxies they predict. Likewise, the validity of these models can be tested by comparing predictions for galaxy clustering from numerical simulations or (semi-)analythic theory (e.g. Kauffmann et al. 1999b; Moustakas & Somerville 2002; Mo & White 2002, and references therein) with the observed clustering of a population of galaxies. In the case of pure monolithic collapse galaxy clustering is dictated by the evolution of galaxy bias under the rules of gravitational perturbation theory, but without the extra non-linear effects arising from galaxy mergers. Such a scenario can be thought of as a baseline model for the clustering of the matter as probed by galaxies situated in average mass haloes. However, in the hierarchical case the evolution of galaxy bias is much more complex, since galaxies are no longer conserved quantitities (Kauffmann et al. 1999b). Comparing their observations to model predictions Daddi et al. (2001) find that such a scenario best explains the clustering evolution of massive ellipticals out to $`z1`$.

Radio surveys can make an important contribution to this study: the use of magnitude-limited surveys for finding high redshift objects is usually a cumbersome task, while any flux density limited sample of radio sources contains objects at redshifts of $`z05`$ (Dunlop & Peacock 1990). Powerful extra-galactic radio sources, or AGN in general, result from the fuelling of a supermassive blackhole (e.g. Rees 1984, 1990), and there is evidence that the host galaxies of these high-redshift AGN are associated with some of the most massive structures in the early Universe (e.g. McCarthy 1988; Crawford & Fabian 1996; Röttgering et al. 1996; Best et al. 1998; Pentericci et al. 1999; Venemans et al. 2002). Moreover, because powerful AGN were far more numerous at $`z12`$ than today, radio surveys can be used to probe a population of massive galaxies in the epoch of galaxy formation.

Despite initial concerns that any cosmological clustering of radio sources may be undetectable due to the relatively broad redshift distribution washing out the signal (e.g. Webster & Pearson 1977; Griffith 1993), Kooiman et al. (1995) detected strong clustering of bright radio sources in the 4.85 GHz 87GB survey. Cress et al. (1996) made a thorough analysis of clustering at the mJy-level. Using the 1.4 GHz FIRST survey (see also Magliocchetti et al. 1998) they obtained the first high-significance measurement of clustering from a deep radio sample, allowing them to investigate the separate contributions of both AGN and starburst galaxies (but see Wilman et al. 2003). Further results on the statistics of radio source clustering have been presented by Loan et al. (1997) and Rengelink (1998), who based their analysis on the 4.85 GHz Parkes-MIT-NRAO survey and the 325 MHz WENNS survey, respectively. In high-resolution surveys such as FIRST, large radio sources can become resolved in several components, thereby spuriously contributing to the cosmological clustering signal. Cress et al. (1996) and Magliocchetti et al. (1998) outlined the basic steps involved in separating the signal due to this effect from the true cosmological clustering, although the angular size distribution of radio sources at the mJy level is still largely unconstrained.

Since the individual redshifts of the radio sources are generally not known, one usually only measures the two-dimensional clustering by means of the angular correlation function, $`w(\theta )`$. However, the redshift distribution of the survey can be used to constrain $`r_0`$. Using this so-called Limber inversion technique (Limber 1953; Rubin 1954; Phillipps et al. 1978; Peebles 1980), radio sources from the above surveys are typically found to have $`r_0515`$ $`h^1`$ Mpc. Rengelink (1998) and Rengelink & Röttgering (1999) pointed out that this broad range in $`r_0`$ measured can be explained by a scenario in which powerful radio sources have a larger $`r_0`$ than less powerful radio sources. This would be highly consistent with the mounting evidence that powerful radio galaxies are the high-redshift progenitors of local cD-galaxies residing in massive environments that are hence strongly clustered. Here, we will further explore the hypothesis of Rengelink et al. by investigating the clustering of radio sources in a number of flux-limited subsamples taken from the 1.4 GHz NRAO VLA Sky Survey (see also Blake & Wall 2002a, b), the largest existing 1.4 GHz survey to date, containing $`1.8\times 10^6`$ radio sources down to a flux density limit of $`2.5`$ mJy at $`45\mathrm{}`$ (FWHM) resolution (Condon et al. 1998). We also present new results on clustering using the latest release of the FIRST survey, carefully taking into account the contribution of multiple-component radio sources, which we found to be severely underestimated in earlier analyses.

The outline of this article is as follows: in Sect. 2 we describe our methods for measuring the angular two-point correlation function. In Sect. 3 we describe the NVSS and FIRST radio surveys, and in Sect. 4 we present measurements of the angular clustering of the sources in these surveys and construct a simple model of the angular size distribution of radio sources. We derive an estimate of $`r_0`$ as a function of flux density limit in Sect. 5. In Sect. 6 we compare our results with the results found for other populations of galaxies taken from literature, and discuss how the combined measurements relate to current theories on galaxy formation and evolution. The main conclusions are summarized in Sect. 7.

## 2 The angular correlation function

The galaxy angular two-point correlation function, $`w(\theta )`$, is defined as the excess probability, over that expected for a Poissonian distribution, of finding a galaxy at an angular distance $`\theta `$ from a given other galaxy (e.g. Peebles 1980):

$$\delta P=n[1+w(\theta )]\delta \mathrm{\Omega },$$
(1)
where $`\delta P`$ is the probablility, $`n`$ is the mean surface density and $`\delta \mathrm{\Omega }`$ a surface area element. The angular two-point correlation function of a given sample of objects can be estimated as follows. For each object, determine the angular distances to all other objects, then count the number of objects in each angular distance interval, denoted by $`DD(\theta )`$. As we want to calculate the excess probability of finding a galaxy at a certain distance from another galaxy due to clustering, we compare the observed distribution, $`DD(\theta )`$, with the expected distribution of distances, $`RR(\theta )`$, calculated from large artificial catalogues of randomly placed sources. We note that several variants of $`w(\theta )`$-estimators exist in literature, of which the methods proposed by Hamilton (1993) and that of Landy & Szalay (1993) (see Blake & Wall 2002b, for application of this estimator to NVSS) are generally considered to be the most robust. We follow Rengelink (1998) and Wilman et al. (2003) and use the Hamilton estimator

$$w(\theta )=\frac{4n_Dn_R}{(n_D1)(n_R1)}\frac{DD(\theta )RR(\theta )}{DR(\theta )DR(\theta )}1,$$
(2)
where $`n_D`$ and $`n_R`$ are the number of sources in the data and random catalogues, respectively, and the numerical factor $`4n_Dn_R/(n_D1)(n_R1)`$ normalizes the pair counts. This estimator additionally makes use of the cross-correlation between data and random catalogues, $`DR(\theta )`$, to minimize effects due to large-scale fluctuations in the mean galaxy density. We estimate $`w(\theta )`$ by averaging over the $`w(\theta )`$ computed using 16 different random catalogues, each containing the same number of sources as the data catalogue to minimize the errors in $`DR(\theta )`$ and $`RR(\theta )`$ (a similar result can be obtained by constructing a single random catalogue that vastly exceeds the size of the data catalogue). Poissonian errors on the binned values of $`w(\theta )`$ are estimated by $`\delta w(\theta )=\sqrt{[1+w(\theta )]/DD(\theta )}`$. Alternatively, errors can be computed using the so-called bootstrap resampling method of Ling et al. (1986). In this method, the standard deviation in $`w(\theta )`$ found among a large number of pseudo-random resamples of the original dataset is used as a measure of the error in $`w(\theta )`$. However, we found that fitting a model to $`w(\theta )`$ (see Sect. 4) using (i) Poissonian errors, and (ii) bootstrap errors gives results that are consistent within the errors of the fitted parameters. Therefore, given the unprecedented volumes of the radio surveys we use the first method instead of the relatively expensive bootstrap technique.

## 3 Survey descriptions and data selection

### 3.1 The NRAO VLA Sky Survey

The NRAO VLA Sky Survey (NVSS) is the largest radio survey that currently exists at 1.4 GHz. It was constructed between 1993 and 1998 (Condon et al. 1998), and covers $`10.3`$ sr of the sky north of $`\delta =40\mathrm{°}`$ ($`82\%`$ of the sky). Fig. 1 indicates the coverage of the NVSS. With a limiting flux density of $`2.5`$ mJy ($`5\sigma _{rms}`$) and an angular resolution of $`45\mathrm{}`$ (FWHM), the NVSS contains about $`1.8\times 10^6`$ sources, and is considered to be $`99\%`$ complete at a flux density limit of 3.4 mJy (Condon et al. 1998). The NVSS is based on 217,446 snapshot observations (of mostly 23 seconds) using the VLA in D- and DnC-configuration. These snapshots were then combined to produce a set of $`4\mathrm{°}\times 4\mathrm{°}`$ datacubes containing Stokes I, Q, and U images. A source catalogue was extracted by fitting the images with multiple elliptical Gaussians. Since the angular resolution of the NVSS ($`\theta 45\mathrm{}`$ FWHM) is well above the median angular size of extra-galactic radio sources ($`\theta 10`$ arcseconds), most sources in the catalogue are unresolved ($`95\%`$ for $`3<S_{1.4}<10`$ mJy). The main NVSS data products have been made publicly available for the use of the astronomical community, and can be obtained from the NRAO website111http://www.cv.nrao.edu/nvss/.

### 3.2 NVSS data selection

To optimize our catalogue for measuring the true cosmological clustering of radio sources, we have carried out a detailed examination of the NVSS source catalogue to identify and correct regions that may spuriously contribute to $`w(\theta )`$ :
(i) The edge of the survey just a few arcminutes south of $`\delta =40\mathrm{°}`$ follows an irregular pattern with right ascension. We select the region $`\delta 40\mathrm{°}`$ to ensure that the boundary of the survey is straight.
(ii) The survey area is known to contain six hexagonal gaps due to missing snapshot observations that we masked from the catalogue by excluding rectangular regions of $`2\mathrm{°}\times 2\mathrm{°}`$ fully covering each gap. The regions are listed in Table 1.
(iii) We constructed a map of the NVSS source density as a function of position on the sky by applying an equal-area projection to the catalogue and plotting filled contours of the number of objects in $`1\mathrm{°}\times 1\mathrm{°}`$ non-overlapping cells covering the survey area. This map is shown in Fig. 1. The scaling of the greyscale was chosen so that underdense regions of $`2\sigma `$ below the mean density are black, and overdense regions of $`2\sigma `$ above the mean are white. Radio emission from the region of the galactic plane, as evidenced by a continuous chain of large white areas in Fig. 1, is dominated by the large population of galactic radio sources that consists mostly of supernova remnants and H$`\mathrm{II}`$ regions. In Fig. 4 we plot the rms-noise level as a function of galactic latitude, where the rms-noise level in each latitude bin is the average of the locally determined rms-noise values listed for every source entry in the NVSS catalogue. The rms-noise level is found to peak at $`b=0\mathrm{°}`$ due to the overcrowding of galactic sources, but falls off to a relatively constant level of $`0.48`$ mJy beam-1 for $`|b|10\mathrm{°}`$. We decided to exclude the region of the galactic plane that is bounded by $`|b|=10\mathrm{°}`$, which was chosen so that the large overdense regions in Fig. 1 are all fully masked and the rms-noise is at a relatively constant level.
(iv) Further inspection of Fig. 1 reveals that some regions are associated with a significant increase in the local source density. From contour maps of these areas it was found that bright and/or extended sources are sometimes accompanied by significant numbers of spurious sources due to a side-effect of the fitting algorithm used to extract the sources, and, in some cases, due to side-lobe contamination. From the catalogue we excluded rectangular regions of mostly $`1\mathrm{°}\times 1\mathrm{°}`$ in size centered on each of these sources (larger regions of up to $`2\mathrm{°}\times 2\mathrm{°}`$ were required in some cases). The excluded regions are listed in Table 1. No regions of $`2\sigma `$ underdensities were found.

(v) Most of the NVSS observations were conducted using the VLA in D-configuration, but the regions $`\delta 10\mathrm{°}`$ and $`\delta +78\mathrm{°}`$ were observed using the hybrid DnC-configuration to counterbalance projection effects which result from foreshortening of the north-south $`uv`$-coverage range. Fig. 2 shows the NVSS source density as a function of declination for various flux-limited sub-samples. Below the flux density limit of 10 mJy, the use of the DnC-configuration has caused a significant decrease in sensitivity leading to a drop in the source density of $`10`$ percent (see also Fig. 1). As this will inevitably cause spurious signal in the angular two-point correlation function, we selected only the regions observed in D-configuration for measuring $`w(\theta )`$ below flux density limits of 10 mJy.

Table 2 lists the final regions and the number of sources in them for various flux density limited subsamples.

### 3.3 The FIRST Survey

The FIRST (Faint Images of the Radio Sky at Twenty centimeters) survey (Becker et al. 1995) is another 1.4 GHz VLA survey, which was started in 1993 and is still under construction. Using the VLA in B-configuration it will ultimately cover $`10,000`$ square degrees of the northern Galactic cap, matching the survey area of the Sloan Digital Sky Survey. Given the large coverage of FIRST, its sensitivity is unprecedented: with a limiting flux density of $`1`$ mJy ($`5\sigma _{rms}`$) and an angular resolution of $`5\stackrel{}{.}4`$ (FWHM) the catalogue contains about 100 sources per square degree with a completeness level of $`95\%`$ at 2 mJy (Becker et al. 1995).

We have obtained the publicly available 2001 October 15 version of the source catalogue222http://sundog.stsci.edu/, which has been derived from the 1993 through 2001 observations, and covers about 8,565 square degrees of the sky. About 4% of the 771,076 sources in the catalogue are flagged as possible side-lobes, which we exclude from the catalogue. We set the lower flux density limit of the catalogue to 3 mJy, the limiting flux density of the NVSS survey. Finally, we select the regions $`+2\mathrm{°}\delta +20\mathrm{°}`$ and $`9^h\alpha 16^h`$, $`+20\mathrm{°}\delta +55\mathrm{°}`$ and $`8^h\alpha 17^h`$ from the catalogue, by requiring a relatively uniform source density and a simple geometric form. This area covers $`5,538`$ square degrees and contains 188,885 sources. As for the NVSS, we construct a map of the FIRST surface density (Fig. 3). and plot the source density as a function of declination (Fig. 5). For the selected region we found no suspicious features in the catalogue. The number of sources in various FIRST subsamples are listed in Table 2.

## 4 The angular clustering of radio sources

### 4.1 The angular correlation function of $`S>10`$ mJy NVSS sources

Following the procedures described in Sect. 2 we compute $`w(\theta )`$ for the $`S>10`$ mJy NVSS subsample. Distances between data and/or random positions are initially measured in bins of $`0\stackrel{}{.}5`$, and rebinned in bins of constant logarithmic spacing to analyse the data. We fit the data using a weighted $`\chi ^2`$-minimization routine, and we determine the $`1\sigma `$ errors from the covariance matrix.

The results are shown in Fig. 6. We find that two power-laws are needed to describe the full range of our measurements. Fitting the data with a power-law angular correlation function $`w(\theta )=A\theta ^{1\gamma }`$ (e.g. Peebles 1980) at angular scales of $`\theta 6\mathrm{}`$ gives a slope of $`\gamma =4.4\pm 0.2`$, while at $`\theta 6\mathrm{}`$ we find a slope of $`\gamma =1.7\pm 0.1`$. The latter value is consistent with the slope of the empirical power-law of $`\gamma 1.8`$ found for the cosmological clustering of objects ranging from normal galaxies to clusters (see Bahcall & Soneira 1983, for a review). However, at small angular scales the power-law is much steeper, presumably caused by the enhancement of $`DD(\theta )`$ due to the decomposition of large radio galaxies into their separate radio components (see Sect. 4.2 & Sect. 4.4; see also Blake & Wall 2002a). If we fit the data simultaneously with a double power-law correlation function of the form $`w(\theta )=B\theta ^{1\gamma _B}+A\theta ^{1\gamma _A}`$ with fixed slopes of $`\gamma _B=4.4`$ and $`\gamma _A=1.8`$, we find amplitudes of $`B=(1.5\pm 0.2)\times 10^6`$ and $`A=(1.0\pm 0.2)\times 10^3`$. The double power-law fit is indicated in Fig. 6.

### 4.2 The effect of multiple component radio sources and the angular correlation function of FIRST

Although the median angular size of radio sources is $`10\mathrm{}`$ (e.g. Condon et al. 1998), radio sources can have sizes of up to several arcminutes. At angular scales comparable to the size of these large radio galaxies, the true cosmological $`w(\theta )`$ can become confused or even dominated by resolving these galaxies into their various radio components, such as lobes, hot spots and cores. The angular scale at which the size distribution of radio galaxies begins to dominate $`w(\theta )`$ is indicated by the clear break around $`6\mathrm{}`$. Earlier studies attempted to correct $`w(\theta )`$ for the contribution of multi-component radio sources by means of component combining algorithms. For example, Cress et al. (1996) calculated the angular correlation function for the FIRST survey considering all sources within $`1\stackrel{}{.}2`$ of each other as a single source. The analysis of the FIRST data was repeated by Magliocchetti et al. (1998), who removed double sources using an algorithm based on the $`\theta \sqrt{S}`$ relation of Oort et al. (1987) and flux ratio statistics of the components of genuine doubles. They found values of $`\gamma =2.5\pm 0.1`$, and $`A=(1.0\pm 0.1)\times 10^3`$ for flux density limits between 3 and 10 mJy. Comparing their results to our measurement for the NVSS presented in Fig. 6, we conclude that despite the efforts of these authors it is likely that a residual contribution from large radio galaxies remained. Fitting the data over the whole range of $`\theta `$ with a single power-law explains the apparently high value of $`\gamma 2.5`$ reported for the clustering of FIRST radio sources.

Here, we present new measurements from the FIRST survey. Our reasons for repeating the work of Cress et al. (1996) and Magliocchetti et al. (1998) are threefold. Firstly, the FIRST catalogue has almost doubled in size, enabling a better statistical measure of $`w(\theta )`$. Secondly, the clear break found in the angular correlation function of the NVSS enabled us to isolate the signal due to true clustering from the signal due to the size distribution of radio galaxies. A similar analysis can be applied to the FIRST data. Thirdly, we found large-scale gradients in the NVSS source density below a flux density limit of 10 mJy (see Sect. 3.2). The FIRST data can be used to verify and complement the results from the NVSS for 3–10 mJy.

In Fig. 7 we present our measurements for the angular correlation function from the $`S>3`$ mJy FIRST subsample. As for the NVSS, we see a clear break in $`w(\theta )`$ due to the presence of multi-component radio sources. Fitting the measurements with our double power-law model yields $`\gamma _B=4.1\pm 0.2`$ and $`B=(2.7\pm 0.3)\times 10^6`$, and $`\gamma _A=1.9\pm 0.2`$ and $`A=(1.0\pm 0.3)\times 10^3`$. Note that the break in $`w(\theta )`$ in this sample occurs at $`\theta 4\mathrm{}`$ compared to $`\theta 6\mathrm{}`$ for $`S>10`$ mJy in NVSS (see Fig. 6). Blake & Wall (2002a) show that this is due to a $`1/\sigma `$ dependency ($`\sigma `$ being the surface density of radio sources) of the amplitude of $`w(\theta )`$ at small angular scales, simply because the weight of pair-counts due to large radio galaxies increases as the surface density decreases (see their equation 4).

We conclude that the cosmological $`w(\theta )`$ of $`S>10`$ mJy NVSS sources and $`S>3`$ mJy FIRST sources, as determined by our analysis, are consistent with having the canonical clustering power-law slope of $`\gamma 1.8`$, and an amplitude of $`A1\times 10^3`$.

### 4.3 $`w(\theta )`$ as a function of flux density limit

To investigate angular clustering as a function of flux density limit, we calculate $`w(\theta )`$ for all NVSS and FIRST subsamples listed in Table 2. We obtain the amplitudes of $`w(\theta )`$ by fitting the data with the double power-law model $`w(\theta )=B\theta ^{1\gamma _B}+A\theta ^{1\gamma _A}`$, fixing the slopes at $`\gamma _B=4.4`$ and $`\gamma _A=1.8`$. However, because the signal for the $`S>200`$ mJy subsample is affected by a ’bump’ at $`\theta 0\stackrel{}{.}2`$ (see Fig. 8), we obtained the amplitudes for this subsample by fitting the small- and large-scale correlation functions separately with power-laws $`w(\theta )=B\theta ^{3.4}`$ for $`\theta 0\stackrel{}{.}1`$ and $`w(\theta )=A\theta ^{0.8}`$ for $`\theta 0\stackrel{}{.}3`$, respectively. The measured amplitudes and their $`1\sigma `$ errors are listed in Table 3. The values of both $`B`$ and $`A`$ are found to increase with increasing flux density limit of the subsamples. The increase in $`B`$ can be explained by the $`1/\sigma `$-dependency of the small-scale correlation function that is dominated by double or multiple component radio sources (see Sect. 4.2). From this point onward, we will be only concerned with the amplitude $`A`$ that is believed to be dominated by the true cosmological clustering. In Fig. 9 we have plotted the amplitude of the cosmological $`w(\theta )`$ as a function of flux density limit. For comparison, we have indicated the results from the 325 MHz WENSS and 4850 MHz GB6 surveys (Rengelink 1998; Rengelink & Röttgering 1999) by extrapolating to 1.4 GHz using a power law spectrum, $`S_\nu \nu ^\alpha `$, with spectral index $`\alpha =0.8`$. Between 3 and 40 mJy the amplitude is approximately constant within the errors and has an (unweighted) average of $`1.2\times 10^3`$. From $`50100`$ mJy the amplitude is $`2\times `$ higher, and it has increased by another factor of $`23`$ at 200 mJy. These measurements indicate a trend of increasing clustering amplitude with increasing flux density limit. However, one has to keep in mind that the sources in the brighter subsamples are also included in the subsamples with lower limiting flux densities. Therefore, we also compute $`w(\theta )`$ for sources that lie in the flux interval $`10<S<40`$ mJy. The results are shown in Fig. 8 together with $`w(\theta )`$ found for $`S>200`$ mJy. The amplitude $`A=(6.6\pm 1.8)\times 10^3`$ that we measure for $`S>200`$ mJy is significantly higher than the amplitudes measured at lower flux densities. This is consistent with Rengelink (1998) and Rengelink & Röttgering (1999) who found $`A=(11.5\pm 3.5)\times 10^3`$ for $`S_{1.4}160`$ mJy in the GB6 survey and Loan et al. (1997) who estimated that $`A`$ has a value between 0.005 and 0.015 for $`S_{1.4}>100270`$ mJy from the combined 87GB and PMN surveys (Fig. 9).

We would like to make the following remarks:
(i) Rengelink (1998) and Rengelink & Röttgering (1999) measured $`w(\theta )`$ from WENSS and GB6 by excluding the first $`5\mathrm{}`$ and $`10\mathrm{}`$, respectively. We have used our routines to measure $`w(\theta )`$ for their catalogues as well (not shown here). The amplitudes and slopes we find are consistent with their values, and we find no evidence for a contribution of multi-component sources at the smallest angular scales allowed by these surveys.
(ii) Below 10 mJy the amplitudes for the NVSS and FIRST data are consistent with $`A1.1\times 10^3`$. However, at 10 mJy the amplitude is $`2\times `$ higher for FIRST than for the NVSS. This is curious since the NVSS and FIRST surveys probe radio sources at exactly the same frequency. Blake & Wall (2002b) give a very nice demonstration (see their Fig. 3) of the most probable cause. The resolution of FIRST is ten times higher than that of NVSS, and therefore the average flux density of a single NVSS source is only equal to the sum of all its possibly resolved components in FIRST. Sources that appear in NVSS with integrated fluxes just above a given flux density limit can thus be missed in FIRST. Therefore, we consider NVSS to be more optimal than FIRST for measuring the clustering of extra-galactic radio sources. Furthermore, if we compute $`w(\theta )`$ for only those NVSS sources that lie in the region covered by FIRST, we find an amplitude of $`A=(1.7\pm 0.3)\times 10^3`$. This is consistent with the results found for the 10 mJy FIRST sample, suggesting that cosmic variance of clustering may be an additional factor contributing to the difference in amplitudes measured for the total NVSS area and FIRST. Future work might show that the region covered by FIRST is especially rich in large-scale structures.
(iii) In the 200 mJy subsample we find an unexpected increase in the correlation signal at $`\theta 0\stackrel{}{.}2`$ (indicated by the connected points in Fig. 8). We investigate two possibilities. (1) Sidelobes: Cress et al. (1996) found a bump in $`w(0\stackrel{}{.}1)`$ for $`S>3`$ mJy sources in FIRST, and found that it was caused by sidelobe contamination. However, if sidelobes are responsible for boosting the correlation function at $`\theta 0\stackrel{}{.}2`$ in the $`S>200`$ mJy NVSS sample, these sidelobes themselves also must have minimum peak fluxes of 200 mJy. It is highly unlikely that such bright sidelobes have found their way into the NVSS catalogue, without being masked in Sect. 3.2. Also, we have visually inspected the contour maps of several tens of source pairs ($`S>500`$ mJy) that contribute to $`w(\theta )`$ at $`\theta 0\stackrel{}{.}2`$. In all cases the pairs consisted of unresolved peaks without signs of diffuse, extended emission or side-lobe contamination. (2) Radio galaxies with large angular sizes: the position of the bump near the break in $`w(\theta )`$ suggests that it may somehow be related to the size distribution. Conveniently, Lara et al. (2001) have constructed a sample of 84 large angular size ($`\theta 4\mathrm{}`$) radio galaxies from the NVSS at $`\delta +60\mathrm{°}`$ and a total integrated flux density of $`100`$ mJy. Candidates were pre-selected by visual inspection of the NVSS maps, and confirmed or rejected following observations at higher resolution. If the bump is caused by $`12\mathrm{}`$-sized radio galaxies, then given the 2-Mpc linear size cutoff of large radio galaxies (see Schoenmakers et al. 2001), these galaxies must lie at $`z0.1`$. It is unlikely that such a large, relatively nearby source with, among other emission, two radio components each with a peak flux of $`200`$ mJy would have been missed by their selection criteria. Lara et al. (2001) determined angular sizes by either measuring the maximum distance between $`3\sigma `$ contours, or by the distance between peaks at the source extremes. Also, sizes were measured along the ’spine’ of a source if significant curvature was present. To investigate how many of these sources could actually contribute to $`w(\theta )`$ at $`12\mathrm{}`$ we redetermine the angular sizes of the sources of Lara et al. (2001). We find that none of these sources consists of $`2`$ components of $`200`$ mJy of $`12\mathrm{}`$ separation. On the other hand, if we extrapolate the clustering power-law derived at larger scales to $`\theta =0\stackrel{}{.}2`$ we find that the bump translates into $`10\times `$ the number of pairs expected. Even allowing for the much larger area of NVSS, the possibility that the bump is caused by large radio galaxies as in the sample of Lara et al. (2001) is therefore unlikely.

Unfortunately, the exact origin of this feature remains unclear. We realize, however, that this bump is situated at a crucial angular scale for our measurements. Therefore, we have obtained the amplitudes $`B`$ and $`A`$ by fitting $`w(\theta )`$ on both sides of the bump with a single power-law. Under the condition that the effect that causes the bump is not responsible for enhancing $`w(\theta )`$ at $`\theta 0\stackrel{}{.}3`$, this will enable us to derive an estimate for the amplitude for the cosmological clustering. At $`\theta 0\stackrel{}{.}3`$ $`w(\theta )`$ is consistent with the classical $`\gamma =1.8`$ power-law clustering model.

### 4.4 Modelling the angular size distribution of radio galaxies

#### 4.4.1 The model

The steepening of the slope of $`w(\theta )`$ at small angular scales is presumably related to multi-component sources spuriously enhancing the true clustering pair counts at small $`\theta `$. To demonstrate the reality of this assumption, we create a simple model for the angular size distribution of radio galaxies in the NVSS, that is able to account for this extra signal contributing to $`w(\theta )`$. We model the physical size distribution of sources in our $`S>10`$ mJy NVSS sample, and use their redshift distribution to obtain the angular size distribution. Because we know the angular resolution of the NVSS, this model can then be used to estimate the fraction of sources likely to be resolved. It is essential to separate sources that are resolved into a single, elongated object from sources that are resolved into a number of components, since only the latter would produce extra pair counts. Here, we assume that the majority of surplus pair counts arise from resolving the two edge-brightened radio lobes of FRII-type radio galaxies (see Fanaroff & Riley 1974), and we estimate that the fraction of FRIIs at 10 mJy is $`40`$% from Wall & Jackson (1997) (assuming a spectral index of $`\alpha =0.8`$ to extrapolate to 1.4 GHz).

Several groups have investigated the median physical sizes of FRII radio galaxies as a function of redshift and radio luminosity by parameterizing the linear size as $`D(1+z)^nP^m`$, where $`P`$ is the radio luminosity (for a review see Blundell et al. 1999). Results using different samples of radio galaxies vary from no size evolution at all (e.g. Nilsson et al. 1993), to size evolution depending only on redshift (e.g. Kapahi et al. 1987), and size evolution depending on both redshift and luminosity with contradictory results (e.g. Oort et al. 1987; Barthel & Miley 1988; Singal 1993). We use the results of Neeser et al. (1995) who found the following linear size–redshift relation from a spectroscopically complete sample of FRII radio galaxies:

$$D(1+z)^{1.7\pm 0.5}\text{(for }\mathrm{\Omega }_M=1\text{ and }\mathrm{\Omega }_\mathrm{\Lambda }=0\text{)},$$
(3)
and remark that no intrinsic correlation was found between $`D`$ and $`P`$ ($`P^m1`$ with $`m=0.06\pm 0.09`$). This observed linear-size evolution may be related to evolution of the confining intergalactic medium, or to evolution of the radio galaxy itself, but the exact underlying physical mechanism is unknown (see Neeser et al. 1995).

For the purpose of our model, we place simulated sources in small redshift intervals ($`\mathrm{\Delta }z=0.01`$) in the range $`0z5`$, and assume that their mean physical size evolves with redshift according to equation 3. We set the total number of input sources equal to the estimated number of $`S>10`$ mJy FRIIs in our NVSS sample ($`40\%of434,000`$), and calculate the number of sources in each redshift interval from the redshift distribution, $`N(z)`$, using the formalism of Dunlop & Peacock (1990) (see Sect. 5 for details). We then assume that in each redshift interval sizes are normally distributed. We take a mean size of 500 kpc and a standard deviation of 250 kpc at $`z0`$, chosen so that the resulting physical size distribution roughly resembles the distribution of projected linear sizes versus redshift as it is given by Blundell et al. (1999) for three complete samples of FRII radio galaxies from the 3C, 6C, and 7C radio surveys. The resulting physical size distribution is shown in Fig. 10, where we plot filled contours of the source density in the linear size-redshift plane to illustrate the underlying redshift distribution. We have also indicated the minimum physical size that is theoretically required for a source to become resolved as a function of redshift, given by the NVSS resolution of 45″ (FWHM). We would like to remark at this point that the distribution of sizes in our model beyond redshifts of $`z3`$ should not be taken too seriously as it is based on a straight extrapolation from measurements made at redshifts $`0z2`$, and does not take into account the fact that at these high redshifts most sources will be extremely young and are thus likely to be very small. However, as can be seen from Fig. 10, our modeled size distribution falls below the NVSS resolution already at $`z1`$. Taking smaller sizes at higher redshifts will have no effect on the modeled size distribution of resolved sources that we want to derive here.

Assuming $`\mathrm{\Omega }_M=1`$ we calculate the angular size distribution associated with our model. We construct 10 such models, and average them to get our final model of the angular size distribution of the sample. This model is presented in Fig. 11. Although the mean angular size is $``$10$`\mathrm{}`$ in agreement with Condon et al. (1998), sizes are found to extend up to several arcminutes beyond the resolution of the NVSS (indicated by the dotted line).

#### 4.4.2 Results

We now compare the number of surplus pairs expected from resolved FRII sources in the model, $`DD_{mod}(\theta )`$, to the actually measured pair counts at angular scales of $`\theta 6\mathrm{}`$. At these scales, the measured pair counts consist of both pair counts due to clustering and pair counts due to doubles, so

$$DD_{tot}(\theta )=DD_{gal}(\theta )+DD_{dbl}(\theta ).$$
(4)
To extract $`DD_{dbl}(\theta )`$ from the total counts, $`DD_{tot}(\theta )`$, we calculate $`DD_{gal}(\theta )`$ by assuming that the galaxy angular correlation function as measured above the break in $`w(\theta )`$ can be extrapolated to angular scales of $`\theta 6\mathrm{}`$:

$$w_{gal}(\theta )=1.0\times 10^3\theta ^{0.8}=DD_{gal}(\theta )F(\theta )1,$$
(5)
where $`F(\theta )=4RR(\theta )/[DR(\theta )]^2`$, the part of the Hamilton estimator that is relatively independent of the presence of doubles. Since we now know both $`DD_{tot}(\theta )`$ and $`DD_{gal}(\theta )`$, we can subtract them to get a measure of the counts arising from the double sources: $`DD_{dbl}(\theta )`$. The final step is to rebin the modeled number of pair separations $`DD_{mod}(\theta )`$ in order to match the binning scheme of $`DD_{dbl}(\theta )`$. Fig. 12 shows the ratio of the observed doubles to the modeled doubles per distance interval. The errors in the observed counts are estimated from the $`1\sigma `$-error in the amplitude of $`w(\theta )`$. The errors in the modeled pair counts are estimated by allowing a 10% error in the estimated fraction of FRIIs in the NVSS. We conclude that: a model in which the small-scale angular correlation function steepens due to resolving FRII radio galaxies into two distinct knots of radio emission is in good agreement with the measurements presented in Fig. 6.

Several remarks that can be made are the following:
(i) The size distribution of radio sources at the mJy level is still largely unconstrained. Recently, however, Lara et al. (2001) presented a new sample of large radio galaxies (LRGs) selected from the NVSS. In the region $`\delta +60\mathrm{°}`$ they found $``$80 radio galaxies with apparent angular sizes larger than $`4\mathrm{}`$ and total flux density greater than 100 mJy. If we roughly extrapolate our model to their sensitivity and correct for the area we successfully predict the number of FRIIs in the range $`4\mathrm{}\theta 6\mathrm{}`$. However, in this interval one third of the sample of Lara et al. (2001) consists of FRIs, while the model only uses FRIIs to estimate the number of surplus pairs expected. The model could be refined by decreasing the fraction of resolved FRIIs to also allow a contribution from large FRIs.
(ii) The model allows objects to be either single or double sources, although visual inspection of NVSS contour maps shows that sources are sometimes split into three or even more components. Therefore, we may expect an extra amount of spurious pair counts on top of the counts due to classical double radio sources. This may become increasingly important with increasing flux density limit.
(iii) The model predicts a fraction of resolved sources in NVSS of $``$0.07, in rough agreement with the value of $``$0.05 predicted by Condon et al. (1998).

The simple model allows us to explore the general relations between the physical size distribution of radio galaxies and $`w(\theta )`$ at small angular scales. Although our crude method is successful in reproducing the observations, it relies on a number of assumptions that are not easily verified from the data currently in literature. Radio sources come in a wide variety of sizes ranging from $`<1`$ kpc for the class of gigahertz peaked spectrum sources (GPS), to 1–20 kpc for the compact steep spectrum sources (CSS), $`>20`$ kpc for FRI- and FRII-type radio galaxies, and $`>1`$ Mpc for giant radio galaxies (Fanti et al. 1990; O’Dea et al. 1991; Blundell et al. 1999; Schoenmakers et al. 2001). Evidently, the distribution of linear sizes of radio sources are very complex, and will remain an important subject for future studies. As we have shown, the angular correlation function can be used to put constraints on the size distribution of large radio galaxies. However, perhaps more ideal would be to make a statistical redshift sample of all radio source pairs within some angular distance interval, and then take high resolution radio observations to constrain the numbers of intrinsic doubles in that sample.

## 5 The spatial clustering of NVSS souces

### 5.1 The redshift distribution

At the mJy level and higher it is standard practice to compute redshift distributions using the Dunlop & Peacock (1990) radio luminosity functions (RLFs). These authors have constructed a range of model luminosity functions using spectroscopically complete samples from several radio surveys at different frequencies. Using a free-form modelling approach they found a number of smooth functions that were consistent with the data. In addition, they attempted two models of a more physical nature by assuming pure luminosity evolution (PLE) and luminosity/density evolution (LDE) to describe the RLF. The total ensemble is expected to agree well at those luminosities and frequencies at which they are best constrained by the data, while uncertainties in the extrapolation of each of these models to those regions that are less constrained by the data may be reduced by taking the ensemble as a whole. We compute redshift distributions, $`N(z)`$, for each flux-limited subsample using the free-form models $`14`$ and the PLE/LDE models for the combined population of flat ($`\alpha =0`$, $`S_\nu \nu ^\alpha `$) and steep ($`\alpha =0.8`$) spectrum radio sources given by Dunlop & Peacock (1990, taking the MEAN-$`z`$ data from their appendix C) from

$`{\displaystyle \frac{dN(z)}{dz}}`$ $`=`$ $`{\displaystyle \frac{dV(z)}{dz}}\times {\displaystyle _{P_{low}(z)}^{\mathrm{}}}\mathrm{\Phi }_i(P,z)𝑑P,`$ (6)
$`P_{low}(z)`$ $`=`$ $`x(z)^2\left({\displaystyle \frac{S}{(1+z)^{1\alpha }}}\right)\left({\displaystyle \frac{2.7GHz}{\nu }}\right)^\alpha ,`$
where $`V(z)`$ is the comoving volume, $`\mathrm{\Phi }_i(P,z)`$ is the model RLF, $`x(z)`$ the comoving distance, $`S`$ the limiting flux density of the subsample, and $`\nu `$ the frequency of FIRST/NVSS. We note that $`N(z)`$ is independent of cosmology as long as the calculations are carried out in the cosmology used to construct the RLFs (i.e. $`\mathrm{\Omega }_M=1.0`$ and $`H_0=50`$ km s-1 Mpc-1).

Figs. 13 and 14 show the redshift distributions for $`S>10`$ mJy and $`S>100`$ mJy, respectively. We calculate the average of the six different models (indicated by the solid curve), which will be our best estimate of $`N(z)`$ use in the analysis below (the same method was used for the $`N(z)`$ applied to the model of the angular size distribution described in Sect. 4.4). It is important to keep in mind that the functional form of $`N(z)`$ remains virtually unchanged from $`3200`$ mJy. Over this range in flux densities the RLFs represent a broad redshift distribution with a peak around $`z1`$, indicating the very large median redshift that is generally probed by radio surveys.

### 5.2 The spatial correlation function

Given the amplitudes of $`w(\theta )`$ determined in Sect. 4 we can use the cosmological Limber equation to estimate the spatial correlation length, $`r_0`$, by deprojecting $`w(\theta )`$ into the spatial correlation function, $`\xi (r)`$ using the redshift distribution and cosmology (e.g. Peebles 1980, §56). We consider two cosmological models: a flat, vacuum dominated, low-density Universe ($`\mathrm{\Lambda }`$CDM; $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$), and an Einstein-de Sitter model Universe ($`\tau `$CDM; $`\mathrm{\Omega }_M=1.0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$). We use $`H_0=100h`$ km s-1 Mpc-1.

We assume an epoch dependent power-law spatial correlation function of the form

$$\xi (r_p,z)=\left(\frac{r_p}{r_0}\right)^\gamma (1+z)^{(3+ϵ)},$$
(7)
where $`r_p`$ is the proper distance, $`r_0`$ is the spatial correlation length333Note that the spatial correlation length is not a physical lengthscale in the space distribution of galaxies. It is just defined as that length at which $`\xi (r)`$ is unity (i.e. the chance of finding a galaxy at the distance $`r_0`$ from another galaxy is twice the Poissonian chance). at $`z=0`$, and $`ϵ`$ parameterizes the redshift evolution of the clustering. To express $`\xi (r_p,z)`$ in terms of comoving coordinates $`r_c=r_p\times (1+z)`$, we write:

$$\xi (r_c,z)=\left(\frac{r_c}{r_0}\right)^\gamma (1+z)^{\gamma (3+ϵ)},$$
(8)
which can be written as

$$\xi (r_c,z)=\left(\frac{r_c}{r_0(z)}\right)^\gamma ,r_0(z)=r_0(1+z)^{1\frac{3+ϵ}{\gamma }},$$
(9)
where $`r_0(z)`$ is the (comoving) correlation length measured at $`z`$. In a flat model Universe, the cosmological Limber equation can be expressed as follows (see e.g. Peebles 1980):

$`w(\theta )`$ $`=`$ $`A\theta ^{1\gamma }=\sqrt{\mathrm{\Omega }_M}\left({\displaystyle \frac{r_0H_0}{c}}\right)^\gamma \theta ^{1\gamma }H_\gamma \times `$
$`{\displaystyle \frac{_0^{\mathrm{}}𝑑zN(z)^2(1+z)^{\gamma 3ϵ}x^{1\gamma }Q(z)}{\left[_0^{\mathrm{}}𝑑zN(z)\right]^2}},`$
with

$`Q(z)`$ $`=`$ $`\left[(1+z)^3+\mathrm{\Omega }_M^11\right]^{0.5},`$ (11)
$`x(z)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\mathrm{\Omega }_M}}}{\displaystyle _0^z}{\displaystyle \frac{dz}{Q(z)}},`$
$`H_\gamma `$ $`=`$ $`\mathrm{\Gamma }\left({\displaystyle \frac{1}{2}}\right)\mathrm{\Gamma }\left({\displaystyle \frac{\gamma 1}{2}}\right)\mathrm{\Gamma }\left({\displaystyle \frac{\gamma }{2}}\right)^1,`$
and using the approximation that angles are small ($`\theta 1`$). We calculate $`N(z)`$ for each subsample.

The evolution parameter $`ϵ`$ can represent a variety of clustering models. Three important cases are the following (see Phillipps et al. 1978; Kundić 1997). (1) The stable clustering model ($`ϵ=0`$): if galaxy clustering is gravitationally bound at small scales, then clusters have fixed physical sizes (i.e. they will neither contract nor expand) and will have a correlation function that decreases with redshift as $`(1+z)^{1.2}`$. (2) The comoving clustering model ($`ϵ=\gamma 3`$): galaxies and clusters expand with the Universe, so their correlation function remains unchanged in comoving coordinates. This case applies well to a low density Universe where there is not enough gravitational pull to counterbalance expansion, and implies that structures have formed very early. (3) The linear growth model ($`ϵ=\gamma 1`$): clustering grows as expected under linear perturbation theory.

Studies of the spatial clustering properties of radio-quiet quasars indicate that the clustering history of active galaxies, unlike that of normal galaxies, is best characterized using a negative value for $`ϵ`$. Kundić (1997) measured the high-redshift quasar-quasar correlation function from the Palomar Transit Grism Survey, and found no evidence for a decrease in the correlation amplitude of quasars with redshift. Moreover, he found that $`\xi _{qq}(z>2)/\xi _{qq}(z<2)1.8`$, suggesting an even higher amplitude at higher redshifts. Similarly, Croom et al. (2001) find almost no evolution in clustering strength for quasars taken from the 2dF QSO Redshift Survey out to $`z2.5`$.

Therefore, we opt for evolution model 2 (i.e. constant clustering in comoving coordinates), which implies $`ϵ=1.2`$ for $`\gamma =1.8`$. In table 4 we list the results obtained using this model for the two different cosmological models. For comparison, we also indicate the results using the stable clustering model ($`ϵ=0`$). For $`ϵ=0`$ the present-day correlation length is $`1.4`$ times higher than for $`ϵ=1.2`$ in both cosmologies. However, given the strong peak in the redshift distribution at $`z1`$, we are effectively measuring clustering at $`z1`$. Calculating $`r_0(z1)`$ in the case of stable clustering using Eq. 9 yields a value that is only $`1.1`$ times lower than $`r_0(z1)=r_0`$ in the case of $`ϵ=1.2`$. Therefore, the value of $`r_0(z1)`$ is relatively independent of the exact value of $`ϵ`$. The results for the $`ϵ=1.2`$ ($`\mathrm{\Lambda }`$CDM) case are presented in Fig. 15. We find an approximately constant spatial correlation length of $`6.0`$ $`h^1`$ Mpc from 3–40 mJy, compared to $`14`$ $`h^1`$ Mpc at 200 mJy.

As we have shown, the possibility that the observed flux-dependency of the clustering is just an effect of projection can be ruled out, since the shape of the redshift distribution is relatively constant with flux over several orders of magnitude (at least above $`1`$ mJy). This automatically implies that the average radio power of the subsamples increases with flux density (indicated by the top axis of Fig. 15). An alternative explanation was therefore suggested by Rengelink (1998) and Rengelink & Röttgering (1999) based on their measurements of the clustering of radio sources in the WENSS and GB6 surveys. They concluded that the clustering signal could change as a function of flux density if relatively low and high power radio galaxies represent different spatial structures at a similar epoch ($`z1`$). Taking the predicted population mix of radio sources from Wall & Jackson (1997), we find that for $`S_{1.4}>10`$ mJy the fractions of FRIs and FRIIs are about equal. However, for $`S_{1.4}>100`$ mJy the fraction of FRIIs is more than $`75`$%. Given the fractional changes of the source populations with flux density limit, the clustering amplitudes measured are very well matched by a scenario in which the clustering of powerful radio sources (mostly FRII) and average power radio sources (FRI/FRII) are intrinsically different, with FRIIs being more strongly clustered at $`z1`$ than the radio galaxy population on average.

As pointed out by Rengelink (1998) and Rengelink & Röttgering (1999) the large difference in observing frequencies and sensitivities of WENSS and GB6 (the limiting 1.4 GHz flux densities probed by these surveys correspond to 10 mJy for WENSS and 70 mJy for Greenbank, respectively) only allowed them to make a comparison between the results, whereas the detection of the inferred flux-dependency of $`r_0`$ within a single survey would be highly desirable. Our analysis of the clustering in the single large-area, intermediate-frequency NVSS survey is in agreement with their conclusions.

## 6 Discussion

### 6.1 Clustering measurements from literature

We start this section by making a survey of other clustering measurements from literature. However, readers may wish to skip directly to Sect. 6.2 for a discussion on these measurements and the results presented in this paper in their cosmological context.

In order to compare results from different studies, all values taken from literature were converted assuming a fixed slope $`\gamma =1.8`$ by setting $`r_{0,1.8}=(r_{0,\gamma })^{\gamma /1.8}`$. All correlation lengths are expressed in comoving units, and we have transformed all values to a $`\mathrm{\Lambda }`$CDM cosmology (see Magliocchetti et al. 2000). Please note that the list given below is not complete, and the reader is kindly invited to consult the individual papers and the references therein for further information.

#### 6.1.1 Clusters

Estimates of the correlation length of rich Abell clusters are given by Bahcall & Soneira (1983) and Postman et al. (1992) who found $`r_0=24\pm 9`$ $`h^1`$ Mpc. Lahav et al. (1989) found $`r_0=21\pm 7`$ $`h^1`$ Mpc from an all-sky sample of the brightest X-ray clusters, and Dalton et al. (1994) and Croft et al. (1997) found $`r_0=19\pm 5`$ $`h^1`$ Mpc and $`r_0=16\pm 4`$ $`h^1`$ Mpc, respectively, for clusters selected from the APM Galaxy Survey. Recently, Gonzalez et al. (2002) measured the correlation length of distant clusters in the Las Campanas Distant Cluster Survey and found a correlation length of $`24.8\pm 4.5`$ $`h^1`$ Mpc at $`\overline{z}=0.42`$.

Different studies may have sampled clusters of different degrees of richness, which can account for most of the scatter in the reported values. In general, however, all results are consistent with clusters being the most strongly clustered objects known in the Universe.

#### 6.1.2 Optically-selected ordinary galaxies and IRAS galaxies

Bright early-type galaxies are found to have a strongly clustered distribution in the local Universe. Willmer et al. (1998) find $`r_0=6.8\pm 0.4`$ $`h^1`$ Mpc for local $`LL_{}`$ ellipticals, and Guzzo et al. (1997) measure a considerably higher $`r_0=11.4\pm 1.3`$ $`h^1`$ Mpc for a sample of similar galaxies. Although these results are only consistent with each other at the $`3\sigma `$ level, the latter sample contains a higher fraction of local clusters, presumably responsible for boosting the $`r_0`$. The dependence of galaxy clustering on luminosity and spectral type has been studied using the ongoing 2 degree Field Galaxy Redshift Survey (2dFGRS). Norberg et al. (2002) find $`r_0=11.8\pm 1.6`$ $`h^1`$ Mpc for the brightest early-type galaxies in the 2dFGRS. Moreover, they find a strong dependence of clustering strength on luminosity, with the amplitude increasing by a factor of $`2.5`$ between $`L_{}`$ and $`4L_{}`$. The ordinary population of galaxies has been found to be less strongly clustered than the population consisting of local (bright) ellipticals: Loveday et al. (1995) find $`r_0=4.7\pm 0.2`$ $`h^1`$ Mpc from the APM survey. At higher redshifts, the clustering strength in a sample of faint $`K`$-selected galaxies with minimum rest-frame luminosities of $`M_K=23.5`$, or about $`0.5L_{}`$, is found to be fairly rapidly declining with redshift: Carlberg et al. (1997) find $`r_0=3.3\pm 0.1`$ $`h^1`$ Mpc, $`r_0=2.3\pm 0.2`$ $`h^1`$ Mpc, $`r_0=1.6\pm 0.2`$ $`h^1`$ Mpc, and $`r_0=1.2\pm 0.2`$ $`h^1`$ Mpc, at $`\overline{z}=0.34`$, $`\overline{z}=0.62`$, $`\overline{z}=0.97`$, and $`\overline{z}=1.39`$, respectively. Carlberg et al. (2000) present measurements on a sample of $`LL_{}`$ galaxies up to $`z0.6`$ and find a much milder decline from $`r_0=5.1\pm 0.1`$ $`h^1`$ Mpc at $`\overline{z}=0.10`$ to $`r_0=4.2\pm 0.4`$ $`h^1`$ Mpc at $`\overline{z}=0.59`$.

Clustering of the local population of IRAS-selected galaxies is best fit by $`r_0=3.4\pm 0.2`$ $`h^1`$ Mpc (Fisher et al. 1994).

#### 6.1.3 Extremely red objects (EROs)

Several recent studies indicate that the comoving correlation length of early-type galaxies undergoes little or no evolution from $`0z1`$. Evidence for this is provided by the clustering of extremely red objects, a population of galaxies having very red optical to infrared colors ($`RK_s>5`$). These red colors are consistent with them being either old, passively evolving elliptical galaxies, or strongly dust-enshrouded starburst galaxies at $`z11.5`$. Indeed, further observations have confirmed that both classes are present in the ERO population (e.g. Dunlop et al. 1996; Cimatti et al. 1998; Dey et al. 1999; Liu et al. 2000). Daddi et al. (2001) have recently embarked on a study of the spatial clustering of a large sample of $`LL_{}`$ EROs at $`z1`$, and found a large correlation length of $`r_0=12\pm 3`$ $`h^1`$ Mpc. In Cimatti et al. (2002) the results are presented involving the EROs that were identified in a large flux limited redshift survey of $`500`$ galaxies with $`K20`$. The derived fraction of early-type EROs from that sample is $`50\pm 20\%`$, while there is an increasing contribution of dusty star-forming EROs at faint magnitudes. Therefore, Daddi et al. (2002) have attempted to analyse separately the spatial clustering of EROs from both categories by studying the frequency of close pairs. They find that the comoving correlation length of the dust-enshrouded starbursts is constrained to be less than $`r_0=2.5`$ $`h^1`$ Mpc, while the old EROs are clustered with $`5.5r_016`$ $`h^1`$ Mpc. This is consistent with the value reported earlier in Daddi et al. (2001), which is still valid as a lower limit for the clustering of early-type EROs based on the argument that the much less clustered dusty star-forming EROs only dilute the clustering signal coming from the ellipticals in this sample (see also Roche et al. 2002). Furthermore, McCarthy et al. (2001) have identified a large sample of such faint red galaxies as being consistent with mildly evolved early-type galaxies at $`z1.2`$. They find a clustering strength of $`r_0=9.5\pm 1`$ $`h^1`$ Mpc.

#### 6.1.4 Radio galaxies

The results on the spatial clustering of radio sources at $`z1`$ presented in this paper indicate that $`r_0`$ depends on radio luminosity in such a way that very luminous (FRII) radio galaxies cluster more strongly than the total population of radio galaxies (both FRI and FRII) on average, reminiscent of a similar luminosity trend found for samples of optically-selected galaxies. We roughly construct two radio luminosity bins from our measurements by comparing the $`r_0`$ found for 3–40 mJy to the $`r_0`$ found for the 200 mJy subsample. We find $`r_06\pm 1`$ $`h^1`$ Mpc-1 for the relatively low power bin ($`P_{1.4}10^{2425}`$ W Hz-1 sr-1), and $`r_014\pm 3`$ $`h^1`$ Mpc-1 for the high power bin ($`P>10^{26}`$ W Hz-1 sr-1).

#### 6.1.5 Optically-selected quasars

Croom et al. (2001) have determined the correlation length of quasars (QSOs) using 10,558 quasars taken from the 2dF QSO Redshift Survey. They find that QSO clustering appears to vary little with redshift, with $`r_0=4.9\pm 0.8`$ $`h^1`$ Mpc at $`\overline{z}=0.69`$, $`r_0=2.9\pm 0.8`$ $`h^1`$ Mpc at $`\overline{z}=1.16`$, $`r_0=4.2\pm 0.7`$ $`h^1`$ Mpc at $`\overline{z}=1.53`$, $`r_0=5.3\pm 0.9`$ $`h^1`$ Mpc at $`\overline{z}=1.89`$, and $`r_0=5.8\pm 1.2`$ $`h^1`$ Mpc at $`\overline{z}=2.36`$.

#### 6.1.6 Lyman-break galaxies

Lyman-break galaxies (LBGs) are found to be associated with star-forming galaxies at $`z3`$, with comoving correlation lengths of $`r_0=3.3\pm 0.3`$ $`h^1`$ Mpc (Adelberger 2000), and $`r_0=3.6\pm 1.2`$ $`h^1`$ Mpc (Porciani & Giavalisco 2002). Ouchi et al. (2001) find $`r_0=2.7\pm 0.6`$ $`h^1`$ Mpc for a sample of LBGs at $`z4`$.

### 6.2 Clustering evolution

#### 6.2.1 The clustering of massive ellipticals at $`z1`$

In Fig. 16 we present an overview of the evolution of galaxy clustering, as it follows from the broad variety of observational results summarized above. The $`r_0`$ that we measure for the brightest radio sources at $`z1`$ is comparable to the $`r_0`$ measured for bright ellipticals locally, and $`2\times `$ higher than the $`r_0`$ measured for relatively faint radio sources and quasars, suggesting that they are considerably more biased and probably probe spatial structures associated with strongly clustered, massive objects. This does not come totally unexpectedly, as there is a range of observational evidence in support of this result. Best et al. (1998) found that powerful 3CR radio galaxies are mostly associated with massive galaxies at $`z1`$, and at high ($`z1`$) and very high ($`z2`$) redshifts the most luminous (i.e. FRII-type) radio sources are found in very dense environments associated with forming clusters. This is based on for example the presence of large X-ray halos (Crawford & Fabian 1996), excesses of companion galaxies (McCarthy 1988; Röttgering et al. 1996; Nakata et al. 2001), and excesses of Ly$`\alpha `$ emitters around powerful radio sources (Kurk et al. 2000; Venemans et al. 2002). Furthermore, most very high redshift radio galaxies ($`z>2`$) are surrounded by giant halos of emission line gas (e.g. Röttgering et al. 1999; De Breuck et al. 2000), and some have very clumpy morphologies suggestive of massive forming systems (e.g. Pentericci et al. 1999, 2000). Using HST/NICMOS observations, Pentericci et al. (2001) have found a number of radio galaxies at $`z2`$ having morphologies that are represented well by a de Vaucouleurs profile, consistent with them being elliptical galaxies or proto-galaxy bulges.

As argued by Best et al. (1999), powerful radio sources must rely on (i) a plentiful supply of gas to fuel a supermassive blackhole that can drive the AGN activity, and (ii) a dense surrounding medium able to contain the radio lobes. These environments are indeed expected to be found in the gas-rich galaxy clusters at high redshift, additionally supporting the conclusion that high redshift FRIIs are associated with strongly clustered, massive objects. One may argue that this conclusion somewhat contradicts the fact that low redshift FRIIs are primarily found to be situated in small, isolated galaxy groups, and not in the centers of large clusters (Butcher & Oemler 1978; Hill & Lilly 1991). This, however, can easily be explained by considering that the local analogs of the gas-rich cluster environments that are suitable for producing powerful FRIIs at high redshifts, are found in relatively small galaxy groups, and not in the gas-depleted centers of local rich clusters (Rengelink 1998).

Interestingly, we find that both EROs and powerful radio galaxies are strongly clustered with $`r_010`$ $`h^1`$ Mpc at $`z1`$. Willott et al. (2001) suggested that high-redshift radio galaxies and EROs could be identical galaxies seen at different stages of their evolution, based on their findings of ERO-like host galaxies for a number of radio galaxies from the 7C Redshift Survey. This, of course, would be highly consistent with the belief that both radio galaxies and EROs may be the progenitors of local bright ellipticals. They conclude that the density of radio sources with minimum radio luminosities of log10 $`P_{151}=24`$ W Hz-1 sr-1 is consistent with a model in which all EROs go through a relatively short period of AGN activity, forming a radio galaxy somewhere between $`z=2`$ and $`z=1`$.

However, if all EROs are radio galaxies at some stage, their highly clustered spatial distribution should be reflected in the spatial distribution of the radio galaxies. Fig. 16 shows that the clustering of EROs and radio galaxies is consistent only for those galaxies with radio luminosities of log10 $`P_{1400}26`$ W Hz-1 sr-1. The surface density of such radio sources in the redshift range $`1<z<2`$ in the NVSS is $`2\times 10^4`$ arcmin-2, while the surface density of EROs having $`K_s19`$ and $`RK_s>5`$ is $`0.5`$ arcmin-2 (Daddi et al. 2001). If we take the fraction of old ellipticals among EROs to be $`70`$% (Cimatti et al. 2002), then only $`0.06\%`$ of these EROs are currently observed in their radio-loud phase. However, because the typically assumed AGN lifetimes are short compared to the cosmological time-scale from $`z=2`$ and $`z=1`$ ($`t_{z=21}3.5`$ Gyr for $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$), the number of EROs that could undergo a radio-loud phase is $`220\%`$ (assuming $`t_{AGN}10^{78}`$ yr.). These fractions can be increased significantly if, for example, we select EROs that are much redder: the density of EROs having $`RK_s>6`$ is a factor of $`10`$ lower compared to $`RK_s>5`$ (Daddi et al. 2001), giving $`14140\%`$. It may be clear from the above that the unification of EROs and radio galaxies, although tempting, relies on a number of issues that have not yet been resolved. Further study of the luminosities, colors and morphologies of radio galaxy hosts, as well as the cluster environments of EROs may be expected to provide important clues for constraining this scenario.

#### 6.2.2 Comparison with theoretical predictions

Linear ($`ϵ1`$, dot-dashed line) or stable ($`ϵ=0`$, dotted line) clustering evolution models have been found to best fit the measurements of ordinary, optically-selected galaxies at $`z1`$ (e.g. Carlberg et al. 1997, 2000; McCracken et al. 2001, and references therein). However, as Fig. 16 shows, these models do not provide a good description for the evolution of massive early type galaxies as inferred from the meassurements of local bright ellipticals and FRII radio galaxies and EROs at $`z1`$. Adjusting these models to the measurements would either require $`z1`$ massive ellipticals to have a correlation length around $`67`$ $`h^1`$ Mpc, or local bright ellipticals to have a correlation length of the order of that of local clusters, far greater than observed. For these galaxies, the current measurements require a model that predicts relatively constant clustering in comoving coordinates, i.e. a negative value of $`ϵ1`$ in the simple $`ϵ`$-model.

Although the parameterization of clustering evolution by means of the $`ϵ`$-model is useful for characterizing the measurements as a function of redshift, it does not provide good physical insight into evolution governed by the clustering of dark matter halos (see Giavalisco et al. 1998; McCracken et al. 2001). Galaxy clustering evolution can be described more precisely by

$$\xi _{gal}(z,r)=D^2(z)b^2(z)\xi _m(0,r),$$
(12)
where $`D(z)`$ is the linear cosmological growth rate (see Carroll et al. 1992), $`b(z)`$ the evolution of the bias, and $`\xi _m(0,r)`$ the correlation function of the underlying matter distribution at $`z=0`$. Since $`b(z)`$ is related to the nature of the mechanism through which the galaxies were formed, measurements of $`\xi _{gal}(z,r)`$ can be used to constrain structure formation models.

In the galaxy conservation model, objects are formed by means of monolithic collapse at arbitrarily high redshift, and their clustering evolution is described solely by the cosmological growth of density perturbations (Fry 1996). In this model, bias evolves as

$$b(z)=1+(b_01)/D(z),$$
(13)
where $`b_0(\sigma _{8,gal}/\sigma _{8,m})`$ and $`\sigma _8`$ is the rms fluctuation amplitude inside a sphere of 8 $`h^1`$ Mpc radius. Taking $`r_{0,m}(0)=5`$ $`h^1`$ Mpc for the present-day correlation length of the dark matter from the GIF/VIRGO $`N`$-body simulations of Jenkins et al. (1998) (thick solid line in Fig. 16) and $`\overline{r}_{0,gal}(0)=8.7`$ $`h^1`$ Mpc for ellipticals, we find $`\sigma _{8,m}=0.9`$ and $`\sigma _{8,gal}=1.5`$ corresponding to $`b_01.65`$. This model is indicated in Fig. 16 (thin solid line). Analogous to the above arguments against simple stable or linear clustering, extrapolating the clustering of local ellipticals to $`z1`$ in the galaxy conservation model does not fit the observed extreme clustering of EROs and powerful radio galaxies. On the other hand, this scenario shows good agreement with the $`r_06`$ $`h^1`$ Mpc measured for lower luminosity radio sources and QSOs at $`z1`$.

Crucial to the picture that is developing may be the recent results of Wilson (2003), who studied the clustering of $`(VI)`$-selected $`L_{}`$ early-type galaxies in the redshift range $`0.2<z<0.9`$. This author found that these galaxies cluster slightly more strongly compared to the field, with a best-fitting $`ϵ`$-model of $`ϵ=0`$ and $`r_0=5.25\pm 0.28`$ $`h^1`$ Mpc. This is in agreement with the correlation length of local $`L_{}`$ early-types in the 2dFGRS. Wilson (2003) remarks that this measurement is inconstent with the large $`r_0`$ found for EROs, which are also believed to be $`LL_{}`$ early-type galaxies. The value of $`r_0`$ for EROs and radio galaxies could be spuriously high due to uncertainties in their redshift distributions which is not included in the quoted errors, although the selection functions of both EROs and powerful radio galaxies are considered to be understood relatively well (e.g. Dunlop & Peacock 1990; Daddi et al. 2001; McCarthy et al. 2001). Alternatively, EROs and radio galaxies at $`z1`$ may be much more strongly clustered because they correspond to a population of massive, bright cluster galaxies in the process of formation. If FRII radio galaxies and EROs are indeed the distant analogs of local $`LL_{}`$ early-types, they are becoming considerably more biased tracers of the underlying galaxy distribution with redshift, while this galaxy distribution itself probably traces the dark matter distribution with relatively constant bias. Interestingly, (semi-) analytic models and $`N`$-body simulations are able to explain this bias evolution and the large inferred $`r_0`$ at $`z1`$ of massive ellipticals, if the assumption that galaxies are conserved quantities (i.e. closed-box systems) is relaxed. These hierarchical merging models (e.g. Mo & White 1996; Matarrese et al. 1997; Moscardini et al. 1998; Kauffmann et al. 1999b; Moustakas & Somerville 2002; Mo & White 2002, and references therein) prescribe that for certain types of objects bias can grow stronger with redshift than the growth of perturbations, resulting in a $`r_0`$ that is constant or even increasing with redshift.

In the (transient) model of Matarrese et al. (1997) it is assumed that the mass of the dark matter halo also determines the physical parameters of the galaxy that it contains. Based on the work of Mo & White (1996) and the formalism of Press & Schechter (1974), Matarrese et al. (1997) derive that the bias in such a model evolves as

$$b(z)=11/\delta _c+[b_0(11/\delta _c)]/D(z)^\beta ,$$
(14)
where $`\delta _c=1.686`$ is the critical linear overdensity for spherical collapse (but see also Lilje 1992). The parameters $`b_0`$ and $`\beta `$ depend on the minimum mass of the halo, and we have used the COBE-normalized ($`\mathrm{\Lambda }`$CDM) values for $`_{\mathrm{min}}=10^{1214}M_{}`$ given by Moscardini et al. (1998) to plot this model in Fig. 16 (thin dashed lines). We find that the $`_{\mathrm{min}}=10^{14}M_{}`$ model is able to fit the measurements at both $`z1`$ and $`z=0`$. Likewise, the model with $`_{\mathrm{min}}=10^{13}M_{}`$ has been found to fit the spatial clustering of QSOs relatively well (Croom et al. 2001), although several serious caveats exist (see Rengelink 1998; Croom et al. 2001). Most importantly, the assumption that there always exists a simple relationship between the mass of the dark matter halo and the property by which a galaxy is selected may not be valid.

In Fig. 16 we have also indicated the predicted evolution of the clustering of early-type galaxies (thick dashed line) from the $`\mathrm{\Lambda }`$CDM-models of Kauffmann et al. (1999a, b) (see also Somerville et al. 2001), normalised to $`r_0`$ found for local ellipticals. An important feature of the models presented in Kauffmann et al. (1999b) is that one naturally expects a dip in $`r_0`$ between $`z=0`$ and $`z1`$, if structure is probed by galaxies of intermediate luminosities residing in haloes of masses $`10^{1112}M_{}`$ that have formed early and are unbiased tracers of the overall mass distribution. However, these simulations also show that this dip is very sensitive to sample selection criteria: massive early-type galaxies exhibit no dip in clustering between $`z=0`$ and $`z1`$, because they occur in rare, very masssive haloes of $`10^{1314}M_{}`$ which are strongly biased locally, and which become even stronger biased with redshift. The agreement of this model with the results presented in this paper and the results of Daddi et al. (2001) and McCarthy et al. (2001) is striking. Although promising, some discrepancies between the model and the observations remain. For instance, Daddi et al. (2001, 2002) find strong disagreement between the model and the high observed space density of EROs, seemingly consistent with the purely passive evolution of local ellipticals. Furthermore, current merging models generally predict that these galaxies should have experienced recent star-formation activity, while this is not observed. It may become possible to still reconcile the observations with the $`\mathrm{\Lambda }`$CDM merging models if, for example, the merging is accompanied by little star-formation (Daddi et al. 2001). Also, the EROs are found to have relatively old stellar populations of $`3`$ Gyr that show no indications of recent formation processes. However, Moustakas & Somerville (2002) point out that the relatively old ages of their stellar populations do not automatically imply similar ages for the host galaxies.

Despite the success of current hierarchical models in predicting the evolution of bias for these massive galaxies, we would like to point out that galaxy conservation or linear/stable clustering evolution could still be able to explain the measurements if EROs and/or powerful radio galaxies are solely found in rich Abell-type clusters with (present-day) $`r_01525`$ $`h^1`$ Mpc. As we have shown there is substantial evidence that this may be the case for, at least, the powerful radio galaxies, and future data may show whether this also holds for (a subset of) the population of EROs.

At the highest redshifts, clustering of LBGs at $`3z4`$ indicate that these objects can be connected to local ellipticals in a galaxy conservation scenario. However, it is now believed that LBGs probably occupy much less massive halos of $`10^{1112}M_{}`$ than those that contain local massive galaxies, suggesting that if these objects are to be the progenitors of local ellipticals, they must have accumulated a considerable amount of mass (Adelberger 2000; Moustakas & Somerville 2002).

#### 6.2.3 Clustering and the occurrence of AGN at high $`z`$

Fig. 16 suggests that the clustering evolution of active galaxies in general is considerably different from that of ordinary galaxies. Albeit at a lower amplitude, the clustering of QSOs also shows a trend of constant or slightly increasing amplitude with redshift, very similar to the trend that we derive for the clustering of the most massive ellipticals. According to the standard paradigm, AGN are powered by the accretion of matter onto a (super-)massive blackhole (e.g. Rees 1984). This fuelling mechanism may very well be associated with the injection and accretion of gas during major merging events, and thus, the occurrence of AGN seems to be logically linked to the hierarchical scenarios for structure formation. Recently, in a series of papers (Haehnelt et al. 1998; Haehnelt & Kauffmann 2000; Kauffmann & Haehnelt 2000, 2002) the simulations of Kauffmann et al. (1999b) were extended to a unified model for the evolution of both galaxies and quasars. In their model, elliptical galaxies, supermassive black holes and starbursts are formed during major merging events, in which a fraction of the available gas is used to trigger quasar activity by accretion for about $`10^7`$ years, and the remaining gas is converted into stars in a single short burst. This model succesfully reproduces the evolution of cold gas that is derived from observations of damped Ly$`\alpha `$ systems, the luminosity functions and clustering properties of QSOs from the 2dF QSO survey, and the relation between bulge velocity dispersion and black hole mass that has been found in demographic studies of black holes in nearby galaxies (e.g. Kormendy & Richstone 1995; Magorrian et al. 1998; Gebhardt et al. 2000).

Although it has yet remained unknown exactly what processes cause the physical differences between radio-quiet and radio-loud AGN, recent results indicate that the hosts of all powerful AGN (both radio-loud and radio-quiet) are almost exclusively $`LL_{}`$ ellipticals (see Dunlop & McLure 2003, and references therein). However, the same studies also indicate that while radio-quiet AGN hosts can have black holes with masses of $`10^{610}M_{}`$, the radio-loud sources are cleanly confined to black hole masses $`M_{bh}5\times 10^8M_{}`$. Furthermore, in the regime of extreme radio luminosities that lie well beyond the FRI/FRII luminosity-break, the power needed can only be achieved by blackholes with $`M_{bh}>10^9M_{}`$, requiring host masses of $`>10^{12}M_{}`$ that imply $`L>L_{}`$ luminosities (Dunlop & McLure 2003). This may explain why the most powerful NVSS sources are extremely clustered compared to the, on average, less massive hosts of QSOs. This is supported by the fact that the radio sources in our lower radio luminosity bin have a correlation length similar to that of QSOs at $`z1`$, while both populations are still clustered more strongly compared to the field at $`z1`$. We conclude that the masses of the haloes, host galaxies, and black holes that are probed by the most powerful radio sources are among the most massive objects in the Universe, possibly formed through massive mergers in hierarchical fashion.

## 7 Summary

The main conclusions that can be drawn from our analysis are the following:
$``$ Below $`6\mathrm{}`$ $`w(\theta )`$ is dominated by the size distribution of multi-component radio sources. A simple model of the physical size distribution of FRII radio galaxies is able to explain the observed enhancement of the cosmological clustering signal.
$``$ The amplitude of the angular two-point correlation function of radio sources increases with increasing radio flux, corresponding to a similar increase in $`r_0`$ with increasing average radio power of the samples. This suggests that powerful FRII radio galaxies are intrinsically more strongly clustered than the average population of radio galaxies at $`z1`$. This is consistent with the extremely rich environments in which high redshift FRIIs are generally found.
$``$ The correlation lengths of powerful radio galaxies and EROs are of comparable magnitude and both are associated with massive ellipticals at $`z1`$. This suggests that we could be looking at identical objects at different stages of their evolution, implying that AGN activity is an important phase in the evolution of massive galaxies in general.
$``$ The evolution that we infer for the clustering of massive ellipticals between $`z1`$ and $`z0`$ is in agreement with predictions from hierarchical models for structure formation, because they can account for the observed lack of evolution in $`r_0`$. However, the large correlation length of powerful radio galaxies at $`z1`$ is also consistent with galaxy conservation models if they are primarily associated with rich, Abell-type clusters.

###### Acknowledgements.

We would like to thank Chris Blake, Emanuele Daddi, Matt Jarvis, Melanie Johnston-Hollitt and Jaron Kurk for productive discussions and reading through the text. We also thank the referee for very helpful comments.