Improved Models for Cosmic Infrared Background Anisotropies: New Constraints on the IR Galaxy Population
Abstract
The power spectrum of cosmic infrared background (CIB) anisotropies is sensitive to the connection between star formation and dark matter halos over the entire cosmic star formation history. Here we develop a model that associates starforming galaxies with dark matter halos and their subhalos. The model is based on a parameterized relation between the dustprocessed infrared luminosity and (sub)halo mass. By adjusting 3 free parameters, we attempt to simultaneously fit the 4 frequency bands of the Planck measurement of the CIB anisotropy power spectrum. To fit the data, we find that the starformation efficiency must peak on a halo mass scale of and the infrared luminosity per unit mass must increase rapidly with redshift. By comparing our predictions with a wellcalibrated phenomenological model for shot noise, and with a direct observation of source counts, we show that the mean duty cycle of the underlying infrared sources must be near unity, indicating that the CIB is dominated by longlived quiescent star formation, rather than intermittent short “star bursts”. Despite the improved flexibility of our model, the best simultaneous fit to all four Planck channels remains relatively poor. We discuss possible further extensions to alleviate the remaining tension with the data. Our model presents a theoretical framework for a future joint analysis of both background anisotropy and source count measurements.
keywords:
submillimetre: diffuse background – submillimetre: galaxies – galaxies: star formation – galaxies: halos1 Introduction
About half of the optical/UV emission from stars is absorbed by dust
in galaxies and reemitted at infrared wavelengths
(Dwek et al., 1998; Fixsen et al., 1998). Measurements of the CIB therefore offer an
important channel to study the starformation history of the
Universe. In particular, CIB measurements are sensitive to the star
formation out to high redshifts, due to the correlation between the
infrared luminosity of a galaxy and its star formation rate
(Kennicutt, 1998). Owing to the relatively poor angular resolution
of farinfrared telescopes, it is challenging to resolve and study
individual infrared sources. Surveys to date are only able to resolve
the brightest sources, which are responsible for less than 15% of the
total CIB (Oliver et al., 2010b).
Although the Atacama Large Millimeter/Submillimeter Array
(ALMA)
Although fluctuations in the CIB were detected over a decade ago
(Lagache & Puget, 2000; Matsuhara et al., 2000), these were limited to probing only the
Poisson contribution to the power spectrum.
Robust measurements of the clustered component in the CIB anisotropy,
probing the spatial correlations of the underlying sources, have
become possible relatively recently, and are now rapidly improving.
The first detection of clustered power came from 160m Spitzer data (Grossan & Smoot, 2007; Lagache et al., 2007). This was followed recently
by the Balloonborne Large Aperture Submillimeter Telescope
(BLAST)
To make full use of this information, and to draw conclusions about the cosmic star formation history, requires the development of a modeling framework. With too complex a model, degeneracies between model parameters threaten to make meaningful conclusions impossible. On the other hand, if the model is too simple, it may be unable to accommodate all relevant information, such as observations of bright individual sources – or even worse: we may risk having conclusions that are artificially driven by inadequate modeling assumptions.
The Planck and Herschel measurements, the most informative to date, have so far been interpreted via halo models of largescale structure (Cooray & Sheth, 2002). These models take advantage of the relatively small uncertainty with which, given a model, statistical properties of the spatial distribution of dark matter can be calculated. They do so by dividing the modeling into two stages. First, the spatial distribution of dark matter is described (in particular, by assuming that all dark matter is confined to lie with spherically symmetric, collapsed halos). Second, the observable of interest is related to the distribution of dark matter. A key element for the latter is the halo occupation distribution (HOD) function of galaxies, which is a statistical description of how galaxies occupy halos, depending on halo mass and redshift.
To avoid degeneracies that might arise in the presence of too many free parameters, the halo model in Ade et al. (2011) and Amblard et al. (2011a), to which we refer as the “standard model” in the following, has been relatively simple – ignoring many details of how infraredluminous galaxies trace dark matter halos. For example, the HODs were assumed to be the same for all spectral types and to not vary with redshift. Even quite coarse details, such as the dependence of the luminosity of a galaxy on halo mass have been neglected.
These models have the virtue of simplicity, and without the data demanding more complicated models, it would perhaps not be worth abandoning such simplicity. The primary motivation of our study, however, is that the existing data are already difficult to understand with these simple models. For example, in the models used to interpret the Planck CIB measurements, there is no single choice of parameters that generates a good fit, simultaneously, to the power spectra at all frequencies (Ade et al., 2011). Also, because these models do not include a dependence of the IR luminosity on halo mass, neither the number counts of individual sources (detectable at the bright end), nor the shotnoise power levels can be calculated reliably.
In this paper, we extend the modeling framework used to date, by allowing galaxies in halos with different masses to have different luminosities. We do so by assuming that satellite galaxies occupy dark matter subhalos, and use a subhalo mass function derived from cosmological N–body simulations. The luminosity of each satellite galaxy is then parametrically related to the mass of the subhalo. Ours is not the first CIB model to allow for IR galaxies to have different luminosities. Amblard & Cooray (2007) discussed such a model, based on the conditional luminosity function (Cooray & Milosavljević, 2005; Yang, Mo & van den Bosch, 2003; Yang et al., 2005). What is new in the present work is the explicit association of IR galaxies with subhalos, and the use of such a model for interpreting the Planck data. As we will show below, the comparison of this model to the data indeed leads to novel conclusions.
The modeling extension necessarily comes at the cost of a few additional parameters, but has several benefits. First, and perhaps most importantly, it allows us to estimate a lower limit to the bright end of the number counts. Second, by introducing a duty cycle, and comparing to source count observations, we can put limits on the duty cycle. In the context of our model, the data favor duty cycles near unity. Physically, this implies that the bulk of the CIB is produced by quiescent star formation, rather than by the intense, short starbursts expected to be triggered in major mergers (Sanders et al., 1988; Barnes & Hernquist, 1991). Third, there exists a tension between the value of the powerlaw index , describing the dependence of the number of satellite galaxies on halo mass, inferred either from cosmological simulations, from optical galaxy data, or from the CIB anisotropy measurements (see §4). Our model naturally resolves this tension.
The remainder of this paper is organized as follows. In §2, we review recent measurements of the CIB anisotropy power spectrum by BLAST, SPT, and the space missions Herschel and Planck. In §3, we first review the basic formulae of the existing halo–based CIB models, and then describe our extensions that include a more realistic relation. In §4, we discuss the ingredients of the extended model in more detail, and illustrate how each ingredient affects the resulting CIB anisotropy power spectra. In §5, we apply our model to the recent measurement by Planck, and derive constrains on both the relation and on the shape of the mean spectral energy distribution (SED) of the sources. With enhanced flexibility, we attempt to simultaneously fit the data at all measured frequencies. In §6, we examine the shot noise levels and source counts computed with the new model, compare them with an observationallycalibrated model and a direct measurement, and place a constraint on the duty cycle of the underlying infrared sources. In §7, we further discuss the details of the fit, as well as several caveats and potential future improvements to the model. Finally, in §8, we summarize our main conclusions. Throughout this paper, we adopt the standard CDM model as our fiducial background cosmology, with parameter values consistent with the seventh year WMAP results (Komatsu et al., 2011), i.e., {,,,,,} = {0.274, 0.726, 0.045, 0.705, 0.810, 0.96}.
2 Measurements of CIB anisotropies
Substantial progress has been made in measuring the power spectrum of the CIB anisotropies in the past few years, thanks to the joint efforts from BLAST, SPT, Herschel and Planck.
In 2009, BLAST measured the power spectra in the GOODSSouth field at 250, 350, and 500 (i.e., 1200, 857, 600 GHz in frequencies, respectively) in the multipole range of (Viero et al., 2009). After subtracting the Galactic cirrus and Poisson noise, they found the variance of the CIB over the scales of 5’ – 25’ is consistent with a constant amplitude of relative to the CIB mean at all observing frequencies. The data could alternatively be well fit by a linear halo bias model, with bias parameters, and at 1200, 857, 600 GHz, respectively. Interestingly, the data could not be fit by the standard halo model based formula, with data points lying below the model curves at small scales and above the model curves at large scales. Later comparison with the Planck measurement shows a general consistency between the two measurements. However, the data points from the BLAST measurement are systematically higher at large scales. This seems to suggest that the BLAST measurement have been contaminated by residual Galactic cirrus at the large scales (Ade et al., 2011). If true, the failure of the halo model fit, likely due to the contamination of the data, does not necessarily mean the model is wrong.
Hall et al. (2010) reported SPT measurements of the auto and crosscorrelation at 150 GHz and 220 GHz in the multipole range of . They also determined that the spectral indices of the Poisson and clustered components between the two observing frequencies are and , respectively. This implies a steep slope of the SED of the dust emission (i.e. a graybody powerlaw index of ; see definition in eq. 24 below). Although the SPT measures the small–scale clustering, a model based on linear halo bias was still able to provide an acceptable fit to the data. This does not necessarily mean the clustering is linear at such small scales, but rather shows that the shape of non–linear power spectrum over the scales probed by SPT can not be distinguished from that of the linear power spectrum. A clear detection of the excess power from non–linear clustering requires the measurement to cover the transition at multiples of . This is exactly the angular coverage provided by Herschel and Planck.
Soon afterwards, Dunkley et al. (2010) published their measurement of the power spectra at 148 GHz and 218 GHz, based on the ACT data collected during their 2008 season over 296 . The clustered component of the IR power spectrum is detected at , assuming an analytic model for its shape. The spectral index of unresolved IR emission between the two frequencies is found to be , consistent with the SPT measurement.
With the Spectral and Photometric Imaging Receiver (SPIRE) on board, the Herschel Space Observatory, during its Science Demonstration Phase, measured the power spectra of CIB anisotropies in the Lockman Hole and GOODS South field at 250, 350, and 500 over the scales 1’ – 40’ ( in multipoles; Amblard et al. 2011b). The coverage of a wide range of angular scales, combined with clean data, yielded a clear detection of the non–linear clustering signal. The standard halo model–based prescription provides satisfactory fits to the data, but with different sets of bestfit parameters at different frequencies. The minimum halo mass for star–formation, , is constrained to be , which the authors interpret as the most efficient mass–scale for star formation. This mass–scale is lower than previous predictions by semianalytical models for galaxy formation (González et al., 2011).
Most recently, the Planck team published their measurement of the angular power spectra of CIB anisotropies, using maps of six regions of low Galactic dust emission (Ade et al., 2011, hereafter A11). The power spectra were determined in 4 frequency bands (217, 353, 545, and 857 GHz) over multipoles between and . Planck found CIB anisotropies at the level compared to the mean CIB, and also that the fluctuating component and the mean CIB have similar frequency spectra, consistent with the BLAST results. Unlike other missions, Planck was unable to independently measure the amplitude of the shot noise, due to its relatively poor angular resolution ( arcmins). Instead, it adopted values calculated using a parametric model by (Béthermin et al., 2011, hereafter B11). Again, different frequency bands yield different bestfit model parameters.
3 Halo model for CIB anisotropies and its extension
3.1 Review of the existing models
The angular power spectrum of CIB anisotropies is defined through,
(1) 
where denotes the observing frequency and is the measured specific intensity at that frequency. In a flat universe as assumed throughout this paper, the specific intensity is related to the comoving specific emission coefficient via
where is the comoving distance to redshift , and is the scale factor. Combining equations (1) and (3.1) and employing the Limber approximation (Limber, 1954), we obtain
(3) 
where is the 3D power spectrum of the emission coefficient, and is defined as follows,
(4) 
The existing models equate with the galaxy power spectrum , assuming the CIB is sourced by galaxies, and that the spatial variations in the emission coefficient trace the galaxy number density,
(5) 
On large scales, in the linear regime, the galaxy power spectrum follows the linear matter power spectrum, , while on small, non–linear scales, is typically computed from the halo model. In the framework of the halo model, the galaxy power spectrum is a sum of three terms
(6) 
where and account for contributions from galaxy pairs in the same halo and in different halos, respectively, and is the shot noise. The analytical expressions for and are (Cooray & Sheth, 2002),
(7)  
and
(8)  
where is the halo mass, is the halo mass function, is the Fourier transform of the halo density profile, and and , specified by the HOD, are the number of central and satellite galaxies inside a halo, with . Motivated by simulations (e.g., Kravtsov et al., 2004), is typically modeled as a simple step function,
(9) 
while is parameterized by a power law,
(10) 
Here, is the arbitrary “pivot” halo mass that hosts one satellite galaxy. Note that in the above equations, , , and could be thought of as the average number of galaxies inside halos with a fixed mass , with negligible scatter. The parameterizations above have been extended to include non–zero stochasticity in the relationship between halo mass and the number of galaxies (e.g., Zheng et al., 2005).
3.2 The improved model
The model above assumes that emissivity density traces galaxy number density (equation 5). This assumption implies that all galaxies contribute equally to the emissivity density, irrespective of the masses of their host halos. In other words, it assumes that all galaxies have the same luminosity. Realistically, however, both the luminosity and the clustering strength are closely related to the mass of the host halo. In particular, galaxies in massive halos are likely both luminous and highly clustered. In general, introducing a monotonic relation should give rise to stronger clustering than the existing models, in which the more strongly clustered sources (massive galaxies) have been assigned artificially low luminosities. This improvement is indeed the main new feature of our model.
In the following, we abandon the assumption of a mass–independent luminosity, and use the standard halo model to directly compute the power spectrum of the emission coefficient, . Note that the emission coefficient is related to the underlying galaxy population as follows,
(11) 
where denotes the infrared luminosity and is the infrared galaxy luminosity function. Neglecting any scatter, the galaxy luminosity is a function of the mass of the host dark matter halo (for central galaxies), or subhalo (for satellite galaxies), and equation (11) can be rewritten as:
where is the subhalo mass, and is the subhalo mass function. Studies show that the luminosity of a satellite galaxy best correlates with the mass or circular velocity of the host subhalo at the time it is accreted into the main halo, i.e. before it looses significant mass due to tidal stripping (Nagai & Kravtsov, 2005; Vale & Ostriker, 2006; Conroy, Wechsler & Kravtsov, 2006; Wang et al., 2006; Wetzel & White, 2010). We therefore use this “unstripped” mass to infer the infrared luminosity.
The computation of the power spectrum of the emission coefficient closely follows that of the galaxy power spectrum in the previous section. The only difference is that the galaxy numbers ( and ) should be replaced by expressions accounting for contributions from galaxies to the overall emissivity. Note that apart from the additional weighting by luminosity, equation (3.2) is essentially identical to the expression for the number density:
(13) 
To compute , we therefore only need to replace by
(14) 
and by
(15) 
in equations (7) and (8). The final results are
(17) 
where
It is straightforward to show that equations (3.2) and (17) reduce to equations (7) and (8) for a flat relation (galaxies of different masses have the same luminosity), as assumed in previous models.
The main ingredients of this model, such as the halo mass function, halo bias, and halo density profiles, have been carefully calibrated using numerical simulations. Throughout this study, we define halos as overdense regions with a mean density equal to 200 times the mean density of the universe. We assume an NFW profile (Navarro, Frenk & White, 1997) for the halo density, and adopt the fitting function of Tinker et al. (2008) for the halo mass function and its associated prescription for the halo bias (Tinker, Wechsler & Zheng, 2010). For the subhalo mass function, we use the fitting function of Tinker & Wetzel (2010, equation (12) in their paper).
We note that a luminosity–weighting scheme similar to the above have been explored by Sheth (2005) and Skibba et al. (2006), and applied to study resolved sources and the environment–dependence of galaxy properties. To the best of our knowledge, this is the first time that such a weighting scheme has been used for an analysis of unresolved brightness fluctuations.
Finally, in addition to the clustering caused by correlated large–scale structures, measurements of the CIB power spectrum include shot noise from random fluctuations in the discrete number of galaxies. In principle, this shot noise can be computed in the above model, through the equation
(19) 
where includes both halos and subhalos, and is the source flux,
(20) 
and is the flux above which individual sources are detected and removed in a given experiment. Unfortunately, this computation will likely remain inaccurate, given the simplicity of our model. In particular, our model neglects any scatter in the relation, for simplicity. This scatter is particularly relevant to the calculation of the shot noise: as is clear from equation (19), shot noise involves the integration of , which increases with scatter for a fixed mean relation. By neglecting scatter, our model therefore underestimates the shot noise, and our result should be understood as a lower limit. For this reason, we adopt the values computed using the parametric model of B11 (their model successfully fits bright source count measurements, and is also adopted by A11 as the shot noise) when we fit the measured angular power spectra. The values can be found in in Table (2) below (and also in Table 3 of A11). However, later in § 6, equation (19) and (20) are employed to constrain the duty cycle of the underlying sources.
4 Modeling the relation and its effect on CIB anisotropies
4.1 Parameterizing the relation
Before applying our model to the Planck data, we here specify and discuss the model ingredients in detail. In the halo model, the galaxy power spectrum is fully determined by the HOD, namely the functions and . In our model for the background fluctuations, the power spectrum depends, additionally, on the function . The latter depends on three variables: the redshift , the mass of the host (sub)halo, and the observing frequency . Limited by the current data quality, we adopt a relatively simple form for in this study. First, we do not make a distinction between halos and subhalos of the same mass, i.e., we assume that the and relations are identical. Second, we assume, for simplicity, that all galaxies have the same SED, independent of their masses and redshifts, and that the relation does not evolve with redshift except for an overall normalization. The dependence of on the three variables is then separable, and can be written as the product
(21) 
where is a normalization factor. In the following, we discuss each of the three components in more detail.
(1) Redshift evolution
For a given (sub)halo mass, the luminosity and the star formation rate (SFR) are expected to increase with redshift, because of the higher gas accretion rates, higher gas fractions, and more compact geometries at early times. Additionally, mergers, which can trigger starbursts, are more frequent at high redshift. In this study, we parameterize as a power law,
(22) 
This is partly motivated by the study of the specific star formation rate (sSFR; SFR per unit stellar mass). If we assume the stellar mass to halo mass ratio evolves only mildly with redshift (supported by semianalytical studies such as Neistein et al. 2011), the sSFR should have a redshift evolution similar to that of due to a correlation between SFR and infrared luminosity (Kennicutt, 1998).
The evidence for a smooth power–law evolution of the sSFR with redshift is somewhat mixed. Semianalytical models of galaxy formation indeed show that the redshift evolution of sSFR, assumed to be primarily driven by the fresh gas supply, could be well fit by a power–law with a slope of (Neistein & Dekel, 2008; Dekel, Sari & Ceverino, 2009; Oliver et al., 2010a). Observations, however, suggest a more complicated shape, with a steep evolution below , and a plateau beyond this redshift (see Bouché et al. 2010 and Weinmann, Neistein & Dekel 2011 for a compilation of measurements). Such a transition is hard to explain from a theoretical perspective (Bouché et al., 2010; Weinmann, Neistein & Dekel, 2011). Due to the different shapes, theoretical models tend to predict a sSFR lower by a factor of a few at and much higher at than inferred from observations, even though these models could more or less reproduce the redshift evolution of the overall SFR (Weinmann, Neistein & Dekel, 2011). On the other hand, there are large uncertainties in current observations, especially at high redshifts. For example, measurements by Yabe et al. (2009) and Schaerer & de Barros (2010) suggest the sSFR continue to rise beyond . Given these uncertainties, we consider two scenarios: has a single value at all redshifts (case 0), and is fixed to 0 at (case 1).

(2) LM relation
A robust conclusion from comparisons between observed galaxy luminosity functions and halo mass functions is that star formation is only effective in a certain range of halo masses. In particular, starformation must be suppressed at both low and high halo masses – for example, by feedback processes such as photoionization heating, supernovae, active galactic nuclei (AGN) and virial shocks (Birnboim & Dekel, 2003; Kereš et al., 2005; Dekel & Birnboim, 2006; Bower et al., 2006; Croton et al., 2006). The quantity of particular interest is the halo mass scale at which the star formation is most efficient, corresponding to a peak in the luminositytomass () ratio. Assuming, for simplicity, that as a function of mass is log–normal, we model the relation as follows,
(23) 
Here describes the peak of the specific IR emissivity per unit mass, and describes the range of halo masses for producing IR luminosity. At the low–mass end, in addition to the exponential drop in equation (23), we impose a minimum mass of , i.e., for . This cutoff is motivated by the work of Bouché et al. (2010). Later, we will see the bestfit is more than 2 away from the cutoff, so effects of applying such a cutoff are mild.
(3) SED shape
We adopt the same form for the SED as Hall et al. (2010). The SED is a graybody at low frequencies, and is a power law at high frequencies,
(24) 
where is the Planck function and is an effective dust temperature. These two functions are connected smoothly at the frequency that satisfies
(25) 
We fix as in Hall et al. (2010), but, unless stated otherwise, we allow and to be free parameters.
4.2 Effects on CIB anisotropies
Next, we illustrate how the CIB anisotropy power spectrum depends on each of the model parameters. For this purpose, we adopt the following set of fiducial model parameter values: , , , K, . The choice of is based on low redshift measurements of the sSFR (Noeske et al., 2007; Dunne et al., 2009; Oliver et al., 2010a; Rodighiero et al., 2010); while is around the peak of the stellarmass to halo mass ratio found by the abundance matching technique (e.g., Guo et al., 2010). The value K is consistent with measurements such as Dunne et al. (2000, K), Chapman et al. (2005, K), and Amblard et al. (2010, K). The assumption follows the common practice in the literature, and is supported by many theoretical considerations (e.g., Draine & Lee, 1984; Mathis & Whiffen, 1989). The choice implies that the “efficient” mass range for star–formation covers about one and half orders of magnitude in halo mass. Though this number is somewhat arbitrary, it is unlikely to severely bias the final results, since, as demonstrated below, the power spectra are relatively insensitive to .
We next vary one parameter at a time, and examine how the power spectrum changes compared to the fiducial case. The overall normalization is fixed using the amplitude of the mean CIB measured by FIRAS (this assumption will be discussed and relaxed in §7 below). More specifically, we choose such that the ratio of the predicted and measured mean intensity, averaged over all frequency bands, is equal to unity,
(26) 
Here is the number of frequency bands in which the power spectra have been measured (for the most recent Planck measurement by A11, ), is the intensity at frequency predicted by our model while is the same quantity measured by the FIRAS instrument. For the latter, we use the values quoted in Gispert, Lagache & Puget (2000).
The results are shown in Figure 1 for the highest (857 GHz) and lowest (217 GHz) frequency bands of Planck. The power spectra in the fiducial model are plotted as solid thick curves, against which the other curves should be compared. There are two distinct ways in which the parameters can affect the power spectra: either by changing the shape of the spectrum () or by changing the effective bias factor (). If the parameter enters primarily through the spectrum, then it must have opposite effects at high and low frequencies, since the average intensity is kept fixed by the normalization (equation 26). On the other hand, if the parameter enters primarily through the bias factor, then it has similar effects at different frequencies. This distinction helps us recognize how the parameters affect the power spectra, and also indicate their degeneracies – parameters which enter the same way are more likely to be degenerate. Figure 1 reveals that three of the parameters, , and , affect the spectrum, while the other two parameters, and , enter through the bias factor. Consequently, degeneracies are likely strong among , and , and between and .
Interestingly, increasing not only changes the overall clustering strength, but also increases the amplitude of the nonlinear (1halo term) contribution, relative to the linear (2halo term) contribution to the power spectrum. This is seen clearly as the upturn in both panels of Figure 1 in the short–dashed [blue] curve at , and is similar to the effect of increasing the value of the parameter in the previous models (see equation 10). In this sense, differs from , since the latter mostly affects the overall clustering. This feature of our model arises from the nonflatness of the relation, and is worth discussing further, because, as we will see below, the data requires significant power on small scales, which can be provided by the 1halo term.
For the intensity fluctuations, what matters for the relative strength of the 1halo and 2halo terms is the total luminosity, rather than total number, of satellite galaxies. To be more concrete, let us suppose that the total luminosity of satellite galaxies inside a halo of mass can be parameterized by a power–law, similar to equation (10),
(27) 
It is then , rather than , that determines the small scale clustering. If the luminosities of galaxies were independent of the masses of their host (sub)halos, would be equal to . However, for a more realistic relation, for which luminosity and mass are positively correlated over the relevant mass range, is larger than . This is because in addition to the number of subhalos increasing with the halo mass, the average subhalo mass increases, as well. Consequently, satellite galaxies are, on average, brighter in more massive halos, rendering larger than .
To illustrate this point explicitly, in Figure 2 we show the total number (solid curve) and the total luminosity (dashed and dotted curves) of satellite galaxies inside a halo of mass . All curves are shown normalized by their values at . As expected, the total luminosity increases more rapidly with halo mass than the number of satellites. Further comparing the dashed () and the dotted () curves shows that the slope of the relation itself increases with . This explains why the amplitude of the 1halo term increases with compared to the 2halo term.
In their recent analysis of the CIB power spectrum, Amblard et al. (2011a) have found that in order to fit the (large) observed power on small angular scales, the number of satellite galaxies had to be increased compared to that expected from numerical simulations and optical surveys. In particular, the expectation is and (e.g., Gao et al., 2004; Kravtsov et al., 2004; Zheng et al., 2005; Hansen et al., 2009), while the analysis of Amblard et al. (2011a) requires either or, effectively, (reducing and increasing have similar effects on the power spectrum since both increase the number of satellites and therefore raise the small scale power). Our model naturally resolves this tension by distinguishing from , and, as we will demonstrate below, is able to fit the data without the presence of any additional low–mass satellite galaxies.

5 Constraints from Planck
reduced  
Case 0  1.82 (33)  
Case 1  2.78 (33)  
Case 2  1.12 (24)  
reduced  
Case 3  2.43 (33)  
Case 4  1.70 (32) 

We are now ready to derive constraints from the Planck data. Given the flexibility of our extended model, we attempt to fit the CIB power spectra in the four Planck frequency bands simultaneously. As mentioned above, however, there are strong degeneracies among our 5 model parameters, which prevent us from simultaneously constraining each. We therefore first proceed by fixing two of the parameters – the dust temperature , and the massrange for significant luminosity, – at their fiducial values. In addition, we also apply a flat prior on , . Hereafter, this will be referred to as “case 0”; we will explore other possibilities below.
To obtain the bestfit values and their confidence levels, we adapted the CosmoMC Monte Carlo Markov Chain (MCMC) code (Lewis & Bridle, 2002), and applied it to the data in the four Planck channels presented in A11. The chains are checked visually and diagnosed using the statistical tests provided by CosmoMC to ensure their convergence. The bestfit values and the marginalized errors are listed in the first row of Table 1 and the power spectrum in the bestfit model is shown explicitly at each of the four frequencies by the solid [red] curves in Figure 3. As this table and the figure demonstrates, with our admittedly crude assumptions, and varying only three parameters, we are able to obtain good fits to the Planck data, with the exception of the high data points at 545 GHz. where we underpredict the observed power.
The bestfit is qualitatively consistent with the rapid evolution of the sSFR indicated by recent measurements. As we will see below, the uncertainties induced by our model assumptions are likely much larger than the statistical errors. It is nevertheless interesting that the data, without any direct measurements of redshifts, requires a fairly rapid evolution of the sSFR.
The halo mass scale for most efficient star formation is constrained to be around . This is consistent with the typical masses of the host halos of submillimetre galaxies, as predicted by semianalytical models (e.g. González et al., 2011). This mass scale is also in general agreement with that inferred from the clustering measurements of resolved bright sources (e.g., Cooray et al., 2010). These studies all converge and indicate that star–formation is most vigorous in halos with masses of a few . It is therefore interesting to note that a recent measurement of the smallscale CIB power by Herschel derived a much lower mass–scale of . Amblard et al. (2011b) There are two possible reasons for this difference. First, the measurements themselves might be inconsistent. A recent comparison between the two measurements indeed shows that the power spectra measured by Planck at 545 GHz and 857 GHz are higher than those measured by Herschel. Second, as mentioned in the Introduction, the meaning of in the old model (adopted by the Herschel team) does not correspond to that of the defined here, invalidating a direct comparison between the two.
The bestfit value for () is larger than the value commonly assumed (). The common assumption is due to the fact that simple models for both insulating and conducting materials naturally give at long wavelengths (Draine & Lee, 1984; Gordon, 1988). However, longrange disorder in the dust grains can lead to (Meny et al., 2007). Meny et al. (2007) also refers to observational and laboratory evidence for a wide range of values, including those greater than 2. Note that others have also inferred a steep spectral index for the background anisotropy spectrum at long wavelengths (Hall et al., 2010; Dunkley et al., 2010; Shirokoff et al., 2011).
Inferring from observations is complicated by contributions from dust grains at varying temperatures. A mixture of cold and hot dust flattens the spectrum in the wavelength range between the two spectral peaks; if analyzed with a singletemperature greybody model, the inferred value of will thus be artificially low as emphasized by Reach et al. (e.g., 1995).
Just as a model that is missing a cold component can lead to an artificially low inference of , missing a hot component can lead to an artificially high inference of . What we infer as a high value of , may be due to a component, not included in our model, with a high effective temperature where . Perhaps we are missing a significant lowredshift component.
More importantly, we emphasize that our best–fit value depends especially sensitively on model assumptions. In particular, as discussed above, is degenerate with ; we will shortly see that is reduced to if we make a different assumptions about the redshiftevolution. Another degeneracy of is with ; a possible cause of a large is that we adopted the wrong dust temperature. The value K is motivated by the temperature measurements in resolved sources, which, on average, are brighter, and reside at lower redshift, than the fainter sources responsible for the unresolved background. Finally, the measurement errors are considerable (K).
For the reasons in the preceding paragraph, we next consider the case where we fix , and allow the dust temperature to vary instead. The results are listed in Table 1 under “case 3”. The values of and are consistent with those from case 0, while the dust temperature is found to be around 45 K, 11 K higher than assumed in case 0. This numerical result confirms a rough estimate based on Figure 1: a reduction of 0.68 in could be compensated by an increase of K in . The bestfit dust temperature is then somewhat higher than the measurements quoted above (although within their 2 errors). One could indeed speculate that at the higher redshifts probed by the CIB measurement, galaxies are more compact, and dust in their interstellar medium is hotter, since dust particles reside closer to the stars.
In Figure 4, we show 2D confidence contours of the parameters in case 0. These the correlations between the 3 free parameters. In particular, is anticorrelated with and , while and are positively correlated with each other. These correlations can be understood by examining the dependence of the power spectrum on individual parameters, shown in Figure 1 above.
Finally, owing to the current uncertainty on whether the sSFR continues to rise beyond redshift (as discussed above), we consider a variation on our fiducial model with a high plateau. In this variant (referred to as “case 1”), we allow the normalization to increase with redshift as before, i.e. as a power–law with a constant (eq.22) at low redshift, but we set at . The results for this case are shown in the second row of Table 1, and by the dashed curves in Figure 3. The bestfit and both change considerably, with increasing to , and reduced close to its expected value of . This indicates that uncertainties in and due to model assumptions are large. Interestingly, the bestfit remains almost unchanged. The minimum in case 1 is somewhat larger than in case 0 (2.78 v.s. 1.82), indicating that the data favor case 0.

6 Shot noise and the dutycycle of the underlying infrared source
As mentioned in §3.2 above, our model can be used to obtain a lower limit on the shot noise. As we show here, this enables us to place a lower limit on the the dutycycle of the underlying infrared sources, , and to thus speculate on their nature (intermittent or quiescent). For simplicity, throughout the discussion below, we will assume a universal, constant duty cycle, independent of redshift and luminosity.
Combining equations (19) and (20) implies that
(28) 
where is the luminosity of “active” galaxies, and
is the number of active galaxies. For a fixed abundance of
galaxies, is proportional to the dutycycle, and if the
total luminosity density () is fixed, then
is inversely proportional to it. Thus, the amplitude of the shot noise
power is inversely proportional to the dutycycle, .
In Table 2, we list the shot noise estimates in our models, in Case 0 and Case 1 (the conclusions from the other cases are very similar) and compare these with results from the empirical model of B11. The values of needed in the calculation can be found in Table 3 of A11. The B11 shot noise estimates have been carefully calibrated with source count measurements. The model employs a double exponential function to parametrize the luminosity function, and double powerlaws to parameterize its redshift evolution. By adjusting 13 free parameters, the model successfully reproduces the counts from the midinfrared to the millimeter wavelengths, as well as the midinfrared luminosity functions. The duty cycle is taken to be either (mimicking the value for quiescent star–formation in long–lived galaxies) or (mimicking short–lived star–bursts, as expected when star–formation is triggered by major mergers). As the table shows, the shot noise levels for are of the same order as the results of B11 (although our shot noise levels for Case 0(1) are somewhat higher(lower) than those of B11). On the other hand, with , the shot noise in our models severely exceed those of B11 (by factors of ). Given this vast difference, we conclude that the average duty cycle of the underlying sources is on the order of unity. This suggests that the infrared background is mainly contributed by normal quiescent galaxies, rather than starbursts. This conclusion is consistent with recent numerical simulations (e.g., Hopkins et al., 2010), which find that of the infrared background is contributed by normal galaxies.
It is interesting to note that between the two cases, our Case 1 has a lower shot noise, due to its flatter redshift evolution. With the total intensity fixed, the shot noise decreases with increasing source number density. In Case 1, a higher fraction of the CIB originates from low redshift, where the galaxy number density is higher.
It might appear worrisome that our lower limits on the shot noise, in Case 0, are slightly higher than those of B11. This, however, could arise because our assumed is too small. Increasing would reduce the average luminosity, and consequently reduce the shot noise levels. If this were the case, then in Table 1 would have been underestimated. On the other hand, comparisons with the BLAST measurement show that the shot noise estimates by the model of B11 could be lower than the true values by a considerable amount (60% at 545 GHz, 20% at 857 GHz; see A11).
217 GHz  353 GHz  545 GHz  857 GHz  

B11    
Case 0  1  14.79  308.39  2305.94  7874.02 
0.01  516.93  5790.21  34813.30  85269.09  
Case 1  1  5.71  114.43  1057.59  5434.20 
0.01  570.56  7884.97  72353.10  89765.42 
In principle, our IR source population model could be compared
directly to measurements of the counts of resolved bright sources, as
well. The distribution of source counts as a function of flux,
especially at the bright end, also depends sensitively on the assumed
the duty cycle. In particular, if the total CIB intensity is kept
fixed, then the number of bright sources increases, while the number
of faint sources decreases with decreasing dutycycle.

7 Discussion
In this section, we discuss a few remaining issues and caveats related to our main results.
(1) Small scale clustering at 545 GHz
Although we are able to fit all four Planck channels with an overall reduced of 1.82, formally this is still a bad fit, and our bestfit model noticeably underpredicts the smallscale power at 545 GHz (by a factor of two). There are at least two possible ways to alleviate this problem.
First, the shot noise, which is almost parallel to the 1halo term over the scales considered in this study, is uncertain. The model could match the data better if the shot noise were raised by a suitable amount. Indeed, as we mentioned previously, the shot noise levels measured by BLAST (Viero et al., 2009) are about and higher than the model predictions by B11 at 545 GHz and 857 GHz, respectively (A11). If we keep the best–fit model parameters in Case 0 as listed in the first row of Table 1, but increase the shot noise at 545 GHz and 857 GHz by these factors, the tension between the model and data at the small scales is significantly reduced.
Second, one can certainly attempt to modify and/or finetune our best–fit models, to increase the small–scale power. As we mentioned in §4, this could be achieved by increasing . However, we have found no simple way to increase the small scale clustering at 545 GHz, without affecting other frequencies, especially at 857 GHz. In particular, this is because there is considerable overlap between different frequencies, caused by the extended range of redshifts over which sources contribute to the CIB. To illustrate this inter–dependence of frequencies, we perform a fit to the three lowest frequency bands, ignoring the 857GHz channel (denoted as “Case 2” in the third row of Table 1, and shown with short dashed curves in Figure 3). The best fit to these three channels yields a of 27 for 24 degrees of freedom. The power at 857 GHz in this case is, however, overpredicted by a factor of two. We conclude that nontrivial modifications to the SED, or other ingredients of our model would be needed in order to improve this situation.
(2) The normalization of the relation
Thus far, the normalization of the has been fixed by equation (26). This makes the MCMC chains converge more easily, and yields tighter constraints on other parameters. However, fixing with infinite precision could raise two possible concerns. First, while equation (26) guarantees that our model is on average consistent with the existing FIRAS measurement of the mean CIB, individual frequencies might still deviate, with a few frequencies far below the measured values, while the others far above. Second, since the uncertainties of the FIRAS measurement are nonnegligible, could our conclusions be altered if the normalization is allowed to vary?
To address the first concern, in Figure 6 we show the spectrum of the mean CIB, as predicted in our model (for Case 0 only, as the spectra in the other cases are very similar). This figure also shows the FIRAS measurement, along with the region (as derived from the errors on the parameters in the FIRAS fitting formula; Gispert, Lagache & Puget 2000). Clearly, the intensities at all frequencies are consistent with the observations.
To address the second concern, we relaxed the normalization in the MonteCarlo fitting. Other settings were kept the same as in Case 3, except that we imposed the prior , in addition to the prior on . This was necessary because otherwise the fits allowed combinations of unphysically high with suitably low normalizations. With the addition of , we then have a total of four free parameters. Only the Planck data was used in the fitting, although in principle, the FIRAS data could be used explicitly for the normalization, as well. The marginalized errors on the parameters are shown in Table 1, marked “Case 4”. The best–fit value for the normalization is found to be , where is the normalization satisfying equation (26). Compared to other cases, is considerably higher and well beyond the halo mass of a typical galaxy. The large raises the small scale clustering, and helps alleviate the small scale tension we discussed previously. It also increases the large–scale clustering (by a smaller factor), which is then compensated by the lower normalization . The negative correlation between and is also clear from the confidence contours shown in Figure 7. Interestingly, the best–fit normalization, , is within the allowed range of the FIRAS measurement, despite the fact that this was not imposed on the fits.


(3) Future improvements of the model
Although successfully describing the Planck data, and an improvement over previous version, our model is still admittedly very simplified. In the following, we give a few examples how the model could be improved and extended to include more physical ingredients.
1) As mentioned above, we neglect any scatter in the relation. This scatter could be modeled, and would impact the estimates of the shot noise and the duty–cycle, in particular. The scatter itself might be mass– and redshift–dependent.
2) We have a single galaxy population in our model. In reality, multiple populations, such as normal blue galaxies, star–bursts and obscured AGNs may all make appreciable contributions to the CIB. These populations might have different SEDs, duty cycles, relations and redshift evolutions. The fraction of starbursts and AGNs could be taken from theoretical merger–tree models, calibrated by existing data on the IR galaxy populations. In addition, the various properties of each population may depend not only on redshift, but also on environment.
3) In the current model, very simple forms have been assumed for the SED and relations, and the SED was not allowed to evolve with redshift. These are features that could be improved in future work. For example, an (asymmetric) double power–law with a suitable set of parameters is probably a better representative of the true relation, though it requires more parameters.
4) Central and satellite galaxies could have been treated with different relations.
With such improvements, the model could be developed so as to properly compute the shot noise, and simultaneously fit number counts of resolved sources and the power spectra of background fluctuations.
8 Conclusions
In this study, we developed a halo–model based formalism to compute the power spectrum of the cosmic infrared background (CIB). Although previous, similar, models provided excellent fits to individual frequency bands, they were based on an unrealistic assumption that the luminosities of galaxies are independent of the host halo masses. We relaxed this assumption by incorporating the subhalo mass function, together with a more realistic relation, into our model. With these improvements, we were able to naturally resolve the tension between the high power in the CIB on small angular scales (requiring , where the number of satellite galaxies in a halo scales as ) and optical data and numerical simulations results (which favor ).
With our improved model, we were also able to simultaneously fit all four frequency bands of the Planck measurement of the CIB anisotropy spectrum, by varying only three parameters: (the describing the redshift evolution of the relation), (the halo mass–scale on which star–formation is most efficient), and (a parameter describing the graybody SED of the unresolved IR sources). We found that star formation is most efficient in halos with relatively high masses of several, in general agreement with semianalytical models of galaxy formation. We also found that the Planck data favor an increase in the normalization with redshift. Finally, by comparing the shot noise and bright source counts predicted by our models with an observationally calibrated empirical model and a direct measurement, respectively, we conclude that the dutycycle of the source is on the order of unity. This implies that the CIB is dominated by long–lived, quiescent starforming galaxies, rather than short–lived star–bursts. This conclusion should be robust, since our estimates of shot noise and bright source counts are already conservative, owing to our neglect of any scatter in the relation.
The model we presented here, with the potential to compute the shot noise and number counts, in addition to the correlated background fluctuations, provides a theoretical framework for a future joint analysis of the unresolved CIB background and resolved source counts. Such a coherent analysis will no doubt afford new insight into the relation between dust–enshrouded star–formation and dark matter halo properties.
Acknowledgments
We thank Bruce Draine, Nicholas Hall, Olivier Dore, Guilanine Lagache, Carlo Giocoli and Dan Marrone for useful discussions. We acknowledge the use of the CNSI Computer Facilities at UC Santa Barbara for majority of the numerical work. LK and SPO acknowledge NSF grants AST 0709498 and AST 0908480 respectively.
Footnotes
 pagerange: Improved Models for Cosmic Infrared Background Anisotropies: New Constraints on the IR Galaxy Population–References
 pubyear: 2009
 https://science.nrao.edu/facilities/alma
 http://blastexperiment.info/
 http://pole.uchicago.edu/
 http://www.physics.princeton.edu/act
 See http://www.rssd.esa.int/index.php?project=Planck and http://herschel.esac.esa.int
 This neglects the fact that the brightest sources have been masked. For Planck, the flux cut is large (a few hundred mJy), and we find that the inverse relation between and holds accurately as long as .
 This can be easily verified by visualizing shifting a loglog diagram “downward” and then toward higher luminosities.
References
 Addison G. E. et al., 2011, ArXiv eprints
 Ade P. A. R. et al., 2011, ArXiv eprints (A11)
 Amblard A., Cooray A., 2007, ApJ, 670, 903
 Amblard A. et al., 2011a, Nature, 470, 510
 —, 2011b, Nature, 470, 510
 —, 2010, A&A, 518, L9+
 Barnes J. E., Hernquist L. E., 1991, ApJL, 370, L65
 Béthermin M., Dole H., Lagache G., Le Borgne D., Penin A., 2011, A&A, 529, A4+ (B11)
 Birnboim Y., Dekel A., 2003, MNRAS, 345, 349
 Bouché N. et al., 2010, ApJ, 718, 1001
 Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
 Chapman S. C., Blain A. W., Smail I., Ivison R. J., 2005, ApJ, 622, 772
 Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201
 Cooray A. et al., 2010, A&A, 518, L22+
 Cooray A., Milosavljević M., 2005, ApJL, 627, L89
 Cooray A., Sheth R., 2002, Physics Reports, 372, 1
 Croton D. J. et al., 2006, MNRAS, 365, 11
 Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
 Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785
 Draine B. T., Lee H. M., 1984, ApJ, 285, 89
 Dunkley J. et al., 2010, ArXiv eprints
 Dunne L., Eales S., Edmunds M., Ivison R., Alexander P., Clements D. L., 2000, MNRAS, 315, 115
 Dunne L. et al., 2009, MNRAS, 394, 3
 Dwek E. et al., 1998, ApJ, 508, 106
 Fixsen D. J., Dwek E., Mather J. C., Bennett C. L., Shafer R. A., 1998, ApJ, 508, 123
 Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS, 355, 819
 Gispert R., Lagache G., Puget J. L., 2000, A&A, 360, 1
 González J. E., Lacey C. G., Baugh C. M., Frenk C. S., 2011, MNRAS, 413, 749
 Gordon M. A., 1988, ApJ, 331, 509
 Grossan B., Smoot G. F., 2007, A&A, 474, 731
 Guo Q., White S., Li C., BoylanKolchin M., 2010, MNRAS, 404, 1111
 Haiman Z., Knox L., 2000, ApJ, 530, 124
 Hall N. R. et al., 2010, ApJ, 718, 632
 Hansen S. M., Sheldon E. S., Wechsler R. H., Koester B. P., 2009, ApJ, 699, 1333
 Hopkins P. F., Younger J. D., Hayward C. C., Narayanan D., Hernquist L., 2010, MNRAS, 402, 1693
 Kennicutt, Jr. R. C., 1998, ARA&A, 36, 189
 Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
 Knox L., Cooray A., Eisenstein D., Haiman Z., 2001, ApJ, 550, 7
 Komatsu E. et al., 2011, ApJS, 192, 18
 Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S., Allgood B., Primack J. R., 2004, ApJ, 609, 35
 Lagache G., Bavouzet N., FernandezConde N., Ponthieu N., Rodet T., Dole H., MivilleDeschênes M.A., Puget J.L., 2007, ApJL, 665, L89
 Lagache G., Puget J. L., 2000, A&A, 355, 17
 Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
 Limber D. N., 1954, ApJ, 119, 655
 Mathis J. S., Whiffen G., 1989, ApJ, 341, 808
 Matsuhara H. et al., 2000, A&A, 361, 407
 Meny C., Gromov V., Boudet N., Bernard J.P., Paradis D., Nayral C., 2007, A&A, 468, 171
 Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557
 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Neistein E., Dekel A., 2008, MNRAS, 383, 615
 Neistein E., Weinmann S. M., Li C., BoylanKolchin M., 2011, MNRAS, 414, 1405
 Noeske K. G. et al., 2007, ApJL, 660, L43
 Oliver S. et al., 2010a, MNRAS, 405, 2279
 Oliver S. J. et al., 2010b, A&A, 518, L21+
 Reach W. T. et al., 1995, ApJ, 451, 188
 Rodighiero G. et al., 2010, A&A, 518, L25+
 Sanders D. B., Soifer B. T., Elias J. H., Madore B. F., Matthews K., Neugebauer G., Scoville N. Z., 1988, ApJ, 325, 74
 Schaerer D., de Barros S., 2010, A&A, 515, A73+
 Sheth R. K., 2005, MNRAS, 364, 796
 Shirokoff E. et al., 2011, ApJ, 736, 61
 Skibba R., Sheth R. K., Connolly A. J., Scranton R., 2006, MNRAS, 369, 68
 Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
 Tinker J. L., Wechsler R. H., Zheng Z., 2010, ApJ, 709, 67
 Tinker J. L., Wetzel A. R., 2010, ApJ, 719, 88
 Vale A., Ostriker J. P., 2006, MNRAS, 371, 1173
 Viero M. P. et al., 2009, ApJ, 707, 1766
 Wang L., Li C., Kauffmann G., De Lucia G., 2006, MNRAS, 371, 537
 Weinmann S. M., Neistein E., Dekel A., 2011, ArXiv eprints
 Wetzel A. R., White M., 2010, MNRAS, 403, 1072
 Yabe K., Ohta K., Iwata I., Sawicki M., Tamura N., Akiyama M., Aoki K., 2009, ApJ, 693, 507
 Yang X., Mo H. J., Jing Y. P., van den Bosch F. C., 2005, MNRAS, 358, 217
 Yang X., Mo H. J., van den Bosch F. C., 2003, MNRAS, 339, 1057
 Zheng Z. et al., 2005, ApJ, 633, 791