1 Introduction

Reconstructing the Near-IR Background Fluctuations from known Galaxy Populations using Multiband Measurements of Luminosity Functions

Abstract

We model fluctuations in the Cosmic Infrared Background (CIB) arising from known galaxy populations using measured UV, optical and NIR luminosity functions (LF) from a variety of surveys spanning a wide range of redshifts. We compare best-fit Schechter parameters across the literature and find clear indication of evolution with redshift. Providing fitting formulae for the multi-band evolution of the LFs out to z5, we calculate the total emission redshifted into the near-IR bands in the observer frame and recover the observed optical and near-IR galaxy counts to a good accuracy. Our empirical approach, in conjunction with a halo model describing the clustering of galaxies, allows us to compute the fluctuations of the unresolved CIB and compare the models to current measurements. We find that fluctuations from known galaxy populations are unable to account for the large scale CIB clustering signal seen by Spitzer/IRAC and AKARI/IRC and continue to diverge out to larger angular scales. This holds true even if the LFs are extrapolated out to faint magnitudes with a steep faint-end slope all the way to =8. We also show that removing resolved sources to progressively fainter magnitude limits, isolates CIB fluctuations to increasingly higher redshifts. Our empirical approach suggests that known galaxy populations are not responsible for the bulk of the fluctuation signal seen in the measurements and favors a very faint population of highly clustered sources.

Subject headings:
cosmology: diffuse radiation — large-scale structure of universe — galaxies: evolution — luminosity function — infrared radiation

1. Introduction

Cosmic infrared background (CIB) includes contributions from emissions over the entire history of the Universe, including from objects inaccessible to the current telescopic studies. Several direct measurements of the total mean levels of the CIB using the wide-beam DIRBE and IRTS instruments claim a significant excess mean flux over the contribution of known galaxies in the near-IR (Dwek & Arendt, 1998; Gorjian et al., 2000; Wright & Reese, 2000; Cambrésy et al., 2001; Matsumoto et al., 2005); also see review by Kashlinsky (2005). The entire excess emission over that from known galaxy populations ( 30  in 1-4, Kashlinsky (2005)) was originally theorized to come from primordial PopIII stars (Santos et al., 2002; Salvaterra & Ferrara, 2003) but this interpretation has been challenged on several grounds since the claimed levels require uncomfortable levels of star formation efficiency (Madau & Silk, 2005; Salvaterra & Ferrara, 2006). It is possible that much of the excess flux seen may be due to inaccurate removal of bright zodiacal emission in the foreground (Dwek et al., 2005; Mattila, 2006). Furthermore, the extragalactic background light (EBL) is a fundamental source of opacity for high energy photons and the -ray attenuation seen in blazar spectra favors low levels of NIR background light (e.g. Aharonian et al., 2006; Mazin & Raue, 2007).

An alternative way to study the CIB, much less sensitive to foreground removal, is to measure background anisotropies after subtracting resolved galaxies down to faint magnitudes (Kashlinsky et al., 1996). Detections of spatial structure in the CIB were initially based on datasets from COBE/DIRBE (Kashlinsky & Odenwald, 2000), the IRTS (Matsumoto et al., 2000) and 2MASS (Kashlinsky et al., 2002; Odenwald et al., 2003). More recently, Kashlinsky et al. (2005, 2007b) using deep exposures from Spitzer/IRAC (3.6-8.0) found significant fluctuations after subtracting galaxies down to 25. The level of these fluctuations, 0.1  at arcminute scales, imply an isotropic CIB flux as low as 1  from the remaining unresolved sources in the IRAC bands (Kashlinsky et al., 2007c). Thompson et al. (2007a) analysis constrained CIB at 1.4-1.8 using HST/NICMOS observations and Matsumoto et al. (2011) measure fluctuations on arcminute scales in the - range using the AKARI satellite. After this paper was submitted, Kashlinsky et al. (2012) have measured the Spitzer/IRAC out to 1 using more extensive datasets from the Spitzer Extended Deep Survey (SEDS), confirming earlier results and extending the fluctuation measurement to much larger angular scales. All the present measurements of CIB-fluctuations are consistent with an extragalactic origin, necessitating an associated unresolved component in the CIB. This component likely requires only a fraction of the CIB excess, which is below limits imposed by -ray photon absorption (Kashlinsky et al., 2007c; Kashlinsky & Band, 2007; Arendt et al., 2010).

There seems to be an emerging consensus that the extragalactic clustering signal is real, but the nature of the sources producing it is still a subject of debate. Plausible candidates for the bulk of the CIB are evolving stellar populations in galaxies, although accreting black holes at high- can also contribute (e.g., Ricotti & Ostriker, 2004). Both Kashlinsky et al. (2005) and Matsumoto et al. (2011) argue that that the clustering is consistent with ”first stars” era objects whereas Cooray et al. (2007); Chary et al. (2008) have posited that the signal originates mostly in the clustering of faint galaxies at redshifts 1-3. Understanding the expected levels of fluctuations from known galactic populations is possible following the establishment of the standard cosmological model for structure formation, the concordance CDM (Komatsu et al., 2011). In order to compute the levels of source-subtracted CIB fluctuations remaining in the Spitzer data, Sullivan et al. (2007) used a halo model combined with conditional luminosity functions and compared it to measurements at 3.6. Their claim is that the fluctuations detected by Kashlinsky et al. (2005) can be explained by ordinary galaxies just beyond the detection threshold of Spitzer/IRAC, although this claim appears to contradict the results of their analysis shown in their Fig. 8.

Kashlinsky (2005) discusses the importance of the shape of the emission history for the resulting fluctuations demonstrating how brief episodes of light production can lead to enhanced fluctuations. In this paper, we construct the entire history of light production produced by known galaxy populations using a novel empirical approach that relies exclusively on observations. We use a compilation of galaxy luminosity functions (LF) in the literature to populate the observed lightcone with galaxies down to faint magnitudes. The many galaxy surveys conducted in recent years provide a wealth of data in multiple bands and cover a wide range of redshifts. Individually, LFs only probe specific rest-frame wavelengths for a limited range of redshifts, while together we can use them to infer the source distribution composing the background light in the 0.1-5.0 range. Our only theoretical assumptions concern the clustering properties of the unresolved sources which are modeled according to the well-established concordance CDM model (see Section 5). We refer to Johnston (2011) for a good review on the properties of luminosity functions and how they are measured.

Modeling the underlying populations of the EBL has been attempted using various mixtures of theory and observations. Backward evolution scenarios take the present galaxy populations and extrapolate them to higher redshift (e.g., Jimenez & Kashlinsky, 1999; Franceschini et al., 2008), while forward evolution follows dark matter merger trees starting from the cosmological initial conditions, using semi-analytical models of galaxy formation (Gilmore et al., 2010; Guo et al., 2011). Domínguez et al. (2011) use directly the measured K-band LFs out to z=4 from Cirasuolo et al. (2010), combined with best-fit SEDs of multi-wavelength galaxy data (AEGIS) to empirically derive the overall EBL spectrum. We however, present an alternative empirical approach by examining the best-fit Schechter parameters (Schechter, 1976) of LFs covering the UV, optical and near-IR out to redshifts z3-8. We provide empirical fitting formulae describing the smooth evolution of multi-band LFs with redshift, and construct lightcones containing all populations seen in the near-IR bands, selected at each redshift such that .

This paper is organized as follows. Section 2 describes the data used and in Section 3 we explain the modeling in detail. In Section 4, we calculate both galaxy number counts and the EBL in the near-IR bands () and compare with existing data. In Section 5 we analyze the source-subtracted CIB-fluctuations implied by our empirical reconstruction and compare with previous work. We discuss the implications of our findings in Section 6. Throughout this paper we adopt the concordance CDM comology with =0.3, =0.7 and =70 . All magnitudes are in the AB system unless stated otherwise (Oke & Gunn, 1983).

2. Measurements of the Galaxy Luminosity Function

Figure 1.— All luminosity functions used in our analysis in Schechter parametrization (see original references in Table 1). The wavelength bins are listed in the panels (lower right) and their effective wavelengths are listed in Table 2 along with other properties. The LFs shown have a range of redshifts.
Reference Rest-frame band Redshift Sample Selection Survey Catalog / Field Symbol / Color


Arnouts et al. (2005)
1500Å 0.2-1.2 1039 NUV24.5 GALEX/VVDS green triangles(up)
1.75-3.4 F450&F60627 HDF
Wyder et al. (2005) 0.055 896,1124 GALEX/2dF blue circles
Oesch et al. (2010) 1500Å 0.5-2.5 284-403 26 HST ERS yellow circles
Oesch et al. (2012) 1500Å 8 70 27.5 CANDLES/HUDF09/ERS pink triangles(up)
Reddy et al. (2008) 1700Å 1.9-3.4 15,000 25.5 blue crosses
Yoshida et al. (2006) 1500Å 4,5 3808,539 26-27 Subaru Deep Field blue squares
McLure et al. (2009) 1500Å 5,6 1500 26 SXDS/UKIDSS purple squres
Ouchi et al. (2009) 1500Å 7 22 26 SDF/GOODS-N
Bouwens et al. (2007) 1600Å,1350Å 4,5,6 4671,1416,627 29 HUDF/GOODS violet triangles(down)
Bouwens et al. (2011) 1600Å,1750Å 7,8 73,59 26-29.4 HUDF09 orange diamonds

Gabasch et al. (2004)
0.45-5 5558 FORS Deep Field Green triangles(down)
Baldry et al. (2005) 0.3 43223 u20.5 SDSS red squares
Faber et al. (2007) 0.2-1.2 34000 DEEP2/COMBO-17 tan squares
Norberg et al. (2002) 0.2 110500 19.45 2dFGRS purple squares
Blanton et al. (2003b) 0.1 147986 16.5-18.3 SDSS blue plus
Montero-Dorta & Prada (2009) 0.2 947053 17-19 SDSS green crosses
Loveday et al. (2012) 0.002-0.5 8647-12860 GAMA yellow squares

Ilbert et al. (2005)
0.05-2.0 11034 VIMOS-VLT Deep Survey pink triangles(up)
Gabasch et al. (2006) 0.45-3.8 5558 FDF green circles
Marchesini et al. (2007) 2.0-3.5 989 MUSYC/FIRES/GOODS/EIS orange circles
Marchesini et al. (2012) 0.4-4.0 19403 27.8,25.6 blue triangles(up)

Hill et al. (2010)
0.0033-0.1 2437-3267 18-21 MGC/UKIDSS/SDSS purple diamonds
1589-1798 17.5-18
Dahlen et al. (2005) 0.1-2 18381 GOODS-HST/CTIO/ESO dark green diamonds
0.1-1 2768
Jones et al. (2006) 0.2 138226 6dFGS/2MASS dark red plus
/SuperCOSMOS
Bell et al. (2003) 22679 SDSS orange circles
6282 2MASS
Kashikawa et al. (2003) 0.6-3.5 439 Subaru Deep Survey red crosses
Stefanon & Marchesini (2011) 1.5-3.5 3496 - MUSYC/FIRES/FIREWORKS green squares
Pozzetti et al. (2003) 0.2-1.3 489 K20 Survey tan plus
Feulner et al. (2003) 0.1-0.6 500 - MUNICS yellow crosses
Eke et al. (2005) 0.01-0.12 16922,15664 2dFGRS/2MASS violet diamonds
Cole et al. (2001) 0.005-0.2 7081,5683 2dFGRS/2MASS blue squares
Smith et al. (2009) 0.01-0.3 40111 , UKIDSS-LAS/SDSS red triangle(up)
Saracco et al. (2006) 0.001-4 285 HDFS/FIRES blue triangles(down)
Kochanek et al. (2001) 0.003-0.03 4192 2MASS/CfA2/UZC magenta circles
Huang et al. (2003) 0.001-0.57 1056 2dF/AAO violet diamonds
Arnouts et al. (2007) 0.2-2 21200 SWIRE/VVDS dark green squares
/UKIDSS/CFHTLS
Cirasuolo et al. (2010) 0.2-4 50000 UKIDSS/SXDS orange plus
Babbedge et al. (2006) 0.01-0.6 34281 SWIRE/INT WFS blue crosses
Dai et al. (2009) 0.01-0.6 4905,5847 IRAC-SS/AGES dark red circles

Note. – The measured LF are shown in Figure 1 and all Schechter parameters are displayed in Figure 3. The compile database of the Schechter parameters is available upon request.

Data taken from multiple surveys/fields

The symbols and color of the corresponding data points in Figure 2 and 3

Table 1Measurements of the Luminosity Function

The total emission seen in the near-IR bands () depends on the contribution of local near-IR galaxies as well as redshifted light radiated at shorter rest-frame wavelengths. To quantify the present day background produced by galaxies, we have utilized measurements of luminosity functions probing all rest-frame wavelengths in the interval 0.15.0 anywhere in the redshift cone. This results in a compilation of LFs from a large variety of surveys which we list in Table 1. Our approach does not depend on stellar population synthesis models (e.g., Bruzual & Charlot, 2003) and we do not need to make an assumption for the IMF. Rather, in this method we predict the levels of CIB fluctuations directly from the available data, assuming only i) standard CDM model of structure formation and ii), the validity of a Schechter-type LF after fitting its parameters to the data. All the LFs we use have been characterized by a Schechter function (Schechter, 1976),

(1)

determined by the normalization, , characteristic absolute magnitude, and the faint-end slope, . By integrating Equation (1), the luminosity density can be shown to be , where is the characteristic luminosity and is the Gamma function. All the Schechter LFs used are shown in Figure 1.

The Schechter LF is usually found to fit the data fairly accurately but deviations are seen, in particular when fitting a wide range of luminosities. At low- for example, Jones et al. (2006) find that the shape does not fit the sharp downturn seen at and both Blanton et al. (2003b) and Montero-Dorta & Prada (2009) find an excess of bright galaxies in the blue SDSS bands. There are also hints of an upturn in the local LF at faint magnitudes where the Schechter fit does a poor job (e.g., Blanton et al., 2005). We address this faint-end issue in Section 3.1, but note that sources at the bright end are efficiently removed from the maps in CIB fluctuations studies. At longer wavelengths (5), a double power-law is found to provide a more adequate fit than the Schechter function (Babbedge et al., 2006; Magnelli et al., 2011). The mutual consistency of measurements is a primary concern when comparing LFs in the literature. Inconsistencies can be caused by field-to-field variations, photometric system, k-corrections, type of LF-estimator, survey depth and completeness, redshift binning, sample statistics, error estimates, etc. These undoubtedly account for differences in shape and amplitude of the measured LF (see Figure 3). We include a discussion of common issues in Appendix B, but these do not affect our results because we let all measurements collectively contribute to our derived LF (see Section 3).

To directly compare flux measurements at different wavelengths, we have adopted the magnitude system which conveniently relates the apparent magnitude, , to the specific flux, , via

(2)

(1 Jy = 10WmHz). Where system conversions are not explicitly given by the authors (Table 1) we have made use of the calculations available at http://mips.as.arizona.edu/cnaw/sun.html. With all magnitudes converted to we do not distinguish between magnitudes of different filter- and photometric variations, e.g. Johnson and SDSS , apart from their center frequencies.

Figure 2.— The local luminosity density according to all available LF measurements at z0.12 in Table 1, with symbols/colors indicated in the same Table. To avoid overcrowding the region of interest we omit error bars. The solid line shows the luminosity density in our fiducial bands as implied by our fits in Figure 3. The sets of gray lines show the contribution from galaxies of different metallicities and ages from synthetic galaxy SED spectra shown in Fig. 14 of (Kashlinsky, 2005). The bottom-gray curves show the early type stellar populations, the upper-dark show late type populations and middle-light lines show the average of the two contributions.

3. Populating the Lightcone with known Galaxy Populations

This section outlines the step-by-step approach leading to the quantification of the galaxy distribution seen on the sky. Using the data in Table 1, we populate the evolving lightcone by placing the rest-frame galaxy distribution at a distance such that the associated emission is shifted into the near-IR bands in the observer frame, defined by . Initially, we bin the LFs according to their rest-frame wavelength in fiducial bands which we call , , , , , , , , , ,  and  (see Table 2). For example, measurements in rest-frame SDSS , Johnson and 2dF are binned together in our -band despite having an offset in center wavelength of about 0.03. The largest offset occurs in our -bin where the centers of SDSS and Johnson is 0.063. The uncertainty associated with the redshift of the population usually dominates these offsets so we do not correct for them. The centers of our fiducial bands, , are taken to be the mean rest-frame wavelength of all measurements in the bin (see Table 2).

Figure 3.— The measured Schechter parameters ,, from the studies in Table 1 including the luminosity density, =, as a function of redshift. The different sybmols/colors are listed along with the corresponding references in Table 1. We have omitted error bars for the sake of clarity. The solid curves show the evolutionary fits according to Equations (3)-(4) with the best-fit parameters listed in Table 2. We have modified to follow the fitting functions of Bouwens et al. (2011) at z3.5 to better match the turnover seen. We note that our fits are only empirically supported for z4, beyond which we extrapolate. The dashed curves in the panels shows the evolution assumed in our default model whereas the dark shades areas encompass the range bracketed by our high faint-end (HFE) and low faint-end (LFE) models. These ranges are ultimately constrained by the observed galaxy counts (see Section 3.1). The shaded areas in the bottom row () is the evolving quantity corresponding to this allowed range in . The dotted curves in in  and -bands are not fits to the data but are instead assumed to have the same form as the -band fits. The light gray shaded areas correspond to the redshift regions for which the rest-frame emission redshifts into the observed NIR wavelengths of interest, defined to encompass the 1.25-4.5 range. We are most concerned with the goodness of fit in these regimes. All the data-points assume =0.7.

By placing the entire population of each LF at the median redshift of the sample, , we examine the evolution of the individual Schechter parameters (,,) in our fiducial bands. In the cases where is not explicitly given by the authors, we choose the midpoint of the redshift bin of the LF measurement. The distances of the galaxies composing the LF is the dominant uncertainty in the resulting counts on the sky and we have therefore examined the effects of placing the LF at the opposite boundaries of the redshift bin (the resulting counts differ by less than a factor of two at the two extremes (see Section 4)). Figure 3 shows the Schechter parameters as a function of redshift from 0.15-4.5. Across the spectrum, we see clear indication of evolution in and and in some cases also in the poorly measured .

Over time, galaxy populations evolve both in brightness and abundance. As small systems merge to form more massive ones, we expect a net increase in the number of bright and massive galaxies with time accompanied by a decrease in fainter ones. This is encoded in the evolution of (the number density of systems), which we expect to increase with time whereas the faint-end slope, should consequently flatten. The difference of the LF among rest-frame bands reflects the tendency of galaxies of different types being preferentially bright/faint at a given wavelength. The decomposition of the LF into red/blue galaxies typically shows an early-type population of individually bright galaxies with a diminishing faint-end whereas a the late-type population is composed of a rising number of faint galaxies (e.g., Faber et al., 2007). The characteristic luminosity, therefore depends heavily on the mixture of spectral types at any given epoch. Much work has been devoted to the K-band LF where the stellar mass-to-light ratio is relatively stable and it can thus be used as an indicator for the stellar mass function (e.g., Cole et al., 2001). It is therefore natural to expect to brighten with cosmic time as more mass becomes locked up in low-mass stars. In the red/NIR bands, the luminosity evolution is typically between redshift 0.1 and 1 whereas it is much stronger in the UV/blue rest-frames indicating higher star formation rates at earlier times. Extensive work has been done on the UV LF which is largely driven by its importance as a tracer of star formation rate and assisted by increasing detection rates of distant Lyman break galaxies in deep surveys. We see brighten with increasing redshift and then turning over, thus roughly exhibiting the same behavior as the derived star formation history (Madau plot). The wide redshift range of available UV LF measurements makes it the only LF in which a non-monotonic evolution is distinctly seen in . In all other bands, the evolution of the Schecter parameters can be fitted with an analytic function to quantify the global evolution, while “washing” out outliers in the process. Several authors have parameterized the evolution in individual bands (e.g., Lin et al., 1999; Cirasuolo et al., 2010), but to our knowledge, our work is the first multi-wavelength parametric study of the evolution of the LF parameters. We find the following forms to fit the data well across our wide range of wavelengths and redshifts:

(3)
(4)

and we assume the following a priori form for the faint-end slope

(5)

These fits are shown in Figure 3. For and we have taken =0.8, but = 0.01 for . The other best-fit parameters are listed in Table 2. Instead of selecting a preferred LF measurement for a given redshift in each band we have chosen to let all measurements contribute equally to the fitting process regardless of depth, area and sample size of the survey. Although there are a few notable discrepancies between the data and the fits we note that our IR-fluctuation results are unaffected as long as the fits remain good in the light shaded areas of Figure 3. These regions correspond to the distance for which the rest-frame emission is redshifted into the observed near-IR wavelengths of interest, defined to encompass the 1.25-4.5 range. In the following sections we will rely on lightcones extrapolated from the highest measured redshift, typically z4, out to z=7 (see Table 2). To account for the turnover observed in , we only use our Equation (3) out to z3 where they intersect the high- fitting formulae given by Bouwens et al. (2011) which we adopt for z3.

Evolution is not easily discerned in the faint-end slope, , which by the very nature of surveys is hard to measure over large distances. For this reason we explore different scenarios for the behavior of which we explain in Section 3.1. In the  and  bands, the redshift range covered by the available measurements is so limited that we can only fit but not the other Schecther parameters. Thus, for these two bands we assume to take on the same form as the neighboring -band. Fortunately, the data available in the -bands covers the redshift range of interest as is indicated by the shaded regions in Figure 3.

Band , , ,
0.15 24 8.0 -19.62,1.1 2.43,0.2 -1.00,0.086
0.36 27 4.5 -20.20,1.0 5.46,0.5 -1.00,0.076
0.45 44 4.5 -21.35,0.6 3.41,0.4 -1.00,0.055
0.55 18 3.6 -22.13,0.5 2.42,0.5 -1.00,0.060
0.65 25 3.0 -22.40,0.5 2.25,0.5 -1.00,0.070
0.79 17 3.0 -22.80,0.4 2.05,0.4 -1.00,0.070
0.91 7 2.9 -22.86,0.4 2.55,0.4 -1.00,0.060
1.27 15 3.2 -23.04,0.4 2.21,0.6 -1.00,0.035
1.63 6 3.2 -23.41,0.5 1.91,0.8 -1.00,0.035
2.20 38 3.8 -22.97,0.4 2.74,0.8 -1.00,0.035
3.60 6 0.7 -22.40,0.2 3.29,0.8 -1.00,0.035
4.50 6 0.7 -21.84,0.3 3.29,0.8 -1.00,0.035

Note. – 1) Fiducial rest-frame band, (2) the effective wavelength in microns, (3) number of LFs used, (4) highest redshift of LF available in band, (5) Best-fit parameters for with =0.8, (6) Best-fit parameters for with =0.8 in units of , (7) The parameters for chosen to reflect the models (HFE&LFE) presented in Section 3.1. assumed to be the same as in

Table 2Properties of the data shown in Figure 3 and the best-fit evolution parameters of Equations (3)-(5)

There is significant degeneracy in the Schechter parameters derived for a given galaxy population which can manifest itself in different values of (,,) depending on the LF-estimator used (see Appendix B). The overall shape of the LF can appear similar despite different Schechter parameters typically resulting in a comparable value for the luminosity density, , which we display in the bottom panels in Figure 3. For example, Ilbert et al. (2005) (VVDS) and Gabasch et al. (2006) (FDF) derive comparable LFs depite giving very different values for the Schechter parameters. The general agreement of the -data and the curves, , indicates that our separate fits do not systematically over- or under-estimate the total luminosity density.

The second step is populating the lightcone seen from the standpoint of the observer. Light from distant galaxies appearing in the observed -band was emitted at wavelength i.e. at all rest-frame wavelengths shortwards of throughout the redshift cone. We extract the Schechter parameters from our fits in Figure 3 at the redshift defined by where corresponds to our fiducial bands () and . Our template LFs then become

(6)

The continuous evolution of the LF seen in the -band is then obtained by interpolating the ’s from to . It should be noted that because of the degeneracy, our separated (,,) fits used in Equation (6), cause some amount of deviation from the original shape of the LF. This is a small effect in comparison with the general disagreement between individual authors on the shape of the LF. We refer to Appendix A where an independent method is used to populate the lightcone, in which the original shapes of the LFs are kept intact. We show that the two different methods produce the same results, confirming the validity of our standard treatment.

As an example we show in Figure 4 the Schechter parameters characterizing the LFs, probing the sky in two different observer-frame bands centered at 2.2, 3.6  respectively.

Figure 4.— Evolution of the Schechter parameters and the luminosity density seen in the observed at 2.2 (light squares) and 3.6 (dark diamonds). The values are extracted from the fits in 3 at the appropriate redshifts. is in units of 10Mpc and in units of 10ergsMpc.

Although the abundance of galaxies diminishes by itself at high- according to our fits, we impose a limit of 7 in our modeling, beyond which we assume that ordinary galaxy populations were not yet established. But due to the steep drop of at high-, our results are not sensitive to this parameter: in fact, using , yields results nearly identical to our fiducial model. We emphasize that our evolution models are empirically supported out to z4 only, beyond which we extrapolate the evolution deduced at lower redshifts.

In order to deduce the rest-frame LF from survey data, absolute magnitudes need to be derived from apparent magnitudes. We derive the flux from galaxies in our lightcone by backtracking the original procedure i.e. going from absolute magnitudes back to apparent magnitudes. This implies undoing any corrections the authors have made in this process

(7)

where is the distance modulus, is the k-correction, is the evolution correction and is the correction due to galactic extinction at the Galactic coordinate . In LF measurements, authors typically use de-reddened magnitudes or correct for extinction using Galactic dust maps (Schlegel et al., 1998). This correction can be large in the UV/optical but becomes less severe towards the infrared where we have 7-10 approaching 15-20 in the IRAC bands (Cardelli et al., 1989; Indebetouw et al., 2005). We are only concerned with emission entering the Milky Way as near-IR where the extinction correction is typically well within 0.1 mag so we neglect it in Equation (7). Correcting for evolution is intended to make a sample drawn from a distribution of redshifts reflect the true luminosity function at a given epoch (usually of the survey/bin) by accounting for changes in luminosity and number density over time (e.g., Blanton et al., 2003b). This has been done for some local surveys where a considerable spread in the redshift distribution leaves more cosmic time for evolution to take place. This typically results in corrections of 0.1 mag (Bell et al., 2003) but since the evolution correction simply acts to make the LF more accurate at a given redshift we do not need to make any adjustments. The only magnitude adjustment in Equation (7) of concern is the k-correction (Hogg et al., 2002) which is needed to transform to the rest-frame by accounting for the redshifted SED of a given source. There are a variety of methods to deal with this SED dependence and we refer to Appendix B for a more complete discussion of two commonly used treatment in the literature. From the k-corrected absolute magnitudes, we simply require the spectral independent term to account for the redshift into the observed frame, . Equation (7) is now reduced to

(8)

which is the conversion we use. In Section 4 we show that we recover the observed number counts to a very good accuracy using this methodology.

3.1. The Faint-End LF Regime

The source subtracted CIB fluctuations are isolated to faint sources. By the nature of galaxy surveys, the faint-end is generally poorly constrained causing large uncertainties and scatter in measurements of , especially at high-. Because of this, many authors prefer to keep fixed in their Schechter fits. Since the data does not show robust evolution in in most bands (unlike and ) we explore variants of the behavior of the faint-end slope to get a feel for the sensitivity of CIB fluctuations to the abundance of faint galaxies. The substantial scatter in measurements of leaves us some freedom in modifying the faint-end regime but we find that deep galaxy counts impose strict limits on the allowed range of faint-end slopes. This is most notable in , where a steep faint-end at =1-3 leads to an overproduction of the observed number counts in the faintest magnitude bins (see Figure 5). We therefore consider the range of allowed scenarios that collectively yield galaxy counts consistent with observations across all bands simultaneously. We leave and unchanged when varying despite degeneracies in the parameters (see appendix A). We consider two models, high faint-end (HFE) and low faint-end (LFE), which, based on the resulting galaxy counts, are likely to bracket the true behavior of the faint-end of ordinary galaxy populations. These are shown in Figure 3 and 5 as the upper and lower boundaries of shaded regions. With the faint-end reasonably well constrained at z=0, ranging from -0.8 to -1.2, we fix at these two values for LFE and HFE respectively and vary later evolution by changing the slope of the power-law in (called in Equation (5))1. Our HFE model is characteristic of strong steepening such as that found by Ilbert et al. (2005) (VVDS) out to z1 whereas the LFE implies a more modest evolution, closer to that of Marchesini et al. (2007, 2012). Our LFE reflects a lack of evolution in the NIR i.e. const., which seems to be favored by some authors (e.g., Cirasuolo et al., 2007). We choose a faint-end cutoff for each template LF at for LFE and for HFE, thereby extrapolating the LF to very low luminosities. For both scenarios we find to be near saturation with flux contribution for fainter magnitude bins always being 0.02 . Our “default” model is the average of HFE and LFE with a cutoff at .

We have chosen our LFE/HFE models so that they remain consistent with number counts data. The LFs dominating the faint counts in Figure 5 are mostly determined by the faint-end slope, , at high and intermediate redshifts and it is important to emphasize that more extreme faint-end evolution models generally yield number counts that are inconsistent with observations. Alternatively, one could in principle imagine an increase in the LF in the faintest magnitudes observed deviating from a Schechter function. In fact, such an upturn has been observed locally, for which a “double” Schechter function provides a better overall fit of the LF (e.g. Blanton et al., 2005; Loveday et al., 2012). Allowing for a much steeper slope at z=0 to accommodate this possibility does not affect the resulting CIB fluctuations because the surface density of sources on the sky tends to be dominated by populations at larger distances. This can be illustrated by examining the underlying LFs of the resulting galaxy number counts in Figure 5, where the gray lines starting at the bright-end (from left) correspond to the local contribution (the thick line being the most local) moving to high- LFs to the right. The rapid redshift evolution of the cosmic volume element prevents a large surface density of low- sources and we find the faint counts always being dominated by populations at intermediate and high redshifts (z1). In order for low- sources to have sufficient densities to dominate the faint galaxy counts, and thereby also the unresolved fluctuations, we would need an extremely steep faint-end at z=0, becoming flatter towards increasing redshift i.e. which is the opposite of the observed evolution trend. Alternatively, a sudden upturn deviating from a Schechter form at faint magnitudes would need to shoot up by roughly two orders of magnitude. We therefore consider our HFE scenario to be sufficiently extreme at low- and making it steeper does not have an effect on our results. On the other hand, if a significant upturn in the LF exists at z0.5 (so far undetected), then this may result in a non-negligible contribution to the unresolved fluctuations. The large number of small halos predicted by the standard CDM model permits such a scenario, especially if the first population of dwarfs with normal stellar populations formed in halos with mass 10 M (Ricotti et al., 2002a, b, 2008). For instance, if the ultra-faint dwarf galaxies recently discovered around the Milky Way can be identified as fossils of the first galaxies formed before reionization, that would imply that we have only discovered a small fraction of a widespread population of dwarfs which were almost certainly brighter in the past (Ricotti & Gnedin, 2005; Bovill & Ricotti, 2009, 2011a, 2011b). However, it is unclear how to make the flux from this population sufficiently large to reproduce the measured fluctuation signal and, furthermore satellite dwarfs are efficiently masked along with their host galaxy in fluctuation measurements as displayed by the masking typically having angular radius of 15 (Arendt et al. (2010), see also Fig. A-3 in Kashlinsky et al. (2012)). In this work we probe whether the known galaxy populations, which we extrapolate to faint magnitudes in our HFE and LFE limits, can account for the observed source-subtracted CIB fluctuations, and the question of the nature of the new populations that can explain these fluctuations is, while important, outside the scope of the current discussion.

4. Number Counts and Background Light from LF Data

Figure 5.— Galaxy number counts in our default description (solid curve) including the regions of bracketed by our two extreme models, HFE and LFE (gray shaded areas). The gray curves show the underlying template LFs in our fiducial bands (Equation (6)) which we interpolate and integrate to obtain the number counts via Equation (9). The low- LF dominate the bright counts whereas high- and intermediate redshift LFs dominate the faint counts (from left to right). The most local available LF is shown as thick gray curves to demonstrate their negligible contribution to the faint counts. For 0.45-0.80 panels the data are from Capak et al. (2004) (red asterisks),Capak et al. (2007) (purple diamonds), McCracken et al. (2003) (green triangles), Yasuda et al. (2001) (turqoise diamonds) and Kashikawa et al. (2004) (blue squares). Data in the 1.25-2.2 panels are taken from Väisänen et al. (2000) (green triangles), Dickinson et al. 1999 (purple squares), Maihara et al. (2001) (green asterisks), Keenan et al. (2010b) (blue triangles), Keenan et al. (2010a) (blue diamonds) Frith et al. (2006) (yellow triangles), Thompson et al. (2005) (yellow asterisks), Metcalfe et al. (2006) (green diamonds), Quadri et al. (2007) (turqiose triangles), Baker et al. (2003) (purple squares), Minowa et al. (2005) (orange squares), Huang et al. (1997) (red crosses) and the 3.6-4.5 data comes from Fazio et al. (2004) (purple symbols).

Galaxy number counts have the advantage of being free of the uncertainties associated with e.g. k-corrections and redshift determination which makes it an important test of both cosmology and galaxy evolution models. We project our lightcones onto the sky to obtain the galaxy number counts in each magnitude bin per unit solid angle:

(9)

where is the comoving volume element per solid angle. In Figure 5 we display the number counts from Equation (9) in the 0.45-4.5 range and compare with existing data in the literature. The good agreement between our modeling and observed counts demonstrates the validity of our method. We also display the range bracketed by out two limiting models for the faint-end slope of the LF, as discussed in Section 3.1 (shaded areas). The gray curves in Figure 5 reflect the underlying template LFs contributing to the number counts in different redshift bins (bright/left to faint/right correspond roughly to low- to high-), elucidating the different populations governing the source surface density on the sky. It is reassuring, although not surprising, that we recover the shape of the galaxy counts using independent observations (the only assumption being the Schechter parametrization of the LF). This explicitly confirms that our multi-wavelength collection of observed LFs provides an accurate description of the photometric properties of resolved galaxies on the sky.

Figure 6.— Left: Our reconstructed number counts in BRIJKHLM compared across the spectrum. The counts have been multiplied by a slope of 10 to bring out features in the shape. In this representation is proportional to the flux contribution from each magnitude bin. Right: The accumulation of integrated background light from galaxies over time. The flux builds-up from high- (right) to low- (left) reaching the present-day observed value listed in Table 3

Figure 6 (left) examines how the shape of the number counts varies across the spectrum (0.45-4.5). Both the shape and amplitude of the counts are governed by the behavior of the (,,)-parameters shown in Figure 4 and some particular features deserve a few remarks. The bright counts all start out with a well known (Euclidian) slope of 0.6 continuing down to 18-20 where it flattens to 0.4. To first order this “knee” is simply caused by the transition from -dominated to -dominated regime. More specifically, a dip appears in the number counts at 18-20 which arises from the lack of very bright galaxies at higher redshifts, i.e. becomes fainter with redshift (see Fig. 4). At higher redshifts (and shorter rest-frame wavelengths) we see a brightening again which is associated with star forming galaxies, bright in UV rest-frames. This brightening causes another feature at  mag revealing a “double-knee” surrounding the dip. This is most pronounced in the -counts but disappears at longer observed wavelengths where the UV rest-frame becomes too distant. Beyond 25-26, the counts are gradually diminished by the CDM volume element. Depending on the exact faint-end model, the logarithmic slope in this regime is 0.2-0.3 in , decreasing as we go to longer wavelengths.

Another clear feature of the number counts seen in Figure 6 (left), is the overall increase per magnitude bin as we go to longer observed wavelengths. We find the reasons for this to be twofold. First, the bright end is typically dominated by galaxies which are more luminous in the red bands such as the case of giant ellipticals. Therefore we see a larger number of them out to greater distances (in Fig. 3 we clearly see becoming overall brighter from blue to red). Second, when we look at the Universe through redder bands, we observe the redshifted light from bluer rest-frames emitted in epochs when the star formation activity was greater and consequently was brighter. We further point out that our reconstructed counts are immune to confusion and agree well with the confusion corrected Spitzer/IRAC counts of Fazio et al. (2004) (confusion enters around 20-22).

We infer the amount of background light from galaxies from our reconstructed counts:

(10)

where of Equation (2) and is the integrated flux in units of . Figure 6 (right) shows how extragalactic background light builds up with cosmic time observed through . This results in present day values of the integrated background light of 9.6, 9.3, 8.1, 4.9 and 3.3  at 1.25, 1.63, 2.2, 3.6 and 4.5 respectively (see Table 3), which agree very well with Table 5 of Kashlinsky (2005) and are also in general agreement with Madau & Pozzetti (2000), but slightly lower than the values found by Keenan et al. (2010a). A subtle underestimation could be due to the smooth fitting of the LF evolution which smears out any abrupt variation of the Schechter parameters which could either be physical. The small deficit with respect to the EBL of Keenan et al. (2010a) arises in the 21-23 mag range where we see a better agreement with Madau & Pozzetti (2000) and Maihara et al. (2001).

5. Near-IR Fluctuations from Unresolved Galaxies

We now turn to evaluating the source-subtracted CIB fluctuations keeping in mind the procedure leading to their detection from raw images. If enough pixels remain in the maps after the masking of resolved sources, the fluctuations can be characterized via their angular power spectrum, which can then be computed more efficiently by using FFTs than the 2-point correlation function. For a detailed description of the process of reducing CIB fluctuation data in the Spitzer/IRAC analysis we refer to Arendt et al. (2010).

Band
All
3.33 2.26 1.17 0.52 4.92
2.95 1.90 0.96 0.42 5.65
2.86 1.75 0.85 0.37 6.56
2.81 1.58 0.72 0.30 7.97
2.59 1.20 0.48 0.18 9.60
2.25 0.96 0.36 0.13 9.34
1.74 0.69 0.24 0.08 8.09
0.98 0.34 0.11 0.03 4.87
0.75 0.24 0.07 0.02 3.28

Note. – The upper and lower values are not error but correspond to the HFE/LFE evolution scenarios of the faint-end slope. All quantities are in .

Table 3Extragalactic Background Light

The measured two-dimensional power spectrum from extragalactic sources consists of two components: i) the shot noise from the fluctuation in the number of unresolved sources entering the instrument beam, and ii) the clustering component arising from the correlation of galaxies on all scales. Additional power arising from local components such as Galactic cirrus and Zodiacal Light has been shown to be comfortably below the measured signal at 1-5 (Kashlinsky et al., 2005; Matsumoto et al., 2011; Kashlinsky et al., 2012). In comparing with observational data we adopt the convention for the power spectrum to approximate the root-mean-square fluctuations as (Kashlinsky, 2005). The angular power spectrum of galaxies projected onto the sky can be related to their evolving 3D power spectrum, , by the Limber approximation (for 1 radian) which we adopt as modified by Fernandez et al. (2010),

(11)

where is the comoving angular diameter distance. The quantity in the square brackets is the flux production rate which is empirically determined by our populated lightcones:

(12)
Figure 7.— Flux production rate (times z) as a function of redshift in the unresolved regime shown for limiting magnitudes of 22, 24, 26, 28 (solid, dashed, dot-dashed, dotted curves respectively). The total unresolved flux under each curve listed in Table 3. The figure illustrates how removal of ever fainter sources isolates the unresolved component to higher redshifts.

It is important to note that the process developed in Kashlinsky et al. (2005, 2007b) removes sources down to a fixed level of the shot-noise power (see Table 4). This is equivalent to removing galaxies down to a limiting magnitude, , so that the remaining unresolved background is given by Equation (12) integrated from to . In Figure 7 we show the unresolved background from our modeling as a function of redshift, which illustrates the process of galaxy removal down to fainter magnitudes isolating the background to progressively higher redshifts. Note, that there is very little contribution (0.1 ) from galaxies at z1 after removing galaxies down to 26 AB mag. We find that for a limiting magnitude brighter than 24 mag, the unresolved flux is mostly dominated by galaxies at intermediate redshifts whereas galaxies at the faint-end takes over once 24. In Table 3 we list the total integrated background in the 0.45-4.5 range including the unresolved background for different limiting magnitudes corresponding to the curves in Figure 7.

5.1. Shot Noise

Figure 8.— Shot-noise power amplitude after integrating the counts as a function of limiting magnitude (connected squares). The gray shaded area corresponds to the allowed range of the faint-end evolution of the LF. The thick gray lines show the levels of reached by Matsumoto et al. (2011) at 2.4, Kashlinsky et al. (2005) (dark) and Kashlinsky et al. (2007b) (light) at 3.6 and 4.5. The intersection corresponds to the limiting magnitude reached in these studies. We tabulate these values in Table 4. We point out that our model counts are immune to the effects of confusion.

The shot-noise level seen in fluctuation measurements is critically important in order to identify the nature of the unresolved populations Kashlinsky et al. (2007c). It can be described as statistical counting noise in the number of unresolved sources within the instrument beam and its power is,

(13)

Shot noise is a directly measurable quantity and is not affected by confusion which may be present. This allows us to evaluate the effective limiting magnitude, , for a given shot noise level using our models which are also immune to confusion. We calculate the shot noise associated with galaxies in our lightcones and display it in Figure 8 as a function of limiting magnitude at the relevant bands. As fainter galaxies are removed the shot noise drops steadily in the same manner as seen in measurements. At 22 mag we have already removed most galaxies at z1 beyond which the shot noise is mostly determined by the faint-end of the LF. The horizontal lines in Figure 8 show the levels reached by the studies listed in Table 4. The intersection with our models agrees well with Kashlinsky et al. (2007b) claiming to have removed galaxies down to 25-26 AB mag but is slightly brighter (24) for the levels reached by Kashlinsky et al. (2005) who claimed to reach 25 mag. Similarly, our shot noise levels agree well with those found by Matsumoto et al. (2011) after removing galaxies down to AB magnitudes 22.9, 23.2 and 23.8 in the AKARI/IRC bands at 2.4, 3.2 and 4.1 respectively. Table 4 lists the limiting magnitude predicted for the shot noise levels reached in several studies.

Reference
       Band [10 nWmsr] (AB)
Thompson et al. (2007a)
       F160W 1.0 27
Thompson et al. (2007b)
       F110W 1.8 27
Kashlinsky et al. (2005)
       IRAC1 5.8 24.4
       IRAC2 6.0 24.0
Kashlinsky et al. (2007b)
       IRAC1 2.0 25.1
       IRAC2 1.0 25.1
Matsumoto et al. (2011)
       IRC 82 23.2
       IRC 33 23.3
       IRC 8.1 23.9

Note. – The upper and lower values are not error but correspond to the HFE/LFE evolution scenarios of the faint-end slope. The values are inferred from Figure 3 of Matsumoto et al. (2011).

Table 4Limiting Magnitudes Implied by Shot Noise Levels

We have defined to separate resolved/removed galaxies from unresolved remaining sources. In practice, however, the accurate value of reached depends on the source detection algorithm and the photometric aperture used to derive magnitudes. Furthermore, source extraction can become limited by confusion, depending on exposure and instrument beam. Since our underlying reconstruction of galaxy counts from LFs is immune to confusion, we assume that the measured shot noise levels serve as a reliable indicator for the faintest sources removed, . This obviously assumes that the source removal is done properly and does not introduce spurious signals in the background fluctuations as discussed at length in Arendt et al. (2010). It also assumes that the (quasi-)flat power seen on small scales is entirely due to shot noise dominating the contribution from non-linear clustering of galaxies which we discuss in the following subsection.

5.2. Galaxy Clustering

The shape and amplitude of the fluctuations produced in each redshift slice is dictated by the two evolving quantities in the Limber equation (eqn. 11), i) the amount of light production given by our reconstructed in a given band, and ii) the clustering pattern of the sources in this epoch, described by their three-dimensional power spectrum, . For the latter quantity we assume that on large scales sources cluster according to the observationally established concordance CDM power spectrum. Prescriptions exist for non-linear evolution that modify the linear power spectrum in the regime where structures have collapsed out of the density field and linear theory breaks down (Peacock & Dodds, 1996; Smith et al., 2003). However, luminous sources are known to be biased tracers of the dark matter distribution particularly in the non-linear regime where the correlations of sources depends on the Halo Occupation Distribution (HOD) of galaxies. We therefore consider a halo model description of the power spectrum which decomposes it into two terms, a two-halo term () on large scales arising from the correlations of isolated halos, and a one-halo () from correlations of particles within the same halo on small scales (e.g. Ma & Fry (2000)). We follow the treatment of Cooray & Sheth (2002) and write,

(14)

where,

(15)
(16)

and is the halo mass function (Sheth et al., 2001, from), is the average number density of galaxies, is the linear CDM power spectrum (computed using the transfer function of Bardeen et al. (1986)), is the normalized Fourier transform of the halo profile (Navarro et al., 1996), and is the halo bias (from Sheth et al., 2001). The occupation number has been separated into central galaxies, , and satellite galaxies, , such that

(17)

We take the mass dependence of our HOD model to follow the four parameter description of Zheng et al. (2005):

(18)
(19)

where is characterized by , the minimum halo mass that can host a central galaxy and , which controls the width of the transition of the step from zero to one central galaxy. The satellite term has a cut-off mass which is twice as large as the one for central galaxies and grows as a power-law with a slope of , normalized by . This form has been explored both numerically and observationally. Since the measurements of HOD-parameters are obtained from samples of resolved galaxies at low-, their validity may not extend into the unresolved regime or, in particular, to higher redshifts. Since we are concerned with the unresolved regime it is important to note that the measured cut-off mass of central galaxies, , is typically set by the lowest luminosity probed by the survey so halos may continue to host central galaxies to lower masses but are excluded due to selection criteria. In Section 4 we showed how the unresolved light is typically dominated by the faint-end of the LF for 25 with most bright central galaxies removed out to 3 in measurements of CIB fluctuations. One would also expect the masking to eliminate most of the surrounding satellite galaxies. We have adopted the following parameters of the HOD-model motivated by SDSS measurements of Zehavi et al. (2011): , , , and where we have deliberately chosen a lower cut-off reflecting low mass halos hosting galaxies well into the unresolved regime, and a lower allowing for large amounts of unresolved satellite galaxies, while keeping . It should be noted, that in the absence of any HOD-assumptions, a simple linear CDM clustering with typical bias, , produces nearly identical fluctuations on large scales. The one-halo term has white-noise power spectrum (const) with its amplitude limited from above by the measurements at small scales and so its modeling is irrelevant to interpreting the clustering signal at scales .

We assume that unresolved sources in our lightcones are uniformly mapped onto the halo distribution i.e. the clustering is independent of galaxy luminosity. In practice however, we expect the most luminous galaxies to be removed in the masking process along with most of the accompanying satellites. This could motivate one to introduce an upper mass limit in the integrals in Equations (15), . However, this would require an additional mass-to-light ratio assumption and since it would always result in a decrease of the clustering amplitude, we do not apply and consider the result to be an upper limit for the resultant power spectrum. This includes the mass-dependent bias which is similarly integrated over the entire range of occupied halos (10). The large scale (linear regime) galaxy bias seen by Zehavi et al. (2011) in the local SDSS sample is 1 when all galaxies are included. At somewhat higher redshifts, Granett et al. (2012) find averaged over 0.5z1.2. Further increase of the linear bias with redshift is expected on theoretical grounds as collapsing density peaks were increasingly rare in the past. The bias prescription used here shows the same general behavior (Sheth et al., 2001). Several CIB studies at far-IR wavelengths claim a linear bias as high as =2-3 for far-IR sources (e.g., Lagache et al., 2007; Viero et al., 2009) but at these redshifts, the samples are already biased towards the most luminous objects due to selection effects. If anything, we expect the bias to be lower in the faint and unresolved regime after the more strongly biased luminous galaxies are masked and removed.

The large scale fluctuations are always dominated by clustering in the linear regime (two-halo term). On the other hand, the non-linear clustering described by the one-halo term in Equation (15) exhibits a =const behavior () making it indistinguishable from shot noise in measurements. Given that we found excellent agreement between the shot noise in our models and the measurements at the same magnitude levels, there does not seem to be any need to invoke non-linear clustering to explain fluctuations on small-scales (unless the data points deviate from a simple white noise spectrum ). In addition, we explored the pure dark-matter treatment of the non-linear clustering of Smith et al. (2003) but find it to be inconsequential in comparison with the shot noise dominated fluctuations on small scales. Although we see the one-halo term contributing somewhat to the HFE fluctuations in Figure 9, it becomes less relevant if one accounts for the more massive halos being masked/removed. In fact, we will see in Section 5.4 that current fluctuation measurements place a limit on the amount of non-linear power in the unresolved regime.

5.3. Comparison with fluctuations from the Millennium Simulation and Semi-Analytic Models

To compare our results with the clustering of halos seen in large scale N-body simulations, we have made use of the theoretical lightcones constructed by Henriques et al. (2012). These mock catalogs are based on semi-analytical models for galaxy evolution (Guo et al., 2011) which are implemented on two very large dark matter simulations, the Millennium Simulation (Springel et al., 2005) and the Millennium-II Simulation (Boylan-Kolchin et al., 2009). The simulations provide a description of the evolving spatial distribution of dark matter halos and subhalos whereas the nature of the baryonic content is described by the latest version of the semi-analytical Munich model (Guo et al., 2011). The Millennium Simulation follows structure formation in a box of side 500Mpc comoving with a resolution limit of M whereas the Millennium-II Simulation focuses on a region of 100Mpc but with complete merger trees down to M2. Henriques et al. (2012) use the Millennium Simulation only in their study, limiting the faint-end of the LF to halos 10M. Even so, the predicted faint near-IR counts are higher than observations suggest due to an unusually high abundance of relatively low mass galaxies (10M) at z1 (Guo et al. (2011) tuned their model to match the local populations). A comparison of the predicted correlation function of these models with local SDSS data shows decent agreement for massive galaxies whereas correlations of low mass systems are overpredicted, particularly at small separations. Henriques et al. (2012) also neglect the effects of dust and PAH emission and consider only starlight using stellar synthesis models of Bruzual & Charlot (2003) and Maraston (2005). For a detailed description of these models we refer to Springel et al. (2005), Guo et al. (2011) and Henriques et al. (2012).

Despite the limitations mentioned above, we find that this study provides a useful comparison to our fluctuation analysis. After constructing images using the publicly available mock data of Henriques et al. (2012), we calculate the projected angular power spectrum, convolved with the instrument beam. We analyze two independent regions observed in H, K, IRAC1 and IRAC2 each covering 1.41.4 degrees on the sky. We extract all galaxies in the magnitude range 30 to produce the unresolved fluctuations which we display alongside our results in Figure 9. Because of the overabundance of faint galaxies at 3.6 and 4.5 in the semi-analytical description of Guo et al. (2011), we need to remove galaxies down to 0.2 mag deeper than the listed in the panels in order to normalize to a common shot noise level. This NIR overabundance (despite the resolution limit of 10M) results in the Millennium fluctuations (dark-gray shades in Figure 9) being in closer agreement with our HFE scenario at 3.6 and 4.5 but can otherwise be considered to be consistent with our main results.

5.4. Results

Figure 9.— Models of the unresolved near-IR fluctuations compared to measurements from authors listed in the panels. We have chosen the limiting magnitude such that the models are normalized to the shot noise levels reached in these studies (including a contribution from a one-halo term). The solid curves show the total contribution from clustering and shot noise whereas the light shaded areas indicate the region bracketed by our HFE and LFE models. These are all suppressed by the instrument beam on small scales. The dotted lines indicate the separate one-halo and two-halo terms of the power spectrum. Shown in each panel is the total unresolved flux associated with the default model (), the values of (in units of nWmsr) and the associated . The dark shaded regions correspond to fluctuations arising from galaxies in the lightcones of Henriques et al. (2012) derived from the Millennium Simulation in the magnitude range 30. Because of their overabundance of faint galaxies at 3.6 and 4.5  we have increased the of the Millennium fluctuations by 0.2 mag to normalize to the correct shot noise levels. In the 3.6 panel we also show the default model from Sullivan et al. (2007) (dashed line). In the 1.6 panel the notation follows Fig. 2 of Thompson et al. (2007a): asterisks correspond to fluctuations with all sources removed whereas the triangles indicate their estimate of the instrumental Gaussian noise.
Figure 10.— The lines show the ratio of the measured source-subtracted power spectrum from the AKARI 2.4 data and the latest Spitzer-based measurements at 3.6 and 4.5 (Kashlinsky et al., 2012) to the HFE and LFE expectations (red/upper and blue/lower respectively). The results show that the measured CIB fluctuations continue to diverge from our models as we go to larger scales and are thus unlikely to result from extra-biasing of these faint populations: in order to explain the measured signal the biasing would have to be 1) scale-dependent, i.e. non-linear, 2) biasing amplification would have to be more non-linear on scales where the amplitude of the underlying correlation function is weaker (larger scales), and 3) the biasing would have to be different at 3.6 and 4.5.

With the emission history reconstructed from LFs and the sources distributed according to the halo model in Section 5.2, we projected onto the sky the clustering pattern in our NIR lightcones using Equation (11) and display the results in Figure 9. The limiting magnitudes have been chosen such as to normalize the shot noise (dot-dashed lines) to the measurements shown in each band. The shot-noise is seen to dominate the fluctuations on small scales whereas the clustering component becomes significant at arcminute scales. In the display we have chosen to focus on 1.6, 2.4, 3.6 and 4.5 where we can compare with measurements from Hubble/NICMOS, AKARI/IRC, and Spitzer/IRAC. Our models have been convolved with the beam profile (or PSF) of these instruments. It is immediately clear from Figure 9 that the contribution from known galaxy populations falls short of the measured clustering signal in every band shown. We briefly discuss each comparison:

Kashlinsky et al. (2007b) find excess fluctuations of 0.05-0.1 at arcminute scales in the Spitzer/IRAC channels after removing sources down to 25 mag or shot-noise levels nW/m/sr. It can be seen from Figure 9 that the known sources remaining at the measured shot-noise levels cannot account for the observed fluctuations for any faint-end modeling of the LF. We have displayed the data of Kashlinsky et al. (2005) and Kashlinsky et al. (2007b) in panels side-by-side illustrating that the discrepancy gets larger as galaxies are removed to deeper levels. The unresolved flux associated with our default model is 0.18  in the deepest 3.6 maps of Kashlinsky et al. (2007b), so in order to explain the observed level of the excess fluctuations the relative levels of the source-subtracted CIB fluctuations would have to be close to non-linear, 1, all the way to . The spatial spectra of the CIB fluctuations from the known galaxy populations is such that the gap increases toward large scales if this behavior of the source-subtracted CIB fluctuations continues as observed (say, ), so these fluctuations would have to be in the same (quasi)non-linear regime at much larger scales making it more difficult to explain them with the known galaxies. The additional linear biasing to amplify the arcminute scale signal to the observed levels but this would require which is highly unlikely for small systems in the range where most of the flux is produced. This, however, can be shown not to be viable in light of the latest Spizter-based results submitted after this paper. Kashlinsky et al. (2012) measure the source-subtracted CIB fluctuations on sub-degree angular scales confirming the earlier results and identifying, for the first time, the fluctuations spectrum to where the discrepancy continues to grow. Figure 10 shows the ratio of the measured power spectrum from the new large scale Spitzer/IRAC data of Kashlinsky et al. (2012) to the power spectra of our HFE and LFE (red and blue), illustrating that the data keeps diverging from our models out to . This shows that if one were to model the measured CIB fluctuations with extra biasing of the known galaxy populations, the biasing would have to be 1) highly scale-dependent, i.e. more prominent on larger scales, where density pattern is in linear, 2) the resultant biasing factors would have to be huge reaching amplifications of over two order of magnitude at the largest scales, and 3) the biasing would have to be wavelength dependent attesting to the different discrepancy ratios in each band. This argues further against the detected CIB fluctuations arising from the a faint-end extension of the known populations.

Following the HOD-model in Section 5.2 we find that small scale power from the one-halo term () is lower than the shot noise (). For the HFE model however, the two are comparable at 3.6 and 4.5 which leads to a slight increase in the upper values of in Table 4 (by 0.3 mag) as normalized by the measured shot noise levels. The fact that our shot noise, calculated from number counts, is in good agreement with the small scale data points of Kashlinsky et al. (2005, 2007b), argues against a significant clustering on small scales. It must therefore remain at or below the measured level of , otherwise sources would have to be removed to , which is not possible in the deepest Spitzer/IRAC maps.

Figure 11.— The contribution of different redshift bins to the unresolved IR-fluctuations shown in Figure 9 for the indicated in the panels. The 3.6 and 4.5 panels correspond to the models at the shot noise levels of Kashlinsky et al. (2007b). The different set of lines correspond to the redshift bins indicated in the legend. This illustrates that depending on the observed band and the depth of source removal, the unresolved fluctuations from known galaxies are dominated by populations at different epochs. The amplitude and shape is governed by 1) the flux production history (see Fig. 7), and 2) the evolving power spectrum, . The non-linear clustering component is important at low- but moves towards small scales for higher . The dependence on the comoving angular diameter distance, (see Eqn. (11)) is easily seen as the peak of the CDM power spectrum shifts towards smaller scales with increasing redshift.

Chary et al. (2008) stack deep Spitzer exposures to detect faint ACS galaxies beyond the detection threshold of the frames used in Kashlinsky et al. (2007b) and explore the sensitivity of the IR-fluctuations to these ACS sources. Their stacked source detections down to 26.0-26.2 mag imply a net flux of 0.12-0.35 . For comparison, the flux associated with our lightcones in the 25-26.2 mag range is 0.04 and 0.2  at 3.6 for LFE and HFE respectively with 0.04-0.35  from still fainter galaxies, 26.2 mag. We note that Kashlinsky et al. (2007a) demonstrate observationally the negligible correlations on arcminute scales between the source-subtracted CIB maps, as constructed by their self-calibration procedure Arendt et al. (2010) and ACS source maps.

Thompson et al. (2007a) measure fluctuations at 1.6 on scales out to 80 using HST/NICMOS (and at 1.1 in Thompson et al. (2007b)) and ascribe the signal to faint galaxies emitting at redshifts 0.5-1.5. Their fluctuations at 80 have amplitudes of 0.4 , which is a factor of 2-7 times higher than the total unresolved component, 0.06-0.20 , for sources fainter than 28 mag, indicating that the clustering of the underlying galaxies must be highly non-linear. For their CIB fluctuation levels to be reconciled with our empirical estimates, the one-halo term would have to be significantly higher, but then its amplitude would overshoot the data at all the other NIR wavelengths. If we take the upper limit on the shot noise at these wavelengths to be at the levels of the estimated instrument noise of Thompson et al. (2007a), then our shot noise already matches at AB magnitude of 27 (see triangles in Figure 9). But even at that level we cannot reproduce the fluctuations (asterisks) with the clustering of known galaxy populations out to . 3 We point out in this context the clearly visible outer halos of the sources removed by Thompson et al (2007a, see their Fig. 4) whose contribution to their CIB fluctuations shown may be significant and should be estimated for more quantitative conclusions at 1.6 .

The Matsumoto et al. (2011) measured fluctuations at 2.4, 3.2, 4.1 using data AKARI satellite and conclude that they are consistent with stars from early epochs confirming the identification proposed in Kashlinsky et al. (2005). The left panel in Fig. 10 confirms that the AKARI signal at 2.4 cannot be explained by the remaining known galaxy populations.

In Figure 9 we also display the default model from Sullivan et al. (2007) (dashed lines) who combined a halo model and conditional luminosity functions to calculate IR-fluctuations at 3.6. Our models have a somewhat lower amplitude considering the fact that we use =24.4 as opposed to the 25.3 mag used by Sullivan et al. (2007) (and quoted in Kashlinsky et al. (2005)) but the two are in rough agreement. For =25.3 our unresolved flux is 0.1  (LFE) which is roughly consistent with the 0.08  found by Sullivan et al. (2007). However, they claim that the fluctuations measured by Kashlinsky et al. (2005) at 3.6 can be explained by galaxies in the magnitude range 25.3 to 28.8 (AB) at 1-3. This is a somewhat puzzling conclusion when comparing their model with the data in Figure 9 as it clearly fails to account for the clustering excess4.

In Figure 11 we show the contribution of different redshift bins to the unresolved IR-fluctuations for the indicated in the panels of Figure 9. This illusrates the different epochs in which unresolved galaxy populations contribute to the fluctuations in different observed NIR bands. The redshift dependence is governed by 1) the flux production history (see Fig. 7), and 2) the evolving power spectrum, . The Figure also reflects the dependence on the comoving angular diameter distance, (see Eqn. (11)) with the overall clustering pattern shifting towards smaller scales with increasing redshift.

6. Summary and Discussion

We have reconstructed the emission histories seen in the near-IR of present-day observers to model the unresolved CIB fluctuations and compared with current measurements. Our compilation of  luminosity functions used to populate lightcones at 7 reproduces the observed number counts remarkably well and accounts for the features shaping them. We assume the Schechter-type LF and model the evolution of its parameters from the available datasets. We then considered high and low faint-end LF limits within the constraints permitted by deep galaxy counts data. Extending these to faint magnitudes and to high- we calculated the range of unresolved background flux in deep images and derived CIB-fluctuations from these galaxy populations predicted by the standard CDM clustering power spectrum. We find good agreement between the predictions of our analysis and semi-analytical galaxy evolution models combined with the large scale Millennium N-body simulation.

By varying the limiting magnitude of source subtraction we normalize our models to the observed shot noise levels, finding good agreement with the depths reached in current fluctuation measurements. We show that the known galaxy populations fail to account for the observed source subtracted CIB clustering signal in either LFE or HFE limits. Although, in principle, by varying one can find a population of brighter galaxies that matches the measured clustering amplitude at some fiducial angular scale, the associated shot noise levels always imply that all such populations have been removed in the source subtraction thereby not contributing to the unresolved fluctuations. Thus it means that the emitters producing the source-subtracted CIB fluctuations on arcminute scales are below the detection limits of current surveys and furthermore, cannot be a part of the known evolving galaxy populations. In other words, the only way to reproduce the clustering excess with extragalactic sources is by introducing a new population of sources that are significantly fainter than the detection threshold of current instruments i.e., a highly clustered population with low shot noise.

The high isotropy of the CIB fluctuation signal measured in Spitzer IRAC data Kashlinsky et al. (2012) argues strongly against the signal originating in Galactic or Solar system foreground emissions as well as very local extragalactic sources. Since the observed galaxy populations (extrapolated to very faint limits) cannot explain the measurements, the CIB fluctuations must originate in new populations so far unobserved in galaxy surveys. Kashlinsky et al. (2007a) show that there are no correlations between the ACS maps with sources down to AB mag of 28, and the source-subtracted CIB maps from Kashlinsky et al. (2007b). This implies that either the CIB fluctuations originate in a large unknown population of very small systems at low/intermediate redshifts, or they are produced by high redshift, z7, populations whose Lyman break (at rest 0.12) is shifted passed the longest ACS channel (at 0.9).

Acknowledgments

This work was supported by NASA Headquarters under the NASA Earth and Space Sciences Fellowship Program - Grant NNX11AO05H. KH is also grateful to The Leifur Eiriksson Foundation for its support. MR acknowledges partial support by NASA Grant NNX10AH10G and NSF CMMI1125285. We also thank to B. Henriques, R. Keenan and T. Matsumoto for useful exchanges and data. Our compiled database of Schechter parameters is available upon request.

Appendix A A. LF Binning and Interpolations

Because of degeneracy in (,,), different sets of Schechter parameters can represent LFs of very similar shapes. The method used in Section 3 disentangles the Schechter parameters to fit their evolution individually. In addition to this, we used an alternative approach in which the shape of each measured LF is kept intact. We took each LF in its rest-frame and redshift the associated emission to the observed wavelength, . We examine the all LFs that meet the criterion where is the center of the NIR band and is roughly the FWHM of the filter. The inserts in Figure 10 show the redshift distribution of available LFs which can be observed through JHKL. In a given band, we place each LF in redshift bins and take the functional average of in common bins so that we have a unique LF at each redshift. We thus have template LFs, , in each of the observed NIR bands and the rest of the analysis is identical to that in Section 3 following from Equation (6) (we interpolate the evolution and project the populations onto the sky). The major shortcoming of this method is the redshift information. Averaging over several LF in a common redshift bins is immune to the effects of Schechter parametrization but comes at the cost of crude evolution i.e. the sampling of is determined by the number of -bins. As seen in Figure 12 there is no guarantee that there exists a LF measurement falling into in each redshift bin. In this case we borrow LFs from neighboring wavelengths scaling them according to synthetic spectra. Figure 10 shows that despite these limitations, we obtain very comparable number counts to the ones in Section 4, agreeing to within 20% in the relevant magnitude range.

Figure 12.— Comparison between our default method (dashed) and the alternative method presented here (solid). The two curves agree to within 20% in the range shown. The data shown in the background is the same as in Figure 5. The insets show the redshift distribution of LFs avalable in for the calculation in each band (i.e. ).

Appendix B B. Consistency notes

K-correction

Calculating the absolute magnitudes of a galaxy sample requires a k-correction to account for the offset in the rest-frame and the observed SED due to the cosmological redshift (e.g. Hogg et al. (2002); Blanton et al. (2003a))

(B1)

where refers the band of interest. The k-correction can be written (in AB magnitudes)

(B2)

where is the observed brightness of a galaxy at redshift and is its rest-frame brightness in -band. The exact value of the k-correction requires knowledge of the spectral energy distribution (SED) of the source and is commonly evaluated by assuming a template SEDs based on the galaxy type/color. This treatment is fairly reliable for low- galaxies but the correction can become large for high- galaxies and dominate the uncertainty in the derived LF, especially in the blue bands. Recent multiband photometric surveys offer a robust way of reducing this SED dependency by utilizing magnitudes in multiple bands to constrain the best-fit SED . Not only does multiband coverage indicate SED shape but when probing the LF in the rest-frame band centered at , the galaxy flux can be sampled in the band which is closest to . In other words, the observed filter () that best matches the redshifted rest-frame band of interest is the one that minimizes . The k-correction needed then becomes the matter of setting this quantity to exactly zero which is typically a small correction. We can rewrite Equation (B1) in this framework

(B3)

where the SED dependence of the k-correction is now small even at high redshifts. Backtracking the original procedure to apparent magnitudes now requires simply which we use in Equation 7.

Photometric Systems

Unfortunately, there is no photometric system which is universally accepted and the different ways used to evaluate the apparent magnitude of galaxies in the survey can introduce biases affecting the derived luminosity functions (see Bessell (2005) for a review of photometric systems). As the flux from a galaxy diminishes from the center it will eventually drop below the background noise to be missed by the aperture. Photometric systems based on total magnitudes, such as Sérsic, are usually preferred since they directly quantify the physical flux while apertures such as Kron and Petrosian will always suffer from missed light to some extent. However, total magnitudes typically assume an extrapolated profile which is model-dependent and has larger measurement errors (Cole et al., 2001). The Petrosian system can be advantageous since it compensates for the effects of seeing by increasing the fraction of the light recovered from a galaxy when its angular size is small (Blanton et al., 2001). Despite this, Petrosian magnitudes are found to underestimate Sérsic by 0.2 mag (Strauss et al., 2002; Blanton et al., 2001). Likewise, 2MASS Kron and isophotal magnitudes may account for only 50-80% of the total flux in the most extreme cases (Andreon, 2002)). For example, Smith et al. (2009) show that their UKIDSS Petrosian magnitudes can be up to 0.5 mag fainter than 2MASS Kron magnitudes. The fraction of the lost flux increases towards fainter galaxies and may cause a systematic underestimation of the faint-end luminosities as well as the luminosity density. Hill et al. (2011) provide a good analysis of the effects of different photometric systems used in surveys. They find an overdensity of faint galaxies when compared with the best-fit Schechter function irrespective of the aperture system used and show that a Schechter function parametrization does not provide a good fit at the faint-end. They also show that the use of a photometric systems based on total magnitudes (e.g. Sérsic extrapolated) have a systematically steeper faint-end slope than photometric systems based on Kron or Petrosian magnitudes. They further show that the r-band Kron & Petrosian photometry underestimates the luminosity density by at least 15% as they do not account for missing light. Blanton et al. (2003b) show that the difference of the luminosity density resulting from Petrosian and Sérsic magnitudes should be within 0.1 mag in the SDSS bands and not worth correcting for given the limitations of both systems. Still many authors apply a correction to estimate the total magnitudes in order to derive quantities such as the luminosity density in physical units (e.g. Kochanek et al. (2001); Bell et al. (2003); Eke et al. (2005)). These can be as high as 0.3 mag in the K-band. It seems that uncertainties in the LF may be dominated by the aperture governing the fraction of flux recovered, especially at the faint end.

Luminosity Function Estimators

In this paper we use LFs derived from a variety of different LF estimators. The choice of LF estimator is unlikely to be a major source of discrepancy between the LFs derived by different authors although it can lead to different combinations of the Schechter parameters. The most commonly used methods are i) the method (Schmidt, 1968), ii) the Sandage-Tammann-Yahil maximum likelihood method (STY) (Sandage et al., 1979) and iii) the StepWise maximum Likelihood Method (SWLM) (Efstathiou et al., 1988). The method is reliable in the sense that it simultaneously gives the shape and normalization of the LF requiring no assumption on the parametric form for the LF. However, it suffers from systematic biases in the presence of density inhomogeneities in the observed field. The STY method is typically preferred when estimating the LF over multiple fields since it has been shown to be unbiased to large scale structure and does not require binning of the data (Efstathiou et al., 1988). It does however require an assumption of a functional form of the luminosity function. The SWML method is widely used since it makes no assumption of the LF shape while still being insensitive to large scale structure. Willmer (1997) compare the properties of each LF estimator and show how different LF estimators tend to be biased towards the faint-end either overestimating or underestimating the slope, depending on the estimator and the underlying catalog. In order to minimize such effects one routinely compares the outputs of more than one method (e.g. Bouwens et al. (2007); Ilbert et al. (2005); Cirasuolo et al. (2010)).

Footnotes

  1. In the rest-frame UV/optical, where the low- contribution does not matter for the observed NIR, we fix the low- slope at -0.9 and -1.1 for LFE and HFE respectively.
  2. The Millennium Simulation and the resulting lightcones of Henriques et al. (2012) assume a WMAP1-based cosmology (Spergel et al., 2003) with parameters = 0.73,