© ESO 2018

1 Introduction

How and when the galaxies assembled their stars are two of the most pressing questions of modern cosmology. In the last decade, it has become clear that the dusty star-forming galaxies are significant contributors to this process as they pin down the episodes of galaxy-scale star formation (e.g. Gispert et al. 2000; Casey et al. 2014). Indeed they are a critical player in the assembly of stellar mass and the evolution of massive galaxies (e.g. Béthermin et al. 2015).

Dusty star-forming galaxies are difficult to detect individually at high redshift because they are so faint and numerous compared to the angular resolution achievable in the far-infrared to millimetre wavelengths that the confusion limits the ultimate sensitivity (e.g. Nguyen et al. 2010). As a result, single-dish experiments, such as Herschel and Planck, can only see the brightest objects that represent the tip of the iceberg in terms of star formation rates and halo masses. These experiments are nevertheless sensitive enough to measure the cosmic infrared background (CIB), the cumulative infrared emission from all galaxies throughout cosmic history (Lagache et al. 2005; Planck Collaboration XVIII 2011). The mean value of the CIB gives the amount of energy released by the star formation that has been reprocessed by dust (Dole et al. 2006), while CIB anisotropies trace the large-scale distribution of dusty star-forming galaxies, and can thus be used to trace the underlying distribution of the dark matter halos in which such galaxies reside (Amblard et al. 2011; Viero et al. 2013; Béthermin et al. 2013; Planck Collaboration XXX 2014).

The cosmic history of star formation is one of the most fundamental observables in astrophysical cosmology. A consistent picture has emerged in recent years, whereby the star-formation rate density peaks at z ~ 2, and declines exponentially at lower redshift. The Universe was much more active in forming stars in the past; the star formation rate density was about ten times higher at z ~ 2 than is seen today (Madau & Dickinson 2014). These highest star formation rate densities are the consequence of the rising contribution from the dusty star formation. At higher redshift (z ≫ 2), the level of dust-obscured star formation density is still a matter of debate; there is no census for the star formation rate density (SFRD) selected from dust emission alone (e.g. Koprowski et al. 2017). With their unmatched redshift depth, CIB anisotropies can be used to improve our understanding of early star formation. The contribution of dusty galaxies to the SFRD can be derived from their mean emissivity per comoving unit volume as derived from CIB anisotropy modelling. This has been investigated by Planck Collaboration XXX (2014). In this paper, we revisit this determination using Planck CIB measurements combined with the latest observational constraints and theoretical developments.

Galaxy clustering is another piece of information embedded into the CIB anisotropies (Knox et al. 2001; Lagache et al. 2007). The distribution of galaxies is known to be biased relative to that of the dark matter. In the simplest model, on large scales, it is assumed that the galaxy and dark matter mass–density fluctuations are related through a constant bias factor b such that ${(\Delta \rho /\rho )}_{\text{light}}=b{(\Delta \rho /\rho )}_{\text{mass}}$. Light traces mass exactly if b = 1 (which is the case in the local Universe, on average and at scales larger than the individual galaxy halos); if b > 1, light is a biased tracer of mass, as expected if the galaxies form in the highest peaks of the density field. The galaxy bias is dependent on the host dark matter halo mass and the redshift (Mo & White 1996). As expected by theory (e.g. Kaiser 1986; Wechsler et al. 1998) and confirmed by observations (e.g. Adelberger et al. 2005), the biasing becomes more pronounced at high redshift. As most of the star formation over the history of the Universe has been obscured bydust (Dole et al. 2006), studying the clustering of dusty galaxies is crucial to exploring the link between dark matter halo mass and star formation. Measurements of correlation function of individually resolved galaxies, whilst not forming the dominant contribution to the CIB, give constraints on the host halo mass of the brightest star-forming objects, the submillimetre galaxies (SMGs). Maddox et al. (2010) and Cooray et al. (2010) measured the correlation function of Herschel/SPIRE galaxies (see also Weiß et al. 2009; Williams et al. 2011 for measurements from ground-based observations), but their results were inconsistent probably because of the difficulty in building well-controlled selections in low angular resolution data (Cowley et al. 2016). One way to get more accurate measurements is to use the cross-correlation with optical/near-IR samples (Hickox et al. 2012; Béthermin et al. 2014; Wilkinson et al. 2017). These analyses found the host halo mass of galaxies with the star formation rate (SFR) > 100 M_⊙ to be ~ 10^12.5−13 M_⊙. Reaching the population that contributes to the bulk of the CIB requires modelling the angular power spectrum of CIB. Making various assumptions such as the form of the relationship between galaxy luminosity and halo mass, a typical halo mass for galaxies that dominate the CIB power spectrum has been found to lie between 10^12.1±0.5 M_⊙ and 10^12.6±0.1 M_⊙ (Viero et al. 2013; Planck Collaboration XXX 2014). In this paper, we exploit the strength of Planck, which is a unique probe of the large-scale anisotropies of the CIB. We use only the linear part of the power spectra, and thus consider a simple model (with few parameters), to measure the galaxy clustering of dusty star-forming galaxies up to z ~ 5.

All previous CIB angular power spectrum analyses have considered a fixed cosmology. Assuming a cosmology is a fair assumption as the model of the CIB galaxy properties has far greater uncertainties than on the cosmological parameters (Planck Collaboration XIII 2016). However, the CIB power spectrum depends on key cosmological ingredients, such as the dark matter power spectrum, cosmological distances, and cosmological energy densities. Therefore, it is worth studying the impact of the cosmology on CIB modelling parameters, which we do here for the first time.

The paper is organised as follows. We describe the CIB anisotropy modelling and the data we use to constrain the model in Sect. 2. In Sect. 3, we present and discuss the redshift distribution of CIB anisotropies and mean level, and the star formation rate density from z = 0 to z = 6 that are obtained by fitting our model to the CIB and CIB × CMB lensing power spectra measurements and by using extra constraints coming from galaxy number counts and luminosity functions. We then discuss the host dark matter halo mass of CIB galaxies in Sect. 4 by first measuring the effective bias b_eff, then converting it to the mean bias of halos hosting the dusty galaxies, before finally computing the host dark matter halo mass of the dusty star formation. In Sect. 5, we study the effect of cosmology on the CIB. We conclude in Sect. 6.

Throughout the paper, we used a Chabrier mass function (Chabrier 2003) and the Planck 2015 flat Λ CDM cosmology (Planck Collaboration XIII 2016) with Ω_m = 0.33 and H₀ = 67.47 km s⁻¹ Mpc⁻¹.

2 CIB modelling on large scales

2.1 CIB power spectrum

The measured CIB intensity I_ν at a given frequency ν, in a flat Universe, is given by $$I_{\nu }={\displaystyle \int \frac{\text{d}\chi }{\text{d}z}}aj(\nu ,z)\text{d}z,$$(1)

where j is the comoving emissivity, χ(z) is the comoving distance to redshift z, and a = 1∕(1 + z) is the scale factor. The angular power spectrum of CIB anisotropies is defined as $$C_l^{\nu \times {\nu }^{\prime }}\times {\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}=\langle \delta I_{lm}^{\nu }\delta I_{l^{\prime }{m}^{\prime }}^{{\nu }^{\prime }}\rangle .$$(2)

Combining these two equations and using the Limber approximation, we get $$C_l^{\nu \times {\nu }^{\prime }}={\displaystyle \int \frac{\text{d}z}{{\chi }^2}}\frac{\text{d}\chi }{\text{d}z}{a}^2\overline{j}(\nu ,z)\overline{j}({\nu }^{\prime },z)P_j^{\nu \times {\nu }^{\prime }}(k=l/\chi ,z),$$(3)

where $P_j^{\nu \times {\nu }^{\prime }}$ is the 3D power spectrum of the emissivities and is defined as $$\langle \delta j(k,\nu )\delta j(k^{\prime },{\nu }^{\prime })\rangle ={(2\pi )}^3\overline{j}(\nu )\overline{j}({\nu }^{\prime })P_j^{\nu \times {\nu }^{\prime }}(k){\delta }^3(k-k^{\prime }),$$(4)

with δj being the emissivity fluctuation of the CIB. In our analysis, we are interested in modelling the CIB anisotropies on large angular scales (ℓ < 800), where the clustering dominated by the correlation between dark matter halos and the non-linear effects can be neglected. We can then model the CIB anisotropies equating P_j with P_gg, which is the galaxy power spectrum. On large angular scales, $P_{gg}(k,z)=b_{\text{eff}}^2(z)P_{\text{lin}}(k,z)$. Here, b_eff is the effective bias factor for dusty galaxies at a given redshift, i.e. the mean bias of dark matter halos hosting dusty galaxies at a given redshift weighted by their contribution to the emissivities (see Planck Collaboration XXX 2014). The fact that more massive halos are more clustered and have host galaxies emitting more far-infrared emission is implicitly taken into account through this term. We consider b_eff as scale independent, which is a good approximation at the scales we are interested in. Here, P_lin (k, z) is the linear theory dark matter power spectrum. Finally, the linear CIB power spectrum is given as $$C_l^{\nu \times {\nu }^{\prime }}={\displaystyle \int \frac{\text{d}z}{{\chi }^2}}\frac{\text{d}\chi }{dz}{a}^2{b}_{\text{eff}}^2\overline{j}(\nu ,z)\overline{j}({\nu }^{\prime },z)P_{\text{lin}}(k=l/\chi ,z).$$(5)

This is known as the linear model of CIB anisotropies. The emissivities $\overline{j}(\nu ,z)$ (in Jy Mpc⁻¹) are derived using the star formation rate density ρ_SFR (M_⊙ yr⁻¹ Mpc⁻³) following $$\overline{j}(\nu ,z)=\frac{{\rho }_{\text{SFR}}(z)(1+z)S_{\nu ,\text{eff}}(z){\chi }^2(z)}{K},$$(6)

where K is the Kennicutt constant (Kennicutt 1998) and S_ν,eff(z) is the mean effective spectral energy distribution (SED) for all the dusty galaxies at a given redshift. We used a Kennicutt constant corresponding to a Chabrier initial mass function (SFR ∕L_IR = 1.0 × 10⁻¹⁰ M_⊙ yr⁻¹). Here, S_ν,eff (z) (in Jy L_⊙) is the average SED of all galaxies at a given redshift weighted by their contribution to the CIB. They have been computed using the method presented in Béthermin et al. (2013), but assuming the new updated SEDs calibrated with the Herschel data and presented in Béthermin et al. (2015, 2017). Compared to the previous SEDs used in Béthermin et al. (2013) and Planck Collaboration XXX (2014), the dust in these new SED templates is warmer at z > 2. We used CAMB¹ to generate cold dark matter power spectra P_lin(k) for given redshifts. For the effective bias, we chose the following parametric form based on redshift evolution of the dark matter halo bias: $$b_{\text{eff}}(z)=b_0+b_{1}z+b_2{z}^2.$$(7)

We also tested alternative parametric forms for the effective bias where the evolution of the bias is slower compared to the above equation. The results are presented in Appendix A.

To describe ρ_SFR, we used the parametric form of the cosmic star formation rate density proposed by Madau & Dickinson (2014) and given by $${\rho }_{\text{SFR}}(z)=\alpha \frac{{(1+z)}^{\beta }}{1+{[(1+z)/\gamma ]}^{\delta }}{M}_{\odot }{\text{yr}}^{-1}{\text{Mpc}}^{-3},$$(8)

where α, β, γ, and δ are free parameters in our CIB model.

2.2 CIB–CMB lensing

Large-scale distribution of the matter in the Universe gravitationally deflects the cosmic microwave background (CMB) photons which are propagating freely from the last scattering surface. This gravitational lensing leaves imprints on the temperature and polarisation anisotropies. These imprints can be used to reconstruct a map of the lensing potential along the line of sight (Okamoto & Hu 2003). Dark matter halos located between us and the last scattering surface are the primary sources for this CMB lensing potential (Lewis & Challinor 2006) and it has been shown (e.g. Song et al. 2003)that a strong correlation between CIB anisotropies and a lensing derived projected mass map is expected.

We calculated the cross-correlation between the CIB and CMB lensing potential (Planck Collaboration XVIII 2014), which is given by $$C_l^{\nu \varphi }={\displaystyle \int b_{\text{eff}}}\overline{j}(\nu ,z)\frac{3}{l^2}{\Omega }_m{H}_0^2\left(\frac{{\chi }_{*}-\chi }{{\chi }_{*}\chi }\right)P_{\text{lin}}(k=l/\chi ,z)\text{d}\chi ,$$(9)

where χ_* is the comoving distance to the CMB last scattering surface, Ω_m is the matter density parameters, and H₀ is the valueof the Hubble parameter today. From Eq. (9), we see that $C_l^{\nu ,\varphi }$ is proportional to b_eff, whereas $C_l^{\text{CIB}}$ is proportional to $b_{\text{eff}}^2$. Therefore, also using the CIB–CMB lensing potential measurement when fitting the CIB model helps us resolve the degeneracy between the evolution of b_eff and ρ_SFR to some extent.

2.3 Constraints on the model through data

2.3.1 Observational constraints on the power spectra

We used the CIB angular power spectra measured by Planck Collaboration XXX (2014). The measurements were obtained by cleaning the Planck frequency maps from the CMB and galactic dust, and they were further corrected for SZ and spurious CIB contamination induced by the CMB template, as discussed in Planck Collaboration XXX (2014). We used the measurements at the four highest frequencies (217, 353, 545, and 857 GHz from the HFI instrument). For the 3000 GHz, far-infrared data from IRAS (IRIS, Miville-Deschênes & Lagache 2005) were used. As we are using the linear model we fit the data points only for l ≤ 600, which are dominated by the 2-halo term (>90%, Béthermin et al. 2013; Planck Collaboration XXX 2014).

CIB–CMB lensing potential cross-correlation values and error bars are available for the six Planck HFI channels (100, 143, 217, 353, 545, and 857 GHz) and are provided in Planck Collaboration XVIII (2014). These values range from ℓ = 163 to ℓ = 1937. As discussed in Planck Collaboration XVIII (2014), the non-linear term can be neglected in this range of multipoles.

We usedthe νI_ν = const photometricconvention. Thus, the power spectra computed by the model need to be colour-corrected from our CIB SEDs to this convention. Colour corrections are 1.076, 1.017, 1.119, 1.097, 1.068, 0.995, and 0.960 at 100, 143, 217, 353, 545, 857, and 3000 GHz, respectively (Planck Collaboration XXX 2014). The CIB power spectra are then corrected as $$C_{l,\nu ,{\nu }^{\prime }}^{\text{model}}\text{\hspace{0.22em}x\hspace{0.22em}}{\text{cc}}_{\nu }\text{\hspace{0.22em}x\hspace{0.22em}}{\text{cc}}_{{\nu }^{\prime }}=C_{l,\nu ,{\nu }^{\prime }}^{\text{measured}}.$$(10)

Calibration uncertainties are not accounted for in the CIB power spectra error bars and are treated differently. Béthermin et al. (2011) introduce a calibrationfactor $f_{\text{cal}}^{\nu }$ for galaxy number counts. We used a similar approach here. We put Gaussian priors on these calibration factors for different frequency channels with an initial value of 1 and the error bars as given in Table 1. The CIB measurements from Planck wereobtained with the PR1 release of the maps. Between the PR1 and PR2 releases, the absolute calibration improved further, and we thus wanted to consider the absolute calibration of the PR2 release for the CIB measurements. Consequently, the CIB × CIB and CIB × CMB lensing power spectra data points were corrected for the absolute calibration difference of the two releases (following the numbers given in Table 1), and we used the absolute calibration errors as given for the PR2 release (Planck Collaboration VIII 2016).

2.3.2 External observational constraints

To put better constraints on ρ_SFR and b_eff parameters, in addition to the CIB × CIB and CIB × CMB lensing angular power spectra, we used some external observational constraints on the star formation rate density at different redshifts, local bias of the dusty galaxies, and mean CIB levels at different frequencies, as detailed below:

• 1.

  We used the ρ_SFR measurements at different redshifts that were obtained by measuring the IR luminosity functions from Gruppioni et al. (2013), Magnelli et al. (2013), and Marchetti et al. (2015) (see the discussion in Sect. 3.2). This helped us to put better constraints on the parameters. The cosmological parameters used in these studies are different from the ones used here. As mentioned before, we also perform a Fisher matrix analysis to study the effect of the CIB parameters on cosmology. For this purpose, we needed to convert all the observational ρ_SFR data points to actual measurements which are cosmology independent. We thus used the observed flux in the range 8–1000 μm (rest frame) per redshift bin per solid angle $\frac{{\text{dB}}_{\text{IR}}}{\text{d}z\text{d}\Omega }$. We perform the conversion as follows:

$$\begin{array}{lll}\frac{\text{d}{B}_{\text{IR}}}{\text{d}z\text{d}\Omega }\hfill & =\hfill & \frac{\text{d}(P_{\text{IR}}/4\pi D_L^2)}{\text{d}z\text{d}\Omega }\hfill \\ \hfill & =\hfill & \frac{1}{4\pi D_L^2}\frac{\text{d}{P}_{\text{IR}}}{\text{d}{V}_c}\frac{\text{d}{V}_c}{\text{d}z\text{d}\Omega }\hfill \\ \hfill & =\hfill & \frac{1}{4\pi D_L^2}\frac{{\rho }_{\text{SFR}}}{K}\frac{D_H{(1+z)}^2{D}_A^2}{E(z)}.\hfill \end{array}$$

Here $\frac{\text{d}{P}_{\text{IR}}}{\text{d}{V}_c}$ is the power emitted in the infrared per unit co-moving volume of space, D_L is the luminosity distance, D_A is the angular diameter distance, D_H is the Hubble distance given by c∕H₀ (with c and H₀ being the speed of light and Hubble’s constant, respectively), and $E(z)=\sqrt{{\Omega }_m{(1+z)}^3+{\Omega }_k{(1+z)}^2+{\Omega }_{\Lambda }}$. This equation can be simplified further using the relation $D_L={(1+z)}^2{D}_A$ and becomes $$\frac{\text{d}{B}_{\text{IR}}}{\text{d}z\text{d}\Omega }=\frac{1}{4\pi }\frac{{\rho }_{\text{SFR}}}{K}\frac{D_H{(1+z)}^2}{{(1+z)}^{2}E(z)}.$$(14)

We also needed to convert the measurements to the same IMF; e.g. Gruppioni et al. (2013), Magnelli et al. (2013), and Marchetti et al. (2015) used the Salpeter IMF, whereas we used the Chabrier IMF. Finally, we used these observed flux values d B_IR∕dzdΩ in our fitting rather than the ρ_SFR values. This conversion to cosmology-independent variables is particularly important for the Fisher matrix analysis presented in Sect. 5, else we would wrongly find that CIB anisotropies can be used to significantly improve the uncertainties on the cosmological parameters.

• 2.

  Saunders et al. (1992) provide the fixed value for the product of the local bias of infrared galaxies b_I and the cosmological parameter σ₈: b_Iσ₈ = 0.69 ± 0.09. We put aconstraint of b = 0.83 ± 0.11 on local bias value which has been converted using the σ₈ measured by Planck Collaboration XIII (2016).

• 3.

  The mean level of the CIB at different frequencies has been deduced from galaxy number counts. The values we used are given in Table 2. Similarly to what is being done for $C_{l,\nu ,{\nu }^{\prime }}$ (Eq. (10)), the CIB computed by the model needs to be colour-corrected, from our CIB SEDs to the photometric convention of PACS, SPIRE, and SCUBA2. The colour corrections are computed using the Béthermin et al. (2012a) CIB model and the different bandpasses and are given in Table 2.

• 4.

  Finally, we also put physical constraints on ρ_SFR and b_eff parameters such that the star formation rate and galaxy bias is positive at all redshifts.

Fig. 1 Posterior confidence ellipses for the CIB model parameters. All the calibration parameters lie within the 1σ range of their prior values. There is no significant degeneracy between any of the calibration parameters and CIB parameters and hence they are not shown here.

2.3.3 Fitting the data

We performed a Markov chain Monte Carlo (MCMC) analysis on the global CIB parameter space. Python package “emcee” (Foreman-Mackey et al. 2013) is used for this purpose. We have a 12-dimensional parameter space: $\{\alpha ,\beta ,\gamma ,\delta ,b_0,b_1,b_2,f_{217}^{\text{cal}},f_{353}^{\text{cal}},f_{545}^{\text{cal}},f_{857}^{\text{cal}},f_{3000}^{\text{cal}}\}$. Our global χ² has a contribution from the CIB × CIB and CIB × CMB lensing angular power spectra measurements, priors on calibration factors, and priors imposed by the external observational constraints mentioned above. We assumed Gaussian uncorrelated error bars for measurement uncertainties. The 1-halo and Poisson terms have very little measurable contribution (less than 10%) at l ≤ 600, where we are fitting the data points. Similar to the procedure used in Planck Collaboration XXX (2014), we added the contributionof these terms derived from the Béthermin et al. (2013) model to our linear model.

2.4 Results

We present the results from our fit in Table 3. We find a good fit for the model and the posterior of all the parameters with a Gaussian prior (local effective bias and calibration factors) are within a 1σ range of the prior values. The 1σ and 2σ confidence regions for the ρ_SFR and b_eff parameters are shown in Fig. 1. As expected, we observe strong degeneracies between ρ_SFR and b_eff, and also within the ρ_SFR and b_eff parameters themselves. Figures 2 and 3 show the fit for the linear model to the observational data points for all CIB × CIB auto- and cross-power spectra, and CIB × CMB lensing power spectra, respectively. We also show the shot noise and one-halo term for all the frequencies in Fig. 2.

In Fig. 3, we show the comparison between our best fit model and the CIB–CMB lensing cross-correlation points. We consider these data points to calculate the best fit of our model. We find a good agreement between the data points and the best fit.

3 CIB redshift distribution and star formation history

3.1 Redshift distribution of CIB anisotropies and mean level

The model we used is based on the SEDs over a range of wavelengths for galaxies which vary with redshift. It is thus interesting to study the redshift distribution of the mean level of the CIB as well as the CIB anisotropies which are based on these SEDs. In Fig. 4, we show the redshift distribution of the CIB anisotropies at ℓ = 300 and CIB mean level intensity for different frequency bands. Both plots have been normalised, i.e. ∫ d(νI_ν)∕dz dz = 1 and ∫ dC_l∕dz dz = 1, so that it is easier to compare the results for different frequencies.

It is observed that as we go to lower frequencies, from 3000 to 217 GHz, the peak of the CIB anisotropy redshift distribution increases in redshift. This is expected as higher frequencies (lower wavelengths) probe lower redshifts and vice versa (e.g. Lagache et al. 2005; Béthermin et al. 2013). The CIB mean level distribution follows the same trend. A similar pattern is observed for the CIB anisotropy and for the mean level distribution in Béthermin et al. (2013) model. We note, however, that the peak of CIB anisotropies predicted in Béthermin et al. (2013) at 353 and 217 GHz is reached around z = 2.4 and the same is reached at lower redshift (around z = 1.7) for our model.

In Fig. 5, we compare our redshift distribution on the CIB mean level with the lower limits from Béthermin et al. (2012b) and Viero et al. (2013) which were derived by stacking the 24 μm-selected and mass-selected sources, respectively. It is observed that the CIB mean level from our model is higher than these two lower limits for most of the values. We also compare our results with redshift distribution of the CIB derived by Schmidt et al. (2015) using the cross-correlation between the Planck HFI maps and SDSS DR7 quasars. It can be seen that although some of the data points from Schmidt et al. (2015) are consistent with our curve, the measurements tend to be higher than our model at Z > 0.5. A similar trend is followed at other frequencies, and so they are not shown here. Koprowski et al. (2017) assumed a Gaussian distribution for the IR and UV ρ_SFR and also provided the best fit values for the distribution. Based on this ρ_SFR form, we calculated the corresponding mean CIB level distribution using Eqs. (6) and (5) using the SED templates from Béthermin et al. (2017) (see Fig. 5, dashed black curve). It can be seen that their mean CIB level distribution is lower than some lower limits (at z ~ 1 and z ~ 3). It is also lower than values from our model between 0.5 < z < 2, and hence in even more tension with the values from Schmidt et al. (2015).

Fig. 2 Measurements of the CIB auto- and cross-power spectra obtained by Planck and IRAS (extracted from Planck Collaboration XXX 2014) and the best fit CIB linear model. Data points on higher ℓ are not shown as they are not used for fitting. We also show the shot noise (black dash-dotted line) and one-halo term (red dashed line) for all the frequencies.

3.2 Star formation history

We show the evolution of star formation density with the redshift in Fig. 6 with the best fit as well as the corresponding 1σ and 2σ confidence regions. We created these confidence regions using the chains from the Monte Carlo sampling of our likelihood. We took random samples from the chains for the star formation density parameters and constructed an array of ρ_SFR with these samples for all the redshifts. The median of these values at every redshift is then the central value, and we took samples within 68.2% and 95.4% around the central value as 1σ and 2σ regions, respectively.

We also show the IR measurements from Gruppioni et al. (2013), Magnelli et al. (2013), Marchetti et al. (2015), and Bourne et al. (2017), and UV measurements from Cucciati et al. (2012), Bouwens et al. (2012), and Reddy & Steidel (2009), in the upper and lower panels, respectively. Our best fit passes through the IR data points considered in the fit. As mentioned in Sect. 2.3.2, these data points were obtained using a Salpeter IMF, whereas in our study, we use a Chabrier IMF. Therefore, both the IR and UV data points have been converted to take into account this change.

As expected, the UV data points (which have not been corrected for dust attenuation) for unobscured star formation density lie below our curves for z < 3 as most of the UV light emitted by young and short-lived stars in this regime is reprocessed by dust. Although the IR contribution is still dominant below z < 4, the contribution from the UV becomes significant and roughly equal to IR for z > 4. Hence IR contribution alone is not a good measure of the total star formation rate density at such high redshifts.

We plot the best fit IR and UV ρ_SFR curves from Koprowski et al. (2017) with dashed black lines in the upper and lower panels of Fig. 6, respectively. It can be seen that their value of the IR SFRD in the local Universe is higher compared to the other measurements. Also, their IR SFRD drops very quickly at higher redshifts and is basically negligible for z > 4. Their results are clearly discrepant with ours in these two regimes. A similar trend is followed by the UV SFRD which drops quickly and is much smaller than current observational constraints at z ~ 6.

To investigate the discrepancy between our values and the Koprowski et al. (2017) ρ_SFR, we performed a MCMC fitting for the CIB anisotropy model (as described in Sect. 2.3.3) and fitting for the effective bias and thecalibration factors, but fixing the SFRD to the best fit of Koprowski et al. (2017). We found the effective bias at redshifts z > 2.5 to be much higher (>50%) than is observed for the SMGs. This effective bias value is too high compared to what is realistically expected, suggesting that their SFRD is underestimated.

Similarly, Fig. 6 shows that the Bourne et al. (2017) points at z ≃ 1 and 2 are a factor >2 and 1.3 below our measurements, respectively. Taking these points as priors for ρ_SFR rather than those from Gruppioni et al. (2013) would give a lower measurement for SFRD and this would affect the measurement of the linear bias. We investigated the consistency of the solution obtained on both ρ_SFR and b_eff either by taking the Bourne et al. (2017) data points as priors or by forcing and fixing the SFRD to go through the Bourne et al. (2017) data points in our MCMC analysis. We found that the linear bias is severely overestimated compared to similar galaxy populations, for example by a factor >3 at z ≃ 1.2 compared to the SMG sample of Wilkinson et al. (2017). This shows that the SFRD determined by Bourne et al. (2017) at low z seems to be underestimated.

Fig. 3 CIB × CMB lensing cross-power spectra. Our best fit model is shown by the blue curves for different frequency channels. Measurements (red data points) are from Planck Collaboration XVIII (2014). They have been included in the likelihood to calculate and get better constraints on the CIB linear model.

Fig. 4 Expected CIB mean level and anisotropy redshift distributions are shown in the lower and upper panels, respectively. Both of the distributions have been normalised for easier comparison, and the different colours show different frequency bands. The CIB anisotropy distribution is shown for ℓ = 300.

Fig. 5 Expected CIB mean level redshift distribution at 857 GHz compared to observational constraints. The lower limits are from Béthermin et al. (2012b) and Viero et al. (2013). They are shown using black and red triangles, respectively. Measurements from Schmidt et al. (2015) are shown in green. The redshift distribution derived from Koprowski et al. (2017) SFRD constraints is shown with the black dashed line.

Fig. 6 Evolution of star formation density with redshift as constrained by the linear CIB model. The ± 1σ and ± 2σ confidence region around the median realisation is shown in yellow and blue, respectively. Measurements of obscured star formation density from Gruppioni et al. (2013), Magnelli et al. (2013), and Marchetti et al. (2015), which were used to fit the CIB model, have been added along with the measurements from Bourne et al. (2017) in the upper panel. The Gaussian form of the ρSFR for the IR obscured star formation rate density used by Koprowski et al. (2017) has been plotted as the black dashed line in the upper panel with a corresponding UV part in the lower panel. Unobscured star formation rate density derived from UV from Cucciati et al. (2012), Bouwens et al. (2012), and Reddy & Steidel (2009) are also shown in the lower panel.

4 Host dark matter halos of the CIB

4.1 Effective bias of CIB galaxies

We also study the evolution of the effective bias and find that it is increasing with redshift, as expected. In Fig. 7, we compare our measurements with the clustering measurements of individual galaxies selected using various criteria.

At z > 2, our measurements and the clustering of SMGs by Wilkinson et al. (2017) are very similar. SMGs are the most star-forming galaxies at z ≳ 2 (e.g. Blain et al. 2002; Chapman et al. 2005), and thus it is not completely surprising to find a similar bias since a large fraction of the SFRD at z > 2 is hosted by galaxies forming more than 100 M_⊙ yr⁻¹ (e.g. Caputi et al. 2007; Magnelli et al. 2009; Béthermin et al. 2011; Gruppioni et al. 2013). At z < 2, Wilkinson et al. (2017) found a bias close to 1, which is 1σ and 2σ below our model in the 1 < z < 1.5 and 1.5 < z < 2 bins, respectively. This might be a statistical fluctuation coming from the small sample in their low-z bins (~60 objects). In addition, the algorithm they use to associate SMGs with optical counterparts and photometric redshift is only reliable at 83%, which could induce biases in the clustering measurements. In the 1.5 < z < 2 bin, if they consider only objects with a solid radio identification of the counterparts, they find a bias of 1.65 ± 1.09, which is compatible at 1σ with our measurements, which means that there might be no real tension in these low-z bins.

We also compare our bias estimates with selections at other wavelengths. Béthermin et al. (2014) measured the clustering of massive star-forming galaxies at z ~ 2 selected using the BzK color criterium (Daddi et al. 2004). Our effective bias agrees with this measurement. Viero et al. (2013) show that majority of the contribution to the CIB is emitted by massive dusty star-forming galaxies with log (M∕M_⊙) 10.0–11.0. We plot the clustering of the mass-selected (M_⋆ = 10^10.6 M_⊙) galaxies measured by Cowley et al. (2018) and find an agreement with our measurements, as expected. Ishikawa et al. (2015) measured the clustering of gzKs-selected galaxies for various depths and found a strong dependance of bias with the depth. Their deep sample (K ≤ 23) has a much smaller bias than the CIB (1.8 vs. 3.1), but the shallowest sample (K ≤ 21) has a stronger bias (4.16). The CIB measurements are indeed dominated by the galaxies, which contribute the most to the star formation budget, while clustering measurements of galaxy populations are dominated by the most numerous objects, which are usually the numerous low-mass objects residing in low-mass halos (M_halo < 10¹² M_⊙) with a lower clustering. It is thus expected that the deepest K-band sample with masses and SFRs well below the knees of the mass and SFR functions have a lower bias than the CIB. In a similar way, the bias of H_α -selected galaxies, which are mainly low-SFR and low-mass galaxies, measured by Cochrane et al. (2017) (1.12, 1.78, 2.52), is significantly lower than our model (1.64, 2.62, 4.06) at redshifts 0.8, 1.47, 2.23 respectively.

4.2 What does the effective bias mean?

Throughoutthe paper, we have used the effective bias term (b_eff) and it should be noted that effective bias might not be equal to the mean bias of halos hosting the dusty galaxies (b(M, z)). As shown in Appendix C of Planck Collaboration XXX (2014), the effective bias is the mean bias of halos hosting the dusty galaxies weighted by their differential contribution to the emissivities and is given by the following equation: $$b_{\text{eff}}=\frac{{\displaystyle \int \frac{\text{d}j}{\text{d}M}}(\nu ,z)b(M,z)\text{d}M}{{\displaystyle \int \frac{\text{d}j}{\text{d}M}}\text{d}M}.$$(15)

In this equation, $\frac{\text{d}j}{\text{d}M}$ is the differential contribution of a range of halo mass to the emissivity. Since the emissivity is directly proportional to ρ_SFR, we can rewrite this equation as $$b_{\text{eff}}=\frac{{\displaystyle \sum _{\text{vol}}\text{SFR}}\times b(M,z)}{{\displaystyle \sum _{\text{vol}}\text{SFR}}},$$(16)

where we sum over all the halos and their host galaxies in a given volume. Here b(M, z), which is the bias as a function of halo mass and redshift, is calculated with the formula which has been calibrated using the numerical simulations from Tinker et al. (2010), $$b(\nu )=1-A\frac{{\nu }^a}{{\nu }^a+{\delta }_c^a}+B{\nu }^b+C{\nu }^c,$$(17)

where ν = δ_c∕σ_M characterises the peak heights of the density field as function of mass and redshift; δ_c is the linear critical density at a given redshift for a given cosmology; and σ_M is the linear matter variance in a top hat filter of width $R={(3M/4\pi \overline{{\rho }_0})}^{1/3}$ with M and $\overline{{\rho }_0}$ being the mass and mean density for a given halo, respectively. Details of these calculations have been given in Appendix A of Coupon et al. (2012). A, a, B, b, C, and c are functions of Δ which is theoverdensity for a given halo defined as the mean interior density for the halo relative to the background. These functions are provided in Table 2 in Tinker et al. (2010). As often used, we took a virial value of Δ ≈ 200. From Eq. (15), we can see that if halos in a given mass range contribute more or less to the emissivity compared to halos of a different mass range, there will be a slight offset between b_eff and b(M, z).

It is interesting to quantify this offset, particularly in the context of the measurement of the mass of the typical host dark matter halos which contribute to the CIB, and hence to the obscured star formation. This is the purpose of the next section.

As a first approximation, we can estimate the mean host halo mass at a given redshift by equating b_eff and b(M, z). At z = 2, we find that the mass of the dark matter halos which contribute the most to the CIB calculated this way is ${\log }_{10}M={12.93}_{-0.081}^{+0.084}$, where 0.084 and 0.081 are the 1σ upper and lower limits, respectively.

Fig. 7 Evolution of the effective bias with redshift derived from the CIB compared to the observational values obtained on selected populations of galaxies: Béthermin et al. (2014), Cochrane et al. (2017), Cowley et al. (2018), and Wilkinson et al. (2017) using a sample of BzK-selected galaxies, H α -selected galaxies, mass-selected galaxies, and SMGs, respectively.

4.3 Host dark matter halo mass of the dusty star formation

As is shown in the previous section, b_eff might not be exactly equal to the mean bias of halos hosting the dusty galaxies (b(M, z)). To understand how this effective bias relates to the bias of the halos contributing the most to the star formation budget, we used the Simulated Infrared Dusty Extragalactic Sky (SIDES) simulation (Béthermin et al. 2017), where the impact of clusteringand angular resolution for far-infrared and millimetre continuum observations has been taken into account.

The SIDES simulation provides the star formation rate for galaxies corresponding to a given halo mass M_halo at a given redshift z. In Fig. 8 we plot the cumulative star formation rate as a function of halo masses for galaxies with 2.0 < z < 2.2. We can see from the red curve in Fig. 8 how the cumulative SFR density grows with the halo masses in SIDES. We define the mass at which the above curve reaches 50% of the total cumulative star formation as the dark matter halo mass hosting the dusty galaxies.

We calculated the bias b(M, z) corresponding to each halo at a given redshift in SIDES using Eq. (17). The blue curve in Fig. 8, shows how the normalised cumulative SFR ×b(M, z) grows with the halo mass. It is seen that 50% of the total SFR or SFR × b is reached around a halo mass of 10¹² M_⊙, but not at the exact same halo mass. Since the bias is higher for massive halos hosting more massive and star-forming galaxies, we reach these 50% for SFR × b at a mass 0.06 dex higher than for just SFR. The effective bias of the CIB, thus corresponds to a bias slightly higher than the bias of the pivot halo mass at which we have 50% of the SFRD. This value of 0.06 dex, however, is still not the offset we should apply to the mass ${\log }_{10}M={12.93}_{-0.081}^{+0.084}$ of the host dark matter halo obtained using b_eff.

The CIB traces only the obscured star formation, but at low stellar mass (and thus halo mass), a significant fraction of the UV from the star formation escapes the galaxies. We used the method of Bernhard et al. (2014) based on the empirical relation between stellar mass and dust attenuation of Heinis et al. (2014) to predict the SFR_IR, i.e. the obscured star formation, from which we deduce an L_IR (obscured IR luminosity) for all the sources in the simulation. In Fig. 8, the green and black curves represent the cumulative contributionof the L_IR and the L_IR × b, respectively. This latest quantity (L_IR × b) is the closest to what we actually measure with the CIB. The halo mass at which 50% of the total quantity is reached for slightly higher halo mass (0.05 dex) than that found for SFR × b because obscured star formation is slightly biased toward massive galaxies (and thus halos). For the dark matter halos between 2 < z < 2.2, we observe a difference of around 0.11 dex between halo masses found using only the SFR and the bias weighted LIR (L_IR × b). Using Eq. (16), we can calculate the b_eff corresponding to these LIR values where we just replace the SFR term in the equation with the LIR obtained here. Using the procedure followed in Sect. 4.2, we convert this b_eff to the corresponding halo mass using Eq. (17). The black vertical line in Fig. 8 represents this mass of the dark matter halo. Finally, the difference between this mass and halo mass found using only the SFR, which is 0.10 dex, is the actual offset we are looking for at this redshift bin.

Following the procedure mentioned above, we calculate this mass offset for all the redshift bins. This offset is shown in Fig. 9. The dashed blue line in the figure shows the mean value of the offset over all the bins which is around 0.16 dex. It is clear from the figure that this offset is not constant at all redshifts. After applying this mean correction to the halo masses at all the redshifts, we also propagate the uncertainty on this offset to the error bars on the mass measurements. We find that the mass of the dark matter halos which contribute the most to the CIB at z = 2 comes out to be ${\log }_{10}M={12.77}_{-0.125}^{+0.128}$, where 0.128 and 0.125 are the 1σ upper and lower limits, respectively (considering additional uncertainty on mass offset), compared to the previous value of ${\log }_{10}M={12.93}_{-0.081}^{+0.084}$.

We show in Fig. 10 the variation of the host halo mass as a function of redshift for z ≤ 4. Contours represent the 1σ and 2σ confidence regions. We obtain an almost constant dark matter halo mass (≈ 10^12.7 M_⊙) contributing to the CIB from 1 < z < 4, and then it starts decreasing from z < 1. Host halos grow over time and the dashed lines show this growth of halos with redshift, as computed by Fakhouri et al. (2010). We see that most of the star formation at z > 2.5 occurred in the progenitors of clusters (M_h(z = 0) > 10^13.5 M_⊙). Then at the lower redshifts from 0.3 < z < 2.5, most of the stars were formed in groups (10^12.5 < M_h(z = 0) < 10^13.5 M_⊙) and later on, finally, inside the Milky Way-like halos (10¹² < M_h(z = 0) < 10^12.5 M_⊙) for z < 0.3. Although a bit high, these results are compatible with the results obtained by Béthermin et al. (2013).

Fig. 8 Normalised cumulative SFR density and LIR density as a function of host dark matter halo mass for galaxies inSIDES simulation are shown in red and green, respectively. We also show the variation of these quantities when they are weighted with the corresponding bias value for the galaxies in blue and black, respectively. The black vertical line represents the mass of the dark matter halo mass calculated converting the beff obtained using the LIR instead of the SFR into the halo mass using Eqs. (17) and (16) (replacing SFR by LIR). The horizontal line represents the 50% cumulative level. This plot has been made for the galaxies in the range 2.0 < z < 2.2.

5 Effect of cosmology on CIB

5.1 Method

All the previous studies on the CIB were performed assuming a fixed fiducial cosmology. We can see the role of the cosmological parameters for the CIB through Eq. (5) where the CIB power spectrum depends upon the dark matter power spectrum, Hubble’s constant, cosmological energy density parameters (through distance measures), etc. It was thus interesting to study the effect of cosmology on the CIB, i.e. whether changing the cosmology makes a significant impact on the CIB parameters. In order to study this effect, we performed a Fisher matrix analysis over all 12 CIB parameters and 6 cosmological parameters.

For an N-variate multivariate normal distribution X ~ N(μ(θ), Σ(θ)), the (m, n) entry of the Fisher matrix is given as $$I_{m,n}=\frac{\partial \mu }{\partial {\theta }_m}{\Sigma }^{-1}\frac{\partial {\mu }^{\text{T}}}{\partial {\theta }_n}+\frac{1}{2}\text{tr}\left({\Sigma }^{-1}\frac{\partial \Sigma }{\partial {\theta }_m}{\Sigma }^{-1}\frac{\partial \Sigma }{\partial {\theta }_m}\right),$$(18)

where ${()}^{\text{T}}$ denotes the transpose of a matrix and tr() denotes the trace of a square matrix; θ = [θ₁, …, θ_K] is a K-dimensional vector of parameters being considered; μ(θ) = [μ₁(θ), …, μ_N(θ)] are the mean values of the N random variables for which the uncertainties on their measurements are available; Σ (θ) is the covariance matrix for all the variables μ(θ); and $$\frac{\partial \mu }{\partial {\theta }_m}=\left[\frac{\partial {\mu }_1}{\partial {\theta }_m},\dots ,\frac{\partial {\mu }_N}{\partial {\theta }_m}\right]$$(19)

is a vector of the partial derivatives of the N variables for a given parameter θ_m. In our case, Σ(θ) = constant, i.e. the covariance matrix is independent of the parameters and hence the second term in Eq. (18) vanishes as $\frac{\partial \Sigma }{\partial {\theta }_m}=0,$ and therefore $$I_{m,n}=\frac{\partial \mu }{\partial {\theta }_m}{\Sigma }^{-1}\frac{\partial {\mu }^{\text{T}}}{\partial {\theta }_n}.$$(20)

In our case, the parameters being considered are $\theta =[\alpha ,\beta ,\gamma ,\delta ,b_0,b_1,b_2,f_{\nu }^{\text{cal}},H_0,{\Omega }_b{h}^2,{\Omega }_c{h}^2,\tau ,n_s,A_s]$. In the following sections, we explain the construction of the $\frac{\partial \mu }{\partial \theta }$ and Σ matrices.

Fig. 9 Difference between the dark halo mass estimated as the halo mass where 50% cumulative SFR (Mh SFR 50%) is achieved and the halo mass calculated by converting the beff, which is obtained using the LIR (replacing SFR with LIR in Eq. (16)) to the corresponding halo mass using Eq. (17) (Mh LIR weighted). This difference has been calculated and is shown here for all the redshift bins. The dashed line shows the mean of these values: around 0.16 dex. This is the offset applied to the halo mass derived from beff to obtain the mean mass of dark matter halos contributing to CIB.

Fig. 10 Mass of the dark matter halos hosting the galaxies contributing to the CIB as a function of redshift. The black dashed lines show the growth of the dark matter halo mass with redshift. The red dashed line shows the mass that would be obtained directly equating beff and b(M, z) (see Sect. 4.2).

5.2 $\frac{\partial \mu }{\partial \theta }$ matrix construction

Figure 11 shows the different components of the $\frac{\partial \mu }{\partial \theta }$ matrix. To construct this matrix, we have to calculate the partial derivatives of all the variables with respect to all the CIB and cosmological parameters. For every parameter θ_m, the partial derivative is calculated as $$\frac{\partial \mu }{\partial {\theta }_m}=\frac{\mu (\Theta ,{\theta }_m+\delta {\theta }_m)-\mu (\Theta ,{\theta }_m-\delta {\theta }_m)}{2\delta {\theta }_m},$$(21)

where while calculating the partial derivative for a parameter θ_m, all the other parameters represented here by Θ are kept constant at their best fit value. As we are performing the Fisher analysis to see the relative effect of cosmological parameters on CIB, we assume a fiducial cosmology and treat all the cosmological parameters as priors. We maintain the value of δθ small enough to ensure that we are within the Gaussian region for a given parameter.

Component I in the upper left part of the matrix contains the observational constraints mentioned in Sect. 2.3.1, i.e. CIB auto- and cross-power spectra, priors on calibration factors $f_{\nu }^{\text{cal}}$, and CIB–CMB lensing cross-correlation points for all the frequencies. Component II in the same upper left part of the figure contains the external observational constraints mentioned in Sect. 2.3.2, which are the ρ_SFR observations, prior on local bias, and mean level of CIB at different frequencies. It should be noted that all the ρ_SFR values have been calculated by different groups using different cosmologies. In order to perform the Fisher analysis, we needed to make them cosmology independentand therefore we converted them to observed flux values using Eq. (14) assuming a single fiducial Planck 2015 cosmology for all the points. Components III and IV in the lower left part of Fig. 11 calculate the effect of the six cosmological parameters on the CIB variables (observational constraints in component III and external observational constraints in component IV).

Calculations in some of the cases are simplified as certain variables are dependent on only one of the parameters, for example the Gaussian prior on the calibration factor for a given frequency $f_{\nu }^{\text{cal}}$ is independent of all the parameters except itself, and its partial derivative with respect to itself is one. Therefore, the column vector corresponding to this variable in the matrix would simply be {0, 0, …, 1, 0, …, 0}.

Component V in the upper right part of the matrix calculates the effect of CIB parameters on six cosmological parameters. As mentioned above, we treat cosmological parameters as priors, hence they are independent of the CIB parameters, and their partial derivatives with respect to CIB parameters is zero. Therefore, this part of a matrix is an array of the size 12 × 6 with all the elements being zero. Similarly, component VI of the matrix is a diagonal matrix with all the diagonal entries being one and all the off-diagonal elements being zero.

Fig. 11 Details of the matrix $\frac{\partial \mu }{\partial {\theta }_m}$ used in the Fisher matrix analysis.

Fig. 12 Details of the matrix Σ used in the Fisher matrix analysis.

5.3 Σ matrix construction

The Σ matrix is the covariance matrix and it contains the information of the error bars for the variables used in the $\frac{\partial \mu }{\partial \theta }$ matrix. As mentioned in Sect. 2.3.3, we assume Gaussian uncorrelated error bars for the data points and this makes it relatively easy to construct this matrix. Its elements are shown in Fig. 12. The upper left panel of the figure contains the covariance matrix from the CIB variables, which contains the observational and the external observational constraints. As there is no correlation between different CIB variables, the upper left panel is just a diagonal matrix with component I containing the square of the error bars for every variable. The off-diagonal components II and III are both zero. We treat the cosmological parameters as priors. Cosmological parameters are very well constrained by Planck (see Planck Collaboration XIII 2016) and the Monte Carlo sampling chains for their likelihood is provided on the Planck Legacy Archive². We used these chains to calculate the covariance matrix between the six cosmological parameters. This matrix goes in component VI of Fig. 12. Components IV and V are matrices with all the elements being zero.

Once we have $\frac{\partial \mu }{\partial \theta }$ and Σ matrices, it is straightforward to calculate the entries of the Fisher matrix with Eq. (20); in our case, has the dimensions of (18,18).

5.4 Results

We show the results from our Fisher matrix analysis in Table 4. The error bars on the CIB parameters change very little when we add the effect of varying the cosmological parameters. The change in the error bars on the parameters varies from 0.04% to 3.4% which is within the uncertainties on the CIB parameters. Therefore, we conclude that the cosmological parameters determined at the current level of precision have a negligible impact on the CIB model.

6 Conclusion

We have developed a linear CIB model and used conjointly the CIB anisotropies and CIB × CMB lensing cross-correlations measured at large scale by Planck to determine the SFRD up to z = 6 and the evolution of the effective bias. Our paper improved upon the analysis of the linear CIB model performed by Planck Collaboration XXX (2014). We used a functional form for the SFRD (Madau & Dickinson 2014) and a polynomial functional form for the effective bias evolution with redshift. We use effective SEDs derived from the latest observationsand modelling (Béthermin et al. 2015, 2017). The inclusion of the CIB–CMB lensing cross-correlation in our fitting helps us to partially break the degeneracy between the effective bias and the SFRD parameters. In order to get better constraints on the SFRD parameters, we also used external observational constraints on the SFRD at different redshifts from different surveys which are converted to observed flux values to account for the different cosmologies used in the analysis. We also used the mean level of the CIB measured at different frequencies to get a better constraint on the CIB model parameters. Gaussian priors have been put on the local value of the effectivebias along with the calibration factors for different frequencies where the improved values of the uncertainties on the calibration factors from Planck PR2 data release compared to the PR1 data have been taken into account.

With these improved constraints on the CIB model parameters, we derived the SFRD of dusty galaxies up to z = 6. We showed that UV SFRD measurements are consistently lower than our IR SFRD measurement for z < 4 and become compatible (at 1σ) with the IR SFRD at z = 5. The effective bias increases steeply with redshift and is compatible with a number of other measurements, and in particular with the bias obtained from the mass-selected M_⋆ = 10^10.6 M_⊙ sample of Cowley et al. (2018), as expected from the mass distribution of the CIB (Viero et al. 2013). The possible deviation in the redshift bin 1.5 < z < 2 with the SMG sample of Wilkinson et al. (2017) might be a statistical fluctuation arising from the small sample size in this redshift bin. Consistency checks performed on the ρ_SFR and effective bias revealed that the IR SFRD from Bourne et al. (2017) at z ≃ 1 and Koprowski et al. (2017) at z > 3.5 are underestimated. We found that the redshift distribution of the CIB mean was consistently above the published lower limits and that it peaks at z ~ 1.7 at 353 and 217 GHz.

Having measured the effective bias value, we have estimated the typical host dark matter halo mass of galaxies contributing to the CIB. As the effective bias of the galaxies we measure is weighted by their SFR (and even more specifically by their L_IR), the value of the host dark matter halo mass obtained using the effective bias value has to be corrected. We have used the SIDES simulations from Béthermin et al. (2017) to quantify the amplitude of the correction, which is found to be 0.1 dex. Using this offset, we find that the typical mass of the host halos contributing to the CIB is ${\log }_{10}M={12.77}_{-0.125}^{+0.128}$ at z = 2, which is slightly on the high side of the range 10^12.1±0.5 M_⊙ to 10^12.6±0.1 M_⊙ found by Viero et al. (2013) and Planck Collaboration XXX (2014) but in very good agreement with Chen et al. (2016) for faint SMGs of ${\log }_{10}M={12.7}_{-0.2}^{+0.1}$.

Finally, we have quantified for the first time the effect of the cosmology on the CIB parameters. All the previous studies on the CIB had been performed using a fiducial background cosmology. We performed a Fisher matrix analysis to study the effect of changing the cosmology on the CIB and found that it is negligible compared to the existing measurement uncertainties on the CIB parameters.

Acknowledgements

We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA, and CNES, France; from the ANR under the contract ANR-15-CE31-0017; and from the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government programme managed by the ANR. Abhishek Maniyar warmly thanks Sylvain De la Torre and Carlo Schimd for the enlightening discussions on the Fisher matrix analysis.

Appendix A Alternative parametric forms for the effective bias

As mentioned in Sect. 2, we tested alternative parametric forms for the effective bias with different rates of evolution with redshift compared to b_eff used in Eq. (7). The first alternative model we used calculates b_eff following $$b_{\text{eff}}(z)=b_0+b_1{z}^{0.5}+b_{2}z+b_3{z}^{1.5}.$$(A.1)

Fig. A.1 Evolution of the effective bias with redshift derived from the CIB using two different parametrisations, as described in Eq. (A.1) (plotted as the blue line with the 1σ and 2σ regions given in yellow and blue, respectively) and Eq. (7) (plotted as the dark line with the 1σ and 2σ regions given with the red dashed and blue dot-dashed lines, respectively).

The evolution of this effective bias with redshift is shown in Fig. A.1. For comparison, we also show the effective bias from Eq. (7) and its 1σ and 2σ regions. We can see that the constraints on b_eff become slightly weaker for the new effective bias compared to the original one. This effect is especially strong at lower redshifts where the 1σ and 2σ regions on the b_eff are much higher with the new parameterisation.

Similar to this case, we also checked the evolution of b_eff by adding one more parameter to Eq. (A.1), i.e. b_eff(z) = b₀ + b₁z^0.5 + b₂z + b₃z^1.5 + b₄z². We get similar results for this case with b_eff evolving faster than the above case and with higher uncertainties.

References

All Tables

All Figures

	Fig. 1 Posterior confidence ellipses for the CIB model parameters. All the calibration parameters lie within the 1σ range of their prior values. There is no significant degeneracy between any of the calibration parameters and CIB parameters and hence they are not shown here.
In the text

	Fig. 2 Measurements of the CIB auto- and cross-power spectra obtained by Planck and IRAS (extracted from Planck Collaboration XXX 2014) and the best fit CIB linear model. Data points on higher ℓ are not shown as they are not used for fitting. We also show the shot noise (black dash-dotted line) and one-halo term (red dashed line) for all the frequencies.
In the text

	Fig. 3 CIB × CMB lensing cross-power spectra. Our best fit model is shown by the blue curves for different frequency channels. Measurements (red data points) are from Planck Collaboration XVIII (2014). They have been included in the likelihood to calculate and get better constraints on the CIB linear model.
In the text

	Fig. 4 Expected CIB mean level and anisotropy redshift distributions are shown in the lower and upper panels, respectively. Both of the distributions have been normalised for easier comparison, and the different colours show different frequency bands. The CIB anisotropy distribution is shown for ℓ = 300.
In the text

Fig. 5 Expected CIB mean level redshift distribution at 857 GHz compared to observational constraints. The lower limits are from Béthermin et al. (2012b) and Viero et al. (2013). They are shown using black and red triangles, respectively. Measurements from Schmidt et al. (2015) are shown in green. The redshift distribution derived from Koprowski et al. (2017) SFRD constraints is shown with the black dashed line.

In the text

Fig. 6 Evolution of star formation density with redshift as constrained by the linear CIB model. The ± 1σ and ± 2σ confidence region around the median realisation is shown in yellow and blue, respectively. Measurements of obscured star formation density from Gruppioni et al. (2013), Magnelli et al. (2013), and Marchetti et al. (2015), which were used to fit the CIB model, have been added along with the measurements from Bourne et al. (2017) in the upper panel. The Gaussian form of the ρSFR for the IR obscured star formation rate density used by Koprowski et al. (2017) has been plotted as the black dashed line in the upper panel with a corresponding UV part in the lower panel. Unobscured star formation rate density derived from UV from Cucciati et al. (2012), Bouwens et al. (2012), and Reddy & Steidel (2009) are also shown in the lower panel.

In the text

	Fig. 7 Evolution of the effective bias with redshift derived from the CIB compared to the observational values obtained on selected populations of galaxies: Béthermin et al. (2014), Cochrane et al. (2017), Cowley et al. (2018), and Wilkinson et al. (2017) using a sample of BzK-selected galaxies, H α -selected galaxies, mass-selected galaxies, and SMGs, respectively.
In the text

Fig. 8 Normalised cumulative SFR density and LIR density as a function of host dark matter halo mass for galaxies inSIDES simulation are shown in red and green, respectively. We also show the variation of these quantities when they are weighted with the corresponding bias value for the galaxies in blue and black, respectively. The black vertical line represents the mass of the dark matter halo mass calculated converting the beff obtained using the LIR instead of the SFR into the halo mass using Eqs. (17) and (16) (replacing SFR by LIR). The horizontal line represents the 50% cumulative level. This plot has been made for the galaxies in the range 2.0 < z < 2.2.

In the text

Fig. 9 Difference between the dark halo mass estimated as the halo mass where 50% cumulative SFR (Mh SFR 50%) is achieved and the halo mass calculated by converting the beff, which is obtained using the LIR (replacing SFR with LIR in Eq. (16)) to the corresponding halo mass using Eq. (17) (Mh LIR weighted). This difference has been calculated and is shown here for all the redshift bins. The dashed line shows the mean of these values: around 0.16 dex. This is the offset applied to the halo mass derived from beff to obtain the mean mass of dark matter halos contributing to CIB.

In the text

	Fig. 10 Mass of the dark matter halos hosting the galaxies contributing to the CIB as a function of redshift. The black dashed lines show the growth of the dark matter halo mass with redshift. The red dashed line shows the mass that would be obtained directly equating beff and b(M, z) (see Sect. 4.2).
In the text

	Fig. 11 Details of the matrix $\frac{\partial \mu }{\partial {\theta }_m}$ used in the Fisher matrix analysis.
In the text

	Fig. 12 Details of the matrix Σ used in the Fisher matrix analysis.
In the text

	Fig. A.1 Evolution of the effective bias with redshift derived from the CIB using two different parametrisations, as described in Eq. (A.1) (plotted as the blue line with the 1σ and 2σ regions given in yellow and blue, respectively) and Eq. (7) (plotted as the dark line with the 1σ and 2σ regions given with the red dashed and blue dot-dashed lines, respectively).
In the text