text:
# Overdensities Of Extremely Red Objects In The Fields Of High-Redshift Radio-Loud Quasars

## 1 Introduction

The age and space density of the oldest and most massive galaxies at each epoch provide the strongest constraints on models of galaxy formation. This is one of main reasons that there has been so much focus on Extremely Red Objects (EROs). It is the high-redshift systems that are most important for testing models, and the ERO selection criterion is designed to find old, passively evolving galaxies at $`z1`$. Because the 4000 Å break in evolved galaxies shifts from optical to near-infrared at $`z1`$, EROs are defined in terms of optical to near-infrared color, with two of the most used definitions being $`RK5.0`$ and $`IK4.0`$ mag.

The co-moving volume density, or even the number density, of massive, evolved galaxies at $`z1`$ can place strong constraints on models of galaxy formation. Whereas monolithic collapse models (Eggen, Lynden-Bell, & Sandage, 1962; Larson, 1974; Tinsley & Gunn, 1976; Loeb & Peebles, 2003) predict pure luminosity evolution of early-type galaxies from the period of formation until the present-day, hierarchical formation models (White & Rees, 1978; Baugh, Cole, & Frenk, 1996; Kauffmann, 1996; Baugh et al., 1998; Somerville, Primack, & Faber, 2001) predict significant density evolution from $`z1`$ to $`z0`$. Testing the models has proven to be difficult, and the results are not yet conclusive. Some surveys have found a deficit of early-type galaxies at $`z1`$ (Barger et al., 1999; Smith et al., 2002; Roche et al., 2002) whereas others have not (Daddi, Cimatti, & Renzini, 2000b; Cimatti et al., 2002; Im et al., 2002). Part of the disagreement may be caused by the strong clustering of EROs on the sky (Daddi et al., 2000a; Firth et al., 2002; Roche et al., 2002), and large fields must be surveyed in order to overcome this problem.

The ERO definition also selects dusty starburst galaxies whose spectral energy distribution can mimic those of old, passively evolving galaxies. It is still unclear how large a fraction of the ERO population consists of dusty starburst galaxies. Several methods are being used to separate the two galaxy types, including photometric (Mannucci et al., 2002; Smail et al., 2002; Miyazaki et al., 2002), spectroscopic (Cimatti et al., 2002; Saracco et al., 2003), and those based on morphology (Moriondo, Cimatti, & Daddi, 2002; Stiavelli & Treu, 2000; Yan & Thompson, 2003). The results so far seem to suggest that roughly two-thirds of the EROs at $`K19`$ mag are old, passively evolving galaxies. However, it is also clear that the situation is more complex than just a simple distinction into old, evolved galaxies and dusty starbursts since systems with old populations and residual star formation also appear in ERO samples (Yan & Thompson, 2003).

In this paper, we are concerned with EROs in the fields of high-redshift quasars. Some early studies serendipitously found EROs in the fields of $`z>2`$ radio galaxies and quasars (McCarthy, Persson, & West, 1992; Hu & Ridgway, 1994; Dey, Spinrad, & Dickinson, 1995). Although some of these EROs were later shown to lie in the foreground of the radio galaxies (Graham & Dey, 1996), probing the fields around high-redshift radio-loud AGN (active galactic nuclei) has been a promising method of finding high-redshift galaxies. Luminous quasars and radio galaxies are often associated with high-density galaxy environments (Ellingson, Yee, & Green, 1991; Hill & Lilly, 1991; Hall & Green, 1998; Wold et al., 2000; Best, 2000; Venemans et al., 2002; Best et al., 2003), and the rich environments are often interpreted as clusters or groups of galaxies hosting the AGN. By combining a targeted search with a color criterion like that of EROs, high-redshift early-type galaxies can be effectively filtered out from the numerous foreground galaxies. This becomes very important at high redshifts, since even for a rich cluster at $`z1`$, most of the galaxies observed using a single broad-band filter are foreground or background galaxies (Dickinson, 1997). A few surveys have found that EROs seem to be more common in high-redshift AGN fields than in the general field (Chapman, McCarthy, & Persson, 2000; Cimatti et al., 2000; Hall et al., 2001; Sánchez & González-Serrano, 2002), and there are individual AGN fields which show striking overdensities of EROs (Thompson, Aftreth, & Soifer, 2000; Liu et al., 2000; Haines et al., 2001).

Even though galaxy excess in the fields of AGN are often being interpreted as groups or clusters physically associated with the AGN (spectroscopically confirmed in a few cases (Dickinson, 1997; Deltorn et al., 1997)), gravitational lensing is also known to cause correlations (and anti-correlations) between high-redshift quasars and optically bright foreground galaxies (Webster et al., 1988; Hammer & LeFévre, 1990; Benítez, Martínez-González, & González-Serrano, 1995; Benítez, Martínez-González, & Martín-Mirones, 1997; Norman & Impey, 2000; Myers et al., 2003). It is possible that a similar situation could arise between high-redshift quasars ($`z2`$) and EROs, because the EROs are likely to be tracers of massive structures which may boost the fluxes of distant background quasars by gravitational lensing.

The aim with our present study is to test whether EROs are more common in the fields of $`z2`$ quasars. To our knowledge, this is the first systematic survey to probe the surface density of EROs in the fields of quasars at $`z>2`$, and preliminary results were reported by Wold et al. (2003). The quasars in our sample are selected from the catalogs of Barthel et al. (1988) and Hewitt & Burbidge (1993), and form a subset of a larger sample of 40 designed to study host galaxies (Armus et al. in preparation). The 40 sample quasars are nearly evenly divided between radio-loud and radio-quiet systems with comparable $`V`$ magnitudes and redshifts. The subsample of 13 analyzed here were selected randomly from the sample of 40 and followed up with wide-field near-infrared and optical imaging. The sample is listed in Table 1. In this table, we have also listed the radio loudness, defined in terms of the ratio of the flux densities at rest-frame 5 GHz and 2500 Å, $`S_{5GHz}/S_{2500}`$ (Sramek & Weedman, 1980). According to the definition by Kellermann et al. (1989), a quasar is radio-loud if $`S_{5GHz}/S_{2500}>10`$, and thus two of the 13 quasars are classified as radio-quiet (one of these might be radio-intermediate instead of radio-quiet) whereas the rest are radio-loud. In calculating the radio loudness, we used the $`K_s`$-band fluxes of the quasars as found from our data, and the radio fluxes at 5 GHz listed in the catalog of Barthel et al. (1988) unless otherwise is noted in the table. The $`K_s`$-band fluxes were converted to rest-frame 2500 Å by assuming a power-law spectrum, $`S_\nu \nu ^\alpha `$, where $`\alpha `$ was found from the median quasar energy distribution by Elvis et al. (1994). The conversion to rest-frame 5 GHz radio flux was done by assuming the same power-law form of the spectrum, but using spectral indices listed in the catalog by Barthel et al. (1988), unless otherwise is noted in the table.

The structure of the paper is as follows. In the next section, we describe the observations and data reduction. Section 3 deals with the construction of object catalogs and addresses their incompleteness. In Section 4 we present the results of the analysis. The surface density of EROs is calculated and the ERO excess evaluated, carefully taking into account their strong clustering. We also investigate the $`K_s`$-band number counts and the radial distribution of EROs around the quasars. In Section 5, we discuss two different interpretations of the results, and the conclusions are drawn in Section 6.

All magnitudes are given in the Vega system, and our assumed cosmology has $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`H_0=70`$ km s-1 Mpc-1.

## 2 Observations and data reduction

The data consist of wide-field $`R`$\- and $`K_s`$-band images taken with the Palomar 200-inch telescope using the prime-focus instruments COSMIC111Carnegie Observatories Spectrograph and Multiple Imaging Camera (Kells et al., 1998) and PFIRCAM222Prime Focus Infrared Camera (Jarrett et al., 1994a). The COSMIC instrument is equipped with a 2048$`\times `$2048 Tektronix CCD with a pixel scale of 0.2856 arcsec, and the PFIRCAM employs a 256$`\times `$256 HgCdTe array with a pixel scale of 0.494 arcsec. The resulting fields of view in $`R`$ and $`K_s`$ are therefore 9.7$`\times `$9.7 and 2.1$`\times `$2.1 arcmin2, respectively.

Images were obtained during several observing runs from September 1995 to July 2001. We used a 9- or 27-point dither pattern for the $`K_s`$ images, and on-chip integration times of 4 sec. The number of co-adds was set so that the total integration time at each dither point was 20–30 sec, and the offsets between individual frames were typically 10–15 arcsec. The $`R`$-band exposures were divided into integrations of 300–600 sec each, with an offset of 10–20 arcsec between each exposure. The total integration time in $`R`$ and $`K_s`$ for each field is listed in Table 1.

Flat-field images in $`K_s`$ were made by taking the median of dark-subtracted quasar images. Typically, a flat-field was made for each dither sequence, and this was found to be sufficient in most cases. After flat-fielding, the images were sky subtracted, corrected for bad pixels and registered using the dimsum package (Stanford, Eisenhardt & Dickinson, 1995) within iraf.333iraf is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.

The sky subtraction in dimsum is performed by subtracting a running sky frame made by taking the median of nine neighboring images. After the sky subtraction, masks are created for objects and cosmic rays and the masks are utilized to make a first pass combined image. The first pass image is then used to make a new and improved mask by including fainter objects. The new object mask is used in the final pass where the sky subtraction is repeated and the images aligned and combined. At this final step, we block replicated the images by a factor of two so that sub-pixel shifts were used for the aligning. The resulting images therefore have a pixel scale of 0.247 arcsec.

The $`R`$-band images were reduced in iraf in a standard manner. The bias level was subtracted either by using a bias frame or by using a 3rd order Legendre to fit the overscan region which was then subtracted from each row. A flat-field was formed for each quasar field by taking the median of the quasar images, or of the quasar images in combination with exposures of other fields taken near in time. Typically 15–30 frames were used for this, and all emission $`2\sigma `$ above the average background was masked out. This method rejected cosmic rays, stars and galaxies visible in the quasar images and yielded a high quality flat-field frame. The flat-fielded images were then aligned using 20–40 stars in each frame, and averaged using imcombine with a high threshold set to remove bad pixels that had been flagged by giving them high, discrepant values. The combined $`R`$ images were aligned with, and put to the pixel scale of, the $`K_s`$ images using the iraf tasks geomap and geotran.

The flux calibration of the $`R`$-band images was accomplished using standard stars from Landolt (1992), and the $`K_s`$-data were calibrated using standard stars from Persson et al. (1998). The conditions were photometric, and the transparency good throughout most of the runs. The final photometric status was checked by the amount of variation within the individual exposures of the same source and by following standard stars throughout the night. For fields that were imaged during non-photometric conditions, calibration images were obtained later during photometric conditions to allow for an accurate flux calibration. During the photometric nights, the zeropoint varied by less than 0.1 mag. For a subset of the images, we checked our absolute photometry against the Two Micron All Sky Survey (2MASS) point source catalog (Skrutskie et al. 1997; Cutri et al.444http://www.ipac.caltech.edu/2mass/releases/second/ doc/explsup.html). A very good agreement, typically $`<0.1`$ mag difference, was found between our $`K_s`$-magnitudes and the 2MASS $`K_s`$-magnitudes for bright stars of $`K_s15`$ mag.

The full width at half maximum (FWHM) of the seeing in the $`K_s`$-images is 0.8–1.0 arcsec, and in the $`R`$-images, typically 1.0–1.5 arcsec. The raw number counts of galaxies with detection significance $`3\sigma `$ turns over at $`K_s20`$ mag, although this varies somewhat from field to field due to the different exposure times. In the $`R`$-band the corresponding turn-over in the number counts is $`R24.5`$ mag. The 3$`\sigma `$ detection limits in the images are $`K_s21.5`$ and $`R26`$ mag. We address the completeness of the images in the next section.

## 3 Construction of object catalogs

### 3.1 Source detection

For object detection and photometry we used the sextractor software package (Bertin & Arnouts, 1996). In order to have the same FWHM of the seeing in both $`R`$ and $`K_s`$ we convolved the $`K_s`$-band images with a Gaussian function using the iraf task gauss prior to object detection. For each field, object catalogs were constructed based on the detections in the $`K_s`$-band image. During the detection process, we used a Gaussian convolution filter with a FWHM matched to the seeing and a detection threshold of 1.5$`\sigma `$. Since the images were block replicated, we set the minimum number of pixels above the threshold to 20 in order to avoid spurious detections and cosmic rays over an area of less than 5 pixels in the original images. The result of the detection process was inspected visually in order to ensure that no obvious objects were missed, and that no false detections were entered into the catalogs. Saturated objects and objects lying close to the image boundaries were rejected from the catalogs.

Because of the dithering technique used in the $`K_s`$-band, the exposure time falls off toward the edges in the final combined images. Since the PFIRCAM field of view is relatively small, we included lower exposure regions, but rejected the regions where the exposure time is less than 25 % of the total, i.e. where the signal-to-noise drops by more than 50 %. This causes the rms background noise to increase toward the edges, resulting in a number of spurious detections. In order to avoid such spurious detections, we used the exposure maps produced by dimsum as weights in sextractor. This allows for an adjustment of the detection threshold according to the noise level.

After the $`K_s`$-catalogs were made, matched $`R`$-band catalogs were created by running sextractor in double image mode. This produces a catalog of fluxes measured in the $`R`$-band image, but using the $`K_s`$-band image for detection. For total magnitudes, we used the sextractor best magnitudes which are evaluated using adaptive apertures, or using isophotal corrections if a neighbor is suspected of biasing the adaptive aperture magnitude (Bertin & Arnouts, 1996). For determining colors, we used aperture magnitudes, and chose an aperture diameter of two times the FWHM of the seeing, corresponding to typically 9–10 pixels, or 2.2–2.5 arcsec.

For each $`K_s`$-selected object, we calculated the detection significance using the uncertainty in the flux as measured by sextractor which takes into account the Poissonian nature of the counts and the standard deviation in the background counts. In addition to this, we added in quadrature the error arising from the fluctuating background level across the image (Best, 2000). In order for a $`K_s`$-selected object to be included in the catalog we required the detection significance within the aperture to be $`5\sigma `$. If the flux measured within the corresponding aperture in the $`R`$-band was significant at the $`3\sigma `$ level, we evaluated the $`RK_s`$ color, and if $`<3\sigma `$, a lower limit to the $`RK_s`$ color was calculated assuming a $`3\sigma `$ detection in $`R`$.

In order to separate stars and galaxies, we used the star-galaxy classifier in sextractor which assigns a probability between 0 and 1 to an object according to how likely it is to be a star. The star-galaxy classifier works well, but becomes less reliable at faint magnitudes. We therefore performed star-galaxy separation only at $`K_s17`$ mag, rejecting objects with a star-galaxy class of $`0.85`$ from the catalogs. The star-galaxy separation was based on the detections in the $`K_s`$ images.

### 3.2 Completeness

To be able to address the completeness of the $`K_s`$-band catalogs, we added artificial galaxies to the images using the iraf task mkobjects. For a given input magnitude, we added galaxies with both deVaucouleurs and exponential disk profiles, in total $`1500`$ galaxies distributed over the 13 fields. The scale-lengths of the artificial galaxies were randomly selected from the FWHM of the objects in the catalogs, and the images convolved with the same Gaussian function as was used to put the $`K_s`$-band images to the seeing of the $`R`$-band images. After this, the images were processed through sextractor using the same detection criteria as described above.

The simulations show that for input objects of $`K_s=18`$ mag, 96 % of the galaxies are recovered, where the missing 4 % is due to blending with other objects. At $`K_s=19`$ mag the completeness drops to 90 %, and at $`K_s=20`$ mag the catalogs are $`40`$ % complete. The numbers were evaluated assuming a 50:50 mix of deVaucouleurs and exponential disk profiles.

## 4 Results

### 4.1 The surface density of EROs

We select EROs in the quasar fields using the criterion $`RK_s5.0`$ mag. This includes old, passively evolving galaxies at $`z0.9`$, but also dusty starbursting galaxies. With only two filters, and at the resolution of our images, it is not possible to distinguish between these two galaxy types. There will also be some contamination by faint, red stars in our ERO sample as discussed below. We have chosen the $`RK_s5.0`$ mag selection criterion because of two reasons. Firstly, because selecting galaxies with red colors enables us to eliminate a large part of the numerous foreground galaxies at $`z<1`$ while at the same time producing a big enough sample to perform meaningful statistics on. Secondly, because several random-field surveys employing the $`RK_s5.0`$ mag criterion exist with which we can compare the ERO surface density. Our areal coverage in $`R`$ and $`K_s`$ is limited by the small field of view of the PFIRCAM instrument, so the field counts cannot be estimated from the quasar images themselves.

The color-magnitude diagram of all the $`K_s`$-selected objects is shown in Fig. 1. Above our $`RK_s5.0`$ mag selection there are 253 objects, and 105 of these have a detection significance in $`R`$ of three or less, i.e. the lower limit on their $`RK_s`$ color is $`5.0`$ mag. The sextractor star-galaxy classifier finds that two of the EROs at $`K_s17`$ mag are stars, and these are excluded from the ERO sample, resulting in a total of 251 galaxies (including the $`3\sigma `$ detections in $`R`$ having a lower limit of $`RK_s5.0`$ mag).

In Fig. 2, we have plotted the average cumulative surface density of EROs in the quasar fields together with ERO surface densities found in the random-field surveys of Daddi et al. (2000a), Scodeggio & Silva (2000), Cimatti et al. (2002), Roche et al. (2002), and Miyazaki et al. (2002). Note that we have not compared with field counts derived from other colors (e.g. $`IK`$ or $`IH`$) since they may bias toward later-type galaxies with prolonged star formation (McCarthy et al., 2001; Yan & Thompson, 2003). Note also that Miyazaki et al. define EROs as having $`RK_s(AB)3.35`$ mag, which, according to Miyazaki et al., corresponds to $`RK_s(Vega)4.95`$ mag.

The literature surveys cover wide, random areas and should therefore be representative of the ERO counts in the general field. The field counts agree well with each other except for the counts in the Chandra Deep Field South by Scodeggio & Silva (2000), which are consistently lower than the other. This is most likely caused by the inhomogeneity of the ERO distribution coupled with the relatively small (and contiguous) area surveyed by Scodeggio & Silva. Fig. 2 demonstrates that the ERO counts in our quasar fields are systematically higher than those in the field. At $`K_s19`$ mag, where our counts are complete, the average surface density of EROs is two–three times larger than the average of the field values (excluding those in the Chandra Deep Field South). According to the average of the random-field surveys, as indicated by the dotted line in Fig. 2, we expect $`34\pm 14`$ EROs at $`K_s19`$ mag over an area equal to the total area of our survey. (The error in this number takes into account the strong clustering of EROs as explained below.) The observed number of EROs in the quasar fields at $`K19`$ mag is 79, hence there is a significant excess above the field counts. At $`K_s20`$ mag we still see a clear excess in the quasar fields, but our counts become strongly affected by incompleteness.

The number of EROs per quasar field and the average cumulative surface density as a function of $`K_s`$-magnitude is listed in Table 2. Some fields are seen to be very rich in EROs, such as 1456$`+`$092 and 2150$`+`$053, with surface densities four–five times larger than the field values. These fields have ERO excesses in every magnitude bin. Other fields have surface densities consistent with the field counts, like 0741$`+`$169 and 1629$`+`$680. The majority of the fields, however, have surface densities in between these two extremes. This is illustrated in Fig. 3, where a histogram of the ERO surface densities at $`K_s19`$ mag is shown. If the quasar fields were random with respect to the ERO distribution, we would expect a normal distribution centered on $`0.546`$ arcmin-2, which is the average surface density in the general field. This is not seen, instead there is a tail toward high surface densities, with as many as four fields (30 per cent of the sample) having a surface density larger than three times the field value, and seven fields where the surface density is twice the field value. The combination of the very rich and the moderately rich fields gives an overall surface density well above the average for the random-field surveys.

### 4.2 Contamination by red stars

Of the 13 quasars fields, 10 lie at high galactic latitudes $`b>35^{}`$, so the contamination by stars is likely to be small. Using a color-cut of $`RK_s5.0`$ mag will eliminate most field stars because of their bluer color, but L- and M-type dwarfs will be red enough to make it into our ERO sample. Based on $`RK^{}`$ and $`JK^{}`$ colors, Mannucci et al. (2002) find that 9 % of their $`K^{}20`$ mag ERO sample (at $`b=23^{}`$) is likely to consist of stars.

In order to estimate the contamination by stars, we used the stellar distribution model developed by Jarrett (1992) (see also Jarrett, Dickman, & Herbst (1994)). The model is an extension of the Bahcall & Soneira (1980) optical star count model to the near-infrared, and its performance with respect to 2MASS near-infrared star counts is demonstrated by Cambrésy et al. (2002). For every quasar position, we used the model to predict the expected number of stars having $`RK_s5.0`$ mag. The predicted average surface density of stars for the whole ERO sample is listed on the bottom line of Table 2. These numbers show that 40–50 % of the brighter EROs at $`K_s18`$ mag might be stars, but that at fainter magnitudes the contamination decreases to 15–20 %. Note that this is an overestimate because the ERO counts suffer from incompleteness at fainter magnitudes. The 9 % contamination estimated by Mannucci et al. (2002) based on two colors therefore seems to agree well with the model predictions. We thus expect roughly 10 % of our EROs to be red stars. Since we cannot use the small quasar fields for estimating the field counts of EROs, we compare the counts in the quasar fields with number densities of EROs from random-field surveys found in the literature. Since the random-field surveys are also contaminated by stars, a contamination of 10–20 % does not alter our conclusions.

The three lowest galactic latitude fields in our sample are 0741$`+`$169, 0758$`+`$120 and 0830$`+`$115 at $`b=2027^{}`$. For these fields the model predicts 2.9, 2.6, and 1.6 stars at $`K_s19`$ mag, respectively. Note also that these are not among the richest fields in our sample.

### 4.3 Significance of the ERO excess

In this section we evaluate the significance of the ERO excess, taking into account the large variance in the ERO field counts caused by their strong clustering. The significance of the excess is given by

$$\sigma _{exc}=\frac{(N_{qso}N_{field})}{\sigma }$$
(1)
(Yee & López-Cruz, 1999), where $`N_{qso}`$ is the total number of EROs in the quasar fields brighter than a given $`K_s`$-magnitude and $`N_{field}`$ is the number of EROs in the random field scaled to the total area of the quasar fields, 61.7 arcmin2. The error in the field counts expected over an area of 61.7 arcmin2 is denoted $`\sigma `$. Ideally $`N_{field}`$ should have been estimated from the edges of the quasar fields, but our fields are too small for this. Instead we use the average of the surface densities found by Daddi et al. (2000a), Roche et al. (2002), and Miyazaki et al. (2002).

Because ERO positions are correlated on the sky, the variance is given by

$$\sigma ^2=N_{field}\left(1+N_{field}A_\omega C\right)$$
(2)
(Daddi et al., 2000a), where $`A_\omega `$ is the amplitude of the ERO angular two-point correlation function. The integral constraint which arises from the fact that the survey area is limited is given by $`A_\omega C`$, and Daddi et al. find that $`C`$ can be approximated by $`58\times Area^{0.4}`$. We use this approximation and find that $`C=11.15`$ for the 13 quasar fields. The approximation is valid for a contiguous area, but serves our purpose well even if the 13 fields are non-contiguous. In any case, $`C=11.15`$ is likely to be an overestimate since our fields are widely separated and not correlated in any way. We use the amplitudes, $`A_\omega `$, found by Daddi et al. and Roche et al. for $`K_s19`$ mag and by Roche et al. for $`K_s19.5`$ mag, and calculate $`\sigma `$.

The result of the calculations for the entire quasar sample on average is shown in Table 3, where it can be seen that the excess at $`K_s18`$, 18.5 and 19 mag is significant at the $`3\sigma `$ level. Because of incompleteness, the excess is not significant at fainter $`K_s`$-levels, but if a correction is made for the incompleteness, the excess appears to be significant also at fainter $`K_s`$ magnitudes, as shown in the last column of Table 3.

### 4.4 Radial distribution

If the EROs are associated with the quasars either through a surrounding cluster or a magnification bias we might expect to see a trend with projected quasar separation. In order to examine the radial distribution of EROs, we evaluated the average surface density in three annuli centered on the quasars, for three different magnitude intervals. The result is shown in Fig. 4. The brighter EROs at $`17.5K_s<19.5`$ mag are seen to be uniformly distributed over the whole field, but the fainter ones at $`K_s19.5`$ mag seem to cluster more in the immediate 20 arcsec region around the quasars. Note that this is not caused by the mosaicing technique in the $`K_s`$ filter. The exposure times across the inner two annuli are uniform, with the third annulus being on the average only slightly less well exposed than the other two. Using completeness simulations we find that the variation in the completeness level across the three annuli is at most a few per cent, i.e. consistent with what we expect from object crowding and blending. Therefore, there appears to be a real trend toward central concentration of the fainter EROs. Cimatti et al. (2000) looked for central concentrations of $`K_s`$-selected galaxies toward $`z1.5`$ radio-loud AGN. They did not find any significant trends, but noted that the regions closest to the AGNs had a higher number density of galaxies than in the general field.

### 4.5 The surface density of $`RK_s6.0`$ mag EROs

For completeness, we also evaluated the surface density of $`RK_s6.0`$ mag EROs. Using the redder color criterion gives a much smaller sample, 14 EROs at $`K_s19`$ mag and 31 EROs at $`K_s20`$ mag. The resulting surface densities are thus $`0.23\pm 0.06`$ and $`0.50\pm 0.09`$ arcmin-2 at $`K_s19`$ and $`20`$ mag, respectively (Poissonian errors). The surface density of EROs at $`K_s19`$ mag is therefore 3–5 times higher than that found in the field surveys by Thompson et al. (1999) and Daddi et al. (2000a) ($`0.039\pm 0.016`$ and $`0.07\pm 0.01`$ arcmin-2, respectively). We also note that the surface density of $`RK_s6.0`$ EROs in our quasar fields agrees with the ERO surface densities found in radio-loud AGN fields at $`z=11.5`$ by Hall & Green (1998) and Cimatti et al. (2000).

### 4.6 $`K_s`$-band counts

Whereas there is a clear excess of EROs in the quasar fields, the overall $`K_s`$-counts show only a slight excess at fainter magnitudes, $`K_s18.5`$ mag. This is seen in Fig. 5 where we have plotted the average $`K_s`$-band number counts together with counts from the surveys by Soifer et al. (1994), Moustakas et al. (1997), and Maihara et al. (2001), as well as the literature averages compiled by Hall, Green, & Cohen (1998) and Best et al. (2003). There are large field-to-field variations in the random-field $`K_s`$-band counts, and even though the counts in the quasar fields appear to lie above the literature averages at fainter $`K_s`$ magnitudes, they are still consistent with the counts from the surveys by e.g. Soifer et al. (1994), Moustakas et al. (1997), and Maihara et al. (2001). We also find that the $`R`$-band number counts are consistent with the counts in the general field (Metcalfe et al., 1995; Hogg et al., 1997). Overplotted in Fig. 5 are also the number counts of EROs in our fields, demonstrating, as expected, that non-ERO galaxies clearly dominate at all magnitudes.

Table 4 shows the fraction of EROs with respect to the whole $`K_s`$-selected sample of galaxies. For comparison, we have also listed the ERO fractions given by Daddi et al. (2000a), Scodeggio & Silva (2000), and Roche et al. (2002), the latter calculated from their table 2. In the quasar fields, the EROs make up typically 20–25 % of all the $`K_s`$-selected galaxies, whereas the field surveys find $`10`$ %. We note that this behavior is what we expect if the EROs are tracing the early-type population in high-redshift clusters of galaxies since the fraction of early-type galaxies in rich clusters is known to be different from the field (Dressler, 1980). This also illustrates the richness of EROs in the quasar fields, and suggests that the EROs are efficiently tracing overdensities not readily apparent from the $`K_s`$ number counts alone. Clearly, the color-selection is very effective in eliminating foreground contamination.

## 5 Discussion

We have demonstrated above that the quasar fields have a significant excess of EROs above what is expected from the general field population of EROs, taking into account their strong clustering. It appears that the quasar fields contain clusters or dense concentrations of red galaxies, and their radial distribution suggests that the fainter EROs may be related to the quasars. However, another possibility is that the EROs are related to structures in the foreground of the quasars. We discuss each of these possibilities below. Note that our sample is dominated by radio-loud quasars, and that we restrict the following discussion to radio-loud quasars only.

### 5.1 ‘Clusters’ of EROs associated with the quasars

An argument in favor of the ERO excess being at the quasar redshifts is that luminous radio-loud quasars and radio galaxies at $`z1`$ are known to reside in regions of enhanced galaxy density, interpreted as clusters or groups of galaxies physically associated with the AGN (Hall & Green, 1998; Best, 2000). Since powerful AGN with supermassive black holes may trace the densest regions, it is natural to assume that this is true also at $`z2`$. For example, Cimatti et al. (2000) find excess EROs in radio-loud AGN fields at $`z1.5`$ and Best et al. (2003) find that red galaxies cluster around powerful radio sources at $`z1.5`$. Our study pushes this out to $`z2`$.

If the ERO excess originates in dense environments associated with the quasars, the majority of the EROs must lie at $`z2`$ since this is the redshift of most of the quasars in our sample. It is therefore appropriate to discuss whether the color and the brightness of the EROs are consistent with being evolved galaxies at $`z2`$. We have used the pégase models of Fioc & Rocca-Volmerange (1997) to examine both the apparent magnitude and the $`RK_s`$ color of passively evolving galaxies at $`z2`$ (assuming no extinction by dust, and an initial metallicity of 1 % solar value). The apparent $`K_s`$ magnitude of $`L^{}`$ galaxies at $`z=2`$ was found by calculating the scale factor between the observed flux of the model and that of a present-day $`L^{}`$ galaxy. The present-day $`L^{}`$ galaxy was assumed to have $`M_B^{}=20.3`$ mag (Loveday et al., 1992) and an age of 14 Gyr.

The standard model in which the stars form in a 1 Gyr long burst at $`z>4`$ has an $`RK_s`$ color at $`z=2`$ which is too red for most of our EROs. In order for this model to produce a color of $`RK_s=56`$ mag at $`z=2`$, the formation redshift has to be $`z3.5`$. The $`K_s^{}`$ magnitude of such a galaxy observed at $`z=2`$ is 19.2 mag. Since the star formation in this model stops abruptly after 1 Gyr, this model represents the galaxy which reddens most quickly. Therefore, $`K_s=19.2`$ mag is the brightest magnitude we can expect for an $`L^{}`$ galaxy at $`z=2`$. If the $`K_s19`$ mag ERO excess is associated with passively evolving galaxies at $`z2`$, they must therefore be $`L^{}`$ galaxies. Other models which have $`RK_s=56`$ mag at $`z=2`$ include those with exponentially decreasing star formation rates. For instance, a galaxy formed at $`z=6`$ with an exponentially decreasing star formation rate of $`\tau =0.5`$ Gyr has $`RK_s=56`$ mag at $`z=2`$, and $`K_s^{}=19.619.7`$ mag. If the EROs are described by this latter model, they must have small amounts of residual star formation in them.

The recent spectroscopic work of van Dokkum et al. (2003) has shown that red galaxies at $`z>2`$ can be very bright. Five of their six galaxies were found to lie at $`2.4z3.5`$ and to have $`K_s`$-magnitudes in the range 19.2–19.9 mag. These galaxies were selected by the $`JK_s>2.3`$ mag criterion (Franx et al., 2003), whereas our galaxies are selected on the basis of having $`RK_s5.0`$ mag. Both criteria are designed to find galaxies with evolved 4000 Å breaks, but at different redshifts. However, these two selection techniques will not necessarily sample the same type of galaxies since the different filters sample different rest-frame wavelength ranges. While the results of van Dokkum et al. prove that high-redshift red galaxies can be bright, their spectroscopic sample (six of eleven $`JK_s`$-selected galaxies) may be biased toward systems with ongoing star formation.

The average ERO fraction in our quasar fields is above 20 %, and for the four richest fields in our sample the fraction is as high as $`50`$ %. This is what we expect if the fields contain clusters or concentrations of galaxies with an already existing relation between morphology and density (Dressler et al., 1997). This is also the suggestion from the work by Best et al. (2003) who find that red galaxies are concentrated toward $`z1.5`$ powerful radio sources.

Whereas the magnitudes and colors of the EROs are consistent with models of passively evolving galaxies, another possibility is that they might be starbursting systems reddened by dust. Hence, an alternative explanation is that the quasars live in dense environments dominated by galaxies with enhanced star formation. In this case, the excess EROs may be dusty starburst galaxies, possibly the progenitors of massive ellipticals in a forming cluster. There are examples of fields around high-redshift radio-loud AGN with excesses of starforming galaxies such as Ly-alpha emitters (Kurk et al., 2000; Pentericci et al., 2000; Venemans et al., 2002). Also, Ivison et al. (2000) find a concentration of submm sources in the field of the $`z=3.8`$ radio galaxy 4C 41.17 raising the possibility that the EROs in this field are dusty starburst systems rather than passively evolving galaxies. This possibility is also discussed by Smail et al. (2003) for a luminous submm source in the field of the radio galaxy 53W002 at $`z=2.4`$.

One might expect a correlation between the radio properties of the quasars and their ERO excesses if there is e.g. a relation between radio power and black hole mass (Wold et al., 2000; Laor, 2000; Lacy et al., 2001). We have checked this for our sample, but do not find that the radio power correlates with the ERO excess.

### 5.2 EROs in foreground lensing structures

Little is known yet of the luminosity function of EROs. In the spectroscopic survey by Cimatti et al. (2002), which is 67 % complete at $`K_s19.2`$ mag, the $`RK_s5.0`$ mag EROs are found to lie at $`z=0.71.5`$. The redshift distribution derived by Miyazaki et al. (2002) based on photometric redshifts is in broad agreement with Cimatti et al’s, but has a tail toward higher redshifts. If we assume that our ERO sample has a redshift distribution similar to that found by Cimatti et al. and Miyazaki et al., we have to explain why the ERO surface density in the quasar fields is so significantly different from that in the large random-field surveys. We are sampling the ERO distribution with several small, but widely separated fields, and should therefore be minimally affected by cosmic variance. Instead, the opposite seems to be the case; the quasar fields sample the densest regions in the ERO distribution. A large number of the quasar fields contain clumps of EROs, similar to what is seen only in the densest areas in the large random-field surveys. If the EROs have a peak redshift of $`z1.0`$, they must therefore be clusters of galaxies, or overdensities in the ERO distribution, lying along the lines of sight to the quasars.

A biasing of quasar lines of sight with respect to $`z1`$ overdensities in the ERO distribution can be explained by gravitational lensing. If we assume that the overdensities in the quasar fields do not extend much beyond the edges of the fields their typical scale is $`1`$ arcmin, corresponding to 2–3 Einstein radii for galaxy clusters at $`z1`$. The most efficient lens redshift for a quasar at $`z=2`$ is 0.5–0.6, but a magnification can still occur if the lens lies at $`z0.81.0`$.

The quasars in our sample are very luminous, eight of them have $`M_V<27.2`$ mag ($`H_0=50`$ km s-1 Mpc-1), and hence classify as highly luminous quasars. They therefore fulfill the three criteria for maximizing the probability of a magnification bias, i.e. they are distant, bright and have a steep number-magnitude slope (Claeskens & Surdej 2002, and references therein). If a magnification bias is the explanation for our results, we might expect that the fields around the most luminous quasars contain the highest overdensities of EROs, but our data do not show any significant correlations between quasar luminosity and ERO surface density to support this.

Furthermore, if the majority of the EROs are associated with overdensities along the line of sight to the quasars, we might expect the spectra of the quasars in the richest fields to have strong (rest-frame equivalent width $`0.7`$ Å) metal absorption line systems (Liu et al., 2000). We have found absorption line data for seven of the 13 quasars in our sample (Junkkarinen, Hewitt, & Burbidge, 1991, 1992). Two of the four richest fields have such data; Q2150$`+`$053 show no strong metal absorption lines and Q2338$`+`$042 have CIV and FeII at $`z_{abs}1.80`$ (Barthel, Tytler, & Thomson, 1990). Other quasars with metal absorption line systems in the foreground is Q0830$`+`$115 with FeII at $`z_{abs}=0.92`$ and MgII at $`z_{abs}=0.80`$ (Sargent, Steidel, & Boksenberg, 1989), Q0941$`+`$261 with MgII systems at $`z_{abs}=0.71`$ and $`z_{abs}=1.09`$, and Q1354$`+`$258 with MgII at $`z_{abs}=0.86`$ and some weaker FeII and MgII systems at $`z_{abs}=1.42`$ (Barthel, Tytler, & Thomson, 1990), but these quasars lie in fields with relatively few EROs compared to the richer fields. The data for Q0758$`+`$120 and Q1629$`+`$680 show no strong metal absorption line systems (Barthel, Tytler, & Thomson, 1990). There is thus no clear trend of increasing number of absorption line systems with ERO richness. Also, the quasars in our sample were not selected on the basis of having strong or numerous absorption line systems.

## 6 Conclusions

On the average, the quasar fields studied here have a cumulative surface density of EROs at $`K_s18`$, 18.5 and 19 mag which is two-three times larger than what is found in the general field. The excess is significant at the $`3\sigma `$ level even after taking into account that the ERO distribution is strongly clustered. There is also an excess of EROs at fainter magnitudes, but at $`K_s20`$ mag our counts become strongly affected by incompleteness and a secure quantitative estimate of the excess for $`K_s>20`$ mag is difficult. However, by making a correction for the incompleteness, we find that the ERO excess at $`K_s19.5`$ and $`K_s20`$ mag may also be significant at the $`3\sigma `$ level.

The quasar fields are small, typically 2.1$`\times `$2.1 arcmin2 (limited by the near-infrared field of view), but widely separated, and therefore sample the ERO population with minimal influence from cosmic variance. The fact that the excess is statistically significant therefore suggests that the quasar fields probe the densest regions in the ERO distribution, implying that the fields around luminous, high-redshift quasars are on the average, overdense in EROs. There are two possible explanations for the observed overdensities. (1) The majority of the EROs may be associated with concentrations of galaxies at the quasar redshifts, or (2) the EROs lie in the foreground, along the line of sight to the quasars.

Based on both spectroscopic observations of $`z>2`$ red galaxies (van Dokkum et al., 2003) and spectral synthesis models (Fioc & Rocca-Volmerange, 1997), we argue that a redshift of $`z2`$ for the EROs is plausible. The $`K_s19`$ mag ERO excess may therefore be physically associated with the quasars, in which case their magnitudes and colors are consistent with $`L^{}`$ passively evolving galaxies.

The EROs at $`K_s1919.5`$ mag are more clustered in the immediate 20 arcsec region around the quasars, hence the radial distribution supports the explanation that the fainter EROs are physically associated with the quasars. Furthermore, the fraction of EROs among the $`K_s`$-selected galaxies in the quasar fields is unusually high, typically 20–25 %, in contrast to $`10`$ % which is found in random-field surveys (Daddi et al., 2000a; Scodeggio & Silva, 2000; Roche et al., 2002). This is consistent with the ERO excess being associated with dense concentrations of red, possibly evolved, galaxies. The high ERO fractions may therefore indicate that a morphology-density relation (Dressler et al., 1997) exists in the quasar fields. The candidate ERO clusters in the quasar fields may therefore represent the oldest and most evolved systems at $`z2`$, and will be crucial for testing models of structure and galaxy formation. If the ERO excess is dominated by dusty starforming galaxies instead of passively evolving galaxies, the EROs might be the progenitors of massive ellipticals in a forming cluster associated with the quasars.

The second explanation, that the ERO excess lies in the foreground, is equally interesting because it implies that the lines of sight toward high-redshift, luminous quasars are biased toward overdensities in the ERO distribution. In this case, the ERO excess is most likely associated with clusters of galaxies at $`z1`$ giving rise to a magnification bias of the quasars. In this interpretation, we would be seeing candidate clusters with passively evolving galaxies at $`z1`$, a redshift where still only less than a handful of clusters are known. However, we see no trend of quasar brightness with ERO overdensity as might be expected in this interpretation. Also, there is no association of the numbers of absorption line systems with ERO richness for seven of the 13 quasars with absorption line data.

With the current data, it is not possible to tell whether the ERO excess is part of a small-scale overdensity centered on, or lying close to, the quasars, or if it is part of a larger structure that extends over a wider area. To address this issue, we have begun an imaging program using $`J`$\- and $`K_s`$-filters at the Palomar 200-inch telescope over much larger fields. With these data we will be able to trace the ERO population around the $`z2`$ quasars on scales up to 4–5 Mpc.

We are grateful to Dave Thompson, Lin Yan and Mark Lacy for helpful discussions, and to the referee for a careful review of the manuscript. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has also 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.


format: text

text:
# Random matrices with external source and multiple orthogonal polynomials

## 1 Random matrices with external source

Following Brézin and Hikami  and P. Zinn-Justin  we consider a random matrix ensemble with an external source,

$$\frac{1}{Z_n}e^{\mathrm{Tr}(V(M)AM)}dM$$
(1.1)
defined on $`n\times n`$ Hermitian matrices $`M`$. The ensemble (1.1) consists of a general unitary invariant part $`V(M)`$ and an extra term $`AM`$ where $`A`$ is a fixed $`n\times n`$ Hermitian matrix, the external source or the external field. Due to the external source, the ensemble (1.1) is not unitary invariant. For the special Gaussian case $`V(x)=\frac{1}{2}x^2`$, we can write $`M`$ in (1.1) as $`M=H+A`$ where $`H`$ is a random matrix from the GUE ensemble, and $`A`$ is deterministic, hence in this case it reduces to the class of deterministic plus random matrices studied in .

Zinn-Justin  showed that the eigenvalue correlations of ensemble (1.1) can be expressed in the determinantal form,

$$R_m(\lambda _1,\mathrm{},\lambda _m)=det(K_n(\lambda _i,\lambda _j))_{i,j=1,\mathrm{},m}$$

for some kernel $`K_n`$. In this paper, we give a different expression for $`K_n`$. We believe that our formulation is useful for asymptotic analysis. Indeed, for the Gaussian case $`V(x)=\frac{1}{2}x^2`$ and for the case where $`A`$ has only two distinct eigenvalues, we have been able to carry out the asymptotic analysis almost completely. This will be reported elsewhere.

Our approach is based on the observation that the average characteristic polynomial

$$P_n(z)=𝔼\left[det(zM)\right]$$

of the ensemble (1.1) can be characterized by the property that

$$_{\mathrm{}}^{\mathrm{}}P_n(x)x^ke^{(V(x)a_jx)}𝑑x=0$$

for every eigenvalue $`a_j`$ of $`A`$ and for $`k=0,\mathrm{},n_j1`$, where $`n_j`$ is the multiplicity of $`a_j`$, see Section 2. We can embed the polynomial $`P_n`$ in a sequence of polynomials $`\{P_k\}_0^n`$ where $`P_k`$ has degree $`k`$. Then our kernel has the form

$$K_n(x,y)=e^{\frac{1}{2}(V(x)+V(y))}\underset{k=0}{\overset{n1}{}}P_k(x)Q_k(y)$$
(1.2)
where the $`Q_k`$ are certain dual functions (not polynomials in general), see Section 3.

When $`A=0`$ (no external source), the polynomials $`P_k`$ are usual monic orthogonal polynomials with respect to the weight $`e^{V(x)}`$ on $``$. In that case, the function $`Q_k`$ is a multiple of $`P_k`$ and the kernel (1.2) reduces to the orthogonal polynomial kernel which is familiar in the theory of random matrices. By the Christoffel-Darboux formula we then have

$$K_n(x,y)=e^{\frac{1}{2}(V(x)+V(y))}\frac{\gamma _{n1}}{\gamma _n}\frac{P_n(x)P_{n1}(y)P_{n1}(x)P_n(y)}{xy},$$
(1.3)
where $`\gamma _k`$ is the leading coefficient of the orthonormal polynomial of degree $`k`$.

In Section 4 we present an analog of the Christoffel-Darboux formula for the kernel (1.2) in the case where $`A`$ has only two eigenvalues. We also relate it to a Riemann-Hilbert problem in Section 5.

## 2 The average characteristic polynomial

We define the monic polynomial

$$P_n(z)=𝔼\left[det(zM)\right]$$

where the expectation is with respect to the ensemble (1.1).

###### Proposition 2.1.

Suppose $`A`$ has eigenvalues $`a_j`$, $`j=1,\mathrm{},n`$ with $`a_ia_j`$ if $`ij`$. Then the following hold.

1. There is a constant $`\stackrel{~}{Z}_n`$ such that

$$P_n(z)=\frac{1}{\stackrel{~}{Z}_n}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{(V(\lambda _j)a_j\lambda _j)}\mathrm{\Delta }(\lambda )d\lambda ,$$
(2.1)
where

$$\mathrm{\Delta }(\lambda )=\underset{i>j}{}(\lambda _i\lambda _j)$$

and $`d\lambda =d\lambda _1d\lambda _2\mathrm{}d\lambda _n`$.
2. Let

$$m_{jk}=_{\mathrm{}}^{\mathrm{}}x^ke^{(V(x)a_jx)}𝑑x.$$

Then we have the determinantal formula

$$P_n(z)=\frac{1}{\stackrel{~}{Z}_n}\left|\begin{array}{cccc}m_{10}& m_{11}& \mathrm{}& m_{1n}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{n0}& m_{n1}& \mathrm{}& m_{nn}\\ 1& z& \mathrm{}& z^n\end{array}\right|.$$
(2.2)
3. For $`j=1,\mathrm{},n`$,

$$_{\mathrm{}}^{\mathrm{}}P_n(x)e^{(V(x)a_jx)}𝑑x=0,$$
(2.3)
and these equations uniquely determine the monic polynomial $`P_n`$.

###### Proof.

Write $`M=U\mathrm{\Lambda }U^{}`$ where $`U`$ is unitary and $`\mathrm{\Lambda }`$ is diagonal, $`\mathrm{\Lambda }=\mathrm{diag}(\lambda _1,\mathrm{},\lambda _n)`$. of $`M`$ on the diagonal. Then by the Weyl integration formula, see e.g. , we have for every integrable function $`f`$ on the space of Hermitian $`n\times n`$ matrices,

$$f(M)𝑑M=\pi ^{n(n1)/2}\left(\underset{j=0}{\overset{n}{}}j!\right)f(U\mathrm{\Lambda }U^{})\mathrm{\Delta }(\lambda )^2𝑑\lambda 𝑑U$$
(2.4)
where $`dU`$ denotes the normalized Haar measure on the unitary group $`U(n)`$. Thus

$$P_n(z)=\frac{\left(_{j=0}^nj!\right)}{Z_n\pi ^{n(n1)/2}}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{V(\lambda _j)}\left(e^{AU\mathrm{\Lambda }U^{}}𝑑U\right)\mathrm{\Delta }(\lambda )^2d\lambda .$$

Because of the Harish-Chandra, Itzykson–Zuber integral

$$e^{AU\mathrm{\Lambda }U^{}}𝑑U=\left(\underset{j=0}{\overset{n1}{}}j!\right)\frac{det\left(e^{a_j\lambda _k}\right)}{\mathrm{\Delta }(a)\mathrm{\Delta }(\lambda )},$$

we obtain that

$$P_n(z)=\frac{n!\left(_{j=0}^{n1}j!\right)^2}{Z_n\pi ^{n(n1)/2}}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{V(\lambda _j)}\frac{det\left(e^{a_j\lambda _k}\right)}{\mathrm{\Delta }(a)}\mathrm{\Delta }(\lambda )d\lambda .$$
(2.5)
We expand the determinant

$$det\left(e^{a_j\lambda _k}\right)=\underset{\sigma S_n}{}(1)^\sigma \underset{j=1}{\overset{n}{}}e^{a_j\lambda _{\sigma (j)}}$$

where $`S_n`$ is the symmetric group. Hence

$$P_n(z)=\frac{n!\left(_{j=0}^{n1}j!\right)^2}{Z_n\pi ^{n(n1)/2}\mathrm{\Delta }(a)}\underset{\sigma S_n}{}(1)^\sigma \underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{V(\lambda _j)}\underset{j=1}{\overset{n}{}}e^{a_j\lambda _{\sigma (j)}}\mathrm{\Delta }(\lambda )d\lambda .$$

We make a change of variables $`\lambda _j^{}=\lambda _{\sigma (j)}`$. Then $`(1)^\sigma \mathrm{\Delta }(\lambda )=\mathrm{\Delta }(\lambda ^{})`$, hence in the sum over $`S_n`$, we have $`n!`$ equal terms and, by dropping the prime, we obtain (2.1) with constant

$$\stackrel{~}{Z}_n=Z_n\frac{\pi ^{n(n1)/2}\mathrm{\Delta }(a)}{\left(_{j=0}^nj!\right)^2}.$$
(2.6)
This proves part (a).

Observe that

$$\underset{j=1}{\overset{n}{}}(z\lambda _j)\mathrm{\Delta }(\lambda )=\mathrm{\Delta }(\lambda ,z)=\left|\begin{array}{cccc}1& \lambda _1& \mathrm{}& \lambda _1^n\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 1& \lambda _n& \mathrm{}& \lambda _n^n\\ 1& z& \mathrm{}& z^n\end{array}\right|$$

and

$$\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{(V(\lambda _j)a_j\lambda _j)}\mathrm{\Delta }(\lambda )=\left|\begin{array}{cccc}e^{(V(\lambda _1)a_1\lambda _1)}& \lambda _1e^{(V(\lambda _1)a_1\lambda _1)}& \mathrm{}& \lambda _1^ne^{(V(\lambda _1)a_1\lambda _1)}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ e^{(V(\lambda _n)a_n\lambda _n)}& \lambda _ne^{(V(\lambda _n)a_n\lambda _n)}& \mathrm{}& \lambda _n^ne^{(V(\lambda _n)a_n\lambda _n)}\\ 1& z& \mathrm{}& z^n\end{array}\right|.$$

Then (2.2) follows immediately from this and (2.1). This proves part (b).

From (2.2) it follows that

$$_{\mathrm{}}^{\mathrm{}}P_n(x)e^{(V(x)a_jx)}𝑑x=\frac{1}{\stackrel{~}{Z}_n}\left|\begin{array}{cccc}m_{10}& m_{11}& \mathrm{}& m_{1n}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{n0}& m_{n1}& \mathrm{}& m_{nn}\\ m_{j0}& m_{j1}& \mathrm{}& m_{jn}\end{array}\right|=0,$$

for every $`j=1,\mathrm{},n`$. This proves (2.3). To prove uniqueness of $`P_n`$ satisfying (2.3), observe that by equating the coefficients of $`x^n`$ in (2.2) we obtain that

$$\left|\begin{array}{cccc}m_{10}& m_{11}& \mathrm{}& m_{1,n1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{n0}& m_{n1}& \mathrm{}& m_{n,n1}\end{array}\right|=\stackrel{~}{Z}_n0.$$
(2.7)
Let $`P_n(x)=x^n+p_{n1}x^{n1}+\mathrm{}+p_0`$ and set $`p=\left(\begin{array}{ccc}p_0& \mathrm{}& p_{n1}\end{array}\right)^T`$. Then the equations (2.3) are written in terms of the vector $`p`$ as

$$Mp=m,M=\left(m_{jk}\right)_{j=1,\mathrm{},n;k=0,\mathrm{},n1},m=(m_{jn})_{j=1,\mathrm{},n}.$$

By (2.7), $`detM0`$, hence $`p`$ and therefore $`P_n`$ is unique. ∎

Proposition 2.1 can be extended to the case of multiple $`a_j`$s as follows.

###### Proposition 2.2.

Suppose $`A`$ has distinct eigenvalues $`a_i`$, $`i=1,\mathrm{},p`$ with respective multiplicities $`n_i`$ so that $`n_1+\mathrm{}+n_p=n`$. Let $`n^{(i)}=n_1+\mathrm{}+n_i`$ and $`n^{(0)}=0`$. Define

$$w_j(x)=x^{d_j1}e^{(V(x)a_ix)},j=1,\mathrm{},n,$$

where $`i=i_j`$ is such that $`n^{(i1)}<jn^{(i)}`$ and $`d_j=jn^{(i1)}`$. Then the following hold.

1. There is a constant $`\stackrel{~}{Z}_n>0`$ such that

$$P_n(z)=\frac{1}{\stackrel{~}{Z}_n}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}w_j(\lambda _j)\mathrm{\Delta }(\lambda )d\lambda .$$
(2.8)
2. Let

$$m_{jk}=_{\mathrm{}}^{\mathrm{}}x^kw_j(x)𝑑x.$$

Then we have the determinantal formula

$$P_n(z)=\frac{1}{\stackrel{~}{Z}_n}\left|\begin{array}{cccc}m_{10}& m_{11}& \mathrm{}& m_{1n}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{n0}& m_{n1}& \mathrm{}& m_{nn}\\ 1& z& \mathrm{}& z^n\end{array}\right|.$$
(2.9)
3. For $`i=1,\mathrm{},p`$,

$$_{\mathrm{}}^{\mathrm{}}P_n(x)x^je^{(V(x)a_ix)}𝑑x=0,j=0,\mathrm{},n_i1,$$
(2.10)
and these equations uniquely determine the monic polynomial $`P_n`$.

###### Proof.

We write

$$\alpha =(\alpha _1,\alpha _2,\mathrm{},\alpha _n)=(\underset{n_1\text{ times}}{\underset{}{a_1,\mathrm{},a_1}},a_2,\mathrm{},a_{p1},\underset{n_p\text{ times}}{\underset{}{a_p,\mathrm{},a_p}})$$

Apply formula (2.5) in the case when all $`a_j=\stackrel{~}{a}_j`$ are different and take a limit to the multiple $`a_j`$’s. In this limit we have that

$$lim\frac{det\left(e^{\stackrel{~}{a}_j\lambda _k}\right)}{\mathrm{\Delta }(\stackrel{~}{a})}=\frac{det\left(\lambda _k^{d_j1}e^{\alpha _j\lambda _k}\right)}{\mathrm{\Delta }_0(a)_{i=1}^p_{k=1}^{n_i1}k!}$$

where $`d_j`$ is as in the statement of the proposition, and

$$\mathrm{\Delta }_0(a)=\underset{i>j}{}(a_ia_j)^{n_in_j}.$$

Thus, formula (2.5) becomes

$$P_n(z)=\frac{n!\left(_{j=0}^{n1}j!\right)^2}{Z_n\pi ^{n(n1)/2}}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{V(\lambda _j)}\frac{det\left(\lambda _k^{d_j1}e^{\alpha _j\lambda _k}\right)}{\mathrm{\Delta }_0(a)_{i=1}^p_{k=1}^{n_i1}k!}\mathrm{\Delta }(\lambda )d\lambda .$$
(2.11)
Then we continue as in the proof of Proposition 2.1, that is, we write

$$det\left(\lambda _k^{d_j1}e^{\alpha _j\lambda _k}\right)=\underset{\sigma S_n}{}(1)^\sigma \underset{j=1}{\overset{n}{}}\lambda _{\sigma (j)}^{d_j1}e^{\alpha _j\lambda _{\sigma (j)}},$$

and insert this into (2.11) to obtain a sum of $`n!`$ equal terms, which leads to (2.8) with

$$\stackrel{~}{Z}_n=Z_n\frac{\pi ^{n(n1)/2}\mathrm{\Delta }_0(a)_{i=1}^p_{k=1}^{n_i1}k!}{\left(_{j=0}^nj!\right)^2}.$$
(2.12)
This proves part (a).

Parts (b) and (c) follow from (2.8) in the same way as parts (b) and (c) of Proposition 2.1 followed from (2.1). Note that in particular we have as in (2.7),

$$\stackrel{~}{Z}_n=\left|\begin{array}{cccc}m_{10}& m_{11}& \mathrm{}& m_{1,n1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{n0}& m_{n1}& \mathrm{}& m_{n,n1}\end{array}\right|0.$$
(2.13)
Remark: Formula (2.8) can be also written in the following form:

$$P_n(z)=\frac{1}{\widehat{Z}_n}\underset{j=1}{\overset{n}{}}(z\lambda _j)\underset{j=1}{\overset{n}{}}e^{(V(\lambda _j)a_{i_j}\lambda _j)}\underset{i=1}{\overset{p}{}}\mathrm{\Delta }(\lambda ^{(i)})\mathrm{\Delta }(\lambda )d\lambda ,$$
(2.14)
where $`\lambda ^{(i)}=(\lambda _{n^{(i1)}+1},\mathrm{},\lambda _{n^{(i)}})`$ and

$$\widehat{Z}_n=\stackrel{~}{Z}_nn_1!\mathrm{}n_p!.$$
(2.15)
When $`A=0`$, (2.14) reduces to the usual formula for $`P_n(z)`$ with respect to the random matrix ensemble without external source.

###### Corollary 2.3.

Under the same assumptions as in Proposition 2.2, we have that

$$_{\mathrm{}}^{\mathrm{}}P_n(x)x^{n_i}e^{(V(x)a_ix)}𝑑x0$$
(2.16)
for $`i=1,\mathrm{},p`$.

###### Proof.

Let $`P_{n+1}`$ be the average characteristic polynomial of an ensemble of $`(n+1)\times (n+1)`$ Hermitian random matrices whose external source has the same eigenvalues as $`A`$ plus an additional eigenvalue $`a_i`$. Then by part (c) of Proposition 2.2 we have that $`P_{n+1}`$ is the unique monic polynomial that satisfies the relations (2.3) with $`n_i`$ replaced by $`n_i+1`$. If $`_{\mathrm{}}^{\mathrm{}}P_n(x)x^{n_i}e^{(V(x)a_ix)}𝑑x`$ would vanish, then $`P_{n+1}+P_n`$ would satisfy these relations as well, which would contradict the uniqueness of $`P_{n+1}`$. ∎

Remark: The relations (2.3) can be viewed as multiple orthogonality conditions for the polynomial $`P_n`$. There are $`p`$ weights $`e^{(V(x)a_jx)}`$, $`j=1,\mathrm{},p`$, and for each weight there are a number of orthogonality conditions, so that the total number of them is $`n`$. This point of view is especially useful in case $`A`$ has only a small number of distinct eigenvalues. We will come back to this in Section 5 when we are considering the case of two distinct eigenvalues in detail.

There is a considerable literature on multiple orthogonal polynomials (also called Hermite-Padé polynomials), see e.g.  and the references therein.

## 3 Determinantal form of joint probability density function

As in Proposition 2.2, we assume that $`A`$ is a fixed Hermitian matrix whose eigenvalues $`a_1,\mathrm{},a_p`$ have respective multiplicities $`n_1,\mathrm{},n_p`$, so that $`_{i=1}^pn_i=n`$. We let $`\mathrm{\Sigma }_n`$ be the collection of functions

$$\mathrm{\Sigma }_n:=\{x^je^{a_ix}i=1,\mathrm{},p,j=0,\mathrm{},n_i1\}.$$
(3.1)
We start with a lemma.

###### Lemma 3.1.

There exists a unique function $`Q_{n1}`$ in the linear span of $`\mathrm{\Sigma }_n`$ such that

$$_{\mathrm{}}^{\mathrm{}}x^jQ_{n1}(x)e^{V(x)}𝑑x=0\text{ for }j=0,\mathrm{},n2,$$
(3.2)
and

$$_{\mathrm{}}^{\mathrm{}}x^{n1}Q_{n1}(x)e^{V(x)}𝑑x=1.$$
(3.3)
###### Proof.

The conditions (3.2) and (3.3) give us $`n`$ linear equations for the $`n`$ coefficients of $`Q_{n1}`$ with respect to the basis $`\mathrm{\Sigma }_n`$ with coefficient matrix

$$\left(\begin{array}{cccc}m_{10}& m_{20}& \mathrm{}& m_{n0}\\ m_{11}& m_{21}& \mathrm{}& m_{n1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ m_{1,n1}& m_{2,n1}& \mathrm{}& m_{n,n1}\end{array}\right)$$

where $`m_{jk}`$ is as in part (b) of Proposition 2.2. This matrix is non-singular by (2.13), so that the linear equations have a unique solution, and therefore $`Q_{n1}`$ exists and is unique. ∎

For the rest of this section, we choose some ordering of the eigenvalues of $`A`$ taking into account the multiplicities, say

$$\alpha _1,\alpha _2,\mathrm{},\alpha _n,$$
(3.4)
so that each $`a_i`$ appears exactly $`n_i`$ times among the $`\alpha `$’s. For each $`k=0,1,\mathrm{},n`$, we can construct $`P_k`$ as in the previous section, but based on $`\alpha _1,\mathrm{},\alpha _k`$. Thus $`P_k`$ is a monic polynomial of degree $`k`$ such that

$$_{\mathrm{}}^{\mathrm{}}P_k(x)x^je^{(V(x)a_ix)}𝑑x=0,i=1,\mathrm{},p,j=0,\mathrm{},k_i1,$$
(3.5)
where $`k_i`$ is the number of times that $`a_i`$ appears among $`\alpha _1,\mathrm{},\alpha _k`$. We also have tht $`P_k`$ is the average characteristic polynomial of the ensemble of $`k\times k`$ Hermitian matrices with external source having eigenvalues $`a_i`$ with multiplicity $`k_i`$.

For each $`k=1,\mathrm{},n`$ we also have by Lemma 3.1 a function $`Q_{k1}`$ from the linear span of the functions

$$\mathrm{\Sigma }_k:=\{x^je^{a_ix}i=1,\mathrm{},p,j=0,\mathrm{},k_i1\}.$$
(3.6)
such that

$$_{\mathrm{}}^{\mathrm{}}x^iQ_{k1}(x)e^{V(x)}𝑑x=0\text{ for }i=0,\mathrm{},k2$$
(3.7)
and

$$_{\mathrm{}}^{\mathrm{}}x^{k1}Q_{k1}(x)e^{V(x)}𝑑x=1.$$
(3.8)
It follows from (3.5), (3.7), and (3.8) that the $`P`$’s and $`Q`$’s are a biorthogonal system in the sense that

$$_{\mathrm{}}^{\mathrm{}}P_j(x)Q_k(x)e^{V(x)}𝑑x=\delta _{jk},\text{ for }j,k=0,\mathrm{},n1.$$
(3.9)
This property explains why we used $`Q_{k1}`$ for the function that satisfies (3.7) and (3.8) (and not $`Q_k`$).

We now introduce the kernel $`K_n`$.

Definition: With the polynomials $`P_k`$ and the functions $`Q_k`$ introduced above, we define

$$K_n(x,y)=e^{\frac{1}{2}(V(x)+V(y))}\underset{k=0}{\overset{n1}{}}P_k(x)Q_k(y).$$
(3.10)
Note that the $`P`$’s and the $`Q`$’s depend on the specific ordering (3.4) that we choose for the eigenvalues of $`A`$. However, it will turn out that $`K_n`$ does not depend on this ordering.

Because of the biorthogonality property (3.9) it is easy to see from the definition (3.10) that we have

$$_{\mathrm{}}^{\mathrm{}}K_n(x,x)𝑑x=n$$
(3.11)
and the reproducing kernel property

$$_{\mathrm{}}^{\mathrm{}}K_n(x,y)K_n(y,z)𝑑y=K_n(x,z).$$
(3.12)
The following is the main theorem of this paper.

###### Theorem 3.2.

The joint probability density function on eigenvalues has the determinantal form

$$\frac{1}{n!}det(K_n(\lambda _j,\lambda _k))_{1j,kn}$$
(3.13)
The $`m`$-point correlation function has the form

$$R_m(\lambda _1,\mathrm{},\lambda _m)=det(K_n(\lambda _j,\lambda _k))_{1j,km}$$
(3.14)
###### Proof.

Any joint probability density function of the form (3.13) with a kernel $`K_n`$ satisfying (3.11) and (3.12) leads to $`m`$-point correlation functions of the form (3.14). So it suffices to prove that (3.13) is the joint probability density function of the eigenvalues.

For each $`j`$, we define

$$w_j(x)=x^{d_j1}e^{a_ix}$$
(3.15)
if $`a_i=\alpha _j`$ and $`a_i`$ appears $`d_j`$ times in the sequence $`\alpha _1,\mathrm{},\alpha _j`$. Note that the functions (3.15) differ from the functions $`w_j`$ used in Proposition 2.2 in two respects. First there is an extra factor $`e^{V(x)}`$ in Proposition 2.2, and second we used a specific ordering of the eigenvalues of $`A`$ in Proposition 2.2 (which only amounts to a renumbering).

A similar calculation as that leading to (2.5) in the proof of Proposition 2.1 shows that the joint probability density of eigenvalues is proportional to

$$\underset{j=1}{\overset{n}{}}e^{V(\lambda _j)}det(w_i(\lambda _j))_{1i,jn}\mathrm{\Delta }(\lambda ).$$

Since $`Q_{i1}`$ is a linear combination of $`w_1,\mathrm{},w_i`$ we can take appropriate row combinations to find that

$$det(w_i(\lambda _j))_{1i,jn}det(Q_{i1}(\lambda _j))_{1i,jn}.$$

We write $`\mathrm{\Delta }(\lambda )`$ as a Vandermonde determinant which we similarly rewrite as

$$\mathrm{\Delta }(\lambda )=det(P_{i1}(\lambda _k))_{1i,kn}.$$

Thus the joint probability density of eigenvalues is proportional to

$$det\left(e^{\frac{1}{2}V(\lambda _j)}Q_{i1}(\lambda _j)\right)_{1i,jn}det\left(e^{\frac{1}{2}V(\lambda _k)}P_{i1}(\lambda _k)\right)_{1i,kn}.$$

Taking the transpose of the matrix in the first determinant, and then using the multiplicative property of determinants, we find that the joint probability density is equal to

$$cdet(K_n(\lambda _j,\lambda _k))_{1j,kn}$$

for some constant $`c`$, which should be such that the integral with respect to $`d\lambda _1\mathrm{}d\lambda _n`$ is $`1`$. Because of the properties (3.11) and (3.12) this is so for $`c=\frac{1}{n!}`$ and the theorem is proved. ∎

Remark: Renumbering the eigenvalues $`a_1,a_2,\mathrm{}a_n`$ leads to the same kernel $`K_n`$ but to different $`P_k`$ and $`Q_k`$.

## 4 Special form of the kernel in case of two eigenvalues

In this section we assume we have only two distinct eigenvalues $`a_1`$ and $`a_2`$ with multiplicities $`n_1`$ and $`n_2`$, respectively. We order the eigenvalues $`\alpha _1,\alpha _2,\mathrm{},\alpha _n`$ in some arbitrary way ($`a_j`$ appear $`n_j`$ times in the sequence), but for convenience we assume that

$$\alpha _{n1}=a_1\alpha _n=a_2.$$
(4.1)
We also put

$$\alpha _{n+1}=a_1\alpha _{n+2}=a_2.$$
(4.2)
As in the preceding section, we have polynomials $`P_k`$ and functions $`Q_k`$ for every $`k=0,\mathrm{},n1`$ such that

$$K_n(x,y)=e^{\frac{1}{2}(V(x)+V(y))}\underset{k=0}{\overset{n1}{}}P_k(x)Q_k(y).$$

It is our aim in this section to simplify this expression. The formula we will find is an analogue of the well-known Christoffel-Darboux formula for orthogonal polynomials.

To present the formula we are going to use multi-index notation. For non-negative integers $`k_1`$ and $`k_2`$, we use $`P_{k_1,k_2}`$ to denote the monic polynomial of degree $`k_1+k_2`$ having $`k_j`$ orthogonality relations with respect to the weight

$$w_j(x)=e^{(V(x)a_jx)},j=1,2.$$

Thus

$$_{\mathrm{}}^{\mathrm{}}P_{k_1,k_2}(x)x^iw_j(x)𝑑x=0,i=0,\mathrm{},k_j1,j=1,2.$$
(4.3)
The polynomial $`P_{k_1,k_2}`$ is called a multiple orthogonal polynomial of type II, see e.g. . We also define

$$Q_{k_1,k_2}(x)=A_{k_1,k_2}(x)e^{a_1x}+B_{k_1,k_2}(x)e^{a_2x}$$
(4.4)
where the degree of $`A_{k_1,k_2}`$ is $`k_11`$, the degree of $`B_{k_1,k_2}`$ is $`k_21`$ and

$$_{\mathrm{}}^{\mathrm{}}x^jQ_{k_1,k_2}(x)e^{V(x)}𝑑x=\{\begin{array}{cc}0,& j=0,\mathrm{},k_1+k_22,\hfill \\ 1& j=k_1+k_21.\hfill \end{array}$$
(4.5)
The polynomials $`A_{k_1,k_2}`$ and $`B_{k_1,k_2}`$ are called multiple orthgonal polynomials of type I, see . For each pair $`(k_1,k_2)`$ of non-negative integers, the polynomials $`P_{k_1,k_2}`$, $`A_{k_1,k_2}`$, and $`B_{k_1,k_2}`$ exist and are uniquely defined by their degree requirements and the relations (4.3), (4.4), and (4.5).

We can express $`P_k`$ and $`Q_k`$ in this new notation as

$$P_k=P_{k_1,k_2},\text{and}Q_{k1}=Q_{k_1,k_2}$$

provided $`a_j`$ appears $`k_j`$ times among the numbers $`\alpha _1,\mathrm{},\alpha _k`$ (for $`j=1,2`$). In particular, we have because of our assumptions (4.1) and (4.2)

$$P_n=P_{n_1,n_2},P_{n1}=P_{n_1,n_21},P_{n2}=P_{n_11,n_21}.$$
(4.6)
and

$$Q_{n1}=Q_{n_1,n_2},Q_n=Q_{n_1+1,n_2},Q_{n+1}=Q_{n_1+1,n_2+1}.$$
(4.7)
We also need the numbers

$$h_{k_1,k_2}^{(j)}=_{\mathrm{}}^{\mathrm{}}P_{k_1,k_2}(x)x^{k_j}w_j(x)𝑑x,j=1,2,$$
(4.8)
which are non-zero, cf. (2.16). For later use we note that

$`1`$ $`=`$ $`{\displaystyle P_{k_1,k_2}(x)Q_{k_1+1,k_2}(x)e^{V(x)}𝑑x}`$
$`=`$ $`{\displaystyle P_{k_1,k_2}(x)(A_{k_1+1,k_2}(x)w_1(x)+B_{k_1+1,k_2}(x)w_2(x))𝑑x}`$
$`=`$ $`{\displaystyle P_{k_1,k_2}(x)A_{k_1+1,k_2}(x)w_1(x)𝑑x}`$
$`=`$ $`\left(\text{leading coefficient of }A_{k_1+1,k_2}\right)\times h_{k_1,k_2}^{(1)}`$
so that

$$\text{leading coefficient of }A_{k_1+1,k_2}=\frac{1}{h_{k_1,k_2}^{(1)}}.$$
(4.9)
Similarly,

$$\text{leading coefficient of }B_{k_1,k_2+1}=\frac{1}{h_{k_1,k_2}^{(2)}}.$$
(4.10)
It also follows from (4.9) and (4.10) that $`h_{k_1,k_2}^{(j)}0`$ for $`j=1,2`$.

Then we can state the following theorem.

###### Theorem 4.1.

With the notation introduced above, we have

$`(xy)e^{\frac{1}{2}(V(x)+V(y))}K_n(x,y)`$ $`=`$ $`P_{n_1,n_2}(x)Q_{n_1,n_2}(y)`$ (4.11)
$`{\displaystyle \frac{h_{n_1,n_2}^{(1)}}{h_{n_11,n_2}^{(1)}}}P_{n_11,n_2}(x)Q_{n_1+1,n_2}(y)`$
$`{\displaystyle \frac{h_{n_1,n_2}^{(2)}}{h_{n_1,n_21}^{(2)}}}P_{n_1,n_21}(x)Q_{n_1,n_2+1}(y)`$
The proof of the theorem needs some preparation. We start working again with the $`P_k`$’s and $`Q_j`$’s (single index) as before. For each $`j`$ and $`k`$, we put

$$c_{jk}=xP_k(x)Q_j(x)e^{V(x)}𝑑x.$$

The coefficients $`c_{jk}`$ appear in the expansion

$$xP_k(x)=\underset{j=0}{\overset{k+1}{}}c_{jk}P_j(x),$$
(4.12)
since by the biorthogonality relation we have indeed

$$c_{jk}=xP_k(x)Q_j(x)e^{V(x)}𝑑x.$$

Similarly, we have for $`j=0,\mathrm{},n1`$,

$$xQ_j(x)=\underset{k=0}{\overset{n+1}{}}c_{jk}Q_k(x).$$
(4.13)
Note that by adding the two values $`\alpha _{n+1}`$ and $`\alpha _{n+2}`$ as we did in (4.2), we have this expansion for every $`jn1`$.

###### Lemma 4.2.

1. If $`jk+2`$ then $`c_{jk}=0`$.
2. If $`kj+3`$ and if both $`a_1`$ and $`a_2`$ appear at least once among $`\alpha _{j+2},\alpha _{j+3},\mathrm{},\alpha _k`$, then $`c_{jk}=0`$.

###### Proof.

(a) We have that

$$P(x)Q_j(x)e^{V(x)}𝑑x=0$$

for every polynomial $`P`$ of degree $`j1`$. Since $`xP_k`$ is a polynomial of degree $`k+1`$, it follows that $`c_{jk}=0`$ if $`k+1j1`$. This proves part (a).

(b) Let $`k`$ and $`j`$ be such that the conditions of part (b) are satisfied. Suppose that $`a_1`$ appears $`k_1`$ times among $`\alpha _1,\mathrm{},\alpha _k`$, and $`j_1`$ times among $`\alpha _1,\mathrm{},\alpha _{j+1}`$. We put $`k_2=kk_1`$ and $`j_2=j+1j_1`$. It follows from the assumptions that $`j_1<k_1`$ and $`j_2<k_2`$. Then $`Q_j(x)=Q_{j_1,j_2}(x)=A_{j_1,j_2}(x)e^{a_1x}+B_{j_1,j_2}(x)e^{a_2x}`$ where $`A_{j_1,j_2}`$ has degree $`j_11`$ and $`B_{j_1,j_2}`$ has degree $`j_21`$. It follows that

$$xQ_j(x)=xA_{j_1,j_2}(x)e^{a_1x}+xB_{j_1,j_2}(x)e^{a_2x}$$

and $`xA_{j_1,j_2}(x)`$ has degree $`j_1k_11`$ and $`xB_{j_1,j_2}(x)`$ has degree $`j_2k_21`$. Thus, by the multiple orthogonality property of $`P_k=P_{k_1,k_2}`$, we have

$$P_k(x)xQ_j(x)e^{V(x)}𝑑x=0.$$

This proves part (b). ∎

We also need the following relations between near-by $`P`$’s and $`Q`$’s.

###### Lemma 4.3.

We have

$`P_{n_11,n_21}`$ $`=`$ $`{\displaystyle \frac{h_{n_11,n_21}^{(1)}}{h_{n_11,n_2}^{(1)}}}\left(P_{n_11,n_2}P_{n_1,n_21}\right)`$ (4.14)
$`=`$ $`{\displaystyle \frac{h_{n_11,n_21}^{(2)}}{h_{n_1,n_21}^{(2)}}}\left(P_{n_11,n_2}P_{n_1,n_21}\right)`$
and

$`Q_{n_1+1,n_2+1}`$ $`=`$ $`{\displaystyle \frac{h_{n_1,n_2}^{(1)}}{h_{n_1,n_2+1}^{(1)}}}\left(Q_{n_1,n_2+1}Q_{n_1+1,n_2}\right)`$ (4.15)
$`=`$ $`{\displaystyle \frac{h_{n_1,n_2}^{(2)}}{h_{n_1+1,n_2}^{(2)}}}\left(Q_{n_1,n_2+1}Q_{n_1+1,n_2}\right)`$
###### Proof.

Since $`P_{n_11,n_2}`$ and $`P_{n_1,n_21}`$ are both monic polynomials of degree $`n`$, their difference is a polynomial of degree $`n1`$. Since this difference has $`n_j1`$ orthogonality conditions with respect to $`w_j`$ for $`j=1,2`$, it must be a multiple of $`P_{n_11,n_21}`$. Thus

$$P_{n_11,n_2}P_{n_1,n_21}=\gamma P_{n_11,n_21}$$

for some $`\gamma `$. Integrating this equation with respect to $`x^{n_11}w_1(x)`$ and $`x^{n_21}w_2(x)`$, we get $`h_{n_11,n_2}^{(1)}=\gamma h_{n_11,n_21}^{(1)}`$, and $`h_{n_1,n_21}^{(2)}=\gamma h_{n_11,n_21}^{(2)}`$, respectively. This gives (4.14).

Next we note that we have

$$x^j(Q_{n_1,n_2+1}Q_{n_1+1,n_2})e^{V(x)}𝑑x=00=0,j=0,\mathrm{},n_1+n_21$$

and also

$$x^{n_1+n_2}(Q_{n_1,n_2+1}Q_{n_1+1,n_2})e^{V(x)}𝑑x=11=0.$$

Since $`Q_{n_1,n_2+1}(x)Q_{n_1+1,n_2}(x)=A(x)e^{a_1x}+B(x)e^{a_2x}`$ where $`A`$ has degree $`n_1`$ and $`B`$ has degree $`n_2`$, it follows that $`Q_{n_1,n_2+1}Q_{n_1+1,n_2}`$ is a multiple of $`Q_{n_1+1,n_2+1}`$, say

$$Q_{n_1,n_2+1}Q_{n_1+1,n_2}=\beta Q_{n_1+1,n_2+1}.$$

This means for the $`A`$-polynomials that

$$A_{n_1,n_2+1}A_{n_1+1,n_2}=\beta A_{n_1+1,n_2+1}$$

and looking at the leading coefficient (= coefficient of $`x^{n_1}`$) we get

$$\beta =\frac{\text{leading coefficient of }A_{n_1+1,n_2}}{\text{leading coefficient of }A_{n_1+1,n_2+1}}=\frac{h_{n_1,n_2+1}^{(1)}}{h_{n_1,n_2}^{(1)}},$$

where we used (4.9). We also get by considering the $`B`$-polynomials that

$$\beta =\frac{\text{leading coefficient of }B_{n_1,n_2+1}}{\text{leading coefficient of }B_{n_1+1,n_2+1}}=\frac{h_{n_1+1,n_2}^{(2)}}{h_{n_1,n_2}^{(2)}}$$

because of (4.10). This proves (4.15). ∎

Now we are ready for the proof of Theorem 4.1.

###### Proof.

We note that $`(xy)e^{\frac{1}{2}(V(x)+V(y))}K_n(x,y)`$ has a telescoping character. Indeed we have by (4.12) and (4.13),

$`(xy){\displaystyle \underset{k=0}{\overset{n1}{}}}P_k(x)Q_k(y)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{n1}{}}}xP_k(x)Q_k(y){\displaystyle \underset{j=0}{\overset{n1}{}}}yP_j(x)Q_j(y)`$
$`=`$ $`{\displaystyle \underset{k=0}{\overset{n1}{}}}{\displaystyle \underset{j=0}{\overset{k+1}{}}}c_{jk}P_j(x)Q_k(y){\displaystyle \underset{j=0}{\overset{n1}{}}}{\displaystyle \underset{k=0}{\overset{n+1}{}}}c_{jk}P_j(x)Q_k(y)`$
$`=`$ $`c_{n,n1}P_n(x)Q_{n1}(y){\displaystyle \underset{j=0}{\overset{n1}{}}}c_{jn}P_j(x)Q_n(y){\displaystyle \underset{j=0}{\overset{n1}{}}}c_{j,n+1}P_j(x)Q_{n+1}(y).`$
Now observe that $`c_{n,n1}=1`$, and that $`c_{jn}=0`$ for $`j=0,\mathrm{},n3`$ and $`c_{j,n+1}=0`$ for $`j=0,\mathrm{},n2`$, which follows from Lemma 4.2. Thus

$`(xy)e^{\frac{1}{2}(V(x)+V(y))}K_n(x,y)`$ $`=`$ $`P_n(x)Q_{n1}(y)`$ (4.16)
$`c_{n2,n}P_{n2}(x)Q_n(y)`$
$`c_{n1,n}P_{n1}(x)Q_n(y)`$
$`c_{n1,n+1}P_{n1}(x)Q_{n+1}(y).`$
In formula (4.16) we have reduced the $`n`$-term expression to four terms, which is already quite nice. However, we want to reduce to three terms only. Changing back to multi-index notation and using (4.6) and (4.7), we see that (4.16) leads to

$`(xy)e^{\frac{1}{2}(V(x)+V(y))}K_n(x,y)`$ $`=`$ $`P_{n_1,n_2}(x)Q_{n_1,n_2}(y)`$ (4.17)
$`c_{n2,n}P_{n_11,n_21}(x)Q_{n_1+1,n_2}(y)`$
$`c_{n1,n}P_{n_1,n_21}(x)Q_{n_1+1,n_2}(y)`$
$`c_{n1,n+1}P_{n_1,n_21}(x)Q_{n_1+1,n_2+1}(y)`$
Comparing (4.17) and (4.11) we see that we need to get rid of $`P_{n_11,n_21}`$ and $`Q_{n_1+1,n_2+1}`$. This can be done using the following relations between near-by $`P`$’s and $`Q`$’s.

Our next task is to express the recurrence coefficients $`c_{n2,n}`$, $`c_{n1,n}`$, and $`c_{n1,n+1}`$ that appear in (4.17) in terms of the $`h`$-numbers. This is rather straightforward from the definition. Indeed, we have

$`c_{n2,n}`$ $`=`$ $`{\displaystyle xP_{n_1,n_2}(x)Q_{n_1,n_21}(x)e^{V(x)}𝑑x}`$ (4.18)
$`=`$ $`{\displaystyle P_{n_1,n_2}(x)\left(xA_{n_1,n_21}(x)w_1(x)+xB_{n_1,n_21}(x)w_2(x)\right)𝑑x}`$
$`=`$ $`{\displaystyle P_{n_1,n_2}(x)xA_{n_1,n_21}(x)w_1(x)𝑑x}`$
$`=`$ $`\left(\text{leading coefficient of }A_{n_1,n_21}\right)\times h_{n_1,n_2}^{(1)}`$
$`=`$ $`{\displaystyle \frac{h_{n_1,n_2}^{(1)}}{h_{n_11,n_21}^{(1)}}}`$
where we used (4.9), and similarly,

$`c_{n1,n}={\displaystyle \frac{h_{n_1,n_2}^{(1)}}{h_{n_11,n_2}^{(1)}}}+{\displaystyle \frac{h_{n_1,n_2}^{(2)}}{h_{n_1,n_21}^{(2)}}},`$ (4.19)
and

$`c_{n1,n+1}={\displaystyle \frac{h_{n_1+1,n_2}^{(2)}}{h_{n_1,n_21}^{(2)}}}.`$ (4.20)
Now we plug all formulas (4.14), (4.15), (4.18), (4.19) and (4.20) into (4.17). Straightforward calculations then lead to (4.11). ∎

## 5 Riemann-Hilbert problem

We use the notation of Section 4.

The Christoffel-Darboux formula (4.11) can be expressed in terms of the solution of a Riemann-Hilbert problem that was given by Van Assche, Geronimo, and Kuijlaars  to characterize the multiple orthogonal polynomials, and which generalizes the Riemann-Hilbert problem for orthogonal polynomials due to Fokas, Its, and Kitaev . The Rieman-Hilbert problem is to find $`Y:^{3\times 3}`$ such that

* $`Y`$ is analytic on $``$,
* for $`x`$, we have

$$Y_+(x)=Y_{}(x)\left(\begin{array}{ccc}1& w_1(x)& w_2(x)\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$$
(5.1)
where $`Y_+(x)`$ ($`Y_{}(x)`$) denotes the limit of $`Y(z)`$ as $`zx`$ from the upper (lower) half-plane,
* as $`z\mathrm{}`$, we have

$$Y(z)=\left(I+O\left(\frac{1}{z}\right)\right)\left(\begin{array}{ccc}z^n& 0& 0\\ 0& z^{n_1}& 0\\ 0& 0& z^{n_2}\end{array}\right)$$
(5.2)
where $`I`$ denotes the $`3\times 3`$ identity matrix.

In  it was shown showed that there is a unique solution

$$Y=\left(\begin{array}{ccc}P_{n_1,n_2}& C(P_{n_1,n_2}w_1)& C(P_{n_1,n_2}w_2)\\ c_1P_{n_11,n_2}& c_1C(P_{n_11,n_2}w_1)& c_1C(P_{n_11,n_2}w_2)\\ c_2P_{n_1,n_21}& c_2C(P_{n_1,n_21}w_1)& c_2C(P_{n_1,n_21}w_2)\end{array}\right)$$
(5.3)
with constants

$$c_1=2\pi i\left(h_{n_11,n_2}^{(1)}\right)^1,\text{and}c_2=2\pi i\left(h_{n_1,n_21}^{(2)}\right)^1,$$
(5.4)
and where $`Cf`$ denotes the Cauchy transform of $`f`$, i.e.,

$$Cf(z)=\frac{1}{2\pi i}_{}\frac{f(s)}{sz}𝑑s.$$

The multiple orthogonal polynomials of type I $`A_{n_1,n_2}`$, $`B_{n_1,n_2}`$ have a Riemann-Hilbert characterization as well. We seek $`X:^{3\times 3}`$ such that

* $`X`$ is analytic on $``$,
* for $`x`$, we have

$$X_+(x)=X_{}(x)\left(\begin{array}{ccc}1& 0& 0\\ w_1(x)& 1& 0\\ w_2(x)& 0& 1\end{array}\right)$$
(5.5)
* as $`z\mathrm{}`$, we have

$$X(z)=\left(I+O\left(\frac{1}{z}\right)\right)\left(\begin{array}{ccc}z^n& 0& 0\\ 0& z^{n_1}& 0\\ 0& 0& z^{n_2}\end{array}\right).$$
(5.6)

The solution to this Riemann-Hilbert problem  is

$$X=\left(\begin{array}{ccc}2\pi iC(A_{n_1,n_2}w_1+B_{n_1,n_2}w_2)& 2\pi iA_{n_1,n_2}& 2\pi iB_{n_1,n_2}\\ k_1C(A_{n_1+1,n_2}w_1+B_{n_1+1,n_2}w_2)& k_1A_{n_1+1,n_2}& k_1B_{n_1+1,n_2}\\ k_2C(A_{n_1,n_2+1}w_2+B_{n_1,n_2+1}w_2)& k_2A_{n_1,n_2+1}& k_2B_{n_1,n_2+1}\end{array}\right)$$
(5.7)
where

$$k_1=\frac{1}{\text{leading coefficient of }A_{n_1+1,n_2}}=h_{n_1,n_2}^{(1)},$$
(5.8)
and

$$k_2=\frac{1}{\text{leading coefficient of }B_{n_1,n_2+1}}=h_{n_1,n_2}^{(2)}.$$
(5.9)
It is easy to see that

$$X=Y^t\text{(inverse transpose)}$$
(5.10)
see also .

Now we form the product $`Y^1(y)Y(x)=X^t(y)Y(x)`$ and we compute the $`21`$-entry using (5.3), (5.4), (5.7), (5.8), (5.9),

$`[Y^1(y)Y(x)]_{21}=\left(\begin{array}{ccc}2\pi iA_{n_1,n_2}(y)& k_1A_{n_1+1,n_2}(y)& k_2A_{n_1,n_2+1}(y)\end{array}\right)\left(\begin{array}{c}P_{n_1,n_2}(x)\\ c_1P_{n_11,n_2}(x)\\ c_2P_{n_1,n_21}(x)\end{array}\right)`$
$`=2\pi i\left(P_{n_1,n_2}(x)A_{n_1,n_2}(y){\displaystyle \frac{h_{n_1,n_2}^{(1)}}{h_{n_11,n_2}^{(1)}}}P_{n_11,n_2}(x)A_{n_1+1,n_2}(y){\displaystyle \frac{h_{n_1,n_2}^{(2)}}{h_{n_1,n_21}^{(2)}}}P_{n_1,n_21}(x)A_{n_1,n_2+1}(y)\right).`$
We get a similar expression for the $`31`$-entry $`[Y^1(y)Y(x)]_{31}`$, but with the $`B`$-polynomials instead of the $`A`$-polynomials. Then it follows that we can rewrite the Christoffel-Darboux formula (4.11) as

$$K_n(x,y)=e^{\frac{1}{2}(V(x)+V(y))}\frac{e^{a_1y}[Y^1(y)Y(x)]_{21}+e^{a_2y}[Y^1(y)Y(x)]_{31}}{2\pi i(xy)}$$
(5.11)
which is a compact form for the kernel in terms of the solution of the Riemann-Hilbert problem. We expect that the Riemann-Hilbert problem for $`Y`$ is tractable to asymptotic analysis using the methods of  or .

## Acknowledgments

The first author was supported in part by NSF Grant DMS-9970625. The second author was supported by FWO research grant G.0176.02, by INTAS project 00-272, and by the Ministry of Science and Technology (MCYT) of Spain, project code BFM2001-3878-C02-02. He is grateful to Peter Dragnev (IUPUI-Fort Wayne) for making a visit to Indiana possible.

P.M. Bleher: Department of Mathematical Sciences, Indiana University-Purdue University Indianapolis, 402 N. Blackford Street, Indianapolis, IN 46202, U.S.A.

E-mail address: bleher@math.iupui.edu

A.B.J. Kuijlaars: Department of Mathematics, Katholieke Universiteit Leuven, Celestijnenlaan 200 B, 3001 Leuven, Belgium

E-mail address: arno@wis.kuleuven.ac.be