© A. Maniyar et al. 2021

1. Introduction

The cosmic infrared background (CIB) is made up of the cumulative emission of the infrared radiation from the dusty star-forming 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 extra-galactic sources. The CIB includes correlated anisotropies that are excellent probes of the large-scale 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 so-called thermal Sunyaev-Zel’dovich (tSZ) effect. A part of the CIB originates in the dusty star-forming 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 cross-correlation 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 large-scale 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: one-halo term (P_1h), which takes the small-scale clustering due to the correlations between the galaxies within the same halo into account; and two-halo 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 semi-empirical 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 so-called Compton parameter (y; Sect. 4). The Planck Collaboration provided an all-sky 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 cross-correlation 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¹. 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₀ = 67.47 km s⁻¹ Mpc⁻¹.

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

$$\begin{array}{c}\hfill ⟨\delta I_{\ell m}^{\nu }\delta I_{{\ell }^{\prime }{m}^{\prime }}^{{\nu }^{\prime }}⟩=C_{\ell }^{\nu \times {\nu }^{\prime }}\times {\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }},\end{array}$$(1)

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

$$\begin{array}{cc}\hfill I^{\nu }& ={\displaystyle \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}aj(\nu ,z)\mathrm{d}z\hfill \\ \hfill & ={\displaystyle \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}a\overline{j}(\nu ,z)(1+\frac{\delta j(\nu ,z)}{\overline{j}(\nu ,z)})\mathrm{d}z,\hfill \end{array}$$(2)

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

$$\begin{array}{c}\hfill C_{\ell }^{\nu \times {\nu }^{\prime }}={\displaystyle \int \frac{\mathrm{d}z}{{\chi }^2}\frac{\mathrm{d}\chi }{\mathrm{d}z}}{a}^2\overline{j}(\nu ,z)\overline{j}({\nu }^{\prime },z)P_j^{\nu \times {\nu }^{\prime }}(k=\ell /\chi ,z),\end{array}$$(3)

where at a given redshift, $P_j^{\nu \times {\nu }^{\mathrm{\prime }}}$ is the 3D power spectrum of the emissivity and is defined as

$$\begin{array}{c}\hfill ⟨\delta j(k,\nu )\delta j(k^{\prime },{\nu }^{\prime })⟩={(2\pi )}^3\overline{j}(\nu )\overline{j}({\nu }^{\prime })P_j^{\nu \times {\nu }^{\prime }}(k){\delta }^3(k-k^{\prime }).\end{array}$$(4)

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 $\frac{\mathrm{d}{j}_{\nu }}{\mathrm{d}log{M}_{\mathrm{h}}}(M_{\mathrm{h}},z)$ 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 galaxy-galaxy 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)

$$\begin{array}{c}\hfill \frac{\mathrm{SFR}}{\mathrm{BAR}}(M_{\mathrm{h}},z)=\eta ={\eta }_{\max }{e}^{-\frac{{(log{M}_{\mathrm{h}}-log{M}_{\max })}^2}{2{\sigma }_{M_{\mathrm{h}}}^2(z)}},\end{array}$$(5)

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¹² − 10¹³ 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¹² 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 free-fall 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

$$\begin{array}{c}\hfill \mathrm{BAR}(M_{\mathrm{h}},z)=⟨\dot{M}(M_{\mathrm{h}},z)⟩\times {\mathrm{\Omega }}_{\mathrm{b}}(z)/{\mathrm{\Omega }}_{\mathrm{m}}(z),\end{array}$$(6)

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

$$\begin{array}{cc}\hfill {⟨\dot{M}⟩}_{\mathrm{mean}}=& 46.1M_{\odot }{\mathrm{yr}}^{-1}{\left(\frac{M_{\mathrm{h}}}{{10}^{12}{M}_{\odot }}\right)}^{1.1}\hfill \\ \hfill & \times (1+1.11z)\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}.\hfill \end{array}$$(7)

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; Dessauges-Zavadsky 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 (semi-analytical 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 high-redshift massive galaxies (often residing in the proto-clusters, 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

$$\begin{array}{c}\hfill {\sigma }_{M_{\mathrm{h}}}(z)={\sigma }_{M_{\mathrm{h}0}}-\tau \times \max (0,z_{\mathrm{c}}-z),\end{array}$$(8)

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,

$$\begin{array}{c}\hfill \mathrm{SFR}(M_{\mathrm{h}},z)=\eta (M_{\mathrm{h}},z)\times \mathrm{BAR}(M_{\mathrm{h}},z).\end{array}$$(9)

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⁵ M_⊙. A change of the minimum mass between 10⁴ M_⊙ and 10⁸ 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

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{sub}}={\mathrm{SFR}}_{\mathrm{c}}\times \frac{m_{\mathrm{sub}}}{M_{\mathrm{h}}},\end{array}$$(10)

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⁵ M_⊙ − 10^12.9 M_⊙. We take the case of a subhalo with mass <10¹¹ 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¹¹ 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¹⁴ 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. X-ray 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¹⁴ 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 one-halo 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 = ℓ/χ),

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}=& {\displaystyle \int \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}{\left(\frac{a}{\chi }\right)}^2[\frac{\mathrm{d}{j}_{\nu ,\mathrm{c}}}{\mathrm{d}log{M}_{\mathrm{h}}}\frac{\mathrm{d}{j}_{{\nu }^{\prime },\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)\hfill \\ \hfill & +\frac{\mathrm{d}{j}_{{\nu }^{\prime },\mathrm{c}}}{\mathrm{d}log{M}_{\mathrm{h}}}\frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)\hfill \\ \hfill & +\frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}\frac{\mathrm{d}{j}_{{\nu }^{\prime },\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}{u}^2(k,M_{\mathrm{h}},z)]\hfill \\ \hfill & \times {\left(\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}\right)}^{-1}\mathrm{d}z\mathrm{d}log{M}_{\mathrm{h}},\hfill \end{array}$$(11)

where $\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}=\frac{\mathrm{d}n}{\mathrm{d}log{M}_{\mathrm{h}}}$ is the halo-mass 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 Navarro-Frenk-White (NFW) profile Navarro et al. 1997), and $\frac{\mathrm{d}j(\nu ,z)}{\mathrm{d}log{M}_{\mathrm{h}}}$ 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 one-halo power spectrum. For simplicity, we omitted the M_h and z dependence from $\frac{\mathrm{d}{j}_{\nu ,\mathrm{c}}}{\mathrm{d}log{M}_{\mathrm{h}}}$ and $\frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}$ terms from all the equations.

For the central galaxies, the differential emissivity is calculated as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{j}_{\nu ,\mathrm{c}}}{\mathrm{d}log{M}_{\mathrm{h}}}(M_{\mathrm{h}},z)=\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}\times {\chi }^2(1+z)\times \frac{{\mathrm{SFR}}_{\mathrm{dc}}}{K}\times S_{\nu }^{\mathrm{eff}}(z),\end{array}$$(12)

where $S_{\nu }^{\mathrm{eff}}(z)$ 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⁻¹⁰ M_⊙ yr⁻¹ $L_{\odot }^{-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),

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}(M_{\mathrm{h}},z)=& \frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}\times {\chi }^2(1+z)\hfill \\ \hfill & \times {\displaystyle \int \frac{\mathrm{d}N}{\mathrm{d}log{m}_{\mathrm{sub}}}}(m_{\mathrm{sub}}|M_{\mathrm{h}})\frac{{\mathrm{SFR}}_{\mathrm{sub}}}{K}\hfill \\ \hfill & \times S_{\nu }^{\mathrm{eff}}(z)\times \mathrm{d}log{m}_{\mathrm{sub}},\hfill \end{array}$$(13)

where $\frac{\mathrm{d}N}{\mathrm{d}log{m}_{\mathrm{sub}}}$ is the subhalo mass function for the satellite galaxies with a subhalo mass m_sub. The effective SEDs $S_{\nu }^{\mathrm{eff}}(z)$ 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). $S_{\nu }^{\mathrm{eff}}(z)$ 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 mid-infrared to milimeter emission of the massive star-forming galaxies up to z = 4. With this technique, we found that for the main-sequence 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 mid-IR 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)

$$\begin{array}{cc}\hfill u(k,M_{\mathrm{h}},z)=& \frac{3{\delta }_{200}}{200c^3}(cos(\mu )[\mathrm{Ci}(\mu +\mu c)-\mathrm{Ci}(\mu )]\hfill \\ \hfill & +sin(\mu )[\mathrm{Si}(\mu +\mu c)-\mathrm{Si}(\mu )]-\frac{sin(\mu c)}{\mu +\mu c}),\hfill \end{array}$$(14)

where it is to be noted that c is not the speed of light; it is the so-called halo concentration parameter defined as c(M_h, z) ≡ r₂₀₀(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₂₀₀ 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 δ₂₀₀ is a dimensionless amplitude, which can be expressed in terms of the halo concentration parameter as

$$\begin{array}{c}\hfill {\delta }_{200}=\frac{200}{3}\frac{c^3}{ln(1+c)-\frac{c}{1+c}}·\end{array}$$(15)

The two-halo term for the power spectrum of the CIB takes the clustering between galaxies in two different haloes of mass M_h and $M_{\mathrm{h}}^{\prime }$ into account and is calculated as (Béthermin et al. 2013) (where k = ℓ/χ)

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}=& {\displaystyle \int \int \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}{\left(\frac{a}{\chi }\right)}^2[\frac{\mathrm{d}{j}_{\nu ,c}}{\mathrm{d}log{M}_{\mathrm{h}}}+\frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)]\hfill \\ \hfill & \times [\frac{\mathrm{d}{j}_{{\nu }^{\prime },c}}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}+\frac{\mathrm{d}{j}_{{\nu }^{\prime },\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}u(k,M_{\mathrm{h}},z)]\hfill \\ \hfill & \times b(M_{\mathrm{h}},z)b(M_{\mathrm{h}}^{\prime },z)P_{\mathrm{lin}}(k,z)\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}log{M}_{\mathrm{h}}^{\prime }\mathrm{d}z,\hfill \end{array}$$(16)

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². 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 $b(M_{\mathrm{h}},z)b(M_{\mathrm{h}}^{\prime },z)P_{\mathrm{lin}}(k,z)$ provides the cross-power spectrum between two different haloes (M_h and $M_{\mathrm{h}}^{\prime }$) under the assumption (Cooray & Sheth 2002) $P(k,M_{\mathrm{h}},M_{\mathrm{h}}^{\prime })=b(M_{\mathrm{h}},z)b(M_{\mathrm{h}}^{\prime },z)P_{\mathrm{lin}}(k,z)$. We can simplify Eq. (16) by defining

$$\begin{array}{c}\hfill D_{\nu }(k,z)={\displaystyle \int b}(M_{\mathrm{h}},z)[\frac{\mathrm{d}{j}_{\nu ,\mathrm{c}}}{\mathrm{d}log{M}_{\mathrm{h}}}+\frac{\mathrm{d}{j}_{\nu ,\mathrm{sub}}}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)]\mathrm{d}log{M}_{\mathrm{h}},\end{array}$$(17)

which is the emissivity of the halo weighted by the corresponding bias. Therefore Eq. (16) becomes (where k = ℓ/χ)

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}={\displaystyle \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}{\left(\frac{a}{\chi }\right)}^2{D}_{\nu }(k,z)D_{\nu }^{\prime }(k,z)P_{\mathrm{lin}}(k,z)\mathrm{d}z.\end{array}$$(18)

This way of calculating the two-halo term significantly reduces the time because it reduces the number of integrals. It is also useful when the CIB−CMB lensing cross-correlation is calclutated (see Sect. 2.5). Finally, we have

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{CIB},\mathrm{clustered}}=C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}+C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}},\end{array}$$(19)

and

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{CIB},\mathrm{tot}}=C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}+C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}+C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{shot}}.\end{array}$$(20)

Here $C_{\ell ,\nu ,{\nu }^{\mathrm{\prime }}}^{\mathrm{CIB},\mathrm{clustered}}$ and $C_{\ell ,\nu ,{\nu }^{\mathrm{\prime }}}^{\mathrm{CIB},\mathrm{tot}}$ describe the clustered and total CIB anisotropy power spectrum, respectively. $C_{\ell ,\nu ,{\nu }^{\mathrm{\prime }}}^{\mathrm{shot}}$ is the shot-noise component, which we cannot predict using our halo model formalism here. The shot-noise 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 star-forming 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

$$\begin{array}{c}\hfill {\mathrm{SFRD}}_{\mathrm{c}}(z)={\displaystyle \int \frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}}\times \mathrm{SFR}(M_{\mathrm{h}},z)\times \mathrm{d}log{M}_{\mathrm{h}},\end{array}$$(21)

and for the satellite galaxies in the subhaloes, it reads

$$\begin{array}{cc}\hfill {\mathrm{SFRD}}_{\mathrm{sub}}(z)=& {\displaystyle \int \frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}}\hfill \\ \hfill & \times {\displaystyle (\int \frac{\mathrm{d}N}{\mathrm{d}log{m}_{\mathrm{sub}}}\times \mathrm{SFR}(m_{\mathrm{sub}},z)\times \mathrm{d}log{m}_{\mathrm{sub}})}\hfill \\ \hfill & \times \mathrm{d}log{M}_{\mathrm{h}},\hfill \end{array}$$(22)

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

$$\begin{array}{c}\hfill \mathrm{SFRD}(z)={\mathrm{SFRD}}_{\mathrm{c}}(z)+{\mathrm{SFRD}}_{\mathrm{sub}}(z).\end{array}$$(23)

2.5. CIB–CMB lensing cross-correlation

The large-scale 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 large-scale 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 cross-correlation as (Béthermin et al. 2013) (where k = ℓ/χ)

$$\begin{array}{c}\hfill C_{\ell }^{\nu \varphi }={\displaystyle \int \frac{\mathrm{d}\chi }{\mathrm{d}z}}\left(\frac{a}{\chi }\right)D_{\nu }(k,z)\mathrm{\Phi }(\ell ,z)P_{\mathrm{lin}}(k,z)\mathrm{d}z,\end{array}$$(24)

where Φ(ℓ, z) is given by (Challinor & Lewis 2005)

$$\begin{array}{c}\hfill \mathrm{\Phi }(\ell ,z)=\frac{3}{{\ell }^2}{\mathrm{\Omega }}_{\mathrm{m}}{\left(\frac{H_0}{c}\right)}^2\frac{\chi }{a}\left({\displaystyle \frac{{\chi }_{\star }-\chi }{{\chi }_{\star }\chi }}\right),\end{array}$$(25)

In this equation, χ_⋆ is the comoving distance to the last scattering surface, Ω_m is the matter density parameter, H₀ is the value of the Hubble parameter today, and a is the scale factor of the Universe. The measurement of this cross-correlation is not used in our likelihood to fit the CIB model parameters, but we show that our best-fit model accurately predicts this cross-correlation.

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 cross-calibration 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 colour-corrected 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

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{model}}\times {\mathit{cc}}_{\nu }\times {\mathit{cc}}_{{\nu }^{\prime }}=C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{measured}}.\end{array}$$(26)

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 $f_{\mathrm{cal}}^{\nu }$ 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 cross-calibrated 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 $f_{\mathrm{cal}}^{\nu }$ 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 best-fit 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 one-halo and the shot-noise 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 one-halo term with these data than using the Planck/HFI data alone. Because Herschel/SPIRE data can distinguish between the one-halo and the shot-noise 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 shot-noise 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 shot-noise 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 shot-noise parameters for only the auto-power spectra and used the correlation matrix to obtain the shot-noise 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 cross-power spectra as parameters instead of only taking the auto-power 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 cross-power 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 over-predicted 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 colour-corrected. 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” (Foreman-Mackey et al. 2013). We have a 14-dimensional parameter space:

• physical model parameters {η_max,M_max,σ_Mh0,τ},

• calibration factors $\{f_{600}^{\mathrm{cal}},f_{857}^{\mathrm{cal}},f_{1200}^{\mathrm{cal}}\}$,

• shot noises $\{{\mathrm{SN}}_{217}^{\mathrm{pl}},{\mathrm{SN}}_{353}^{\mathrm{pl}},{\mathrm{SN}}_{545}^{\mathrm{pl}},{\mathrm{SN}}_{857}^{\mathrm{pl}},{\mathrm{SN}}_{600}^{\mathrm{sp}},{\mathrm{SN}}_{857}^{\mathrm{sp}},{\mathrm{SN}}_{1200}^{\mathrm{sp}}\}$.

The global χ² 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 cross-power spectra from Planck/HFI and Herschel/SPIRE, respectively. For most of the cases, the one-halo term is smaller than the shot noise and can potentially be ignored on smaller scales. Based on the best-fit parameters and Eq. (24), we calculated the CIBxCMB lensing power spectra. Figure 4 shows these CIBxCMB lensing power spectra. The cross-correlation values and error bars are available for the six Planck HFI channels (100, 143, 217, 353, 545, and 857 GHz) and are provided in Planck Collaboration XVIII (2014). These values range from ℓ = 163 to ℓ = 1937, and as discussed in Planck Collaboration XVIII (2014), the non-linear 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 best-fit model agrees very well with the data points.

Fig. 2. Measurements of the CIB auto- and cross-power spectra obtained by Planck/HFI (Planck Collaboration XXX 2014) and the best-fit CIB halo model with its different components.

Fig. 3. Measurements of the CIB auto- and cross-power spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the best-fit CIB halo model with its different components.

Fig. 4. CIB × CMB lensing cross-power spectra measured by Planck Collaboration XVIII (2014) in orange and our best-fit model in blue curves.

We present the results from our fit for our model parameters in Table 1. The shot-noise 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 χ² of 113 for 80 data points for Planck power spectra. Similarly, we derive a χ² 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 χ² of 85 for Planck for 80 data points. Similarly, fitting only for the Herschel power spectra with the external priors results in the χ² of 221 for 102 data points. This shows that when we individually fit for either the Planck or Herschel data, we derive a better χ² value than when we fit for both of them together.

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 χ² for Herschel is reduced to 222, but the χ² for Planck data remains the same.

In none of the approaches does the χ² 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 χ² value (${\chi }_{\mathrm{reduced}}^2\sim 1.6$). 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 best-fit values of the one-halo, 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 best-fit 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 shot-noise level³, a factor of ∼1.9 is measured between the two. This discrepancy cannot be reconciled considering the difference in amplitude of the one-halo term. This discrepancy of the Herschel and Planck shot-noise measurements, and with the shot-noise 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 χ² 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 χ² of 221 in Table C.4 for 102 data points. However, the best-fit value for the calibration factors $f_{\mathrm{cal}}^{\nu }$ as shown in Tables C.3 and C.4 is too high. As mentioned before, the cross-calibration 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 best-fit values obtained for the $f_{\mathrm{cal}}^{\nu }$ parameters are unrealistically high and cannot be accepted. In the second approach, we keep the original priors for the Herschel $f_{\mathrm{cal}}^{\nu }$ 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 χ² (36 for 42 data points), and the shot-noise 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₁₀M_max $=12.94\begin{array}{c}+0.02\\ -0.02\end{array}{M}_{\odot }$. The highest efficiency mass found here is slightly on the high side of the range ${log}_{10}M=12.1\begin{array}{c}+0.50\\ -0.50\end{array}$ M_⊙ to $12.6\begin{array}{c}+0.10\\ -0.10\end{array}$ 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 ${log}_{10}M=12.7\begin{array}{c}+0.1\\ -0.2\end{array}$ and ${log}_{10}M=12.77\begin{array}{c}+0.128\\ -0.125\end{array}$ 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 intermediate-age 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 non-zero 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 intra-cluster 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)

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}=& f(\nu )f({\nu }^{\prime }){\displaystyle \int \frac{\mathrm{d}V}{\mathrm{d}z}}\mathrm{d}z\hfill \\ \hfill & \times {\displaystyle {\int }_{m_{\min }}^{m_{\max }}}\mathrm{d}log{M}_{\mathrm{h}}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}logV}{|y_{\ell }(M_{\mathrm{h}},z)|}^2,\hfill \end{array}$$(27)

where f(ν) gives the frequency dependence of the tSZ. It is given by

$$\begin{array}{c}\hfill f(x)=x{\displaystyle \frac{e^x+1}{e^x-1}}-4\mathrm{where}x={\displaystyle \frac{h_{\mathrm{p}}{\nu }_{\mathrm{obs}}}{k_{\mathrm{B}}{T}_{\mathrm{CMB}}}},\end{array}$$(28)

where h_p is the Planck constant, and k_B is the Boltzmann constant. y_ℓ(M_h, z) represents the two-dimensional 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)

$$\begin{array}{c}\hfill y_{\ell }(M_{\mathrm{h}},z)=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}\frac{4\pi r_{500}}{{\ell }_{500}^2}{\displaystyle {\int }_{x_{\min }}^{x_{\max }}}\mathrm{d}x{x}^2\frac{sin(\ell x/{\ell }_{500})}{\ell x/{\ell }_{500}}{P}_{\mathrm{e}}(x).\end{array}$$(29)

In this equation, c is the speed of the light, σ_T is the Thomson cross-section, m_e is the electron mass, x ≡ r/r₅₀₀, with r being the radial distance from the centre of the halo, r₅₀₀ is the radius of the sphere that contains the over-density mass M_500c of 500 times the critical density of the universe, and ℓ₅₀₀ ≡ d_A/r₅₀₀, 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

$$\begin{array}{c}\hfill P_{\mathrm{e}}(x)=C\times P_0{(c_{500}x)}^{-\gamma }{[1+{(c_{500}x)}^{\alpha }]}^{(\gamma -\beta )/\alpha },\end{array}$$(30)

where parameters (P₀, c₅₀₀, γ, α, β) were set to their best-fitting 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

$$\begin{array}{c}\hfill C=1.65{\left(\frac{h}{0.7}\right)}^2{\left(\frac{H}{H_0}\right)}^{\frac{8}{3}}{\left[\frac{(h/0.7){\stackrel{\sim }{M}}_{500c}}{3\times {10}^{14}{M}_{\odot }}\right]}^{\frac{2}{3}+0.12}\mathrm{eV}{\mathrm{cm}}^{-3},\end{array}$$(31)

where H is the Hubble constant at redshift z, with H₀ being its local value and h the reduced Hubble constant (h = H₀/100). The mass used here M̃_500c is not necessarily the true mass, but can contain a bias due to observational effects and non-thermal 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⁻⁶ and x_max = 10 in Eq. (29).

Because massive clusters are rare, the contribution of the two-halo term to the total tSZ angular power spectrum is far lower than that of the one-halo 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 two-halo term plays a significant role in the two-halo term of the CIB × tSZ cross-power spectra and cannot be ignored. In our analysis, we therefore do not neglect this component.

It is straightforward to calculate the two-halo component of the power spectrum as (e.g. Salvati et al. 2018)

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}=& f(\nu )f({\nu }^{\prime }){\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z{P}_{\mathrm{lin}}(k=\ell /\chi ,z)\hfill \\ \hfill & \times {\{{\displaystyle {\int }_{m_{\min }}^{m_{\max }}}\mathrm{d}log{M}_{\mathrm{h}}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}logV}b(M_{\mathrm{h}},z)y_{\ell }(M_{\mathrm{h}},z)\}}^2.\hfill \end{array}$$(32)

The total tSZ power spectrum is the sum of the one-halo and two-halo components:

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{tSZ}-\mathrm{tot}}=C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}+C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}.\end{array}$$(33)

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₀, c₅₀₀, γ, α, β) are not varied, this allows us to tabulate y_ℓ, which speeds up the process considerably. Figure 6 shows the prediction for the tSZ auto-power 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 two-halo term contributes about 2% to the total power spectrum and the one-halo 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 best-fit model is consistent with these measurements. At ℓ = 257.5, using the best-fit parameters mentioned before, we obtain $\ell (\ell +1)C_{\ell }^{\mathrm{tSZ}}/2\pi =0.178$. For this ℓ, Planck Collaboration XXII (2016) provided the measured y power spectrum (10¹²y²) value as 0.217 ± 0.049, which is consistent with our predictions (the best-fit 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 best-fit 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 cross-correlation

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 cross-correlation 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 one-halo term provides the correlation between the CIB sources and the tSZ within the same halo. On the other hand, the two-halo 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 one-halo term would be zero, but the two-halo 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 two-halo 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 two-halo terms are given as (for k = ℓ/χ)

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}=& {\displaystyle \int \mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z}\int \mathrm{d}log{M}_{\mathrm{h}}}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}logV}{y}_{\ell }\hfill \\ \hfill & \times [\{\frac{d{j}_{{\nu }_1,c}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}}+\frac{\mathrm{d}{j}_{{\nu }_1,\mathrm{sub}}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)\}f({\nu }_2)\hfill \\ \hfill & +\{\frac{d{j}_{{\nu }_2,c}^{\prime }}{dlog{M}_{\mathrm{h}}}+\frac{\mathrm{d}{j}_{{\nu }_2,\mathrm{sub}}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}}u(k,M_{\mathrm{h}},z)\}f({\nu }_1)],\hfill \end{array}$$(34)

and

$$\begin{array}{cc}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}=& {\displaystyle \int \mathrm{d}z}\frac{\mathrm{d}V}{\mathrm{d}z}{P}_{\mathrm{lin}}(k,z)\hfill \\ \hfill & \times {\displaystyle \int \mathrm{d}}log{M}_{\mathrm{h}}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}logV}{y}_{\ell }b(M_{\mathrm{h}},z)\hfill \\ \hfill & \times {\displaystyle \int \mathrm{d}log{M}_{\mathrm{h}}^{\prime }}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }\mathrm{d}logV}b(M_{\mathrm{h}}^{\prime },z)\hfill \\ \hfill & \times [\{\frac{\mathrm{d}{j}_{{\nu }_1,c}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}+\frac{\mathrm{d}{j}_{{\nu }_1,\mathrm{sub}}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}u(k,M_{\mathrm{h}}^{\prime },z)\}f({\nu }_2)\hfill \\ \hfill & +\{\frac{\mathrm{d}{j}_{{\nu }_2,c}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}+\frac{\mathrm{d}{j}_{{\nu }_2,\mathrm{sub}}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}^{\prime }}u(k,M_{\mathrm{h}}^{\prime },z)\}f({\nu }_1)],\hfill \end{array}$$(35)

where

$$\begin{array}{c}\hfill \frac{\mathrm{d}{j}_{\nu ,c(\mathrm{sub})}^{\prime }}{\mathrm{d}log{M}_{\mathrm{h}}}(M_{\mathrm{h}},z)=\frac{\mathrm{d}{j}_{\nu ,c(\mathrm{sub})}}{\mathrm{d}log{M}_{\mathrm{h}}}(M_{\mathrm{h}},z)\frac{a}{{\chi }^2}{\left(\frac{{\mathrm{d}}^{2}N}{\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}V}\right)}^{-1}.\end{array}$$(36)

M_h and $M_{\mathrm{h}}^{\prime }$ represent the mass of the two different haloes between which the correlation is computed when the two-halo term is computed. The tSZ spectral variation does not dependent on the redshift or halo mass in the non-relativistic limit, therefore the cross spectra are simpler than they would be otherwise.

The total CIB × tSZ power is then

$$\begin{array}{c}\hfill C_{\ell ,\nu ,{\nu }^{\prime }}^{\mathrm{CIB}\times \mathrm{tSZ}-\mathrm{tot}}=C_{\ell ,\nu ,{\nu }^{\prime }}^{1\mathrm{h}}+C_{\ell ,\nu ,{\nu }^{\prime }}^{2\mathrm{h}}.\end{array}$$(37)

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 template-based approach for the CIB−tSZ correlation was used by introducing a parameter ξ to act as a correlation coefficient. This template-based 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 non-negligible 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₅₀₀, 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₂₀₀ 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₂₀₀, then converting the M₂₀₀ into corresponding M₅₀₀ for the given redshift range and cosmology using the procedure given in Komatsu & Seljak (2001), and then using this M₅₀₀ to calculate the corresponding tSZ terms. We presented the CIB × tSZ results here, we worked with a single definition of M₅₀₀ for the halo masses for the CIB and tSZ terms. We verified the effect of using M₅₀₀ instead of M₂₀₀ 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₅₀₀ instead of M₂₀₀: 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₂₀₀ 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₅₀₀ 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 cross-correlation 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 cross-power spectra. The tSZ term dominates at lower frequency channels and drives the cross-correlation. 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 one-halo term of the power spectrum to the total power spectrum to relatively decrease compared to the two-halo term at higher redshifts. This is indeed what we observe in Fig. 8, where the multipole where the one-halo and two-halo term contribute equally moves slightly to higher values in higher redshift bins, that is, the range in which the two-halo term dominates the one-halo 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 two-halo 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 cross-correlation, 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 power-spectrum 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 scale-independent 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 cross-correlation. 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 best-fit 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 signal-to-noise 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

$$\begin{array}{c}\hfill {\xi }_{\ell }(\nu )=\frac{C_{\ell ,\nu }^{\mathrm{CIB}\times \mathrm{tSZ}}}{\sqrt{C_{\ell ,\nu }^{\mathrm{CIB}}\times C_{\ell ,\nu }^{\mathrm{tSZ}}}},\end{array}$$(38)

where the correlation coefficient ξ_ℓ(ν) depends on frequency and angular scale. The CIB power spectra $C_{\ell ,\nu }^{\mathrm{CIB}}$ 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 best-fit 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² sr⁻¹. This was added to the one- and two-halo 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 cross-correlation 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 cross-correlation 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 x-axis 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 template-based 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 cross-correlation 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 $0.20\begin{array}{c}+0.14\\ -0.11\end{array}$ and $0.16\begin{array}{c}+0.10\\ -0.10\end{array}$, respectively, at ℓ = 3000. Our best-fit 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 dark-matter haloes for the highest efficiency is ${log}_{10}{M}_{\max }=12.94\begin{array}{c}+0.02\\ -0.02\end{array}$ M_⊙, while the maximum efficiency at this mass is $\eta =0.42\begin{array}{c}+0.03\\ -0.02\end{array}$. The mass of the highest efficiency found here is slightly on the high side of the range ${log}_{10}M=12.1\begin{array}{c}+0.50\\ -0.50\end{array}$ M_⊙ to $12.6\begin{array}{c}+0.10\\ -0.10\end{array}$ 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 ${log}_{10}M=12.7\begin{array}{c}+0.1\\ -0.2\end{array}$ and ${log}_{10}M=12.77\begin{array}{c}+0.128\\ -0.125\end{array}$ 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 cross-correlation using the best-fit value of the CIB parameters and find that it agrees well with the measurements (Fig. 4). While we obtain a decent χ² of 113 for 80 Planck data points, the χ² 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 two-halo 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 two-halo 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 two-halo term has a significant contribution to the two-halo 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 cross-correlation 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 two-halo term for the CIB−tSZ correlation is also not negligible with respect to the one-halo term and should be considered in the total power spectrum calculation. This two-halo 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 template-based 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.

Acknowledgments

We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France and support from the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government 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 Aix-Marseille University-A*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

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 best-fit 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 best-fit parameters of the CIB model are defined. In this case, the χ² value obtained is 90 for 80 Planck data points and 156 for 102 Herschel data points. However, the best-fit value for the calibration factor for 1200 GHz $f_{1200}^{\mathrm{cal}}$, which is negatively correlated with the shot-noise 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 χ² 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 χ² 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 power-law parametrisation. For our case, we used a double power law to describe the relation between η and M_h, which reads

$$\begin{array}{c}\hfill \frac{\mathrm{SFR}}{\mathrm{BAR}}=\eta ={\eta }_{\max }{[{\left(\frac{M_{\mathrm{h}}}{M_{\max }}\right)}^{-\beta }+{\left(\frac{M_{\mathrm{h}}}{M_{\max }}\right)}^{\gamma }]}^{-1},\end{array}$$(B.1)

where β, and γ describe the η at the low- and high-mass 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 low-mass 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 best-fit values for the halo model obtained under different conditions. The shot-noise 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 auto-power 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 shot-noise 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 best-fit 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 cross-calibration 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, $C_{\ell }^{Planck,857\mathrm{GHz}}=1.008\times C_{\ell }^{\mathrm{SPIRE},857\mathrm{GHz}}$.

In Fig. C.2 we show the comparison of the measured CIB auto-power 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 best-fit 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χ² 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.

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 χ² value in this case. The best-fit 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.

In order to verify the effect of the cross-calibration 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 cross-calibration 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 χ² 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 cross-calibration measurement between Planck and Herschel (Bertincourt et al. 2016). This test therefore does not solve the compatibility problem of the two data sets.

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 χ² value, and the shot-noise 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.

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

All Tables

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 cross-power spectra obtained by Planck/HFI (Planck Collaboration XXX 2014) and the best-fit CIB halo model with its different components.
In the text

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

	Fig. 4. CIB × CMB lensing cross-power spectra measured by Planck Collaboration XVIII (2014) in orange and our best-fit 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 cross-correlation 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 x-axis 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 best-fit 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 cross-calibration 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, $C_{\ell }^{Planck,857\mathrm{GHz}}=1.008\times C_{\ell }^{\mathrm{SPIRE},857\mathrm{GHz}}$.

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 cross-power spectra obtained by Herschel/SPIRE (Viero et al. 2013) and the best-fit CIB halo model with its different components when Herschel/SPIRE data alone were fit for ℓ <  3000.
In the text