Impact of polarised extragalactic sources on the measurement of CMB B-mode anisotropies

One of the main goals of Cosmology is to search for the imprint of primordial gravitational waves in the CMB polarisation field, to probe inflationary theories. One of the obstacles toward the detection of the primordial signal is to extract the B-mode polarisation from astrophysical contaminations. We present a complete analysis of extragalactic foreground contamination due to polarised emission of radio and dusty star-forming galaxies. We update or use up-to-date models that are validated using the most recent measurements. We predict the flux limit (confusion noise) for the future CMB space or balloon experiments (IDS, PIPER, SPIDER, LiteBIRD, PICO), as well as ground-based experiments (C-BASS, NEXT-BASS, QUIJOTE, AdvACTPOL, BICEP3+Keck, BICEPArray, CLASS, SO, SPT3G, S4). Telescope aperture size (and frequency) is the main characteristic impacting the level of confusion noise. Using the flux limits and assuming constant polarisation fractions for radio and dusty galaxies, we compute the B-mode power spectra of the three extragalactic foregrounds (radio source shot noise, dusty galaxy shot noise and clustering), discuss their relative levels and compare their amplitudes to that of the primordial tensor modes parametrized by the tensor-to-scalar ratio r. At the reionization bump (l=5), contamination by extragalactic foregrounds is negligible. At the recombination peak (l=80), while the contamination is much lower than the targeted sensitivity on r for large-aperture telescopes, it is at comparable level for some of the medium- and small-aperture telescope experiments. For example, the contamination is at the level of the 68 per cent confidence level uncertainty on the primordial r for the LiteBIRD and PICO space experiments. Finally we also provide some useful unit conversion factors and give some predictions for the SPICA B-BOP experiment. Abridged


Introduction
The ΛCDM model is the standard model of cosmology. It is the simplest parametrization of the Hot Big Bang model, with two principal ingredients: Λ refers to a cosmological constant (i.e. the energy density of the vacuum) and CDM stands for cold dark matter, i.e. dark matter particules moving slowly. Being very successful in predicting a wide variety of observations related to the cosmic microwave background (CMB), the large-scale structure, and gravitational lensing, the ΛCDM model has reached the status of a paradigm. In this paradigm, an era of early exponential expansion of the Universe, dubbed cosmic inflation, has been proposed to explain why the Universe as revealed by the CMB radiation is almost exactly Euclidean and so nearly uniform in all directions. While the basic ΛCDM model fits all the data (with parameters known at the per cent level), the physics of inflation is still unknown. Thus, one of the central goals of modern cosmology is to determine the nature of inflation. One generic prediction is the existence of a background of gravitational waves, which produces a distinct, curl-like, signature in the polarisation of the CMB, referred to as primordial B-mode polarisation (which are due to tensor perturbations). The detection of this pri-mordial B-mode polarisation would provide a definitive proof that inflation occurred in the Early Universe. Unfortunately, cosmic inflation does not provide a unique prediction for the amplitude of the primordial tensor modes parametrized by the "tensorto-scalar ratio" r. We are in a situation where there is no natural range for r, in particular there is no relevant lower bound. The natural goal is to be able to measure r beyond doubt for the Higgs inflation (which is an inflationary scenario where the inflaton field is the Higgs boson), i.e. r≥∼ 2×10 −3 at 5σ. If this does not lead to a detection, it will discard the whole class of "large field" models. If the inflaton field was nothing but the Higgs field, this would have tremendous consequences for physics. Thus a precise measurement of (or upper bound on) r is essential for constraining inflationary physics. Present 95% CL upper limit on r as measured by Planck 1 combined with ground-based CMB experiments is r < 0.056 (Planck Collaboration et al. 2018c) at a pivot scale of k= 0.002/Mpc. The search for primordial B-mode is an outstanding challenge motivating a number of experiments designed to measure the anisotropies of the CMB in polarisation with an ever increasing precision.
B-modes are also generated by gravitational lensing of E-mode polarisation, providing a unique window into the physics of the evolved Universe and invaluable insights into late-time physics, such as the influence of dark energy and the damping of structure formation by massive neutrinos. These lensing B-modes are a nuisance for the primordial Bmodes. Several approaches have been studied for the CMB B-mode delensing, using large-scale structure surveys (galaxies or the Cosmic Infrared Background, e.g. Smith et al. 2012;Sherwin & Schmittfull 2015;Manzotti et al. 2017), or assuming that the lensing potential can be estimated internally from CMB data (e.g. Carron et al. 2017;Sehgal et al. 2017).
In addition to instrumental challenges, future experiments targeting r ∼10 −3 will have to face the critical problem of component separation. In addition to lensing, polarised Galactic foreground contamination dominates by several orders of magnitude the amplitude of the large-scale CMB B-modes. The capabilities of future experiments to clean the contamination due to polarised Galactic emissions have been investigated by e.g., Errard et al. (2016); Remazeilles et al. (2016); Philcox et al. (2018). Here, we investigate the polarisation fluctuations caused by the extragalactic contaminants, namely radio galaxies and dusty starforming galaxies (DSFG). While polarised compact extragalactic sources are expected to be a negligible foreground for CMB B-modes near the reionization peak (ℓ <10), they are expected to be the dominant foreground for r = 10 −3 once delensing has been applied to the data, from the recombination peak to higher multipoles, ℓ >50 (Curto et al. 2013).
Extragalactic radio sources are typically assumed to be Poisson distributed in the sky. The clustering of radio sources is strongly diluted by the broad distribution in redshifts of objects that contribute at any flux density. The contribution of clustering to the angular power spectrum is therefore small and can be neglected if sources are not subtracted down to very faint flux limits, S ≪ 10 mJy (González-Nuevo et al. 2005).
For DSFG, we have to consider polarisation fluctuations not only for the Poisson distribution of point sources but also for the clustering, i.e., the cosmic infrared background (CIB) anisotropies (Planck Collaboration et al. 2014b;Viero et al. 2013). The CIB power spectrum can be represented as the sum of two contributions usually called the 1-halo and 2-halo terms. The 1-halo represents the correlation of galaxies in the same dark-matter halo (pairs of galaxies inside the same halo); the 2-halo, capturing the galaxy correlations in different dark-matter haloes, describes the large-scale clustering. While we expect some polarisation fluctuations from the 1-halo (which is close to Poisson fluctuations), polarised 2-halo fluctuations are expected to be null, provided there is no correlation of the polarisation of galaxies within distinct halos. We could have a contribution from the large-scale clustering because of galaxy spin alignments in the filamentary dark-matter structure (e.g. Codis et al. 2018;Piras et al. 2018, and references therein). However, as recently shown by Feng & Holder (2019), this contribution is >100 and 1000 times lower than the shot noise of DSFG at ℓ=100 and ℓ=1000, respectively. Thus we consider that it has a negligible impact given its extreme weakness. and funded by Denmark, and additional contributions from NASA (USA).
In this paper we compute the expected level of polarised fluctuations from the shot noise of radio galaxies and DSFG and from the CIB 1-halo using up-to-date or updated models for a large set of future CMB space or balloon experiments (IDS, PIPER, SPIDER, LiteBIRD, PICO), as well as ground-based experiments (C-BASS, NEXT-BASS, QUIJOTE, AdvACTPOL, BICEP3+Keck, BICEPArray, CLASS, SO, SPT3G, S4). Our predictions use a point source detection limit self-consistently computed for each experiment (taking into account the sensitivities and determining confusion noises using our number count models). We also include some predictions for SPICA B-POP. An accurate computation of the flux detection limit is mandatory for the prediction of the shot noise of radio sources, as changing the flux cut by 30% affects the shot noise by 30%, while it is less important for DSFG as a small variation in the flux cut leads to only a small variation in shot-noise power (Planck Collaboration et al. 2011c).
Our work extends over previous studies that concentrate either on a single experiment or a restricted frequency area (e.g. Bonavera et al. 2017b;Remazeilles et al. 2018), or on a given galaxy population (e.g. radio galaxies Puglisi et al. 2018) or on high multipoles (e.g. Gupta et al. 2019 for ℓ 2000; e.g. Datta et al. 2019 for CMB EE). One of the originalities resides in using our radio and DSFG models, in combination with the CIB and CMB contamination and instrument noise, to iteratively predict the confusion noise due to extragalactic sources for all experiments.
The paper is organised as follows. We present the evolutionary models for radio sources and DSFG and discuss their polarised emission in Sect. 2 and 3. In Sect. 4, we give the formalism to compute polarised shot noise from galaxy number counts in intensity. We then describe our halo model of CIB anisotropies that is used to compute the polarisation power spectra due to the clustering of DSFG (Sect. 5). We use these models to compute the flux limit (caused by the fluctuations of the background sky brightness below which sources cannot be detected individually, i.e. the confusion noise) for a large number of future CMB experiments and for SPICA B-POP (Sect. 6). The flux limits allow us to compute the expected level of radio and dusty galaxy polarised shot noises that we discuss (together with the polarised 1-halo) in Sect. 7.1, and that we compare to the CMB primordial B-mode power spectrum in Sect. 7.2, for all experiments. We finally conclude in Sect. 8.

Number counts at cm to mm wavelengths
Number counts of extragalactic radio sources are well determined at radio frequencies ν < ∼ 10 GHz down to flux densities of S < ∼ 1 mJy (and even S < ∼ 0.03 mJy at 1.4 GHz) thanks to deep and large area surveys (see de Zotti et al. 2010, for a review). Moving to higher frequencies, i.e. from tens of GHz to mm wavelengths, observational data on radio sources are mainly provided by CMB experiments. Space missions like WMAP and Planck, which cover the full sky, was able to detect only bright sources, with flux densities larger than few hundreds of mJy at best. On the other hand, ground-based experiments, thanks to a better angular resolution, go deeper in flux density but only in small areas of the sky. Uncertainties on number counts are therefore still large, especially in the frequency range where CMB dominates, i.e. between 70 and 300 GHz.
Evolutionary models for extragalactic radio sources (e.g., Toffolatti et al. 1998;de Zotti et al. 2005;Massardi et al. 2010) are able to provide a good fit to data on luminosity functions and multi-frequency source counts from ∼ 100 MHz and to > ∼ 5 GHz. They adopt a schematic description of radio source populations, divided in steep-and flat-spectrum (or blazars) sources, according to the spectral index of the power-law spectrum, S (ν) ∝ ν α , at GHz frequencies that is lower or larger than −0.5. A simple power law is also used for extrapolating spectra to high frequencies, ν ≫ 5 GHz. However, especially for blazars, real source spectra are generally more complex than a power law that can hold only for limited frequency ranges. As a consequence, these models tend to over-predict the number counts of radio sources at ν > ∼ 100 GHz, as e.g. measured by the Atacama Cosmology Telescope (ACT) at 148 GHz (Marriage et al. 2011) or by Planck in all the High Frequency Instrument (HFI) channels (Planck Collaboration et al. 2011b. The main reason for this disagreement is the spectral steepening observed in Planck radio source catalogues above ∼ 70 GHz (Planck Collaboration et al. 2011b,a, 2016b and already suggested by other data sets (González-Nuevo et al. 2008;Sadler et al. 2008).
A first attempt of taking such steepening in blazar spectra into account has been done by Tucci et al. (2011). In that work the spectral behaviour of blazars at cm-mm wavelengths was statistically described by considering the main physical mechanisms responsible for the emission. In agreement with classical models of the synchrotron emission in the inner jets of blazars (Blandford & Königl 1979;Konigl 1981;Marscher & Gear 1985), the spectral high-frequency steepening was interpreted as caused, at least partially, by the transition from the optically-thick to the optically-thin regime. The frequency ν M at which the spectral break occurs depends on the relevant physical parameters of AGNs: the redshift, the Doppler factor (δ), and the linear dimension of the region (approximated as homogeneous and spherical) that is mainly responsible for the emission at the break frequency. In particular, Tucci et al. (2011) showed that the break frequency can be written in an approximated form as where D L is the luminosity distance of the sources and C is a function of the spectral indices before and after the break frequency (α f l and α st respectively) and of the flux density S ν 0 at a reference frequency (typically 5 GHz; see their Appendix B). Finally, the parameter r M is the distance from the AGN core of the jet region that dominates the emission at the frequency ν M (for a conical jet model, this parameter can be easily related to the dimension of the emitting jet region). It defines the dimension and thus the compactness of the emitting region at that frequency. This is the most critical parameter for determining ν M because of the large uncertainty on its actual value. Based on 5 GHz number counts and on information of spectral properties of radio sources at GHz frequencies, the Tucci et al. (2011) model provided predictions of number counts at cm/mm wavelengths by extrapolating flux densities of radio sources from low (1-5 GHz) to high frequencies. The model considered three populations of radio sources (steep-, invertedand flat-spectrum sources), and a different high-frequency spectral behaviour for each of them. Here we focus on blazars, which are the dominant class at ν > ∼ 70 GHz. The most successful model studied in the paper (referred as "C2Ex") assumes different distributions of the break frequency for BL Lac objects and Flat Spectrum Radio Quasars (FSRQs). According to it, most of FSRQs should bend their otherwise flat spectra between 10 and 100 GHz, whereas in BL Lac spectral breaks should be typically observed at ν > ∼ 100 GHz (implying that the observed synchrotron radiation comes from more compact emitting regions with respect to FSRQs). This dichotomy has been indeed found in the Planck radio catalogues (Planck Collaboration et al. 2011b, 2016b. This model provides a very good fit to all the data of bright (S > ∼ 100 mJy) radio sources for number counts and spectral index distributions up to ∼500-600 GHz (Planck Collaboration et al. 2011b. However, the agreement is only partial when other surveys, deeper in flux than Planck, are considered. In Fig. 1 we compare number counts from the model with observational data at frequencies between 70 and 220 GHz. Beyond Planck, data are from ACT (150, 218 GHz;Marsden et al. 2014) and SPT (95, 150, 220 GHz;Mocanu et al. 2013b). We see that the model tends to underestimate SPT/ACT counts in the flux-density range [20,60] mJy. More recently, Datta et al. (2019) measured new differential source number counts from the ACTPol two-season catalogs. Their new measurements confirm the underestimate of the C2Ex model in the flux density range [20,300] mJy.

Updated model for number counts
The recent data from ACT and SPT experiments give us the opportunity to better constrain the model parameters for blazars. We have seen above that the break frequency depends on a set physical parameters related to AGNs. In Tucci et al. (2011), most of them are imposed on the basis of observational constraints (as the redshift distribution of the different radio source populations; the doppler factor; spectral indices) and on the basis of typical assumptions for AGN model (equipartition condition, narrow conical jets, etc.). The only free parameter used in the model is the distance r M to the AGN core of the emitting jet region at the break frequency. In the best model of Tucci et al. (2011) r M is taken to be log-uniformly distributed in the range of [0.3, 10] pc for FSRQs and of [0.01, 0.3] pc for BL Lacs.
We now find the best estimate of the r M range by fitting number counts from Planck, ACT (Marsden et al. 2014) and SPT (Mocanu et al. 2013b) between 70 and 220 GHz. This is done only for BL Lacs, while for FSRQs we maintain the same range of r M values as before. We have verified in fact that a change in the r M interval for this class of objects does not give any significant improvement of the fit of number counts at sub-Jy level (i.e., for ACT/SPT data). This is not surprising because FSRQs provide the dominant contribution to number counts of bright sources, with S ≫ 100 mJy (see Fig. 1). At these flux levels the strong constraints come from Planck measurements that are already well described by the model. On the other hand, at fainter fluxes we see that the relevance of BL Lacs increases, and we expect they become the dominant population at few tens of mJy. This is exactly the range of fluxes where the model underestimates observational number counts. By increasing the contribution of BL Lacs we should remove or reduce the discrepancy between model and SPT/ACT data, without affecting the predictions for the very bright sources.
Jointly with r M , we consider the spectral index α st of blazars after the break frequency (i.e., in the optically thin regime) as an extra free parameter in the fit. In Tucci et al. (2011), it was taken Gaussianly distributed around α st = −0.8 with a dispersion of 0.2, in agreement with the canonical values for the optically-thin synchrotron spectral index. No differences between the two classes of blazars were considered. However, Planck Collaboration et al. (2011b, 2016b found that the aver- age spectral index of blazars after the spectral break is somewhat flatter than −0.8. The results of the fit give a more "compact" radio-emission regions in BL Lacs with respect to previous values, with 0.0025 ≤ r M ≤ 0.05 pc, i.e. about a factor 5 smaller than before. In addition, the average high-frequency spectral index is flatter, α st = −0.7, consistent with the trend observed in Planck data. Number counts predicted by the updated model differ mainly at low/intermediate flux densities, S < 0.1 mJy, and provide a clearly improved fit to observational data (see Fig. 1). The reduced χ 2 is now very close to 1. SPT data at 95 GHz are still slightly higher between 20 and 60 mJy, but the discrepancy is reduced and is not significant. Instead, the model provides a very good fit to the data at 150 and 220 GHz. To be noted that the change in the average value of α st produces a small increase also in number counts of FSRQs at ν > 100 GHz.

Statistical properties of polarised emission
Polarisation in radio sources is typically observed to be a few percent of total intensity at cm or mm wavelengths (e.g., (Murphy et al. 2010;Battye et al. 2011;Sajina et al. 2011;Massardi et al. 2013;Galluzzi et al. 2019), with only very few objects showing a fractional polarisation, Π = P/S , as high as ∼10%. Steep-spectrum radio sources are on average more polarised than flat-spectrum sources at ν < ∼ 20 GHz (Tucci et al. 2004;Klein et al. 2003), with a fractional polarisation that strongly depends on the frequency, from ∼ 2.5% at 1.4 GHz to ∼ 5.5% at 10.5 GHz (Klein et al. 2003). At low frequencies, flatspectrum sources instead are characterised by an almost constant and low degree of polarisation (∼ 2.5%).
Extensive studies of high-frequency polarisation properties have been done by Tucci & Toffolatti (2012) and Massardi et al. (2013) using the Australia Telescope 20 GHz (AT20G) survey (Murphy et al. 2010). This is a quite deep survey in intensity (with a completeness level of 91% at S ≥ 100 mJy and 79% at S ≥ 50 mJy in regions south of declination −15 • ) with a high detection rate in polarisation. Moreover, simultaneous measurements at 5 and 8 GHz are also available for a consistent fraction of objects. These analysis found that the distribution of the polarisation degree (in blazars) is well described by a log-normal function (see also Battye et al. 2011) with an average fractional polarisation of ∼ 3%. No clear correlation between the fractional polarisation and the flux density was observed, with a slight dependence on the frequency of the polarisation degree.
At frequencies ν > 20 GHz, polarisation measurements of very bright sources (S > ∼ 1 Jy) seem to indicate an increase of the fractional polarisation with frequency. Using the VLA for polarisation measurements of a complete sample of the WMAP catalogue, Battye et al. (2011) found that Π rad = 2.9, 3.0 and 3.5% at 8.4, 22, and 43 GHz respectively, and a fractional polarisation that is typically higher at 86 GHz than at 43 GHz. This was confirmed by measurements at 86 GHz from Agudo et al. (2010), obtained with the IRAM 30 m Telescope. They found that, for those sources with detected polarisation at 15 GHz, fractional polarisation at 86 GHz is larger than at 15 GHz by a mean factor of ∼ 2. However, these results were not confirmed using new data and/or improved data analysis procedures (Hales et al. 2014;Bonavera et al. 2017a;Galluzzi et al. 2017;Puglisi et al. 2018;Trombetti et al. 2018;Datta et al. 2019;Gupta et al. 2019). No significant trends of the polarisation degree with flux density or with frequency are found at the frequencies of interest for CMB B-mode search. Latest measurements of fractional polarisation at ν > 50 GHz vary from ∼1.5 to 3.5%, and are obtained either using log-normal fits to the distribution of observed polarisation fractions, or using stacking or statistical approaches.
To compute the radio source contamination in polarisation to the CMB B-mode (Sect. 7), we will assume a constant Π rad =2.8%, in agreement with the recent Planck, SPT and ACT measurements and radio source follow ups from 90 to 220 GHz.

Shot noise predictions
In this section, we compare the shot-noise level from residual radio sources found in observational data with values expected from our reference model, to check model validity. As the radio shot noise level is highly sensitive to the flux limit, we also provide some useful empirical relations that allow to compute the shot noise level as a function of the flux limit.

Shot noise levels in current CMB experiments
We report the residual shot-noise level in ACT and SPT data estimated by Dunkley et al. (2013); George et al. (2015), and compare them with predictions from Tucci et al. (2011) model, before and after our update in Table 1. The agreement is quite good for both cases, although the updated model provides a shot-noise level closer to the observational estimates.
In Table 2 we report auto-and cross-power spectra (shot noise only) due to residual radio sources in Planck data according to the updated model. We also compute the error of these predictions due to an uncertainty in the flux cut of 20 and 30 per cent. Moreover, we give a tentative estimate of the error associated to the uncertainty on the model, that is computed as the difference between results from the old and the updated model. The uncertainties we find are probably quite conservative, but, nevertheless, they are smaller than errors due to a 20% uncertainty in S cut at frequencies where radio sources are dominant (i.e. ν ≤ 217 GHz).
The consistency between the measured Poisson amplitude in the Planck auto-and cross-power spectra at 100, 143 and 217 GHz with the updated model discussed here was already investigated in Planck Collaboration et al. (2015, see their Table 20). The agreement is good, with the exception at 100 GHz where the predicted amplitude is significantly lower than the observed value. However, this discrepancy was attributed by the authors to a residual unmodelled systematic effect in the data rather than to a foreground modelling error. Moreover, the Poisson power at 100 GHz was found to be smaller in Planck Collaboration et al. (2014a), in better agreement with the model prediction.

Shot-noise level as a function of flux limits
It can be useful to know the dependence of the shot-noise level from residual radio sources with the flux cut S lim . We start considering the Planck frequencies, and a range of flux limits between 1 mJy and 1 Jy, i.e. more or less the range covered by CMB experiments.
We first start by auto-power spectra. We know that differential number counts for radio sources scale approximately as n(S ) ∝ S −2 , and power spectra as C ℓ ∝ S lim . Therefore, it is convenient to consider the quantity D S N = C ℓ /S lim . At a given frequency, we fit D S N ≡ D S N (S lim ) as a double power law: D S N (S lim ) from the updated model and the best fits given by Eq. 2 are shown in Fig. 2. The parameters of the fits are provided in Table 3. Cross-power spectra depend on the flux cuts at the two frequencies considered. In order to describe C ν 1 , ν 2 ℓ ≡ C ν 1 , ν 2 ℓ (S ν 1 lim , S ν 2 lim ) we choose to use a sixth-degree polynomial function. After computing cross-power spectra in an uniform grid of log(S lim /Jy) between −3 and 0, we find the polynomial fit using the IDL routine SFIT. For arbitrary flux limits (but always between 1 mJy and 1 Jy) at frequencies ν 1 and ν 2 , cross-power spectra can be estimated by means of where K i, j are the coefficients of the fit. We have verified that the fit has a typical error of 2-3%, with maximum errors of about 10-15% (usually at the borders of the grid). Figure 2 also shows examples of cross-power spectra and the corresponding fits fixing S ν 1 lim .

Model for number counts
For DSFG, we use the model of Béthermin et al. (2012b) which is based on the main assumption of two star formation modes in star forming galaxies: main sequence (MS) and starburst (SB). Main sequence galaxies are secularly evolving galaxies with a tight correlation between stellar mass (M ⋆ ) and star formation rate (SFR) at a given redshift. The evolution of MS and SB galaxies are based on Sargent et al. (2012) formalism, which used jointly the mass function of star forming galaxies, the redshift evolution of the sSFR (specific star formation rate, sSFR = SFR/M ⋆ ), and its distribution at fixed M ⋆ , with a separate contribution from MS and SB galaxies to reproduce infrared luminosity functions. The model uses a characteristic spectral energy distribution (SED) template for MS and SB based on fits of Draine & Li (2007) models to Herschel observations of distant galaxies as presented in Magdis et al. (2012). Despite its simplicity, the model provides one of the best fits achieved so far to the number counts from the mid-IR to radio wavelengths, including counts per redshift slice in the SPIRE bands, which were poorly reproduced by the previous generations of models.  Table 1. Shot-noise power of residual radio sources, D ℓ = ℓ(ℓ + 1)C ℓ /2π [µK 2 CMB ], at ℓ = 3000, estimated in ACT and SPT data, and predicted by models. Figure 2. (Left panel) Power spectra (divided by the flux limit) of residual radio sources as a function of the flux limit from 30 to 857 GHz (from top to bottom). Points are from the updated T11 model; solid lines are fits using Eq. 2 with parameters given in Table 3. (Right panel) Cross-power spectra at the frequencies indicated in the figure as a function of the flux limit S ν 2 lim (S ν 1 lim is equal to 0.4 Jy for ν 1 = 30, 70 GHz and 0.1 Jy for ν 1 = 100, 143 GHz). Solid lines are obtained from Eq.3.
We show on Fig. 3 the comparison of the model with the far-IR/sub-millimeter counts. The model agrees with the faint source counts measured by Herschel and SCUBA2 and the Euclidian plateau at high flux measured by Planck (Planck Collaboration et al. 2013a). For bright fluxes ( 1 Jy), the redshift grid of the model is too coarse to estimate properly the Euclidian plateau. We thus computed directly the value of the plateau using Eq. 6 of Planck Collaboration et al. (2013a). Note that an updated model has been published in Béthermin et al. (2017). This model has been used to build a 2 deg 2 simulation, and does not give any analytical predictions. It shows that Herschel measured number counts are overestimated in the range [10,100] mJy. This brings Béthermin et al. (2012b) model in better agreement with the data than can be seen in Fig. 3 at 250 and 350 µm. The volume of the simulation is too small to derive accurate predictions for the large volume surveys discussed here. We thus decided to use the 2012 analytical model, which agrees with the updated model at the ∼20% level.

Polarised emission
Little is known on the polarisation emission of dusty galaxies. Dust enshrouding star-forming galaxies absorbs UV radiation from stars, and re-emits light at longer wavelengths, being responsible for the far-IR SED of CIB galaxies. Thermal emission from interstellar dust in CIB galaxies, as in our Galaxy, is polarized due to the alignment of the dust grains with interstellar magnetic fields. The degree of polarisation is not very well known; it is likely to be low because the complex structure of galactic magnetic fields with reversals along the line of sight and the disordered alignment of dust grains reduce the global   Table 2. Auto-and cross-power spectra due to residual radio sources for Planck according to the updated model, for the flux cuts reported in the Table. polarisation signal from DSFG, although median values from Bonavera et al. (2017b) are consistent with their upper limits. Trombetti et al. (2018) derived a 90% confidence upper limit at 353 GHz, Π 2.2%. The upper limit at the same confidence level is looser at 217 GHz, Π 3.9%, where dusty galaxies are substantially fainter. Accordingly, to study the contamination from polarised emission of DSFG to the CMB B-modes (Sect. 7), we will adopt Π IR =1%.

Shot noise predictions
Béthermin et al. (2017) show that counts obtained from single dish antenna observations in the far-IR to the millimetre are biased high due to source multiplicity and clustering in the large beams (10 to 30 arcsec). This may cause important discrepancies between shot noises measured from the integral of the observed number counts and shot noises measured from CIB power spectra. For Herschel/SPIRE, there is an other complexity in the comparison introduced by the variation of the beam profile and aperture efficiency across the passband, giving relative spectral response function (RSRF) different for point sources and extended emission. To compare model predictions to shot noise measurements from CIB power spectra, we thus also run the model with extended RSRF. Comparison between model and observations are given in Table 4 and 5 for Herschel/SPIRE and Planck/HFI, respectively. The shot noise levels from observations are obtained either by fitting the CIB power spectra using the halo model (Viero et al. 2013;Planck Collaboration et al. 2014b) or by fitting the total power spectra using a parametric model and assuming a power law for the CIB (Mak et al. 2017).
In the first case, there is a strong degeneracy between the 1-halo term and the shot noise, especially at the Planck angular resolution. It is very difficult to derive any conclusion from Tables 4 and 5 as, i) some measured values are incompatible (i.e. when the shot noise derived with a higher flux limit is smaller than that derived with a lower flux limit). This is the case for Planck at 545 and 353 GHz and for Herschel at the three wavelengths and ii) the model is not systematically higher or lower than the measurements. In the frequencies of interest (ν 500 GHz), observations and model predictions agree by 20%, which we assume to be the uncertainty in our prediction. We stress that contrary to the radio, a small variation in the flux limit S lim leads to only a small variation in shot-noise power. For example, changing S lim by 30% leads to variation of the shot-noise level seen by Planck by less than 1% at 217 GHz (Planck Collaboration et al. 2011c)

Polarised shot noise from point sources: formalism
One can wonder why we should expect a polarisation term if galaxies have random orientations. To answer that question, we can define the complex linear polarisation of a source with flux S , where Π is the fractional polarisation and ψ the polarisation angle. If the polarisation angles of different sources are uncorrelated then < P s >= 0 , but the variance is non zero: We derive the shot-noise fluctuations of polarised point sources following Tucci et al. (2004). For Poisson distributed sources, the temperature power spectrum follows We can consider a similar expression for the polarisation power spectrum,  Table 4. Herschel/SPIRE shot noise levels as measured from CIB anisotropies and predicted using the integral of the number counts as modelled by Béthermin et al. (2012b). Values for the shot noise are given with the photometric convention νI ν =cst, obtained using either the "Point Source" or "Extended" RSRF (see text for more details).   Table 9, corrected to νI ν =constant and corrected from the calibration difference between PR1 and PR2 data releases (at 545 and 857 GHz). At 217 GHz, the contribution from radio sources has also been removed. where P = Q 2 + U 2 and C P Since the emission will contribute equally to EE and BB on average, we can consider The power spectrum due to sources with a given fractional polarisation is assuming that Π is not varying with S. Considering the distribution of fractional polarisation for all sources, the power spectrum becomes where is P(Π) is the probability density function of fractional polarisation.
This formulation is very convenient, as C P ℓ is defined as a function of a flux cut derived in total intensity. Thus it assumes that sources are removed from polarisation maps using total intensity data. This is the case with current CMB experiments and will probably also be most likely the case with future CMB data with the use of higher angular resolution and sensitivity surveys to remove the source contamination. With this formulation, one can also consider different source populations with different fractional polarisations.
The probability density function can be constrained from the observed distributions of fractional polarisations. However, due to the lack of constraints at CMB frequencies (∼90-200 GHz) for radio and dusty galaxies, we will consider a fix polarisation fraction for each population (see Sect.7).

Clustering of dusty star-forming galaxies
For the computation of polarisation power spectra due to the clustering of CIB galaxies we use the halo model, which provides a phenomenological description of the galaxy clustering at all relevant angular scales (Cooray & Sheth 2002). Assuming that all galaxies are located in virialized dark-matter halos, the CIB clustering power spectrum is expressed as the sum of two components: a 1-halo term, accounting for correlations between galaxies in the same halo, and a 2-halo term, due to correlations between galaxies belonging to separated dark matter halos. The former term, together with the shot-noise power spectrum, dominate the small-scale clustering, and the latter is prominent at large angular scales. Thus, the total CIB angular power spectrum at frequencies ν and ν ′ can be written as: In the following section, after briefly introducing the model and its main parameters (we refer the reader to Shang et al. (2012); Viero et al. (2013); Planck Collaboration et al. (2014b) for a detailed discussion), we show that the amplitudes of CIB polarisation power spectra are, at most, a small fraction of the 1-halo term of the clustering spectra, and we will derive upper limits on these amplitudes by fitting the model to current measurements of CIB angular power spectra from Herschel/SPIRE (Viero et al. 2013).

A halo model with luminosity dependence
In Limber approximation (Limber 1954), the CIB clustering power spectrum at frequencies ν and ν ′ is: where the term χ(z) denotes the comoving distance at redshift z, while a(z) is the scale factor. The total emissivity from all CIB galaxiesj ν (z) is computed from the luminosity function dn/dL as: where the galaxy luminosity L ν(1+z) is linked to the observed flux S ν as: Finally, the term P νν ′ (k, z) is the 3D power spectrum of the emission coefficient, expressed as: This term includes the 2-halo and 1-halo term. Expressing the luminosity of central and satellite galaxies as L cen,ν(1+z) (M H , z) and L sat,ν(1+z) (m S H , z) (where M H and m S H denote the halo and sub-halo masses, respectively), Eq. 14 can be written as the sum of the contributions from central and satellite galaxies as: Here dN/dm and dn/dm denote the halo and sub-halo mass function, from Tinker et al. (2008) and Tinker et al. (2010), respectively, and N cen is the number of central galaxies inside a halo, which will be assumed equal to zero if the mass of the host halo is lower than M min = 10 11 M ⊙ and one otherwise. Introducing f cen ν and f sat ν as the number of central and satellite galaxies weighted by their luminosity as the 3-dimensional CIB power spectrum at the observed frequencies ν, ν ′ in Eq. 16 can be expressed as the sum of 1-halo term and 2-halo term as: where The term u(k, M, z) is the Fourier transform of the halo density profile (Navarro et al. 1997) with a concentration parameter from Duffy et al. (2010), and b(M, z) denotes the halo bias (Tinker et al. 2010). The linear dark matter power spectrum P lin (k) in Eq. 21 is computed using CAMB (http://camb.info/). The parametrization of the term L (1+z)ν (M, z) is the key ingredient of the model. Following Shang et al. (2012), we assume a simple parametric function to describe the link between galaxy luminosity and its host dark matter halo, where the dependence of the galaxy luminosity on frequency, redshift, and halo mass is factorized in three terms as: The free normalization parameter L 0 is constrained by the data and has no physical meaning. The galaxy Spectral Energy Distribution (SED) is modelled as (see Blain et al. 2003, and reference therein): where the Planck function B ν has an emissivity index β = 1.5, (Planck Collaboration 2014; Serra et al. 2016). The power-law functional form at frequencies ν ≥ ν 0 has been already used in a number of similar analyses (Hall et al. 2010;Viero et al. 2013;Shang et al. 2012;Planck Collaboration et al. 2014b), and is more in agreement with observations than the exponential Wien tail. The free parameter T d is the mean temperature of the dust in CIB galaxies, averaged over the redshift considered range. We assume a redshift-dependent, global normalization of the L-M relation of the form and we consider a log-normal function to describe the luminosity-mass relation, as: The term σ L/m (fixed to σ L/m = 0.5, as in Shang et al. 2012;Viero et al. 2013;Planck Collaboration et al. 2014b;Serra et al. 2016) accounts for the range of halo masses mostly contributing to the infrared luminosity. The parameter M e f f describes the existence of a narrow range of halo masses around M e f f ∼ 10 12 M ⊙ associated with a peak in the star-formation efficiency, and due to various mechanisms suppressing star formation at both high and low halo masses (Benson et al. 2003;Silk 2003;Bertone et al. 2005;Croton et al. 2006;Dekel & Birnboim 2006;Béthermin et al. 2012b;Behroozi et al. 2013).

Results
We constrain the main parameters of our halo model using six measurements of CIB angular auto-and cross-power spectra at 250, 350, 500 µm from Herschel/SPIRE (Viero et al. 2013) in the multipole range 200 < l < 23000, and assuming the "extended" flux limit case. To further constrain the model, we also compute the star formation rate density in the range 0 < z < 6, and we fit to the latest compilation of star formation rate density measurements from Madau & Dickinson (2014). We perform a Monte Carlo Markov Chain (MCMC) analysis of the parameter space using a modification of the publicly available code CosmoMC (Lewis & Bridle 2002), and varying the following set of four halo model parameters: together with 6 free parameters A i=1,...6 for the amplitudes of the shot-noise power spectra. We obtain a good fit to the data, with a total χ 2 of 104.9 for 97 degrees of freedom. Mean values and marginalized limits for all free parameters used in the fit and comparison between Herschel/SPIRE measurements of the CIB power spectra with our best estimates of the 1-halo, 2-halo and shot-noise, are shown in Serra et al. (2016). Shot noises derived from this model are very close to those found for Béthermin et al. (2017) simulations. This gives us confidence on the level of the 1-halo term.

CIB power spectrum in polarisation
The polarisation fraction Π can be expressed, for a given intensity of dust emission I, in terms of the Stokes parameters Q and U as: where Q and U are related to the polarisation angle ψ, through: Polarisation power spectra can be computed with the same formalism used to compute the CIB intensity power spectrum, by substituting the galaxy luminosity L (1+z)ν (M, z) for Q and U as, respectively: (32) It is easy to see that, if the polarisation among different sources is uncorrelated (see Sect. 1), the 2-halo term cannot produce any polarisation power spectrum, since its computation involves an average over the polarisation angle of all sources, which is zero.
The contribution from the 1-halo term is slightly more complicated. In fact, the dark-matter halos mostly contributing to the CIB power spectra have mass in the range 10 12 < M H < 10 13 M ⊙ , and the typical number of satellite galaxies in this range is too small to average to zero the quantities F Q (1+z)ν (M, z) and F U (1+z)ν (M, z). As a result, when computing the 1-halo contribution, it is possible that terms proportional to give a positive contribution to the polarisation power spectra. Note that we do not consider here the terms proportional to f sat f cent since it has been shown, both in simulations and observationally, that the tidal field of a large central galaxy can torque its satellites such that the major axis of satellite galaxies points towards their hosts (see e.g., Fig. 8 in Pereira et al. 2008 or Fig. 6 in Joachimi et al. 2015) and we thus do not expect any polarised signal. While accurate estimates of the amplitude of the polarisation power spectrum would require numerical simulations, in this paper we estimate the maximum contribution from the 1halo term, and we show that it is almost negligible with respect to the contribution from the shot noise (see Sect.7.1). The maximum amplitude of polarisation can be obtained assuming the (unphysical) case where the polarisation angle ψ of all sources is perfectly correlated and equal to zero (for Q) or π/2 (for U).
Assuming Π IR the mean fractional polarisation of all DSFG, it is easy to see that the maximum amplitude of the polarisation power spectra is simply Π IR 2 times the amplitude of the 1-halo contribution to the CIB intensity power spectrum, keeping only the term proportional to f sat ν (M, z) 2 . Thus, the EE of BB CIB power spectra are computed following: Note that scaling the total CIB power spectrum using a fractional polarisation Π to derive the polarised CIB power spectrum (as done in Trombetti et al. 2018) overestimates the CIB contamination.

Future CMB experiments
We considered all future CMB experiments, either already selected, funded or in advanced discussion. Their name, frequency, angular resolution, sky coverage and instrument noise (in intensity) are given in Table 6 for balloon-and space-based experiments and Table 7 for ground-based experiments, respectively. We also considered Planck for reference and for cross-checks   6. CMB space and balloon experiments. From left to right: experiment name, frequency, angular resolution, sky fraction, and instrument noise (σ P inst , in polarisation). The standard deviations (σ) in mJy give the contributions of instrument noise, radio and dusty (IR) galaxies, CIB clustering and CMB to the total noise σ tot when measuring a point source flux (in intensity). They have been corrected for the flux lost by the aperture photometry procedure. S lim is the point source flux limit. SN radio and SN IR are the radio and dusty galaxy shot noises, respectively, corresponding to a flux cut equal to S lim .

Unit conversions and bandpass corrections
In the milimetre wavelength domain, two different units are often used. While for studies of Galactic emission or extragalac- tic sources, Janskies (Jy) are used, K CMB is the natural unit for CMB. Transforming Jy to K CMB is not only a unit conversion but also requires a colour correction (to account for the different spectral energy distribution that is implicitly assumed in the two units). This transformation is detailed in Appendix A. The conversion factors that are given in Tables B.1 and B.2 assume a square bandpass, with a δν and a central frequency ν given in the tables. Colour corrections are not computed for each experiment as it requires to know precisely the bandpasses (e.g. for Planck, assuming a square bandpass rather than the true bandpass leads to error in the colour corrections that are of the same order of the correction). Consequently all the numbers given in the Tables in Jy are given for the true spectra (but σ inst and σ CMB which are given for the convention νIν =constant, using the square bandpasses).
For current experiments with known bandpass, accurate unit conversions are given in Appendix A. Also for current experiments, comparison of foreground levels (CIB and SZ especially) necessitate their extrapolation between nearby frequencies of different experiments. To this end, useful conversion factors are given in Appendix A.

Confusion noise and flux limit
As seen in Eq. 11, we chose to use a flux cut in total intensity rather than in polarised intensity, for mainly two reasons: i) we assume that sources are removed or masked from polarisation maps using total intensity data, for which we could have a high resolution survey complete to some level in total intensity, as opposed to the equivalent in polarised intensity (e.g. Battye et al. 2011;Datta et al. 2019) ; ii) source number counts in polarisation are very scarce and more polarisation data are required to constrain dN/dP. By contrast, thanks to the numerous data in intensity obtained in the last decade, accurate modelling is available for number counts in intensity.
Consequently we compute in this section the confusion noise and flux limit in intensity, for each CMB experiment listed in Sect. 6.1.

Method and validation
The confusion noise 2 is usually defined as the fluctuations of the background sky brightness below which sources cannot be detected individually. These fluctuations are caused by intrinsically discrete extragalactic sources. In the far-IR, sub-millimetre and millimetre, due to the limited size of the telescopes compared to the wavelength, the confusion noise is an important part of the total noise budget. In fact, the confusion noise is often greater than the instrument noise, and is thus severely limiting the surveys depth.
When measuring a point source flux, the root mean square (rms) fluctuations due to extragalactic point sources is the sum of three components: where σ S Nrad , σ S Nir , and σ Clus are the rms fluctuations associated with the radio shot noise, dusty galaxy shot noise, and dusty galaxy clustering, respectively (we recall that clustering from radio sources is neglected, see Sect. 1). They are related to the power spectrum P k following where W k is the power spectrum of the beam (we assume Gaussian beams) and i stands for S Nrad, S Nir, and Clus, respectively. T k is the transfer function linked to the flux measurement of the sources. We assume that fluxes are measured using aperture photometry, where is the rectangular function, and R 1 and R 2 are the radii of the two circular apertures (with R 2 > R 1 ) and: The Fourier transform of f (r) is which gives the following power spectrum for our aperture photometry filter: The confusion noise can be determined using two criteria, the so-called photometric and source density criteria Lagache et al. 2003). The photometric case is derived from the fluctuations of the signal due to the sources below the detection threshold S lim in the beam. The source density case is derived from a completeness limit and evaluates the density of the sources detected above the detection threshold S lim , such that only a small fraction of sources is missed because they cannot be separated from their nearest neighbour. The choice of the criterion depends on the shape of the source counts and the solid angle of the beam . The transition between the two is around 200 µm, depending on telescope diameters . In this paper, we have thus to use the photometric criterion.
The photometric criterion is related to the quality of the photometry of detected sources, the flux measured near S lim being severely affected by the presence of fainter sources in the beam. It is defined by the implicit equation, where q phot measures the photometric accuracy (we assume q phot =5) and S lim is the confusion limit. σ tot is defined as where σ 2 con f is given in Eq. 36 and σ inst is the instrument noise per beam (given in Tables 6 and 7). We also added the noise introduced by CMB fluctuations, σ CMB , which is given by Eq. 37 where we replaced P k by the power spectrum of the CMB. F is a correction factor that accounts for the flux lost by the aperture photometry procedure (which is not covering the entire beam size).
In the range of confusion limits of CMB experiments, only P S Nrad k and P S Nir k depend on S lim . They are derived following where dN/dS are given by the models described in Sect. 2.1.1 and Sect. 3.1 for radio and dusty galaxies, respectively.
Confusion noises and flux limits are given in Tables 6 and 7. They have been obtained using R 1 =FWHM/2 and R 2 = 2×R 1 .
We checked that our predictions give confusion noises in very good agreement with those measured by ISO/ISOPHOT, Herschel/SPIRE and Planck. For SPIRE, we obtain σ con f =6.4, 6.6, and 5.3 mJy/beam while Nguyen et al. (2010) measure 5.8±0.3, 6.3±0.4, and 6.8±0.4 mJy/beam at 250, 350 and 500 µm, respectively. For Planck, we compare our flux limit to the flux cuts given in the PCCS2 source catalog for 90% completeness (in the extragalactic zone) in Table 8. We have an overall agreement at better than ∼2σ. We can however notice that our flux cut is systematically below the PCCS2 flux limit for the highest frequencies (217, 353, 545 and 857 GHz). We checked that this underestimate can be easily explained by the cirrus contamination, which may be quite high in the extragalactic zone (covering |b| > 30 o ) and which is ignored in the present paper. Finally, we also checked our results for SPT, by substituting σ inst from SPT-3G to the SPT-SZ survey. Considering σ S PT −S Z inst = 2, 1.2 and 4 mJy, we obtain S lim = 11, 7.1 and 20.5 mJy, at 95, 150 and 220 GHz respectively, in very good agreement with Mocanu et al. (2013a) (see their Table 3, for 95% completeness limit).
The very good agreement with previous far-IR, submillimetre and millimetre experiments gives us confidence in our computations.

Contributions to the point source sensitivity
Ground-based experiments have a maximum frequency of 280 GHz. The contribution of the different components to the point source sensitivity mostly depends on the frequency and size of the telescope apertures.
The smallest telescopes, with size <1m (BICEP, CLASS, SO-SAT, CMB-S4-SAT) or the low-frequency telescopes (C-BASS, NEXT-BASS, QUIJOTE, with ν < 40 GHz) have quite poor angular resolutions. The contribution of radio sources dominates up to ∼10-15GHz, then the confusion noise from the CMB becomes dominant. If we can get ride of the CMB, the CIB clustering dominates the noise budget at the higher frequencies (ν >200 GHz). Instrument noise is always much lower than the astrophysical components.
As expected, larger aperture telescope gives lower flux limits, because the confusion noise is much smaller (and generally the instrument noise is smaller too). For larger aperture telescopes (AdvACTPOL, SO-LAT, SPT-3G, CMB-S4-LAT), the instrument noise is at the same order of magnitude as confusion noises. For ν >145 GHz, the dominant contribution to the σ tot comes from the shot noise of DSFG.
In space, telescopes have in general smaller apertures and instrument noise is always negligible compared to confusion noise. Confusion from CMB is always dominating, except at the highest frequencies (ν 300 GHz). Appart from the CMB, we have a large contribution from galaxy clustering above ∼150-200 GHz. PIPER, SPIDER and LiteBIRD have large S lim (>1Jy) that will consequently lead to a large contamination to the CMB-B mode measurements.

The case of B-POP
We also considered the SPICA B-POP polarised experiment, which is at shorter wavelength. B-POP will provide 100-350 µm images of linearly polarised dust emission with an angular resolution, signal-to-noise ratio, and dynamic ranges comparable to those achieved by Herschel images of the cold ISM in total intensity. The angular resolution of B-BOP at 200 µm will also be a factor 30 better than Planck polarisation data.
At those wavelengths, and with such a high angular resolution, only the shot noise of DSFG contributes to the confusion noise (σ con f ). Flux limits are about 0.4, 19.6 and 35.3 mJy at 100, 200 and 350 µm, respectively (see Table 9). This is sightly above the SPIRE/Herschel 350 µm flux limit due to the smaller telescope aperture. For a 5'×5' survey, confusion noise levels are reached in 40, 0.07, and 0.08 seconds at 100, 200 and 350 µm, respectively 3 . For a 1 Sq. Deg. survey, they are reached in 1.6 hours, 10 seconds, and 12 seconds. This shows that the 200 and 350 µm maps, even on large areas, will be severely limited in depth due to extragalactic confusion, contrary to shorter wavelengths, such as the 100 µm.
In polarisation, confusion from galaxies will seriously limit the sensitivity of the high-latitude Galactic polarimetric surveys. Assuming a fractional polarisation for DSFG Π IR =1% (see Sect. 7), we obtain a confusion noise in polarisation σ P con f = 2.8, 138, 250 µJy after masking all the sources detected in intensity at 100, 200 and 350 µm, respectively. Those σ P con f levels are reached in 114, 0.18, and 0.23 hours for a single pointing and 67 500, 111, and 135 hours for a 1 Sq. Deg. survey, at 100, 200 and 350 µm, respectively. Therefore, in polarization, confusion should never be reached at 100 µm, but could be reached for the deepest integrations at longer wavelengths.

Contamination of the CMB B-modes
In order to provide reliable predictions of the radio source and DSFG contamination to CMB anisotropy polarisation measurements, we now need to assume a fractional polarisation for each population of galaxies. For radio sources, at CMB frequencies, there are still few polarisation measurements and very scarce polarisation fraction measurements for the different types of radio  Table 9. Confusion noise, flux limit and DSFG shot noise level for the SPICA B-POP experiment. Figure 4. Ratio of shot noise and clustering (1-halo CIB anisotropies) for dusty galaxies at ℓ=80 for all CMB experiments (ℓ=80 corresponds to the recombination B-peak).
sources (see Sect. 2.2). Thus, we use a constant Π rad =2.8%, in agreement with the recent Planck, SPT and ACT measurements and radio source follow ups from 90 to 220 GHz. For DSFG, the situation is even worse (see Sect. 3.2) and polarisation properties are almost completely unexplored. We adopt Π IR =1%.

Polarised power spectra of the extragalactic components
We give in Tables B.1 and B.2 the level of radio (C rad ℓ ) and DSFG (C IR ℓ ) shot noise power spectra, and the clustering power spectra (C CIB ℓ ) for three multipoles (ℓ=80, 1000, 4000). We first compare on Fig. 4 the relative level of DSFG shot noise and clustering power spectra at ℓ=80. We recall that the clustering power spectra is somehow an upper limit as we estimated the maximum contribution of the 1-halo term (see Sect.5.3). We see that the ratio C IR ℓ / C CIB ℓ is mostly constant, and between 2 and 3 for 120 < ν < 700 GHz. At lower frequencies, it is much higher (from 4 to 30) and thus C CIB ℓ can be neglected. Consequently, we do not compute the clustering power spectra for frequencies ν ≤90 GHz. The ratio increases very slowly with ℓ, by up to ∼30% at ℓ=4000 and ν <400 GHz.
We then compare on Fig. 5 the level of the radio power spectra and DSFG+clustering power spectra as a function of frequency. As expected, the general trend is an increase of  Thus, the DSFG power spectra level is higher than that of radio galaxies at a frequency which decreases with telescope size: ∼280, 200 and 130 GHz, from small to large apertures. These results do not depend on the multipole (as C IRsn ℓ / C IRclus ℓ is weakly varying with ℓ).

Comparison with the CMB B-modes
We first illustrate the contaminations of extragalactic components to the CMB B-mode power spectrum at two frequencies, ∼220 GHz (Fig. 6) and 145 GHz (Fig. 7). At each frequency, we plot the power spectra for two different aperture telescopes to illustrate the turnover between radio/DSFG dominant contaminations. The CMB B-mode power spectrum has been calculated for the Planck 2018 cosmology (using TT,TE,EE+lowE+lensing+BAO and a pivot scale for r of 0.002 Mpc −1 , Planck Collaboration et al. 2018a).
We compare on Fig. 6 Planck at 217 GHz with LiteBIRD at 235 GHz. While we have an equal contamination of DSFG and radio galaxies for Planck, the power spectrum of radio galaxies is 10 times larger than that of DSFG for LiteBIRD (even if the frequency of 235 GHz is higher). It is at the same level of the r = 0.01 (r = 0.001) B-mode power spectrum for ℓ=160 (ℓ= 83). For Planck, the contamination is negligible compared to the last 95% CL upper limit r 0.002 < 0.056 (Planck Collaboration et al. 2018c). On Fig. 7, we show the  level of the extragalactic components for the ground-based S4-SAT and S4-LAT experiments. Contamination by radio sources is dominating for S4-SAT, at a level of r=1.7×10 −3 at ℓ=80. For S4-LAT, the dominant contamination comes from DSFG shot noise, at a level of r=6.2×10 −6 at ℓ=80.
We finally compute the equivalent tensor-to-scalar ratio (r eq ) of the total extragalactic contamination (radio galaxy shot noise, DSFG shot noise and clustering) for each individual frequency at given multipoles. We show on Fig. 8 the variation of r eq as a function of frequencies at the recombination B-peak, ℓ=80. Minimum r eq is reached for 90 ν 300 GHz depending on the experiment. Similarly to Fig. 5 (and see Sect 7.1), we can distinguish three cases according to the telescope aperture size: -Large-aperture telescopes: minimum contamination is at the level of r eq =5.7×10 −6 for SPT-3G at 95 GHz. For SO-LAT, AdvACT and S4-LAT, r eq is about 1.1-1.4×10 −5 at both 90-93 and 145-150 GHz. These levels are well below the targeted σ r of these experiments (by a factor of 30-400). -Medium-aperture telescopes: minimum contamination is at the level of r eq ≃10 −4 and is reached at ν ≃200 GHz. While this is ∼40 times higher than σ r for IDS alone, it is at the same level as σ r for PICO (Hanany et al. 2019). -Small-aperture telescopes: the contamination reaches a level of ∼4×10 −4 for S4-SAT, SO-SAT and BICEPArray at ∼220 GHz. It increases to 7-8×10 −4 for CLASS, PIPER and LiteBIRD at 217, 270 and 280 GHz, respectively. Finally, it is about 2.5×10 −3 for SPIDER at 150 GHz. The level of contamination is below the targeted σ r for the Bicep/Keck experiment, for which they project 0.002 < σ r < 0.006 by the end of the planned BICEP Array program, assuming current modelling of polarised Galactic foregrounds and depending on the level of delensing that can be achieved with higher-angular resolution maps from the South Pole Telescope (Hui et al. 2018). It is on the contrary at the level of the 68% confidence level uncertainty of σ r < 10 −3 for LiteBIRD (this σ r includes statistical, instrumental systematic, and Galactic foreground uncertainties, Matsumura et al. 2016).
Note that the above comparison between r eq and σ r is done considering independently each frequency for r eq , while σ r has been usually estimated for each experiment combining the whole set of available bands, and under specific assumptions (e.g., taking into account systematic effects or foregrounds residual impacts). Multi-frequency component separations should be able to decrease the level of extragalactic foreground contamination.
To get a complementary view, rather than comparing r eq with σ r , we could compare r eq with the equivalent instrument noise r inst eq computed independently at each frequency. We calculate r inst eq following: where σ P inst is the instrument noise in polarisation (given in Tables 6 and 7). We show in Fig. 9 the ratio of r eq and r inst eq for all frequencies and experiments. At least ten percent contamination (r eq /r inst eq ≥ 0.1), for 70 ≤ ν ≤ 250 GHz, is reached for BICEP at 95 and 150 GHz, CLASS, SO-SAT and SPIDER at 93 GHz, LiteBIRD from 78 to 140 GHz, S4-SAT from 85 to 155 GHz, and PICO from 75 to 129 GHz. Combining higher and lower frequencies to decrease the Galactic foreground residuals may also add more contamination. For example, for PICO, 0.9 ≤ r eq /r inst eq ≤ 2.7 for 21 ≤ ν ≤ 52 GHz and for S4-SAT, it can reach few tens for ν=30-40 GHz. Figure 8. Equivalent tensor-to-scalar ratio (r eq ) of the sum of the extragalactic foregrounds at the recombination B-peak, ℓ=80, for the different CMB experiments (r eq is computed for each individual frequency). Figure 9. Ratio of equivalent r of the extragalactic foregrounds (r eq ) and instrument noise (r inst eq ), at ℓ=80. Figure 10. Equivalent tensor-to-scalar ratio (r eq ) of the sum of the extragalactic foregrounds at ℓ=5, corresponding to the reionisation B-bump. Also shown are the σ r for LiteBIRD, PICO and CLASS (dashed lines). For Planck, the current 1σ upper limit is r < 0.028 and is thus not visible on the figure.
The scale dependency of extragalactic foregrounds compared to the CMB makes the ratio of the primordial CMB signal over foregrounds more favorable at larger scale, in particular at the reionization B-bump (ℓ=5). Only nearly full-sky ( f sky ≥70%) experiments can provide some measurements at such low multipoles. The r equivalent in this case are very small (1.8×10 −6 for PICO at 186 GHz, 3.6-3.1×10 −5 for LiteBIRD at 195-235 GHz, and 2.3×10 −5 for CLASS at 217 GHz; see Fig. 10). They are much smaller than the targeted limits on the primordial r for PICO and LiteBIRD, and σ r = 8.5 × 10 −3 for CLASS (including diffuse Galactic thermal dust and synchrotron foregrounds, Watts et al. 2015). For Planck, the contamination is much smaller than the current upper limit (Planck Collaboration et al. 2018a).
Finally, we can look at the ratio of the extragalactic foreground and CMB lensing BB power spectra (at ℓ=1000). This ratio is ∼120 times larger than the equivalent tensor-to-sclar ratio r eq at ℓ=80. It goes from ∼ 10 −3 for large-aperture experiments to ∼ 10 −1 for small-aperture experiments. As already known, ground-based large-aperture telescopes will provide the ability to delens the maps from future satellite CMB missions, such as LiteBIRD (e.g., Namikawa & Nagata 2014).

Conclusion
We have computed the expected level of polarised fluctuations from the shot noise of radio galaxies and DSFG and from the CIB clustering using up-to-date or updated models. Using those models, we predicted the point source detection limits (confusion noises, in intensity) for future CMB space or balloon experiments (IDS, PIPER, SPIDER, LiteBIRD, PICO), as well as ground-based experiments (C-BASS, NEXT-BASS, QUIJOTE, AdvACTPOL, BICEP3+Keck, BICEPArray, CLASS, SO, SPT3G, S4). Those limits have been computed by taking into account the instrument noise, the three extragalactic foregrounds, as well as the CMB. The models, as well as the point-source detection flux limits, have been validated using most recent measurements on number counts, CIB power spectra, confusion noises and shot noise levels. As expected, we found that the confusion noise levels are mostly driven by the telescope-aperture sizes and frequency. We also predict the confusion noise for SPICA B-BOP and showed that confusion will limit the sensitivity of polarised surveys (with the confusion noise in polarisation reached in 0.18, and 0.23 hours for a single pointing at 200 and 350 µm, respectively).
Assuming a constant polarisation fraction for the radio and dusty sources of Π rad =2.8% and Π IR =1%, respectively, we then predicted the shot noises and CIB 1-halo clustering B-mode power spectra. We compared the amplitude of the different extragalactic foregrounds as a function of frequency and telescopeaperture size. We found that CIB clustering is almost negligible. The relative levels of radio and DSFG shot noises are mainly driven by the telescope sizes, which can be classified in three categories: large-aperture (≥ 6m, i.e. SPT-3G, S4-LAT, SO-LAT, AdvActPol), medium-aperture (∼1.5m, i.e. Planck, IDS, PICO) and small-aperture (≤0.6m, i.e. LiteBIRD, SPIDER, CLASS, SO-SAT, S4-SAT, BICEP-Keck) telescopes. While we have an equal contribution between radio shot noise and DSFG shot noises (+ clustering) at ν ≃ 130 GHz for large-aperture telescopes, it reaches ν ≃280 GHz for small-aperture telescopes, which are thus dominated by the radio shot noise at the frequencies dedicated to the CMB measurement. González-Nuevo et al. (2005) show that the contribution of radio source clustering to the temperature angular power spectrum is small and can be neglected if sources are not subtracted down to very faint flux limits, S ≪ 10 mJy. However, future ground-based experiments such as S4-LAT will be able to reach flux limits of the order of 2-3 mJy. At these levels, the clustering of radio sources might be not negligible for ℓ < 30 but its expected level in polarisation will likely be much smaller than the polarised shot noise level, as it is the case for DSFGs.
We computed the equivalent tensor-to-scalar ratio (r eq ) of the total extragalactic contamination (radio galaxy shot noise, DSFG shot noise and clustering) for given multipoles. At the reionization B-bump (ℓ=5), the extragalactic contamination will not limit the measurements. At the recombination B-peak (ℓ=80), the contamination for large-aperture telescope experiments is much below the targeted primordial r, but this is not the case for some of the small-and medium-aperture telescopes. For example for the LiteBIRD and PICO space experiments, the contamination is at the level of the 68% confidence level uncertainty on the primordial r (not considering a multi-frequency component separation that should globally decrease r eq ). On the other side of the multipole range, extragalactic components represent 10% of the CMB lensing BB power spectrum at ℓ=1000 for LiteBIRD. Moreover, a similar slope is observed between the extragalactic components and the CMB lensing BB power spectrum up to ℓ=200 and between the extragalactic components and the primordial B-mode power spectrum for 15 ℓ 50, leading to degeneracies in any model fitting. Cleaning the data from this extragalactic contamination is thus mandatory for some of the small-and medium-aperture telescope experiments.
Foreground mitigation has been studied for the Galactic components. It has been shown that it requires a multi-frequency coverage (but see Philcox et al. 2018 for a method based on anisotropy statistics, or Aylor et al. 2019 for the use of Neural Network). This multi-frequency approach will be difficult to apply for extragalactic foregrounds, as the three extragalactic components are degenerated (i.e. same power spectra at the multi-pole of interest) and the sum of the three does not have a welldefined frequency dependency. Moreover, even if more precise polarised source counts for radio galaxies will be obtained in the near future, the variation of the radio shot noise with the flux limit (changing the flux cut by 30% affects the shot noise by 30%, see Table 2), together with the variability of radio sources, may prevent using more accurate modeling to precisely predict the shot noise level.
Polarised Galactic foregrounds are dominated by dust and synchrotron emissions with spatial variation of their SEDs. Using a parametric maximum-likelihood approach, Errard et al. (2016) found that combinations from ground, balloon and space experiments can significantly improve component separation performance, delensing, and cosmological constraints over individual datasets. In particular, they find that a combination of post-2020 ground-and space-based experiments could achieve constraints such as σ r ∼1.3×10 −4 after component separation and iterative delensing. However, such results (see also e.g. Stompor et al. 2016) are often derived ignoring complexities in the Galactic foreground emission due to synchrotron and dust, and neglecting potential other contaminants such as anomalous microwave emission and extragalactic foregrounds. Also, they are adopting component separation methods that essentially assume a model that well-matches the simulated foregrounds under study. Remazeilles et al. (2016) test some of these assumptions explicitly, and find biases in the derived value of r of more than 1σ by e.g. neglecting the curvature of the synchrotron emission law. Given their levels for some of mid-and small-aperture telescopes, extragalactic foregrounds have definitively to be considered in the component separation methods dedicated to the extraction of the CMB B-modes. For that purpose, our detailed computation of flux limits and shot-noise levels will allow to precisely include these foregrounds into the input sky models. 2) for dusty star forming galaxies are given for two different CIB spectral energy distributions ("model" refers to the model of Béthermin et al. (2012a) while "measure" refers to Gispert et al. (2000) fit of FIRAS measurements). For SPIRE, we give the colour corrections for the two spectral responses (extended or point-source RSRF). To convert an intensity in K CMB to an equivalent specific intensity MJy sr −1 , the original intensity has to be multiplied by the factors given in the Table. A.5. Converting CIB power spectra between HFI, ACT, SPT and SPIRE The purpose here is to convert the measurement through one bandpass into a measurement as it would be obtained through another bandpass (often close in frequency, e.g., HFI at 143 GHz versus SPT at 150 GHz). It means we want to find K such that For ease of reading, lets write I 1 and I 2 the fiducial monochromatic flux densities from spectral response 1 and 2  Table A.3. y S Z to K CMB unit conversion. To convert an intensity in K CMB to y S Z , the original intensity has to be multiplied by the factors given in the Table. (with the convention νI ν = constant) at their respective reference frequencies ν 1 and ν 2 . Combining Eq. A.1 and A.2 gives It then follows, where R 1 and R 2 are the normalised spectral responses 1 and 2, respectively. Values for K for HFI, ACT and SPT are given in Table A Table A.4. Factors to convert the CIB intensity (in Jy/sr with the convention νI ν =constant) in the HFI, ACT, and SPT bandpasses (see Eqs. A.3 and A.5). ν 1 and ν 2 are given in the first column and first line, respectively (e.g., K(ν 1 , ν 2 ) = K(857, 545) = 1.989. The factors have been computed using Béthermin et al. (2012a) CIB SED. Some factors can be deduced from combinations of other ones. We give all of them for convenience. Appendix B: B-mode power spectrum of the extragalactic foreground components We give in the following tables the values of the B-mode power spectrum of the extragalactic foreground components for all ex-periments, together with the unit conversion factor. Table B.1. C BB ℓ of the extragalactic foreground components for space and balloon experiments: radio galaxies, dusty galaxies (IR) and CIB 1-halo (completely negligible for ν ≤ 90 GHz and thus not computed). They are given in Jy 2 /sr. The unit conversion factor is also given (C = MJy sr −1 [νI ν = constant] K −1 CMB ). The power spectra in Jy 2 /sr has to be divided by C 2 to obtain power spectra in µK CMB .