Issue 
A&A
Volume 645, January 2021



Article Number  A40  
Number of page(s)  21  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202038790  
Published online  05 January 2021 
Simple halo model formalism for the cosmic infrared background and its correlation with the thermal SunyaevZel’dovich effect
^{1}
Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France
^{2}
Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA
email: abhishek.maniyar@nyu.edu
Received:
29
June
2020
Accepted:
16
October
2020
Modelling the anisotropies in the cosmic infrared background (CIB) on all the scales is a challenging task because the nature of the galaxy evolution is complex and too many parameters are therefore often required to fit the observational data. We present a new halo model for the anisotropies of the CIB using only four parameters. Our model connects the mass accretion on the dark matter haloes to the star formation rate. Despite its relative simplicity, it is able to fit both the Planck and Herschel CIB power spectra and is consistent with the external constraints for the obscured star formation history derived from infrared deep surveys used as priors for the fit. Using this model, we find that the halo mass with the maximum efficiency for converting the accreted baryons into stars is log_{10}M_{max} = 12.94_{0.02}^{+0.02} M_{⊙}, consistent with other studies. Accounting for the mass loss through stellar evolution, we find for an intermediateage galaxy that the star formation efficiency defined as M_{⋆}(z)/M_{b}(z) is equal to 0.19 and 0.21 at redshift 0.1 and 2, respectively, which agrees well with the values obtained by previous studies. A CIB model is used for the first time to simultaneously fit Planck and Herschel CIB power spectra. The high angular resolution of Herschel allows us to reach very small scales, making it possible to constrain the shot noise and the onehalo term separately, which is difficult to do using the Planck data alone. However, we find that large angular scale Planck and Herschel data are not fully compatible with the smallscale Herschel data (for ℓ > 3000). The CIB is expected to be correlated with the thermal SunyaevZel’dovich (tSZ) signal of galaxy clusters. Using this halo model for the CIB and a halo model for the tSZ with a single parameter, we also provide a consistent framework for calculating the CIB × tSZ cross correlation, which requires no additional parameter. To a certain extent, the CIB at high frequencies traces galaxies at low redshifts that reside in the clusters contributing to the tSZ, giving rise to the onehalo term of this correlation, while the twohalo term comes from the overlap in the redshift distribution of the tSZ clusters and CIB galaxies. The CIB × tSZ correlation is thus found to be higher when inferred with a combination of two widely spaced frequency channels (e.g. 143 × 857 GHz). We also find that even at ℓ ∼ 2000, the twohalo term of this correlation is still comparable to the onehalo term and has to be accounted for in the total crosscorrelation. The CIB, tSZ, and CIB × tSZ act as foregrounds when the kinematic SZ (kSZ) power spectrum is measured from the cosmic microwave background power spectrum and need to be removed. Because of its simplistic nature and the low number of parameters, the halo model formalism presented here for these foregrounds is quite useful for such an analysis to measure the kSZ power spectrum accurately.
Key words: infrared: diffuse background / cosmic background radiation / submillimeter: galaxies / galaxies: clusters: general / cosmology: observations / methods: data analysis
© A. Maniyar et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The cosmic infrared background (CIB) is made up of the cumulative emission of the infrared radiation from the dusty starforming galaxies throughout the Universe. It traces the star formation history of the Universe, which spans a wide range of redshift 0 ≤ z ∼ 6. Measurements of the CIB can thus be used as a powerful tool to map the star formation at high redshifts (Knox et al. 2001). Although the CIB was first detected by Puget et al. (1996), Lagache & Puget (2000) and Matsuhara et al. (2000) were the first to detect and discuss the anisotropies in the CIB that are due to unresolved extragalactic sources. The CIB includes correlated anisotropies that are excellent probes of the largescale structure of the Universe (e.g. Hanson et al. 2013). These were first discovered by Spitzer (Lagache et al. 2007) and were then subsequently accurately measured by Planck and Herschel (Planck Collaboration XXX 2014; Viero et al. 2013).
Another such tracer of the underlying dark matter distribution are massive galaxy clusters. Hot electrons in these galaxy clusters Compton scatter the CMB photons and give rise to the socalled thermal SunyaevZel’dovich (tSZ) effect. A part of the CIB originates in the dusty starforming galaxies residing in the galaxy clusters. Thus, the tSZ and the CIB are expected to be correlated to a certain extent. This correlation has indeed been indirectly measured and shown to be positive bypg
Reichardt et al. (2012, 2020), Dunkley et al. (2013), George et al. (2015), and Choi et al. (2020) and thus has to be considered in the CMB power spectrum data analysis. Planck Collaboration XXIII (2016) also reported a measurement of the crosscorrelation between the tSZ and CIB using Planck data. In order to model the CIB × tSZ signal, we in turn need accurate models of the CIB anisotropies and the tSZ.
On large angular scales, we can use the fact that the clustering of the CIB traces the largescale distribution of matter in the Universe up to some bias factor. This makes modelling the CIB anisotropies on large angular scales quite straightforward, as reported in for example Planck Collaboration XXX (2014) and Maniyar et al. (2018). However, in order to describe the anisotropies on both the large and small angular scales coherently, a “halo model” approach as developed by Cooray & Sheth (2002) is generally used. With the assumption inside the halo model that all the galaxies reside in the dark matter haloes, the clustering can be considered as the sum of two components: onehalo term (P_{1h}), which takes the smallscale clustering due to the correlations between the galaxies within the same halo into account; and twohalo term (P_{2h}), which accounts for the clustering on large scales due to the correlations between galaxies in different haloes. Along with the assumption that all the dark matter lies within the collapsed and symmetric haloes, four ingredients are required to characterise the galaxy power spectrum within the halo modelling context: the number density of the haloes per unit mass given by the halo mass function; the halo bias between the haloes and the dark matter; the spatial distribution of the dark matter inside a halo given by the halo density profile; and the halo occupation distribution (HOD), which is a prescription for filling the dark matter haloes with galaxies.
The first generation of the models built to interpret the CIB anisotropies was based either on a HOD model or a combination of models of emissivities of the infrared galaxies and a linear bias (Knox et al. 2001; Lagache et al. 2007; Amblard & Cooray 2007; Viero et al. 2009; Planck Collaboration XVIII 2011; Amblard et al. 2011; Xia et al. 2012). These approaches assumed that the emissivity density is traced by the galaxy number density, implying that all galaxies contribute equally to the emissivity, regardless of their host halo masses. This would mean that all the galaxies have the same luminosity. However, as has been pointed out by Shang et al. (2012), both the luminosity and clustering of the galaxies are closely related to the host halo mass. In general, galaxies situated in more massive haloes are more luminous as a result of a higher stellar mass, and they are also more clustered. When this effect is neglected, the clustering signal on smaller angular scales might be interpreted as being due to a very high number of satellite haloes (which was the case for Amblard et al. 2011) compared to what is found in numerical simulations (discussion in Viero et al. 2013).
Subsequently, several studies such as Shang et al. (2012), Viero et al. (2013), and Planck Collaboration XXX (2014) have improved upon the previous halo models by considering a link between the galaxy luminosity (L) and the host halo mass (M_{h}) in their model (through a L − M_{h} relation). Although their approach is able to fit the CIB power spectra, their description of the infrared galaxies is quite simple (e.g. a single spectral energy distribution, SED, for all galaxy types, but with evolving dust temperature or without scatter on the L − M_{h} relation). These models are useful to derive quantities such as the halo mass for the most efficient star formation, but it is hard to test their validity because a good fit can be obtained easily because of the high number of free parameters in these models. Without considering any priors, the predictions of these models for physical quantities such as the star formation rate density (SFRD) moreover do not match the corresponding constraints from the linear model or galaxy surveys (e.g. Planck Collaboration XXX 2014). This shows the need for physically motivated models that in addition to the power spectra can provide a good fit or prediction for other physical quantities such as the SFRD. Béthermin et al. (2013) used a semiempirical model based on the observed relation between the stellar mass M_{*} and the SFR, which they linked to the corresponding halo mass using abundance matching. This model gives CIB power spectra that are consistent with the measurements. Inspired by their findings of the SFR/BAR relation with respect to the halo mass (where BAR represents the baryonic accretion rate), we develop a simpler halo model for the CIB anisotropies with just four parameters. Our model connects the mass accretion onto the dark matter haloes to the corresponding SFR.
The tSZ is measured through the socalled Compton parameter (y; Sect. 4). The Planck Collaboration provided an allsky map of this Compton y parameter and an estimate of the tSZ angular power spectrum up to ℓ ≈ 1300 (Planck Collaboration XXI 2014; Planck Collaboration XXII 2016). Bolliet et al. (2018) used these data to constrain the cosmological parameters (equation of state of the dark energy w in particular) along with the tSZ parameter. We use this halo model of the tSZ to calculate the tSZ power spectra.
Addison et al. (2012) calculated the CIB−tSZ correlation within the halo model framework. They first presented a formalism to calculate this correlation using a CIB halo model that does not account for the dependence of the source flux on the halo mass, and then expanded their formalism to account for this effect. However, the CIB halo model they finally considered (from Xia et al. 2012) to calculate the CIB−tSZ correlation does not account for the aforementioned L − M_{h} dependency. We present a halo model formalism to calculate the CIB × tSZ crosscorrelation using our new halo model for the CIB and the halo model for the tSZ from Bolliet et al. (2018).
This paper is structured as follows. We begin in Sect. 2 by presenting our halo model for the CIB anisotropies and subsequently present the CIB × CMB lensing correlation within this framework. Section 3 then presents the constraints on the CIB model parameters through the data and corresponding results. In Sect. 4 we provide the halo model for the tSZ. Finally, in Sect. 5 we present the halo model formalism for the CIB × tSZ correlation and the predictions for its power spectra and angular scale dependence^{1}. Conclusions are given in Sect. 6. Appendices A–C provide details of the CIB halo model formalism and the comparison between the Planck and Herschel CIB power spectrum data.
When required, 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_{0} = 67.47 km s^{−1} Mpc^{−1}.
2. New halo model for the CIB power spectrum
The starting point for our model is the accretion of matter onto the dark matter haloes. We then connect the accretion of the baryonic gas onto the dark matter haloes to the SFR corresponding to these haloes. The SFR is defined separately for the central and satellite galaxies in a halo. Using the SFR from a given halo, we then calculate the emissivity of all the haloes of mass M_{h} at a given redshift, which then is used to calculate the angular power spectrum of the CIB anisotropies.
The angular power spectrum of the CIB anisotropies is defined as
where ν is the frequency of the observation and I^{ν} is the intensity of the CIB measured at that frequency. The intensity is a function of the comoving emissivity j through
where χ(z) is the comoving distance to redshift z, and a = 1/(1 + z) is the scale factor of the Universe, and δj(ν, z) are the emissivity fluctuations of the CIB. We expand the Eq. (2) this way between the mean value and its fluctuations because Eq. (1) shows that the power spectrum is calculated using the fluctuations around the mean value. Combining Eqs. (1) and (2), and using the Limber approximation (Limber 1954) for the flat sky, which helps us avoid the spherical Bessel function calculations and makes the computation easier, we therefore obtain
where at a given redshift, is the 3D power spectrum of the emissivity and is defined as
Thus we have to calculate δj to obtain the CIB angular power spectrum. For this purpose we connect (Sect. 2.3) the SFR of the haloes with the specific emissivity and integrate it over all the halo masses and redshift range to obtain the CIB power spectra.
2.1. From accretion onto the dark matter haloes to SFR
The dark matter haloes grow in mass over time through diffuse accretion and mergers with other lower mass haloes (e.g. Fakhouri et al. 2010). Accretion and merger processes are also responsible for the growth of stellar mass in galaxies through galaxygalaxy mergers and through accretion of the gas. The baryonic gas accreted by a given dark matter halo would form stars depending upon certain factors. This is the starting point of our model. As we mentioned earlier, previous studies used a parametric L − M_{h} relation to derive the power spectra. Instead of assuming an L − M_{h} relation with an evolution in redshift, we connect the accretion rate onto a dark matter halo described above with the corresponding SFR. This gives us an SFR−M_{h} relation that in substance is similar to an L − M_{h} relation. The difference between our approach and that of others is a more physical starting point of the parametrisation.
This approach is physically motivated. This link between the accreted baryons and SFR is quite natural. The stars form out of the gas reservoirs within their host galaxies that reside in the dark matter haloes. The amount of gas at a given time depends upon the amount of gas accreted by the host dark matter halo. We assumed that this accreted gas is converted into stars with an efficiency that is a function of the mass of the halo and redshift. We used a lognormal parametrisation between the halo mass and the ratio of the SFR and the baryonic accretion rate BAR for a halo (i.e. SFR/BAR)
where M_{h} is the halo mass, M_{max} represents the mass with the highest star formation efficiency η_{max}, and σ_{Mh}(z) is the variance of the lognormal, which here represents the range of masses over which the star formation is efficient. SFR/BAR represents the efficiency (η) of the dark matter halo of a given mass (M_{h}) at a redshift (z) to convert the accreted baryonic mass into stars.
The choice of the lognormal shape is quite logical. Several studies (e.g. Viero et al. 2013; Planck Collaboration XXX 2014; Maniyar et al. 2018) have found that the dark matter haloes within the mass range 10^{12} − 10^{13} M_{⊙} form the stars most efficiently. With an empirical model, Béthermin et al. (2013) showed in their Fig. 17 that the star formation efficiency as a function of instantaneous halo mass is highest in haloes with masses ∼10^{12} M_{⊙}. This mass does not change considerably over a range of redshifts, whereas the efficiency falls off drastically for masses above and below the most efficient mass. This effect can be understood physically. Below the mass for which star formation is most efficient, the gravitational potential of the dark matter halo is lower and the supernovae feedback is strong enough to remove the gas from the galaxy (e.g. Silk 2003) and thereby decreases the star formation. On the other side of the spectrum, at higher masses, the cooling time of the gas becomes much longer than the freefall time (e.g. Kereš et al. 2005). This suppression of the isotropic cooling of the gas could be due to the energy injection in the halo atmosphere by active galactic nuclei (AGN), which in turn suppresses the star formation (e.g. Somerville et al. 2008). Therefore we assumed a lognormal shape for the SFR efficiency whereby the star formation is highest for halo mass M_{max}, and a significant contribution to the SFR comes from a range of masses around M_{max} driven by σ_{Mh}, and the SFR falls off considerably on very high and very low masses.
For a given halo mass at a given redshift, the BAR is given as
where Ω_{b} and Ω_{m} are the dimensionless cosmological baryonic density and total matter density, respectively (thus, the ratio of the two gives the baryonic matter fraction at a given redshift, which is in fact constant with redshift because they have the same evolution with redshift). Ṁ(M_{h},z) is the mass growth rate. Fakhouri et al. (2010) provided the mean and median mass growth rates of haloes of mass M_{h} at redshift z. We used the mean mass growth rate, given as
This approach assumes that there is no “gas reservoir effect”, that is, the accreted gas is immediately converted into stars and is not collected over time to form reservoirs in Daddi et al. (2010) and around galaxies (Cantalupo et al. 2012). It has been shown (e.g. Saintonge et al. 2013; Béthermin et al. 2015; DessaugesZavadsky et al. 2015) that the depletion timescale, which is the ratio of the mass of the molecular gas to the SFR, ranges from ∼0.5−1 Gyr, which is much shorter than the typical galaxy evolution timescales (several billion years). Thus, the gas reservoir effect is not expected to affect the results over the long timescale that we considered here. Along the same lines, it is also assumed that no recycled gas (i.e. the gas expelled by the supernovae) contributes to the star formation (semianalytical models of e.g. Cousin et al. 2015, 2019 showed that feedback from supernovae can play a part in regulating star formation). In spite of these assumptions, we show that this simple physical model describes the CIB power spectra well.
In our model, we consider the maximum efficiency η_{max} and mass of maximum efficiency M_{max} to be independent of redshift, that is, they do not evolve with redshift. However, we let σ_{Mh} evolve with redshift. The motivation for letting σ_{Mh} evolve with redshift is that we wished to accommodate the star formation from massive haloes at higher redshifts and reduce it at lower redshifts. It has been observed (Popesso et al. 2015) that at lower redshifts (z ≤ 1.5−2), star formation is quite inefficient in massive haloes (typical galaxy cluster environments), that is, at low redshifts, massive haloes contain mostly passive galaxies. In contrast, it has been shown that at highredshift massive galaxies (often residing in the protoclusters, i.e., the progenitors of the clusters at redshift zero) can have efficient star formation (e.g. Miller et al. 2018; Wang et al. 2018). Because the lognormal parametrisation leaves a tail on the high mass end, it might mimic this effect, and the choice of the this shape is therefore justified.
Thus we let σ_{Mh} evolve with redshift as
where z_{c} is a redshift below which σ_{Mh} evolves with redshift, σ_{Mh0} is the value of σ_{Mh} above z_{c}, and τ is the parameter driving this evolution with redshift (τ here should not be confused with the optical depth parameter from the CMB analysis). Following the reasoning mentioned before, this evolution was applied only for haloes with masses greater than the mass of maximum efficiency M_{max} and below redshift z_{c}, that is, the parametrisation is not a symmetrical lognormal below redshift z_{c}. The width of the lognormal is smaller at the side of the curve with haloes higher in mass than M_{max} below redshift z_{c}. However, above redshift z_{c}, the parametrisation is a symmetrical lognormal with no evolution in the width of the lognormal σ_{Mh0}. We fixed z_{c} = 1.5. Other values for z_{c} were tried and gave approximately the same results, but the model with z_{c} = 1.5 provided the best fit for the SFRD history.
2.2. SFR for the haloes and subhaloes
For a given value of the halo mass and redshift, we can calculate η using Eq. (5) and multiply it by the corresponding BAR calculated using Eq. (6) to obtain the SFR, that is,
This is the procedure with which the SFR can be obtained for the haloes. To calculate the SFR for the subhaloes residing within these haloes, the procedure is slightly modified. We first assumed that for a given halo with mass M_{h}, the subhalo masses (m_{sub}) range from M_{min} to M_{h}. In this analysis, we fixed M_{min} = 10^{5} M_{⊙}. A change of the minimum mass between 10^{4} M_{⊙} and 10^{8} M_{⊙} changes the calculation of the power spectra only negligibly. The SFR for the subhaloes can be estimated in two ways. The first way is an approach similar to the one for the haloes, which is calculating the efficiency η and then multiplying with the BAR value to obtain the SFR, that is, replacing M_{h} by m_{sub} in Eq. (9). This assumes the same lognormal parametrisation of η for subhaloes as of the central haloes. The other way to estimate the SFR in subhaloes is
that is, the SFR for the subhalo is obtained by weighing the halo SFR (SFR_{c}) by the ratio of subhalo to halo mass. For every subhalo of a given halo, we estimated the SFR with both these approaches and took the smaller of the two as representative of the SFR for the subhalo.
The reasoning for this is explained in Fig. 1. We first consider case 1 in the figure. In this case, the main halo has a mass (∼10^{12.9} M_{⊙}) very near to the efficiency peak of star formation (for this particular choice of parameters for the lognormal), that is, the central galaxy forms stars very efficiently (see Eq. (9)). This halo has subhaloes ranging from 10^{5} M_{⊙} − 10^{12.9} M_{⊙}. We take the case of a subhalo with mass <10^{11} M_{⊙}. As was pointed out before, subhaloes with very low mass have a low gravitational potential, and it is hard for them to hold on to the gas inside against the pressure from supernova feedback, for instance. Thus, they are expected to have low star formation. This is satisfied in case 1 because at lower masses, the efficiency is indeed very low and is not expected contribute significantly to the total SFR of the halo. In this case, the SFR for the subhaloes can therefore be directly estimated by substituting the halo mass (M_{h}) by subhalo mass (m_{sub}) in Eq. (9). As an example, for a subhalo of mass 10^{11} M_{⊙} belonging to a central halo of mass 10^{12.9} M_{⊙}, the SFR calculated using Eq. (9) is 2% of the value that we derive for the SFR obtained using Eq. (10), and we therefore take the former as the SFR value.
Fig. 1. Lognormal parametrisation (Eq. (5)) between the halo mass (M_{⊙}) and ratio between the SFR and the baryonic accretion rate, η. We show two extreme cases: haloes near the efficiency peak contain subhaloes with very low mass (case 1), and very massive haloes that contain subhaloes near the efficiency peak (case 2). If the same recipe were used to calculate the satellite galaxy SFR in these two cases, unphysical values might result within the assumptions of our model, and therefore we suggest two different ways to calculate the SFR for satellite galaxies (Sect. 2.2). 
In the second case, the main halo is quite massive (>10^{14} M_{⊙}) and far away from the efficiency peak. It therefore does not have strong star formation. However, in this case, this halo can contain a subhalo with a mass of about the efficiency peak (∼10^{12.9} M_{⊙}). According to Eq. (9), when we substitute M_{h} by m_{sub}, this subhalo will have a significant amount of star formation. Again, as pointed out earlier, the gas inside massive haloes is quite hot, and there are mechanisms at play (e.g. Xray heating and AGN feedback) that suppress the gas cooling and hence make it difficult to form stars. Moreover, the central massive galaxies in these haloes can strip the gas out of the satellite galaxies in the subhaloes and thus decrease star formation. These subhaloes are therefore not expected to contribute significantly to the SFR, in contrast to what would be obtained with Eq. (9). This inherently assumes instantaneous quenching, that is, the satellite galaxies are quenched at the same time as the central galaxies in a given parent halo. If we were to explicitly avoid this assumption, we would need to introduce an additional quenching parameter as a function of subhalo mass and redshift. This would result in one or two additional parameters for the model. However, as we mentioned earlier, our main purpose here is to build a very simple halo model of the CIB with as few parameters as possible. One way to correct for this in these cases therefore is to weight the SFR of the main halo by the mass fraction of the corresponding subhalo, that is, use Eq. (10) to obtain the SFR of the subhalo. The SFR in this case would be lower than the rate obtained using Eq. (9) (substituting M_{h} by m_{sub}, of course). Again as an example, for a subhalo of mass 10^{12.9} M_{⊙} belonging to a central halo of mass 10^{14} M_{⊙}, the SFR calculated using Eq. (9) is twice higher than the value we obtain for the SFR with Eq. (10) at z = 2, and we therefore take the latter as the SFR value.
Although Eqs. (9) and (10) seem useful to estimate the SFR for subhaloes when we have cases similar to cases 1 and 2, for every halo we therefore estimate the SFR for the corresponding subhaloes using both these methods at every redshift and select the SFR with the lower value. This automatically takes care of the extreme cases and helps us avoid adding more parameters to the model.
2.3. SFR to CIB power spectra
The onehalo term for the CIB power spectrum takes the clustering of the galaxies within a halo of mass M_{h} into account and was calculated following Béthermin et al. (2013) (where k = ℓ/χ),
where is the halomass function, u(k, M_{h}, z) is the Fourier transform of the density profile describing the density distribution inside the halo (here we consider the density distribution to be a NavarroFrenkWhite (NFW) profile Navarro et al. 1997), and is the specific emissivity of the central and satellite subhaloes at a given frequency and redshift for a given halo mass as defined in Béthermin et al. (2013). After we calculate the specific emissivity term for the central and satellite terms, it is therefore straightforward to calculate the onehalo power spectrum. For simplicity, we omitted the M_{h} and z dependence from and terms from all the equations.
For the central galaxies, the differential emissivity is calculated as
where is the effective SED of the infrared galaxies at a given redshift for a given frequency. SFR_{c} is the SFR for the central galaxies with a given halo mass (Eq. (9)). K is the Kennicutt constant (K = SFR/L_{IR}), which has a value of 1 × 10^{−10} M_{⊙} yr^{−1} for a Chabrier IMF, and L_{IR} is the infrared luminosity (8−1000 μm).
For the satellite galaxies in the subhaloes (Béthermin et al. 2013),
where is the subhalo mass function for the satellite galaxies with a subhalo mass m_{sub}. The effective SEDs for the satellite galaxies are assumed to be the same as those of the central galaxies. The SFR_{sub} is calculated using Eqs. (9) and (10), and the smaller of the two values is taken as the SFR value for those galaxies.
In our analysis, we assumed the halo mass function from Tinker et al. (2008) and the subhalo mass function from Tinker & Wetzel (2010). are the same as we used for the linear clustering model of the CIB anisotropies from Maniyar et al. (2018). They were computed using the method presented in Béthermin et al. (2013), but assuming the new updated SEDs calibrated with Herschel data presented in Béthermin et al. (2015, 2017). A stacking analysis was used to measure the evolution of the average midinfrared to milimeter emission of the massive starforming galaxies up to z = 4. With this technique, we found that for the mainsequence galaxies we used in the analysis, the mean intensity of the radiation field, which is strongly correlated with the dust temperature, rises with redshift. Thus the dust in these new SED templates is warmer at z > 2 than in the previous templates used in Béthermin et al. (2013). We prefer these templates over the other templates (e.g. from Gispert et al. 2000 using FIRAS measurements) because they reproduce all recent measurements of galaxy counts from the midIR to the radio wavelength range, including counts per redshift slice.
The Fourier transform of the NFW profile is given as (e.g. van den Bosch et al. 2013)
where it is to be noted that c is not the speed of light; it is the socalled halo concentration parameter defined as c(M_{h}, z) ≡ r_{200}(M_{h}, z)/r_{⋆}(M_{h}, z), with r_{⋆}(M_{h}, z) being a characteristic radius (see e.g. Navarro et al. 1997) and r_{200} being the radius containing the mass giving 200 times the critical density of the Universe at a redshift (some studies use the average density instead of the critical density; we used the critical density), Ci(x) and Si(x) are standard cosine and sine integrals, respectively, μ ≡ kr_{⋆}, where δ_{200} is a dimensionless amplitude, which can be expressed in terms of the halo concentration parameter as
The twohalo term for the power spectrum of the CIB takes the clustering between galaxies in two different haloes of mass M_{h} and into account and is calculated as (Béthermin et al. 2013) (where k = ℓ/χ)
where b(M_{h}, z) is the halo bias prescription given by Tinker et al. (2010), P_{lin}(k, z) is the linear matter power spectrum, which we calculated using CAMB^{2}. We used the u(k, M_{h}, z) term only for the subhaloes. This means that in our analysis, the central galaxy is assumed to be at the centre of the halo and the satellite galaxies in the subhaloes are distributed according to the NFW profile. The clustering term given by provides the crosspower spectrum between two different haloes (M_{h} and ) under the assumption (Cooray & Sheth 2002) . We can simplify Eq. (16) by defining
which is the emissivity of the halo weighted by the corresponding bias. Therefore Eq. (16) becomes (where k = ℓ/χ)
This way of calculating the twohalo term significantly reduces the time because it reduces the number of integrals. It is also useful when the CIB−CMB lensing crosscorrelation is calclutated (see Sect. 2.5). Finally, we have
and
Here and describe the clustered and total CIB anisotropy power spectrum, respectively. is the shotnoise component, which we cannot predict using our halo model formalism here. The shotnoise component is scale independent and thus has a constant flat power spectrum for a given frequency combination. We thus fit for this constant directly at every frequency channel as described in Sect. 3.1.
2.4. ρ_{SFR}
When the formalism for calculating the CIB angular power spectrum in the halo modelling context is defined, it is desirable that the halo model is able to fit or reproduce statistical properties of dusty starforming galaxies, such as the infrared (IR) SFRD of the Universe. In the context of our halo model, the SFRD for the central galaxies is calculated as
and for the satellite galaxies in the subhaloes, it reads
where SFR(m_{sub}, z) is the satellite galaxy SFR calculated as explained in Sect. 2.2.
Thus, the total SFRD is calculated by adding Eqs. (21) and (22) as
2.5. CIB–CMB lensing crosscorrelation
The largescale distribution of the matter in the Universe gravitationally deflects the CMB photons that freely propagate toward us from the last scattering surface. This phenomenon is called gravitational lensing, and it leaves imprints on the temperature and polarisation anisotropies of the CMB. Because the CIB is an excellent tracer of the largescale structure of the Universe (Maniyar et al. 2019), its anisotropies are expected to be strongly correlated with the CMB lensing, which has indeed been observed and measured (Planck Collaboration XVIII 2014). Within the context of our halo model, we can calculate this CIB−CMB lensing crosscorrelation as (Béthermin et al. 2013) (where k = ℓ/χ)
where Φ(ℓ, z) is given by (Challinor & Lewis 2005)
In this equation, χ_{⋆} is the comoving distance to the last scattering surface, Ω_{m} is the matter density parameter, H_{0} is the value of the Hubble parameter today, and a is the scale factor of the Universe. The measurement of this crosscorrelation is not used in our likelihood to fit the CIB model parameters, but we show that our bestfit model accurately predicts this crosscorrelation.
3. Constraints on the model through data
3.1. Observational constraints on the power spectra
We used the CIB angular power spectra as measured by Planck Collaboration XXX (2014). The measurements were obtained by cleaning the CMB and Galactic dust from the Planck frequency maps. They were further corrected for the 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). The CIB × CIB power spectra were corrected for absolute calibration difference between the PR1 and PR2 data release of Planck. Absolute calibration uncertainties for PR2 are equal to 6.1% and 6.4% at 545 and 857 GHz (5% of which comes from planet models), respectively (Planck Collaboration VIII 2016). However, it has been shown that the planet calibration agrees with CMB calibrations within 1.5% for 545 GHz (Planck Collaboration Int. XLVI 2016).
We also used the CIB power spectra measurement from Herschel/SPIRE data at 600, 857, and 1200 GHz provided by Viero et al. (2013). Bertincourt et al. (2016) calculated the crosscalibration factors between the power spectra measured by Viero et al. (2013) at 600 and 857 GHz and the power spectra measured by Planck Collaboration XXX (2014) at 545 and 857 GHz. They are 1.047 ± 0.0069 and 1.003 ± 0.0080, respectively (for the Planck PR2 data release). Absolute calibration uncertainties for these SPIRE data are 9.5%.
Planck and Herschel measurements are given with the νI_{ν} = constant photometric convention. As a result, the power spectra computed by the model need to be colourcorrected from the CIB SEDs to this convention. These colour corrections are 1.119, 1.097, 1.068, and 0.995 at 217, 353, 545, and 857 GHz, respectively, for Planck (Planck Collaboration XXX 2014). They are equal to 0.974, 0.989, and 0.988 for 600, 857, and 1200 GHz channels of Herschel/SPIRE, respectively, using the extended relative spectral response function (Lagache et al. 2020). The CIB power spectra were then corrected as
The CIB power spectra error bars do not account for the absolute calibration uncertainties. To account for these uncertainties, Béthermin et al. (2011) introduced a calibration factor for galaxy number counts. Using a similar approach here, we set a Gaussian prior on these calibration factors for every frequency channel. The Herschel power spectra were crosscalibrated with respect to the Planck power spectra, therefore we fixed the calibration factors at the Planck frequencies at 1 and set Gaussian priors on calibration factors for 600 and 857, and 1200 GHz channels from Herschel/SPIRE with an initial value of 1.0470, 1.0030, and 1.0000 with 1σ error bars at 0.0069, 0.0080, and 0.0500, respectively (Bertincourt et al. 2016). In Table C.1 we show our bestfit values when we fit our model to the Planck data alone and allowed the individual calibration parameters for all the Planck channels to vary with a 5% prior on them. The posterior values show that the error bars on the calibration parameters for Planck (which are very close to 1) are indeed very small. For our purposes here, the assumption of taking a perfect calibration for Planck therefore does not strongly affect the results.
When only Planck/HFI data are used, it is hard to differentiate between the onehalo and the shotnoise term, and therefore, Planck Collaboration XXX (2014) placed strong flat priors on the shot noise at all the frequencies. They took values of the shot noise based on the model of Béthermin et al. (2012), and was the width of the prior given by their 1σ error bars. However, using Herschel/SPIRE data, the CIB power spectra were measured to very small scales (ℓ_{max} ∼ 30 360). At these highest multipoles, the power spectra are dominated by the shot noise, and it is therefore better possible to constrain the shot noise together with the onehalo term with these data than using the Planck/HFI data alone. Because Herschel/SPIRE data can distinguish between the onehalo and the shotnoise term, this paves the way to include shot noise as a parameter for every pair of the power spectra with broad priors, that is, 10 shotnoise parameters for HFI, as done by Planck Collaboration XXX (2014), and 6 for SPIRE data. This would be 16 additional parameters. In order to circumvent this problem, we calculated the correlation of the 10 shotnoise parameters for HFI and 6 for SPIRE using the values predicted by the model of Béthermin et al. (2012). We then fit for the shotnoise parameters for only the autopower spectra and used the correlation matrix to obtain the shotnoise level for every frequency pair. Thus, in the end we have only 1 shot noise parameter per frequency, that is, 4 for HFI and 3 for SPIRE. We placed broad flat priors on them with sufficient width ([0.1−2] times the value estimated by the Béthermin et al. 2012 model on either side) to avoid biasing the parameter estimation. We also tested the model by placing the shot noise at all the SPIRE auto and crosspower spectra as parameters instead of only taking the autopower as the parameters because SPIRE might be able to differentiate between the clustering and the shot noise because its angular resolution is very high. We find similar values for the shot noise in both cases and therefore selected the case with fewer parameters.
3.2. External observational constraints
In addition to the CIB auto and crosspower spectra, we used external constraints from the mean CIB intensity values at different frequencies, and the SFRD measurements at different redshifts:

1.
As mentioned before, the previously used halo models overpredicted the SFRD of the Universe compared to the SFRD measured by external groups using the infrared luminosity functions. We would like the halo model to be able to fit the CIB power spectra and also produce the correct SFRD history. Thus, we used the ρ_{SFR} measurements at different redshifts that were obtained by measuring the infrared luminosity functions from Gruppioni et al. (2013), Magnelli et al. (2013), Marchetti et al. (2016) as priors while performing the fit. In order to account for the different sets of cosmological parameters used by them, we converted the SFRD values into the observed flux between 8 and 1000 μm per redshift bin per solid angle, as done in Maniyar et al. (2018).
In Appendix A we show the SFRD history produced by our halo model when this external prior is not considered in the fit. There is indeed a significant change in the SFRD history predicted by the model when we do not include any priors at all. This means that including this prior is quite important to obtain physical results from our model.

2.
The mean level of the CIB was deduced at different frequencies using the galaxy number counts. We used these measurements as constraints on our model. Similar to the case of the power spectra, the mean level of the CIB computed by the model needs to be colourcorrected. The values of the mean level with their corresponding frequencies and the colour corrections we used are given in Table 2 of Maniyar et al. (2018).
3.3. Fitting the data
We performed a Markov chain Monte Carlo (MCMC) analysis in the global CIB parameter space using the Python package “emcee” (ForemanMackey et al. 2013). We have a 14dimensional parameter space:

physical model parameters {η_{max},M_{max},σ_{Mh0},τ},

calibration factors ,

shot noises .
The global χ^{2} has a contribution from the CIBxCIB, priors on calibration factors, and priors imposed by the external observational constraints mentioned above. We assumed Gaussian uncorrelated error bars for measurement uncertainties, which simplifies the covariance matrix that ultimately contains only diagonal terms.
3.4. Results
Figures 2 and 3 show the best fit for the halo model to the observational data points for all CIBxCIB auto and crosspower spectra from Planck/HFI and Herschel/SPIRE, respectively. For most of the cases, the onehalo term is smaller than the shot noise and can potentially be ignored on smaller scales. Based on the bestfit parameters and Eq. (24), we calculated the CIBxCMB lensing power spectra. Figure 4 shows these CIBxCMB lensing power spectra. The crosscorrelation 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, and as discussed in Planck Collaboration XVIII (2014), the nonlinear term can be neglected in this range of multipoles. As was done for the CIB, these power spectra were corrected for absolute calibration difference between the PR1 and PR2 data release of Planck. The bestfit model agrees very well with the data points.
Fig. 2. Measurements of the CIB auto and crosspower spectra obtained by Planck/HFI (Planck Collaboration XXX 2014) and the bestfit CIB halo model with its different components. 
Fig. 3. Measurements of the CIB auto and crosspower spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the bestfit CIB halo model with its different components. 
Fig. 4. CIB × CMB lensing crosspower spectra measured by Planck Collaboration XVIII (2014) in orange and our bestfit model in blue curves. 
We present the results from our fit for our model parameters in Table 1. The shotnoise values provided here have to be multiplied with the corresponding colour corrections at each frequency to obtain them in the νI_{ν} = constant convention. The posterior of all the parameters with a Gaussian prior (calibration factors) are within a 1σ range of the prior values. We derive a χ^{2} of 113 for 80 data points for Planck power spectra. Similarly, we derive a χ^{2} of 247 for Herschel power spectra for 102 data points. When we only fit for the Planck data without the Herschel power spectrum data but with the external priors, we derive a χ^{2} of 85 for Planck for 80 data points. Similarly, fitting only for the Herschel power spectra with the external priors results in the χ^{2} of 221 for 102 data points. This shows that when we individually fit for either the Planck or Herschel data, we derive a better χ^{2} value than when we fit for both of them together.
Marginalised values of all the model parameters given at a 68% confidence level.
We also performed this analysis with the updated values of the mean level of the CIB at 350, 350, and 500 μm Herschel/SPIRE bands used as priors in our work (Duivenvoorden et al. 2020). With the updated priors, the χ^{2} for Herschel is reduced to 222, but the χ^{2} for Planck data remains the same.
In none of the approaches does the χ^{2} value for the Herschel power spectra appear to fare as well as the Planck value. Similarly, using a different halo model with a larger number of parameters, Viero et al. (2013) did not obtain a good χ^{2} value (). We investigated this problem by testing the compatibility of the Planck and Herschel measurements. The details of the tests we performed are given in Appendix C.
In Fig. C.1 we show the measurements of the CIB power spectra from Planck/HFI (Planck Collaboration XXX 2014), Herschel/SPIRE (Viero et al. 2013), and Lenz et al. (2019) at 857 GHz. We also show the bestfit values of the onehalo, shot noise, and the total power spectrum when the model was fit to the Planck or Herschel data alone. The Herschel/SPIRE data points and the bestfit curves were scaled to obtain them as they would be measured using the Planck filters at 857 GHz. The procedure for this is provided in Appendix A.1 of Lagache et al. (2020). A similar trend is observed for the power spectra at 545 GHz. The Herschel and Planck CIB measurements agree very well at large scales and up to ℓ = 2000. The extrapolation of the Planck best fit at higher ℓ largely overestimates the Herschel data points, however. This mostly comes from the shot noise, which is not compatible for Planck and Herschel. While the difference of flux cuts (710 mJy for Planck and 300 mJy for Herschel) would give a variation of at most 12% of the shotnoise level^{3}, a factor of ∼1.9 is measured between the two. This discrepancy cannot be reconciled considering the difference in amplitude of the onehalo term. This discrepancy of the Herschel and Planck shotnoise measurements, and with the shotnoise values derived from models of galaxy number counts, has previously been pointed out by Lagache et al. (2020). Inconsistencies also exist inside a single experiment. For example, we show in Fig. C.2 the comparison of the CIB power spectra measured by Thacker et al. (2013) and Viero et al. (2013) for the same flux cut of 50 mJy at 857 GHz. The two measurements are clearly different at very high multipoles.
We explored two ways of reconciling the Planck and Herschel measurements using our halo model. In the first approach, we broadened the error bars on the calibration factors for the 600, 857, and 1200 GHz Herschel channels to 0.05 each with their central values at 1.00 instead of using the values obtained for the relative calibration between Planck and Herschel by Bertincourt et al. (2016), as done in our original approach. In this setting, we indeed find a very good fit to the Herschel power spectra, with a χ^{2} of 127 when fit together with Planck data (Table C.3), and 137 when the Herschel data are fit alone (Table C.2) compared to the χ^{2} of 221 in Table C.4 for 102 data points. However, the bestfit value for the calibration factors as shown in Tables C.3 and C.4 is too high. As mentioned before, the crosscalibration between Herschel/SPIRE and Planck/HFI has been measured very precisely at 545 and 857 GHz by Bertincourt et al. (2016) (the 545 GHz Planck channel has been calibrated to better than 1.5%, see Sect. 3.1). Thus the bestfit values obtained for the parameters are unrealistically high and cannot be accepted. In the second approach, we keep the original priors for the Herschel parameters from Bertincourt et al. (2016) and fit the data for low multipoles only, that is, ℓ < 3000. Results are provided in Table C.5 and Fig. C.3, which show that we obtain a good χ^{2} (36 for 42 data points), and the shotnoise levels for Planck and Herschel are quite similar (see Tables C.1 and C.5), which is what we expect. This shows that the Herschel data at ℓ > 3000 are not compatible with the values expected from the ℓ < 3000 Planck and Herschel data.
The mass of the dark matter haloes for converting the accreted baryons into stars most efficiently is log_{10}M_{max} . The highest efficiency mass found here is slightly on the high side of the range M_{⊙} to M_{⊙} found by Viero et al. (2013) and Planck Collaboration XXX (2014), but it agrees well with Chen et al. (2016) for faint SMGs of and using the linear clustering model for the CIB anisotropies (Maniyar et al. 2018). The error bars on the M_{max} are quite tight. Based on our best fit, we find that remaining three physical parameters of the model (η_{max}, σ_{Mh0}, and τ) are quite correlated amongst themselves (∼90%), while M_{max} is much less correlated with these parameters (∼10%). Furthermore, we did not vary M_{max} with redshift. This, combined with the fact that we have just four parameters in the model, appears to be the reason that the constraints on M_{max} are much tighter than those derived from other studies.
With this simple model, η_{max}, that is, the highest efficiency with which the accreted baryonic gas forms into stars, is ∼40%. We recall that we did not account for the effects such as quenching, AGN and supernovae feedback, and their evolution with redshift through explicit parameters. This efficiency parameter η_{max} can therefore be considered as the efficiency of converting the accreted baryons into stars after marginalising over these effects. Behroozi et al. (2013) reported that the baryonic accretion efficiency as defined by us, that is, SFR/BAR, varied between 20 and 40% across all redshifts, which is consistent with our result. Moreover, Moster et al. (2018) calculated the star formation efficiency as M_{*}(z)/M_{b}(z)=M_{*}(z)/(M_{h}(z) × f_{b}) = M_{*}(z)/M_{h}(z) × Ω_{m}/Ω_{b}, where M_{*}(z) is the stellar mass in a given halo at a given redshift. They reported that this value is 0.17 at redshift z = 0.1 and 0.2 at z = 2.0.
We derived a similar quantity from our model by integrating the mean SFR and BAR over cosmic times for various halo mass seeds at z = 10. We find a maximum efficiency of 0.34 at z = 0.1 and 0.37 at z = 2.0, which is higher than the values reported by Moster et al. (2018). However, in this approach, we did not take the mass loss from stellar evolution into account. Similar to Zahid et al. (2014), and assuming 45% mass loss, corresponding to what was found by Leitner & Kravtsov (2011) for an intermediateage galaxy, we find M_{*}(z)/M_{b}(z) = 0.19 at redshift z = 0.1 and 0.21 at z = 2. These values agree very well with the previous values of Moster et al. (2018).
Our model provides a nonzero value for τ, which means that it supports an evolution of the width of the lognormal parametrisation over time. This agrees with the observed suppression of star formation in massive haloes at low redshifts.
We tried different parametrisations for η instead of the lognormal. The details of one such parametrisation are provided in Appendix B. Because of its simplicity and because it provides physical results, we continue with the lognormal parametrisation.
Finally, as pointed out in the introduction, one of the shortcomings of the halo model used by Planck Collaboration XXX (2014) is that it is does not match the infrared SFRD derived from galaxies. A primary reason for this is that they did not include the independent measurements of the SFRD from galaxy surveys as priors in their likelihood analysis. Consequently, the result from their linear clustering model is not consistent with the halo model either. It is therefore a good consistency check of the model to verify whether the prediction from the linear model and the infrared SFRD measurements from galaxies (which are added as priors in our likelihood) match the SFRD derived from the best fit of the model. Figure 5 shows the SFRD constrained by our CIB halo model along with the SFRD constrained by the linear model presented in Maniyar et al. (2018). Although at low redshifts the SFRD constrained by the halo model is slightly higher than the measurements from galaxies or the constraints from the linear model, it is overall consistent with the latter two. Thus the SFRD obtained from the halo model is able to pass through the SFRD priors derived from galaxy surveys. This is again impressive, considering the fact that we modelled the CIB using only four parameters.
Fig. 5. Measurements of the infrared SFRD using galaxies (Madau & Dickinson 2014). The solid black (blue) line shows the SFRD as constrained by the linear (halo) model with the dotted black (dotted blue) lines showing the 1σ regions around it. We also show the SFRD value determined by Khusanova et al. completely independently using the data from the ALMA ALPINE large program (priv. comm.). 
4. tSZ halo model
While at lower multipoles (ℓ < 1000) the tSZ angular power spectrum is sensitive to the amplitude of the matter fluctuations, at higher multipoles (ℓ > 1000) the power spectrum also depends on the details of the pressure profile of the intracluster gas within the haloes. This means that along with the halo mass function, the model also has to account for the pressure profile of the gas.
The power spectrum of the tSZ is calculated as (e.g. Bolliet et al. 2018)
where f(ν) gives the frequency dependence of the tSZ. It is given by
where h_{p} is the Planck constant, and k_{B} is the Boltzmann constant. y_{ℓ}(M_{h}, z) represents the twodimensional Fourier transform of the electron pressure profile P_{e} for a halo of mass M_{h}. It is given as (e.g. Komatsu & Seljak 2001)
In this equation, c is the speed of the light, σ_{T} is the Thomson crosssection, m_{e} is the electron mass, x ≡ r/r_{500}, with r being the radial distance from the centre of the halo, r_{500} is the radius of the sphere that contains the overdensity mass M_{500c} of 500 times the critical density of the universe, and ℓ_{500} ≡ d_{A}/r_{500}, with d_{A} being the angular diameter distance.
As was done for the CIB modelling, we used the NFW profile to represent the density distribution inside the dark matter halo. Using this distribution, we have
where parameters (P_{0}, c_{500}, γ, α, β) were set to their bestfitting values obtained by Planck Collaboration Int. V (2013). They are equal to 6.41, 1.81, 0.31, 1.33, and 4.13, respectively. The coefficient C goes with the mass as
where H is the Hubble constant at redshift z, with H_{0} being its local value and h the reduced Hubble constant (h = H_{0}/100). The mass used here M̃_{500c} is not necessarily the true mass, but can contain a bias due to observational effects and nonthermal pressure, that is, there can be a difference between the true mass of the cluster and that obtained assuming hydrostatic equilibrium. In order to account for this possible bias, a variable B is often used, which relates the true mass M_{500c} to M̃_{500c} as M̃_{500c} = M_{500c}/B. In our analysis, we took the mass and redshift range to be the same as was used for the CIB model with x_{min} = 10^{−6} and x_{max} = 10 in Eq. (29).
Because massive clusters are rare, the contribution of the twohalo term to the total tSZ angular power spectrum is far lower than that of the onehalo term and can therefore be neglected (Komatsu & Kitayama 1999) for ℓ ≥ 300. It is important only for ℓ ≤ 300, but on these scales, the primary CMB anisotropies are clearly dominant. However, as we show in Sect. 5.3, the twohalo term plays a significant role in the twohalo term of the CIB × tSZ crosspower spectra and cannot be ignored. In our analysis, we therefore do not neglect this component.
It is straightforward to calculate the twohalo component of the power spectrum as (e.g. Salvati et al. 2018)
The total tSZ power spectrum is the sum of the onehalo and twohalo components:
f(ν) is negative for ν ≤ 217 GHz. This means that the tSZ power spectrum for certain choices of frequency combinations will be negative. The case for the CIB−tSZ correlation described in Sect. 5 is similar. To calculate the f(ν) parameter correctly for a broadband filter, we have to convolve the tSZ frequency dependence with the bandpass filter at that particular frequency. The f(ν) parameter is used to convert the dimensionless y parameter into CMB temperature units. Although the plots presented here for the tSZ power spectra are made independent of frequency (dimensionless) by dividing out the factors of f(ν), plots for the CIB−tSZ correlation are calculated at the value of the f(ν) at that particular frequency channel and are convolved with the filters. For 100, 143, 217, 343, 545, and 857 GHz corresponding to the Planck/HFI frequencies, f(ν) is −1.51, −1.04, −0.01, 2.24, 5.60, and 11.09, respectively, when it is not convolved with the respective bandpass filters. The effective values become −4.03, −2.79, 0.19, 6.21, 14.46, and 26.34, respectively, after convolving with the Planck/HFI bandpass filters at these respective frequencies, as given in Table 1 of Planck Collaboration XXII (2016).
Because the values of the parameters (P_{0}, c_{500}, γ, α, β) are not varied, this allows us to tabulate y_{ℓ}, which speeds up the process considerably. Figure 6 shows the prediction for the tSZ autopower spectrum with this model and with the given values of the parameters along with B = 1.41 as obtained by Bolliet et al. (2018) after fitting their model to the data from Planck Collaboration XXII (2016). The only frequency dependence of the tSZ power spectra comes from the f(ν) factors, and this has been divided out in the plot. Aat about ℓ ∼ 1000, the twohalo term contributes about 2% to the total power spectrum and the onehalo term dominates. In their Figs. 1 and 2, Bolliet et al. (2018) showed that the tSZ power spectra are very sensitive to the choice of the halo mass function, to the c − M_{h} relation (where c is the halo concentration parameter), and to cosmological parameters. This means taht several sources of uncertainties arise when the tSZ power spectra are calculated. In Fig. 6 we show the marginalised y power spectrum after subtraction of the foreground residuals as measured by Planck Collaboration XXII (2016). Bolliet et al. (2018) improved upon the Planck analysis by accounting for the trispectrum in the covariance matrix, and they provided new values for the y power spectrum after marginalising over the foreground residuals, which are shown in the figure as well. The bestfit model is consistent with these measurements. At ℓ = 257.5, using the bestfit parameters mentioned before, we obtain . For this ℓ, Planck Collaboration XXII (2016) provided the measured y power spectrum (10^{12}y^{2}) value as 0.217 ± 0.049, which is consistent with our predictions (the bestfit value from their model is 0.203). We placed the error bars on the Planck values as the quadrature sum of the statistical and the foreground errors provided by them. At the same time, Bolliet et al. (2018) provided a value of 0.179 ± 0.034, with a bestfit value of 0.131. This illustrates the range of values that the tSZ power spectra can have based on the different assumptions in the modelling.
Fig. 6. Predictions for the tSZ power spectrum based on our model for the given parameter values. The frequency dependence of the power spectrum through f(ν) has been divided out and the power spectrum is dimensionless. 
5. CIB−tSZ crosscorrelation
As pointed out in the introduction, CIB and tSZ are expected to be correlated to a certain degree. Thus in order to interpret this correlation, we need a crosscorrelation model that predicts this signal. In the sections before, we presented the halo models for calculating the CIB and tSZ power spectra, and we use them below to calculate the CIB−tSZ correlation in a consistent halo model formalism.
5.1. Halo model formalism
The onehalo term provides the correlation between the CIB sources and the tSZ within the same halo. On the other hand, the twohalo term arises from the correlation between the CIB sources in one halo with the tSZ in other halo. This formulation has a very interesting consequence. In an extreme scenario where a very massive halo contributes significantly to the tSZ without star formation, it hence has no CIB sources. In this case, the onehalo term would be zero, but the twohalo term would still contribute provided that there is some overlap in the redshift distribution of the CIB sources and the tSZ haloes. This means that the twohalo term, which dominates on large angular scales, does not depend significantly on the astrophysical processes governing the star formation in massive haloes contributing to the tSZ.
The corresponding one and twohalo terms are given as (for k = ℓ/χ)
and
where
M_{h} and represent the mass of the two different haloes between which the correlation is computed when the twohalo term is computed. The tSZ spectral variation does not dependent on the redshift or halo mass in the nonrelativistic limit, therefore the cross spectra are simpler than they would be otherwise.
The total CIB × tSZ power is then
When the CIB and tSZ model parameters are known, it is straightforward to calculate the CIB−tSZ power spectra and there is no need of additional parameters. This is one of the advantages of the halo model approach over previous studies (e.g. Dunkley et al. 2013; George et al. 2015; Reichardt et al. 2020; Choi et al. 2020), where a templatebased approach for the CIB−tSZ correlation was used by introducing a parameter ξ to act as a correlation coefficient. This templatebased approach is also limited by the fact that ξ might be angular and scale dependent, and fitting for a single parameter ξ might therefore not produce a true picture of the CIB−tSZ correlation. This has indeed been shown to be the case by Addison et al. (2012), who found nonnegligible variation in ξ with the angular scales and the considered frequencies.
5.2. Halo mass definition
Our models for CIB and tSZ (and thus CIB × tSZ) rely on the halo mass M_{h}. For the tSZ effect, the pressure profile from Eq. (30) is in general given in terms of M_{500}, which is defined as the mass contained within a radius in which the mean overdensity is 500 times the background critical density. In contrast, the CIB studies are generally carried out with M_{200} as the definition of the halo mass. We have to be consistent with the definition of the halo mass while considering the CIB and tSZ together. One way of doing this is calculating the CIB terms with M_{200}, then converting the M_{200} into corresponding M_{500} for the given redshift range and cosmology using the procedure given in Komatsu & Seljak (2001), and then using this M_{500} to calculate the corresponding tSZ terms. We presented the CIB × tSZ results here, we worked with a single definition of M_{500} for the halo masses for the CIB and tSZ terms. We verified the effect of using M_{500} instead of M_{200} on the CIB model parameters, and we found that the values of the CIB parameters slightly change. The most notable change was observed in the value of η_{max}, which was found to be higher when we used M_{500} instead of M_{200}: because with this definition of the mass, the CIB contribution is calculated very close to the centre of the cluster and does not include the star formation that occurs far away. The model accordingly tries to increase the conversion efficiency of the accreted baryons into stars to produce the same level of CIB as observed. The CIB model parameters given in Table 1 were calculated using the commonly used M_{200} definition for the halo masses.
5.3. Predictions for the CIB × tSZ power spectra
Figure 7 shows the predictions for the CIB, tSZ, and CIB × tSZ power spectra correlating the Planck 143 GHz frequency with other frequencies based on the models of the CIB, tSZ, and CIB × tSZ given above. We note again that our results were derived using the M_{500} definition for the halo masses, and for the tSZ, the f(ν) factor was convolved with the Planck/HFI frequencies. This provided values of −4.03, −2.79, 0.19, 6.21, 14.46, and 26.34 for 100, 143, 217, 343, 545, and 857 GHz, respectively, as mentioned in Sect. 4. This means that f(ν) is negative for ν < 217 GHz and the tSZ power spectra are negative for certain choices of the frequency pairs. In the figure, the tSZ power spectra are negative for a combination of 143 GHz with all the frequencies above 143 GHz, that is, 217, 353, 545, and 857 GHz (we show their absolute values). Similarly, CIB × tSZ power spectra shown here are negative in all the cases from 100 to 857 GHz.
Fig. 7. Predictions for the different components of the CIB, tSZ, and CIB × tSZ power spectra correlating the Planck 143 GHz frequency with other frequencies based on the models given above. 
A special case to consider when Eqs. (34) and (35) are used is that of the two widely separated frequency channels. Figure 7 shows that the relative power of the CIB−tSZ to CIB increases with the frequency separation for frequencies ν ≥ 217 GHz. In this case, one of the terms in the square parentheses is generally far smaller than the other and can be neglected because the CIB and tSZ depend differently on frequency. The crosscorrelation is then driven by the dominating term. The CIB at higher frequencies traces the haloes at lower redshifts, and this increases the overlap with the tSZ clusters. This results in the increase in the CIB−tSZ power compared to the CIB power spectra as the frequency separation increases. This scenario can provide a good opportunity to constrain the CIB−tSZ crosspower spectra. The tSZ term dominates at lower frequency channels and drives the crosscorrelation. The tSZ has a null point at ν ≈ 217 GHz. The figure shows this, and the location in which the tSZ power spectra in 143 × 217 GHz are far smaller than for 143 × 143 and 143 × 353 GHz.
The observational constraints for the tSZ and especially the CIB have improved dramatically since results from Addison et al. (2012). This has resulted in an improved understanding and modelling of both the tSZ and the CIB, which we take advantage of. Our predictions for the CIB × tSZ correlation for a given frequency channel pair are therefore not exactly the same as those obtained by Addison et al. (2012) in terms of their amplitude and/or shape.
5.4. Redshift contributions to the power spectra
Figure 8 shows the contributions to the CIB−tSZ power spectra from different redshift bins. The power spectra are shown for 100−100 GHz and 100−857 GHz channel pairs.
Fig. 8. Redshift contribution to the power spectra of the CIB−tSZ correlation at for 100−100 GHz channel pair (left panel) and 100−857 GHz channel pair (right panel). All the values shown here are absolute values. 
The clusters hosting hot gas with energetic electrons that cause the tSZ effect reside at relatively low redshifts (∼z < 1). We therefore expect the contribution of the onehalo term of the power spectrum to the total power spectrum to relatively decrease compared to the twohalo term at higher redshifts. This is indeed what we observe in Fig. 8, where the multipole where the onehalo and twohalo term contribute equally moves slightly to higher values in higher redshift bins, that is, the range in which the twohalo term dominates the onehalo term increases with redshift. This of course depends upon the width of the redshift bins considered and on the pair of frequency channels considered, but the overall trend is the same.
Figure 8 shows that although majority of the tSZ clusters are expected to reside at lower redshifts, the CIB−tSZ contribution coming from the lowest redshift bin (0 < z < 0.5) is smaller than or equal to that from other redshift bins. For the 100−100 GHz channel pair, our model predicts that the CIB−tSZ contribution coming from 0 < z < 0.5 is of the same order as the z > 3 bin where we do not expect to find many tSZ clusters, while for the 100−857 GHz channel pair, the contribution from the 0 < z < 0.5 redshift bin is higher than the z > 3 bin. The CIB at lower frequencies traces the galaxies at higher redshift to a certain extent, and vice versa. In the 100−100 GHz case, the CIB contribution mainly comes from galaxies at relatively higher redshift (e.g. see Fig. 4 of Maniyar et al. 2018), and thus the twohalo term caused by the overlap in the redshift distribution of the CIB sources and the tSZ haloes gives a slightly higher power in the z > 3 redshift bin than the 0 < z < 0.5 bin in the 100−100 GHz case. The case for the 100−857 GHz pair is exactly the reverse: the CIB is mainly fed by galaxies at lower redshifts, thereby increasing the power in the lowest redshift bin.
Most of the contribution to the CIB−tSZ power spectra appears to stem from the intermediate redshift bins, that is, 0.5 < z < 1.5 and 1.5 < z < 3.0. Going back to Fig. 4 of Maniyar et al. (2018) again, this is expected because most of the CIB power is coming from the dusty star forming galaxies between 0.5 < z < 4. Thus the CIB part driving the CIB−tSZ correlation has most of its contribution coming from these redshifts. Again, we can see that going from 100−100 GHz pair to 100−857 GHz pair, the contribution of the 0.5 < z < 1.5 slightly increases in comparison to 1.5 < z < 3 bin as even lower redshifts are probed by the CIB at 857 GHz than at 100 GHz.
5.5. Angular scale dependence of the CIB × tSZ power spectra
Measuring the kSZ power spectrum from the CMB power spectrum is one of the challenging goals of current CMB cosmology data analysis. It is important to measure the kSZ power spectrum because it can help us to understand the reionisation history of the Universe better. On very small angular scales (ℓ > 2000) in the CMB power spectrum, the kSZ power spectrum starts to show a similar amplitude as other components such as the CIB, tSZ, and their crosscorrelation, which act as foregrounds to correctly measure the kSZ. Figures 5 and 6 from George et al. (2015) show that the kSZ is degenerate with the tSZ and CIB × tSZ correlation. A consistent modelling of these foregrounds is highly desirable.
The kSZ has been constrained by Dunkley et al. (2013), George et al. (2015), Reichardt et al. (2020), and Choi et al. (2020) using the CMB powerspectrum measurements from the Atacama Cosmology Telescope (ACT) or the South Pole Telescope (SPT), while Planck Collaboration XXIII (2016) reported a measurement of the CIB × tSZ correlation using Planck data. George et al. (2015) and Reichardt et al. (2020) calculated the CIB × tSZ power spectra at every step of their analysis using the CIB and tSZ power spectra obtained with templates at this step, whereas Dunkley et al. (2013) and Planck Collaboration XXIII (2016) used the CIB × tSZ template provided by Addison et al. (2012). All of them fit for a single scaleindependent amplitude parameter for the CIB × tSZ correlation across all the frequency channels. Addison et al. (2012) used the CIB model from Xia et al. (2012) to calculate the CIB × tSZ crosscorrelation. One of the shortcomings of this model is that it assumes the spectral properties of the CIB, that is, SEDs, luminosity to be independent of the host halo mass. Our understanding of the CIB and its properties has improved significantly in the last few years, and we know that the spectral properties of the CIB are highly dependent on the host halo mass. We therefore compare the CIB × tSZ correlation calculated with our newly developed CIB halo model with that developed by Addison et al. (2012) using the CIB model from Xia et al. (2012).
Figure 9 shows the ratio of the CIB × tSZ template from Addison et al. (2012) used in the Planck Collaboration V (2020) likelihood analysis and the corresponding power spectra calculated within our halo model framework at 217 × 217 GHz in red. The two power spectra were normalised such that they have equal value at ℓ = 3000. The ratio of the CIB × tSZ template used in the Planck Collaboration V (2020) with our model is different at all scales (e.g. 1.1 at ℓ ∼ 2000 and 0.9 at ℓ ∼ 3700). The results for the kSZ power spectra derived using these templates might therefore be different than if the halo model developed in this paper were used. In Fig. 9 we also plot the ratio of the CIB × tSZ power spectra for different pairs of Planck frequency channels with that at 217 × 217 GHz calculated using our halo model, again normalised such that they have the same value at ℓ = 3000. Although at very high multipoles the ratio does not vary strongly because the CIB at these nearby frequency channels is highly correlated, it does vary by up to 10% for ℓ < 2000. Addison et al. (2012) also showed that for widely separated frequency channels, the CIB × tSZ correlation is different at different multipoles. Fitting for a single amplitude parameter independent of the multipoles for the CIB × tSZ power spectra across all frequency channels is therefore a crude approximation. The bestfit templates used in the previous analysis should be replaced by such physically motivated halo models. This will matter even more for future CMB data that will have a very high signaltonoise ratio on the small angular scales of interest for the kSZ power spectrum.
Fig. 9. Ratio of the CIB × tSZ power spectra at different frequencies and at 217 × 217 GHz. The red curve shows the ratio with the CIB × tSZ template from Addison et al. (2012) used in the Planck Collaboration V (2020) likelihood analysis. All the power spectra are normalised such that they have the same value at ℓ = 3000. 
5.6. Correlation coefficient for CIB × tSZ
We define the correlation coefficient for the CIB and tSZ correlation as
where the correlation coefficient ξ_{ℓ}(ν) depends on frequency and angular scale. The CIB power spectra also include the contribution from the shot noise at that particular frequency. In Fig. 10 we show ξ_{ℓ}(ν) for 143 GHz Planck frequency channel based on the CIB, tSZ, and the CIB × tSZ bestfit models in this paper. We did not fit for the shot noise of the CIB at 143 GHz. We took this value as calculated by the model from Béthermin et al. (2012), which is 0.87 Jy^{2} sr^{−1}. This was added to the one and twohalo terms at 143 GHz. The CIB × tSZ correlation at 143 GHz is negative, so that we took its absolute value. Based on the halo model formalism in Addison et al. (2012), they provide the value of ξ_{ℓ} for crosscorrelation of 150 GHz SPT frequency channel with other frequencies. They also provide the factors to convert these values into those that would be obtained for the corresponding Planck 143 GHz channel that we used here. George et al. (2015) and Reichardt et al. (2020) measured the value for ξ_{ℓ} at ℓ = 3000 using the 95, 150, and 220 GHz SPT frequency channels. We show these points in Fig. 10.
Fig. 10. Correlation coefficient for the CIB × tSZ crosscorrelation calculated with the halo model for the 143 GHz Planck channel. We also show the corresponding measurements for the correlation coefficient from Addison et al. (2012), George et al. (2015), and Reichardt et al. (2020) for SPT frequency channels at ℓ = 3000. They are slightly offset on the xaxis for clarity. 
The correlation coefficient ξ_{ℓ} is defined by Addison et al. (2012), George et al. (2015), and Reichardt et al. (2020) in a slightly different manner. The values shown in the plot here are twice as high as the values provided by them because of this change in definition. Moreover, as mentioned in Sect. 5.1, George et al. (2015) and Reichardt et al. (2020) used a templatebased approach to calculate the CIB and the tSZ power spectra. They then used these values to calculate the CIB × tSZ correlation by fitting for the crosscorrelation coefficient ξ. They therefore obtained a single value of ξ_{ℓ} for the power spectra across all the frequencies, that is, 95, 150, and 220 GHz across all multipoles. These values were additionally calculated using a different telescope (SPT), which corresponds to a different flux cut (∼6.4 mJy compared to 710 mJy for Planck). Although we can compare the ξ_{ℓ} provided by Addison et al. (2012) with our predictions, it is therefore not a direct comparison with the values from George et al. (2015) and Reichardt et al. (2020). We show these points just to give an idea of the value they derived for ξ_{ℓ}.
George et al. (2015) and Reichardt et al. (2020) reported a value of and , respectively, at ℓ = 3000. Our bestfit model gives a value of ξ_{ℓ} of 0.45, which is higher than the value at ℓ = 3000 reported by Addison et al. (2012), which is 0.39. There should be an error bar on ξ_{ℓ} calculated here and by Addison et al. (2012) because of the uncertainties corresponding to the CIB and the tSZ halo models used. The calculation of these uncertainties is left for future work.
6. Conclusions
One of the main motivations of this work was designing a consistent framework for calculating the CIB, tSZ, and CIB × tSZ power spectra in a halo model setting. For this purpose, we developed a simple and physically motivated halo model of the CIB with only four parameters describing the relationship between the mass of the dark matter haloes and their efficiency to convert the accreted baryons into stars using a lognormal parametrisation. Because previous evidence showed that massive haloes do not contribute significantly at lower redshift to the total star formation budget but are efficient at high redshift, we allowed the width of the lognormal to evolve with redshift. We find that the mass of the darkmatter haloes for the highest efficiency is M_{⊙}, while the maximum efficiency at this mass is . The mass of the highest efficiency found here is slightly on the high side of the range M_{⊙} to M_{⊙} found by Viero et al. (2013) and Planck Collaboration XXX (2014), but it agrees well with Chen et al. (2016) for faint SMGs of and using the linear clustering model for the CIB anisotropies (Maniyar et al. 2018). This agreement is quite motivating considering the simple nature of our model. The model is also able to fit the SFRD measured using the galaxies and is consistent with the SFRD history obtained using the linear clustering model from Maniyar et al. (2018). The SFRD constraints derived from the galaxy surveys are added for the first time as priors in the likelihood for the halo model. This helps us achieve the significant result of reproducing the CIB power spectra and the SFRD history at the same time with this model. We also calculated the CIB−CMB lensing crosscorrelation using the bestfit value of the CIB parameters and find that it agrees well with the measurements (Fig. 4). While we obtain a decent χ^{2} of 113 for 80 Planck data points, the χ^{2} for Herschel power spectra comes is 247 for 102 data points, which is not good. We investigated this discrepancy using different methods and found that the CIB power spectrum measurements are not fully compatible with each other for ℓ > 3000.
For the tSZ, we followed the approach of Bolliet et al. (2018) and developed a halo model to calculate the one and twohalo power spectra. Because we kept the parameters for the pressure profile of the gas constant, the power spectra calculations were sped up by tabulating the values of y_{ℓ}. The only parameter for the tSZ power spectra is the mass bias B, which takes the difference between the true halo mass and the mass derived assuming the hydrostatic equilibrium, for example, into account. Although we find that the contribution of the twohalo term is not significant ∼2% at ℓ ∼ 1000 to the total power spectra, we still considered it in our calculations for a complete analysis and because the twohalo term has a significant contribution to the twohalo term of the CIB × tSZ, which is significant and cannot be ignored.
Following our two halo models for the CIB and tSZ, we computed the crosscorrelation between the CIB and tSZ within a halo model framework. In addition to being consistent with the CIB and the tSZ halo models, another advantage of this approach is that we do not require an additional parameter to calculate the correlation. When the CIB and tSZ model parameters are known, it is straightforward to calculate the CIB−tSZ correlation. This is quite advantageous compared with other studies that used the template approach for the power spectra and a fit for the global amplitude. Based on our model and the best fit of the CIB and tSZ parameters, we find that the relative power of the CIB−tSZ correlation with respect to the CIB and tSZ increases with frequency separation of the maps. The twohalo term for the CIB−tSZ correlation is also not negligible with respect to the onehalo term and should be considered in the total power spectrum calculation. This twohalo term arises from the correlation between the tSZ, which is fed from one halo, and the CIB galaxies, which are located in the other halo. Moreover, most of the contribution to the CIB−tSZ power spectra comes from 0.5 < z < 4 because most of the power of the CIB part that contributes to the CIB−tSZ comes from these redshifts.
CIB and tSZ act as foregrounds in CMB anisotropy measurements. Especially on small scales, the CIB, tSZ, and CIB × tSZ correlation acts as a hindrance for measuring the kinematic SZ (kSZ) signal, which is the SZ effect caused by the peculiar velocities of the galaxy clusters after reionisation and by the ionised bubbles during the reionisation. In order to measure the kSZ from the CMB power spectrum, it is therefore important to correctly model and remove these foregrounds. Previous analyses that measured the kSZ power spectra have used a templatebased approach to model the foregrounds and fit for a single amplitude parameter across all the frequency channels independent of the multipole range considered. We showed that this approach is not sufficient to capture the scale dependence of the CIB × tSZ power spectra, which is moreover model dependent. A similar argument is also applicable to modelling the other foregrounds, that is, the CIB and tSZ as well. The kSZ constraints obtained using the template approach have thus to be revised, replacing the templates by power spectra computed by the physically developed halo models. Because our halo model is effective, it might be used for this purpose.
Code for calculating the CIB, tSZ, and CIB × tSZ power spectra is made available online at https://github.com/abhimaniyar/halomodel_cib_tsz_cibxtsz
Acknowledgments
We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSUIN2P3INP, CEA and CNES, France and support from the OCEVU Labex (ANR11LABX0060) and the A*MIDEX project (ANR11IDEX000102) funded by the “Investissements d’Avenir” French government program managed by the ANR. GL has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 788212) and funding from the Excellence Initiative of AixMarseille UniversityA*Midex, a French “Investissements d’Avenir” programme. AM warmly thanks Marco Tucci for useful discussions and inputs on the halo modelling of the CIB power spectra.
References
 Addison, G. E., Dunkley, J., & Spergel, D. N. 2012, MNRAS, 427, 1741 [NASA ADS] [CrossRef] [Google Scholar]
 Amblard, A., & Cooray, A. 2007, ApJ, 670, 903 [NASA ADS] [CrossRef] [Google Scholar]
 Amblard, A., Cooray, A., Serra, P., et al. 2011, Nature, 470, 510 [NASA ADS] [CrossRef] [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Bertincourt, B., Lagache, G., Martin, P. G., et al. 2016, A&A, 588, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Dole, H., Lagache, G., Le Borgne, D., & Penin, A. 2011, A&A, 529, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthermin, M., Wu, H.Y., Lagache, G., et al. 2017, A&A, 607, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bolliet, B., Comis, B., Komatsu, E., & MacíasPérez, J. F. 2018, MNRAS, 477, 4957 [NASA ADS] [CrossRef] [Google Scholar]
 Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
 Challinor, A., & Lewis, A. 2005, Phys. Rev. D, 71, 103010 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C.C., Smail, I., Swinbank, A. M., et al. 2016, ApJ, 831, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, S. K., Hasselfield, M., Ho, S. P. P., et al. 2020, ArXiv eprints [arXiv:2007.07289] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Cousin, M., Lagache, G., Bethermin, M., & Guiderdoni, B. 2015, A&A, 575, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cousin, M., Guillard, P., & Lehnert, M. D. 2019, A&A, 627, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686 [NASA ADS] [CrossRef] [Google Scholar]
 DessaugesZavadsky, M., Zamojski, M., Schaerer, D., et al. 2015, A&A, 577, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duivenvoorden, S., Oliver, S., Béthermin, M., et al. 2020, MNRAS, 491, 1355 [Google Scholar]
 Dunkley, J., Calabrese, E., Sievers, J., et al. 2013, JCAP, 2013, 025 [NASA ADS] [CrossRef] [Google Scholar]
 Fakhouri, O., Ma, C.P., & BoylanKolchin, M. 2010, MNRAS, 406, 2267 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 George, E. M., Reichardt, C. L., Aird, K. A., et al. 2015, ApJ, 799, 177 [Google Scholar]
 Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1 [NASA ADS] [Google Scholar]
 Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 436, 2875 [NASA ADS] [CrossRef] [Google Scholar]
 Hanson, D., Hoover, S., Crites, A., et al. 2013, Phys. Rev. Lett., 111, 141301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]
 Knox, L., Cooray, A., Eisenstein, D., & Haiman, Z. 2001, ApJ, 550, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., & Kitayama, T. 1999, ApJ, 526, L1 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Komatsu, E., & Seljak, U. 2001, MNRAS, 327, 1353 [Google Scholar]
 Lagache, G., & Puget, J. L. 2000, A&A, 355, 17 [NASA ADS] [Google Scholar]
 Lagache, G., Bavouzet, N., FernandezConde, N., et al. 2007, ApJ, 665, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Lagache, G., Bethermin, M., Montier, L., Serra, P., & Tucci, M. 2020, A&A, 642, A232 [CrossRef] [EDP Sciences] [Google Scholar]
 Leitner, S. N., & Kravtsov, A. V. 2011, ApJ, 734, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Lenz, D., Doré, O., & Lagache, G. 2019, ApJ, 883, 75 [CrossRef] [Google Scholar]
 Limber, D. N. 1954, ApJ, 119, 655 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maniyar, A. S., Béthermin, M., & Lagache, G. 2018, A&A, 614, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maniyar, A., Lagache, G., Béthermin, M., & Ilić, S. 2019, A&A, 621, A32 [CrossRef] [EDP Sciences] [Google Scholar]
 Marchetti, L., Vaccari, M., Franceschini, A., et al. 2016, MNRAS, 456, 1999 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuhara, H., Kawara, K., Sato, Y., et al. 2000, A&A, 361, 407 [NASA ADS] [Google Scholar]
 Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Nature, 556, 469 [NASA ADS] [CrossRef] [Google Scholar]
 Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
 Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2016, A&A, 594, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. V. 2013, A&A, 550, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVI. 2016, A&A, 596, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Popesso, P., Biviano, A., Finoguenov, A., et al. 2015, A&A, 579, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puget, J.L., Abergel, A., Bernard, J.P., et al. 1996, A&A, 308, L5 [NASA ADS] [Google Scholar]
 Reichardt, C. L., Shaw, L., Zahn, O., et al. 2012, ApJ, 755, 70 [Google Scholar]
 Reichardt, C. L., Patil, S., Ade, P. A. R., et al. 2020, ApJ, submitted [arXiv:2002.06197] [Google Scholar]
 Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Salvati, L., Douspis, M., & Aghanim, N. 2018, A&A, 614, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shang, C., Haiman, Z., Knox, L., & Oh, S. P. 2012, MNRAS, 421, 2832 [NASA ADS] [CrossRef] [Google Scholar]
 Silk, J. 2003, MNRAS, 343, 249 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Thacker, C., Cooray, A., Smidt, J., et al. 2013, ApJ, 768, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J. L., & Wetzel, A. R. 2010, ApJ, 719, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 [NASA ADS] [CrossRef] [Google Scholar]
 van den Bosch, F. C., More, S., Cacciato, M., Mo, H., & Yang, X. 2013, MNRAS, 430, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Viero, M. P., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1766 [NASA ADS] [CrossRef] [Google Scholar]
 Viero, M. P., Wang, L., Zemcov, M., et al. 2013, ApJ, 772, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, T., Elbaz, D., Daddi, E., et al. 2018, ApJ, 867, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Xia, J.Q., Negrello, M., Lapi, A., et al. 2012, MNRAS, 422, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Zahid, H. J., Dima, G. I., Kudritzki, R.P., et al. 2014, ApJ, 791, 130 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: The SFRD prior importance
Along with providing a good fit to the CIB power spectra, the halo model should be able to reproduce the SFRD history measured by extrapolating the galaxy luminosity functions. For this reason, we used the SFRD measurements from galaxies at different redshifts as priors while fitting for the halo model parameters. The bestfit parameters then result in the SFRD history, as shown in Fig. 5, which is consistent with the external measurements.
Figure A.1 shows the results when these external measurements are not considered as priors when the bestfit parameters of the CIB model are defined. In this case, the χ^{2} value obtained is 90 for 80 Planck data points and 156 for 102 Herschel data points. However, the bestfit value for the calibration factor for 1200 GHz , which is negatively correlated with the shotnoise at 1200 GHz, is 0.74 against the Gaussian prior set upon it centred at 1.00 with 1 σ error bar of 0.05. This χ^{2} is much better especially for Herschel data than the case when we considered the external SFRD measurements as priors in our likelihood. Although we are able to fit the CIB power spectra much better with these parameters, it is therefore evident that they predict an excessive SFRD, at least at low redshifts. It is therefore quite important to include these measurements as priors in our model, which comes at the expense of a poor χ^{2} value.
Fig. A.1. Measurements of the SFRD using galaxy surveys (Madau & Dickinson 2014). The solid black line shows the SFRD as measured by the CIB linear model of Maniyar et al. (2018). The solid blue line shows the corresponding constraints when the SFRD values from the galaxies are not considered as priors while performing the fit of the CIB halo model. The dotted black and red lines show the 1σ regions. 
Appendix B: Alternate parametrisations
In addition to the lognormal parametrisation considered for η, we studied several different parametrisations. One was inspired by the approach taken by Moster et al. (2013), who connected the galaxy mass with the corresponding host halo mass through a double powerlaw parametrisation. For our case, we used a double power law to describe the relation between η and M_{h}, which reads
where β, and γ describe the η at the low and highmass end, respectively. These two parameters allow us to have asymmetrical distribution around the mass of maximum efficiency. We fit this parametrisation with all the data sets and external constraints as mentioned in Sect. 3. The data used are not sensitive and thus cannot constrain the lowmass end slope β. Therefore we fixed the value of β and let γ evolve with redshift (for z ≤ 1.5) to avoid significant contribution by very massive haloes to the SFR at low redshift, in line with the reasoning presented in Sect. 2.1. Although this parametrisation gave a good fit to the data (even better than the fiducial model using lognormal) and was able to produce a decent SFRD history, the contribution from the haloes above the mass of maximum efficiency at high redshift was unrealistically low. Even though this parametrisation provided a better fit to the data than the fiducial model, we therefore continued with the lognormal because the results were more physical.
Appendix C: Comparison of the Planck and Herschel CIB power spectra
Here we show some figures and plots that compare the measurements of the CIB power spectra from Planck and Herschel as well as the bestfit values for the halo model obtained under different conditions. The shotnoise values provided here have to be multiplied with the corresponding colour corrections at each frequency to obtain them in the νI_{ν} = constant convention.
Figure C.1 shows the measured Planck 545 GHz and SPIRE 600 GHz autopower spectra in the same plot. On large angular scales, both data sets appear to agree well. However, on small scales, the two data sets differ, which appears to be mostly a result of the inconsistent shotnoise values of Planck and Herschel. A similar trend appears for the measurements from Lenz et al. (2019). A more detailed discussion of this can be found in Lagache et al. (2020) and is also presented in Sect. 3.4.
Fig. C.1. Measurements of the CIB power spectra from Planck Collaboration XXX (2014), Viero et al. (2013), and Lenz et al. (2019). The bestfit models shown in this figure are obtained when Planck and Herschel data are fit separately. The SPIRE data have been rescaled to Planck data using the crosscalibration factor derived in Bertincourt et al. (2016) and the factor that allows converting the measurement through the SPIRE bandpass into a measurement as it would be obtained through the Planck bandpass from Lagache et al. (2020). Overall, . 
In Fig. C.2 we show the comparison of the measured CIB autopower spectra at 857 GHz for SPIRE by Thacker et al. (2013) and Viero et al. (2013) for the same flux cut. Similar to the comparison with Planck, the two measurements agree well on the large angular scales. However, they are clearly different on small angular scales, but we expect them to be the same. This shows that there are inconsistencies within the analysis of a single experiment.
Fig. C.2. Measurements of the CIB power spectra from Thacker et al. (2013) and Viero et al. (2013) for a flux cut of 50 mJy. 
In order to determine whether the Planck and Herschel data can be reconciled within our model, we performed some tests, the results of which are provided in the following tables. Instead of fitting the Planck and Herschel data together with our model, we only fit for the Planck data and present the bestfit values in Table C.1. In this case, instead of assuming a fixed perfect calibration for the Planck experiment, we let the calibration factors at all the frequencies vary with a Gaussian prior centred at 1.00 with 1σ error of 0.05. The Planckχ^{2} value improves compared to when we fit for both the Planck and Herschel data together. The posterior values of calibration parameters are very close to one, and it is therefore a justified assumption to keep the Planck calibration fixed while fitting for both the data together.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Planck data considering the parameters at all Planck frequency channels centred at one with error bars of 0.05.
Similar to the previous test, we then fit our model to Herschel data alone. The results are shown in Table C.2. Unlike the Planck data, we obtain a poor χ^{2} value in this case. The bestfit values for the shot noise as well as the calibration parameters are consistent with the case when we fit for both the Planck and Herschel data together. There is some shift in the physical parameters of the halo model, but there is no particular trend as these parameters are correlated with each other.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data considering the parameters at all Herschel frequency channels, as done in our original analysis.
In order to verify the effect of the crosscalibration between Planck and Herschel, we performed several tests and show the results in Tables C.3 and C.4. In this case, instead of setting tight crosscalibration priors on the Herschel calibration parameters, we allowed them to vary with larger Gaussian priors with 1σ error bars of 0.05 centred at 1.00. In the first test, we fit for both the Planck and Herschel data, and in the second case, we fit only for the Herschel data. This test indeed provided very good χ^{2} values for both Planck and Herschel. However, it comes at the cost of unrealistically high values of the calibration parameters at 545 and 857 GHz that are completely inconsistent with the very precise crosscalibration measurement between Planck and Herschel (Bertincourt et al. 2016). This test therefore does not solve the compatibility problem of the two data sets.
Marginalised values of all the model parameters given at a 68% confidence level when we fit for both the Planck and Herschel data considering the parameters at all Herschel frequency channels centred at one with error bars of 0.05.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data considering the parameters at all Herschel frequency channels centred at one with error bars of 0.05.
We determined the consistency of the Planck and Herschel data on the large angular scales by fitting for the Herschel data at ℓ < 3000. The results are shown in Table C.5 and Fig. C.3. In this case, we obtain a decent χ^{2} value, and the shotnoise values are similar to the levels obtained for Planck; this was expected. This tests shows that the Herschel and Planck data appear to be compatible with each other on large angular scales, but they differ on small scales.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data upto ℓ < 3000 considering the parameters at all Herschel frequency channels as done in our original approach.
Fig. C.3. Measurements of the CIB auto and crosspower spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the bestfit CIB halo model with its different components when Herschel/SPIRE data alone were fit for ℓ < 3000. 
All Tables
Marginalised values of all the model parameters given at a 68% confidence level.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Planck data considering the parameters at all Planck frequency channels centred at one with error bars of 0.05.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data considering the parameters at all Herschel frequency channels, as done in our original analysis.
Marginalised values of all the model parameters given at a 68% confidence level when we fit for both the Planck and Herschel data considering the parameters at all Herschel frequency channels centred at one with error bars of 0.05.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data considering the parameters at all Herschel frequency channels centred at one with error bars of 0.05.
Marginalised values of all the model parameters given at a 68% confidence level when we fit only for the Herschel data upto ℓ < 3000 considering the parameters at all Herschel frequency channels as done in our original approach.
All Figures
Fig. 1. Lognormal parametrisation (Eq. (5)) between the halo mass (M_{⊙}) and ratio between the SFR and the baryonic accretion rate, η. We show two extreme cases: haloes near the efficiency peak contain subhaloes with very low mass (case 1), and very massive haloes that contain subhaloes near the efficiency peak (case 2). If the same recipe were used to calculate the satellite galaxy SFR in these two cases, unphysical values might result within the assumptions of our model, and therefore we suggest two different ways to calculate the SFR for satellite galaxies (Sect. 2.2). 

In the text 
Fig. 2. Measurements of the CIB auto and crosspower spectra obtained by Planck/HFI (Planck Collaboration XXX 2014) and the bestfit CIB halo model with its different components. 

In the text 
Fig. 3. Measurements of the CIB auto and crosspower spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the bestfit CIB halo model with its different components. 

In the text 
Fig. 4. CIB × CMB lensing crosspower spectra measured by Planck Collaboration XVIII (2014) in orange and our bestfit model in blue curves. 

In the text 
Fig. 5. Measurements of the infrared SFRD using galaxies (Madau & Dickinson 2014). The solid black (blue) line shows the SFRD as constrained by the linear (halo) model with the dotted black (dotted blue) lines showing the 1σ regions around it. We also show the SFRD value determined by Khusanova et al. completely independently using the data from the ALMA ALPINE large program (priv. comm.). 

In the text 
Fig. 6. Predictions for the tSZ power spectrum based on our model for the given parameter values. The frequency dependence of the power spectrum through f(ν) has been divided out and the power spectrum is dimensionless. 

In the text 
Fig. 7. Predictions for the different components of the CIB, tSZ, and CIB × tSZ power spectra correlating the Planck 143 GHz frequency with other frequencies based on the models given above. 

In the text 
Fig. 8. Redshift contribution to the power spectra of the CIB−tSZ correlation at for 100−100 GHz channel pair (left panel) and 100−857 GHz channel pair (right panel). All the values shown here are absolute values. 

In the text 
Fig. 9. Ratio of the CIB × tSZ power spectra at different frequencies and at 217 × 217 GHz. The red curve shows the ratio with the CIB × tSZ template from Addison et al. (2012) used in the Planck Collaboration V (2020) likelihood analysis. All the power spectra are normalised such that they have the same value at ℓ = 3000. 

In the text 
Fig. 10. Correlation coefficient for the CIB × tSZ crosscorrelation calculated with the halo model for the 143 GHz Planck channel. We also show the corresponding measurements for the correlation coefficient from Addison et al. (2012), George et al. (2015), and Reichardt et al. (2020) for SPT frequency channels at ℓ = 3000. They are slightly offset on the xaxis for clarity. 

In the text 
Fig. A.1. Measurements of the SFRD using galaxy surveys (Madau & Dickinson 2014). The solid black line shows the SFRD as measured by the CIB linear model of Maniyar et al. (2018). The solid blue line shows the corresponding constraints when the SFRD values from the galaxies are not considered as priors while performing the fit of the CIB halo model. The dotted black and red lines show the 1σ regions. 

In the text 
Fig. C.1. Measurements of the CIB power spectra from Planck Collaboration XXX (2014), Viero et al. (2013), and Lenz et al. (2019). The bestfit models shown in this figure are obtained when Planck and Herschel data are fit separately. The SPIRE data have been rescaled to Planck data using the crosscalibration factor derived in Bertincourt et al. (2016) and the factor that allows converting the measurement through the SPIRE bandpass into a measurement as it would be obtained through the Planck bandpass from Lagache et al. (2020). Overall, . 

In the text 
Fig. C.2. Measurements of the CIB power spectra from Thacker et al. (2013) and Viero et al. (2013) for a flux cut of 50 mJy. 

In the text 
Fig. C.3. Measurements of the CIB auto and crosspower spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the bestfit CIB halo model with its different components when Herschel/SPIRE data alone were fit for ℓ < 3000. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.