Simple halo model formalism for the cosmic infrared background and its correlation with the thermal Sunyaev Zel'dovich effect

Modelling the anisotropies in the CIB on all the scales is a challenging task due to the complex nature of the galaxy evolution and thus often requires too many parameters in order to fit the observational data. In this paper, we present a new halo model for the anisotropies of the CIB using only four parameters. Our model connects the mass accretion on the dark matter halos to the star formation rate. Using this model, we find that the halo mass with the maximum efficiency for converting the accreted baryons into stars is $\log_{10}M_\mathrm{max} = {12.94}^{+0.02}_{-0.02} \: M_\odot$, consistent with other studies. Accounting for the mass loss through stellar evolution, we find, for an intermediate age galaxy, that the star formation efficiency defined as $M_\star(z)/M_b(z)$ is equal to 0.19 and 0.21 at redshift 0.1 and 2 respectively, in good agreement with the values obtained by previous studies. It is the first time that a CIB model is used to fit Planck and Herschel CIB power spectra simultaneously. However we find that the large angular scale Planck and Herschel data are not fully compatible with the small scale Herschel data (for $\ell>3000$). The CIB is expected to be correlated with the tSZ signal of the galaxy clusters. Using this halo model for the CIB, and a halo model for the tSZ with a single parameter, we also provide a consistent framework to calculate this cross-correlation, which requires no extra parameter. The CIB$\times$tSZ correlation is found to be higher when inferred with a combination of two widely spaced frequency channels (e.g. 143x857 GHz). The CIB, tSZ, and CIB$\times$tSZ act as foregrounds while measuring the kSZ power spectrum from the CMB power spectrum and need to be removed. Due to its simplistic nature and low number of parameters, the halo model formalism presented here is quite useful for such an analysis to measure the kSZ power spectrum accurately.


Introduction
The cosmic infrared background (CIB) is made up of the cumulative emission of the infrared radiation from the dusty starforming galaxies throughout the Universe. It traces the star formation history of the Universe, which spans a wide range of redshift 0 ≤ z ∼ 6. Measurements of the CIB can thus be used as a powerful tool to map the star formation at high redshifts (Knox et al., 2001). Although the CIB was first detected by Puget et al. (1996),  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 et al., 2014c, Viero et al., 2013a. Another such tracer of the underlying dark matter distribution are massive galaxy clusters. Hot electrons in these galaxy clus-ters Compton scatter the CMB photons and give rise to the socalled 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 by Reichardt et al. (2012), Dunkley et al. (2013), George et al. (2015), Reichardt et al. (2020), and{Choi et al. (2020) and thus has to be considered in the CMB power spectrum data analysis. Planck Collaboration et al. (2016c) also reported a measurement of the crosscorrelation between the tSZ and CIB using Planck data. In order to model the CIB×tSZ signal, we in turn need accurate models of the CIB anisotropies and the tSZ.
On large angular scales, we can use the fact that the clustering of the CIB traces the 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 et al. (2014c) and Maniyar et al. (2018). However, in order to describe the Article number, page 1 of 21 arXiv:2006.16329v2 [astro-ph.CO] 17 Nov 2020 A&A proofs: manuscript no. halo_model_cib anisotropies on both the large and small angular scales coherently, a 'halo model' approach as developed by Cooray & Sheth (2002) is generally used. With the assumption inside the halo model that all the galaxies reside in the dark matter haloes, the clustering can be considered as the sum of two components: onehalo term (P 1h ), which takes the 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 et al., 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. 2013b).
Subsequently, several studies such as Shang et al. (2012), Viero et al. (2013b), andPlanck Collaboration et al. (2014c) 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 et al., 2014c). 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; Sec. 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 et al., 2014a, 2016d. Bolliet et al. (2018) used these data to constrain the cosmological parameters (equation of state of the dark energy w in particular) along with the tSZ parameter. We use this halo model of the tSZ to calculate the tSZ power spectra. Addison et al. (2012) calculated the CIB-tSZ correlation within the halo model framework. They first presented a formalism to calculate this correlation using a CIB halo model that does not account for the dependence of the source flux on the halo mass, and then expanded their formalism to account for this effect. However, the CIB halo model they finally considered (from Xia et al. (2012)) to calculate the CIB-tSZ correlation does not account for the aforementioned L − M h dependency. We present a halo model formalism to calculate the CIB×tSZ crosscorrelation using our new halo model for the CIB and the halo model for the tSZ from Bolliet et al. (2018).
This paper is structured as follows. We begin in Sect. 2 by presenting our halo model for the CIB anisotropies and subsequently present the CIBxCMB lensing correlation within this framework. Sect. 3 then presents the constraints on the CIB model parameters through the data and corresponding results. In Sect. 4 we provide the halo model for the tSZ. Finally, in Sect. 5 we present the halo model formalism for the CIB×tSZ correlation and the predictions for its power spectra and angular scale dependence 1 . Conclusions are given in Sect. 6. Appendices A, B, and 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 et al., 2016b) with Ω m = 0.33 and H 0 = 67.47 km s −1 Mpc −1 .

New halo model for the CIB power spectrum
The starting point for our model is the accretion of matter onto the dark matter haloes. We then connect the accretion of the baryonic gas onto the dark matter haloes to the SFR corresponding to these haloes. The SFR is defined separately for the central and satellite galaxies in a halo. Using the SFR from a given halo, we then calculate the emissivity of all the haloes of mass M h at a given redshift, which then is used to calculate the angular power spectrum of the CIB anisotropies.
The angular power spectrum of the CIB anisotropies is defined as where ν is the frequency of the observation and I ν is the intensity of the CIB measured at that frequency. The intensity is a function of the comoving emissivity j through where χ(z) is the comoving distance to redshift z, and a = 1/(1 + z) is the scale factor of the Universe, and δ j(ν, z) are the emissivity fluctuations of the CIB. We expand the Eq. 2 this way between the mean value and its fluctuations because Eq. 1 shows that the power spectrum is calculated using the fluctuations around the mean value. Combining Eqs. 1 and 2, and using the Limber approximation (Limber, 1954) for the flat sky, which helps us avoid the spherical Bessel function calculations and makes the computation easier, we therefore obtain where at a given redshift, P ν×ν j is the 3D power spectrum of the emissivity and is defined as Thus we have to calculate δ j to obtain the CIB angular power spectrum. For this purpose we connect (Sec. 2.3) the SFR of the haloes with the specific emissivity d j ν d log M h (M h , z) and integrate it over all the halo masses and redshift range to obtain the CIB power spectra.

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) where M h is the halo mass, M max represents the mass with the highest star formation efficiency η max , and σ M h (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., 2013a, Planck Collaboration et al., 2014c, Maniyar et al., 2018 have found that the dark matter haloes within the mass range 10 12 − 10 13 M form the stars most efficiently. With an empirical model, Béthermin et al. (2013) showed in their Fig. 17 that the star formation efficiency as a function of instantaneous halo mass is highest in haloes with masses ∼ 10 12 M . This mass does not change considerably over a range of redshifts, whereas the efficiency falls off drastically for masses above and below the most efficient mass. This effect can be understood physically. Below the mass for which star formation is most efficient, the gravitational potential of the dark matter halo is lower and the supernovae feedback is strong enough to remove the gas from the galaxy (e.g. Silk, 2003) and thereby decreases the star formation. On the other side of the spectrum, at higher masses, the cooling time of the gas becomes much longer than the 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 σ M h , and the SFR falls off considerably on very high and very low masses.
For a given halo mass at a given redshift, the BAR is given as where Ω b and Ω m are the dimensionless cosmological baryonic density and total matter density, respectively (thus, the ratio of the two gives the baryonic matter fraction at a given redshift, which is in fact constant with redshift because they have the same evolution with redshift).Ṁ(M h , z) is the mass growth rate. Fakhouri et al. (2010) provided the mean and median mass growth rates of haloes of mass M h at redshift z. We used the mean mass growth rate, given as This approach assumes that there is no 'gas reservoir effect', that is, the accreted gas is immediately converted into stars and is not collected over time to form reservoirs in (Daddi et al., 2010) and around galaxies (Cantalupo et al., 2012). It has been shown (e.g. Saintonge et al., 2013, Béthermin et al., 2015, 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. (2015Cousin et al. ( , 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 σ M h evolve with redshift. The motivation for letting σ M h 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 σ M h evolve with redshift as where z c is a redshift below which σ M h evolves with redshift, σ M h0 is the value of σ M h 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 σ M h0 . 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.

SFR for the haloes and subhaloes
For a given value of the halo mass and redshift, we can calculate η using Eq. 5 and multiply it by the corresponding BAR calculated using Eq. 6 to obtain the SFR, that is, This is the procedure with which the SFR can be obtained for the haloes. To calculate the SFR for the subhaloes residing within these haloes, the procedure is slightly modified. We first assumed that for a given halo with mass M h , the subhalo masses (m sub ) range from M min to M h . In this analysis, we fixed M min = 10 5 M . A change of the minimum mass between 10 4 M and 10 8 M changes the calculation of the power spectra only negligibly. The SFR for the subhaloes can be estimated in two ways. The first way is an approach similar to the one for the haloes, which is calculating the efficiency η and then multiplying with the BAR value to obtain the SFR, that is, replacing M h by m sub in Eq. 9. This assumes the same lognormal parametrisation of η for subhaloes as of the central haloes. The other way to estimate the SFR in subhaloes is that is, the SFR for the subhalo is obtained by weighing the halo SFR (SFR c ) by the ratio of subhalo to halo mass. For every subhalo of a given halo, we estimated the SFR with both these approaches and took the smaller of the two as representative of the SFR for the subhalo.
The reasoning for this is explained in Fig. 1. We first consider case 1 in the figure. In this case, the main halo has a mass (∼ 10 12.9 M ) very near to the efficiency peak of star formation (for this particular choice of parameters for the lognormal), that is, the central galaxy forms stars very efficiently (see Eq. 9). This 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 (Sec. 2.2). halo has subhaloes ranging from 10 5 M − 10 12.9 M . We take the case of a subhalo with mass < 10 11 M . As was pointed out before, subhaloes with very low mass have a low gravitational potential, and it is hard for them to hold on to the gas inside against the pressure from supernova feedback, for instance. Thus, they are expected to have low star formation. This is satisfied in case 1 because at lower masses, the efficiency is indeed very low and is not expected contribute significantly to the total SFR of the halo. In this case, the SFR for the subhaloes can therefore be directly estimated by substituting the halo mass (M h ) by subhalo mass (m sub ) in Eq. 9. As an example, for a subhalo of mass 10 11 M belonging to a central halo of mass 10 12.9 M , the SFR calculated using Eq. 9 is 2% of the value that we derive for the SFR obtained using Eq. 10, and we therefore take the former as the SFR value.
In the second case, the main halo is quite massive (> 10 14 M ) and far away from the efficiency peak. It therefore does not have strong star formation. However, in this case, this halo can contain a subhalo with a mass of about the efficiency peak (∼ 10 12.9 M ). According to Eq. 9, when we substitute M h by m sub , this subhalo will have a significant amount of star formation. Again, as pointed out earlier, the gas inside massive haloes is quite hot, and there are mechanisms at play (e.g. Xray heating and AGN feedback) that suppress the gas cooling and hence make it difficult to form stars. Moreover, the central massive galaxies in these haloes can strip the gas out of the satellite galaxies in the subhaloes and thus decrease star formation. These subhaloes are therefore not expected to contribute significantly to the SFR, in contrast to what would be obtained with Eq. 9. This inherently assumes instantaneous quenching, that is, the satellite galaxies are quenched at the same time as the central galaxies in a given parent halo. If we were to explicitly avoid this assumption, we would need to introduce an additional quenching parameter as a function of subhalo mass and redshift. This would result in one or two additional parameters for the model. How-ever, as we mentioned earlier, our main purpose here is to build a very simple halo model of the CIB with as few parameters as possible. One way to correct for this in these cases therefore is to weight the SFR of the main halo by the mass fraction of the corresponding subhalo, that is, use Eq. 10 to obtain the SFR of the subhalo. The SFR in this case would be lower than the rate obtained using Eq. 9 (substituting M h by m sub , of course). Again as an example, for a subhalo of mass 10 12.9 M belonging to a central halo of mass 10 14 M , the SFR calculated using Eq. 9 is twice higher than the value we obtain for the SFR with Eq. 10 at z = 2, and we therefore take the latter as the SFR value.
Although Eq. 9 and Eq. 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.

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) 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 d j(ν,z) d log M 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 For the central galaxies, the differential emissivity is calculated as where S e f f ν (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 −10 M yr −1 L −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), (13) where dN d log m sub is the subhalo mass function for the satellite galaxies with a subhalo mass m sub . The effective SEDs S e f f ν (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 e f f ν (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) and Béthermin et al. (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) where it is to be noted that c is not the speed of light; it is the so-called halo concentration parameter defined as being a characteristic radius (see e.g. Navarro et al., 1997) and r 200 being the radius containing the mass giving 200 times the critical density of the Universe at a redshift (some studies use the average density instead of the critical density; we used the critical density), Ci(x) and Si(x) are standard cosine and sine integrals, respectively, µ ≡ kr , where δ 200 is a dimensionless amplitude, which can be expressed in terms of the halo concentration parameter as The 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 h into account and is calculated as (Béthermin et al., 2013) Article number, page 5 of 21 A&A proofs: manuscript no. halo_model_cib is the halo bias prescription given by , P lin (k, z) is the linear matter power spectrum, which we calculated using CAMB 2 . We used the u(k, M h , z) term only for the subhaloes. This means that in our analysis, the central galaxy is assumed to be at the centre of the halo and the satellite galaxies in the subhaloes are distributed according to the NFW profile. The clustering term given by b(M h , z)b(M h , z)P lin (k, z) provides the cross-power spectrum between two different haloes (M h and M h ) under the assumption (Cooray & Sheth, 2002) . We can simplify Eq. 16 by defining which is the emissivity of the halo weighted by the corresponding bias. Therefore Eq. 16 becomes (where k = /χ) This way of calculating the 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 Sec. 2.5). Finally, we have and Here C CIB,clustered ,ν,ν and C CIB,tot ,ν,ν describe the clustered and total CIB anisotropy power spectrum, respectively. C 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 Sec. 3.1.

ρ 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 and for the satellite galaxies in the subhaloes, it reads

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 et al., 2014b). Within the context of our halo model, we can calculate this CIB-CMB lensing cross-correlation as (Béthermin et al., 2013) where Φ( , z) is given by (Challinor & Lewis, 2005) In this equation, χ is the comoving distance to the last scattering surface, Ω m is the matter density parameter, H 0 is the value of the Hubble parameter today, and a is the scale factor of the Universe. The measurement of this 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. We also used the CIB power spectra measurement from Herschel/SPIRE data at 600, 857, and 1200 GHz provided by Viero et al. (2013a). Bertincourt et al. (2016) calculated the crosscalibration factors between the power spectra measured by Viero et al. (2013a) at 600 and 857 GHz and the power spectra measured by Planck Collaboration et al. (2014c) 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%.

Constraints on the model through data
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 et al., 2014c). 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 . The CIB power spectra were then corrected as The CIB power spectra error bars do not account for the absolute calibration uncertainties. To account for these uncertainties, Béthermin et al. (2011) introduced a calibration factor f ν cal 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 ν cal 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 Tab. 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 et al. (2014c) 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), andwas 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 ∼ 30360). 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 et al. (2014c), 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 autoand 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.

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 overpredicted the SFRD of the Universe compared to the SFRD measured by external groups using the infrared luminosity functions. We would like the halo model to be able to fit the CIB power spectra and also produce the correct SFRD history. Thus, we used the ρ SFR measurements at different redshifts that were obtained by measuring the infrared luminosity functions from Gruppioni et al. (2013), Magnelli et al. (2013), Marchetti et al. (2016) as priors while performing the fit. In order to account for the different sets of cosmological parameters used by them, we converted the SFRD values into the observed flux between 8 − 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  The global χ 2 has a contribution from the CIBxCIB, priors on calibration factors, and priors imposed by the external observational constraints mentioned above. We assumed Gaussian uncorrelated error bars for measurement uncertainties, which simplifies the covariance matrix that ultimately contains only diagonal terms. 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 bestfit 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)  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.

Results
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 χ 2 of 113 for 80 data points for Planck power spectra. Similarly, we derive a χ 2 of 247 for Herschel power spectra for 102 data points. When we only fit for the Planck data without the Herschel power spectrum data but with the external priors, we derive a χ 2 of 85 for Planck for 80 data points. Similarly, fitting only for the Herschel power spectra with the external priors results in the χ 2 of 221 for 102 data points. This shows that when we individually fit for either the Planck or Herschel data, we derive a better χ 2 value than when we fit for both of them together.
We also performed this analysis with the updated values of the mean level of the CIB at 350, 350, and 500 µm Herschel/SPIRE bands used as priors in our work (Duivenvoorden et al., 2020). With the updated priors, the χ 2 for Herschel is reduced to 222, but the χ 2 for Planck data remains the same.
In none of the approaches does the χ 2 value for the Herschel power spectra appear to fare as well as the Planck value. Similarly, using a different halo model with a larger number of parameters, Viero et al. (2013a) did not obtain a good χ 2 value (χ 2 reduced ∼ 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 et al., 2014c), Herschel/SPIRE (Viero et al., 2013a), 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 A1 of Lagache et al. (2019). 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 3 , 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. (2019). Inconsistencies also exist inside a single experiment. For example, we show in We explored two ways of reconciling the Planck and Herschel measurements using our halo model. In the first approach, we broadened the error bars on the calibration factors for the 600, 857, and 1200 GHz Herschel channels to 0.05 each with their central values at 1.00 instead of using the values obtained for the relative calibration between Planck and Herschel by Bertincourt et al. (2016), as done in our original approach. In this setting, we indeed find a very good fit to the Herschel power spectra, with a χ 2 of 127 when fit together with Planck data (Tab. C.3), and 137 when the Herschel data are fit alone (Tab C.2) compared to the χ 2 of 221 in Tab C.4 for 102 data points. However, the bestfit value for the calibration factors f ν cal as shown in Tab. 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 Sec.3.1). Thus the best-fit values obtained for the f ν cal parameters are unrealistically high and cannot be accepted. In the second approach, we keep the original priors for the Herschel f ν cal parameters from Bertincourt et al. (2016) and fit the data for low multipoles only, that is, < 3000. Results are provided in Tab. C.5 and Fig. C.3, which show that we obtain a good χ 2 (36 for 42 data points), and the shot-noise levels for Planck and Herschel are quite similar (see Tab. C.1 and C.5), which is what we expect. This shows that the Herschel data at > 3000 are not compatible with the values expected from the < 3000 Planck and Herschel data. The mass of the dark matter haloes for converting the accreted baryons into stars most efficiently is log 10 M max = 12.94 +0.02 −0.02 M . The highest efficiency mass found here is slightly on the high side of the range log 10 M = 12.1 +0.50 −0.50 M to 12.6 +0.10  (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 , σ M h0 , 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-40% across all redshifts, which is consistent with our result. Moreover, Moster et al. (2018) calculated the star formation efficiency as 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,   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 et al. (2014c) 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 Fig. 3. Measurements of the CIB auto-and cross-power spectra obtained by Herschel/SPIRE (Viero et al., 2013a) and the best-fit CIB halo model with its different components. 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.

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) where f (ν) gives the frequency dependence of the tSZ. It is given by where h p is the Planck constant, and k B is the Boltzmann constant. y (M h , z) represents the 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) y (M h , z) = σ T m e c 2 4πr 500 2 500 x max x min dxx 2 sin( x/ 500 ) x/ 500 P e (x) . (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 500 , with r being the radial distance from the centre of the halo, r 500 is the radius of the sphere that contains the over-density mass M 500c of 500 times the critical density of the universe, and 500 ≡ d A /r 500 , with d A being the angular diameter distance. As was done for the CIB modelling, we used the NFW profile to represent the density distribution inside the dark matter halo. Using this distribution, we have where parameters (P 0 , c 500 , γ, α, β) were set to their best-fitting values obtained by Planck Collaboration et al. (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   5. Measurements of the infrared SFRD using galaxies (Madau & Dickinson, 2014 where H is the Hubble constant at redshift z, with H 0 being its local value and h the reduced Hubble constant (h = H 0 /100).
The mass used hereM 500 c 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 500 c toM 500 c as M 500 c = M 500 c /B. In our analysis, we took the mass and redshift range to be the same as was used for the CIB model with x min = 10 −6 and x max = 10 in Eq. 29. Because massive clusters are rare, the contribution of the 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 Sec. 5.3, the two-halo term plays a significant role in the twohalo 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) The total tSZ power spectrum is the sum of the one-halo and two-halo components: f (ν) is negative for ν ≤ 217 GHz. This means that the tSZ power spectrum for certain choices of frequency combinations will be negative. The case for the CIB-tSZ correlation described in Sec. 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 Tab.1 of Planck Collaboration et al. (2016d).
Because the values of the parameters (P 0 , c 500 , γ, α, β) are not varied, this allows us to tabulate y , which speeds up the process considerably. Figure 6 shows the prediction for the tSZ 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 et al. (2016d). 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 Article number, page 11 of 21 A&A proofs: manuscript no. halo_model_cib  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 ( + 1) C tSZ /2π = 0.178. For this , Planck Collaboration et al. (2016d) provided the measured y power spectrum (10 12 y 2 ) 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.

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.

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 and where M h and M h 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 CIBxtSZ power is then When the CIB and tSZ model parameters are knonw, 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 nonnegligible variation in ξ with the angular scales and the considered frequencies.

Halo mass definition
Our models for CIB and tSZ (and thus CIB×tSZ) rely on the halo mass M h . For the tSZ effect, the pressure profile from Eq. 30 is in general given in terms of M 500 , which is defined as the mass contained within a radius in which the mean overdensity is 500 times the background critical density. In contrast, the CIB studies are generally carried out with M 200 as the definition of the halo mass. We have to be consistent with the definition of the halo mass while considering the CIB and tSZ together. One way of doing this is calculating the CIB terms with M 200 , then converting the M 200 into corresponding M 500 for the given redshift range and cosmology using the procedure given in Komatsu & Seljak (2001), and then using this M 500 to calculate the corresponding tSZ terms. We presented the CIB×tSZ results here, we worked with a single definition of M 500 for the halo masses for the CIB and tSZ terms. We verified the effect of using M 500 instead of M 200 on the CIB model parameters, and we found that the values of the CIB parameters slightly change. The most notable change was observed in the value of η max , which was found to be higher when we used M 500 instead of M 200 : because with this definition of the mass, the CIB contribution is calculated very close to the centre of the cluster and does not include the star formation that occurs far away. The model accordingly tries to increase the conversion efficiency of the accreted baryons into stars to produce the same level of CIB as observed. The CIB model parameters given in Tab. 1 were calculated using the commonly used M 200 definition for the halo masses. Figure 7 shows the predictions for the CIB, tSZ, and CIBxtSZ power spectra correlating the Planck 143 GHz frequency with other frequencies based on the models of the CIB, tSZ, and CIBxtSZ given above. We note again that our results were derived using the M 500 definition for the halo masses, and for the tSZ, the f (ν) factor was convolved with the Planck/HFI frequencies. This provided values of -4.03, -2. 79,0.19,6.21,14.46,and 26.34 for 100,143,217,343,545,and 857 GHz,respectively,as mentioned in Sec. 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, CIBxtSZ power spectra shown here are negative in all the cases from 100 to 857 GHz. A special case to consider when Eqs. 34 and 35 are used is that of the two widely separated frequency channels. Fig. 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 143x217 GHz are far smaller than for 143x143 and 143x353 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 CIBxtSZ 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. 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.

Redshift contributions to the power spectra
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. Fig. 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.

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.
A&A proofs: manuscript no. halo_model_cib The kSZ has been constrained by Dunkley et al. (2013), George et al. (2015), Reichardt et al. (2020), andChoi et al. (2020) using the CMB power-spectrum measurements from the Atacama Cosmology Telescope (ACT) or the South Pole Telescope (SPT), while Planck Collaboration et al. (2016c) 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) andPlanck Collaboration et al. (2016c) 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 et al. (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 et al. (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-tonoise ratio on the small angular scales of interest for the kSZ power spectrum.

Correlation coefficient for CIB×tSZ
We define the correlation coefficient for the CIB and tSZ correlation as where the correlation coefficient ξ (ν) depends on frequency and angular scale. The CIB power spectra C 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 2 /sr. This was added to the oneand two-halo terms at 143 GHz. The CIB×tSZ correlation at 143 GHz is negative, so that we took its absolute value. Based on Article number, page 15 of 21 A&A proofs: manuscript no. halo_model_cib 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. 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 Sec. 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 ob-tained 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 +0.14 −0.11 and 0.16 +0.10 −0.10 , 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.

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 +0.02 −0.02 M , while the maximum efficiency at this mass is η = 0.42 +0.03 −0.02 . The mass of the highest efficiency found here is slightly on the high side of the range log 10 M = 12.1 +0.50 −0.50 M to 12.6 +0.10 −0.10 M found by Viero et al. (2013a) and Planck Collaboration et al. (2014c), but it agrees well with Chen et al. (2016) for faint SMGs of log 10 M = 12.7 +0.1 −0.2 and log 10 M = 12.77 +0.128 −0.125 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 χ 2 of 113 for 80 Planck data points, the χ 2 for Herschel power spectra comes is 247 for 102 data points, which is not good. We investigated this discrepancy using different methods and found that the CIB power spectrum measurements are not fully compatible with each other for >3000.
For the tSZ, we followed the approach of Bolliet et al. (2018) and developed a halo model to calculate the one-and 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. 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. (2019). Overall, C Planck,857GHz = 1.008 × C SPIRE,857GHz . 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. (2019) and is also presented in Sec. 3.4. 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. (2013a) 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.
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 Tab. C.1. In this case, instead of assuming a fixed perfect calibration for the Planck experiment, we let the calibration factors at all the frequencies vary with a Gaussian prior centred at 1.00 with 1σ error of 0.05. The Planck χ 2 value improves compared to when we fit for both the Planck and Herschel data together. The posterior values of calibration parameters are very close to one, and it is therefore a justified assumption to keep the Planck calibration fixed while fitting for both the data together.
Similar to the previous test, we then fit our model to Herschel data alone. The results are shown in Tab. C.2. Unlike the Planck data, we obtain a poor χ 2 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 Tab. 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 χ 2 values for both Planck and Herschel. However, it comes at the cost of unrealistically high values of the calibration parameters at 545 and 857 GHz that are completely inconsistent with the very precise 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 Tab. C.5 and Fig. C.3. In this case, we obtain a decent χ 2 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.