Clustering of Extremely Red Objects in Elais-N1 from the UKIDSS DXS with optical photometry from Pan-STARRS 1 and Subaru

The UKIRT Infrared Deep Sky Survey (UKIDSS) began in 2005 and consists of 5 sub-surveys (Lawrence et al. 2007). The Wide Field Camera (WFCAM, Casali et al. 2007) mounted on the UK Infrared Telescope (UKIRT) has been used to obtain UKIDSS images. The Deep eXtragalactic Survey (DXS) is a deep, wide survey mapping 35 deg² with a 5$\sigma$ point-source sensitivity of $J_{AB}\sim 23.2$ and $K_{AB}\sim 22.7$ as one of the sub-surveys. It covers 4 different fields and aims to detect photometric samples of $z\sim 1-2$ galaxies.

WFCAM is composed of four Rockwell Hawaii-II 2K$\times$2K array detectors (Casali et al. 2007). The sky coverage of each detector is 13.7$\times$13.7 arcmin² with 0.4 arcsec/pixel. In order to avoid an undersampled point spread function caused by the relatively large pixel scale, a microstepping technique is applied, i.e. 0.2 arcsec/pixel for the final science image. Since there are gaps between each detector, 4 exposures are necessary to generate a contiguous image covering 0.8 deg².

In this study we deal with the Elais-N1 field centred on $\alpha=$ 16^h 11^m 09.7^s, $\delta=$ +55^d 0^m 47.0^s (J2000). The datasets from UKIDSS data releases 7 and 8 were used for this work¹¹1http://surveys.roe.ac.uk/wsa/. In these releases, $K$-band data covers the whole region but $J$-band coverage is currently only $\sim$56 per cent. The typical seeing is $\sim 0.9^{\prime\prime}$ at $J$ and $\sim 0.8^{\prime\prime}$ at $K$. Although the $K$-band dataset has mapped 6.5 deg² after masking unreliable regions, the actual area for this work depends on the optical datasets (see next subsection).

There are some known issues about the WFCAM images and database catalogues, such as cross-talk and non-optimal galaxy photometry (Dye et al. 2006). To avoid these, we therefore created our own photometric catalogues from the stacked images. Full details are described in Kim et al. (2011a) so we simply summarise the main steps in this paper. Stacked images from the UKIDSS standard pipeline were combined into individual images for each pointing using the Swarp software package (Bertin et al. 2002). Then astronomical objects were extracted with SExtractor (Bertin & Arnouts 1996) and a 2-arcsec aperture magnitude for colour and the AUTO magnitude for total magnitude were also measured. Finally spurious objects such as cross-talk images, diffraction spikes and duplicated objects in overlapping regions were removed. We found 670 214 objects and determined the completeness from an artificial star test to be $>90$ per cent at the DXS magnitude goals of $J_{AB}=23.2$ and $K_{AB}=22.7$. The Vega-AB offsets for these bands are: $-0.938$ mag. for $J$ and $-1.900$ mag. for $K$.

In order to identify EROs we require deep optical imaging. We use two different datasets for this paper, one that best matches the DXS area (Pan-STARRS) and one that is deeper but over a smaller area (Subaru). The combination of the two allows us to determine the clustering properties of EROs as a function of depth and area.

In this study EROs are selected using the $i-K$ colour from DXS/PS1 and DXS/Subaru. Firstly, Galactic stars must be removed to avoid contamination. In the case of DXS/PS1 various schemes were applied. Bright stars ($K_{AB}<16.3$) were removed using the magnitude difference between $K$-band aperture and total magnitudes. Then stellar sequences in $(i-K)$ vs. $(g-i)$ and $(r-[3.6])$ vs. $(r-i)$ colour-colour diagrams were extracted. These criteria are $(i-K)_{AB}<0.76(g-i)_{AB}-0.85$ and $(r-[3.6])_{AB}<2.29(r-i)_{AB}-0.66$. Stars in the DXS/Subaru data were selected by comparing with those in the DXS/PS1. The faintest candidate stars in DXS/Subaru were not removed, but this does not affect our analysis because very few of these objects meet the ERO colour cut. We note that the fraction of faint stars ($i_{AB}>25$) selected as EROs is less than 1 per cent, based on the UKIDSS Ultra Deep Survey DR3 catalogue (Simpson et al. 2006).

From the catalogues with Galactic stars removed, colour criteria were applied to select EROs. For DXS/PS1 $(i-K)_{AB}>2.45,2.95$, $i_{AB}<25$ and $K_{AB}<22.7$ limits were applied to satisfy the classical colour cut ($I-K)_{vega}>4$) for EROs and match the observed magnitudes. For the DXS/Subaru dataset, $(i-K)_{AB}>2.55,3.05$, $i_{AB}<25.5$ and $K_{AB}<22.7$ limits were used. The magnitude limits applied for each band are set to ensure we detect objects uniformly across the whole area (see the previous section). The difference in the applied colour cuts between the two catalogues arises purely from the magnitude difference between PS1 $i$-band and Subaru $i$-band, derived by comparing common objects in both catalogues. Hereafter we quote $(i-K)_{AB}=2.45$ and $2.95$ for EROs from both datasets instead of $2.55$ and $3.05$. Using these criteria we selected 17 250 and 23 916 EROs with $(i-K)_{AB}>2.45$ and 5 039 and 7 959 with $(i-K)_{AB}>2.95$ from the DXS/PS1 and DXS/Subaru datasets respectively.

Fig. 1 shows the number counts for all galaxies and for just EROs. The circles indicate the number counts of all galaxies, and triangles and squares are for EROs with bluer ($(i-K)_{AB}=2.55$ for DXS/Subaru and $2.45$ for DXS/PS1) and redder ($(i-K)_{AB}=3.05$ for DXS/PS1 and $2.95$ for DXS/PS1) cuts from this work, respectively. The number counts of galaxies in Lane et al. (2007, dashed line) and Kong et al. (2006, upper dotted line), and EROs in Kong et al. (2006, lower dotted lines) are also displayed. The counts from this work are consistent with previous studies. However the number counts decrease at fainter magnitudes ($K_{AB}>21$) due to the depth of the optical datasets. The shallow depth of the optical datasets prevents the detection of the reddest galaxies. This effect is more significant for the DXS/PS1 sample than DXS/Subaru because of the more restricted depth of the PS1 dataset.

The main purpose of this paper is to compare the properties of haloes which host EROs at different redshifts by measuring their angular clustering. For this purpose the photometric redshifts of EROs were measured using the $g$, $r$, $i$, $z$, $J$, $K$, $3.6$ and $4.5\mu$m photometric data from DXS-PS1-SWIRE. The EAZY photometric redshift code (Brammer, van Dokkum & Coppi 2008) was run to measure the photometric redshift of all objects using the default parameters of the EAZY code. To test the redshift accuracy of this method we used spectroscopic redshifts from Rowan-Robinson et al. (2008), which included those in Berta et al. (2007) and Trichas et al. (2010). The normalised median absolute deviation (NMAD) in $\Delta z/(1+z_{\rm spec})$ was found to be $\sim 0.066$.

However, there are only a small number of spectroscopic redshifts at $z>1$, where most EROs are located. Therefore we also applied the empirical method of Quadri & Williams (2010) to constrain the photometric redshift uncertainty for EROs. This method assumes that close pairs of galaxies should show a significant probability of being located at the same redshift. We counted pairs of EROs having an angular separations $0.04^{\prime}<\theta<0.25^{\prime}$ and those with randomised positions. Then, the difference between the two sets in $\Delta z/(1+z_{\rm mean})$ was used to measure the photometric redshift uncertainty of EROs. This gave a dispersion of $\sigma_{z}\sim 0.059$ which is consistent with the NMAD value from spectroscopic samples. Fig. 2 shows the redshift distributions of EROs in DXS/PS1 using a best-fit photometric redshift. The solid histogram is for $(i-K)_{AB}>2.45$ EROs, and the dashed one is for $(i-K)_{AB}>2.95$ EROs. It is apparent that most EROs are located at $z>1$. However, the $(i-K)_{AB}>2.45$ ERO selection also contains galaxies at $z<1$. The median redshifts are 1.176 and 1.291 for $(i-K)_{AB}>2.45$ and $2.95$ EROs respectively, and displayed in Fig. 2 as arrows. The trend that redder EROs are to be found at higher redshift was also predicted by Gonzalez-Perez et al. (2009).

The angular two-point correlation function is the excess probability of finding a galaxy pair at a given angular separation compared to a random distribution (Peebles 1980). We used the estimator from Landy & Szalay (1993) to estimate the angular two-point correlation function:

where DD is the number of observed ERO pairs with separations $[\theta,\theta+\Delta\theta]$. For this study we used $\Delta\log\theta=$0.15. DR and RR are data-random and random-random pairs in the same interval, respectively. The random catalogue was generated with 30 times more unclustered points than the observed EROs, and had the same angular mask as the EROs. All pair counts were normalised to have the same total numbers.

One of the aims of this study is to investigate the properties of haloes hosting EROs at different redshifts. However, our photometric redshift measurement may not be accurate enough to split samples cleanly into redshift bins. Therefore we used the probability distribution function of our photometric redshifts to measure the angular correlation function of EROs in different redshift bins and estimate the redshift distribution function of EROs for the halo modeling. The details are described in Wake et al. (2011). Briefly, the weight of each ERO is defined as the fractional probability of a particular ERO being within a given redshift interval, using the probability distribution function from the EAZY code. For the angular correlation function this weight was used to count pairs. Also, the weighted probability distribution function of photometric redshifts was used to estimate the redshift distribution. This strategy is similar to that introduced by Myers et al. (2009). In addition, the estimated redshift distribution was also used to measure the number density of EROs in each redshift bin using the same scheme as Ross & Brunner (2009).

The error on the correlation function was estimated using the Jackknife resampling method to compute the deviation of the correlation functions between subfields (for a description see Norberg et al. 2009). We divided the whole area into 25 subfields for DXS/Subaru and 30 subfields for DXS/PS1, then repeated the measurement of the correlation function. From each set of resamplings we can estimate the error using

where $w_{i}$ and DR_i are the correlation function and data-random pairs excluding the $i^{th}$ subfield, and $N$ is the total number of times that the data is resampled. In the Jackknife, $N$ corresponds to the number of subfields into which the dataset is divided. Then the covariance matrix is calculated with

where $\overline{w(\theta_{i})}$ is the mean correlation function of jackknife subsamples in the $i$th bin. The covariance matrix was used to fit the halo model.

The restricted survey area leads to a negative offset of the observed correlation function from the actual one, which is known as the integral constraint (IC, Groth & Peebles 1977). In order to correct for this bias, we applied the same method as Kim et al. (2011a) using the equation in Roche et al. (1999),

Since it is known that the correlation function of EROs is not well described by a single power-law (Gonzalez-Perez et al. 2011; Kim et al. 2011a), the functional form of $w({\theta})=\alpha_{1}\theta^{-\beta_{1}}+\alpha_{2}\exp(-\beta_{2}\theta)$ was used to describe the correlation function. This form is a close match to the angular correlation function measured for ERO samples defined by magnitude and colour cut, and is adopted as the true underlying correlation function to compute the integral constraint in Eq. 4. Then Eq. 4 was used to estimate the integral constraint, which is then added to our estimate of the angular correlation function. Another approach is to use the correlation function obtained from the halo model as the actual correlation function in Eq. 4 (Wake et al. 2011). For the redshift limited samples of EROs, we use the HOD modeled correlation function to calculate the integral constraint with the Eq. 4. We note that the integral constraint ranges from 0.004 to 0.008 with the functional form of the correlation function for magnitude or colour limited EROs, and that brighter or redder EROs have larger integral constraints due to their enhanced clustering strength. Similarly the value obtained from the halo model for redshift limited samples is between 0.004 and 0.007. We are now probing sufficient area that the integral constraint does not have a major effect on our results.

The halo model (see Cooray & Sheth 2002 for a review) is widely used to estimate the mass of the host dark matter haloes of observed galaxies (Blake et al. 2008; Wake et al. 2008a; Ross & Brunner 2009; Zehavi et al. 2011). We apply it to describe the angular correlation functions of EROs and hence to measure the properties of the dark matter haloes which host EROs.

The halo occupation distribution (HOD) describes the mean number of galaxies being hosted by a halo of a given mass ($M$). In the halo model, galaxies are separated into centrals and satellites. The mean number of galaxies, $N(M)$, is the combination of the mean number of central galaxies, $N_{\rm c}(M)$, and satellites, $N_{\rm s}(M)$ (Zheng et al. 2005; Blake et al. 2008; Wake et al. 2008a; Ross & Brunner 2009), i.e.

where the mean number of the central galaxies and satellites are assumed to be described by

and

where $M_{\rm cut}$, $\sigma_{\rm cut}$, $M_{0},M_{1}$ and $\alpha$ are parameters defining the shape of HOD. If $M$ is smaller than $M_{1}$, $N_{s}$ is set to 0.

In order to generate the real-space correlation function we followed the scheme set out in Ross & Brunner (2009). Firstly, we modeled the power spectrum contributed by galaxies in a single halo (1-halo term) and those in separate haloes (2-halo term, $P_{\rm 2h}$). Also the power spectrum for the 1-halo term is split into that of central-satellite pairs ($P_{\rm cs}(k)$) and satellite-satellite pairs ($P_{\rm ss}(k)$). The equations for each term are

and

where $n(M)$ is the halo mass function as parameterised in Tinker et al. (2010), $u(k|M)$ is the Fourier transform of the halo density profile of Navarro, Frenk & White (1997) and $P_{\rm mat}(k)$ indicates the matter power spectrum at the redshift of the sample. The term $g(k,r)$ can be thought of as an asymptotic bias. The dependence on $r$ arises because only haloes with virial radii less than half of the pair separation $r$ of interest are considered; more massive haloes would experience an exclusion effect at separation $r$ (see the discussion in Appendix B of Tinker et al. 2005). In addition the virial mass ($M_{vir}(r)$) is defined as

where $\overline{\rho}$ is the mean comoving background density. To generate $P_{\rm mat}(k)$ we used the ‘CAMB’ software package (Lewis, Challinor & Lasenby 2000) including the fitting formulae of Smith et al. (2003) to model non-linear growth. The average number density of galaxies ($n_{\rm g}$) is expressed as

A halo model with three free parameters ($\sigma_{\rm cut},M_{0}$ and $\alpha$) was used to produce the angular correlation function. In this case, $M_{\rm cut}$ was fixed by matching the observed number density using the other given parameters, and $M_{1}$ was set as $M_{\rm cut}$ which was found to be suitable in previous studies (Zehavi et al. 2011). The modeled correlation function was projected to angular space using the Limber equation (Limber 1954). Then the covariance matrix was used to find the best fitting parameters having the minimum $\chi^{2}$ value. The fitting range used was $0.001^{\circ}<\theta<0.33^{\circ}$, where the influence of the integral constraint is minimal.

From the fitted parameters, the effective mass ($M_{\rm eff}$), the effective bias ($b_{\rm g}$) and the satellite fraction can be estimated using

and

In order to determine the properties of haloes hosting EROs at different redshifts, we compare all the fitted and estimated values for EROs in three redshift bins (see § 4.3).

It is known that the clustering properties of EROs depend on magnitude and colour (Daddi et al. 2000; Roche et al. 2002, 2003; Brown et al. 2005; Georgakakis et al. 2005; Kong et al. 2006, 2009; Gonzalez-Perez et al. 2011; Kim et al. 2011a). In this section we discuss the properties of the angular two-point correlation function of EROs from the DXS/Subaru and the DXS/PS1 samples.

Fig. 3 shows the correlation functions of EROs selected using various criteria. All correlation functions in this plot show a clear break at $\sim 0.02^{\circ}$, which corresponds to $\sim 1h^{-1}$ Mpc at $z\sim 1$ in comoving coordinates. This break was already reported in Kim et al. (2011a) and implies that a single power-law cannot properly describe the correlation function of EROs. The presence of significant larger scale clustering was confirmed by Kim et al. (2011a) with a detailed analysis of multiple sub-areas of the field and over different ranges in magnitude. Therefore we tried to fit a power-law ($w(\theta)=A_{w}\theta^{-\delta}$) to these correlation functions on small, $0.001^{\circ}<\theta<0.02^{\circ}$, and large, $0.02^{\circ}<\theta<0.33^{\circ}$, scales separately. The boundaries for this fitting were selected to minimise the influence of the integral constraint on the largest scales. The values measured are listed in the top four rows of Table 1. The amplitudes measured for redder or brighter EROs are larger than found for bluer or fainter samples. These features can also be seen in Fig. 3. The top two panels display the dependence of clustering on limiting magnitude, and the third panel from the top shows the colour dependence, which are based on the DXS/Subaru dataset.

In order to check the consistency of these measurements we compared the results with those from the DXS SA22 field in Kim et al. (2011a) which showed good agreement with previously published results. The slopes ($\delta_{small},\delta_{large}$) in Kim et al. were ($0.99\pm 0.09,0.40\pm 0.03$) for $K_{AB}<20.7,(i-K)_{AB}>2.95$ EROs and ($1.00\pm 0.05,0.51\pm 0.02$) for $K_{AB}<20.7,(i-K)_{AB}>2.45$ EROs. The values using the same criteria in this work are ($1.05\pm 0.15,0.72\pm 0.15$) and ($1.10\pm 0.07,0.68\pm 0.12$). On small scales these values are in agreement within the uncertainty range. However, the correlation functions measured in this work are slightly steeper than previous results, particularly on the largest scales. Since Kim et al. used a smaller area, these results might be more affected by cosmic variance, explaining the differences on the large scales. The most pertinent point is that the correlation function is steeper on small scales and flatter on large scales than the single power-law with $\delta=0.8$ assumed in most previous studies. Furthermore, to compare amplitudes directly, we measured the amplitudes again with fixed slopes of $\delta=0.99$ and $0.40$ for small and large scales respectively. For $K_{AB}<20.7,(I-K)_{AB}>3.05$ EROs, the amplitudes ($A_{w}^{small},A_{w}^{large}$) were ($4.14\pm 0.3,42.05\pm 0.9$)$\times 10^{-3}$ in Kim et al., and ($3.66\pm 0.5,33.26\pm 3.7$)$\times 10^{-3}$ in this work. These are also consistent on small scales but not on the largest scales. The reason why the uncertainty is larger in this work than the previous one is that we used the Jackknife resampling method in this study. It is known that the Poissonian error underestimates the uncertainty on the correlation function, especially on the largest scales (Ross et al. 2007; Sawangwit et al. 2011; Nikoloudakis, Shanks & Sawangwit 2013). We also note that the Jackknife resampling method works well on large scales (Scranton et al. 2002; Zehavi et al. 2005). Based on the scale of inflection and the values of amplitude and slope, we conclude that our measurements are consistent with previous results and improve upon them by studying a larger area of sky.

In the bottom panel of Fig. 3, we compare the angular correlation functions from the DXS/Subaru and the DXS/PS1 samples. As mentioned above the area for DXS/PS1 is larger than that of DXS/Subaru, but DXS/Subaru is deeper. The correlation functions of EROs from different optical datasets are well matched, suggesting that there is no significant bias for DXS/PS1 caused by the optical dataset. We use DXS/PS1 EROs to constrain the halo properties in the next sections. Overall, the measurements in this work are consistent with previous work, although even wider data are necessary to measure the angular clustering on large scales more accurately and to test for field to field variations due to cosmic variance. In the near future, PS1 will cover a deeper magnitude range and provide more reliable measurements for the clustering of high redshift galaxies in all four DXS fields.

Photometric redshifts provide an opportunity to compare the clustering properties of EROs in different redshift bins. From the DXS/PS1 catalogue we classified EROs based on their redshift and absolute magnitude. Firstly, we applied a $K$-band absolute magnitude cut ($M_{K}<-23$) with $(i-K)_{AB}>2.45$ to select EROs with a similar $K$-band luminosity in different redshift bins. The $K$-band absolute magnitude was calculated using the k-correction value obtained from the Bruzual & Charlot (2003) GALAXEV code and a luminosity distance calculated using the Javascript Cosmology Calculator (Wright 2006). We assumed a formation redshift of $4<z_{f}<5$ and considered a range of metallicities around solar to find a best fit to $izJK$ colours of the ERO. The simple stellar population templates with a Chabrier (2003) initial mass function were used. Then, the EROs were split into three redshift bins : $1.0<z<1.2,1.15<z<1.45$ and $1.4<z<1.8$. The bin size was chosen to select sufficiently large samples to enable accurate clustering measurements. The angular two-point correlation functions for each sample were measured using the probability distribution function of photometric redshift as described in § 3.1. Finally a power-law fit was performed on small and large scales separately. The same fitting range as above was used. Fig. 4 shows the angular two-point correlation functions of DXS/PS1 EROs at $1.0<z<1.2$ (top), $1.15<z<1.45$ (middle) and $1.4<z<1.8$ (bottom). Dotted lines indicate the power-law fits on small and large scales separately. The bottom three rows in Table 1 list the fitted results. The number of EROs in each bin is the sum of the weights estimated from the probability distribution function of photometric redshift. Although the correlation function shows a slightly flatter shape in the highest redshift bin than the others, all the estimated power-law slopes have similar values within the uncertainty range.

We also note that the amplitudes on large scales with a fixed power-law slope ($\delta=0.7$) are 0.008$\pm$0.0007 at $1.15<z<1.45$ and 0.009$\pm$0.0008 at $1.4<z<1.8$. The power-law slope of $0.7$ was derived for EROs at $1.0<z<1.2$. These similar amplitudes on large scales indicate a higher bias at higher redshift. This will be discussed in the next section.

In this section we study the properties of those haloes hosting EROs at different redshifts. The angular correlation functions for DXS/PS1 EROs in different redshift bins with $M_{K}<-23$ and $(i-K)_{AB}>2.45$ were used for the halo modeling. The halo models were generated at the median redshift of the bins, $z=1.1,1.3$ and $1.5$. The halo model with the three free parameters ($\sigma_{\rm cut}$, $M_{0}$ and $\alpha$) mentioned in § 3.2 was applied at each redshift. We note that this assumed HOD frame work is basically appropriate for mass or luminosity limited samples. In fact EROs are not mass limited samples. Therefore some central or satellite galaxies may be missed, which affects the shape of HODs. In this section, we simply apply the standard HOD model to the clustering of EROs. Then a modified HOD will be discussed in the next section.

Fig. 5 shows the angular correlation function estimated from the observed $(i-K)_{AB}>2.45$ EROs, brighter than $M_{K}=-23$ in the different redshift bins (circles), and the best fit halo models (lines). The dotted and dashed lines represent the 1- and 2-halo terms, respectively. The features of the correlation functions of all ERO subsets are relatively well fitted by the standard halo model. The HOD fit parameters are listed in Table 2. The errors on the HOD parameters in the fits were by determined finding the minimum and maximum parameter values with $\Delta\chi^{2}\leq 1$ from the best fit solution. For other derived parameters ($b_{g}$, $M_{eff}$ and $f_{sat}$), $\Delta\chi^{2}\leq 3.53$ was used.

The top panel in Fig. 6 displays the effective bias which is estimated using the equation in § 3.2. The EROs at higher redshift show a higher bias which is similar to the trend found in previous studies of various populations (Blake et al. 2008; Wake et al. 2008a; Sawangwit et al. 2011 for Luminous Red Galaxies (LRG), Matsuoka et al. 2011; Wake et al. 2011 for stellar mass limited samples and Ross, Percival & Brunner 2010 for absolute magnitude cut samples). The biases of EROs reported previously were $\sim 3$ at $z=2.1$ for $(R-K)_{AB}>3.3$ and $K_{AB}<22.1$ EROs from the semi-analytical model in Gonzalez-Perez et al. (2011) and $2.7\pm 0.1$ in Moustakas & Somerville (2002), which are similar to our measurements, although Moustakas & Somerville (2002) applied a single power-law and used a different criteria, $(I-H)_{AB}>2$ and $H_{AB}<21.9$, for the Las Campanas Infrared Survey data (McCarthy et al. 2001; Firth et al. 2002). These values are similar to those for low redshift LRGs which have $b_{\rm g}\cong 2-3$ (Blake et al. 2008; Wake et al. 2008a; Sawangwit et al. 2011). If EROs have similar stellar masses to LRGs, EROs should be more biased than the lower redshift LRGs. However, we note that the median stellar mass of an SDSS LRG is $\sim 10^{11.5}M_{\odot}$ with a narrow distribution (Barber, Meiksin & Murphy 2007), but in the case of EROs at $K_{AB}<21.6$, the distribution shows a peak at $\sim 10^{11.3}M_{\odot}$ and a sharp cut-off at $\sim 10^{11.5}M_{\odot}$ (Conselice et al. 2008). Therefore our samples are probably marginally less massive than LRGs at lower redshifts.

The middle panel in Fig. 6 shows the effective mass of dark matter haloes hosting EROs. We found that the average halo mass hosting EROs is over $10^{13}h^{-1}M_{\odot}$ and that EROs at higher redshift are in slightly less massive haloes than those at lower redshift. Gonzalez-Perez et al. (2011) reported that the median mass of haloes hosting $(R-K)_{AB}>3.3$ and $K_{AB}<20.9$ EROs at $z=1.1$ is $4.4\times 10^{12}h^{-1}M_{\odot}$ from a semi-analytical model, and Moustakas & Somerville (2002) estimated the average halo mass hosting EROs as $5\times 10^{13}h^{-1}M_{\odot}$ at $z\sim 1.2$. Recently Palamara et al. (2013) published the halo mass of $K_{AB}<18.9$ EROs as $10^{13.09}$ $M_{\odot}$ based on the clustering strength with a power law slope of 1.73. Since we have dealt with a full halo model in this work rather than a simplified conversion, the result may be the better measurement than previous work for EROs.

Other studies use a variety of definitions of “passive” galaxies at $1<z<2$, which have differing degrees of overlap with the ERO sample considered in this paper. Hartley et al. (2010) used star formation history modelling to define galaxies as passive, selecting those in which they inferred that the current star formation rate is less than 10% of the initial star formation rate and using a red cut on the rest-frame $U-B$ colour. These authors estimated that passive galaxies defined in this way with $M_{K}<-23$ at $z<2$ are located in haloes ranging from $10^{13}M_{\odot}$ to $5\times 10^{13}M_{\odot}$. We infer a host halo mass that is slightly higher than the prediction of Gonzalez-Perez et al. (2011), but slightly lower than that for the passive galaxies as defined by Hartley et al. (2010). The comparison with the theoretical work is discussed further in § 5. However, overall, EROs may reside in slightly less massive haloes than passive galaxies of comparable luminosity. It is known that EROs can be split on the basis of their colours into passive galaxies with old stellar populations or dusty, star-forming galaxies (Pozzetti & Mannucci 2000; Smail et al. 2002; Roche et al. 2002; Cimatti et al. 2002, 2003; Moustakas et al. 2004; Sawicki et al. 2005; Simpson et al. 2006; Conselice et al. 2008; Kong et al. 2009). In Kim et al. (2011a), the fraction of old, passive EROs defined in this way was found to be more than $\sim 60$ per cent. In this work it is not possible to split the whole sample into these two sub-populations by using $(i-K)$ and $(J-K)$ colours due to the lack of $J$-band imaging over the full area. So, we simply apply the criterion to the two-colour diagram for the region where $J$-band imaging exists. The fraction of old, passive EROs is 45.4 per cent at $z=1.1$ and 54.4 per cent at $z=1.5$. The significant fraction of dusty, star-forming galaxies may dilute the clustering of EROs which might explain why we find lower halo masses than pure passive galaxy samples would predict.

The bottom panel in Fig. 6 shows the evolution of the satellite fraction with redshift derived from the HOD model. EROs are made up of a larger fraction of satellite galaxies at lower redshifts, although the best fitting power-law slope of satellites ($\alpha$) at $z=1.5$ is much larger than unity, which is the value typically seen in simulations of galaxy formation (e.g. Almeida et al. 2008) and as recovered in many previous analyses of observational measurements of clustering. However, in Wake et al. (2008a), the slope for satellite LRGs at $z=0.55$ was $\sim$2.0, and Matsuoka et al. (2011) also reported a similar value for massive galaxies at $0.8<z<1.0$. So this is not the first result showing a large slope. However, we also note there are other components of the HOD that could lead to a low satellite fraction other than a steep satellite slope. The mass thresholds ($M_{\rm cut}$ and $M_{0}$) for the HOD may be factors causing the result. The threshold for the HOD of central EROs ($M_{\rm cut}$) is smaller than that for LRGs or massive galaxies, which are a few times $10^{13}h^{-1}M_{\odot}$. Moreover, $M_{0}$ is larger than SDSS galaxies at $z\sim 0.3$ (Ross et al. 2010). This means that EROs are in less massive haloes than LRGs when they are the central galaxy and in more massive haloes than typical galaxies when they are a satellite. Therefore, these effects may lead to the lower satellite fraction than the previous results for ordinary galaxies.

As mentioned above, EROs can be split into old, passive galaxies (OG) or dusty, star-forming galaxies (DG) on the basis of their $i-K$ and $J-K$ colours. Kim et al. (2011a) show that these two sub-populations have very different clustering properties with OGs being more strongly clustered than DGs. Unfortunately, the $J$-band coverage of the EN1 field is not complete so we cannot perform the same analysis as Kim et al. (2011a). Gonzalez-Perez et al. (2011) determine the difference in clustering of these sub-populations in semi-analytical simulations and find a comparable difference in the clustering to that in Kim et al. (2011a). Therefore it is likely that OGs are more biased and/or in more massive haloes than DGs but the clustering of DGs is sufficiently similar to that of OGs that any halo modeling of the combined population is representative. We will perform a more detailed halo modeling of each sub-population with the full DXS dataset in a future paper.

Additionally we note the results of halo modeling for stellar mass limited samples. Wake et al. (2011) used data from the NEWFIRM medium band survey (NMBS; van Dokkum et al. 2009; Brammer et al. 2009; van Dokkum et al. 2010; Whitaker et al. 2011) to measure the clustering of stellar mass limited samples of galaxies at $1<z<2$. The highest mass limits of their samples were $M_{*}=10^{10.7}M_{\odot}$ at $z=1.1$ and $10^{10.78}M_{\odot}$ at $z=1.5$. Comparing our results to theirs, the effective mass and bias for EROs are higher than those in Wake et al.(2011). This can be easily explained by the higher stellar mass of EROs, as more massive or brighter galaxies reside in more massive haloes (Zehavi et al. 2005, 2011; Foucaud et al. 2010; Hartley et al. 2010; Matsuoka et al. 2011; Furusawa et al. 2011) and show a higher bias (Coil et al. 2006; Ross & Brunner 2009; Zehavi 2011), if we assume EROs have stellar masses greater than $10^{11}M_{\odot}$.

Foucaud et al. (2010) measured the mass of haloes hosting $10^{11}M_{\odot}<M_{*}<10^{11.5}M_{\odot}$ Palomar/DEEP2 galaxies at $1.2<z<1.6$ based on the model in Mo & White (2002). They reported the halo mass and bias as $10^{13.17}h^{-1}M_{\odot}$ and 2.8$\pm$0.6, which are similar to our result for EROs at $z=1.3$. Therefore EROs may have similar properties with those galaxies.

We start with Fig. 7 in which we compare our new observational estimate of the angular correlation function of EROs with $(i-K)_{AB}>2.45$ and $M_{K}<-23$ (filled circles) with the clustering predicted by GALFORM for such galaxies (open circles). The model galaxies were extracted with the same colour cut and an apparent magnitude³³3The magnitudes are $K_{AB}=20.7$ and $21.0$ for $z=1.1$ and $1.5$, respectively. corresponding $M_{K}=-23$ at $z=1.1$ and 1.5. The theoretical clustering is obtained from the ERO HOD output directly by GALFORM at a given redshift, along with the predicted redshift distributions. We note that this may be different from previous works about the clustering based on GALFORM and Gonzalez-Perez et al. (2011). The lines in Fig. 7 show the angular correlation function obtained by fitting different parametric forms for the HOD to the clustering predicted by GALFORM and will be discussed further in § 5.3.

This method of deriving the angular correlation function from GALFORM does not readily yield an appropriate estimate of the error on the angular clustering, as the effective volume considered, the volume of the Millennium N-body simulation, is much larger than that covered by the observational data. Instead, we take the fractional errors inferred on the observational estimate of the correlation function, and apply these to the GALFORM predictions (note this are not plotted in Fig. 7). In the future it will be possible to extract exactly the same volume sample from the models, using the lightcone mock catalogue building capability developed by Merson et al. (2013).

The upper panel of Fig. 7 shows that at $z=1.1$ the GALFORM predictions show slightly stronger angular clustering than is observed, with some tension at the 2–3$\sigma$ level (see § 5.2). On large scales, the observational estimate may be affected by sampling variance. On the other hand, the discrepancy at small angular separations ($\theta<0.02$) can be interpreted as the model predicting a larger 1-halo term than is suggested by the observations. This in turn implies that the model predicts that too many EROs are satellite galaxies at this redshift.

The lower panel of Fig. 7 shows the angular clustering at $z=1.5$. In this case, the predictions agree remarkably well with the new observational estimate on large scales. The larger volume covered at $z=1.5$ than at $z=1.1$ reduces the effect of cosmic variance. However, the discrepancy on small scales also implies a larger satellite fraction in the models compared with the observations, in the same sense as that suggested at $z=1.1$.

There are several factors which could be responsible for the discrepancy between the observed and predicted clustering. Firstly, as EROs are unusual objects drawn from the tail of luminosity and colour distributions, the precise colour criterion applied to construct the ERO samples can have a substantial impact on the predictions. Small differences between the predicted and observed colour distributions can result in large variations in the number of galaxies selected and their clustering. Gonzalez-Perez et al. (2011) noted that the agreement between theory and observations was greatly improved on perturbing the colour cut used to select EROs in the model. Secondly, the semi-analytical modeling predicts too many red satellites resulting in a higher clustering amplitude in the 1-halo term. As already mentioned, most EROs predicted by GALFORM are quiescent galaxies (Gonzalez-Perez et al. 2009). In addition, the clustering amplitude of quiescent EROs is signicantly higher than dusty, star-forming EROs on small scales in the model (Gonzalez-Perez et al. 2011). However observational results show a smaller discrepancy than this model prediction (Miyazaki et al. 2003; Kim et al. 2011a). Hence, it may be necessary to amend the prescription for satellite formation in the models (Kim et al. 2009; Kimm et al. 2009). This issue will be addressed in more detail in § 5.4. Thirdly, even for the large number of EROs used in our analysis, the relatively small area of the fields (a combined areas of just over $5$ square degrees) means that the clustering estimates are susceptible to sampling variance. This effect is taken into account to some extent in the Jackkinfe errors plotted on our measurements. However, Norberg et al. (2009) demonstrate that an internal estimate of the error such as that obtained using Jackknife resampling can still vary in amplitude between different realisations of the data, particularly if the intrinsic clustering is strong, as is the case with EROs.

We now compare the HODs derived by fitting the parametric form of Zheng et al. (2005) to the observed clustering and to the clustering predicted by the model. To allow a meaningful comparison, we perform the HOD analysis of the observed clustering again adopting the background cosmology used in the Bower et al. (2006) model, which matches that used in the Millennium N-body simulation of Springel et al. (2005)⁴⁴4The cosmological parameters used in the Millennium Simulation are $\Omega_{\rm m}=$0.25, $\Omega_{\Lambda}=$ 0.75, $\sigma_{8}=0.9$ and $H_{0}=$ 73 km s^-1 Mpc^-1..

The results of this exercise are shown by the dotted lines in Fig. 7, which show the angular clustering obtained using the parametric form for the HOD given by Eqs. 5, 6 & 7, with parameters chosen to give the best fit to the theoretical angular clustering obtained from the GALFORM predicted HOD. Table 3 lists the best fitting HOD parameters and some derived quantities. The effective bias factors and host halo masses deduced from the best fitting HOD to the observed clustering are higher than those derived from the fits to the predicted clustering. The effective masses are a factor of $\sim 1.9$ higher at $z=1.1$ and a factor $\sim 2.7$ higher at $z=1.5$ in the fits to the observations compared with the fits to the model. Another discrepancy between models and observations regards the fraction of EROs that are satellites, which is $\approx 2.5$ times higher in the model fits than in the fits to the measured clustering.

The main panels of Fig. 8 display the best fit central HODs to the angular clustering of EROs in the observations (solid line) and the GALFORM predictions (dotted line) at $z=1.1$ (upper) and $z=1.5$ (lower). The insets of Fig. 8 show the HODs for all galaxies, combining the HODs of central and satellite EROs. The best fit HOD to the predicted clustering extends to lower mass haloes than the fit to the observed clustering. In addition, satellites in the model also reside in less massive haloes than the observations. These departures may explain the higher halo masses and lower satellite fractions in the observations than in the model. The form of the HOD output by GALFORM will be discussed in next subsection.

The galaxy HOD is a prediction of GALFORM and not an input. GALFORM models the physics of the baryonic component of the Universe to predict the number of galaxies per halo and their properties. The HOD is extracted by applying the observational selection to the model output and simply counting the number of galaxies which are retained in each dark matter halo, distinguishing between the central galaxy and its satellites. The HOD and the quantities derived from it can be extracted directly from the model, without having to go through the intermediate step of fitting a parametric form for the HOD to the predicted clustering. As we will see later on in this subsection, the form of the actual HOD output by GALFORM can be different from the standard parametrization we have adopted so far.

We begin by listing in Table 4 some basic properties predicted by GALFORM for the EROs and their host haloes. The effective host halo mass defined by Eq. 12 with the direct output HOD of GALFORM is in very close agreement with that derived from the HOD fit to the predicted angular clustering listed in Table 3. This means that EROs predicted by GALFORM are in less massive haloes than the observations. The major differences between the model prediction and observations are number density and satellite fraction. The predicted number density of EROs is much larger than the observational estimate in the previous section. Also the fraction of EROs that are satellites from the direct model prediction is around 4–5 times higher than the fraction obtained from the best fitting HOD to the observation.

We plot the HOD predicted directly by GALFORM in Fig. 8 (dashed line). We compare this with the parameterised form for the HOD proposed by Zheng et al. (2005), which was motivated by earlier simulations of galaxy formation which did not include AGN feedback. In the Zheng et al. (2005) framework, the HOD of central galaxies is assumed to reach unity, i.e. above some halo mass, every central galaxy is assumed to meet the selection criteria of the sample. The parametric HOD (dotted line) plotted in Fig. 8 is the best fit to the angular clustering predicted by GALFORM.

The predicted HOD can differ substantially from the canonical form typically assumed for the HOD as shown in Fig. 8 and remarked upon by Contreras et al. (2013) in their comparison of the HOD predicted in different semi-analytical models. The deviation is largely driven by AGN feedback, which shuts down gas cooling in massive haloes in GALFORM. In general, the onset of AGN feedback above some halo mass alters how galaxy properties vary with halo mass. This could be manifest as a dramatic break in the relation between the galaxy property and halo mass, as in the case of cold gas mass (Kim et al. 2011b), or as a change in the slope and scatter of the relation, as in the case of $K$-band luminosity (Gonzalez-Perez et al. 2011). With AGN feedback, the most massive galaxy in a sample may no longer be in the most massive dark matter halo due to the increased scatter in the correlation between galaxy properties and halo mass.

Fig. 8 shows the HOD of central galaxies that are EROs at $z=1.1$ (upper) and $1.5$ (lower), respectively. The predicted HOD of central EROs differs from the observationally inferred in terms of the transition from zero galaxies per halo and also in the number of central galaxies which are EROs. At both redshifts, the mean number of centrals is below unity for the predicted HOD. In addition the fitted HOD shows a smoother shape than the GALFORM prediction. These features emphasize the importance of the form adopted for the HOD when interpreting the results of halo modeling.

The modified HOD fits the HOD of central EROs in GALFORM much better at both redshifts than was possible using the Zheng et al. (2005) HOD. This emphasizes the importance of using a realistic form for the HOD form to derive robust halo properties from the galaxy distribution.

The comparisons above point to GALFORM predicting that, overall, EROs are found in less massive haloes and that more EROs are satellites compared with the conclusions reached by fitting HOD models to the observed clustering.

Both problems could be solved by changing the treatment of cooling gas in GALFORM. Gonzalez-Perez et al. (2011) looked at the predictions of a different semi-analytical model, that of Font et al. (2008). The Font et al. model includes partial stripping of the hot gas halo from satellite galaxies, whereas the default assumption is that all of the hot gas is stripped from a galaxy once it becomes a satellite within a larger halo. In the Font et al. model, depending on the orbit of the satellite and the ram pressure that it experiences, some hot gas may be retained and can cool onto the galaxy even after it has become a satellite. This change to the model makes satellites bluer, by permitting more star formation to take place. However, Font et al. also invoked a factor of two increase in the stellar yield without changing the IMF, which leads to redder galaxy colours. Gonzalez-Perez et al. report that the Font et al model produces stronger clustering for EROs than the model of Bower et al. but also leads to more EROs than are observed.

Additionally, Gonzalez-Perez et al. (2009) mentioned that the Bower et al. (2006) model quenched the star formation of massive galaxies too efficiently by the AGN feedback, based on the predicted redshift distribution of passive galaxies. Zheng et al. (2009) also pointed out that this model predicts more red galaxies in a few times $10^{12}h^{-1}M_{\odot}$ haloes than observations but less in more massive haloes, comparing the HODs of observed and predicted LRGs. Therefore it may be possible that GALFORM predicts more red central galaxies such as EROs and LRGs in lower mass haloes.

These comparison illustrate the potential of our new clustering results to constrain the modeling of different elements of the physics of galaxy formation.