format: text

text:
# Red galaxy overdensities and the varied cluster environments of powerful radio sources with 𝐳∼1.6

## 1 Introduction

Understanding the formation of large scale structure and the related formation and evolution of galaxies are two of the most important issues in modern–day astronomy. Rich clusters of galaxies at high redshifts ($`z\mathrm{}>1`$) provide an opportunity to investigate both of these problems directly: first, since clusters are the largest, most massive, collapsed structures in the Universe, discerning the epoch and process of their assembly places strong constraints on models of large scale structure formation (e.g. Bahcall & Fan 1998); second, because they contain large numbers of galaxies at the same distance, clusters provide a unique resource for investigating galaxy evolution through the evolution of galaxy scaling relations such as the fundamental plane (e.g. van Dokkum et al 1998) and colour–magnitude relations (e.g. Stanford et al. 1998); third, clusters contain the oldest galaxies known, and these can strongly constrain the first epoch of formation of ellipticals (cf Dunlop et al. 1996, Spinrad et al. 1997).

The identification of significant samples of clusters much above redshift one is both challenging and inefficient using optical or X–ray selection techniques. Existing X–ray surveys suffer from sensitivity limits (the ROSAT Deep Cluster Survey, for example, found very few clusters above $`z1`$, and none above $`z=1.3`$; e.g. Stanford et al. 2002) whilst the smaller field of view of the current generation of X-ray telescopes make them relatively inefficient for wide–area surveys. At optical wavelengths the contrast of a cluster above the background counts is minimal above redshift one, requiring very deep wide–area multi–colour surveys including near–infrared wavebands to identify a significant sample of $`z1`$ clusters; whilst certainly feasible, this is very expensive on telescope time. For this reason, an alternative method has frequently been adopted, using targetted studies towards powerful AGN. At low redshifts AGN are frequently found in relatively rich environments (e.g. Yates, Miller and Peacock 1989; Hill & Lilly 1991), and the host galaxies of high redshift AGN are amongst the most massive galaxies known in the early Universe (e.g. de Breuck et al. 2002), making these promising candidates for residing in galaxy overdensities.

The environments of powerful radio sources at redshifts $`z1`$ have been investigated by many authors over recent years using a variety of different techniques, including searches for luminous extended X–ray emission (e.g. Crawford & Fabian 1996; Fabian 2001), searches for over-densities of galaxies in infrared colours (e.g. Dickinson 1997, Hall & Green 1998) or emission line images (e.g. McCarthy et al. 1995), investigations of the colour–magnitude relations for red cluster ellipticals (e.g. Best 2000, hereafter B00), cross–correlation analyses (e.g. B00; Wold et al 2000) and direct spectroscopic studies of individual sources (e.g. Dickinson 1997, Deltorn et al. 1997). The evidence that at least some powerful $`z1`$ radio sources are located at the centres of galaxy overdensities is overwhelming (see B00 for a review). What remains unclear is the ubiquity, scale, and nature of these (proto?) clusters.

At higher redshifts, radio sources have been detected out to $`z=5.2`$ , and some well–studied sources have been spectroscopically confirmed to lie in overdense environments (e.g. 1138-215 at $`z=2.2`$; Kurk et al. 2000; Kurk et al. in preparation). An on–going VLT large project (PI: Miley) has found order–of–magnitude overdensities of Ly-$`\alpha `$ emitters around all 4 deeply studied $`z>2`$ radio galaxies to date (Venemans et al. in preparation), including the detection of the most distant known large structure around the $`z=4.1`$ radio galaxy TN J1338-1942 . The $`20`$ line emitting galaxies in the structures around each of these high redshift radio galaxies are found over regions extending to $`\mathrm{}>3`$ Mpc, and have velocity dispersions of 300–1000 km s -1(Venemans et al., private communication). Thus, it does indeed seem that powerful radio galaxies offer unique111Of course, within a few years deeper X–ray surveys and Sunyaev–Zel’dovich effect surveys will also find clusters at redshifts $`z1`$, if such clusters already contain a significant hot gas component. prospects for tracking galaxy clusters over cosmic time. The detection of overdensities at these redshifts, however, is based predominantly on observing emission–line galaxies near the radio galaxy redshift; the strong red sequences characteristic of nearby clusters are not observed towards these $`z>2`$ overdensities, restricting the evolutionary studies of cluster ellipticals using the fundamental plane and colour–magnitude relations from being extended to these high redshifts.

The redshift range $`1.2\mathrm{}<z\mathrm{}<2.0`$ is therefore a particularly interesting redshift range to study. This is the redshift range where cluster studies go beyond those currently being carried out using optical and X–ray selected clusters, whilst still being of low–enough redshift that red sequences of cluster ellipticals are likely to exist, and can be studied in a practical amount of observing time. Such clusters will also allow a link to be made between the ‘proto’-clusters being found around $`z>2`$ radio galaxies, and the rich virialised clusters seen in the nearby Universe. The elliptical galaxy population in clusters may well be undergoing a fundamental change within this redshift range

Hall & Green  have previously investigated the small–scale environments of a sample of 31 quasars with redshifts $`1\mathrm{}<z\mathrm{}<2`$, variably selected from previous surveys. They found an excess of K–band galaxy counts on two scales around the quasars; a pronounced peak in the galaxy counts within the central 40 arcsec radius surrounding the quasars was accompanied by a uniform galaxy overdensity across their entire fields (to about $`1.5`$ arcmins radius). Comparing R and K–band data they showed that the excess counts correspond to redder galaxies than the typical field population, with a significant overdensity of extremely red objects. This is indicative that a red elliptical sequence may be present in these fields. Hall et al.  investigated the physical extent of the large–scale overdensity in the fields of two of these AGN, and found a net excess of galaxy counts out to 140 arcsec ($`1.3`$ Mpc) in one of the two fields.

We have embarked upon a project to investigate in detail the environments of a complete sample of the most luminous radio sources (radio galaxies and quasars) within the narrow redshift range $`1.44z1.7`$. The aim of this project is to use deep, wide–area multi–colour imaging to study in detail any galaxy overdensity, in particular the nature of the excess galaxies, the presence of a red sequence of cluster ellipticals, and the radial extent of the overdensity. By studying a complete AGN sample it will be possible to investigate the variety of environments found, and whether the environment of a radio–loud AGN is in any way correlated with the radio source properties. The candidate cluster galaxies in the richest overdensities can then be followed–up using multi–object spectroscopy in order to extend the galaxy evolution studies.

This paper provides the initial results of this project, focussing on the imaging data and results from the first 6 clusters in the sample. Specifically, the lay-out of the paper is as follows. In Section 2 the sample selection, observations, and data reduction are described. The resulting multi–colour images of each field are provided and described in Section 3, and galaxy catalogues are constructed from these. The galaxy number counts and the colour distribution of the galaxies are investigated in Section 4, and a cross–correlation analysis carried out. The properties of the different fields are related to the properties of their AGN in Section 5. The results are discussed and conclusions drawn in Section 6. Throughout the paper, the values adopted for the cosmological parameters are $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`H_0=65`$ km s-1Mpc-1.

## 2 Observations

### 2.1 Sample Selection and Observations

The radio sources were drawn from the equatorial sample of 178 powerful radio sources defined by Best, Röttgering and Lehnert (1999; the BRL sample). This sample was designed to be equivalent to the northern 3CR sample, and was drawn from the Molonglo radio survey at 408 MHz , according to three criteria: (i) the sources must have declinations $`30^{}<\delta <+10^{}`$; (ii) they must lie away from the galactic plane, $`|b|>10`$; (iii) they must be brighter than $`S_{408\mathrm{M}\mathrm{H}\mathrm{z}}=5`$ Jy (roughly equivalent to the definition of the 3CR sample). Spectroscopic redshifts are available for 177 of these 178 sources .

From the BRL sample, the 9 radio sources with redshifts $`1.44z1.7`$ and galactic latitudes $`|b|>25`$ were selected. This sample comprised 4 radio galaxies and 5 radio loud quasars. Of these, the environments of 6 radio sources have been studied to date, and are presented in this paper. These 6 radio sources, limited by the constraints of telescope time, correspond to all of those with $`15<\mathrm{RA}<02`$ and therefore constitute a complete subsample of the most powerful radio sources at these redshifts.

The aim of the current observations was to deeply image these fields in the $`R`$, $`J`$ and $`K`$ wavebands. This filter combination was chosen to provide maximal sensitivity to old cluster ellipticals: the strongest continuum features of such ellipticals is the 4000Å break, which at $`z1.5`$ falls between the $`R`$ and $`J`$ bands, making this colour very sensitive to old cluster galaxies. The addition of the $`K`$–band data provides a long red colour baseline to sample the old stellar populations of these ellipticals, whilst also facilitating identification of bluer cluster galaxies with on-going star formation (and hence less pronounced 4000Å breaks), and helping to filter out non–cluster objects with very red $`RJ`$ colours, such as dusty starbursting galaxies. At $`z1.5`$ passively evolving old ellipticals would have colours of $`JK1.8`$ and $`RJ3.7`$, and brightest cluster galaxies at these redshifts are expected to have magnitudes of $`K17.5`$ to 18 (cf. those of radio galaxies at these redshifts, e.g. Best et al. 1998). The goal of detecting elliptical galaxies in all colours to $`3`$ magnitudes fainter than this set the required depths of the different observations. Full details of the observations are provided in Table 1.

The six fields were observed in the $`J`$ and $`K`$–bands using the infrared imager SOFI on the ESO New Technology Telescope (NTT) during September 2000 and August 2001. SOFI has a 1024 by 1024 HgCdTe array with 0.29 arcsec pixels, providing a field of view of 5 by 5 arcmins. Images were taken using a jittering technique, with the pointing centre moved randomly within a 40 arcsec square box between each exposure. At each pointing position, 1–minute of data was taken, split into six 10s co-adds in the K–band and three 20s co-adds in the J–band. In the $`K`$–band the $`K_s`$ filter was used to reduce the sky background; the $`K_s`$ magnitude is related to the $`K`$-magnitude by $`KK_s0.005(JK)`$, and all quoted magnitudes have been converted to the $`K`$–band magnitude. Between 1.7 and 2.5 hours of data were obtained for each source (see Table 1 for details), reaching a ‘limiting magnitude’ of $`K20.6`$; this value corresponds to the 50% completeness limit for galaxies with sizes typical of those at redshifts $`z1.6`$. In the $`J`$–band the total exposure time varied between 1 and 2 hours, reaching a limiting depth of $`J22.4`$.

For 5 of the 6 fields (2128$``$208 was randomly excluded due to telescope time constraints), $`R`$-band data for the field were taken through the Bessel $`R`$ filter using SUSI-2 on the NTT. SUSI-2 comprises two 2048 by 4096 pixel EEV CCDs, with a pixel scale of 0.08 arcsec, which combined provide a 5.5 by 5.5 arcminute field of view. The two CCDs are separated by about 10 arcseconds on the sky; the total exposure time was therefore separated into 5-min blocks, with the pointing position moved through a 9-point pattern of offsets \[(0,0), (10,10), (20,-10), (30,0), (40,10), (-40,-10), (-30,0), (-20,10), (-10,-10) arcsec\] in order to attain equal sensitivity over the entire SOFI field. This pattern was repeated until between 3.75 and 6 hours of data had been obtained on each field, reaching a limiting depth of $`R25.9`$.

Standard stars were observed approximately every 2 hours throughout the nights to allow the images to be photometrically calibrated. For the SUSI-2 images, once each night the image was rotated by 180 degrees in order to place the standard stars on the opposite CCD; no significant differences in calibration were apparent between the two CCDs. All of the nights were photometric, and the instrument zero-points varied by less than 0.05 magnitudes within each of the Sept 2000 and Aug 2001 runs, although the $`R`$–band zero point changed by 0.1 magnitudes between the two runs.

### 2.2 Data reduction

The infrared observations were reduced in the standard manner, using the dimsum222dimsum is the ‘Deep Infrared Mosaicing Software’ package developed by Eisenhardt, Dickinson, Stanford and Ward. package within iraf. After subtraction of a dark frame, and removal of known bad pixels, a median of the images (after masking of bright objects) was constructed to obtain an accurate sky flat–field. The flat–fielded images were then sky subtracted, cosmic ray events were removed, and then the individual frames were block replicated by a factor of 2 in each dimension to allow more accurate alignment, resulting in a final image scale of 0.144 arcsec per pixel. The frames were accurately registered using the peak positions of several bright unresolved objects visible on all images, and were combined. Because of the long exposure times, the seeing sometimes varied throughout the dataset, and so during the combination process the images were weighted according to their seeing in order to provide a final image with the optimal combination of seeing and sensitivity. Stars from the US Naval Observatory catalogue A2.0  were used to provide an astrometric solution for the images. The absolute astrometry should be accurate to better than 0.5 arcsec. The fields were photometrically calibrated onto the standard Vega photometric system.

The $`R`$–band images were reduced within iraf. The images were bias subtracted, and then flat–fielded using a flat–field constructed by from the median of the various images (after masking bright objects). An illumination correction was then applied to remove large–scale variations in the flat–fielded images. The images were registered using the peak positions of several bright unresolved objects, and then they were combined, weighted according to their seeing as for the $`J`$ and $`K`$ band data. Cosmic rays were rejected during the combination process. The resultant image was then astrometrically transformed using a 5 parameter fit (RA and Dec of the central pixel, pixel scale in x and y directions, and rotation angle) to align it with the infrared data, converting it to 0.144 arcsec pixels whilst conserving the object magnitudes. The relative astrometry of the IR and optical frames is accurate to $`<0.1`$ arcsec.

## 3 The radio source fields

K-band images of the six fields are shown in Figures 1 to 3, with the radio source host galaxy marked by a cross in each case. Three–colour images of the 100 by 100 arcsec region around the radio source host in each field are provided in Figures 4 to 6. A brief discussion of each individual field follows.

0000–177 ($`𝐳=\mathbf{1.47}`$): This radio–loud quasar lies in a relatively sparse region of sky. A faint object is seen 1-2 arcsec to the north–east of the quasar host, and a further seven objects with $`RK>4`$ lie within 30 arcsec ($`275`$ kpc at $`z1.5`$), indicating a possible compact group of galaxies at the quasar redshift. There is no obvious large–scale galaxy overdensity.

0016–129 ($`𝐳=\mathbf{1.59}`$): This radio galaxy is the westernmost of a string of three luminous red galaxies in the centre of this field, and shows extended blue emission along the radio axis direction, characteristic of the alignment effect in high redshift radio galaxies (cf. McCarthy et al. 1987, Chambers et al. 1987). The redder galaxies in the field are more concentrated towards the central regions of the field, within 1–2 arcmins of the AGN.

0139–273 ($`𝐳=\mathbf{1.44}`$): This is the richest field of those imaged in the current observations, and contains a large number of very red galaxies. These galaxies are found predominantly within the inner regions of the field, although no red galaxies are seen close to the AGN. Very luminous elongated blue emission extends from the radio galaxy host, along the direction of the radio source axis.

1524–136 ($`𝐳=\mathbf{1.69}`$): This radio loud quasar is the highest redshift AGN in the current sample, and its surrounding field is relatively rich in the K–band. A moderate excess of red galaxies is seen on scales out to about 120 arcsec ($``$ 1 Mpc). No sharp central peak (such as that seen for 2025$``$155) is seen, suggesting that any cluster environment may still be forming.

2025–155 ($`𝐳=\mathbf{1.50}`$): A very compact group of several red galaxies is found within 15–20 arcsec (150 kpc) radius of this radio–loud quasar. All of these objects have colours consistent with old passive ellipticals at the radio source redshift. A significant overdensity of red galaxies is found extending over the entire 5 arcmin field, although predominantly within the inner regions.

2128–208 ($`𝐳=\mathbf{1.62}`$): This field does not contain a particular excess of K-band galaxies and, except for a companion object to the quasar, no other galaxy is found within 15 arcsec of the AGN host. With only J and K–band data available, there is little information about the colour distribution of the galaxies across the field.

It is clear from this direct examination of the different images that a wide variety of environments in found amongst the sample. At one extreme, the field of 2128$``$208 shows little evidence for any net overdensity of red galaxies on either small ($`\mathrm{}<150`$ kpc) or large ($``$ Mpc) scales, whilst in contrast 2025$``$155 has a very rich environment on small scales coupled with an apparent excess of red galaxies across the entire field. Between these extremes, the other fields show evidence for structure on differing scales: a small scale excess of red galaxies is found around 0000$``$177, but not much large scale structure; the 0016$``$129 field has three central red galaxies and a weak larger–scale excess; above average numbers of red galaxiess are seen on both moderate and large scales around 1524$``$136; and the field of 0139$``$273 shows a strong excess of red galaxies on larger scales, but with few red objects close to the AGN. In the discussions of the statistical properties of the fields that follow, it is important to bear in mind that global results quoted correspond to the average values across 6 rather different environments.

### 3.1 The source catalogues

Source extraction and catalogue construction followed the method described in detail by B00, and the reader is referred to that paper for full details. This method is summarised briefly below.

Throughout these analyses, the K–band data were used as the primary dataset from which the catalogue is defined. Source detection and photometry were carried out on the K–band frames using sextractor version 2.2.1 , using the output exposure map produced by dimsum as a weight map for sextractor to compensate for the varying noise levels associated with the jittering procedure. 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 15 (equivalent to 4 pixels prior to the block replication during the mosaicing procedure). A search for negative holes using the same extraction parameters (cf B00) found only one negative feature over all 6 frames brighter than the 50% completeness limit, implying that essentially all of the extracted objects down to that level are real. The sextractor catalogues were examined carefully, and obvious problems, such as objects associated with the diffraction spikes of bright stars, were corrected. Those objects within the central 2000 by 2000 pixels (290 by 290 arcsec) of each image were retained for the subsequent analysis; more distant objects lie in the regions of the image where the noise levels are higher due to the jittering process, and were discarded; at the edges of the retained regions of the fields, the noise level is only $`10`$% higher than in the field centre. Small regions of sky (of 10, 20 and 10 arcsecond radius) were also discarded around the three bright stars in the images of 0016$``$129, 1524$``$136 and 2025$``$155. The six fields combined provide a sky area of 140 square arcminutes.

sextractor’s mag\_best estimator was used to determine the magnitudes of the extracted sources; this yields an estimate for the ‘total’ magnitude using Kron’s  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. sextractor isophotal apertures smaller than 1 arcsecond were replaced by circular 1 arcsecond apertures. The determined magnitudes were corrected for galactic extinction using the $`E(BV)`$ values provided in the NASA / IPAC Extragalactic Database (NED). The accuracy of these total magnitudes and the completeness level of the source extraction were investigated as a function of position on each image through a Monte–Carlo simulation, as fully described by B00. Briefly, a total of 10000 stellar objects and 25000 galaxies (see B00 for the prescription by which these were constructed) were added, 25 at a time, to each image with a magnitude between 15 and 22.5, and were then re-extracted using sextractor with the same input parameters as for the original source extraction, to see if the added objects were detected. For each recovered source, the difference between the total magnitude measured by sextractor and the input total magnitude was determined.

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 mean 50% completeness limits are $`K=20.55`$ for galaxies and $`K=21.05`$ for stars (see Figure 7); the galactic completeness limits for each individual field are provided in Table 1. 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 8. This difference was used to correct the observed magnitudes and obtain true total magnitudes.

The $`J`$ and $`R`$ band magnitudes for the catalogue sources were then measured by running sextractor in its double image mode, where one image (the $`K`$ band) is used to detect the objects and the fluxes and magnitudes are then measured from the second ($`J`$ or $`R`$ band) image through the same aperture. Those objects whose $`J`$ or $`R`$ band fluxes were measured to have a significance of less than 3$`\sigma `$, where $`\sigma `$ is the uncertainty on the flux measurement provided by sextractor, were replaced by 3$`\sigma `$ upper limits. The sextractor 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 blank regions of 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 in all three colours.

### 3.2 Star–galaxy separation

sextractor provides a ‘stellaricity index’ for each object, where in the ideal case a galaxy has a stellaricity index of 0.0 and a star has 1.0; low signal–to–noise, and galaxies more compact than the seeing, lead to an overlap in the calculated stellaricity indices for the two types of object at the faintest magnitudes, and hence uncertain star–galaxy separation. Stellaricity indices were taken from the K–band data. At magnitudes $`K\mathrm{}<18`$, star–galaxy separation could be carried out directly from these stellaricity values taking values above 0.8 (the generally adopted cut-off) to correspond to stars. At fainter magnitudes, separation is less certain, but was aided by the addition of the multi–colour information; by plotting the R-K versus J-K colour–colour diagram as a function of stellaricity index, a region of colour space typically occupied by stars was defined (see Figure 9), and all objects in this colour space region with stellaricities above 0.5 were also classified as stars. This classification is somewhat conservative, and a small proportion of stars may be included along with the galaxies, but this was considered preferable to falsely classifying galaxies as stars.

Figure 10 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.201K0.128`$. This relatively good fit to the differential star counts using a single power–law distribution demonstrates that the star–galaxy separation is working well to at least $`K20`$. All objects classified as stars were removed from further analysis.

## 4 Galaxy Counts and Colours

### 4.1 Galaxy counts

The galaxy counts were determined by once again following the procedure of B00. In brief, the measured magnitudes of the galaxies were converted to true ‘total’ magnitudes using the offsets determined from the simulations above, and binned in 0.5 magnitude bins over the magnitude range $`14<K<21`$. These raw galaxy counts, $`n_{\mathrm{raw}}`$, were then corrected for completeness and for magnitude biasing effects (again as calculated from the results of the simulations) to produce corrected galaxy counts ($`n_\mathrm{c}`$). Scaling by the observed sky area produced final counts per magnitude per square degree, $`N_\mathrm{c}`$. In calculating the uncertainty on this value, in addition to the Poissonian error on the raw galaxy counts, an extra error term was included to account for the uncertainty of the completeness correction procedure; this was estimated at 20% of the number of galaxies added in the correction. The resulting galaxy counts are tabulated in Table 2 and are plotted in Figure 11 compared with counts from various K–band blank–field surveys .

The derived counts are fully consistent with the literature counts at magnitudes $`K\mathrm{}<17.5`$. This magnitude corresponds to the expected magnitude of brightest cluster galaxies at $`z1.5`$. Fainter than this magnitude an surplus of galaxies is seen; the galaxy counts still lie within the scatter of the other observations, that scatter being largely dominated by cosmic variance (cf Daddi et al. 2000a), but the counts in every magnitude bin from $`17.5<K<20.5`$ exceed the literature average by at least 1$`\sigma `$.

In Table 3 the total galaxy counts in the magnitude range $`17.5<K<20.5`$ are provided for each individual field, and compared to the mean value of the blank–field surveys. To carry out such a comparison it is important to bear in mind that the uncertainty on the counts in the AGN fields will be significantly larger than solely the Poissonian error, due to cosmic variance and the effects of galaxy clustering. The true uncertainty on the counts can be estimated from the equation $`\sigma _{\mathrm{true}}^2=N(1+NAC)`$, where $`N`$ is the total counts, $`A`$ is the angular cross–correlation amplitude, and $`C`$ is the integral constraint which depends upon the size and geometry of the field. As derived in Section 4.6, for these fields the values of $`A`$ and $`C`$ are 0.0005 and 16.7 respectively. Hence, for $`N300`$ galaxies detected in each field, the uncertainty on the counts is approximately double that which would be expected from Poissonian statistics alone. The errors tabulated in Table 3 are these true errors, including this contribution of cosmic variance. For the literature counts, given that these are derived by combining many independent surveys over different areas of sky, the combination should average out the cosmic variance effects to a large extent.

It is noteworthy that five of the six fields have counts above the literature average, with three of the fields being $`\mathrm{}>2\sigma `$ above. Combining the fields, the cumulative counts over this magnitude range are 3$`\sigma `$ above the average of blank–field surveys; it should be emphasised that these quoted significance levels including the effects of cosmic variance, and are twice this size considering Poisson statistics alone. Indeed, since the cosmic variance is caused by the clustering of galaxies, fields which contain clusters of galaxies, like these are proposed to, are in fact the type of field that contributes significantly to this cosmic variance, and so the Poissonian estimate may be more appropriate.

The cumulative excess of galaxy counts corresponds to about 290 ($`\pm 100`$) additional galaxies, that is, of order 50 galaxies per field. Although it is possible that some fraction of these correspond to faint stars that were misclassified as galaxies (cf. the fall in star counts in the final bin of Figure 10), it is clear that these fields are, on average, richer than average environments. For comparison, in Abell’s original cluster definition  he considered the excess counts of galaxies within a radius of 1.5$`h^1`$ Mpc of the cluster centre, with magnitudes between $`m_3`$ and $`m_3+2`$, where $`m_3`$ is the magnitude of the third ranked cluster member: an excess of 30 to 49 galaxies corresponded to Abell Class 0 clusters, and an excess of 50 to 79 galaxies was Abell Class 1. Assuming that the radio host galaxies are of order a magnitude brighter than a third–ranked cluster member333If the radio galaxies were dramatically brighter in the K-band than other cluster members, as is the case for some cD galaxies in nearby clusters, then that would require the fields to be very rich, probably implausibly so. the average excess number counts correspond to moderate richness clusters.

To obtain a better indication of the cluster richnesses, the parameter $`N_{0.5}`$ defined by Hill and Lilly  can be used. This is an Abell–type measurement defined as the net excess number of galaxies within a projected radius of 0.5 Mpc (at the cluster redshift) of the central galaxy with magnitudes between $`m_1`$ and $`m_1+3`$, $`m_1`$ being the magnitude of the brightest cluster galaxy444For the quasars, the value of $`m_1`$ was estimated from the redshift using the Hubble K$`z`$ relation for radio galaxies of the same radio power (e.g. Lilly & Longair 1984, Best et al. 1998); this is appropriate due to the similarity of radio galaxy and quasar hosts at high redshifts (e.g. Lehnert et al. 1992, Kukula et al. 2001).. The value of $`N_{0.5}`$ was calculated for each field, in each case using the region of the image at projected distances greater than 1 Mpc distant from the AGN to estimate the background counts, thus accounting to a large an extent as possible for the cosmic variance in the background sources.

The values of $`N_{0.5}`$ derived are provided in Table 3. Although the errors on each individual value are large, it is notable that the three values significantly above zero correspond to the three fields with large excesses in their number counts. This implies that the number count excesses are predominantly in the inner regions around the AGN, and that these are therefore consistent with cluster environments. Hill and Lilly (and see also Wold et al. 2000) calculated a conversion between the value of $`N_{0.5}`$ and the traditional measure of Abell richness. $`N_{0.5}10`$ corresponds to the an Abell richness 0 cluster, whilst $`N_{0.5}20`$ corresponds to Abell richness 1. Three of the 6 fields studied are compatible with being relatively rich clusters, whilst the other three seem to be at most of Abell class 0 richness. It should be noted, however, that conclusions based upon galaxy counts alone are not robust.

### 4.2 Galaxy colours

The presence of red galaxies associated with a cluster ought to redden the average galaxy colours with respect to blank–field expectations, and indeed this is seen in the mean or median colours of each field (see Table 3). Interestingly, the average colour in the radio galaxy fields is typically reddest in the fields with the lowest number counts, whereas naively the opposite result would be expected if the excess number counts were solely associated with typically red cluster galaxies. This result may imply that there are excess red galaxies around all of the AGN, and in the cases where the non-cluster (foreground and background) galaxy counts are randomly lower, the red cluster galaxies are therefore a larger proportion of the total counts and influence the average colours more. This should offer further caution about interpreting the cluster environments based on counts in a single colour only.

To investigate in more detail the nature of the excess galaxies, the R$``$K colour distribution of the galaxies was derived as a function of projected distance from the AGN for all galaxies with K–magnitudes in the range $`16<K<21`$ across the five fields (excluding 2128$``$208 for which $`R`$–band data is lacking). Four radial bins were considered (0–100, 100–250, 250–500 and 500–1000 kpc, as evaluated at the redshifts of the AGN), and the colour distribution of the galaxies in these radial ranges are shown in Figure 12. Also shown as a dotted line in these plots is the colour distribution of those sources beyond 1 Mpc radius from the AGN, constructed as a comparison sample. There is a clear surplus of red ($`RK\mathrm{}>4`$) galaxies at all radii $`<1`$ Mpc compared to the outer ‘blank–field’ regions. This is most pronounced in the inner 100 kpc radius, but is very significant all the way out to 1 Mpc.

This colour distribution was also derived as a function of K–band magnitude, both over the total combined area of the five images and over just the inner quadrant of the images, corresponding to distances out to about 700 kpc from the AGN. These distributions are displayed in Figure 13. The median colours of the galaxy distributions in each magnitude range can be found in Table 4. To compare these, a blank–field distribution was also determined by combining the data of the CADIS 16hr field (Thompson et al. 1999; data kindly supplied by Dave Thompson) and the K20 survey (Cimatti et al. 2002; data kindly supplied by Andrea Cimatti). The data from the CADIS survey were colour corrected to account for differences in the $`R`$ and $`K`$ passbands according to the relation $`RK=(RK^{})_{\mathrm{CADIS}}+0.07`$. The distributions of the CADIS and K20 data were then well matched. These two fields cover a combined sky area of 206 square arcminutes, and are complete for $`K<20`$. For $`20<K<21`$, a blank–field galaxy $`RK`$ distribution was constructed using the compilation of Hall & Green , colour correcting the passbands according to $`RK=(rK)_{\mathrm{HG}}0.26`$ . These expected blank–field distributions are indicated by the shaded regions in Figure 13 and the median colours tabulated in Table 4.

At magnitudes $`K<18`$, the colour distribution of the total fields and the central regions are both well matched to the field distribution, with comparable median colours. By $`18<K<19`$ the median colour begins to become redder in the AGN fields than in the blank fields, particularly in the central regions, due to a small excess of galaxies with $`RK\mathrm{}>5`$. This red excess becomes even more pronounced in the two fainter bins. The combined colour distribution for all galaxies with $`16<K<21`$ clearly shows that the higher number counts in these fields are associated with red galaxies: the blank–field predictions are matched very well for galaxies bluer than $`RK=3`$, but the AGN fields show a strong surplus of redder galaxies, particularly those with $`RK>5`$. These results are even more pronounced when only the central regions of the images, around the AGN, are considered.

Combining the results for the colour distributions of the galaxies as a function of both radius and magnitude, it is clear that the predominant excess galaxy population corresponds to a population of red ($`RK\mathrm{}>45`$) galaxies, with magnitudes $`K\mathrm{}>18`$, lying at projected distances within $`1`$ Mpc from the AGN. These are exactly the typical magnitudes, colours and locations that would be expected of elliptical galaxies in a cluster environment at the AGN redshifts, if such ellipticals form at early cosmic epoch and evolve relatively passively. By contrast, in semi–analytic models of galaxy formation, very few ellipticals this red and this luminous would be expected in the field at redshifts $`z1.6`$ (cf. Daddi et al. 2000b).

### 4.3 Radial distribution of individual fields

Knowing the magnitudes and colours of the galaxies producing the excess counts, it is possible to investigate their radial distribution in each of the galaxy fields separately. Considering only those galaxies with magnitudes $`17.5<K<20.7`$ and colours $`RK>4`$, the radial distribution of the galaxy counts was derived for each field. These results are shown in Figure 14, together with a combined histogram produced by averaging across all five fields. Also included on these plots is the expected number density for blank-field sources of these magnitudes and colours; this was calculated by taking the K20 and CADIS surveys, making the same cuts in colour and magnitude as for the AGN fields, and then convolving the results with the completeness versus magnitude function (Figure 7).

The results from combining the five AGN fields are in agreement with those derived by Hall & Green , namely that a pronounced peak of red galaxies is seen within the inner $``$150 kpc, together with a weaker large–scale excess. With the larger field–of–view of the data provided here, it is now possible to quantify the radial extent of this large–scale overdensity: the galaxy counts remain above field expectations out to between 1 and 1.5 Mpc radius.

Comparing the five individual fields in Figure 14, however, it is immediately apparent that this combined result is the average of five fields with very different properties. 0016$``$129 and 2025$``$155 do both show radial distributions similar to the average result, with pronounced peaks in the inner bin followed by weak red galaxy overdensities to $`1`$ Mpc. 0000$``$177 has a central pronounced peak of red galaxies, but shows no large–scale excess (indeed quite the opposite); this may be indicative of a compact group of galaxies. 0139$``$273 and 1524$``$136 are both much less centrally concentrated, with smoother central peaks and much more gradual radial declines in red galaxy surface density out to the blank–field expectations at radial distances of between 1 and 1.5 Mpc. There is clearly a large range in the nature, scale and significance of any galaxy overdensities around these AGN.

That the radial distribution of galaxies within these large–scale overdensities is not more sharply peaked should not be surprising. A galaxy cluster with a mass of a few $`\times 10^{14}M_{}`$ and an initial scale of 3–5 Mpc has an overall density of 1–4 $`\times 10^{25}`$ kg m-3. The time that such a structure will take to separate out from Hubble flow and collapse is $`t(3\pi /32G\rho )^{1/2}`$ 3–6 Gyr . The Universe at redshifts $`z1.5`$ is about 4.5 Gyr old, and so such structures will still be collapsing into a virialised cluster such as those we see today; the large–scale structures seen around the AGN might even not yet be gravitationally bound.

### 4.4 Bright red galaxies and sub-structure

At high redshifts most clusters observed to date have several giant galaxies within them, presumably due to merging subclusters. Further, there is some evidence that AGN host galaxies are not always the only bright galaxies in their environment, for example due to sub-cluster mergers or being part of a group falling into a larger structure (e.g. Bremer, Baker & Lehnert 2002; Simpson & Rawlings 2002). Therefore, the fields have been searched for any galaxies other than the AGN that are brighter than $`K=18`$ and satisfy the colour criteria of the ‘cluster elliptical candidates’ above. The images around these bright red galaxies have then been examined; the majority of these galaxies appear to be fairly isolated, and so are likely to be foreground extremely red objects with redshifts of $`z1`$, since such foreground objects will tend to dominate the extremely red galaxy population at brighter magnitudes (e.g. Cimatti et al. 2002). However, in two cases these bright red galaxies have associated fainter red galaxies, and may represent substructures at the AGN redshift.

The first of these cases is in the field of 0000$``$177, where three galaxies with $`K=17.41`$, $`K=17.77`$ and $`K=17.85`$ and colours $`RK>5`$ are found close together, along with fainter red galaxies, about an arcminute east of the AGN. Indeed, these galaxies represent the small peak in the red galaxy radial distribution (Figure 14) at 600 to 700 kpc radius around this source. This region of the field is of comparable richness in red galaxies to the regions around the AGN. The second case is a group found around a $`K=17.09`$ galaxy about 1.5 arcminutes to the south–west of 0016$``$129, partially causing the 800–1000 kpc bump in the radial distribution around that source. This appears even richer than the red galaxy distribution around the radio galaxy itself, which might explain why the value of $`N_{0.5}`$ for this source is fairly low even though the number counts are well above the literature average.

It is not clear whether these extra structures are indeed at the same redshifts as the AGN, but if so then this would indicate that any large–scale structure around the AGN is still in the process of building up. It is interesting that this evidence for substructure is found in the two cases with the weakest evidence for red galaxy excesses in the 200–500 kpc radius range.

### 4.5 The morphology–density relation

In nearby clusters, a well–known morphology–density relation exists , whereby in regions of high galaxy density, such as the centre of the cluster, higher proportions of early–type galaxies are found. Dressler et al.  investigated the evolution of this morphology–density relation with redshift, out to $`z0.5`$, and find that this is still present at these redshifts but only in high-concentration, regular clusters; no correlation between galaxy–type and local galaxy density is found for low-concentration, irregular clusters. They interpret this as meaning that the mechanisms that produce morphological segregation (e.g. dynamical friction, or conversion of late–types into early types by gas stripping and halting star formation) work most rapidly in high mass systems.

The nature of the small–scale galaxy excess around these distant radio sources (and in the two potential substructures) is particularly interesting because, where present, it is comprised almost entirely of red galaxies (cf. Figs 4 to 6). This indicates that morphological segregation may already have taken place in the very inner regions around the AGN, implying that the morphology–density relation is imprinted into cluster centres at a very early epoch.

### 4.6 Angular cross–correlation analyses

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 ($`\delta P`$) 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.

The value of $`w(\theta )`$ can be evaluated for all galaxy–galaxy pairs, for just AGN–galaxy pairs, or restricting the galaxies included to certain ranges of colours and magnitudes. Here, four different set of galaxy magnitudes and colours were considered: (1) ‘all galaxies’; (2) ‘blue galaxies’; (3) ‘red galaxies’; and (4) ‘cluster elliptical candidates’, which are the subset of the red galaxies that have K magnitudes, $`RK`$ and $`JK`$ colours in the range that would be expected for old passive cluster galaxies. The exact colour and magnitude definitions of each class of objects, and the number of galaxies contained in that class, are provided in Table 5.

For each of these subsamples of galaxies the value of $`w(\theta )`$ was estimated for a variety of bins in $`\theta `$, following the method described by B00, considering both all galaxy–galaxy pairs and only AGN–galaxy pairs. The results obtained over all 5 fields with multi–colour data were averaged to produce overall results, with an estimate of the uncertainty in the value of the combined $`w(\theta )`$ in each bin provided by the scatter in the values between the different fields. $`w(\theta )`$ usually shows a power–law form, $`w(\theta )=A(\theta /\mathrm{deg})^\delta `$, with a canonical value for the exponent of $`\delta =0.8`$; this value appears also to be appropriate at high redshifts , and so was adopted here for all fits here. This does not always provide an excellent fit to the data, but the current data do not have sufficient signal–to–noise to warrant an attempt to independently fit this parameter; note that because this parameter is fixed, the formal errors on the cross-correlation amplitudes may be slightly higher than those quoted.

The finite area of the images leads to a small correction, known as the integral constraint ($`C`$), such that $`w(\theta )_{\mathrm{obs}}=A[(\theta /\mathrm{deg})^{0.8}C]`$. The value of $`C`$ can be estimated by integrating $`w(\theta )`$ over the area of each field, and for the current data corresponds to $`C_{\mathrm{gg}}=16.7`$ and $`C_{\mathrm{ag}}19.1`$ for all galaxy-galaxy pairs and just AGN-galaxy pairs respectively. Using these expressions, the amplitudes of the angular cross–correlation function were derived for both all galaxy–galaxy pairs and just AGN–galaxy pairs for each of the selected colour–magnitude subsamples of galaxies. The cross–correlation functions are displayed in Figure 15, and Table 5 provides the derived cross–correlation amplitudes.

The angular cross–correlation amplitude of ‘all galaxies’ in these fields is within the range of the values recently derived to $`K=20.75`$ for ‘blank–fields’ (e.g. Carlberg et al. 1997; Roche et al. 2002; Daddi et al. 2000a). However, comparing the values of $`A_{\mathrm{gg}}`$ and $`A_{\mathrm{ag}}`$ in Table 5 shows that the galaxies are much more highly clustered around the AGN than around other galaxies in general, indicating that a significant proportion of the galaxies are connected with the AGN. Evidence for such clustering around the AGN is particularly apparent when comparing the results for galaxies of different colours. The blue galaxies show no evidence for clustering, not even around the AGN, whilst red galaxies are known to be more highly clustered (e.g. Roche et al. 1996, Shepherd et al. 2001 and references therein) and this is reflected in the value of $`A_{\mathrm{gg}}`$ derived for this class; their higher value of $`A_{\mathrm{ag}}`$ shows that these are preferentially clustered around the AGN. The ‘cluster elliptical candidates’ have even larger values of $`A_{\mathrm{gg}}`$ and $`A_{\mathrm{ag}}`$, suggesting that this selection criteria picks out a highly clustered class of object, many of which are indeed likely to be old passive cluster ellipticals. It should be noted that this clustering signal is dominated by the points at small radii, $`\theta \mathrm{}<10`$ arcsec, corresponding to the excess in the upper left panel of Figure 12; however, even if the inner two bins of the cross–correlation analysis are ignored, the clustering amplitude around the AGN is still higher than that of the field in general, indicating that it is not only these nearby galaxies which produce the clustering signal.

If a population of $`N_\mathrm{c}`$ cluster galaxies, with a high intrinsic angular cross–correlation amplitude ($`A_\mathrm{c}`$), is observed together with a sample of $`N_\mathrm{i}`$ foreground or background interloper galaxies with zero clustering amplitude, then the resultant combined sample will have an observed cross–correlation amplitude of $`\frac{N_\mathrm{c}}{N_\mathrm{c}+N_\mathrm{i}}A_\mathrm{c}`$. In other words, the observed amplitude will be the amplitude that would be found for the cluster galaxies alone, scaled down by the proportion of the total sample that are cluster galaxies. This situation is a good approximation to those of the ‘red galaxy’ and ‘cluster elliptical candidate’ samples, since the observed value of $`A_{\mathrm{gg}}`$ for these samples are almost an order of magnitude below those of $`A_{\mathrm{ag}}`$. Comparing the values of $`A_{\mathrm{ag}}`$ then indicates that the colour definition for the ‘cluster elliptical candidate’ sample approximately halves the fraction of interlopers relative to the red galaxies.

If the results from the number count analysis are taken at face value, then there are about 270 clusters galaxies across the five fields studied in the cross–correlation analysis (290 excess galaxies determined in Section 4.1, of which 20 are from the 2128$``$208 field, not considered here); using this, an estimate can be made for the number of cluster galaxies in each of the colour–selected samples. Assuming that about 80% of the cluster galaxies are ‘red’ (from Figure 13) then within the ‘red galaxy’ sample there are about 220 cluster galaxies, which corresponds to about 25% of the overall ‘red galaxy’ sample (see Table 5). This would suggest that about 40% of the ‘cluster elliptical candidate’ sample, that is $``$100 galaxies, are indeed cluster galaxies. If these results are correct, the $``$100 early–type cluster galaxies at redshifts $`z1.6`$ would provide a substantial sample for cluster galaxy evolution studies.

As has been described by many authors (e.g. Longair & Seldner 1979, Prestage and Peacock 1988), if the form of the galaxy luminosity function at the redshift of the AGN is assumed, it is possible to convert the amplitude of the angular cross–correlation function for galaxies around the AGN into a spatial cross–correlation amplitude, $`B_{\mathrm{ag}}`$, derived from the spatial cross–correlation function $`\xi (r)=B_{\mathrm{ag}}\left(\frac{r}{\mathrm{Mpc}}\right)^\gamma `$. Following the method of B00, the value of $`A_{\mathrm{ag}}=0.0061\pm 0.0022`$ for the ‘all galaxies’ sample corresponds to a spatial cross–correlation amplitude of $`B_{\mathrm{ag}}=660\pm 240`$. This amplitude is comparable to the value of $`B_{\mathrm{ag}}=510\pm 120`$ derived by B00 around 3CR radio galaxies with redshifts $`z1`$. These values can be interpreted physically by comparing with the equivalent values ($`B_{\mathrm{cg}}`$) for Abell clusters calculated between the central galaxy and the surrounding galaxies; these have been derived by several authors and average to $`B_{\mathrm{cg}}350`$ for Abell class 0 and $`B_{\mathrm{cg}}710`$ for Abell class 1 . The environments surrounding these $`z1.6`$ AGN seem to be, on average, comparable to those of local clusters of richness between Abell classes 0 and 1.

## 5 Environmental variations compared with radio source properties

The results have clearly shown that the environments of these high redshift AGN are heterogeneous, with considerable variation in both the amplitude and the physical extent of any galaxy excess. In this section the variations in the galaxy overdensities are compared with the other properties of the radio source, to see if there are any indications of why the environments differ so much. The properties of the radio sources studied are provided in Table 6.

The first thing to note is that all of these radio sources have essentially the same redshifts and radio powers, and so the variations are not generally correlated with either of these properties. Equally, all 6 studied sources are steep–spectrum radio sources according to the standard definition of $`\alpha >0.5`$ (where $`\alpha `$ is defined as $`S_\nu \nu ^\alpha `$), although 1524$``$136 has a spectral index close to this cut–off; the weaker small–scale excess of this source could conceivable be related to this in some way. No fundamental difference is seen between the quasar and radio galaxy populations, and nor is any correlation with the size of the radio source seen for any of the extended radio sources.

The only issue which this comparison does bring to light is that the one source for which no overdensity is seen on either large or small scales, 2128$``$208, is an unresolved radio source; this could potentially explain this result. Extended radio sources must have been fueled continuously for timescales of at least $`10^6`$ years, and possibly up to $`10^8`$ years, which requires a plentiful and steady supply of fueling gas, as well as a relatively dense surrounding medium to confine the expanding radio lobes. Perhaps at high redshifts these conditions are only satisfied for galaxies within larger–scale structures. On the other hand, compact radio sources are relatively young, $`\mathrm{}<10^5`$ years old (e.g. Owsianik and Conway 1998), and it has been suggested that some proportion of these may fizzle out before evolving into larger radio sources (e.g. Alexander 2000). The reduced gas supply required for these short–lived compact sources may then be available to galaxies in any environment, and could explain the lack of any overdensities around 2128$``$208.

This situation has been modelled by Kauffmann & Haehnelt  who have investigated the quasar–galaxy cross–correlation function around quasars at different redshifts within hierarchical clustering models of galaxy formation. In their model they assume that quasars are triggered by galaxy mergers, with some percentage of the gas brought in by the merger then fuelling the black hole. They show that the cross–correlation function around the quasars is significantly higher for quasars with long lifetimes than for short–lived sources, in agreement with the discussion above. This difference is most pronounced for the most powerful AGN, such as those studied here.

An alternative may be that this is a lower power radio source with its flux density Doppler boosted by beaming; in this case a correlation between radio power and environmental density would explain the result. However, a beaming scenario is considered unlikely due to the steep spectral index of the source.

## 6 Discussion and Conclusions

The results of this work can be summarised as follows.

* On average, a significant overdensity of K–band selected galaxies with $`17.5<K<20.5`$ is found across the AGN fields.
* These excess counts are predominantly associated with galaxies with magnitudes $`K\mathrm{}>17.5`$ and colours $`RK\mathrm{}>4`$, which are (although not uniquely) exactly those expected for old passive cluster galaxies at redshifts $`z1.6`$.
* These galaxies are found on two different spatial scales around the AGN: pronounced concentrations of galaxies at radial distances $`\mathrm{}<150`$ kpc, and weaker large–scale overdensities extending out to between 1 and 1.5 Mpc.
* The presence or absence of galaxy excesses on these two scales differs greatly from AGN field to AGN field. All of the extended radio sources show overdensities on some scale: two show overdensities on both scales, two predominantly on large scales and one only on small scales. The one field which shows little evidence for a galaxy excess on any scale is associated with an unresolved radio source.
* Where overdensities are present on $`\mathrm{}<150`$ kpc scales, these are composed almost entirely of red galaxies. This suggests that the morphology–density relation is imprinted into the centres of clusters at a very early stage of cluster formation.
* The amplitude of the angular cross–correlation function around the AGN is a strong function of galaxy colour. The selection of only galaxies with magnitudes and colours consistent with being old passive elliptical galaxies at the AGN redshifts maximises this amplitude, indicating that old ellipticals do indeed exist in these environments.
* The spatial cross–correlation amplitude for all galaxies around the AGN is similar to that around the 3CR radio galaxies at redshifts $`z1`$; both are consistent with having the same richness as local Abell clusters of richness class 0 to 1.

It is interesting to consider what these results imply for the nature of the AGN at these redshifts, and also for cluster evolution.

Powerful radio galaxies are known to be hosted by giant elliptical galaxies at all redshifts $`z\mathrm{}<1.3`$ (e.g. Best et al. 1998, McLure & Dunlop 2000), and even out to redshifts $`z<2.4`$ a significant proportion still have radial profiles well–matched by de Vaucouleurs’ law . Little difference is found between the host galaxies of radio galaxies and radio–loud quasars at these redshifts , which suggests that the hosts of all of the $`z1.5`$ AGN studied here are likely to be giant elliptical galaxies.

The variety of different environments found for these AGN may simply reflect this fact: that the primary requirement for forming a powerful radio–loud AGN is a giant elliptical galaxy host. Such galaxies can be found in a range of environments from galaxy clusters, through groups of galaxies, to relatively isolated environments, but are generally only formed by redshifts $`z1.6`$ when they are found at the highest peaks of the primordial density perturbations, within larger–scale structures. This explains the high average richness, but still significant variation, of the large–scale galactic environments around these AGN.

To understand the pronounced smaller scale galaxy overdensity, the key question is whether these central galaxies represent a stable virialised structure, or whether they are still infalling and will merge with the central radio galaxy. If these are already virialised then the small scale galaxy overdensity is clearly related to the larger–scale structure, and indicates that at least some of the AGN lie at a special location in their environments: at the centres of their clusters or groups. This again is not too surprising: the kinetic energy of the relativistic radio jets of these AGN corresponds to the Eddington luminosity of a black hole with $`M10^9M_{}`$ , and so even assuming that these sources are fueled close to the Eddington limit they must be powered by very massive central engines and the radio lobes must be confined by dense gas. It is not unreasonable to expect these most powerful radio sources to be found within the most massive elliptical galaxies, at the centre of the cluster or group of galaxies.

In this interpretation, the four different scenarios of the presence or absence of galaxy overdensities on small and large scales correspond to four different possible locations of the AGN host galaxy: it may be a central cluster galaxy in cluster with a virialised core (excesses on both scales), an off–centre giant elliptical in a cluster, or the central galaxy in a cluster at an early stage of formation (only large–scale overdensity), the central galaxy of a compact group (only small–scale excess), or an isolated giant elliptical (no excess on any scale). With respect to this model, it is interesting that all of the extended radio sources lie in group or cluster environments.

Alternatively, if the galaxies producing the small–scale excess are not yet virialised, then they may instead be related to a hierarchical build up of the central AGN host galaxy. Pentericci et al.  found that the optical characteristic sizes of powerful radio galaxies are typically a factor of two smaller at $`z2`$ than those at $`z1`$. Since the galaxy luminosity and the characteristic size of ellipticals are related by $`L_{\mathrm{int}}r_\mathrm{e}^{0.7}`$ , this means that the radio galaxies must grow in mass by about a factor of two over this redshift interval, a large part of which may be due to mergers with companion galaxies. If the 150 kpc scale galaxy excess has a velocity dispersion $`\sigma 500`$ km s-1, then the galaxy crossing time for this structure is $`0.3`$ Gyr. A galaxy merger timescale of a few crossing times is therefore comparable to the interval of cosmic time between redshifts $`z2`$ and $`z1`$, and some of the small–scale excesses of galaxies may be destined to merge with the AGN host galaxy soon.

Whichever of these interpretations for the small–scale galaxy excess around the AGN is correct (and undoubtedly a combination of both effects is occurring), the Mpc–scale galaxy overdensities appear be associated with (proto?) cluster environments. The wealth of galaxies with magnitudes and colours consistent with being old elliptical galaxies at $`z1.6`$ indicate that the epoch of elliptical galaxy formation must be earlier than this redshift; the colours of these objects further suggest that the bulk of their stellar populations must have formed at redshifts $`z\mathrm{}>5`$ (this will be discussed further in Best et al. in preparation). By using these radio–loud AGN as signposts, it is therefore possible to define significant samples of $`z1.6`$ cluster early–type galaxies with which to extend studies of galaxy evolution back to earlier cosmic epochs.

Finally, it is interesting to compare the comoving number density of these powerful radio sources with that of nearby rich clusters. According to the pure luminosity evolution model for steep–spectrum radio sources of Dunlop and Peacock , at redshifts $`z=1`$ to 2, there are about $`10^7`$ radio sources per comoving Mpc3 within the upper decade of radio luminosities containing the sources studied here. By comparison, the comoving number density of nearby clusters of Abell richness class 0 or greater is $`4\times 10^6`$ Mpc-3, and for Abell class 1 and richer it is $`1.7\times 10^6`$ Mpc-3 (e.g. Bahcall & Soniera 1983; converted to the cosmology adopted in this paper; cf. West 1994). However, radio sources only live for between $`10^7`$ and $`10^8`$ years (e.g. Kaiser et al. 1997), and so the observed number density of powerful radio sources underestimates the number density of environments which will contain a powerful radio source at some point between redshifts 2 and 1 by a factor $`50`$ ($`=`$ cosmic time between redshift 2 and 1 divided by the lifetime of a typical radio source). Therefore, the comoving number density of environments which will contain a powerful radio source between redshifts 2 and 1 is comparable to the comoving number density of nearby Abell clusters of richness class 0 or greater. This may simply be an interesting coincidence, or alternatively may imply that the onset of a powerful radio source is a key stage in the formation of a rich cluster.

## Acknowledgements

PNB would like to thank the Royal Society for generous financial support through its University Research Fellowship scheme. The authors are very grateful to Dave Thompson and Andrea Cimatti for providing data on the R-K distributions of the CADIS and K20 surveys in electronic form. We thank the referee, Malcolm Bremer, for a careful reading of the original manuscript and helpful comments. This work is based upon observations made at the European Southern Observatory, La Silla, Chile, proposals ESO:65.O-0590(A,B) and ESO:67.A-0510(A,B). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of the USNOFS Image and Catalogue Archive operated by the United States Naval Observatory, Flagstaff Station. This work was supported in part by the European Community Research and Training Network ‘The Physics of the Intergalactic Medium’.


text:
# 1 Introduction

## 1 Introduction

In the recent past, some interest has been shown in applying methods of statistical linguistics and information theory for the analysis of DNA sequences, in particular in investigating whether the frequency distribution of nucleotides or sequences of nucleotides follows Zipf’s law , and using the Shannon entropy to identify the redundancy or the bias of a nucleotides sequence. Let us recall that, at the end of the forties, Zipf remarked that, in natural languages and in many other domains, the distribution function follows an inverse power law, which can be described, denoting by rank $`n=1`$ the most used word, by $`n=2`$ the next one and so on, and with $`a>0`$, by

$$f_n=\frac{f_1}{n^a}.$$
(1)
In 1992, it was shown  that the distribution of the nucleotides in DNA follows a Zipf law of the type given by eq. (1). However, the opinions on this statement are divided, see . Soon after , it was shown that noncoding sequences of DNA are more similar to natural languages than the coding ones, and the Shannon entropy has been used to quantify the redundancy of word. This work also raised a debate in the literature . Also, Zipf’s law seems well adapted to represent the abundance of expressed genes, with an exponent $`a1`$, as can be seen in the recent work of ref. . The origin of Zipf’s law is thought to arise from stochastic processes , especially when they can be modelled as random walks in log scale ; for further analysis, see . However, an analysis of the rank distribution for codons, performed in many genes for several biological species, led the authors of  to fit experimental data with an exponential function. In particular, by considering separately different coding DNA sequences, they studied the relation between the parameter in the exponential, the frequency of the rank one and the length of the sequence for different genes. From this very short overview, it follows that the determination of the kind of law the codon rank distribution follows is extremely interesting, in investigations of the nature of the evolutionary process, which has acted upon the codon distribution, i.e. the eventual presence of a bias.

In the last few years, the number of available data for coding sequences has considerably increased, but apparently no analysis using the whole set of data has been performed. The aim of this paper is to perform such a study. As the result of our investigation, we will point out that the rank of codon usage probabilities follows a universal law, the frequency function of the rank-$`n`$ codon showing up as a sum of an exponential part and a linear part. A quantum-mechanics-inspired model is proposed, which provides for the form of the observed law and gives the correct sign in the different terms. Such a universal behaviour suggests the presence of general biases, one of which is identified with the total *exonic* $`GC`$ content. Indeed, the values of the parameters appearing in the fitting expression are plotted versus the total percentage of exonic $`GC`$ content of the biological species and are reasonably well fitted by a parabola. Finally, from the obtained expression, we derive the theoretical prediction that the usage probability for *rank-ordered* codons is the same in any genic region having the same exonic $`GC`$ content for any biological species.

We compute the Shannon entropy  for amino-acids and find that its behaviour in function of the exonic $`GC`$ content is also a parabola, whose apex is around the value 0.50 of the $`GC`$ content.

## 2 Codon usage probabilities distribution

Let us define the usage probability for the codon $`XZN`$ ($`X,Z,N\{A,C,G,U\}`$) as

$$P(XZN)=\underset{n_{tot}\mathrm{}}{lim}\frac{n_{XZN}}{N_{tot}},$$
(2)
where $`n_{XZN}`$ is the number of times the codon $`XZN`$ has been used in the analysed biosynthesis process for a given biological species and $`N_{tot}`$ is the total number of codons used in all considered processes. It follows that our analysis and predictions hold for biological species with sufficiently large statistics of codons. For each biological species, codons are ordered following the decreasing order of the values of their usage probabilities, i.e. codon number 1 corresponds to the highest value, codon number 2 is the next highest, and so on. We denote by $`f(n)`$ the probability $`P(XZN)`$ of finding $`XZN`$ is in the $`n`$-th position. Of course the same codon occupies in general two different positions in the rank distribution function for two different species. We plot $`f(n)`$ versus the rank and we determine that the best fit to the data can be reached by the sum of an exponential function, a linear function in the rank and a constant, i.e.

$$f(n)=\alpha e^{\eta n}\beta n+\gamma ,$$
(3)
where $`0.0187\alpha 0.0570`$, $`0.050\eta 0.136`$, $`\mathrm{0.82\hspace{0.33em}10}^4\beta \mathrm{3.63\hspace{0.33em}10}^4`$, $`\gamma =0,016`$ are constant depending on the biological species. These four constants have to satisfy the normalization condition

$$\underset{n}{}f(n)=1.$$
(4)
In table 1 we list the 40 biological species (6 vertebrates, 4 plants, 3 invertebrates, 2 fungi and 25 bacteriae) – with a number of codons ranging between 800 000 and 20 000 000 in decreasing order (data from GenBank release 129.0 ), whose codon usage has been fitted, specifying for each biological species the value of the parameters, computed by a best-fit procedure and the corresponding $`\chi ^2`$. Here and in the following, the $`\chi ^2`$ coefficient is defined by

$$\chi ^2=\underset{i}{}\frac{(y_iy(x_i))^2}{y(x_i)},$$
(5)
where $`x_i`$ are the experimental abscissae, $`y_i`$ the experimental values and $`y(x_i)`$ the fitted ones. In some cases, $`y(x_i)`$ takes vanishing or negative values for a few points, hence the $`\chi ^2`$ is not reported. In figs. 14, we report the plot of $`f(n)`$ as a function of $`n`$ for a few biological species (Homo sapiens, Drosophila melanogaster, Arabidopsis thaliana, and Escherichia coli). The plot has been cut to $`n=61`$ to take into account the fact that in standard code there are three Stop codons (to end the biosynthesis process), whose function is very peculiar. For the same reason, the $`\chi ^2`$ has been computed by taking into account the 61 coding codons only. In table 6, we report the type of the twenty-first most used codons of the observed rank distribution $`f(n)`$. A similar study, for a sample of twenty vertebrates with a codon statistics larger than 100 000, reveals that, for almost all biological species, the four most used codons are $`GAG`$, $`CUG`$, $`AAG`$ and $`CAG`$. All these codons have a $`G`$ nucleotide in third position and three of them encode doublets. An analysis performed on the chloroplast codon usage for a sample of five plants gives the same result for the rank distribution $`f(n)`$, see table 2 and fig. 5 (Chloroplast Arabidopsis thaliana). We also report, in table 3, the values of the parameters and the $`\chi ^2`$ for a sample of nine mitochondriae with a codon statistics larger than 15 000. The fits for Homo sapiens and Arabidopsis thaliana are presented in figs. 6 and 7. We point out, however, that for mitochondriae the codon usage frequencies distribution for several species (e.g. Arabidopsis thaliana or Drosophila melanogaster) is ill fitted by eq. (3). It may be an indication that mitochondriae do not follow the universal law (3). Note that the mitochondrial codes have a few differences with the eukaryotic code and vary slightly between species, see e.g. . In these cases, the $`\chi ^2`$ has been computed over the corresponding coding codons. The value of the constant $`\gamma `$ is approximately equal to $`1/61=0.0164`$ or to $`1/64=0.0156`$, i.e. the value of the codon usage probability in the case of uniform and not biased codon distribution. Therefore the other two terms in eq. (3) can be viewed as the effect of the bias mechanism. The appearance of the linear term is more intriguing. Let us remark that in , where an exponential function is used to fit the rank of usage in genes (not the rank of usage probability), the linear term has not been observed, as its contribution becomes noticeable for approximately $`n20`$. Owing to the analysis of genes (with at most a few hundred codons), the fits in that paper end before this value of the rank. It is believed that the main causes of codon usage bias are the translational efficiency, the selection pressure and the spontaneous mutations. From the smallness of parameter $`\beta `$ in (3), it is tempting to identify the latter as the consequence of mutation effect and the first term in (3) as the effect of selection pressure, i.e. the interaction with the environment.

Let us try to build up a simple model to explain the universal structure of eq. (3). In a physical quantum system, with discrete energy levels, the occupation number of the $`i`$-th energetic level is proportional to $`\mathrm{exp}(E_i/kT)`$, where $`E_i`$ is the energy of the level, $`T`$ is identified with the temperature and $`k`$ is the Boltzmann constant. Inspired by this analogy, let us assume that the rank of the probability distribution is a function of some kind of energy. In this case the codon usage probability is, for any codon, equal to $`f(n)=\mathrm{exp}(aE_n)`$, where $`a`$ is a positive constant. We assume that the system satisfies periodic conditions in the $`y`$ and $`z`$ directions, with a very small spacing in the $`z`$ direction, so that the excited modes along this direction can be neglected, and that the interaction with the environment can be modelled as the switching on of a constant uniform magnetic field in the direction of the $`z`$ axis. As a consequence, the Landau levels appear, the energetic levels now being in natural units

$$E_nE_0=\omega (n+1/2),$$
(6)
where $`\omega `$ is the so-called cyclotron frequency. The wave function is

$$\psi (x,y,z)=\mathrm{\Phi }_n(xlx_0)e^{ik_yy}e^{ik_zz},$$
(7)
where

$$k_y=\frac{2\pi l}{L_y},k_z=\frac{2\pi }{L_z},l𝐙,l0,$$
(8)
and $`\mathrm{\Phi }_n(xlx_0)`$ is the wave function of the $`n`$-th energy level of a one-dimensional harmonic oscillator shifted by $`lx_0`$ ($`x_0=2\pi /qBL_y`$). There is a degeneracy due to the arbitrariness of $`k_y`$. A quantun charged physical system in an excited state spontaneously falls down to a lower state, the equilibrium configuration being the system in the fundamental state. In the case of codons, the equilibrium condition is the configuration where they occupy different “energy” levels, the situation being the analogue of a quantum system in the presence of an external pumping-back source. In order to model this situation in a simple way, we assume that the equilibrium probability distribution is a function $`f(n)f(n)_{B=0}=F(E_n,\gamma _n)`$, where $`\gamma _n`$ is the “spontaneous emission coefficient” of the $`n`$-th level, which is taken to be very small because of the assumed equilibrium condition. We can make a first-order development in $`\gamma _n`$:

$$f(n)f(n)_{B=0}=F(E_n,\gamma _n=0)+\frac{dF}{d\gamma _n},$$
(9)
where $`dF/d\gamma _n`$ is negative, the occupation probability decreasing with the increasing emission. We compute $`\gamma _n`$ in the standard dipole approximation  between the $`n`$-th and the $`k`$-th levels ($`k<n`$). In the case of our model, we get

$$\gamma _nn\delta _{k,n1}.$$
(10)
So, replacing $`F(E_n,\gamma _n=0)`$ by the Maxwell–Boltzmann expression, we obtain eq. (3). The values of the parameters depend on the biological species, as the selection pressure (in the model, the value of the charge and of the magnetic field) is different for different biological species. Even if our simple model is able to explain the universal behaviour and the right sign of the linear term, one should rather consider it as a toy model. However, the fact that one can build a simple mathematical scenario, able to reproduce the observed distribution, provides us with an indication of the existence of some strong physico-chemical constraints that add to the random effects.

Since it is well known that the $`GC`$ content plays a strong role in the evolutionary process, we expect the parameters to depend on the total $`GC`$ content of the genes region (here the total exonic $`GC`$ content), which is indeed correlated with the evolution of the system (see  and references therein). We have investigated this dependence and report, in table 4, the fits of $`\alpha `$ and $`\beta `$ to the total exonic $`GC`$ content $`Y_{GC}`$ of the biological species. One finds that the values of $`\alpha `$ and $`\beta `$ are well fitted by polynomial functions (with $`0Y_{GC}100`$ in percentage):

$`\alpha `$ $`=`$ $`0.211450.00776Y_{GC}+\mathrm{7.92\hspace{0.33em}10}^5Y_{GC}^2,\chi ^2=0.0262,`$ (11)
$`10^2\beta `$ $`=`$ $`0.100960.00345Y_{GC}+\mathrm{3.50\hspace{0.33em}10}^5Y_{GC}^2,\chi ^2=0.0170.`$ (12)
The two parameters $`\alpha `$ and $`\beta `$ appear to be correlated. Indeed the plot representing $`\beta `$ as a function of $`\alpha `$ is satisfactorily fitted by a regression line (see fig. 5):

$`10^2\beta `$ $`=`$ $`0.00851+0.375\alpha ,\chi ^2=0.0218.`$ (13)
The value of the $`\eta `$ parameter is largely uncorrelated with the total exonic $`GC`$ content. Let us recall that, however, $`\eta `$ is a function of $`\alpha `$ and $`\beta `$ due to the normalization condition of eq. (4). Indeed we have (assuming $`e^{65\eta }0`$) 111Note that the result is almost unchanged if the data are normalized on the 61 coding codons.

$$1=\frac{\alpha e^\eta }{1e^\eta }+2080\beta +64\gamma .$$
(14)
Using the fits for $`\alpha `$ and $`\beta `$, we can write the probability distribution function for any biological species, whose total $`GC`$ content in per cent in the exonic regions is $`Y_{GC}`$, as

$$f(n)=(\alpha _0+\alpha _1Y_{GC}+\alpha _2Y_{GC}^2)e^{\eta n}n(\beta _0+\beta _1Y_{GC}+\beta _2Y_{GC}^2)+\gamma ,$$
(15)
where $`\eta `$ is obtained by solving eq. (14). Of course we are not able to predict which codon occupies the $`n`$-th rank. Finally, let us remark that the total exonic $`GC`$ content $`Y_{GC}`$ has to satisfy the consistency condition

$$Y_{GC}=\frac{1}{3}\underset{iI}{}d_if(i),$$
(16)
where the sum is over the set $`I`$ of integers to which the 56 codons containing $`G`$ and/or $`C`$ nucleotides belong and $`d_i`$ is the multiplicity of these nucleotides inside the $`i`$-th codon.

## 3 Amino-acids rank distribution

It is natural to wonder if some kind of universality is also present in the rank distribution of amino-acids. From the available data for codon usage, we can immediately compute (using the eukaryotic code) the frequency of appearance of any amino-acid $`F(n)`$ ($`1n20`$) in the whole set of coding sequences. The calculated values as a function of the rank are satisfactorily fitted by a straight line,

$$F(n)=F_0Bn.$$
(17)
The parameters $`F_0`$, $`B`$ and the corresponding $`\chi ^2`$ for the fits are reported in table 7. A better fit can be obtained in general by using a third-degree polynomial; however, the range of the four parameters for this fit is larger than the range of the 2-parameter fit. For a few biological species, we give below the parameters for the two fits, see also the figures of table 9. The plots of the linear fits for few biological species are given in table 8. Note that the 21st point is just the contribution of the Stop codons, which of course has not been taken into account for the fits. One can remark that the most frequent amino-acid is always above the line. This can be easily understood in the light of eq. (3). Indeed, the most frequent amino-acids get, in general, a contribution of the exponential term of (3) with a low value of $`n`$.

Of course, the frequency of an amino-acid is given by the sum of the frequencies of its encoding codons given by (3). If the ranks of the encoding codons were completely random, we do not expect that their sum should take equally spaced values, as is the case in a regression line. Therefore, we can infer, for the biological species whose amino-acids frequency is very well fitted by a line, the existence of some functional constraints on the codon usage.

We report in table 10 the distribution of the amino-acids for the different biological species.

However, the behaviour predicted by eq. (3) fits the experimental data very well, while the shape of the distribution of amino-acids seems more sensible to the biological species. In fact, one can remark on many plots of the amino-acid distributions (see e.g. table 8), the existence of one or two plateaux, for which we do not have any explanation yet.

## 4 Consequences of probability distribution

We now derive a few consequences of eq. (3). In the following, we denote by $`y`$ the *local* exonic $`GC`$ content (i.e. for coding sequences of genes) for a given biological species. Let us assume that the exonic $`GC`$ content of a biological species is essentially comprised in the interval $`y_1y_0=\mathrm{\Delta }`$ (e.g. for Homo sapiens $`y_0=35\%`$ and $`y_1=70\%`$). We can write

$$f(n)=\frac{1}{\mathrm{\Delta }}_{y_0}^{y_1}f(y,n)𝑑y.$$
(18)
Since the l.h.s. of the above equation has the form given by eq. (3) for any $`n`$ and for any biological species, if we do not want to invoke some “fine tuning” in the integrand function $`f(y,n)`$, we have to assume that

$$f(y,n)=a(y)e^{\eta n}b(y)n+\gamma $$
(19)
with the condition

$$\alpha =\frac{1}{\mathrm{\Delta }}_{y_0}^{y_1}a(y)𝑑y,\beta =\frac{1}{\mathrm{\Delta }}_{y_0}^{y_1}b(y)𝑑y.$$
(20)
As a consequence, we predict that the codon usage probability is the same for any codon in any exonic genic region with the same $`GC`$ content. The form of the $`a(y)`$ and $`b(y)`$ functions is yet undetermined.
For Homo sapiens, we remark that the total exonic $`GC`$ content $`Y_{GC}`$ is, in a very good approximation, equal to the mean value of the interval $`[y_0,y_1]`$. Therefore, inserting (18) and (19) into (16), we derive that the functions $`a(y)`$ and $`b(y)`$ have to be *linear* functions of $`y`$. This theoretical derivation is in accordance with the conclusions of Zeeberg  obtained by an analysis for 7357 genes. On a quantitative level, using the numerical linear fits of Zeeberg, we find a very good agreement with our calculations. Note that this result is not in contradiction with eq. (15), since the previous analysis is valid for the fixed value of the exonic $`GC`$ content for Homo sapiens. For bacteriae, the range of variation $`\mathrm{\Delta }`$ of the local exonic $`GC`$ content is very small. Therefore we expect the functions $`a(y)`$ and $`b(y)`$ to have the same shape as the functions $`\alpha `$ and $`\beta `$ given in eqs. (11) and (12). Hence the functions $`\alpha `$ and $`\beta `$ depend on the biological species.

We compute the Shannon entropy, given by

$$S=\underset{n}{}f(n)\mathrm{log}_2f(n)$$
(21)
for the codons of a biological species and plot it versus the total exonic $`GC`$ content; see figure in table 5. The Shannon entropy is rather well fitted by a parabola:

$$S=2.2186+0.144Y_{GC}0.00146Y_{GC}^2,\chi ^2=0.0315.$$
(22)
Note that the parabola has its apex for $`y0.50`$, which is expected for the behaviour of the Shannon entropy for two variables (here $`GC`$ and its complementary $`AU`$).

The same behaviour have been put in evidence by analogous computations made by Zeeberg  for Homo sapiens. So it seems that the entropy in the gene coding sequences and the total exonic region entropy in function of the exonic $`GC`$ content show the same pattern.

In conclusion, the distribution of the experimental codon probabilities for a large total exonic region of several biological species has been very well fitted by the law of eq. (3). The spectrum of the distribution is universal, but the codon, which occupies a fixed level, depends on the biological species. Indeed, a more detailed analysis shows that, for close biological species, e.g. vertebrates, a fixed codon occupies almost the same position in $`f(n)`$, while for distant biological species the codons occupy very different positions in the rank distribution. We have also derived that the codon frequency for any genic region is, for fixed biological species and fixed $`GC`$ content, the same. Entropy analysis has shown that the behaviour observed in the genes with different $`GC`$ content for the same biological species is very similar to the one shown by the total exonic region with different $`GC`$ content for the different biological species.

Acknowledgements: The authors would like to thank A. Arneodo for fruitful discussions. A.S. is indebted to the Université de Savoie for financial support and LAPTH for its kind hospitality.


text:
# Study of the formation and decay of electron-hole plasma clusters in a direct-gap semiconductor CuCl

thanks: Author to whom correspondence should be addressed

## Abstract

The master equation for the cluster-size distribution function is solved numerically to investigate the electron-hole droplet formation claimed to be discovered in the direct-gap CuCl excited by picosecond laser pulses \[Nagai et al., Phys. Rev. Lett. 86, 5795 (2001); J. Lumin. 100, 233 (2002)\]. Our result shows that for the excitation in the experiment, the average number of pairs per cluster (ANPC) is only around 5.2, much smaller than that (106 typically for Ge) of the well studied electron-hole droplet in indirect-gap semiconductors such as Ge and Si.

In photoexcited semiconductors, electrons and holes are generated by the photon absorption and interact with each other via the Coulomb interaction. When the exchange- and correlation-induced attractions among the high density photogenerated carriers are strong enough to overcome the fermionic repulsive pressure and to reduce the chemical potential below the exciton energy, electron-hole droplets (EHDs) appear as first order transition.lvk  Investigation of the formation kinetics of EHD provides an opportunity to understand the quantum many-body effects. In semiconductor with an indirect gap and correspondingly long optical lifetimes, such as Ge and Si, EHDs have been investigated both theoretically and experimentally in great detail.sige  Within the coexistence regime, relatively large EHD clusters are formed and the kinetics of the EHD formation is well understood in terms of the master equationsil ; haug  where the ANPC $`n`$ is very large (for Ge $`n`$ is about 106).

Unlike indirect gap semiconductors, in direct gap semiconductors the situation is very complicated because the system cannot reach quasiequilibrium during the carrier lifetime (typically in the order of 1 ns). Despite extensive experimental efforts, no decisive evidence of the EHD formation has yet been reported for a very long time.haug2  EHD formation has been predicted theoretically in polar CdS and CdSe.rice1  Also, EHD-like behavior has been observed in them experimentally in which time-resolved emission measurements are performed.vbtim1 ; ssion1  This method has also been applied to GaAs. It reveals that electron-hole plasma created by band-to-band excitation has extremely high temperature due to the overheating caused by the reduction of band gap. Very recently, Nagai et al. used CuCl which has large binding energy and band gap (binding energy of 1s-Z3 exciton Eex=213 meV and band-gap energy Eg=3.396 eV) to overcome the problem. They therefore are able to study the formation of EHD in CuCl via time-resolved emission and pump-probe reflection measurements by femtosecond excitation in the mid-infrared region under resonant exciton excitation. They reported that clear EHD-like behavior in CuCl has been found.ku1 ; ku2

However, due to the fact of the short lifetime of the carrier, one expects only small clusters can be formed and sustained. Two decades ago, Haug and Abraham performed a theoretical investigation of the EHD formation in GaAs.haug1  They reported for the typical picosecond excitation the ANPC is at most around 2 where there should be no EHD formation can be found. In this paper we use the master equation theory developed by Haug and Abrahamhaug1  to investigate the femtosecond excitation in CuCl in the experiments of Nagai et al..ku1 ; ku2  This theory has been proven to be very sucessful in dealing with the EHD formation in indirect gap semiconductors.sil ; haug  In this study, we focus on how big the ANPC can be in this material and how this number is affected by different parameters.

The concentration $`f_n(t)`$ to find at a given time t a cluster with $`n`$ e-h pairs per unit volume is given by the following master equation equations:haug1  For $`n>1`$,

$`{\displaystyle \frac{f_n}{t}}=j_{n1}j_n`$ (1)
$`j_n=g_nf_nl_{n+1}f_{n+1}`$ (2)
where $`j_n`$ is the net probability current between the clusters with n and n+1 e-h pairs, $`g_n`$ is the gain rate of a cluster with n e-h pairs and is approximated by

$`g_n`$ $`=`$ $`4\pi R_n^2v_0f_1`$ (3)
$`v_0`$ $`=`$ $`\sqrt{kT/2\pi m^{}}`$ (4)
$`n`$ $`=`$ $`{\displaystyle \frac{4\pi }{3}}R_n^3\rho _0`$ (5)
$`R_n`$ is the radius of the cluster, $`v_0`$ is the thermal exciton velocity, $`\rho _0`$ is the liquid density. The loss rate $`l_n`$ is the sum of the evaporation rate $`e_n`$ and the recombination rate $`n/\tau _n`$:

$$l_n=e_n+n/\tau _n.$$
(6)
$`e_n`$ is given by the Richardson-Dushman currentsil

$$e_n=\gamma \left(\frac{m^{}kT}{2\pi \mathrm{}^2}\right)^{3/2}\frac{g_{n1}}{f_1}\mathrm{exp}[(\varphi +c\sigma n^{1/3})\beta ]$$
(7)
$`\gamma `$ is the degeneracy of the exciton level, $`\varphi `$ is the binding energy of the liquid with respect to the exciton. $`c\sigma n^{1/3}`$ is the correction of the binding energy with $`\sigma `$ denoting the surface energy of the liquid. Typically we select $`\gamma `$ =1. $`\sigma =C\phi \rho _0r_0`$ with $`C1`$, $`\phi `$ standing for the condensate energy per electron-hole pair and $`r_0=(3/4\pi \rho _0)^{1/3}`$.lvk  The change of exciton (corresponding to $`n=1`$) is treated separately:

$$\frac{f_1}{t}=G(t)\underset{n=1}{\overset{\mathrm{}}{}}\frac{f_n}{\tau _n}2j_1\underset{n=2}{\overset{\mathrm{}}{}}j_n$$
(8)
where G(t) is the generation rate of excitons.

We solve above equations numerically for the example of CuCl. The parameters for CuCl at 4.2 K are $`m^{}`$=2.0$`m_0`$ and $`\varphi =5.12\times 10^{14}`$ erg.d1  The value of the binding energy $`\varphi `$ is taken to be the biexciton binding energy which is the lower limit of the binding energy. Our computation indicates that the result is almost unchanged if one increases $`\varphi `$ by one order of magnitude. We obtain $`\rho _0`$ from $`r_s`$=$`\sqrt[3]{3/4\pi na_0^3}`$, where $`a_0`$ is the Bohr radius and $`r_s`$=1.8. ku1  $`\sigma =0.54`$ erg$``$cm-2. In our calculation we assume $`\tau _n\tau `$. According to the experiment, we choose a Gaussian pulse:

$$G(t)=G_0e^{(tt_0)^2/t_r^2}$$
(9)
where the pulse width $`t_r=0.2`$ ps, and the peak generation rate is $`G_0=6.1\times 10^{33}`$ cm-3s-1, corresponding to a excitation density of 4.3 mJ$``$cm-2. The main results of our calculation are plotted in Figs. 1-3.

Figure 1 shows the development of the distribution of cluster distribution. For the short pulse, the exciton concentration $`f_1(t)`$ rapidly reaches its maximum near the pulse center, $`t_0=2`$ ps. Then it decreases quickly due to the formation of the clusters. It is also noted from the figure that within 1 ps after the center of the pulse, clusters of all sizes reaches to their maximum values. As any cluster with size $`n2`$ comes from the exciton, $`f_1(t)`$ sharply drops 5 orders of magnitude from its peak within that period of time. After that, because of the recombination and the evaporation, almost all clusters start to decay, except for some small clusters (say $`n=1`$ and 2) which increase a little bit due to the evaporation from bigger clusters. The lifetime of large cluster (say $`n=50`$) is much shorter than that of the small one as the later gets compensation from the decay of the larger clusters and also the recombination rate $`n/\tau _n`$ is propotional to $`n`$. After 10 ns almost all cluster concentrations are below $`10^{10}`$ cm-3. This can be seen more quantitatively in Fig. 2 where the concentration of some selected clusters are plotted as functions of time.

In order to see in average how large the cluster can be under the experiments by Nagai et al., we plot in Fig. 3 the ANPC $`n`$ as function of time for three different lifetimes $`\tau `$. $`n`$ can be calculated from

$$n=\frac{\underset{n=1}{\overset{\mathrm{}}{}}nf_n}{_{n=1}^{\mathrm{}}f_n}.$$
(10)
One notices from the figure that for the excitation in the experiment, ANPC is only 5.2, orders of magnitude smaller than the number in the EHD of indirect semiconductor such as Ge and Si. This result is consistent with the regular expectation of the excitation in the direct gap semiconductors. Moreover, if one assumes that the larger clusters are stable and therefore $`\tau _n`$ ($`n2`$) is larger than $`\tau _1`$. We demonstrate in the same figure the result when $`\tau _1=1`$ ns but $`\tau _n=10`$ ns when $`n2`$. One finds that this does not affect the ANPC, although the metastable time is much longer. We further find that if we increase the excitation $`G_0`$ by one order of magnitude, we can only double the ANPC. These results indicate that the feature of low ANPC is very robust. These results indicate that what measured in CuCl by Nagai et al.ku1 ; ku2  may not come from the EHD formed from exciton gas, but instead possibly come from some bubbles of excitons in metallic liquid.ku2

MWW is supported by the “100 Person Project” of Chinese Academy of Sciences and Natural Science Foundation of China under Grant No. 10247002. He would also like to thank S. T. Chui at Bartol Research Institute, University of Delaware for hospitality. LJ would like to thank M. Q. Weng for his help during this work. MWW, MN and MK is supported in part by the Solution Oriented Research for Science and Technology (SORST) of the Japan Science and Technology Corporation (JST).