THE GALAXY CLUSTER MID-INFRARED LUMINOSITY FUNCTION AT 1.3 < z < 3.2

As part of the Spitzer snapshot program CARLA, we observed the fields of 209 high-redshift radio galaxies (HzRGs) and 211 RLQs at 1.3 < z < 3.2. A description of the sample selection, observation strategy, data reduction, source extraction, determination of completeness limits, and initial scientific results are given in Wylezalek et al. (2013).

Briefly, the fields, covering 52 × 52 each, corresponding to a physical size of ∼2.5 × 2.5 Mpc for the redshift range of the targeted RLAGN, were mapped with the Infrared Array Camera (IRAC; Fazio et al. 2004) on board the Spitzer Space Telescope at 3.6 and 4.5 μm (referred to as IRAC1 and IRAC2). The total exposure times give similar depths for both channels with the IRAC2 observations having been tailored to be slightly deeper for our intended science. The HzRG and RLQ samples have been matched with respect to their redshift and radio-luminosity distributions. The L_500MHz–z plane is covered relatively homogeneously in order to be able to study the environments around the AGN as a function of both redshift and radio luminosity.

The data were reduced using the MOPEX (Makovoz & Khan 2005) package and source extraction was performed using SExtractor (Bertin & Arnouts 1996) in dual image mode with the IRAC2 image serving as the detection image. The CARLA completeness was measured by comparing the CARLA number counts with number counts from the Spitzer UKIDSS Ultra Deep Survey (SpUDS, PI: J. Dunlop) survey. The SpUDS survey is a Spitzer Cycle 4 legacy program that observed ∼1 deg² in the UKIDDS UDS field with IRAC and the Multiband Imaging Spectrometer (Rieke et al. 2004) aboard Spitzer. SpUDS reaches 5σ depths of ∼1 μJy (mag ≃ 24). We determine 95% completeness of CARLA at magnitudes of [3.6] = 22.6 (=3.45 μJy) and [4.5] = 22.9 (=2.55 μJy).

The cluster candidates studied in this paper were identified as IRAC color-selected galaxy overdensities in the fields of CARLA RLAGN. In order to isolate galaxy cluster candidates, we first identify high-redshift sources (z > 1.3) by applying the color cut [3.6]−[4.5] ⩾−0.1 (e.g., Papovich 2008). This color selection has proven to be very efficient at identifying high-redshift galaxies independent of their evolutionary stage. A negative k-correction, caused by the 1.6 μm bump that enters the IRAC bands at z ∼ 1, leads to an almost constant IRAC2 apparent magnitude at z > 1.3 and a red [3.6]−[4.5] color (Stern et al. 2005; Eisenhardt et al. 2008; Wylezalek et al. 2013). We apply a counts-in-cells analysis to such color-selected IRAC sources, which we simply refer to as IRAC-selected sources, to identify overdense fields.

Specifically, we define a source to be an IRAC-selected source if (1) it is detected above the IRAC2 95% completeness limit and above an IRAC1 flux of 2.8 μJy (3.5σ detection limit) and has a color of [3.6]−[4.5] >−0.1 or (2) it is detected above the IRAC2 95% completeness limit and has an IRAC1 flux <2.8 μJy but has a color > − 0.1 at the 3.5σ detection limit of the IRAC1 observation (e.g., Figure 1). This refined criterion means that all IRAC-selected sources are 95% complete in IRAC2 down to [4.5] = 22.9 (=2.55 μJy) but are not necessarily well-detected in IRAC1.

Zoom In Zoom Out Reset image size

We measure the density of IRAC-selected sources in a radius of 1 arcmin centered on the RLAGN and compare it to the mean blank-field density. A radius of 1 arcmin corresponds to ∼500 kpc over the targeted redshift range and matches typical sizes for z > 1.3 mid-IR selected clusters with log(M₂₀₀/M_☉) ∼ 14.2 (e.g., Brodwin et al. 2011). The typical blank-field density of IRAC-selected sources was measured by placing roughly 500 independent apertures of 1 arcmin radius onto the SpUDS field with its catalogs having been cut at the CARLA depth. Since the publication of Wylezalek et al. (2013), 33 new CARLA fields have been observed, and we provide the full table of CARLA fields and their overdensities in the online journal. A short, example version is presented in Table 1.

In Wylezalek et al. (2013), we used the IRAC1 95% completeness limit of 3.45 μJy as the IRAC1 limit in the color selection process, rather than the 3.5σ detection limit of 2.8 μJy used above. As the IRAC1 95% completeness limit is brighter than the IRAC2 95% completeness limit, this IRAC1 limit together with the color criterion introduced an artificially brighter IRAC2 limit ([4.5] = 22.7, 3.0 μJy). In this work, we therefore raise the IRAC1 limiting magnitude to [3.6] = 22.8 (=2.8 μJy), i.e., we lower the flux density cut, to include all IRAC2 sources down to the formal completeness limit of [4.5] = 22.9. This new IRAC1 flux density cut still corresponds to a >3.5σ detection. This is necessary as we are aiming to include as many faint sources as possible in the analysis. We apply this refined color criterion to the CARLA and SpUDS fields and plot the distribution of their densities in Figure 2.

Zoom In Zoom Out Reset image size

Promising galaxy cluster fields are defined as fields that are overdense by at least 2σ (Σ_CARLA > Σ_SpUDS + 2σ_SpUDS; Wylezalek et al. 2013). This criterion is met by 46% of the CARLA fields: 27% are ⩾3σ overdense and 11% are overdense at the 4σ level. Excluding bad fields (e.g., fields that are contaminated by nearby bright stars), this provides us with 178, 101, and 42 high-redshift galaxy cluster fields at the 2σ, 3σ, and 4σ level, respectively.

Note that this selection is a measure of the overdensity signal compared to a blank field, i.e., to a distribution of cell densities centered on random positions. Since the sample is selected using specific galaxies, RLAGN, there will always be at least one source in the overdensity. As all galaxies are clustered to some extent (Coil 2013), the exact measurement of the overdensity signal of RLAGN would be to compare it to a distribution of background cells centered on random galaxies not random positions.

However, the average number density of IRAC-selected sources in SpUDS is very high, on average 34 sources per aperture with radius r = 1 arcmin. This means that if the IRAC-selected sources were extended and filled the aperture maximally, their radius would only be r_ext = 10.3 arcsec and the average distance between two sources would be 2 × r_ext = 20.6 arcsec. We repeated the blank-field analysis by measuring the densities in roughly 500 apertures centered on random SpUDS IRAC-selected sources. Because of the large aperture radius of r = 1 arcmin and small r_ext, the difference between the two background measurements is not significant (Σ_SpUDS = 10.3 ± 2.6 arcmin⁻² compared to Σ_SpUDS = 9.6 ± 2.1 arcmin⁻², see Section 3.2). Since we will work with the blank field background for the luminosity function (LF) analysis and to be consistent with Wylezalek et al. (2013), we show the SpUDS blank field distribution in Figure 2.

Wylezalek et al. (2013) shows that the radial density distribution of the IRAC-selected sources is centered on the RLAGN, implying that the excess IRAC-selected sources are associated with the RLAGN. For the following analysis, we assign the redshift of the targeted RLAGN to the IRAC-selected sources in the cell, and we study the evolution of these galaxy cluster member candidates as a function of redshift. To exclude problematic cluster candidates, we also checked that the median [3.6]−[4.5] color of the IRAC-selected sources per field is in agreement with the [3.6]−[4.5] color expected for a source with the redshift of the targeted RLAGN. No obvious problematic fields were found. For the rest of the manuscript, we will use the term "galaxy clusters" to refer to these cluster and protocluster candidates.

Luminosity functions, Φ(L), provide a powerful tool to study the distribution of galaxies over cosmological time. They measure the comoving number density of galaxies per luminosity bin, such that

$$\begin{equation} dN = \Phi (L)dLdV, \end{equation} \tag{ 1 }$$

where dN is the number of observed galaxies in volume dV within the luminosity range [L, L + dL].

There are many ways in which Φ(L) can be estimated and parameterized, but the most common of these models is the Schechter function (Schechter 1976)

$$\begin{equation} \Phi (L)dL = \Phi ^{*}\left(\frac{L}{L^{*}}\right)^{\alpha }\exp \left(\frac{-L}{L^{*}}\right)\frac{dL}{L^{*}}, \end{equation} \tag{ 2 }$$

where Φ* is a normalization factor defining the overall density of galaxies, usually quoted in units of h³ Mpc⁻³, and L* is the characteristic luminosity. The quantity α defines the faint-end slope of the luminosity function and is typically negative, implying large numbers of galaxies with faint luminosities. The luminosity function can be converted from absolute luminosities L to apparent magnitudes m and can be written as

$$\begin{equation} \Phi (m)= 0.4\ln (10)\Phi ^{*}\frac{(10^{0.4(m^{*}-m)})^{(\alpha +1)}}{e^{10^{0.4(m^{*}-m)}}}, \end{equation} \tag{ 3 }$$

where m* is the characteristic magnitude.

In this work, we study the evolution of m* in galaxy clusters as a function of redshift for both α as a free fitting parameter and fixed α = −1. We include all CARLA fields that are overdense at the 2σ level or more.

We measure the IRAC2 luminosity function of the galaxy clusters from the CARLA survey as a function of redshift, defined as the redshift of the AGN, and galaxy cluster richness. The CARLA galaxy cluster richness is defined in terms of the significance of the overdensity of IRAC-selected sources within the cell centered on the RLAGN (Figure 2). For the largest sample of CARLA clusters, i.e., those overdense at the >2σ level, we measure the luminosity function in six redshift bins chosen in a way that all redshift bins contain the same number of galaxy clusters. We also consider CARLA clusters overdense at the 2.5–3.5σ and >3.5σ level. For these smaller subsamples, we measure the luminosity function in three redshift bins. IRAC1 luminosity functions are presented in the Appendix.

For each CARLA cluster candidate j, we compute the k- and evolutionary corrections (k_j- and e_j) that are required to shift the galaxy cluster members to the center of the redshift bin and apply them to the apparent magnitudes of the galaxy cluster members. The corrections were computed using the publicly available model calculator EzGal (Mancone & Gonzalez 2012) using the predictions of passively evolving stellar populations from Bruzual & Charlot (2003) with a single exponentially decaying burst of star formation with τ = 0.1 Gyr and a Salpeter initial mass function. The k_j+e_j corrections are of order 0.04 mag and do not have a significant impact on the LF analysis; they are simply applied for completeness. The typically bright targeted RLAGNs have been removed from the analysis.

For each CARLA cluster candidate j, we measure the number of cluster members $n_{m_{i}}$ (in a cell with radius r = 1 arcmin centered on the RLAGN) in the ith magnitude bin with m_i = [m_i, m_i + δm_i]:

$$\begin{equation} \Phi _{i_{\rm {CARLA}}, j, k_j+e_j } = n_{m_{i}}. \end{equation} \tag{ 4 }$$

This number density is a superposition of the cluster luminosity function and the luminosity function of background/foreground galaxies. In the following, we describe the statistical background determination and how the true galaxy cluster luminosity function is determined.

The CARLA fields cover an area of ∼52 × 52, roughly corresponding to a region with a radius of 1–1.5 Mpc for the typical redshift of the RLAGN. This radius is in good agreement with sizes of typical mid-IR selected clusters (e.g., Brodwin et al. 2011). Therefore, a local background subtraction in each field is not possible. Instead, we determine a global background in a statistical way using the SpUDS survey.

As described in Section 2.2, we placed roughly 500 random, independent (i.e., non-overlapping) apertures with radius r = 1 arcmin onto the SpUDS survey to estimate the typical blank field density of IRAC-selected sources. Fitting a Gaussian to the low-density half of the SpUDS density distribution finds a mean surface density of Σ_SpUDS = 9.6 arcmin⁻² with a width of σ_SpUDS = 2.1 arcmin⁻². The tail at larger densities arises, as even a 1 deg² survey has large-scale structure and contains clusters. To determine the mean background, we consider cells in the SpUDS survey with surface densities of IRAC-selected sources in the range 9.6 ± 2.1 arcmin⁻² (see Figure 2).

The average background luminosity function per cell is given by

$$\begin{equation} \Phi _{i_{\rm {BG}}} = n_{m_{i}, \rm {BG}}/N_{\rm {BG}}, \end{equation} \tag{ 5 }$$

where $n_{m_{i}, \rm {BG}}$ is the number of galaxies with magnitudes m_i = [m_i, m_i + δm_i] in the SpUDS background cells and N_BG is the number of SpUDS cells used for the background determination.

We then subtract this average blank field luminosity function from the luminosity function of each CARLA cluster field j after having applied the same k_j- and e_j-correction as for the corresponding CARLA field, such that

$$\begin{equation} \Phi _{i,j} = \Phi _{i_{\rm {CARLA}},j, k_j+e_j } - \Phi _{i_{\rm {BG}}, k_{j}+e_{j}. } \end{equation} \tag{ 6 }$$

After this background subtraction the signal of all CARLA cluster fields per redshift bin is stacked to obtain the background-subtracted luminosity function:

$$\begin{equation} \Phi _{i} = \left(\sum \limits _{j=0}^{N}{\Phi _{i,j}}\right)\times \frac{1}{N}\times \frac{1}{A}, \end{equation} \tag{ 7 }$$

where N is the number of CARLA clusters in the redshift bin and A = 1 arcmin² × π is the area of a cell with a radius of 1 arcmin. Figure 3 shows the luminosity functions for CARLA fields of richness ⩾2σ divided into six redshift bins.

Zoom In Zoom Out Reset image size

The accurate measurement of the background number counts is essential for measuring the true cluster luminosity function. Because the SpUDS survey, used in this paper, covers only ∼1 deg², it needs to be confirmed that it represents a typical blank field and is not significantly affected by cosmic variance. At the depth of the CARLA observations, however, it is the largest contiguous survey accessible.

The 18 deg² Spitzer Extragalactic Representative Volume Survey (SERVS; Mauduit et al. 2012) reaches almost the same depth as CARLA and allows a test of the goodness of SpUDS as a blank field. SERVS maps five well-observed astronomical fields (ELAIS-N1, ELAIS-S1, Lockman Hole, Chandra Deep Field South, and XMM-LSS) with IRAC1 and IRAC2. Coverage is not completely uniform across the fields but averages ∼1400 s of exposure time. We extracted sources from the SERVS images (M. Lacy 2013, private communication) in the same way as we did for CARLA and SpUDS to allow for a consistent comparison. In Figure 4, we show the IRAC2 number counts in SpUDS and SERVS observations of the XMM-LSS field, illustrating the difference in depths of the two surveys. Comparing the SERVS number counts with the SpUDS number counts gives an IRAC2 95% completeness limit of 2.85 μJy for the SERVS observations. As we aim to go as deep as possible for the CARLA analysis, the SERVS 95% completeness is slightly shallow compared to the corresponding depth of 2.55 μJy for the CARLA observations. We therefore use the smaller area SpUDS survey for the background determinations.

Zoom In Zoom Out Reset image size

To test the validity of our background subtraction, we placed 200 random apertures onto the 18 deg² SERVS survey and measured the density of IRAC-selected sources in those random cells. There is a chance that a CARLA cluster field and a non-associated large-scale structure are found in projection. In this case, our background subtraction would underestimate the actual background in that field. By placing 200 random apertures onto the 18 deg² SERVS survey, we can get a qualitative upper limit of this probability. Figure 5 shows the distribution of densities of IRAC-selected sources in those 200 random cells at the magnitude limit of the SERVS survey and compares it to the distribution in SpUDS. SpUDS only covers ∼1 deg² and is known to be biased to contain clusters and large-scale structure. The apertures placed on SpUDS cover almost the full area and will thus pick up any large-scale structure or clusters in that field. That is why a prominent high-density tail arises. By placing random apertures onto a larger survey that is less biased by large-scale structure, the high-density tail becomes less prominent, as expected.

To make sure, however, that we do not underestimate the background luminosity function by not taking enough of the high-density tail into account, we test here what impact our choice of background has on the overall results in this work. We repeat the luminosity function fitting analysis by choosing different density intervals to estimate the background, namely, Σ_SpUDS ± 0.8 × σ_SpUDS and Σ_SpUDS ± 1.6 × σ_SpUDS. Taking an even narrower interval than Σ_SpUDS ± 0.8 × σ_SpUDS gives too few apertures and results in a very noisy background luminosity function. In Figure 6, we show that the results for the luminosity function fits with a fixed α for the original background subtraction (1σ) and for the test background subtractions (0.8σ, 1.6σ). The results for the different runs agree remarkably well and no systematic effect is seen. At Σ_SpUDS = Σ_SpUDS ± 1.6 × σ_SpUDS, we are already sampling part of the high-density tail. Table 2 shows that this has a minimal impact on the Schechter function fits. The exact choice of the density interval of the background subtraction is not critical and illustrates that choosing SpUDS cells with Σ_SpUDS ± σ_SpUDS is a sensible measure for the blank field density of IRAC-selected sources.

Table 2. Schechter Fit Results for Both α Free and α Fixed to −1 to the 4.5 μm Luminosity Function

ΣCARLA	〈z〉	$m^*_{4.5\,\mu {\rm m}}$	α	$m^*_{4.5\, \mu {\rm m}, \alpha = -1}$	N
>2σ	1.45	$19.84 ^{+ 0.30 }_{- 0.20 }$	$-1.01 ^{+ 0.10 }_{- 0.10 }$	$19.85 ^{+ 0.15 }_{- 0.15 }$	28
>2σ	1.77	$20.23 ^{+ 0.30 }_{- 0.20 }$	$-1.07 ^{+ 0.15 }_{- 0.15 }$	$20.32 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.05	$20.19 ^{+ 0.30 }_{- 0.20 }$	$-0.92 ^{+ 0.15 }_{- 0.15 }$	$20.08 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.26	$19.91 ^{+ 0.25 }_{- 0.25 }$	$-1.28 ^{+ 0.15 }_{- 0.15 }$	$20.41 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.51	$19.96 ^{+ 0.35 }_{- 0.20 }$	$-1.10 ^{+ 0.15 }_{- 0.15 }$	$20.10 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.92	$19.60 ^{+ 0.30 }_{- 0.20 }$	$-1.29 ^{+ 0.10 }_{- 0.15 }$	$20.10 ^{+ 0.20 }_{- 0.15 }$	30
2.5 < σ < 3.5	1.65	$20.38 ^{+ 0.20 }_{- 0.20 }$	$-0.75 ^{+ 0.20 }_{- 0.10 }$	$20.10 ^{+ 0.15 }_{- 0.15 }$	25
2.5 < σ < 3.5	2.23	$19.99 ^{+ 0.20 }_{- 0.20 }$	$-1.13 ^{+ 0.20 }_{- 0.15 }$	$20.21 ^{+ 0.15 }_{- 0.15 }$	27
2.5 < σ < 3.5	2.81	$19.74 ^{+ 0.20 }_{- 0.20 }$	$-1.17 ^{+ 0.15 }_{- 0.15 }$	$19.99 ^{+ 0.15 }_{- 0.15 }$	27
>3.5σ	1.49	$19.88 ^{+ 0.20 }_{- 0.25 }$	$-0.89 ^{+ 0.10 }_{- 0.10 }$	$19.72 ^{+ 0.15 }_{- 0.15 }$	18
>3.5σ	1.88	$20.18 ^{+ 0.30 }_{- 0.20 }$	$-1.01 ^{+ 0.20 }_{- 0.15 }$	$20.20 ^{+ 0.25 }_{- 0.15 }$	19
>3.5σ	2.49	$20.23 ^{+ 0.30 }_{- 0.20 }$	$-0.95 ^{+ 0.20 }_{- 0.15 }$	$20.16 ^{+ 0.20 }_{- 0.15 }$	20
>2σ	1.45	$19.81 ^{+ 0.30 }_{- 0.20 }$	$-1.02 ^{+ 0.10 }_{- 0.15 }$	$19.84 ^{+ 0.20 }_{- 0.15 }$	28
>2σ	1.77	$20.25 ^{+ 0.30 }_{- 0.20 }$	$-1.07 ^{+ 0.15 }_{- 0.15 }$	$20.34 ^{+ 0.15 }_{- 0.15 }$	30
>2σ	2.05	$20.16 ^{+ 0.30 }_{- 0.20 }$	$-0.93 ^{+ 0.15 }_{- 0.15 }$	$20.07 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.26	$19.86 ^{+ 0.25 }_{- 0.25 }$	$-1.32 ^{+ 0.15 }_{- 0.10 }$	$20.44 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.51	$19.91 ^{+ 0.20 }_{- 0.30 }$	$-1.12 ^{+ 0.10 }_{- 0.15 }$	$20.09 ^{+ 0.20 }_{- 0.15 }$	30
>2σ	2.92	$19.49 ^{+ 0.30 }_{- 0.25 }$	$-1.33 ^{+ 0.15 }_{- 0.10 }$	$20.09 ^{+ 0.20 }_{- 0.15 }$	30
2.5 < σ < 3.5	1.65	$20.39 ^{+ 0.20 }_{- 0.20 }$	$-0.73 ^{+ 0.15 }_{- 0.15 }$	$20.10 ^{+ 0.15 }_{- 0.15 }$	25
2.5 < σ < 3.5	2.23	$19.94 ^{+ 0.20 }_{- 0.20 }$	$-1.16 ^{+ 0.15 }_{- 0.15 }$	$20.22 ^{+ 0.15 }_{- 0.10 }$	27
2.5 < σ < 3.5	2.81	$19.69 ^{+ 0.25 }_{- 0.15 }$	$-1.19 ^{+ 0.15 }_{- 0.10 }$	$19.97 ^{+ 0.15 }_{- 0.15 }$	27
>3.5σ	1.49	$19.86 ^{+ 0.25 }_{- 0.20 }$	$-0.89 ^{+ 0.10 }_{- 0.10 }$	$19.71 ^{+ 0.15 }_{- 0.15 }$	18
>3.5σ	1.88	$20.19 ^{+ 0.20 }_{- 0.20 }$	$-1.01 ^{+ 0.10 }_{- 0.15 }$	$20.21 ^{+ 0.15 }_{- 0.15 }$	19
>3.5σ	2.49	$20.21 ^{+ 0.20 }_{- 0.20 }$	$-0.95 ^{+ 0.15 }_{- 0.10 }$	$20.15 ^{+ 0.15 }_{- 0.15 }$	20

Notes. Below the horizontal line, we also show Schechter fit results for the same analysis but taking SpUDS cells Σ_SpUDS = Σ_SpUDS ± 1.6 × σ_SpUDS to estimate the background.

Download table as:  ASCIITypeset image

We make use of the Levenberg–Marquardt technique to solve the least-squares problem and to find the best solution for m*, α, and Φ* using the parameterization given in Equation (3). This method is known to be very robust and to converge even when poor initial parameters are given. However, it only finds local minima. In order to find the true global minimum and the true best Schechter fit to our data, we vary the starting parameters and choose the best fit solution. We show the Schechter fits for the CARLA fields of richness ⩾2σ in Figure 3.

Our data are deep enough to not just fit for m* and Φ* and assume a fixed α as had to be done in previous studies but to fit for all three quantities, α, m* and Φ* simultaneously. Mancone et al. (2012) shows that α = −1 fits the galaxy cluster luminosity function at 1 < z < 1.5 well and that it stays relatively constant down to z ∼ 0. We therefore also repeat the Schechter fits with fixed α = −1 (Table 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We estimate the uncertainties of the fitted parameters and the confidence regions of our fits using a Markov Chain Monte Carlo (MCMC) simulation. Figure 7 shows the 1σ and 2σ contours for α and m* and the best Schechter fit.

Zoom In Zoom Out Reset image size

We test our MCMC simulation by choosing a random set of Φ*, α, and m* and by computing Φ_i at m_i following Equation (3). We then add normally distributed uncertainties ΔΦ_i so that Φ_{i, err} = Φ_i + ΔΦ_i, and we fit a Schechter function to Φ_{i, err}. The original Schechter function Φ(Φ*, α, m*), the disturbed data points Φ_{i, err} and the Schechter fit agree very well with the 1 and 2σ confidence regions derived from the MCMC simulation and prove that the MCMC simulation provides us with a proper description of the uncertainties. Figure 3 shows the results of the Schechter fits to the CARLA clusters with Σ_CARLA > 2σ for all redshift bins and the confidence regions derived from MCMC simulations.

The evolution of the [3.6]−[4.5] color is very steep at 1.3 < z < 1.7; cluster members in our lowest redshift bin are expected to have bluer colors, i.e., closer to −0.1, than cluster members at higher redshift. Faint sources which have a larger color uncertainty would therefore be expected to pass the criterion in some cases and not pass the criterion in other cases. For the contaminating background/foreground sources, this scattering is expected to be the same at all cluster redshifts. We therefore test if this effect has a statistically significant influence on the faint end of the luminosity function in the lowest redshift bin. Figure 8 shows the normalized number difference between sources where [3.6]−[4.5] >−0.1 but [3.6]−[4.5]-σ_{[3.6]–[4.5]} < −0.1 and sources where [3.6]−[4.5] <−0.1 but [3.6]−[4.5]+σ_{[3.6]–[4.5]} > −0.1. In the following, we refer to these sources as "scatter-in sources" and "scatter-out sources,'' respectively. As the sum of the background and foreground color distribution is not expected to be dependent on cluster redshift and will therefore be the same for all redshift bins, any evolution in the difference between "scatter-in" and "scatter-out" sources will be caused by the cluster members. As Figure 8 shows, the

evolution of the difference is consistent with being flat with redshift. The mean of the distribution is 0.06 (N_fields)⁻¹ and a Spearman rank correlation analysis only gives a 33% chance that there is an evolution of the difference with redshift. This test confirms that the color selection is very efficient and stable with redshift and that statistically the same portion of cluster members is selected at all redshifts.

We use the ISCS (Eisenhardt et al. 2008) to validate the methodology used in this paper and to compare to previously obtained results. Mancone et al. (2010) uses the ISCS to derive luminosity functions, though he restricted the analysis to candidate cluster members based on photometric redshifts, computed as in Brodwin et al. (2006) and with deeper photometry from the Spitzer Deep, Wide-Field Survey (SDWFS; Ashby et al. 2009). We measure the luminosity function around ISCS clusters at 1.3 < z < 1.6 and 1.6 < z < 2, i.e., overlapping in redshift with the CARLA cluster sample. We again first apply the IRAC color criterion to the sources in a cell with radius r = 1 arcmin around the cluster center to determine cluster member candidates and carry out a background subtraction using SpUDS as described above. We then fit a Schechter function to the resulting data with α fixed to −1. The SDWFS data is shallower than CARLA, and we can therefore only measure the luminosity function down to [4.5] =21.4. The fitted values for m* at 4.5 μm are $20.25^{+0.30}_{-0.30}$ and $20.57^{+0.35}_{-0.35}$ for clusters at 1.3 < z < 1.6 and 1.6 < z < 2, respectively. The m* for the same redshift bins and α = −1 found by Mancone et al. (2010) are $20.48^{+0.12}_{-0.09}$ and $20.71^{+0.18}_{-0.12}$. Our results agree very well with these values, though please note the larger error bars inherent to our simple color selection compared to the photometric redshift selection of Mancone et al. (2010), which minimizes background contamination by incorporating extensive multi-wavelength supporting data. This test confirms that the method used in this paper, i.e., color-selecting galaxy cluster members and carrying out a statistical background subtraction, provides robust results. With this test, we also confirm the trend toward fainter magnitudes for m* in ISCS clusters in the highest redshift bins reported by Mancone et al. (2010).

Zoom In Zoom Out Reset image size

Figure 9 shows the evolution of m* with redshift in the context of passively evolving stellar population models. We also include results from Mancone et al. (2010). The measured m* values are compared to previous work and to model predictions of passively evolving stellar populations from Bruzual & Charlot (2003) normalized to match the observed m* of galaxy clusters at z ∼ 0.82, [4.5]$^{*}= 19.82^{+0.08}_{-0.08}$ (AB) for α = −1 (Mancone et al. 2012). Note that the results and implications of our work are independent of the model normalization. Although Mancone et al. (2010) used a normalization obtained from lower redshift cluster analyses, we use a normalization at relatively high redshift to match the redshift of the CARLA data as best as possible.

Zoom In Zoom Out Reset image size

Mancone et al. (2010) carried out a similar study and measured the IRAC1 and IRAC2 luminosity function for 296 galaxy clusters from SDWFS over 0.3 < z < 2. The clusters were identified as peaks in wavelet maps generated in narrow bins of photometric redshift probabilities. At z < 1.3, the evolution of m* agreed well with predictions from passively evolving stellar populations and no mass assembly, with an inferred formation redshift of z_f = 2.4. For the two highest redshift bins at z > 1.3, the results disagreed significantly from the continuation of the passive evolution model. Mancone et al. (2010) interpreted this deviation as possible ongoing mass assembly at these epochs. The results here for the CARLA clusters do not agree with nor continue the trend Mancone et al. (2010) report at z > 1.3. In the redshift bins in which our study and Mancone et al. (2010) overlap, we find m* to be ∼0.5 mag brighter and to continue the trend found at lower redshifts by Mancone et al. (2010). At lower redshift (1.3 < z < 1.8), the models do not differ much and m* is therefore consistent with a range of formation redshifts. At higher redshift, the models diverge and show that clusters may have formed early with formation redshifts in the range 3 < z_f < 4.

At all redshifts, our results are consistent with passive evolution models with 3 < z_f < 4. As CARLA clusters at z ∼ 3 will not necessarily be progenitors of CARLA clusters at z ∼ 1.5, it does not necessarily imply that they will remain passive subsequently. However, even at the highest redshifts probed here (z ∼ 3), we do not measure a prominent deviation from passive evolution models nor do we see signs of significant mass assembly.

We also study the evolution of m* for two CARLA subsamples of different richnesses (2.5σ < Σ_CARLA < 3.5σ and Σ_CARLA > 3.5σ) and plot the results in Figure 9. Although it is not fully known how the richness of IRAC-selected sources scales with mass, we assume here that statistically these two subsamples will represent clusters of increasing mass. For the lower richness sample, the evolution of m* appears to be similar to the full sample with formation redshifts of 3 < z_f < 4. Excluding the highest and lowest richness clusters therefore does not seem to have a significant effect on the luminosity functions. At lower redshift (i.e., 1.3 < z < 2.0) the result for m* and free α is consistent with results from Mancone et al. (2010) at z ⩽ 1.3 which might suggest that the CARLA lower richness clusters at this redshift are similar to the ISCS clusters studied Mancone et al. (2010).

Only considering the highest richness sample of CARLA clusters results in a monotonically increasing m*, with m* ∼ 0.7 mag fainter in the highest redshift bin compared to the lower richness samples. This might imply that the high richness sample includes slightly older stellar populations. This subsample is thus of particular interest for future follow-up as this population of rich, high-redshift clusters could provide a powerful probe to study the early formation of massive galaxies in the richest environments and—from a cosmological point of view—test the predictions of ΛCDM (Brodwin et al. 2010).

As can be seen in Figures 7 and 9, the Schechter fit results for α as a free parameter and α fixed to −1 generally agree well within the 1σ uncertainties. For the few exceptions, the fits do agree within the 2σ uncertainties. We conclude that a fixed α of −1 describes the luminosity function functions well at all redshifts and all densities probed in this paper. With deeper IRAC observations (1400 s of exposure time), albeit not as deep as the CARLA observations (2000 s of exposure time) and deep multi-wavelength supporting data, Mancone et al. (2012) study a subset of the original Böotes sample and measure the faint end slope α of the mid-IR galaxy cluster luminosity function to be ∼ − 1 at 1 < z < 1.5. Similar studies at lower redshift measure similar slopes. This suggests that the shape of the cluster luminosity function is mainly in place at z = 1.3.

The luminosity functions measured in this paper are consistent with α = −1 at all redshifts and all richnesses probed. Combined with results from Mancone et al. (2010) and Mancone et al. (2012), this result suggests that galaxy clusters studied in this paper have already started to assemble low-mass galaxies at early epochs. Further processes that govern the build up of the cluster and that are discussed in more detail in Section 6 then have probably no net effect on the shape of the luminosity function. This is consistent with the results found by Mancone et al. (2010).

The above results suggest an early formation epoch for galaxy clusters that are passively evolving and an early build-up of the low-mass galaxy population. These measurements, however, seem to be at odds with results investigating lower redshift clusters. Thomas et al. (2010) shows that the age distribution in high-density environments is bimodal with a strong peak at old ages and a secondary peak composed of young, ∼2.5 Gyr old galaxies. This secondary peak contains about ∼10% of the objects. Similarly, Nelan et al. (2005) derives a mean age of low-mass objects of low redshift galaxy clusters to be about 4 Gyr in low-redshift clusters. Although their observation suggests a decline of star-forming galaxies and a trend of downsized galaxy formation, low mass galaxies are still assembling until relatively recent times. Measurements of the star-formation activity for higher redshift cluster galaxies also provide evidence for continuous star-formation activity, albeit evolving more rapidly than the star-formation activity in field galaxies (Alberts et al. 2014, e.g.,).

We therefore investigate the extent to which our results can be explained by a sum of various stellar populations to estimate the maximum fraction of a star-forming cluster population that is still consistent with the data. We divide the cluster population into two stellar populations, a simple stellar population (SSP) with a delta-burst of star formation at high redshift and a composite stellar population (CSP) with a continuous, only slowly decaying star-formation rate (SFR) of the form ∝exp(− t/τ) and large τ. The details of the three different model sets we derive are as follows.

• 1.  
  Model 1. Sum of SSP with z_f = 3 and CSP with z_f = 3 and τ = 10 Gyr with ratios (SSP:CSP) ranging between 90:10 and 0:100¹² by mass.
• 2.  
  Model 2. Sum of SSP with z_f = 3 and CSP with z_f = 3 and τ = 1 Gyr with ratios (SSP:CSP) ranging between 90:10 and 0:100 by mass.
• 3.  
  Model 3. Sum of SSP with z_f = 5 and CSP with z_f = 3 and τ = 10 Gyr with ratios (SSP:CSP) ranging between 90:10 and 0:100 by mass.

A prolonged mass assembly means that at high redshift the observed 4.5 μm magnitude of galaxies is fainter because the stellar population is still forming. Therefore, accounting for mass assembly, i.e., allowing for star-forming population to contribute to the observed m*, causes m* to become fainter at high redshift depending on the contribution of this star-forming population.

In Figure 10, we show these models in the context of our results. They allow us to set an upper limit on the star-forming fraction, P, in our candidate cluster sample. Model 1, which allows for a significant star formation that is only very slowly decaying with cosmic time, shows that at all redshifts the mass fraction of the star-forming population cannot be larger than 40% (with one outlying exception of up to 60% at z ∼ 1.7). In Model 2, the SFR of the CSP decays faster and the contribution of the SSP becomes dominant much earlier than in Model 1. It therefore predicts m* to be brighter at high redshifts and to resemble the prediction of the passive evolution model. Consequently, our empirical results are in agreement with large contributions of up to 90% of the CSP in Model 2, with the CSP starting to passively evolve ∼2.3 Gyr earlier (at z ∼ 1.7) than in Model 1 (at z ∼ 0.9). Model 3 shows the evolution of a mixed population with a delta burst of star formation at z_f = 5 (SSP) and a long burst of star formation at z_f = 3 (CSP with τ = 10 Gyr). The SSP with z_f = 5 is fainter in IRAC2 at 1.5 < z < 3 than an SSP with z_f = 3, so that an additional starburst at z = 3 leads to an even fainter m* at 1 < z < 3 for Model 3. Model 3 does not reproduce the results of the LF analysis, and therefore, this scenario can be ruled out by the data.

Zoom In Zoom Out Reset image size

This shows that the results for the evolution of m* obtained in this analysis—although consistent with passive evolution models—also allow for a limited contribution of a star-forming population in galaxy clusters. Our models show that this contribution is small (up to 40% by mass, but probably on average around ∼20%) for a population with a high and slowly decaying SFR, or that this contribution is large (up to 80%) for a population with a fast decaying SFR and an evolution that resembles passive evolution ∼2.3 Gyr earlier.

For a sample of 10 rich clusters at 0.86 < z < 1.34, van der Burg et al. (2013) compares the contributions of quiescent and star-forming populations to the total mass function. We integrate the published mass functions for the quiescent and star-forming population (van der Burg et al. 2013) over galaxy masses with 10^10.1 < M_*/M_☉ < 10^11.5. We find a mass fraction of the quiescent population of ∼80% compared to the total stellar mass of the clusters. This is also in agreement with the upper limit for the fraction of the star-forming population derived for CARLA clusters in the lowest redshift bin.

As CARLA clusters at z ∼ 3 will not necessarily be the progenitors of CARLA clusters at z ∼ 1.5, we unfortunately cannot constrain the evolution of the maximum fraction of a star-forming population and cannot make conclusions about the quenching timescales and processes.

Our analysis shows that the evolution of the CARLA clusters seems to be significantly different from the cluster sample analyzed by Mancone et al. (2010). Using Spitzer 24 μm imaging for the same cluster sample, Brodwin et al. (2013) analyzed the obscured star formation as a function of redshift, stellar mass, and clustercentric radius. They find that the transition period between the era where cluster galaxies are significantly quenched relative to the field and the era where the SFR is similar to that of field galaxies occurs at z ∼ 1.4. Combining these measurements with other independent results on that sample (Snyder et al. 2012; Martini et al. 2013; Alberts et al. 2014), the authors conclude that major mergers contribute significantly to the observed star formation history and that merger-fueled AGN feedback may be responsible for the rapid truncation between z = 1.5 and z = 1.

While the ISCS clusters studied in Mancone et al. (2010) and Brodwin et al. (2013) were selected from a field survey as three-dimensional overdensities using a photometric redshift wavelet analysis, the clusters studied here are found in the vicinity of RLAGN. With RLAGN belonging to the most massive galaxies in the universe (m ∼ 10^11.5 M_☉; Seymour et al. 2007; De Breuck et al. 2010), these clusters could reside in the largest dark matter halos, deepest potential wells, and densest environments.

Indeed, Mandelbaum et al. (2009) derives halo masses for 5700 RLAGN from the Data Release 4 of the Sloan Digital Sky Survey and finds the halo masses of these RLAGN to be about twice as massive as those of control galaxies of the same stellar mass. Previous work (e.g., Best et al. 2005)has shown that more massive black holes seem to trigger radio jets more easily, but as this boost in halo mass is independent of radio luminosity, the authors conclude that the larger-scale environment of the RLAGN must play a crucial role for the RLAGN phenomenon. Similarly, albeit at higher redshift, N. A. Hatch et al. (in preparation) finds that the environments of CARLA RLAGN are significantly denser than similarly massive quiescent galaxies. They detect a weak positive correlation between the black-hole mass and the environmental density on Mpc-scales, suggesting that even at high redshift the growth of the black hole is also linked to collapse of the surrounding cluster.

This peculiar interplay between radio jet triggering, stellar mass, black hole, and halo mass of the RLAGN and the larger-scale environment suggest that (proto-)clusters and the large-scale environments of RLAGN are distinct from clusters found in field surveys.

If mergers are significantly contributing to the transition of clusters from unquenched to quenched systems, as suggested in Brodwin et al. (2013), this transition redshift will be dependent on cluster halo mass. If the environments and dark matter halos around RLAGN are indeed more massive this would explain the conflict of our results with those from Mancone et al. (2010). CARLA clusters have probably undergone this transition period much earlier than the ISCS clusters. At 1.4 < z < 1.8 where the star-forming fraction of ISCS clusters analyzed by Mancone et al. (2010) still seems to be very high (∼80%), this contribution is already much smaller in CARLA clusters. As mentioned earlier, our measurements allow us to derive upper limits on the contribution of a star-forming population but do not allow for a definite constraint on the transition redshift.